浅圆仓和立筒仓粮堆气流分布的影响因素及模拟研究进展

李兴军 张洪清

摘要：浅圆仓和立筒仓占我国现代化粮仓总量的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.

（作者单位：国家粮食和物资储备局科学研究院）