储粮生态系统数学模型和数值模拟研究进展

王远成 吴子丹 李福君 曹 阳 张忠杰 赵会义

(山东建筑大学热能工程学院1,济南 250101)(国家粮食局科学研究院2,北京 100037)

储粮生态系统数学模型和数值模拟研究进展

王远成1吴子丹2李福君2曹 阳2张忠杰2赵会义2

(山东建筑大学热能工程学院1,济南 250101)(国家粮食局科学研究院2,北京 100037)

分析了储粮生态系统的各种因子及影响因素，重点介绍了国内外储粮生态系统的主流数学模型、数值模拟方法及常用的软件。通过数值模拟实例和结果，对国内外储粮生态系统中空气流动、热湿迁移、生物现象及其演化规律的数值模拟研究进展和现状进行了综述，并探讨了采用数值模拟技术对储粮生态系统研究的优势以及未来发展趋势。

粮食储藏 生态系统 数学模型 数值模拟

储粮生态系统是一个封闭或半封闭生态系统，主要由粮仓围护结构、粮堆、物理和生物因子四部分组成。该系统不仅涉及粮仓内部各种生物因子(虫霉等)、非生物因子(温湿度、水分和气体浓度等)，而且还受到仓外的局地大气环境(风速、大气温湿度和太阳辐射等)的影响。一般来说，储粮生态系统的各种物理、化学和生物现象都不是单独存在的，各种生物和非生物因素及其参量之间相互耦合相互影响。要控制和调节储粮生态系统，显然需要对吸湿性多孔介质流体力学、传热传质学、害虫和霉菌的生物学、生物化学以及化学动力学有深刻的理解。因此，储粮生态系统是一门跨学科、涉及多物理和生物场的复杂系统，储粮生态系统的研究要比单独去分析研究某一物理场要复杂得多，单纯依靠传统的试验方法很难解决这样复杂系统的科学和技术问题。

数值模拟也叫计算机模拟或计算机仿真，它是基于数学建模，结合有限元或有限容积的概念，依靠电子计算机，通过数值计算和图像显示的方法，达到对工程问题和物理问题乃至自然界各类现象进行研究的目的。因此，数值模拟又称为计算机上的实验。基于计算流体力学(Computational Fluid Dynamics，简称CFD)、数值传热学(Numerical Heat Transfer，简称NHT)的数值模拟技术是国外近年发展起来的一种研究流动、传热传质等物理过程的新方法，它可以形象地再现流动、热湿传递过程的情景，也为分析和预测储粮生态系统的演替规律提供了一个有效的手段。

1 储粮生态系统的模型

储粮生态系统模型是由多个子模型组成。主要包括：1)非通风条件下(非人工干预)仓储粮堆内部自然对流和热湿耦合传递模型。该模型针对特定地点的局地天气条件下(空气干球温度，相对湿度，风速，雨雪和太阳辐射等)及粮堆内部谷物和虫霉呼吸影响下的密闭粮仓内部温度和水分变化规律、结露以及通风时机的选择等问题进行预测分析，并确定边界条件(粮仓外围护结构)的影响。2)通风条件下(人工干预)仓储粮堆内部空气流动和传热传质模型。该模型主要预测和分析通风时粮堆内部流动阻力和通风系统的功耗、粮堆温度和水分与通风量、入口空气的温度和相对湿度以及它们之间的关系。3)通风时粮堆内部各向异性的流动阻力模型。4)基于温度和水分的虫霉种群数量模型。5)气调和熏蒸模型等。

1.1非通风情况下仓储粮堆内部自然对流和热湿传递模型

在非通风情况下，由于仓外气温的季节性变化和昼夜性变化，使得粮堆内部容易产生温度梯度。在温度梯度的作用下，粮堆内部会形成一个与温度梯度方向相同的水蒸气分压力梯度，导致水蒸气在粮堆内部的扩散，并通过粮堆的吸湿和解吸湿作用引起水分的迁移。同时，由于温度梯度而产生的粮堆内部空气的密度差，形成自然对流运动，它也会促进粮堆内部的水分迁移，从而导致仓储粮堆的水分发生迁移和再分配。

