淺圓倉和立筒倉糧堆氣流分布的影響因素及模擬研究進展

李興軍 張洪清

摘要：淺圓倉和立筒倉占我國現代化糧倉總量的12%，節約用地。從倉頂進糧時由于糧食靜止角的作用總是倉中心富集雜質導致孔隙度低于倉壁附近，倉中心糧層阻力往往大于倉壁附近，國內對通風期間這兩種倉型的糧堆氣流分布和壓強降缺乏數學模擬研究。高大的圓柱形糧堆上層糧食對底層糧食產生靜態壓強，可導致糧粒變形、破壞，破碎率和生化活性成分隨著局部糧堆含水率和溫度增加而顯著變化。從發展角度看需要研究淺圓倉和立筒倉糧堆氣流分布的影響因素和數學模擬，減少機械通風期間風機噪音和糧食水分損失，保持糧食品質。

關鍵詞：淺圓倉 立筒倉 氣流分布 孔隙度 曲折度 糧食品質

我國糧庫建設用地上，堅持節約用地原則，當前我國現代化糧倉中，高大平房倉占比例85%，立筒倉占7%，淺圓倉占5%。從倉儲和貿易角度看，淺圓倉和立筒倉是發展方向。國外近年研究顯示，農場筒倉中谷物糧食的堆密度，在筒倉底部能夠增加8%-10%，是由于上面糧食產生的靜態壓強，增加了氣流通道的曲折度[1]。圓柱形筒倉倉頂中心入糧方法，糧粒從糧倉中心位點進入，徑向向外流出，導致絕大多數的粉塵材料保持在糧堆的中間核心。沿著筒倉直徑向外，粉塵的數量減少[2]。這種中心入糧的方法導致在糧堆的中心位置較低的孔隙度，而在倉壁處則逐漸轉變為較高的孔隙度。在高大的筒倉中油菜籽儲藏期間，糧堆上層對底層施加靜態壓強，可引起底層油菜籽變形和破壞，在模擬試驗中，20-60kPa壓強導致13%-16%含水率的油菜籽在25-30℃范圍隨著壓強、含水率及溫度增加，總植物甾醇、β-谷甾醇、蕓苔素甾醇的含量降低，而7%-9%含水率的油菜籽隨著壓強和溫度增加植物甾醇含量變異小[3]。值得分析淺圓倉和立筒倉中糧堆氣流分布和品質影響因素，為優糧優儲提供技術保障及評價指標。

一、我國現代化糧倉建設歷程及糧堆氣流分布研究的目的

自從20世紀80年代改革開發以來，我國調動了農民種植糧食的積極性，糧食總產量由1978年的3億噸持續增加到1998年的5.1億噸。到1998年由于中央和地方儲備糧食的倉容不足，導致農民買糧難，影響了農民種植糧食的積極性，糧食總產量下降到2003年的4.3億噸，在2003-2016年期間我國加大了糧庫倉容建設，同時對農民種糧積極引導，到2007年糧食產量又突破了5億噸，當前糧食年產量保持在6.3億噸左右[4-5]。上世紀50年代到1998年之前，我國以蘇式倉、基建平房倉、土圓倉及地下倉等倉型為主，儲糧基礎參數來自前蘇聯，推廣自然低溫、低溫密閉及地下低溫儲糧技術。在80年代開始示范推廣自然低溫輔助機械通風、谷物冷卻機、制冷機及氣調等技術，均能夠延緩糧食品質劣變[6]。在1998年我國開始建設以高大平房倉型為主的現代化儲備糧倉，單倉倉容為5千-3萬噸[7]。2003-2016年為現代化糧食儲備倉建設的重要時期。目前各種倉型中，高大平房倉占比例85%，適合長儲藏周期。在第二次世界大戰之后，美國國家儲備糧采用高大平房倉儲藏，在20世紀50-90年代，美國農業工程學會雜志（Transactions of the ASAE）發表了大量關于高大平房倉糧堆微生態學的參數及變化規律[8-9]。由于高大平房倉受大氣濕度影響而倉內局部糧情不穩定，我國在“十五”到“十二五”期間對高大平房倉儲藏的糧堆溫度及水分變化規律研發投入很大。但是，對糧堆通風期間氣流分布和壓強降一直缺乏深入研究。

