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

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

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

0 引言

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

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

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

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

1.1 儲(chǔ)糧通風(fēng)倉(cāng)原型

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

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

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

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

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

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

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

1.2.1 動(dòng)量方程

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

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

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

式中：αp為一個(gè)瞬變因子，Pa。

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

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

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

1.2.2 能量守恒方程

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

1.2.3 水分守恒方程

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

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

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

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

式中：Qr為糧堆累積耗氧量，kg；T為糧堆中溫度分布，℃；t為儲(chǔ)藏時(shí)間，s；m為糧食的干物質(zhì)損失的重量，kg；M為糧食的濕基水分，%；b1、b2、b3、b4、b5、b6、b7、b8為糧食的特定經(jīng)驗(yàn)常數(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 谷蠹生長(zhǎng)模型

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

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

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

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

圖4 程序運(yùn)行界面圖

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

2 結(jié)果與分析

2.1 儲(chǔ)糧通風(fēng)實(shí)驗(yàn)

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

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

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

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

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

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

表1 糧倉(cāng)模擬平均濕基水分?jǐn)?shù)據(jù)與實(shí)測(cè)數(shù)據(jù)對(duì)比表

2.3 糧堆溫度實(shí)驗(yàn)與模擬結(jié)果對(duì)比分析

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

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

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

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

圖7 糧倉(cāng)中心各溫度測(cè)點(diǎn)的實(shí)驗(yàn)數(shù)據(jù)與模擬數(shù)據(jù)對(duì)比圖

圖8 靠近壁面各個(gè)溫度測(cè)點(diǎn)實(shí)驗(yàn)數(shù)據(jù)與模擬數(shù)據(jù)的對(duì)比圖

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

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

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

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

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

3 結(jié)論

通過(guò)上述研究可知：

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

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

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

高。

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

參考文獻(xiàn)：

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

［2］王遠(yuǎn)成，張忠杰，吳子丹，等.計(jì)算流體力學(xué)技術(shù)在糧食儲(chǔ)藏中的應(yīng)用［J］.中國(guó)糧油學(xué)報(bào)，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è)機(jī)械化學(xué)院學(xué)報(bào)，1984，3：79-94.

［5］張忠杰，李瓊，楊德勇，等.準(zhǔn)靜態(tài)倉(cāng)儲(chǔ)糧堆溫度場(chǎng)的CFD模擬［J］.中國(guó)糧油學(xué)報(bào)，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］姜俊伊，李倩倩，沈邦灶等.平房倉(cāng)橫向與豎向通風(fēng)系統(tǒng)儲(chǔ)糧溫度變化對(duì)比研究［J］.糧食儲(chǔ)藏技術(shù)，2017，46（6）：1-7.

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

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

［12］尉堯方，王遠(yuǎn)成，王興周，等.儲(chǔ)糧通風(fēng)模型的構(gòu)建及其應(yīng)用分析［J］.山東建筑大學(xué)學(xué)報(bào)，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.