粮食在自然储藏条件下，其内部的自然对流和热湿传递方程主要基于多孔介质内部质量、动量和能量守恒关系而得到的。其中，被大家普遍接受的是Khankari等[1]以及王远成等[2]建立的控制方程，该控制方程的特点就是借助于谷物吸湿/解吸湿等温曲线(即平衡水分方程)建立了描述粮堆内部水分迁移的方程。

动量方程：

(1)

能量方程：

(2)

水分迁移方程：

(3)

方程(1)描述的是在浮升力驱动下粮堆内部自然对流流动及其阻力的张量形式的达西扩展方程，其中方程(1)右边第3项即为由于温度梯度而引起的浮升力。方程(2)左边第1项代表粮堆热量随时间的变化率，第2项代表由于自然对流作用而发生的热量传递，右边第1项是由于传导作用而引起的热量扩散(傅里叶热传导)，右边第2项表示粮堆吸湿或解吸湿的热量，右边第3项表示粮堆内部由于谷物呼吸以及虫霉的生物活动而产生的热量。方程(3)左边第1项代表粮堆水分随时间的变化率，第2项代表由于自然对流作用而发生的水分迁移，右边第1项是由于水分梯度而引起的水分扩散，右边第2和第3项是由于温度梯度而引起的水分扩散和对流迁移(即热的湿效应)，右边第四项表示粮堆内部由于谷物呼吸以及虫霉的生物活动而产生的水分。

方程(1)～方程(3)中：K为渗透率；p为压力；uj为张量形式的达西速度(表观速度)；xj为坐标；μ为空气动力黏度；ρ0为温度为T0时空气的密度；T0为参考温度；β为空气的体积膨胀系数；ρa和ρb分别是空气密度和粮堆密度；ca和cb分别是空气比热容和粮堆比热容；keff为粮堆导热系数；T为温度；Wg为粮堆干基水分；μ粮堆的谷物吸湿或解吸湿热；YH2O是粮堆呼吸24小时释放的水分；YCO2是粮堆呼吸24 h的CO2释放率；qh是呼吸热值；qw是呼吸水分值。基于温度梯度的水蒸气扩散系数DT=Deffω；基于水分梯度的水蒸气扩散系数DM=Deffσ；Deff是水蒸气在粮堆中的有效扩散系数；σ是一定温度下含水量改变引起的分压变化量；ω是一定含水量下温度改变引起的分压变化量。

1.2 通风过程中粮堆流动和传热质模型

在通风条件下，粮堆内部的流动和热湿传递方程大多采用Thorpe[3]以及王远成等[4-5]建立的控制方程。该控制方程假设粮堆是连续性的、均匀分布、各向同性的多孔介质，而且粮堆内部满足局部热平衡原理。同时，考虑粮食颗粒的吸湿和解吸湿特性，并通过吸湿/解吸湿速率方程使得粮堆内部热湿传递相互耦合。由于通风时间相对较短，谷物和虫霉呼吸作用所产生的热量和水分可以忽略不计。通风时粮堆内部流动及热湿耦合传递的控制方程如下：

动量方程：

(4)

能量方程：

(5)

水分平衡方程：

(6)

吸湿/解吸湿速率方程：

(7)

