陳桂香 張宏偉 王海濤 劉超賽 尹 君

(河南工業(yè)大學(xué)土木建筑學(xué)院1, 鄭州 450001)(國(guó)家糧食局科學(xué)研究院2, 北京 100037)

我國(guó)每年糧食總產(chǎn)量高達(dá)5億t，倉(cāng)內(nèi)害蟲、糧堆發(fā)熱和谷物霉變是引起儲(chǔ)藏過(guò)程中糧食損失的主要原因[1-2]。糧堆內(nèi)的溫度和水分是影響谷物霉變、糧堆發(fā)熱和倉(cāng)內(nèi)害蟲的重要因素[3]。糧倉(cāng)冷卻干燥通風(fēng)可以把糧堆溫度降低至安全溫度以下，把糧堆水分減小到安全水分以下，是一種有效減少糧食在儲(chǔ)存過(guò)程中損失的技術(shù)措施。研究冷卻干燥通風(fēng)過(guò)程中糧堆內(nèi)熱濕耦合傳遞規(guī)律，對(duì)指導(dǎo)糧倉(cāng)的通風(fēng)冷卻和改善糧食儲(chǔ)存環(huán)境具有重要意義。

在冷卻干燥通風(fēng)過(guò)程中，糧堆內(nèi)部熱濕遷移是一個(gè)復(fù)雜的多場(chǎng)耦合過(guò)程[4]，糧堆內(nèi)部熱濕遷移過(guò)程的影響因素很多，包括糧倉(cāng)所在地氣象條件、糧堆物性參數(shù)和糧食生物特性等。冷卻干燥通風(fēng)過(guò)程中糧堆熱濕耦合遷移規(guī)律研究涉及湍流、濕空氣傳熱、多孔介質(zhì)傳熱、多孔介質(zhì)傳質(zhì)和環(huán)境大氣等多個(gè)物理場(chǎng)。

數(shù)值模擬作為一種經(jīng)濟(jì)有效的研究手段，在探索糧堆內(nèi)部的熱濕傳遞規(guī)律過(guò)程中日益受到學(xué)者的重視和關(guān)注[5]。Chang等[6-7]提出了一個(gè)預(yù)測(cè)小麥通風(fēng)狀態(tài)下溫度和水分變化的數(shù)學(xué)模型，并據(jù)此對(duì)糧堆熱濕傳遞規(guī)律進(jìn)行了數(shù)值模擬研究。尹君等[8-9]利用MATLAB模擬軟件重現(xiàn)了不同倉(cāng)型的小麥糧堆溫度場(chǎng)和水氣分壓場(chǎng)的變化過(guò)程，并研究了不同季節(jié)糧堆的“冷芯”“熱芯”和結(jié)露現(xiàn)象。陳桂香等[10]采用CFD軟件數(shù)值模擬研究了通風(fēng)過(guò)程中糧堆內(nèi)熱濕傳遞及霉變預(yù)測(cè)。李祥利等[11]通過(guò)數(shù)值模擬方法研究了通風(fēng)狀態(tài)下“圭”字形風(fēng)道的高大平房倉(cāng)糧堆內(nèi)部溫度以及水分的變化規(guī)律。王雪等[12]模擬研究了糧食溫度場(chǎng)在倉(cāng)內(nèi)部熱源作用下的分布及變化情況。張前等[13]采用回歸模型方法預(yù)測(cè)了高大平房倉(cāng)內(nèi)部糧堆的溫度變化。Daniela de Carvalho Lopes等[14]采用數(shù)值模擬軟件AERO研究了糧堆內(nèi)發(fā)熱點(diǎn)對(duì)糧堆一維傳熱的影響。Carrera-Rodríguez等[15]采用Fortran軟件研究了周圍環(huán)境溫度對(duì)高粱糧堆空隙內(nèi)熱環(huán)境的影響。

