基于COMSOL的平房倉冷卻通風過程中糧堆熱濕耦合傳遞研究

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

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

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

在冷卻干燥通風過程中，糧堆內部熱濕遷移是一個復雜的多場耦合過程[4]，糧堆內部熱濕遷移過程的影響因素很多，包括糧倉所在地氣象條件、糧堆物性參數和糧食生物特性等。冷卻干燥通風過程中糧堆熱濕耦合遷移規律研究涉及湍流、濕空氣傳熱、多孔介質傳熱、多孔介質傳質和環境大氣等多個物理場。

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

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

1 糧堆內熱濕耦合傳遞的COMSOL數值模擬

1.1 高大平房倉三維物理模型

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

圖1 高大平房倉的COMSOL三維物理模型

1.2 COMSOL軟件開發接口與參數設置

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

以實測的送風濕空氣的溫濕度擬合函數作為模擬研究的入口邊界條件。參數的正確設置對糧堆內熱濕耦合傳遞COMSOL 數值模擬十分重要。表1給出了COMSOL數值模擬時的相關條件設置和參數設置。

表1 COMSOL數值模擬的相關條件和參數設置

1.3 糧堆內熱濕耦合傳遞COMSOL數值模擬控制方程

冷卻干燥通風過程中糧堆內熱濕耦合傳遞COMSOL數值模擬需要濕空氣和多孔介質兩種材料，利用COMSOL內置材料方程將小麥糧堆物理參數賦值給多孔介質材料。COMSOL數值模擬涉及湍流、濕空氣傳熱、濕空氣傳質、多孔介質傳熱和多孔介質傳質等多個物理場，還需要熱濕耦合、溫度耦合和流動耦合多個耦合物理場。多孔介質熱濕耦合傳遞數值模擬需要描述谷粒吸附解析水分時的熱量交換和水分遷移，該過程的實現需將COMSOL 內置方程修改為糧堆內熱濕耦合傳遞控制方程。基于非穩態的多孔介質熱濕傳遞偏微分方程數值解，以含濕量、溫度和壓力為驅動勢，建立了冷卻干燥通風過程中糧堆內熱濕耦合傳遞動態數學模型。

1.3.1 糧堆內水分遷移控制方程

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

(1)

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

根據菲克定律和達西定律：

jv=-δpPv+jaxa

(2)

jl=KlPk

(3)

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

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

(4)

假設溫度對糧堆平衡含濕量的影響忽略不計，則有：

(5)

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

(6)

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

Pv=φPs

(7)

xa=6.2×10-6Pv

(8)

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

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

(9)

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

(10)

(11)

1.3.2 糧堆內熱量平衡控制方程

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

(12)

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

qcond=-kT

(13)

qconv=Cp,ajaT+jvhlv

(14)

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

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

(15)

1.3.3 糧堆內動量守恒控制方程

COMSOL數值模擬根據泊肅葉定律[20]，將糧堆內的空氣平均流動速度表示為壓力梯度與速度的關系式。

ja=-kaPa

(16)

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

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

(17)

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

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

(18)

1.4 計算機軟件和硬件配置

本研究采用AutoCAD軟件繪制高大平房倉三維幾何模型，然后把高大平房倉三維幾何模型導入COMSOL Multiphysics 5.3軟件，構建形成聯合體的高大平房倉COMSOL三維物理模型。采用COMSOL Multiphysics 5.3軟件模擬計算冷卻干燥通風過程中糧堆內熱濕耦合傳遞過程，COMSOL Multiphysics 5.3和Origin8.5軟件作為數值模擬的繪圖及數據分析軟件。

2 模擬結果與討論

數值模擬了冷卻干燥通風38 h糧堆內熱濕耦合傳遞過程，分析了糧堆溫度場、濕空氣速度場和濕空氣壓力場的變化情況，通過對比預測值與實測值，驗證了COMSOL數值模擬結果的正確性，進而得出了冷卻干燥通風過程中糧堆內熱濕耦合傳遞的規律。

2.1 數值模擬結果的準確性驗證

圖2為冷卻干燥通風過程中的糧堆平均溫度模擬預測值和實測值的變化關系。分析圖2可知，冷卻干燥通風過程中糧堆平均溫度模擬預測值與實測值吻合較好，兩者的最大溫度偏差為0.211 ℃，表明COMSOL數值模擬的溫度預測精度較高。

圖2 糧堆實測和模擬預測溫度平均值