方程(4)描述的是通风时粮堆内部强迫对流流动及其阻力的矢量形式的动量方程，其中方程(4)右边第3项为黏性阻力，右边第四项为惯性阻力。方程(5)左边第1项代表通风时粮堆热量随时间的变化率，第2项代表由于强迫对流作用而发生的热量迁移，右边第1项是由于传导而引起的热量扩散，右边第2项表示粮堆的吸湿或解吸湿热量，俗称水分蒸发消耗的热量。方程(6)根据整个粮堆内部粮粒间水蒸气绝对含湿量与粮粒中含水量满足质量守恒定律而得到的，左边第1项代表粮粒间水蒸气随时间的变化率，第2项代表由于强迫对流作用而发生的粮粒间水蒸气的迁移，右边第1项是由于粮粒间水蒸气的扩散，右边第2项是粮堆中水分的时间变化率。方程(7)为仿照牛顿冷却公式型式的吸湿/解吸湿速率方程，俗称干燥速率方程。方程(4)～方程(7)中：ε为空隙率；ρa为空气密度；ρb为粮堆的容重；dp为谷物颗粒的等效直径；u为粮堆内部空气的表观速度或达西速度；p为压力；t为时间；为微分算子；ca、cb分别为空气和粮堆的比热；T为粮堆绝对温度；W为粮堆的水分；keff为粮堆的有效导热系数；μ为空气的动力黏度；w为粮粒间空气中的绝对含湿量；Deff为湿空气在粮堆中的有效扩散系数；hs为谷粒吸湿或解吸湿热；k为吸湿/解吸湿经验常数；We为平衡水分。

值得注意的是，通过数量级分析方程(5)可以发现，方程(5)右边第2项不能忽略，即粮粒的吸湿/解吸湿热量不能忽略。这是因为在通风过程中，粮堆水分会因粮粒的解吸湿而降低，而解吸湿过程中势必消耗粮堆的热力学能量(即蒸发潜热)，从而加速了粮堆温度的降低。因此，在计算粮堆内部热量传递过程时，必须要考虑吸湿或解吸湿热量，否则会产生较大的误差，甚至产生谬误。

1.3 通风时粮堆内部流动和各向异性阻力模型

由于粮堆堆积过程中自溜分级以及深层粮堆重力作用，粮堆内部的孔隙率分布是不均匀的，因此，粮堆内部流动阻力是各向异性的，即水平方向阻力与垂直方向的阻力是各不相同的。Jayas等[6]采用数值模拟的方法调查了通风过程中粮堆内的气流分布和各向异性阻力。流动和各向异性阻力模型为：

(8)

(9)

(10)

B=(Bx+By)/2

(11)

式中：p为压力；x、y为在水平和垂直坐标值；Ax、Bx、Ay、By和B为实验数据回归得到的常数。Sinicio等[7]使用了有限元方法求解方程(8)，并用使用实验室测试数据验证了模型准确性。

Hood等[8]通过研究发现，当粮粒入粮仓时，在重力场的作用下粮粒在自然堆积时，由于要满足自稳定性，对于大多数粮种来说，粮粒长轴往往接近水平方向。因此，粮堆水平方向空隙的迂曲率小于垂直方向，水平方向流动阻力通常低于垂直方向上的阻力。对于粮堆来说，在水平和垂直方向的压力梯度(每米粮堆的阻力)可以表示如下：

(12)

(13)

式中：v为粮堆内部空气的表观速度或达西速度矢量；u和v为v在x、方向上分速度值；Rx、Ry、Sx和Sy是由实验确定的经验常数[8]。

1.4 虫霉发展和分布的预测模型

1.4.1 害虫种群数量预测模型

基于粮堆温度和水分的适用于各种谷物的害虫发育和生长模型已经被初步开发用于害虫种群数量的预测和模拟。因此，将粮堆温度、水分分布与害虫种群数量模型相结合，可以预测粮堆内部害虫种群的数量。目前，害虫种群数量模型主要分为2类：静态解析模型和动态的微分模型，前者忽略了害虫的爬行移动行为，后者在模型中增加了害虫爬行移动的影响因子，充分考虑了害虫的爬行移动行为。

1.4.1.1 静态解析模型

Thorpe[9]根据Desmarchelier[10]的研究结果，分析了储粮害虫种群数量与谷物的水分含量和温度的关系，建立了储粮害虫的增长率与粮粒间的空气湿球温度成线性比例关系式，并给出了害虫种群数量模型。害虫种群数量模型[10]：

P=P0exp(rst)

(14)

式中：rs=c1(Tw-c2)

(15)

(T-Tw)+w(2502+1.809T-4.186Tw)-(2502-2.377Tw)w*=0