本研究建立了冷卻干燥通風(fēng)過(guò)程中高大平房倉(cāng)COMSOL三維物理模型，以實(shí)測(cè)的送風(fēng)空氣溫度和濕度為入口邊界條件，通過(guò)修改COMSOL內(nèi)置材料方程、耦合滲流、能量守恒、動(dòng)量守恒和水分遷移控制方程，數(shù)值模擬糧堆溫度場(chǎng)、濕空氣速度場(chǎng)、濕空氣壓力場(chǎng)和水分變化過(guò)程，得出了冷卻干燥通風(fēng)過(guò)程中糧堆內(nèi)熱濕耦合傳遞規(guī)律。

1 糧堆內(nèi)熱濕耦合傳遞的COMSOL數(shù)值模擬

1.1 高大平房倉(cāng)三維物理模型

以位于華東地區(qū)長(zhǎng)寬高為42 m×24 m×11 m的高大平房倉(cāng)為研究對(duì)象。該糧倉(cāng)采用六個(gè)“一機(jī)兩道”的壓入式通風(fēng)系統(tǒng)，采用“U”字型地籠通風(fēng)道。糧倉(cāng)裝備了一套無(wú)線傳感器的糧情監(jiān)控系統(tǒng)。糧倉(cāng)布置了83根測(cè)溫電纜，共安裝了四層332個(gè)溫度傳感器。高大平房倉(cāng)進(jìn)行冷卻干燥通風(fēng)時(shí)，糧倉(cāng)內(nèi)小麥堆糧高度為8 m。糧倉(cāng)冷卻干燥通風(fēng)時(shí)間為70 h。在COMSOL軟件中，建立了如圖1所示的高大平房倉(cāng)三維物理模型。整個(gè)糧倉(cāng)物理模型采用非結(jié)構(gòu)網(wǎng)格形式共劃分為350多萬(wàn)個(gè)網(wǎng)格。

圖1 高大平房倉(cāng)的COMSOL三維物理模型

1.2 COMSOL軟件開發(fā)接口與參數(shù)設(shè)置

COMSOL 為用戶提供了多種開發(fā)接口方式，其主要開發(fā)接口方式有：利用MATLAB 進(jìn)行二次開發(fā)、利用應(yīng)用程序接口API 進(jìn)行二次開發(fā)、利用物理模型創(chuàng)建器進(jìn)行二次開發(fā)和利用修改COMSOL 內(nèi)置方程進(jìn)行二次開發(fā)。用戶可以利用MATLAB 的數(shù)據(jù)處理功能、利用COMSOL 提供的API 函數(shù)、利用COMSOL 提供的創(chuàng)建器開發(fā)出自定義物理場(chǎng)模型或修改COMSOL 內(nèi)置方程，根據(jù)研究需要進(jìn)行COMSOL 軟件二次開發(fā)。在本研究中，COMSOL 數(shù)值模擬涉及的糧食物性參數(shù)和邊界條件可為任意的常數(shù)、變量或函數(shù)表達(dá)式等多種形式，因此本文主要采用兩種COMSOL 軟件開發(fā)接口形式：通過(guò)修改COMSOL 內(nèi)置方程進(jìn)行二次開發(fā)和利用應(yīng)用程序接口API 進(jìn)行二次開發(fā)。

以實(shí)測(cè)的送風(fēng)濕空氣的溫濕度擬合函數(shù)作為模擬研究的入口邊界條件。參數(shù)的正確設(shè)置對(duì)糧堆內(nèi)熱濕耦合傳遞COMSOL 數(shù)值模擬十分重要。表1給出了COMSOL數(shù)值模擬時(shí)的相關(guān)條件設(shè)置和參數(shù)設(shè)置。

表1 COMSOL數(shù)值模擬的相關(guān)條件和參數(shù)設(shè)置

1.3 糧堆內(nèi)熱濕耦合傳遞COMSOL數(shù)值模擬控制方程