糧堆中存在的粉塵、雜質，以及糧堆構造、糧食的物理特性，均可能引起倉內局部區域糧堆通風暴露率差，通風暴露率高的局部區域通風過度而造成糧食過度干燥。通風效果差相關于糧堆局部區域水分升高，霉菌生長，昆蟲進攻和糧食敗壞[10-11]。研究不同糧食床條件的氣流分布和模式，高效化通風作業，以阻止糧食品質劣變、真菌毒素污染發生及減少經濟損失。

二、我國現有糧食行業標準分析

在糧倉設計和機械通風作業風量、風壓選擇上，在我國儲糧機械通風技術規程LS/T 1202-2002[12]中，風機壓強的計算方法如下：

LS/T 1202-2002對單位通風量的規定是，緩速通風時，可采用排風扇通風系統，q應該小于8 m3/（ht）;房式倉或淺圓倉選用離心風機或軸流風機的通風系統，q應該小于20 m3/（ht）;立筒倉選用離心風機的通風系統，q應該小于10 m3/（ht）。

全球徑流數據中心（GRDC）[13]建議，干燥高水分糧食到儲存安全水分，對谷物含水率小于14%，需要的氣流速率在54-90 m3/（ht）;通風冷卻糧食，對谷物含水率12%-14%，需要的氣流速率7.2-14.4 m3/（ht）。對100噸筒倉，這些值等價于表觀風速（空氣體積流速除以筒倉橫切面積），對糧食干燥表觀風速是0.07-0.11 m/s，對糧食冷卻表觀風速則是0.01-0.02m/s。

國內在糧堆降溫通風作業中選擇風機功率和風壓時，在計算中壓緊系數C1取為1.3-1.5，風道阻力取為200-300 Pa;總風機壓強為壓緊壓強與風道阻力之和再乘以系數1.2。這些計算基于經驗，過于粗糙。由于糧食粉塵、雜質含量的影響，通常選擇大風量風機，導致降溫通風作業時，糧堆水分丟失嚴重;或者選擇的風機風量過小，通風不均勻，局部糧堆仍然存在高溫區，沒有達到降溫通風的預期效果，導致儲糧害蟲生長繁殖。還有設計立筒倉或淺圓倉時，當前的經驗是對有的糧種如大豆，倉房高度不能夠超過24米。這些問題的產生主要是由于國內對糧堆內壓強分布缺乏研究。

三、糧堆內壓強分布研究進展

糧堆主要的非生物因子有溫度、水分及籽粒間隙氣體等。糧堆溫度和水分的變異主要導致水分遷移和熱斑點產生，而籽粒間隙空氣充當熱和水分遷移的載體。研究糧堆中氣流分布，有助于設計科學的通風系統和保持糧食品質[14-16]。

任何結構的糧倉中，糧堆內氣流發生自然對流和強力對流兩種類型。自然對流主要由溫度變異產生，而強力對流則是采用風機將空氣壓入或推出糧堆。在強力對流作業中，低氣流速率3.6-7.2 m3/（ht）用于冷卻糧堆，去除糧堆中的溫度梯度，而高氣流速率36-90 m3/（ht）用于干燥糧堆。強力對流作業表現效果依賴于糧堆中空氣的均勻性，而時常存在機械通風的“死區”。在大倉容糧倉進行實倉試驗勞動強度過大，于是依據物理學原理的數學模型用于研究糧堆中氣流的分布[17-18]，這是最佳的選擇。糧堆氣流分布通常認為是流體流過多孔的介質，糧堆氣流模型則用于描述氣流壓強降與氣流速度的關系。氣流分布模型可分為簡單和高級模型。在簡單模型中糧堆被認為是同質的或各向同性，在高級模型中糧堆則被考慮為各向異性或非均質的特征。

