圓筒倉儲糧通風(fēng)微環(huán)境模擬研究

魯子楓，王遠(yuǎn)成，尉堯方，曲安迪，杜傳致

（山東建筑大學(xué)熱能工程學(xué)院，山東濟南250101）

0 引言

糧食以及糧食的安全儲藏對于社會的穩(wěn)定甚至人類文明的發(fā)展有著非常重要的作用。倉儲糧堆內(nèi)部是一個復(fù)雜的多場耦合問題，其中包括溫度場、濕度場、水分場、生物場等。糧倉外部會受到大氣溫濕度和太陽輻射的作用，而糧食本身是一種有活性的物質(zhì)，種種因素相互耦合構(gòu)成復(fù)雜的生態(tài)系統(tǒng)，其中溫度和水分是糧食安全儲藏的重要因素。經(jīng)過研究發(fā)現(xiàn)，對于小麥儲藏問題，糧堆溫度＜15℃，且水分在12%～12.5%之間時有利于安全儲藏，可以有效避免蟲害的發(fā)生，抑制糧堆中生命體的活動，延緩糧食的劣變［1］。因此，將糧食的溫度和水分控制在合理的范圍之內(nèi)，可以保證糧食的安全儲藏。

基于計算流體力學(xué)的數(shù)值模擬方法是國內(nèi)外近年發(fā)展起來的一種研究流動、傳熱傳質(zhì)等現(xiàn)象的新方法，可以形象的再現(xiàn)流動、熱濕傳遞過程的情景，為解決儲糧通風(fēng)問題提供了一個良好的數(shù)值分析和優(yōu)化的工具［2］。借助于數(shù)值模擬方法國內(nèi)外研究人員取得了豐碩的成果。Cerconse等采用實驗和數(shù)值模擬相結(jié)合的方法，對圓筒倉內(nèi)的溫度和水分進行了研究［3］。曹崇文對各物的干燥過程做了詳細(xì)的綜述［4］。張忠杰等采用CFD模擬方法模擬了倉儲食堆內(nèi)溫度的變化［5］。Wang等和Gao等利用數(shù)值模擬方法對谷物的干燥和降溫過程進行了分析［6-7］。Li等基于壓力傳感器對水分的測量進行了探究［8］。姜俊伊等就糧倉的通風(fēng)方式對糧食的降溫效果進行了對比分析［9］。王利敏等對引起糧食霉菌的原因進行了探索，并提出了有效的防治措施［10］。但是他們只考慮了糧堆內(nèi)部單區(qū)域的溫濕水的變化，沒有考慮外部太陽輻射以及糧倉頂部空氣區(qū)域和糧食區(qū)域之間的耦合過程。王遠(yuǎn)成等對雙區(qū)域模型中溫度和水分進行了很好地分析和驗證［11-12］。但在計算流體力學(xué)模擬時，具有軟件對電腦硬件設(shè)施要求比較高，及耗費時間較長等問題。

文章基于倉儲糧堆內(nèi)部自然對流、熱濕耦合傳遞以及局部熱平衡原理，建立了熱濕傳遞和流動的數(shù)學(xué)模型，并引入干物質(zhì)損耗、谷蠹數(shù)量關(guān)聯(lián)式，基于Frotran語言編寫模擬計算程序，對通風(fēng)過程中圓筒倉內(nèi)部的小麥溫度、水分、干物質(zhì)損耗、谷蠹數(shù)量的變化進行了模擬，并且與實驗數(shù)據(jù)進行了比較，驗證了這個程序的準(zhǔn)確性和實用性。

1 儲糧通風(fēng)微環(huán)境模型及程序?qū)崿F(xiàn)

1.1 儲糧通風(fēng)倉原型

儲糧通風(fēng)倉型結(jié)構(gòu)及尺寸如圖1所示，圓筒倉為鋼板倉。筒倉高度為4.5 m，直徑為4.2 m。圖1、2分別為糧倉的正視圖和側(cè)視圖，圖3是糧倉的俯視圖，其中陰影部分是4個通風(fēng)口，分布在對稱的位置，其中每個通風(fēng)口長度為1.2 m。

圖1 倉型結(jié)構(gòu)及尺寸正視圖／mm

圖2 傳感器布置位置圖／mm

1.2 數(shù)學(xué)模型

基于多孔介質(zhì)傳熱傳質(zhì)理論和局部熱平衡的方法，考慮糧堆內(nèi)部的熱濕耦合和吸濕解吸濕的影響，建立儲糧機械通風(fēng)過程中糧堆內(nèi)部的流動和熱濕傳遞的數(shù)學(xué)模型［13］。

用圓柱坐標(biāo)（r，φ，z）表示動量方程、能量方程、水分方程、干物質(zhì)損耗、以及谷蠹和霉菌的數(shù)學(xué)模型。其中，r為半徑方向；φ為角度方向；z為高度方向。