冷卻干燥通風(fēng)過(guò)程中糧堆內(nèi)熱濕耦合傳遞COMSOL數(shù)值模擬需要濕空氣和多孔介質(zhì)兩種材料，利用COMSOL內(nèi)置材料方程將小麥糧堆物理參數(shù)賦值給多孔介質(zhì)材料。COMSOL數(shù)值模擬涉及湍流、濕空氣傳熱、濕空氣傳質(zhì)、多孔介質(zhì)傳熱和多孔介質(zhì)傳質(zhì)等多個(gè)物理場(chǎng)，還需要熱濕耦合、溫度耦合和流動(dòng)耦合多個(gè)耦合物理場(chǎng)。多孔介質(zhì)熱濕耦合傳遞數(shù)值模擬需要描述谷粒吸附解析水分時(shí)的熱量交換和水分遷移，該過(guò)程的實(shí)現(xiàn)需將COMSOL 內(nèi)置方程修改為糧堆內(nèi)熱濕耦合傳遞控制方程?；诜欠€(wěn)態(tài)的多孔介質(zhì)熱濕傳遞偏微分方程數(shù)值解，以含濕量、溫度和壓力為驅(qū)動(dòng)勢(shì)，建立了冷卻干燥通風(fēng)過(guò)程中糧堆內(nèi)熱濕耦合傳遞動(dòng)態(tài)數(shù)學(xué)模型。

1.3.1 糧堆內(nèi)水分遷移控制方程

COMSOL數(shù)值模擬將糧堆視為各項(xiàng)同性的連續(xù)多孔介質(zhì)。根據(jù)單元體質(zhì)量守恒，式(1)給出了糧堆內(nèi)水分遷移控制方程，用于描述糧堆內(nèi)熱濕耦合傳遞的水分遷移情況。

(1)

式中：ω為體積含濕量/kg/m3；t為時(shí)間/s；jv為水蒸氣傳遞速率/kg/(m2·s)；jl為液態(tài)水傳遞速率/kg/(m2·s)；Sw為水分源項(xiàng)[16]，Sw利用Thorpe[17]介紹的方法計(jì)算確定。

根據(jù)菲克定律和達(dá)西定律：

jv=-δpPv+jaxa

(2)

jl=KlPk

(3)

式中：δp為水蒸氣滲透率/kg/(m·s·Pa)；Pv為水蒸氣分壓力/Pa；ja為空氣流動(dòng)速率/kg/(m2·s)；xa為糧食顆粒間濕空氣含濕量/kg/kg；Kl為液態(tài)水滲透率/kg/(m·s·Pa);Pk為毛細(xì)水壓力/Pa。

將式(2)和式(3)代入式(1) 得：

(4)

假設(shè)溫度對(duì)糧堆平衡含濕量的影響忽略不計(jì)，則有：

(5)

式中：ζ是等溫吸放濕曲線的斜率/kg/m3；φ是糧食顆粒間濕空氣相對(duì)濕度，φ的開爾文關(guān)系式可以表示為：

(6)

糧食顆粒間濕空氣按理想氣體處理，則有：

Pv=φPs

(7)

xa=6.2×10-6Pv

(8)

式中：ρl為水的密度/kg/m3；RD為水蒸氣氣體常數(shù)/J/(kg·K)；T為開爾文溫度/K；Ps為飽和水蒸氣壓力(Pa)。

將式(5)～式(8)代入式(4)，可以得到：

(9)

式中：Dφ和DT可以分別利用式(10)和式(11)計(jì)算確定。

(10)

(11)

1.3.2 糧堆內(nèi)熱量平衡控制方程

COMSOL數(shù)值模擬假設(shè)糧堆內(nèi)的水分只有汽、液兩相[18]，根據(jù)單元體能量守恒，式(12)給出了對(duì)流傳熱過(guò)程熱量守恒偏微分方程。

(12)

式中：ρm是小麥顆粒的干基密度/kg/m3；ρm=790.3 kg/m3。Cp,m是糧堆的恒壓熱容/J/(kg·K)，Cp,m=1 934.2 (J/(kg·K))；qcond是導(dǎo)熱熱流密度/W/m2；qconv是對(duì)流熱流密度/W/m2；Sh是熱量源項(xiàng)[19]，Sh利用Thorpe[17]介紹的方法計(jì)算確定。