圖3 糧堆4 m高度處實測溫度平均值與預測溫度平均值

利用糧堆4 m高度處溫度傳感器測量值驗證COMSOL數值模擬結果的準確性。圖3為糧堆4 m高度處實測溫度與預測溫度平均值。由圖3可知，冷卻干燥通風開始階段糧堆4 m高度處實測溫度平均值與預測溫度平均值之間有較大的溫差，隨著冷卻通風時間的增加實測溫度平均值與預測溫度平均值之間溫差逐漸減小，冷卻通風結束時糧堆4 m高度處預測溫度平均值已十分接近實測溫度平均值。這是因為數值模擬假定糧堆初始溫度為16.7 ℃的等溫體，導致冷卻通風開始時糧堆4 m高度處預測溫度平均值與實測溫度平均值之間有較大的溫差。隨著冷卻通風時間的增加，糧堆4 m高度處預測溫度平均值逐漸接近實測溫度平均值。這進一步證明COMSOL數值模擬具有較高的預測精度。

2.2 數值模擬結果分析與討論

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

圖4 通風70 h糧堆截面y=6 m的濕空氣壓力場云圖

圖5 通風70 h糧堆截面y=6 m的速度場云圖

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

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

圖6是通風過程中糧堆溫度場云圖。分析圖6可知，糧堆溫度場分布存在明顯的分層現象，而且基本上是對稱分布的。與糧堆的初始溫度相比，溫度下降較為明顯。糧堆平均溫度由初始平均溫度32 ℃降為14.74 ℃。糧堆底部的送風道周圍區域溫度最低，自糧堆底部向上糧堆溫度逐漸升高，未被冷卻的糧層溫度最高。相鄰兩個風道之間糧堆溫度下降速度相對較慢，糧堆底部兩個風道之間的中心區域比其周圍區域的溫度高，主要原因是該溫度較高區域處于通風死角，通過該區域的冷卻干燥通風濕空氣量相對較少，該區域的體積大小跟糧層厚度、送風風速、送風道間距和送風道類型等因素有關，可以通過糧倉通風系統優化來減小通風死角區域的體積。

圖6 通風70 h糧堆截面y=6 m的溫度場云圖

圖7 冷卻干燥通風過程中糧堆干基水分的預測值

圖7給出了冷卻干燥通風過程中糧堆干基水分變化圖。分析圖7可知，糧堆干基水分由初始水分16.7%降為14.5%，表明冷卻干燥通風可以干燥小麥糧堆。在冷卻干燥通風前期階段水分下降速度較快，隨著通風時間的增加糧堆干基水分下降速度不斷減小，分析水分遷移控制方程可知，小麥顆粒的水分與顆粒空隙間濕空氣的相對濕度和溫度之間存在對應關系。隨著通風時間的增加，糧堆干基水分和溫度不斷降低，而糧食顆粒內水蒸氣分壓力與濕空氣水蒸氣分壓力之間的壓力差不斷減少，因此糧堆干基水分的下降速度不斷減小，直至小麥糧堆水分與濕空氣水分達到近似平衡水分狀態。小麥糧堆達到平衡水分狀態后，繼續進行冷卻干燥通風，糧堆干基水分將不再發生明顯變化。

3 結論

3.1 通過多場耦合有限元軟件COMSOL 的二次開發接口，可以實現冷卻干燥通風過程中糧堆內熱濕耦合傳遞的數值模擬。通過修改COMSOL 材料內置方程，使多孔介質材料具有了小麥糧堆的物性參數。通過將COMSOL 內置方程修改為糧堆內水分遷移、熱量守恒和動量守恒控制方程，可以實現冷卻干燥通風過程中糧堆內熱濕耦合傳遞過程的數值模擬。糧堆孔隙率、滲透系數和密度等物理參數會影響數值模擬的預測精度，數值模擬時應該根據實際情況合理設置這些參數。

3.2 在糧倉冷卻干燥通風過程中，糧堆降溫和降水是同時存在的。糧堆干基水分的變化方向是逐漸趨于糧食的平衡水分狀態，糧堆干基水分的變化方向和變化速率主要跟糧堆干基水分值和通風的濕空氣參數有關。通風正上方區域糧堆降溫較為明顯，相鄰通風道之間的區域糧堆降溫效果較差，冷卻干燥通風過程中存在通風死角區域，通過糧倉通風系統優化可以消除或減小通風死角區域的體積，可以更好地保障儲糧安全。