基于Fortran程序对储粮通风温度水分和害虫演替的模拟研究

王 柯 王远成 刘济洲 尹 君 尉尧方 余 海

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

世界各国因储粮害虫对粮食造成的损失非常严重，一般损失率为10%～40%[1]，所以研究储粮通风过程中害虫数量的变化至为重要。小麦在储藏期间的重要害虫主要有谷蠹、米象和玉米象等，温度、湿度、水分等生态因素会影响储粮害虫种群生长及分布。目前在实际储粮过程中控制害虫增长的方法主要有采用杀虫剂、熏蒸剂、通风降温等方式，其中，最常见的方法是采用机械通风，即利用通风机将冷空气送入粮堆内部，与粮堆进行热湿传递，降低粮堆的温度，从而抑制害虫的繁衍，降低粮食在储藏期间的损耗。

利用计算机编程来预测实仓通风过程中的温度、湿度、水分以及生物场等变化规律，是近年来兴起的一种既科学又高效的研究方法。国内外研究人员借助计算流体力学模拟软件对储粮通风过程的研究已经取得丰硕的成果。Hunter等[2]设计了一种通风降温系统，可以形成抑制害虫增长的储粮环境。Thorpe等[3]提出了储粮通风中有关粮粒内部的热湿耦合与传递、杀虫剂衰减的有限差分数学模型，得出粮食初始温度、水分含量以及通风时长等因素对杀虫剂的衰减率的影响。朱振刚等[4]研究了圆形粮仓的通风过程，对比得出立式通风形式更有利于粮食储藏。王远成等[5-7]对国内外储粮系统的数学模型、数值模拟方法进行了综述，并模拟预测了高大圆筒仓内温度场、水分场等。

利用商业CFD(Computational Fluid Dynamics)软件模拟储粮通风过程的文献已经屡见不鲜。由于各种条件的限制，粮食仓储企业使用商用CFD软件还有诸多不便。本研究基于多孔介质热湿耦合传递理论，建立了浅圆仓的粮堆内部热湿传递和流动的数学模型，采用Fortran语言编程模拟了储粮通风后粮堆内部温度、水分含量、杀虫剂的浓度衰减规律以及谷蠹的分布情况，为优化储粮通风工艺，降低粮食在储藏期间的损耗，提高经济效益提供参考。

1 储粮通风多场耦合模型及其实验验证

1.1 实验仓和实验工况

储粮通风实验仓结构及尺寸如图1所示，其中图1a、图1c分别为粮仓的正视图和俯视图。筒仓高度为4.5，直径为4.2 m。通风道采用十字型通风道，其中十字型通风口长度为1.2 m。

注：单位：mm。图1 仓型结构图及温度传感器布置图

实验粮种为小麦，初始温度为40 ℃，初始湿基状态下含水量为12%，吨粮通风量为49.68 m3/h· t。每天通风6 h，共通风10 d，利用布置在粮堆内部的探头测出粮堆温度的变化，水分变化采用埋袋取样进行水分化验。

1.2 数学模型

定性和定量分析储粮生态系统内的各种生物和非生物因素之间的相互作用是研究储粮生态系统的关键。数学模型是定性和定量分析必不可少的工具。基于动量，质量和能量方守恒方程以及多孔介质热湿耦合传递理论，构建了浅圆仓储粮通风过程中物理和生物参数变化的数学模型。通过坐标转换将锥形底筒仓映射到柱坐标下的空间变量(ξ,φ,η)中，其中ξ代表半径方向，φ代表角度方向，η代表高度方向。

1.2.1 流动方程

圆柱坐标下通风时粮仓内部流动如式(1)表示：

(1)

式中：Vr、Vφ、Vz分别为半径方向、圆周方向和高度方向的速度场；p为粮堆的压力/Pa；ξ为半径方向；φ为角度方向；η为高度方向；α、β为假设的瞬态因子；α1、α2、α3、α4为等比的经验常数；RP为瞬变因子。

1.2.2 能量守恒方程

根据局部热平衡理论，粮堆内部的能量守恒方程可由式(2)表示：

(2)

式中：T为粮堆中温度分布/K；W是粮堆中水分的含量/kg；H为空气中的绝对湿度/kg/kg；Keff为有效热扩散量/J/kg·k；Hw为粮堆吸湿而放出的热量/J/kg；ρb、ρa分别为粮堆、干空气的密度/kg/m3；hvap为自由水放出的潜热/J/kg；hs为粮堆吸收的热量/J/kg；Qr为粮堆累积耗氧量/kg。ca、(c1)r、(c1)σ和(c2)σ分别为干空气、水蒸气、液态水和粮食的比热容/J/kg·K；c2v为潜热方程中的特定常数/J/kg·K；εr为粮堆的孔隙率；m为干物质质量/kg。