qcond=-kT

(13)

qconv=Cp,ajaT+jvhlv

(14)

式中：k是糧堆的導(dǎo)熱系數(shù)/W/(m·K)；Cp,a是干空氣的恒壓熱容/J/(kg·K)；hlv是水蒸氣汽化潛熱/J/kg。

將式(13)和式(14)代入式(12)得：

(15)

1.3.3 糧堆內(nèi)動(dòng)量守恒控制方程

COMSOL數(shù)值模擬根據(jù)泊肅葉定律[20]，將糧堆內(nèi)的空氣平均流動(dòng)速度表示為壓力梯度與速度的關(guān)系式。

ja=-kaPa

(16)

式中：ka是糧堆內(nèi)濕空氣的滲透率/m2；ka是指在空氣流動(dòng)方向上空氣流動(dòng)速率與壓力梯度的比值，其物理意義是指在一定壓差下糧堆允許濕空氣通過(guò)的能力；Pa是濕空氣壓力/Pa。

根據(jù)連續(xù)性方程，式(16)可以變形為

(17)

式中：ρa(bǔ)是濕空氣密度/kg/m3；ε是糧堆孔隙率，ε=0.4。

在糧堆內(nèi)熱濕耦合傳遞過(guò)程的COMSOL數(shù)值模擬研究中，糧堆內(nèi)的濕空氣流速低，濕空氣流動(dòng)的壓力低，濕空氣流動(dòng)過(guò)程中溫度變化小，因此糧堆內(nèi)濕空氣可以作為不可壓縮氣體處理，則式(17)可簡(jiǎn)化為：

(18)

1.4 計(jì)算機(jī)軟件和硬件配置

本研究采用AutoCAD軟件繪制高大平房倉(cāng)三維幾何模型，然后把高大平房倉(cāng)三維幾何模型導(dǎo)入COMSOL Multiphysics 5.3軟件，構(gòu)建形成聯(lián)合體的高大平房倉(cāng)COMSOL三維物理模型。采用COMSOL Multiphysics 5.3軟件模擬計(jì)算冷卻干燥通風(fēng)過(guò)程中糧堆內(nèi)熱濕耦合傳遞過(guò)程，COMSOL Multiphysics 5.3和Origin8.5軟件作為數(shù)值模擬的繪圖及數(shù)據(jù)分析軟件。

2 模擬結(jié)果與討論

數(shù)值模擬了冷卻干燥通風(fēng)38 h糧堆內(nèi)熱濕耦合傳遞過(guò)程，分析了糧堆溫度場(chǎng)、濕空氣速度場(chǎng)和濕空氣壓力場(chǎng)的變化情況，通過(guò)對(duì)比預(yù)測(cè)值與實(shí)測(cè)值，驗(yàn)證了COMSOL數(shù)值模擬結(jié)果的正確性，進(jìn)而得出了冷卻干燥通風(fēng)過(guò)程中糧堆內(nèi)熱濕耦合傳遞的規(guī)律。

2.1 數(shù)值模擬結(jié)果的準(zhǔn)確性驗(yàn)證

圖2為冷卻干燥通風(fēng)過(guò)程中的糧堆平均溫度模擬預(yù)測(cè)值和實(shí)測(cè)值的變化關(guān)系。分析圖2可知，冷卻干燥通風(fēng)過(guò)程中糧堆平均溫度模擬預(yù)測(cè)值與實(shí)測(cè)值吻合較好，兩者的最大溫度偏差為0.211 ℃，表明COMSOL數(shù)值模擬的溫度預(yù)測(cè)精度較高。

圖2 糧堆實(shí)測(cè)和模擬預(yù)測(cè)溫度平均值

圖3 糧堆4 m高度處實(shí)測(cè)溫度平均值與預(yù)測(cè)溫度平均值