(16)

方程(14)中P0为初始时刻害虫虫群数量；t为时间；P为t时刻害虫虫群数量。方程(15)中：rs为害虫种群的内禀增长率；c1和c2为不同虫群特定的常数；Tw为粮堆中空气的湿球温度。方程(16)中：w为粮堆中空气的绝对含湿量；w*为相应湿球温度下的饱和绝对湿度；T为粮堆中空气的温度。通过方程(15)和方程(16)可以看出，害虫虫群数量与粮堆内部的温湿度是密切相关的。

Driscolla[11]等基于4种害虫种群(谷蠹、米象、锯谷盗和赤拟谷盗)生长率的数据，建立了一个新的害虫种群的内禀增长率经验公式：

rs=(ka+kb·r+kc·r2)exp(k1·T)+ln[k2·(Tm-T)]

(17)

式中：ka、kb和kc为虫群增长率系数；k1和k2为经验常数；Tm为害虫由生长转为死亡的极限温度；T为粮堆中空气的温度。该模型与干燥和就仓通风模拟程序相结合，可以预测了在干燥和就仓通风时粮仓内部害虫种群数量的变化规律。在热带条件下该模型预测的害虫种群数量与实验数据吻合良好，但在低温和低湿度时误差稍大一些。

需要注意的是上述模型都是假设害虫是不爬行移动的，而在现实中，他们可能会转移到更有利的生存环境中。此外，通风对害虫移动的影响也被忽略。另外，模型还忽略了不同品系的害虫的行为差异。造成这些缺陷的原因是由于缺乏相关的数据而造成的，这也是储粮害虫学家值得研究的一个领域。

1.4.1.2 动态的微分模型

Jian等[12]在实验的基础上建立了害虫种群数量分布的偏微分型式的扩散输运方程，确定了与输运方程的相关系数，并通过测定有无温度梯度的小麦储藏箱内锈赤扁谷盗种群分布数据验证了所建立的模型。

(18)

式中：Dx、Dy和Dz分别为x、y、z方向害虫的扩散系数；Bx、By和Bz分别为x、y、z方向害虫的移动速度；α为害虫繁殖和死亡的数量变化量；x、y、z为空间坐标；t为时间；n害虫种群数量。

1.4.2 霉菌生长发展的预测模型

储粮最初生长的危害真菌主要是局限曲霉和灰绿曲霉[13-14]。真菌的活性依赖于很多因素，涉及粮食的贮藏条件，如温度、水分和气体浓度等。而且不同类型粮种的主要有害霉菌并非完全相同，而且贮藏于不同环境下的同一种粮种的主要有害霉菌也不同。因此，描述霉菌活动的模型的建立都是艰难的。由于储粮真菌生长影响因素的复杂性，加之缺乏相关实验数据及大量试验样本的支撑，目前针对曲霉和青霉生长的影响因素(生物的和非生物的)进行模拟研究的文章还较少发现，这也是今后的储粮生态系统研究任务之一。

1.5 气调和熏蒸模型

气调储粮模型是描述气体在粮堆内部扩散输运过程的控制方程，通过求解气调储粮的控制方程，可以预测气调和熏蒸(氮气或二氧化碳，磷化氢气体等)过程中各种气体的浓度分布规律，并用于指导气调和熏蒸工艺的设计和优化。常见的气调模型属于组分输运方程[15-16]：

(19)

式中：方程左边第1项代表某种气体组分浓度的随时间的变化率；第2项代表由于对流运动而产生的某种气体组分迁移；方程右边第1项为某种气体的分子扩散项；第2项为粮粒对某种气体组分的吸附和虫霉活动产生的气体源项。方程(19)中：C为某种气体组分的浓度；t为时间；uj为jth方向气流速度分量；xj为坐标；Deff粮堆内部气体组分的有效扩散系数。

2 储粮生态系统常用的数值模拟软件

