糧堆熱學特性的數(shù)值模擬

吳才章，程明福

（河南工業(yè)大學 電氣工程學院，河南 鄭州 450001）

0 引言

糧食是民生的根本，是關系到社會經(jīng)濟穩(wěn)定發(fā)展的核心物資.糧食儲藏是糧食產(chǎn)業(yè)的核心環(huán)節(jié)，我國糧食儲藏體系龐大，糧庫企業(yè)數(shù)量高達2萬多個，儲備糧總共有1.5億 t左右.長期以來，糧庫儲糧中容易出現(xiàn)倉房內(nèi)積熱多、糧堆局部發(fā)熱等問題，從而導致局部乃至整體糧堆中微生物和害蟲的大量繁殖，造成糧食品質(zhì)劣化以及數(shù)量損耗，由此給國家造成巨大的損失.因此，對糧堆熱物理學特性進行研究，進而通過非生物因素（傳熱、調(diào)質(zhì)等方法）干預控制糧堆的溫濕度變化，確保儲糧安全，成為當前糧食儲藏技術研究的一個熱點[1-3].

對糧食熱學特性的研究，國外開展得相對較早，系統(tǒng)地測定并研究了各種糧食的比熱和導熱系數(shù)及變化規(guī)律；國內(nèi)相關的研究起步較晚，科研基礎較薄弱[4-5]，研究主要針對糧庫糧堆溫度場的演化，采用宏觀的研究方法，基于現(xiàn)場直接測量的數(shù)據(jù)，代入散粒體熱學方程，通過擬合的方式獲得糧堆溫度場的變化，但由于缺乏對我國主要糧食品種的熱物理學特性的系統(tǒng)研究，對糧庫儲藏技術的指導有限.作者從微觀的角度出發(fā)，研究糧食熱物理學基本參數(shù)，測定糧食的比熱和導熱系數(shù)，進而模擬糧庫糧堆溫度場的演化，為糧庫通過傳熱、調(diào)質(zhì)等方法干預控制糧堆局部發(fā)熱奠定基礎.同時，該課題的研究也為我國糧食濕熱物理學參數(shù)測試標準的制定提供可靠的依據(jù)，為糧食儲藏技術改進提供重要的理論依據(jù).

1 糧堆導熱系數(shù)測定

FD-TC-B是基于穩(wěn)態(tài)法來測量不良導體的導熱系數(shù).其基本過程是，利用熱源對樣品加熱，樣品內(nèi)部的溫差使熱量從高溫向低溫處傳導，適當控制試驗條件和參數(shù)使加熱和傳熱的過程達到平衡狀態(tài)，待測樣品內(nèi)部可形成穩(wěn)定的溫度分布，根據(jù)溫度分布可以計算出導熱系數(shù)[6-7].

待測樣品小麥為散粒體，在試驗中需要將小麥均勻的包在很薄的纖維網(wǎng)（網(wǎng)孔直徑1.9 mm）中，并形成要求的形狀.注意使測試樣品直徑與散熱盤一致，而且樣品厚度不應過大，否則會造成較大誤差.試驗中準備4組樣品：豫麥49（A）、鄭麥366（B）、周麥16（C）和矮抗58（D）.單次測量的試驗步驟為：

（1）準備樣品，測量其厚度及直徑.將樣品放入加熱盤與散熱盤中間，然后固定，調(diào)節(jié)微調(diào)螺絲使其與加熱盤和散熱盤接觸良好.