1.2.3 水分守恒方程

根据质量守恒原理，粮堆内部的质量守恒方程可由式(3)表示：

(3)

式中：Deff为有效扩散率/m2/s；dm/dt为干物质损耗率。

1.2.4 杀虫剂浓度衰减

在储粮通风中加入杀虫剂，杀虫剂会附着在粮粒上从而抵抗害虫造成的损耗。由于粮食颗粒属于吸湿性多孔介质，当杀虫剂作用于粮粒时会在较短时间内产生分解，致使自身药效减小。杀虫剂浓度衰减过程由式(4)表示[3]。

Cpes=Cie(-1.386 3(rh)time×10B(T-30)/t*)

(4)

式中：Ci为杀虫剂的初始浓度；B和t*为杀虫剂的特性参数，对于谷蠹来说，分别是0.036和8.467E06。rh和T分别为粮堆内部的温湿度，从方程(4)可以看出杀虫剂浓度的衰减是与粮堆内部的温湿度相关的，因此，可以看出粮堆内部温度、湿度或水分是杀虫剂浓度(单向)耦合的。

1.2.5 储粮害虫生长模型

储粮害虫的生长及分布与谷物的含水量和温度有关。Thorpe[3]研究了昆虫种群每周的增长速率。谷蠹生长率见式(5)。

(5)

式中：Rs为谷蠹的生长率；Cr、tr分别为谷蠹的经验常数，其值分别为0.043 5、13.0，Tweb为湿球温度/K。

从热力学角度，湿球温度是与干球温度和空气湿度相关的，即粮堆温湿度是害虫生长密切相关的、相互耦合的。

1.3 数学模型的验证

1.3.1 储粮通风实验

实验粮种为小麦，通风时间为每天0∶00—6∶00，共通风10 d。实验送风温度为20 ℃，相对湿度55.5%。在通风期间，采集设置在粮堆内部的温湿度传感器检测粮温的数据，通风过程水分采用埋袋取样进行水分化验。同时利用程序进行通风工况的模拟，然后比较两者结果，验证模型的准确性。本次实验共设置11个测点，Z1、Z2、Z3布置在筒仓轴心处，在东南、东北、西南、西北不同高度各设有2个测点，NE1、SE1、NW1、SW1设置在下部；NE2、SE2、NW2、SW2设置在上部，具体传感器布置见图1b。粮库内配备一套精确的粮情监测软件，可以将某时刻测点的温湿度输入电脑端。本次实验每隔4 h采集一次温度，但由于粮食水分含量变化缓慢，所以每天取样一次测定水分含量。

1.3.2 模拟的初始条件及程序运行

模拟状态下稻谷粮堆温度为40 ℃，初始湿基状态下含水量为12%，送风温度为20 ℃，相对湿度55.5%，吨粮通风量为49.68 m3/h·t，均与实验保持一致。首先对三维筒仓进行网格划分，在ξ方向，筒仓间距是不均匀的，在其他两个方向上，网格间距等分。最终选用知21×20×21的网格划分方案，可以在满足精确度的同时减少不必要的运算量。模拟中通风方案采用间歇式，每天通风6 h共通风10 d，通风状态下循环时间步长为90 s，非通风状态下循环时间步长为900 s。将筒仓尺寸、粮食的各项物性参数以及害虫增长的经验常数等定义为全局变量，设立在子程序中，主程序运行时调用子程序即可运行。

1.3.3 数值模拟与实验结果比较

将实验实测的温度数据与模拟数据导入Origin软件进行处理后，选取其中三个测点呈现结果见图2。

图2 部分测点温度的实验数据与模拟数据对比

模拟水分数据与实测水分数据的对比见表1。

表1 模拟水分数据与实测水分数据的对比

综合图2及表1可知，探头所测出的粮堆温度与模拟结果近乎相似，水分的实测数据与模拟数据误差不超过0.5%，相对较小，验证了数学模型和这套程序的准确性。

2 不同通风工况的模拟结果与对比分析

2.1 不同通风工况的初始条件

为了进一步优化储粮通风工艺，将筒仓高度设为2.7 m，直径为2.23 m。通风道采用十字型加环型组合式通风道，其中十字型通风口长度及环型通风口半径均为0.8 m，优化后的储粮仓俯视图如图3所示。模拟初始温度为30 ℃，初始湿基状态下水分含量为12%，吨粮通风量为5.65 m3/h·t，其他参数设置同实验工况。分别模拟工况一即每天通风10 h、工况二即全天24 h通风两种工况，通风总时间均为10 d，对比分析温度场、水分含量、储粮害虫的分布以及杀虫剂浓度衰减的结果，为以后粮食的安全储藏和通风设计提供参考。

注：单位：mm。图3 优化后的筒仓俯视图

2.2 不同工况下粮堆温度及水分含量的变化与分析