（一）糧堆氣流與壓強降之間的關系

為了準確模擬氣流通過筒倉中糧堆，模擬了各向異性的多孔介質。許多學者采用有限元（FE）和有限體積（FV）方法[21-23]，以及Ergun方程（方程7）模擬了筒倉中一個元件的每單位長度的壓強降。基于流體流動阻力的雷諾爾德（Reynolds）理論，方程7中壓強降是粘性阻力（第一項）和惰性阻力（第二項）的函數。

粘性阻力項解釋了粘性流體穿過多孔介質的壓強損失，它與表觀風速成比例;該項顯示，較大的顆粒尺寸，發生較少的阻力，而孔隙度增加則減少粘性阻力。惰性阻力項是穿過多孔介質的空氣密度的函數，它相關于當氣流遷移通過糧堆曲折途徑空隙期間的惰性損失。它隨表觀風速的平方增加，不依賴于顆粒尺寸和孔隙度。Lawrence and Maier [21]發現，粘性阻力和惰性阻力成分需要分別乘以因子3和1.7，以匹配儲存玉米的壓強降。需要修正這些因子包含進入阻力成分的糧堆曲折度參數。因此，當采用Ergun方程，通過糧堆的任何氣流研究必須確定合適的修正因子。以前的數值研究采用多孔介質元件內具有表觀風速的質量和動量守恒定律。但是，在實際中，多孔介質的空氣穿過孔空間更加快速。

（二）簡化假說和邊界條件

通過綜合流體流動方程和連同體方程，以Laplacian形式描述了氣流分布，報道的模型滲透率（K）不同，K被計算。一些科學家根據他們的方便性、可利用性及編程技術，以不同的計算機語言寫出他們自己的計算機編碼。一些學者修正可利用的程序以滿足研究和應用的目的。對壓強P解方程，畫出等壓線和蒸汽線，以研究他們的氣流模型。

邊界條件幾乎對所有模型是一樣的，（1）糧食倉房壁被認為是不透氣的，采用Neumann條件，即？犖P=0;（2）空氣入口和出口點被認為是等壓值，采用Dirichlet條件，即P=Pe。 在大多數模型中，空氣存在的糧堆外部表面被認為是大氣壓強（Pe=0），在少數情況，一些常數值被采用。

（三）影響糧堆氣流分布的因素

1. 糧堆孔隙度和曲折度

在筒倉地板之上，糧堆內氣流分布依賴儲糧的孔隙度和曲折度。孔隙度是糧堆內空氣比例的體積測定，而曲折度是測定空隙的連通性。筒倉中心入糧的方法導致在糧倉的中心較低的孔隙度，在倉壁處則逐漸轉變為較高的孔隙度。Lawrence and Maier[21]分析直徑10m、高度8m 的筒倉內玉米糧堆，核心處孔隙度為0.34，而近倉壁處孔隙度為0.38。Olatunde等[22]分析直徑15.63m、高度4m 內稻谷糧堆，糧堆核心處孔隙度為0.45，而近倉壁處孔隙度為0.55。結果顯示，較低孔隙度的糧堆對穿過它的氣流提供較大的阻力。

Atungulu等[25]研究新收獲的四個品種稻谷在含水率18%-27%濕基范圍，其中扣除物含量范圍是0.2%-2.0%，扣除物含量依賴于收割機和天氣條件。當稻谷送入儲存倉，在倉中心稻谷糧堆低孔隙度，而倉周邊較高孔隙度。這導致通風期間倉中心顯示較低的氣流速度，而倉周邊糧堆較高的氣流速度。Siebenmorgen and Jindal[26]認為氣流阻力來自粉塵的含量，他們采用含水率12%、18%、24%的長粒稻谷品種，粉塵濃度是0%、5%、10%、15%、20%、25%、30%， 氣流速度范圍是0.0135-0.387 m/s， 發現粉塵濃度增加1%，氣流阻力增加0.87%。Chung等[27]獲得相似的結果，采用的稻谷含水率12%、13%、15%、16%、18%，粉塵含量0%、1%、3%、5%，氣流速率0.05-0.38 m/s。