研究储粮生态系统的方法有现场观测方法、物理模拟方法和数值模拟方法。现场观测所得到的结果无疑是可靠的，但其需要投入较大的人力物力，成本较高，而且结果不具有可重复性。物理模拟方法是借助模拟仓进行研究的方法，但它的时空尺度较小，其结果与实际情况有较大的差距。基于计算流体动力学的数值模拟方法恰好克服了前面两种方法的弱点，它具有成本低、适用范围广、结果形象逼真的特点[17]。

目前，常用于储粮生态系统的CFD商用软件主要是基于有限体积法的FLUENT和基于有限元法的COMSOL。前者在流动、传热与传质、化学反应与燃烧、生物、食品等工程领域有着广泛的应用，但其对于变物性和复杂边界条件的问题的处理稍显麻烦，需要编制用户自定义程序。后者通过求解多物理场来实现对真实物理现象的仿真，从而到达模拟工程领域的各种物理过程的目的，而且其中的模型定义非常灵活，材料属性、源项、以及边界条件等可以是常数、任意变量的函数或者直接是一个代表实测数据的插值函数等。

3 未来储粮生态系统模拟研究的方向

尽管国内外学者不懈地探究储粮生态系统的数学模型和数值模拟方法，但是由于储粮生态系统的复杂性，目前建立的模型还不能完全适用储粮生态研究的需要，虫霉发展和分布的预测模型还有待于完善和修正。另外，储粮生态系统许多其他问题的数值模拟研究尚需开展。

3.1 虫霉种群数量预测模型的完善和改进

由于储粮中虫霉生长和发展的复杂性，现有的虫霉种群数量预测模型大多是在实验室实验数据基础上建立起来的，还缺乏大量试验样本的支撑，更缺乏实仓数据的验证，距离实际应用还有一定距离。因此，虫霉种群数量模型的完善和验证工作还有待于加强。

3.2 储粮品质的变化模型及其模拟研究

考虑到储粮最终用途是食用或者其他工业用途，因此，需要研究分析储粮生态系统中生物因子和非生物因子对储粮品质(加工和食用品质)的影响关系，建立储粮品质变化的数学模型，如储粮的干物质损耗模型、储粮的生物活性变化、储粮发芽率模型等，开展储粮品质变化规律的数值模拟研究。

3.3 粮仓维护结构和仓型结构对储粮生态系统中物理和生物因子影响的数值模拟研究

3.4 各种仓外大气环境对储粮生态中物理和生物因子影响的数值模拟研究。

[1]Khankari K K, Patankar S V, Morey R V. A mathematical model for natural convection moisture migration in stored grain [J]. Transaction of ASAE, 1995, 38(6): 1777-1787

[2]王远成,亓伟,张中涛. 圆筒仓内自然对流对粮堆热湿传递的影响研究[J]. 水动力学研究与进展, 2014,29(4):487-496

Wang Y C, Qi W, Zhang Z T. The effect of natural convection on heat and moisture transfer of stored grain in a silo [J]. Chinese Journal of Hydrodynamics, 2014, 29(4):487-496

[3]Thorpe G R. The application of computational fluid dynamics codes to simulate heat and moisture transfer in stored grains [J]. Journal of Stored Products Research, 2008, 44, 21-31

[4]王远成,段海峰,张来林, 等. 就仓通风时粮堆内部热湿耦合传递过程的数值预测[J]. 河南工业大学学报,2009,30(6):75-79

Wang Y C, Duan H F, Zhang L L. Numerical simulation of heat and moisture transfer in stored grains at near-ambient drying [J]. Journal of Henan University of Technology, 2009, 30(6): 75-79

[5]Wang Yuancheng, Duan Haifeng, Zhang Hao, et al. Modeling on heat and mass Transfer in stored wheat during forced cooling ventilation[J]. Journal of Thermal Science, 2010, 19(2):167-172

[6]Jayas D S, Sokhansanj S, Moysey E B, et al. Predicting airflow pressure patterns in canola (rapeseed) bins[J]. Canadian Agricultural Engineering, 1990, 32: 249-254

