倉(cāng)儲(chǔ)稻谷黃變的數(shù)值模擬研究

余海，王遠(yuǎn)成，*，陳朝暉，尹君，楊開(kāi)敏

(1.山東建筑大學(xué)熱能工程學(xué)院，山東濟(jì)南 250101；2.山東省公共資源交易中心，山東濟(jì)南 250014；3.國(guó)家糧食和物資儲(chǔ)備局科學(xué)研究院，北京 100037)

0 引言

稻谷的顏色是衡量其品質(zhì)和價(jià)格的重要標(biāo)準(zhǔn)之一，但稻谷儲(chǔ)藏一段時(shí)間后會(huì)產(chǎn)生黃變，其影響因素主要為溫度、濕度和儲(chǔ)藏時(shí)間。 新收獲的稻谷水分含量＞20%，該水分含量下稻谷易生長(zhǎng)霉菌，因此對(duì)新收獲的稻谷采用干燥處理，將水分含量降至＜15%，方可保證稻谷安全儲(chǔ)藏。 稻谷儲(chǔ)藏可分為自然儲(chǔ)藏階段(即非人工干預(yù)階段)和人工干預(yù)階段，自然儲(chǔ)藏時(shí)圓筒倉(cāng)是封閉的，與外界不可滲透，又稱為密閉儲(chǔ)藏。

GRAS 等[1]和 BASON 等[2]分別通過(guò)實(shí)驗(yàn)研究了溫度、相對(duì)濕度和空氣成分對(duì)兩個(gè)品種白米和稻谷顏色的影響，發(fā)現(xiàn)溫度和相對(duì)濕度對(duì)白米和稻谷黃變速率有影響，氧氣對(duì)黃變速率有輕微影響，而二氧化碳對(duì)黃變速率無(wú)影響。 PHILLIPS 等[3]研究了稻谷在干燥前的黃變現(xiàn)象，指出在一定時(shí)間內(nèi)，散裝高水分的稻谷比袋裝的低水分稻谷黃變速度更快。戚禹康等[4]、俞曉靜等[5]、王小萌等[6-7]和魯子楓等[8]采用計(jì)算流體動(dòng)力學(xué)(Computational Fluid Dynamics，CFD)的方法模擬了糧堆溫度和水分的變化，發(fā)現(xiàn)兩種參數(shù)互有影響。 歐陽(yáng)毅等[9]研究了糧堆溫度和水分的變化與真菌生長(zhǎng)的關(guān)系，得到了糧堆儲(chǔ)藏溫度和水分與真菌生長(zhǎng)的預(yù)測(cè)關(guān)系曲線。 劉進(jìn)吉等[10]探索了高大淺圓倉(cāng)儲(chǔ)糧溫度的變化，發(fā)現(xiàn)在壁面附近溫度梯度較大。 唐為民[11]分析了黃粒米的成因并提出了一些預(yù)防方法。 萬(wàn)立昊等[12]在探究稻谷儲(chǔ)藏質(zhì)量影響因素的過(guò)程中發(fā)現(xiàn)，黃粒米隨溫度和水分含量升高而增多。

雖然學(xué)者們對(duì)稻谷黃變進(jìn)行過(guò)大量的實(shí)驗(yàn)研究[1-3，12]，結(jié)果較為準(zhǔn)確且合乎實(shí)際，但是稻谷的黃變周期很長(zhǎng)，且在高溫高濕的情況下才容易黃變，隨著時(shí)間的積累，黃變將越來(lái)越嚴(yán)重。 國(guó)內(nèi)對(duì)稻谷黃變的研究與國(guó)外相比相對(duì)較少，對(duì)安全儲(chǔ)糧的研究主要停留于溫度和水分的變化規(guī)律，少數(shù)研究稻谷黃變還只是通過(guò)實(shí)驗(yàn)進(jìn)行的，但是通過(guò)實(shí)驗(yàn)研究將耗費(fèi)大量的時(shí)間和人力物力。 近年來(lái)，隨著電腦技術(shù)的高速發(fā)展，涌現(xiàn)出一種用電腦模擬流體流動(dòng)、傳熱傳質(zhì)的新技術(shù)，即CFD 的方法，為糧食儲(chǔ)藏問(wèn)題提供了一個(gè)良好的工具[13-14]。對(duì)比大量的研究現(xiàn)狀發(fā)現(xiàn)，數(shù)值模擬時(shí)主要為設(shè)定常壁溫邊界條件，然而大氣溫度是變化的，因此采用此條件模擬的結(jié)果與真實(shí)糧情有偏差。