2. 籽粒形狀

不同種類糧食的形狀范圍從近似球形如油菜籽到縱長型，且長度大于寬度和厚度，如小麥。對球形和縱長型籽粒，立筒倉內糧堆對穿過它的氣流產生各向異性的阻力，水平阻力小于垂直方向的阻力（如表2）。當孔隙度不變時曲折度隨著方向有差異[28]。

3. 糧堆高度

影響立筒倉內氣流均勻性的最后一個因子是糧堆的高度。倉頂中心填充糧食的圓柱筒倉，產生中心峰，以糧食的靜止角（稻谷的自然靜止角30°）升高，空氣總是采取最小阻力的途徑。這意味著空氣趨向流動到糧堆周邊部分，避開了通過頂峰的抬高[10]。解決這個問題的一個方法是，通過人工弄平糧堆或者采用去核心工藝設備最小化頂峰，從筒倉騰空一些糧食，減少糧堆頂峰的高度。除去糧堆頂峰改善了氣流通過上部糧堆的均勻性[21-22]。

Lai[29]在圓柱體糧食床中采用Ergun方程，給出三維非線性偏微分方程描述穿過多孔介質的軸對稱的氣流分布，模擬糧堆核心和周邊的氣流分布采用的空隙度分布是0.4和0.6，Ergun方程（方程7）作為氣流阻力的源項添加到控制方程中。Singh and Thorpe[30]采用有限差異方法對帶有圓錐頂的糧堆提出了三維熱、質量、動量（自由對流）傳遞模型，模擬儲存在澳大利亞類型沙坑糧倉（bunkers）中的糧食，復雜的帶有圓錐頂的糧倉幾何學轉化為立方體計算結構域，他們采用Darcy定律模擬糧堆中氣流阻力，采用矢量勢能概念解非線性方程。

Smith[31]采用基于Ergun方程的壓強和速度關系預測穿過糧食介質的氣流速度，將非線性動量方程簡化為曲線方程。Garg （2005）采用有限體積方法研究糧堆中非均勻氣流分布，Lai的可變孔隙度概念被用于發展儲糧的非均勻氣流模型，在糧堆中心核區域（高雜質含量）和周邊區域（低雜質含量）采用了兩種孔隙度。

Bartosik and Maier[10]試驗測定帶有圓錐頂、平整、去核心的玉米糧堆中氣流分布，對每種糧堆測定糧面中心附近和周邊的氣流速率。基于非均勻因子（NUF）確定氣流分布的非均勻性，在玉米筒倉中，NUF定義為， （糧堆周邊風速-中心風速）*100/（糧堆周邊風速+中心風速），發現將帶有圓錐頂的糧堆弄成平整的糧面，非均勻度因子由89%減少到36%。

Lawrence and Maier[21]采用非均勻度因子驗證Bartosik and Maier[10]的工作，假定從糧堆中心到周邊孔隙度是線性變化，采用基于多孔介質理論的fluent 計算流體力學軟件解3D非均勻氣流模型，針對玉米筒倉中糧面帶有圓錐頂、平整、反轉的情況。他們對玉米糧堆采用兩種恒定的孔隙度（0.38和0.40）和三種可變的孔隙度（0.34-0.38、0.36-0.38、0.38-0.40）驗證氣流分布，對可變的孔隙度0.34-0.38，模型預測非常接近試驗結果。采用線性插值函數估計的糧堆中心孔隙度0.34，糧堆周邊孔隙度0.38。