[7]Simcio R, Jayas D S, Muir W E, et al. Finite element prediction of non-uniform airflow in fixed beds of wheat [J]. Postharvest Biology and Technology, 1992, 2:51-59

[8]Hood T J A, Thorpe G R. The effects of the anisotropic resistance to airflow on the design of aeration systems for bulk stored grains [J].Australia Agricultural Engineering, 1992, 21(1 and 2): 18-23

[9]Thorpe G R. Modeling ecosystems in ventilated conical bottomed farm grain silos [J]. Ecological Modelling, 1997, 94, 255-286

[10]Desmarchelier J M. 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, 1988, 24(3): 107-113

[11]Driscolla R, Longstaff B C, Beckett S. Prediction of insect populations in grain storage [J]. Journal of Stored Products Research, 2000, 36,131-151

[12]Jian F, Jayas D S, White N D G, et al. Two dimensional diffusion of Cryptolestes ferrugineus (Stephens) (Coleoptera: Laemophloeidae) populations in stored wheat under constant environmental conditions [J]. Journal of Stored Products Research, 2007, 43: 342-348

[13]陈畅,丁伟,唐芳,等. 小麦储存水分、温度和真菌生长危害进程预测[J]. 中国粮油学报, 2012, 27(5):5-9

Chen C, Ding W, Tang F, et al. Study on the changes of moisture content and temperature and forecast of growth and spoilage process of fungi in stored wheat [J]. Journal of the Chinese Cereals and Oils Association, 2012, 27(5):5-9

[14]岳晓禹,张恒业,辛婷,等. 储粮预测微生物模型的研究进展[J]. 中国粮油学报, 2012,27(5):118-123

Yue X Y, Zhang H Y, Xin T, et al. Research progress in predictive microbiology models of stored grain [J]. Journal of the Chinese Cereals and Oils Association, 2012, 27(5):118-123

[15]Alagusundaram K, Jayas D S, Muir N W E, et al. Convective-diffusive transport of carbon dioxide through stored-grain bulks [J]. Transaction of ASAE, 1996, 39(4): 1505-1510

[16]王远成,Graham R Thorpe,赵会义,等. 储粮熏蒸过程中磷化氢浓度的数值模拟研究[J]. 中国粮油学报, 2015, 30(7):81-84

Wang Y C, Thorpe G R, Zhao H Y, et al. Study on the model and validation of phosphine fumigation in stored grains in silos [J]. Journal of the Chinese Cereals and Oils Association, 2015, 30(7):81-84

[17]王远成,张忠杰,吴子丹,等. 计算流体力学技术在粮食储藏中的应用[J]. 中国粮油学报, 2012, 27(5):86-91

Wang Y C, Zhang Z J, Wu Z D,et al. Application of computational fluid dynamics technology in grain storage systems[J]. Journal of the Chinese Cereals and Oils Association, 2012, 27(5): 86-91.

Overview of Mathematical Models and NumericalSimulations in the Stored Grain Ecosystem

Wang Yuancheng1Wu Zidan2Li Fujun2Cao Yang2Zhang Zhongjie2Zhao Huiyi2

(School of Thermal Energy Engineering, Shandong Jianzhu University1, Jinan 250101) (Academy of State Administration of Grain2, Beijing 100037)

In present paper, the various factors in the stored grain ecosystem are briefly analyzed, and then main mathematical models, the numerical simulation methods and the software used commonly for the domestic stored grain ecosystem are also introduced. The advance and current situations in numerical research on air flow, heat and moisture migration, and biological phenomenon as well as evolution laws in the stored grain ecosystem are reviewed through examples and results of the numerical simulation . The advantages of numerical simulation for the stored grain ecosystem and the future development trend are also discussed.

grain storage, ecosystem, mathematical models, numerical simulations

TK121; TS207.7

A

1003-0174(2016)10-0157-06

国家自然科学基金(51276102),国家粮食公益专项(201513001)，国家重点专项(2016YFD0400100，2016YFD0401002)

2015-02-01

王远成，男，1963出生，教授，博士，复杂系统流动和多孔介质传热传质