文章建立了倉(cāng)儲(chǔ)稻谷自然儲(chǔ)藏過(guò)程中糧堆內(nèi)部流動(dòng)和熱濕耦合傳遞的數(shù)學(xué)模型，并結(jié)合已有的稻谷黃變動(dòng)力學(xué)模型，通過(guò)多物理場(chǎng)仿真軟件COMSOL Multiphysics 數(shù)值模擬了稻谷黃變過(guò)程，溫度邊界條件設(shè)置為動(dòng)態(tài)變化的大氣環(huán)境溫度，模擬環(huán)境更加接近真實(shí)糧情，研究分析了局地氣候條件下倉(cāng)內(nèi)糧堆局部溫度、水分的分布規(guī)律以及稻谷黃變規(guī)律，以期為安全儲(chǔ)糧提供精確的理論指導(dǎo)。

1 稻谷黃變模型的建立與條件設(shè)置

1.1 稻谷黃變物理模型

以中試圓筒倉(cāng)為研究對(duì)象，由于倉(cāng)壁厚度不大，通過(guò)倉(cāng)壁的傳熱影響較小，為了簡(jiǎn)化模型未考慮圓筒倉(cāng)壁厚和其對(duì)糧食的熱傳導(dǎo)。 圓筒倉(cāng)的高度為6 m、筒倉(cāng)半徑為1.5 m。 糧食在儲(chǔ)藏過(guò)程中，內(nèi)部糧情參數(shù)是動(dòng)態(tài)變化的，因此在糧堆內(nèi)部設(shè)定監(jiān)測(cè)點(diǎn)來(lái)記錄參數(shù)變化情況，其中監(jiān)測(cè)點(diǎn)1 位于糧堆左下角(-1.25 m，0.25 m)處，監(jiān)測(cè)點(diǎn)2 位于糧堆中心(0 m，3 m)處，而監(jiān)測(cè)點(diǎn)3 位于糧堆右上角(1.25 m，5.75 m)處，位置示意圖如圖1 所示。 由于圓筒倉(cāng)為軸對(duì)稱圖形，沿直徑各截面都一樣，因此選取沿筒倉(cāng)直徑的截面作為研究對(duì)象。 數(shù)值計(jì)算區(qū)域的網(wǎng)格劃分如圖2 所示，并且對(duì)糧倉(cāng)壁面進(jìn)行網(wǎng)格加密處理。

圖1 監(jiān)測(cè)點(diǎn)位置示意圖

圖2 計(jì)算區(qū)域網(wǎng)格圖

1.2 稻谷黃變數(shù)學(xué)模型

1.2.1 連續(xù)性方程

密閉儲(chǔ)藏過(guò)程中，由于淺圓倉(cāng)是封閉的，與外界不可滲透，并假定空氣不可壓縮，根據(jù)質(zhì)量守恒，連續(xù)性方程[15]由式(1)表示為

式中 ρa(bǔ)為空氣密度，kg/m3；ε 為糧堆孔隙率；t 為時(shí)間，s；?為哈密頓算子；u 為空氣的表觀流速，m/s。

1.2.2 動(dòng)量方程

考慮浮力效應(yīng)，空氣密度隨溫度的變化近似于布森內(nèi)斯克方程(Boussinesq equation)[16]，根據(jù)動(dòng)量守恒，其方程由式(2)表示為