利用糧堆4 m高度處溫度傳感器測(cè)量值驗(yàn)證COMSOL數(shù)值模擬結(jié)果的準(zhǔn)確性。圖3為糧堆4 m高度處實(shí)測(cè)溫度與預(yù)測(cè)溫度平均值。由圖3可知，冷卻干燥通風(fēng)開始階段糧堆4 m高度處實(shí)測(cè)溫度平均值與預(yù)測(cè)溫度平均值之間有較大的溫差，隨著冷卻通風(fēng)時(shí)間的增加實(shí)測(cè)溫度平均值與預(yù)測(cè)溫度平均值之間溫差逐漸減小，冷卻通風(fēng)結(jié)束時(shí)糧堆4 m高度處預(yù)測(cè)溫度平均值已十分接近實(shí)測(cè)溫度平均值。這是因?yàn)閿?shù)值模擬假定糧堆初始溫度為16.7 ℃的等溫體，導(dǎo)致冷卻通風(fēng)開始時(shí)糧堆4 m高度處預(yù)測(cè)溫度平均值與實(shí)測(cè)溫度平均值之間有較大的溫差。隨著冷卻通風(fēng)時(shí)間的增加，糧堆4 m高度處預(yù)測(cè)溫度平均值逐漸接近實(shí)測(cè)溫度平均值。這進(jìn)一步證明COMSOL數(shù)值模擬具有較高的預(yù)測(cè)精度。

2.2 數(shù)值模擬結(jié)果分析與討論

利用建立的高大平房倉(cāng)COMSOL三維物理模型，通過(guò)改變初始條件，模擬了糧堆高度相同的新糧糧堆冷卻干燥通風(fēng)過(guò)程。新糧糧堆的初始溫度為32 ℃，糧堆的初始干基水分為16.7%。冷卻干燥通風(fēng)的數(shù)值模擬時(shí)間為70 h。數(shù)值模擬研究冷卻干燥通風(fēng)過(guò)程中糧堆熱濕耦合傳遞過(guò)程，給出糧堆溫度場(chǎng)、濕空氣壓力場(chǎng)和濕空氣速度場(chǎng)變化情況。

圖4 通風(fēng)70 h糧堆截面y=6 m的濕空氣壓力場(chǎng)云圖

圖5 通風(fēng)70 h糧堆截面y=6 m的速度場(chǎng)云圖

圖4是冷卻干燥通風(fēng)過(guò)程中糧堆濕空氣壓力場(chǎng)云圖。由圖4 可知, 通風(fēng)地籠內(nèi)的濕空氣壓力最大，濕空氣壓力從糧堆底部往上逐漸減小。糧層阻力是造成這一現(xiàn)象的主要原因。濕空氣壓力分布存在明顯的分層現(xiàn)象，而且基本上是對(duì)稱分布的。

圖5 是通風(fēng)過(guò)程中糧堆濕空氣速度場(chǎng)云圖。由圖5 可知，通風(fēng)地籠內(nèi)的濕空氣速度最大，糧堆中濕空氣速度最小，糧堆上部濕空氣速度介于兩者之間，這一現(xiàn)象出現(xiàn)的原因也主要是糧層阻力的存在。在冷卻干燥通風(fēng)過(guò)程中，糧堆溫度場(chǎng)變化會(huì)受濕空氣速度的影響，而濕空氣速度分布又取決于其壓力的分布。

圖6是通風(fēng)過(guò)程中糧堆溫度場(chǎng)云圖。分析圖6可知，糧堆溫度場(chǎng)分布存在明顯的分層現(xiàn)象，而且基本上是對(duì)稱分布的。與糧堆的初始溫度相比，溫度下降較為明顯。糧堆平均溫度由初始平均溫度32 ℃降為14.74 ℃。糧堆底部的送風(fēng)道周圍區(qū)域溫度最低，自糧堆底部向上糧堆溫度逐漸升高，未被冷卻的糧層溫度最高。相鄰兩個(gè)風(fēng)道之間糧堆溫度下降速度相對(duì)較慢，糧堆底部?jī)蓚€(gè)風(fēng)道之間的中心區(qū)域比其周圍區(qū)域的溫度高，主要原因是該溫度較高區(qū)域處于通風(fēng)死角，通過(guò)該區(qū)域的冷卻干燥通風(fēng)濕空氣量相對(duì)較少，該區(qū)域的體積大小跟糧層厚度、送風(fēng)風(fēng)速、送風(fēng)道間距和送風(fēng)道類型等因素有關(guān)，可以通過(guò)糧倉(cāng)通風(fēng)系統(tǒng)優(yōu)化來(lái)減小通風(fēng)死角區(qū)域的體積。