Olatunde et al 2016[22]采用有限體積方法模擬評價筒倉中籽粒2.94mm、孔隙度0.55的長粒稻谷糧堆，帶有圓錐頂、反轉、平整糧面情況下的氣流分布特征，對氣流速率0.55、0.825、1.1 m3/（min·t）進行模擬，采用試驗尺度壓強降系統和倉容700噸的帶有圓錐頂糧面的稻谷筒倉進行驗證，長粒稻谷品種的粘性阻力和惰性阻力系數分別是9.72E+06和36185。帶有圓錐頂或反轉糧面的糧堆氣流呈現非均勻分布，帶有圓錐頂糧面的糧堆氣流阻力最大，圓錐頂處的氣流速率顯著低于糧堆周邊部位。從稻谷筒倉中測定的平均非均勻因子是34%，對恒定和可變孔隙度從模型獲得的非均勻因子分別是19%和71%。為了保證氣流分布的非均勻因子小于50%，對反轉糧面的情況需要移走最大50噸稻谷。

另外，通風期間氣流分布影響糧堆熱傳遞。大多數研究假定通風期間氣流均勻分布，而Garg （2005）研究了糧堆中由于非均勻氣流分布引起了的兩維熱、質量、動量傳遞。他采用基于軟件fluent的有限體積方法模擬非均勻動量傳遞，采用2D PHAST-FEM 編碼（Post-Harvest Aeration & Storage Simulation Tool）作為非均勻氣流分布的函數模擬熱和質量傳遞。

四、展望

（一）研究淺圓倉和立筒倉內糧堆氣流分布影響因素

有幾個相互關聯的物理因子糧堆孔隙度和曲折度、糧粒形狀、糧堆高度等影響儲糧的氣流均勻性，能夠導致差的氣流區域。為了減少儲存損失和成本，需要了解淺圓倉和立筒倉內氣流分布，最佳化通風期間的操作。這樣，數值模擬很重要，顯示如何最佳化糧食筒倉設計和操作最小化成本。

（二）機械通風期間實時顯示糧堆各點的壓強和氣流速率

由于受大都市的溫室效應影響，儲糧存儲周期內在降溫通風作業中水分損失嚴重，風機噪音大，根本原因是對糧堆氣流壓強場缺乏研究。本研究根據我國淺圓倉通風系統和糧食儲藏周期長的特點，實時顯示糧堆各點的壓強和氣流速率，有助于在降溫過程中合理選擇風機，減少風機噪音，減少糧食水分損失，節約電能，保持糧食的品質，為民眾提供優質的原糧食品和美麗的宜居環境。

參考文獻：

[1]Bhadra R， Casada ME， Turner AP， et al. Stored grain pack factor measurements for soybeans， grain sorghum， oats， barley， and wheat. Transactions of the ASABE， 2018， 61： 747-757.

[2]Jian F， Narendran RB， Jayas DS. Segregation in stored grain bulks： Kinematics， dynamics， mechanisms， and minimization- A review. Journal of Stored Products Research， 2019， 81： 11-21.

[3]Gawrysiak-Witulska M， Rudzinska M， Wawryzyniak J， Siger A. The effect of temperature and moisture content of stored rapeseed on the phytosterol degradation rate. Journal of the American Oil Chemists Society 2012， 89： 1673-1679.

[5]顧莉麗. 中國糧食主產區的演變與發展研究. 吉林農業大學博士學位論文2012， pp 22-47.

[6]王帥，趙秀梅.中國糧食流通與糧食安全關鍵節點的風險識別. 西北農林科技大學學報（社會科學板）2019，19（2）： 124-132.

[6]國家糧食儲備局儲運管理司.中國糧食儲藏大全. 重慶： 重慶大學出版社， 1994.

[7]國家糧食局糧食行政管理司. 儲糧新技術教程. 北京： 中國商業出版社， 2001.

[8]李興軍， 王雙林. 我國糧食平衡水分測定及應用研究進展. 中國科技成果2017（04）：55-57.

[9]李興軍. 谷物糧堆通風的理論依據與目標. 糧食加工2017，42（1）： 28-38.

[10]Bartosik R， Maier D. Effect of airflow distribution on the performance of NA/LT in-bin drying of corn. Transactions of the ASABE 2006， 49（4）： 1095-1104.