式中P 為多孔介質(zhì)中的壓力，Pa；μ 為流體的動(dòng)力黏度，Pa·s；g 為重力加速度，m/s2；β 為空氣的體積膨脹系數(shù)，β ＝ 1/T0，K-1；T0和 T 分別為糧堆初始溫度和實(shí)時(shí)溫度，K；d 為稻谷顆粒當(dāng)量直徑，mm；K 為糧堆的滲透率；ui(i ＝ 1，2) 為空氣的流動(dòng)速度，m/s；δij為糧粒空隙間距，m；xi為坐標(biāo)方向，x1＝x 時(shí)，為水平方向，x2＝ y 時(shí)，為垂直方向。

1.2.3 能量方程

密閉儲(chǔ)藏過(guò)程中，由于淺圓倉(cāng)是封閉的，與外界不可滲透，根據(jù)能量守恒定律，能量方程由式(3)表示為

式中 ρs為糧堆的密度，kg/m3；Cs和 Ca分別為稻谷和空氣的比熱容，J/(kg·K)；keff為糧堆的有效導(dǎo)熱系數(shù)，W/(m·K)；Wg為稻谷糧堆的干基水分，為吸濕或解吸濕熱，J。

1.2.4 水分守恒方程

自然儲(chǔ)藏時(shí)糧堆內(nèi)的水分遷移滿足質(zhì)量守恒，由式(4)表示為

式中ω 為稻谷顆粒間空氣的絕對(duì)含濕量，kg/kg；Deff吸濕產(chǎn)生的干基水分質(zhì)量，kg。

1.2.5 稻谷黃變模型

稻谷在倉(cāng)內(nèi)儲(chǔ)藏時(shí)，在溫、濕度的共同影響下，稻谷將慢慢變黃，黃度由b 值計(jì)量[17]，由式(5)表示為

其中，

式中φ 為相對(duì)濕度，%。

1.3 數(shù)值模擬參數(shù)設(shè)置

模型采用有限元法研究稻谷，糧堆的孔隙率ε 為0.48、密度ρs為580 kg/m3、滲透率K ＝ 7.27 × 10-9m3、比熱容Cs為(1269 ＋ 34.89M) W/m·k-1，其中M 為糧堆的濕基水分，M ＝Wg/ (1 ＋Wg) × 100。 文章對(duì)不同溫度和相對(duì)濕度的9 種工況進(jìn)行模擬研究，具體參數(shù)見(jiàn)表1，其中恒溫恒濕工況1、2、3 的模擬為對(duì)文獻(xiàn)[17] 實(shí)驗(yàn)結(jié)果的驗(yàn)證(該工況數(shù)據(jù)源自SOPONRONNARIT 的實(shí)驗(yàn)數(shù)據(jù))，其他工況為秋冬季(2016 年 10 月～2017 年 2 月)與春夏季(2017 年4 月～2017 年 8 月)各 5 個(gè)月的工況模擬，秋冬季與春夏季溫度為大氣溫度(簡(jiǎn)稱氣溫)(數(shù)據(jù)取自浙江省某糧庫(kù)記錄的溫度)如圖3 所示。

圖3 秋冬季與春夏季溫度圖

表1 9 種模擬工況參數(shù)表

2 模擬結(jié)果與分析

2.1 工況 1、2 和 3 模擬結(jié)果分析

SOPONRONNARIT 等[17]采用實(shí)驗(yàn)的方法測(cè)量了稻谷黃變的黃度，步驟為(1) 將飽和鹽溶液裝在一個(gè)玻璃材質(zhì)的器皿中，依次在器皿上方鋪上一層塑料網(wǎng)和初始濕基水分為19.35%的稻谷；(2) 將該玻璃器皿放置在塑料箱內(nèi)，將塑料箱進(jìn)行密封后放于烤箱中。 實(shí)驗(yàn)時(shí)，將烤箱溫度、飽和鹽溶液相對(duì)濕度分別調(diào)為 35 ℃與 80%、60 ℃與 80%、60 ℃與95.2%，實(shí)驗(yàn)總時(shí)長(zhǎng)為6 d，用型號(hào)為JP7100 的色差儀測(cè)量稻谷的黃度，初始黃度為11.50。