自然储藏时，粮仓内粮堆的温度变化较缓慢，但对粮仓进行通风时，温度会发生较为明显的变化，以筒仓半径为横坐标，筒仓高度为纵坐标绘制不同工况的粮堆温度分布如图4a、图4c所示。粮食初温为30 ℃，两种工况下均通风10 d后，通风道附近最低温度降至13.4 ℃，相比初温将低16.6 ℃，下降较多，说明通风降温效果较好。东西两侧壁面的温度分布并不对称，这是因为在不同方向的筒仓壁面所接受到的太阳辐射不同，东、北、西和南四个方向的壁面温度径向分布如图4e、图4f所示。可以看出由于接受了高强度的辐射，南向壁面温度最高，两种工况下壁面温度径向分布规律相似，但由于工况二是连续通风，距筒仓中心较近处，壁面温度下降幅度高于工况一。在筒仓中心区域，温度从下向上逐渐增大，平均温度约为17 ℃。对比图4a与图4c可知工况二条件下粮堆温度分布更加均匀，在筒仓中上部区域温度略高于工况一，下部区域低于工况一，但相差不大。且底部通风笼附近温度均升高，可能与粮粒呼吸作用产热有关。

由图4b、图4d可知，粮堆的水分含量近似对称分布，但工况二条件下水分含量整体低于工况一，在锥形仓底部，工况二粮堆水分含量降至11.3%，比初始含水量低约0.7%，比工况一低约1%。同时，由于种子的呼吸作用的影响，底部通风道周围粮堆水分含量升高。由于组合式通风道通风面积较大，在整个通风过程中，水分整体变化较小，说明该风道形式通风保水效果较好。

图4 不同工况粮堆温度与水分含量分布等高线图

2.3 不同工况下储粮害虫的分布与对比分析

影响储粮害虫生长繁殖的因素很多，包括温度、湿度、光线等环境因素和粮粒本身的生物因素，各环境因素共同作用对储粮生物场造成综合性影响，且储粮害虫正常的生长、繁殖必须依赖一定范围的温度、湿度条件。对储粮仓进行通风能有效控制昆虫种群生长，本次主要模拟了谷蠹在粮仓内的增长量及其分布，如图5所示。对比图5a、图5b可知，两种工况下谷蠹在筒仓内的分布规律相似，筒仓中心区域均有分层现象，但工况二条件下谷蠹分层更明显。由于谷蠹具有喜温性，所以主要分布在壁面处及上层区域。在筒仓中心区域，工况二条件下同高度处谷蠹的增长量高于工况一，但相差不是很大。当筒仓高度小于1 m时，工况二的谷蠹增长量小于工况一，这是由于连续通风导致该区域粮堆温度及水分含量都低于工况一。

图5 不同工况下谷蠹数量分布图

2.4 不同工况下杀虫剂浓度衰减的变化与分析

对粮仓通风降温能够有效缓解储粮害虫的生长，但简单的机械通风很难杀死害虫，防治储粮害虫的最直接有效的方法便是使用杀虫剂。图6a、图6b分别是两种工况下杀虫剂浓度衰减等高线图，可以看出在通风道附近，靠近风口位置，杀虫剂浓度较高，最高可达97.6%和97.7%。在壁面附近及筒仓中心区域，从下向上，杀虫剂浓度逐渐降低，这是由于筒仓中心区域的气流温度较高，影响杀虫剂的降解，所以在这个区域的杀虫剂浓度较低，且工况二条件下杀虫剂浓度衰减分层更明显，整体浓度水平低于工况一。另外，太阳辐射也会影响杀虫剂浓度的衰减，由图6可知，壁面处杀虫剂浓度下降较快，这是由于太阳辐射主要作用于筒仓壁面，致使壁面温度较高，加速了杀虫剂的降解。

图6 不同工况下杀虫剂浓度衰减等高线图

3 结论

本研究基于多孔介质热湿耦合传递理论，建立了浅圆仓的粮堆内部热湿传递和流动的数学模型，采用Fortran语言编程模拟了储粮通风后粮堆内部温度、水分含量、杀虫剂的浓度衰减规律以及谷蠹的分布情况。

基于本文建立的储粮通风数学模型预测的温度和水分值与实验测定数据最大相差分别为1.1 ℃与0.5%，初步证明数学模型的可靠性。由于受到太阳辐射的影响，浅圆仓内部粮堆的水分含量近似对称分布，粮仓不同方向壁面的温度分布并不对称，朝阳一面贴壁处粮温高于被阴面。储粮害虫在粮仓内的数量分布与温度、水分等因素有关。储粮害虫在壁面处分布较多，在筒仓中心区域，害虫数量分布出现分层现象。杀虫剂浓度衰减受温度的影响，温度高会影响杀虫剂的降解，导致杀虫剂浓度较低。