[11a]Khatchatourian OA， Binelo MO. Simulation of three dimensional air flow in grain storage bins. Biosystem Engineering 2008， 101（2）：225-238.

[12]LS/T 1202-2002.機械通風儲糧技術規程.2002.

[13]GRDC. GRDC grain storage grow notes national. Grain Res Dev Corp. Retrieved from https//grdc.com.au/_data/assests/pdf_file/0019/244045/GRDC-Grain-Storage-GrowNotes-National.pdf. 2017.

[14] Gayathri P， Jayas DS. Mathematical modeling of airflow distribution in grain bulks- a review. An ASABE Meeting Presentation， Paper number： 076226. Minnesota， USA： ASABE， 2007-06-17

[15]Crozz DE， Pagano AM. Modeling resistance to airflow to airflow through beds of agropyron and corn. Estimation of power ventilation. Latin Am Appl Res 2006， 36： 7-14.

[16]Khatchatourian OA， Savicki DL. Mathematical modelling of airflow in an aerated soya bean stored under non-uniform conditions. Biosystem Engineering 2004， 88（2）：201-211.

[17]Molenda M， Montross MD， McNeil SG， Horabik J. Airflow resistance of seeds at different bulk densities using Erguns equation. Transactions of the ASAE 2005， 48（3）： 1137-1145.

[18]Neethirajan S， Karunakaran C， Jayas DS， White NDG. X-ray computed tomography image analysis to explain the airflow resistance differences in grain bulks. Biosystem Engineering 2006， 94（4）：545-555.

[19]張來林， 趙思孟.機械通風糧層阻力測定. 糧食儲藏1993，22（5）：6-10.

[20]尉堯方，王遠成，潘 鈺，張曉靜. 糧堆通風阻力的研究方法及阻力模型研究進展. 糧食儲藏2016，45（4）： 9-15.

[21]Lawrence J， Maier DE. Three-dimensional airflow distribution in a maize silo with peaked， levelled and cored grain mass configurations. Biosystem Engineering， 2011， 110： 321-329.

[22]Olatunde G， Atungulu GG， Sadaka S. CFD modeling of air flow distribution in rice bin storage system with different grain mass configurations. Biosystem Engineering， 2016， 151： 286-297.

[23]Binelo MO， Faoro V， Kathatourian OA， Ziganshin B. Airflow simulation and inlet pressure profile optimization of a grain storage bin aeration system. Computers and Electronics in Agriculture， 2019， 104923.

[24]Hood TJA， Thorpe GR. The effects of the anisotropic resistance to airflow on the design of aeration systems for bulk stored grains. Agric Eng Aust 1992， 21 （1 and 2）： 18-23.

[25]Atungulu G， Prakash B， Wang X， Wang T， Fu R， Khir R， et al. Determination of dockage for accurate rough rice quality assessment. Applied Engineering in Agriculture， 2013， 29（2）： 253-261.

[26]Siebenmorgen TJ， Jindal VK. Airflow resistance of rough rice as affected by moisture content， fines， concentration and bulk density. Paper/ASAE， 1987， 86： 3036.

[27]Chung D， Maghirang R， Kim Y， Kim M. Effects of moisture and fine material on static pressure drops in a bed of grain sorghum and rough rice. Transactions of the ASAE， 2001， 44（2）：331

[28]Nwaizu C， Zhang Q. Characterizing tortuous airflow paths in a grain bulk using smoke visualization. Canada Biosystem Engineering， 2015， 57： 13-22.

[29]Lai F. Three-dimensional flow of air through nonuniform grain beds. Transactions of the ASAE， 1980， 23（3）： 729-934.

[30]Singh AK， Thorpe GR. A solution procedure for three-dimensional free convection flow in peaked bulks of grain. J Stored Prod Res， 1993， 29（3）： 221-235.

[31]Smith EA. Pressure and velocity of air during drying and storage of cereal grains. Transport in Porous Media， 1996，23： 197-218.

（作者單位：國家糧食和物資儲備局科學研究院）