文章基于建立的數(shù)學(xué)模型，采用 COMSOL Multiphy 軟件模擬上述3 種實(shí)驗(yàn)工況，將模擬結(jié)果與實(shí)驗(yàn)測(cè)量的結(jié)果對(duì)比如圖4 所示，兩者的變化趨勢(shì)基本相同，最大誤差＜3%，說(shuō)明通過(guò)數(shù)值模擬的方法研究稻谷黃變是可行的。

圖4 模擬結(jié)果與實(shí)驗(yàn)結(jié)果的比較圖

3 種工況下稻谷黃度的分布結(jié)果如圖5 所示?？梢钥闯觯銣睾銤駰l件下，倉(cāng)內(nèi)稻谷黃變分布大致均勻；在相對(duì)濕度恒定的情況下，溫度越高，黃度值越大；在恒溫情況下，相對(duì)濕度越大，稻谷黃度值也越大，結(jié)果與圖4 一致，說(shuō)明稻谷在高溫高濕的環(huán)境中黃變?cè)絿?yán)重，稻谷的黃變與時(shí)間呈線性關(guān)系，且隨著時(shí)間的增加黃變?cè)絹?lái)越嚴(yán)重。

圖 5 6 d 時(shí)工況 1、2、3 的黃度分布圖

2.2 秋冬季與春夏季各工況模擬結(jié)果分析

2.2.1 自然對(duì)流速度場(chǎng)分布

秋冬季和春夏季工況下90 d 時(shí)的速度場(chǎng)分布結(jié)果如圖6 所示。 自然對(duì)流流動(dòng)速度大小的數(shù)量級(jí)大概為10-5～10-4m/s，秋冬季工況下，在對(duì)稱軸(x ＝0)右側(cè)，微氣流為順時(shí)針流動(dòng)，對(duì)稱軸左側(cè)，微氣流為逆時(shí)針流動(dòng)，且壁面附近與糧堆內(nèi)部流動(dòng)相對(duì)更強(qiáng)；春夏季工況下，氣體流動(dòng)方向與秋冬季工況相反。

圖6 秋冬季與春夏季工況90 d 時(shí)速度場(chǎng)圖

2.2.2 糧堆溫度分布

秋冬季1 與春夏季1 工況下溫度場(chǎng)及溫度變化圖如圖7 所示。

圖7 秋冬季1 與春夏季1 工況下溫度場(chǎng)及溫度變化圖

由圖7(a)和(b)可以看出，整個(gè)儲(chǔ)藏周期內(nèi)，秋冬季1 工況下糧堆平均溫度由20 ℃下降至14.5 ℃，其降幅為5.5 ℃；春夏季1 工況下糧堆平均溫度由 0 ℃上升至 28.2 ℃，其升幅為 28.2 ℃。150 d時(shí)，秋冬季工況糧堆處于“熱芯糧”狀態(tài)，春夏季工況時(shí)糧堆為“冷芯糧”狀態(tài)，這是由于秋冬季時(shí)外界環(huán)境溫度較低，春夏季時(shí)其溫度較高，而倉(cāng)壁附近的糧食與倉(cāng)頂糧面受外界環(huán)境的影響較大，因此秋冬季時(shí)倉(cāng)壁與倉(cāng)頂附近的糧溫較低，糧堆內(nèi)部溫度相對(duì)較高；春夏季時(shí)倉(cāng)壁與倉(cāng)頂附近的糧溫較高，糧堆內(nèi)部溫度相對(duì)較低。 由于外界溫度呈波動(dòng)狀，溫度的傳遞有一定的延滯性，且糧食本身導(dǎo)熱系數(shù)比較小，所以外界溫度對(duì)糧堆影響區(qū)域較小。 由圖7(c)和 (d)可以看出，秋冬季1 工況下3 個(gè)監(jiān)測(cè)點(diǎn)溫度均呈下降趨勢(shì)，春夏季1 則呈上升趨勢(shì)，由于監(jiān)測(cè)點(diǎn)1 和3 靠近倉(cāng)壁，而倉(cāng)壁附近受外界環(huán)境影響很大，因此監(jiān)測(cè)點(diǎn)1 和3 的溫度變化趨勢(shì)與大氣溫度的變化趨勢(shì)大致一致。 同時(shí)還可看出，無(wú)論是秋冬季還是春夏季，監(jiān)測(cè)點(diǎn)3 處的溫度變化最大，在春夏季時(shí)，監(jiān)測(cè)點(diǎn)3 處溫度最高。