（2）加熱樣品，打開電扇開關.加熱盤的溫度上升到設定溫度值（一般設定為75～80℃）時，開始記錄散熱盤的溫度，每隔1 min記錄一次，直到加熱盤和散熱盤的溫度在一個足夠長的時間段內(nèi)基本保持不變，即認為傳熱已經(jīng)達到穩(wěn)定狀態(tài)，記錄樣品上下表面的溫度分別為θ1和θ2.相關的參數(shù)測試為：環(huán)境溫度為16℃；散熱盤的質(zhì)量m為891.4 g，比熱容c為385 J/(kg·K），半徑RP為65.0 mm，厚度hP為7.66 mm；樣品的半徑dB為65.0 mm，厚度hB為18.0 mm.4種樣品A、B、C、D水分含量分別為12.1%、11.8%、11.6%、12.2%.各樣品的測試結果如表1所示.

表1 樣品的上下表面穩(wěn)態(tài)溫度

（3）取走樣品后加熱散熱盤，使其溫度在原溫度基礎上增加20℃左右，然后移去加熱盤，開啟風扇讓散熱盤冷卻，每隔10 s記錄該盤溫度.作散熱曲線，計算散熱盤的冷卻速率，散熱盤的冷卻曲線，如圖1所示.當傳熱達到穩(wěn)定狀態(tài)時，可以認為加熱盤通過樣品傳遞的熱流量與散熱盤向周圍環(huán)境的散熱量相等，通過散熱盤在穩(wěn)態(tài)溫度θ2時的散熱速率來求熱流量ΔQ/Δt.

圖1 散熱盤冷卻曲線

（4）將試驗測得樣品的上、下表面穩(wěn)態(tài)溫度θ1和θ2，熱流量ΔQ/Δt以及樣品和散熱盤的相關參數(shù)代入導熱系數(shù)計算公式，計算所測樣品的導熱系數(shù)λ.

式中：λ為樣品的導熱系數(shù)；hB為樣品的厚度，mm；dB為樣品的半徑，mm；m為散熱盤的質(zhì)量，g；c為散熱盤的比熱容，J/(kg·K）；Rp為散熱盤的半徑，mm；hp為散熱盤厚度，mm.

由圖1可知，豫麥49的樣品下表面穩(wěn)態(tài)溫度θ2為38.6℃時，散熱盤的冷卻速率將其代入導熱系數(shù)計算公式，可得豫麥49的導熱系數(shù)為λA=0.161 8 J/（m·s·K）.同理，λB=0.159 5 J/（m·s·K），λC=0.159 6 J/（m·s·K），λD=0.158 3 J/（m·s·K）.

2 糧堆熱學方程

2.1 糧堆的導熱方程

在圓柱坐標下，根據(jù)傅里葉定理，糧堆的導熱方程：

式中：ρ為糧食的密度；C為糧食的比熱；λ為糧食的導熱系數(shù)，J（/m·s·K）；T為糧食的溫度，℃；r為柱體的半徑，m；t是時間，s；φ是柱體內(nèi)任意點在其截面上的角度；Z是柱體的軸向尺寸；qν為內(nèi)熱源.

糧堆是由有生命、非剛性、膠體顆粒組成的非均質(zhì)散粒體，糧倉內(nèi)糧堆的傳熱是一個復雜的非穩(wěn)態(tài)傳熱過程，影響糧堆傳熱過程的因素較多，譬如環(huán)境的溫濕度、糧食物性、糧倉的形狀和結構等.為方便研究，將做適當簡化，假設：

（1）糧食為各向同性的均勻介質(zhì)；

（2）空氣及糧食熱物性參數(shù)近似不變；

（3）在初始時刻，糧堆內(nèi)溫度分布一致.

根據(jù)假設，在圓柱坐標下糧堆的導熱方程可以簡化為：

2.2 邊界條件與初始條件的選取

在圓筒倉中，糧堆有3個導熱邊界，如圖2所示.

圖2 圓筒倉中糧堆的導熱邊界

糧堆頂部邊界S1，由于糧堆與頂部空氣有對流換熱，它屬于第三類邊界條件：

式中：h為空氣和糧堆表層的對流換熱系數(shù)；Ta為頂層空氣的溫度；T為糧堆表層的溫度；n為外法線方向.

糧堆四周邊界S2，由于倉壁較厚，又考慮隔熱設計，它屬于第一類邊界條件：

式中：Tb為糧倉內(nèi)壁面的溫度.

糧堆底部邊界S3，糧倉底部一般深入地下，一年四季保持恒溫，其絕熱邊界條件為：

初始時刻，設糧堆整體的溫度一致，即

3 糧堆溫度場的數(shù)值模擬

3.1 模擬對象及參數(shù)選取

以河南省鄭州市某糧食儲備庫中的立筒倉為研究對象，其立筒倉倉高12 m，堆糧線11 m，內(nèi)徑8 m，糧倉內(nèi)儲藏小麥，小麥品種為鄭麥366.糧庫所在地具有年溫差較大，空氣干燥等氣候特點.相關的熱物性參數(shù)如表2所示，鄭麥366的導熱系數(shù)為0.159 5 J（/m·s·K）.

表2 小麥相關的熱物性參數(shù)

3.2 溫度場的MATLAB模擬及求解

偏微分方程數(shù)值解方法有拉格朗日元法、有限元法、離散元法、有限差分法等.文中采用有限元法，即利用部分插值把區(qū)域連續(xù)求解的微分方程離散成解線性代數(shù)方程組，然后進行求解.對于前面糧堆非穩(wěn)態(tài)導熱偏微分方程，利用MATLAB的PDE工具箱進行求解，求解過程如圖3所示，以圖示的方式表達糧堆的溫度分布.

4 結果及分析

按照求解偏微分方程的步驟，對糧堆的溫度場求解需要進行幾何模型設定、方程設定、初始條件和邊界條件代入、網(wǎng)格劃分及求解方程.糧倉上層空間氣溫設為40℃，糧倉墻壁內(nèi)側(cè)和糧堆溫度均為15℃，經(jīng)過求解得到糧堆的溫度分布圖.計算結果表明，由于糧倉墻壁具有較強的隔熱效果，糧堆受四周溫度影響變化不明顯.糧堆頂部由于存在空氣和糧堆對流換熱，糧堆頂部熱量交換較顯著，造成糧堆上層溫度自上而下呈梯度分布.1 h和10 h后糧堆頂部溫度隨縱深距離（Z軸）變化情況，如圖4所示.

圖3 利用MATLAB的PDE工具箱求解偏微分方程步驟

圖4 經(jīng)過1 h和10 h后糧堆頂部溫度隨縱深距離變化曲線

由圖4可知，1 h后，糧堆溫度從40℃到15℃自上而下的縱深僅為2 cm；經(jīng)過10 h后，糧堆溫度從40℃到15℃自上而下的縱深為12 cm，這一結論也適用于熱量從四周向糧堆中間傳遞時溫度隨傳遞深度的變化特性.數(shù)值計算結果與糧庫實測的溫度場隨時間變化情況一致，證明了測得的糧食熱學參數(shù)值與實際是相符的.

［1］ 張忠杰，李瓊，楊德勇,等．準靜態(tài)倉儲糧堆溫度場的CFD模擬［J］．中國糧油學報，2010，25（4）：46-50．

［2］ 李燦，高彥棟，黃素逸.熱傳導問題MATLAB數(shù)值計算［J］．華中科技大學學報，2002，30（9）：91-93．

［3］ 賈力，方肇洪，錢興華．高等傳熱學［M］．北京：高等教育出版社，2003．

［4］ 閆艷霞，曹玲芝．糧倉谷物溫度場數(shù)學模型的研究與應用［J］．糧食流通技術，2007，15（4）：40-42．

［5］ Iguaz A，Arroqui C，Esnoz A，et al．Modeling and validation of heat transfer in stored rough rice without aeration［J］．Biosystems Engineering，2004，88（4）：429-439．

［6］ 張德豐.MATLAB數(shù)值分析與應用［J］．北京：國防科技出版社，2007．

［7］ 賈燦純，曹崇文．倉儲小麥溫度場的數(shù)學模擬和實驗研究［J］．農(nóng)業(yè)機械學報，1998，29（1）：70-74.