圖3 倉型結(jié)構(gòu)及尺寸俯視圖 ／mm

1.2.1 動量方程

糧堆中壓強的分布符合拉普拉斯方程，由式（1）表示為

式中：P為糧堆的壓力，而當(dāng)時間趨近于無窮的時，近似認(rèn)為糧倉里壓力分布為穩(wěn)態(tài)，由式（2）表示為

由式（1）、（2）相等可得式（3）為

式中：αp為一個瞬變因子，Pa。

由達(dá)西定律可知，速度場和壓力場的分布關(guān)系用式（4）表示為

將式（2）、（3）帶入式（4）進而可得出3個方向上速度場的分布關(guān)系分別由式（5）～（7）表示為

式中：Rp為空氣有效阻力系數(shù)［14］，Pa·s／m。

1.2.2 能量守恒方程

根據(jù)局部熱平衡理論，熱平衡公式可以由式（8）表示為

1.2.3 水分守恒方程

根據(jù)質(zhì)量守恒原理，水分平衡公式由式（9）表示為

式中：Deff是由Thorpe提出來的表示有效擴散率［15］，m2／s；為干物質(zhì)損耗率；ρb、ρɑ分別為糧堆的密度和干空氣的密度，kg／m3。

1.2.4 干物質(zhì)損耗模型

Lacey等提出的干物質(zhì)損耗模型可由式（10）、（11）［16］表示為

式中：Qr為糧堆累積耗氧量，kg；T為糧堆中溫度分布，℃；t為儲藏時間，s；m為糧食的干物質(zhì)損失的重量，kg；M為糧食的濕基水分，%；b1、b2、b3、b4、b5、b6、b7、b8為糧食的特定經(jīng)驗常數(shù)，其值分別為3.4583×10-4、3.478×10-8、0.1737、20.33、0.9143、-2.878×10-7、-2.878×10-7、-0.013634、24.38。

1.2.5 谷蠹生長模型

Desmarchelier在報告中提出了谷蠹生長率的數(shù)學(xué)模型由（12）、（13）［17］表示為

式中：Rs為害蟲數(shù)量增長率；Tweb為濕球溫度，℃；谷蠹的Cc1，Cc2分別為經(jīng)驗常數(shù)0.0435、13.0；T為糧堆溫度，℃；RH為糧堆小數(shù)級的相對濕度。

1.3 程序?qū)崿F(xiàn)

模擬采用Codeblocks軟件中的Fortran語言編程。Fortran運行界面如圖4所示。模擬邊界條件設(shè)置采用實驗數(shù)據(jù)，模擬了10 d通風(fēng)實驗。每天固定在0：00 am～6：00 am進行通風(fēng)，程序中最初設(shè)定把糧倉的網(wǎng)格劃分為19個×19個×19個、39個×39個×39個、79個×79個×79個，通過模擬發(fā)現(xiàn)19個×19個×19個的模擬結(jié)果與另外2個方案相差無幾，為了節(jié)約模擬時間最終選用19個×19個×19個的網(wǎng)格劃分方案。程序初始需要調(diào)入一個狀態(tài)參數(shù)文件，文件中包括外部空氣的逐時溫度和相對濕度、倉型尺寸、糧堆的初始溫度、初始水分等，最后運行程序即可。

圖4 程序運行界面圖

選取的具體參數(shù)為：糧種為小麥，容重為800 kg／m3，初始平均糧食溫度為40℃，初始平均濕基水分為 12%，噸糧通風(fēng)量為 13.8 L／（s·t），通風(fēng)階段日最大平均干球溫度為28.7℃，平均日最小干球溫度為13℃，9：00平均相對濕度為47%，而15：00平均相對濕度為28%，每天通風(fēng)6 h，默認(rèn)通風(fēng)參數(shù)和大氣環(huán)境中的參數(shù)一致。

2 結(jié)果與分析

2.1 儲糧通風(fēng)實驗

選取入倉之后的小麥進行通風(fēng)，通風(fēng)時間從2016年11月9日至19日。糧倉中一共布置了11個溫度和濕度傳感器，其中3個在糧倉中心線上，分別布置在距最下端1.5、3、4.5 m處，另外8個測點分別分布在東南（SE）、東北（NE）、西南（SW）、西北（NW）4個方向上，每個方向有2個測點，測點名稱為 SEA、SEB、NEA、NEB、SWA、SWB、NWA、NWB，測點距倉壁52 cm。

實驗采集系統(tǒng)是由傳感器（vaisala HMP45A vo520014）傳來的信號到達(dá)數(shù)據(jù)記錄儀（HP E1326B 5 1／2＠14083）轉(zhuǎn)換之后傳遞給PC記錄儀記錄數(shù)據(jù)。數(shù)據(jù)記錄間隔為31 min。將實驗獲得數(shù)據(jù)和模擬數(shù)據(jù)進行對比分析。