2.2.3 糧堆濕基水分分布

秋冬季1 與春夏季 1 糧堆初始濕基水分為14%，150 d 濕基水分場(chǎng)及水分變化圖如圖8 所示。由圖8(a)和(b)可知，秋冬季時(shí)，倉(cāng)頂糧面及倉(cāng)壁附近濕基水分較高，糧堆內(nèi)部和倉(cāng)底濕基水分較低；在春夏季，倉(cāng)頂糧面及倉(cāng)壁附近、糧堆底角處濕基水分較低，糧堆內(nèi)上半部水濕基水分較高；由圖8(c)和 (d)可知，秋冬季1 工況下，監(jiān)測(cè)點(diǎn)1 和2 處濕基水分均呈降低趨勢(shì)，其中監(jiān)測(cè)點(diǎn)1 處濕基水分下降最快，監(jiān)測(cè)點(diǎn)3 處濕基水分略微升高；春夏季1 工況時(shí)，監(jiān)測(cè)點(diǎn)1 處濕基水分降低，監(jiān)測(cè)點(diǎn)2 和3 處濕基水分略微升高，這是因?yàn)樵谧匀粚?duì)流的作用下，秋冬季時(shí)濕基水分從糧堆內(nèi)部和底部向倉(cāng)壁和倉(cāng)頂糧面遷移，其中倉(cāng)底附近濕基水分最低；春夏季時(shí)濕基水分從倉(cāng)壁和倉(cāng)頂糧面向糧堆上半部遷移，其中糧堆左下角與右下角濕基水分最低。 同時(shí)還可知，無(wú)論是秋冬季還是春夏季，監(jiān)測(cè)點(diǎn)3 處都是高水分區(qū)。

圖8 秋冬季1 與春夏季1 工況下濕基水分場(chǎng)及水分變化圖

2.2.4 糧堆內(nèi)部黃度變化

不同初始濕基水分時(shí)，各工況不同時(shí)間點(diǎn)下計(jì)算得到的糧堆平均黃度值，見(jiàn)表2。

表2 秋冬季與春夏季各工況不同時(shí)間點(diǎn)的模擬黃度表

由表2 可知，當(dāng)?shù)竟瘸跏紳窕譃?4%時(shí)，秋冬季1 工況和春夏季1 工況的稻谷平均黃度值均變化較小，但春夏季和秋冬季相比，春夏季環(huán)境下稻谷黃變略微嚴(yán)重。 當(dāng)?shù)竟瘸跏紳窕譃?6%時(shí)，在秋冬季2 工況下，稻谷的平均黃度值由11.50 上升至11.54，在春夏季 2 工況下則由 11.50 上升至11.62。 當(dāng)?shù)竟瘸跏紳窕譃?9%時(shí)，秋冬季3 工況下的稻谷平均黃度值由11.50 上升至11.64，和秋冬季1、2 工況相比，稻谷黃變程度相對(duì)更加嚴(yán)重；春夏季3 工況下，稻谷的平均黃度值由11.50 上升至11.87。 表明初始濕基水分和氣溫越高，則稻谷黃變程度越嚴(yán)重。 另外，從表2 還可知，相同初始濕基水分時(shí)，前期時(shí)間段內(nèi)秋冬季工況下比春夏季工況下的稻谷平均黃度值大，這是因?yàn)榍捌跁r(shí)間段內(nèi)，秋冬季工況下稻谷的平均溫度高于春夏季，但到了后期，春夏季工況下稻谷的平均溫度高于秋冬季，因此到150 d 時(shí)，春夏季工況下的稻谷平均黃度值更大。無(wú)論是秋冬季工況還是春夏季工況，稻谷的初始濕基水分越高，稻谷的黃度速度越快，與初始濕基水分為14%和16%相比，當(dāng)?shù)竟鹊某跏紳窕譃?9%時(shí)，稻谷的黃變速度明顯更快。