圖6 通風(fēng)70 h糧堆截面y=6 m的溫度場(chǎng)云圖

圖7 冷卻干燥通風(fēng)過(guò)程中糧堆干基水分的預(yù)測(cè)值

圖7給出了冷卻干燥通風(fēng)過(guò)程中糧堆干基水分變化圖。分析圖7可知，糧堆干基水分由初始水分16.7%降為14.5%，表明冷卻干燥通風(fēng)可以干燥小麥糧堆。在冷卻干燥通風(fēng)前期階段水分下降速度較快，隨著通風(fēng)時(shí)間的增加糧堆干基水分下降速度不斷減小，分析水分遷移控制方程可知，小麥顆粒的水分與顆?？障堕g濕空氣的相對(duì)濕度和溫度之間存在對(duì)應(yīng)關(guān)系。隨著通風(fēng)時(shí)間的增加，糧堆干基水分和溫度不斷降低，而糧食顆粒內(nèi)水蒸氣分壓力與濕空氣水蒸氣分壓力之間的壓力差不斷減少，因此糧堆干基水分的下降速度不斷減小，直至小麥糧堆水分與濕空氣水分達(dá)到近似平衡水分狀態(tài)。小麥糧堆達(dá)到平衡水分狀態(tài)后，繼續(xù)進(jìn)行冷卻干燥通風(fēng)，糧堆干基水分將不再發(fā)生明顯變化。

3 結(jié)論

3.1 通過(guò)多場(chǎng)耦合有限元軟件COMSOL 的二次開發(fā)接口，可以實(shí)現(xiàn)冷卻干燥通風(fēng)過(guò)程中糧堆內(nèi)熱濕耦合傳遞的數(shù)值模擬。通過(guò)修改COMSOL 材料內(nèi)置方程，使多孔介質(zhì)材料具有了小麥糧堆的物性參數(shù)。通過(guò)將COMSOL 內(nèi)置方程修改為糧堆內(nèi)水分遷移、熱量守恒和動(dòng)量守恒控制方程，可以實(shí)現(xiàn)冷卻干燥通風(fēng)過(guò)程中糧堆內(nèi)熱濕耦合傳遞過(guò)程的數(shù)值模擬。糧堆孔隙率、滲透系數(shù)和密度等物理參數(shù)會(huì)影響數(shù)值模擬的預(yù)測(cè)精度，數(shù)值模擬時(shí)應(yīng)該根據(jù)實(shí)際情況合理設(shè)置這些參數(shù)。

3.2 在糧倉(cāng)冷卻干燥通風(fēng)過(guò)程中，糧堆降溫和降水是同時(shí)存在的。糧堆干基水分的變化方向是逐漸趨于糧食的平衡水分狀態(tài)，糧堆干基水分的變化方向和變化速率主要跟糧堆干基水分值和通風(fēng)的濕空氣參數(shù)有關(guān)。通風(fēng)正上方區(qū)域糧堆降溫較為明顯，相鄰?fù)L(fēng)道之間的區(qū)域糧堆降溫效果較差，冷卻干燥通風(fēng)過(guò)程中存在通風(fēng)死角區(qū)域，通過(guò)糧倉(cāng)通風(fēng)系統(tǒng)優(yōu)化可以消除或減小通風(fēng)死角區(qū)域的體積，可以更好地保障儲(chǔ)糧安全。