2.2 糧堆水分模擬結(jié)果對比分析

程序一共模擬了通風(fēng)10 d的結(jié)果，其中輸出了3、6和9 d的溫度和水分?jǐn)?shù)據(jù)水分的模擬結(jié)果如圖5所示。水分分布圖近似對稱分布，由于溫度場和微氣流的影響靠近通風(fēng)道和倉壁的水分偏大，隨著通風(fēng)天數(shù)的增加，水分略有降低，但對于降低的幅度不如溫度降低的大。糧倉模擬平均濕基水分?jǐn)?shù)據(jù)與實測數(shù)據(jù)的對比見表1。

根據(jù)表1分析可知，隨著通風(fēng)時間的延長，糧堆水分略有增加，這與大氣中的濕度及糧堆中的溫度場有直接性的關(guān)系。水分的模擬數(shù)據(jù)和實測數(shù)據(jù)誤差最大為5.9%，驗證了文章中的數(shù)學(xué)模型和編寫的Fortran程序的準(zhǔn)確性和實用性，對以后糧食的安 全儲藏有著重要指導(dǎo)意義。

圖5 通風(fēng)過程中水分的模擬結(jié)果圖

表1 糧倉模擬平均濕基水分?jǐn)?shù)據(jù)與實測數(shù)據(jù)對比表

2.3 糧堆溫度實驗與模擬結(jié)果對比分析

通風(fēng)過程中溫度模擬結(jié)果如圖6所示。由于風(fēng)道處于糧倉的對稱位置，可以明顯看出溫度場的分布近似為對稱分布，而且由于風(fēng)道的位置處于兩側(cè)，糧倉兩側(cè)的溫度比中心的溫度低?？拷Z倉四周部分由于受到太陽輻射的影響，溫度比糧倉中心的溫度高。圖6（a）為通風(fēng)3 d后的溫度分布，糧堆上半部分與下半部分的溫差大概為7℃。圖6（b）為通風(fēng)6 d后的溫度分布，上、中、下層溫差近似為4℃。圖6（c）為通風(fēng)9 d后的溫度分布，溫度場比較均勻。分析可知糧堆由于通風(fēng)的原因，糧堆的由下而上溫度逐漸降低，最后溫度場變得均勻。

同時，用模擬出的逐時溫度數(shù)據(jù)與實驗實測的溫度數(shù)據(jù)進行對比，結(jié)果如圖7～8所示。

通過比較圖7糧倉中心部分3個測點和圖8中靠近倉壁的8個測點的溫度數(shù)據(jù)可知：隨著通風(fēng)時間的延長，糧堆各個位置的溫度明顯降低，且靠近倉壁的溫度降低幅度比糧倉中心降低的幅度大。這是因為通風(fēng)道的位置和倉型的影響；在通風(fēng)100和180 h時，溫度上升較明顯，由于受外界大氣溫濕度的影響較大；模擬數(shù)據(jù)和實驗數(shù)據(jù)較為吻合，溫度誤差最大為2℃

圖6 通風(fēng)過程中溫度的模擬結(jié)果圖

圖7 糧倉中心各溫度測點的實驗數(shù)據(jù)與模擬數(shù)據(jù)對比圖

圖8 靠近壁面各個溫度測點實驗數(shù)據(jù)與模擬數(shù)據(jù)的對比圖

2.4 干物質(zhì)損耗、谷蠹數(shù)量變化的模擬結(jié)果對比分析

干物質(zhì)損耗、谷蠹數(shù)量變化的模擬結(jié)果如圖9～10所示。

通過分析圖9、10可知，由于通風(fēng)樓處于對稱位置，干物質(zhì)損耗及谷蠹的分布近似對稱分布，干物質(zhì)損耗一般集中在靠近倉壁和通風(fēng)樓的區(qū)域。通風(fēng)3 d時，通風(fēng)樓附近的干物質(zhì)損耗為2.5×10-5kg；通風(fēng)6 d時，通風(fēng)樓附近的干物質(zhì)損耗為7.9×10-4kg；通風(fēng)9 d時，通風(fēng)樓附近的干物質(zhì)損耗為3.5×10-3kg?？傮w來看干物質(zhì)損耗不是很大，但是隨著通風(fēng)時間的延長，干物質(zhì)損耗增加。而從谷蠹的分布圖可知，通風(fēng)3 d時，谷蠹為1.108個／kg；通風(fēng)6 d的時候，谷蠹為1.416個／kg；通風(fēng)9 d時，谷蠹為1.533個／kg。綜合糧倉中溫度和水分的分布可知，隨著通風(fēng)時間的延長，溫度和水分變化梯度大的區(qū)域，干物質(zhì)損耗和谷蠹數(shù)量較大。可見，干物質(zhì)損耗和谷蠹數(shù)量的變化和溫度、水分和通風(fēng)的時間有關(guān)。