秋冬季與春夏季工況下3 個(gè)監(jiān)測(cè)點(diǎn)處的黃度變化如圖9～11 所示，由圖9～11 的(a)可知，秋冬季工況下，3 個(gè)監(jiān)測(cè)點(diǎn)的黃度值最終相差較小。 由圖9 ～11 的(b)可知，春夏季工況下，監(jiān)測(cè)點(diǎn)1 和2 的黃度值變化也較小，監(jiān)測(cè)點(diǎn)3 處的黃度變化較大，在春夏季1 工況(初始濕基水分為14%)，監(jiān)測(cè)點(diǎn)3 處的黃度由11.50 上升至11.67，升幅為0.17，在春夏季2工況(初始濕基水分為16%)，監(jiān)測(cè)點(diǎn)3 處的黃度由11.50 上升至 11.95，升幅為 0.45，在春夏季 3 工況(初始濕基水分為19%)，監(jiān)測(cè)點(diǎn)3 處的黃度由11.50上升至12.8，升幅為1.3，說(shuō)明在春夏季工況下，監(jiān)測(cè)點(diǎn)3 處即靠近倉(cāng)頂和倉(cāng)壁交界附近區(qū)域黃度受影響較大，結(jié)合圖7(c)和 (d)、圖8(c)和 (d)也可知，監(jiān)測(cè)點(diǎn)3 處溫度和水分變化比較大，而黃變與溫濕度密切相關(guān)，在高溫度和高水分的雙重影響下，該區(qū)域的稻谷黃變更加嚴(yán)重，因此，靠近倉(cāng)頂和倉(cāng)壁交界附近區(qū)域?yàn)槲ｋU(xiǎn)區(qū)域。

圖9 初始濕基水分為14%時(shí)各監(jiān)測(cè)點(diǎn)黃度變化情況圖

圖10 初始濕基水分為16%時(shí)各監(jiān)測(cè)點(diǎn)黃度變化情況圖

圖11 初始濕基水分為19%時(shí)各監(jiān)測(cè)點(diǎn)黃度變化情況圖

3 結(jié)論

文章采用有限元法數(shù)值模擬了圓筒倉(cāng)內(nèi)稻谷自然儲(chǔ)藏過(guò)程中糧堆溫度、水分和黃變的變化規(guī)律，分析了3 種恒溫恒濕工況(溫度35 ℃和相對(duì)濕度80%、溫度60 ℃和相對(duì)濕度80%、溫度60 ℃和相對(duì)濕度95.2%)，以及春夏季與秋冬季大氣溫度環(huán)境6種工況，得出以下結(jié)論:

(1) 稻谷的黃變情況受儲(chǔ)藏環(huán)境的影響較大，恒溫恒濕工況下，稻谷的黃變與時(shí)間呈線性關(guān)系，隨著時(shí)間的增加黃變?cè)絹?lái)越嚴(yán)重；低溫環(huán)境下稻谷黃變緩慢，高溫環(huán)境中稻谷黃變較快，高溫高濕環(huán)境下稻谷極易黃變。

(2) 在秋冬季，由于糧溫較低，稻谷黃變緩慢；在春夏季，隨著糧溫的升高，稻谷黃變較快；在相同的儲(chǔ)藏環(huán)境中，和稻谷的初始水分為14%和16%相比，當(dāng)?shù)竟鹊某跏妓譃?9%時(shí)，稻谷的黃變速度明顯更快。

(3) 稻谷儲(chǔ)藏過(guò)程中，靠近倉(cāng)頂和倉(cāng)壁交界附近區(qū)域?yàn)槲ｋU(xiǎn)區(qū)域，該區(qū)域溫度和水分變化較大且黃變嚴(yán)重。