圖9 干物質(zhì)損耗隨通風(fēng)時間的模擬分布圖

圖10 谷蠹隨通風(fēng)時間延長的模擬分布圖

3 結(jié)論

通過上述研究可知：

（1）溫度和水分的模擬數(shù)據(jù)和實驗數(shù)據(jù)較為吻合。其中，其最大誤差分別為2℃和5.9%，但是誤差在糧食安全儲存可接受范圍之內(nèi)。同時，通過與實驗真實測得的數(shù)據(jù)對比結(jié)果驗證了由Fortran程序設(shè)計的數(shù)學(xué)模型的實用性。

（2）相對其他 CFD模擬軟件而言，設(shè)計的Fortran程序模擬時間極短（47 s）。這套程序不限糧種，只需調(diào)整物性參數(shù)。

（3）溫度和水分近似呈對稱分布，靠近壁面和通風(fēng)樓的區(qū)域，溫度和水分變化梯度較大，且由于太陽輻射的影響，靠近壁面的溫度比糧倉中心的溫度

高。

（4）干物質(zhì)損耗和谷蠹在糧倉中分布近似呈對稱分布，并且與溫度、水分以及通風(fēng)時間有關(guān)。

參考文獻(xiàn)：

［1］王遠(yuǎn)成，段海峰，張來林.就倉通風(fēng)時糧堆內(nèi)部熱濕耦合傳遞過程的數(shù)值模擬［J］.河南工業(yè)大學(xué)學(xué)報，2009，30（6）：76-77.

［2］王遠(yuǎn)成，張忠杰，吳子丹，等.計算流體力學(xué)技術(shù)在糧食儲藏中的應(yīng)用［J］.中國糧油學(xué)報，2012，27（5）：86-91.

［3］Converse H H，Graces A H，Chung D S.Transient heat transfer within wheat stored in a cylindrical bin［J］.Transaction of American Society Agriculture Engineering，1973，16（1）：129-133.

［4］曹崇文.谷物干燥的數(shù)學(xué)模擬［J］.北京農(nóng)業(yè)機械化學(xué)院學(xué)報，1984，3：79-94.

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

［6］Wang R L，Liu S Q，Huan L I，etɑl..Study of temperature changes in wheatgrain heap during cooling process［J］.Cereals＆Oils，2017，30（5）：23-27.

［7］Gao S，Wang Y，Qiu H，etɑl..Numerical simulation for anisotropic resistance of wheat under the condition ofmechanical ventilation［J］.Journal of the Chinese Cereals＆Oils Association，2017，32（2）：116-119.

［8］Li X J，Zhang Y Q，Wu FW，et.ɑl.Research of grain moisture measurement method based on pressure sensor［J］.Advanced Materials Research，2012（472）：1518-1523.

［9］姜俊伊，李倩倩，沈邦灶等.平房倉橫向與豎向通風(fēng)系統(tǒng)儲糧溫度變化對比研究［J］.糧食儲藏技術(shù)，2017，46（6）：1-7.

［10］王利敏，邢福國，呂聰，等.復(fù)合植物精油防霉劑對玉米霉菌及真菌霉素的控制效果［J］.核農(nóng)學(xué)報，2018，32（4）：0732-0739.

［11］王遠(yuǎn)成，潘鈺，尉堯方，等.倉儲糧堆內(nèi)部自然對流和熱濕傳遞的數(shù)學(xué)分析和驗證［J］.中國糧油學(xué)報，2017，32（9）：120-130.

［12］尉堯方，王遠(yuǎn)成，王興周，等.儲糧通風(fēng)模型的構(gòu)建及其應(yīng)用分析［J］.山東建筑大學(xué)學(xué)報，2017，32（3）：251-257.

［13］Thorpe GR.The application of computational fluid dynamics codes to simulate heat and moisture transfer in stored grains［J］.Journal of Stored Products Research，2008，44（1）：21-31.

［14］Hunter A J.An isostere equation for some common seeds［J］.Journal of Agricultural Engineering Research，1987，37（2）：93-105.

［15］Thorpe G R，Ochoa JA，Whitaker S.The diffusion ofmoisture in food grains［J］.Journal of Stored Products Research，1991，27（3）：1-9.

［16］Lacey J，Hamer A，Magan N.Respiration and loss in stored wheat grain under different environmental condition［C］.Proceedings of the 6th international Working Conference on Stored-product Protection，1994，2：1007-1013.

［17］Desmarchelier M J.The relationship between wet-bulb temperature and the intrinsic rate of increase of eight species of stored-product coleoptera［J］.Journal of Stored Products Research，1998，24（3）：107-113.