基于多物理场分析稻谷保质储藏的研究

刘家琦，杨开敏，余 海，王远成，杨 泰，李加斌

（山东建筑大学 热能工程学院，山东 济南 250101）

储藏阶段，粮食品质的好坏取决于粮堆的温度和水分，如果温湿度过高，将导致储粮发热霉变和害虫滋生繁衍[1]。杭州地区稻谷的安全储藏条件为：平均粮温不超过20 ℃；水分在12%～15%之间，稻谷有其原有的正常颜色，若出现淡黄色，则称为黄变[2]，当水分超过 15%时稻谷易黄变。研究稻谷的温湿度及黄变规律，可以为储粮减损保质、实现绿色储粮提供理论支撑与指导。

Arias等[3]通过数值模拟的方法研究了不同粮种在自然储藏时的温度分布和水分迁移规律。Hammamia等[4]研究了机械通风对袋装粮堆温度的影响，发现机械通风可以很好的控制粮温。尹君等[5]采用检测粮温拟合算法，通过Matlab编程模拟了不同仓型结构储粮在春、夏和冬季期间自然储藏时的温度场分布规律，未考虑水分在储藏过程中的变化，缺少在整个储藏周期内变化的研究。葛蒙蒙等[6]通过多物理场数值模拟软件COMSOL模拟了静态稻谷粮堆储藏一年的温度变化，研究发现粮温与仓外环境温度变化关系密切，仓壁处稻谷温度变化尤为明显。任芳等[7]在实验过程中发现，堆放的稻谷因呼吸放热，粮堆内部温度可达 40 ℃以上，且随着水分的升高稻谷黄变现象严重。

综上所述，关于粮堆温湿度变化的研究大多为阶段性储藏研究，针对大跨度时间周期内的演变特征研究相对较少，而粮食自入仓后开始储藏，该周期为1～3年，涉及非人工干预和人工干预两种典型储藏阶段。文章通过多物理场仿真软件模拟稻谷仓储过程，研究全年自然储藏（无通风条件）与分阶段自然储藏（有通风条件）的仓内稻谷温度、水分和黄度变化规律，为稻谷长期安全储藏提出预警策略。

1 模型的建立

1.1 物理模型

以杭州某浅圆仓为研究对象，材质为钢筋混凝土，其导热系数为1.54 W/（m·K），仓高32 m，内径25 m，装粮高度25 m，檐高27 m，仓壁厚0.27 m，仓顶厚0.14 m。因为浅圆仓为轴对称结构，可以做简化处理，沿正南正北方向的直径选取二维截面作为数值模拟研究对象，为了便于观察浅圆仓各处的温度、水分和黄度参数值的变化，在物理模型上设置了如图1（a）所示的9个监测点，两侧的监测点1、2、3、7、8、9距离墙壁0.5 m，4、5、6三点距离两侧墙壁均为12.5 m，上层监测点距离装粮面0.5 m，下层监测点距离仓底0.5 m，中间层监测点距离仓底12.5 m。图1（b）为物理模型的网格划分示意图，为了提高网格质量和计算精度，对粮仓内部壁面进行网格加密处理，网格的平均单元质量为0.94，网格总数为10 970。

图1 物理模型Fig.1 The physical model

1.2 数学模型

文章通过连续性方程、能量方程、动量方程、水分守恒方程[8-11]和稻谷黄变模型[12]来研究大时间跨度、多物理场的全年自然储藏与分阶段自然储藏过程，谷物在储藏过程中会进行有氧呼吸，所以在能量方程中考虑稻谷呼吸作用放出的热量，在水分守恒方程中考虑稻谷呼吸作用释放的水分[13]。

1.2.1 连续性方程

根据质量守恒原理，假定空气是不可压缩的，粮堆内的连续性方程由式（1）表示：

式中：ρa为空气密度，kg/m3；ε为粮堆孔隙率；t为时间，s；∇为哈密顿算子；u为空气的表观流速，m/s。

1.2.2 动量方程

粮堆内流体的动量方程由式（2）表示：

式中：ui为张量形式的达西速度；t为时间；δij为谷物颗粒间空隙距离；T为储藏期间粮堆温度，T0为粮堆初始温度，K；ρ0为温度T0时的空气密度，kg/m3；φ为稻谷颗粒当量直径，mm。

1.2.3 能量方程

粮堆内的热量传递满足能量守恒，由式（3）表示：

1.2.4 水分守恒方程

水分守恒方程由式（4）表示：

1.2.5 稻谷黄变模型

在储粮过程中，粮堆内部的环境条件，即温度和相对湿度条件下，稻谷的颜色由白色变为黄色，颜色的变化由b值衡量[12]，由式（5）表示：

其中，

式中：k为稻谷的黄变率；T为温度，℃；φ为小数级的粮堆相对湿度，%。

1.3 数据分析

通过 COMSOL软件对仓储稻谷传热传质及黄变过程进行数值模拟，得到粮堆在一年储藏周期内的温度云图、水分云图和黄度云图；利用Origin软件将各监测点的温度、水分和黄度变化数据绘制成折线图，对比分析全年自然储藏与分阶段自然储藏工况下的变化规律。

2 数值模拟参数设置

采用有限元法研究稻谷粮堆，粮堆的孔隙率为 0.48、密度为 580 kg/m3、渗透率为 7.27×10–9m3、比热容为(1269+34.89M) J/(kg·k)，其中M为粮堆的湿基水分。稻谷粮堆的初始温度为 28 ℃，初始水分为14%，黄度值为11.50，文章将对2020年10月1日至2021年9月30日实行全年自然储藏与分阶段自然储藏的两种工况进行数值模拟，具体储藏周期如表1所示。

表1 工况储藏时间设置Table 1 Setting of storage time under working conditions

数值模拟时，各边界设置为无滑移边界条件，且仓壁气密性良好，北侧仓壁外侧温湿度设置为当地大气逐时温湿度，南侧仓壁的边界条件设置为当地大气温度与太阳辐射的综合作用，2020年10月1日至2021年9月30日杭州地区的气象参数如图2所示。工况二的通风阶段，进风温湿度为外界空气温湿度，通风方式为全孔板垂直通风，进风口在仓底，仓底为绝热边界条件，2个出风口在仓顶，单位吨粮通风量为10 m3/（h·t）。

图2 杭州地区气象参数Fig.2 Meteorological parameters of Hangzhou

3 模拟结果与分析

对全年自然储藏和分阶段自然储藏的工况进行数值模拟，模拟结果温度、水分和黄度变化分析如下。

3.1 工况一模拟结果分析

图3为工况一自然储藏365 d时的仓内温度和稻谷水分、黄度云图。由图3（a）可知，粮仓内部温度高达37.9 ℃，原因是粮食发生呼吸作用释放热量；由于南侧仓壁受太阳辐射的影响较大，该侧仓壁附近的稻谷温度高于北侧仓壁处的稻谷温度。由图3（b）发现，装粮面顶部中间区域有一薄层粮食水分偏高，高达 15.1%，原因是粮仓内部存在温度梯度，形成自然对流微气流，气流携带水分自下往上迁移，导致自然储藏365 d后，粮仓底部稻谷所含水分下降，粮仓上部装粮面处水分升高。稻谷粮堆的初始黄度值为11.50，图3（c）显示，自然储藏一年后粮堆内部黄度值升高很多，高达 13.16；仓壁附近稻谷黄度也有所增加，但增加幅度明显小于内部，且南侧仓壁附近黄度值增加远大于北侧仓壁附近，综合图3中各云图特征，可推测黄度值受温度的影响大于相对湿度。

图3 自然储藏365 d的模拟结果Fig.3 Simulation results after 365 days of natural storage

3.2 工况二模拟结果分析

分阶段自然储藏的工况二在2020年10月1日至12月8日进行为期69 d的第一阶段自然储藏；在2020年12月9日至11日进行为期3 d的第二阶段通风，通风温湿度为仓外空气的温湿度，通风平均温度约为 5 ℃，平均相对湿度约为 75%；在2020年12月12日至2021年9月30日进行为期293 d的第三阶段自然储藏，各阶段的温度、水分、黄度变化分析如下。

3.2.1 第一阶段

图4为工况二在第一阶段自然储藏69 d时的各物理量云图变化。图4（a）为温度云图，发现粮仓内部粮温约为 29.7 ℃，比稻谷粮堆的初始平均温度升高1.9 ℃；受太阳辐射的影响，南侧仓壁附近粮温高于北侧仓壁附近的粮温；粮仓顶部空气区域的温度低于粮堆内部稻谷温度，原因是此时为冬季，外界环境温度低，空气热阻很小，上部空气区域易受其仓外环境温度影响，且装粮面附近的粮食也易受环境温度影响，所以其温度比粮仓内部稻谷的温度低。图4（b）为水分云图，粮仓内部平均水分仍为14%，等于稻谷粮堆所含初始水分，仓壁附近的粮食水分有所升高，且北侧仓壁附近的稻谷水分高于南侧仓壁附近，分析原因是粮仓内部存在温度梯度，存在自然对流微气流，气流携带水分自下往上，在仓顶处分别向左右两侧沿着仓壁两侧向下流动，因南侧仓壁附近温度高于北侧仓壁温度，导致北侧微气流循环大于南侧，水分迁移程度大，所以北侧仓壁附近所含水分高于南侧。图4（c）为黄度云图，平均黄度值约为11.58，比初始黄度之值增加了0.08，原因是同阶段稻谷粮堆温度有所上升，黄度值随着粮温升高而变大；南侧仓壁附近稻谷的黄度高于北侧墙壁，原因是受太阳辐射影响，南侧仓壁附近粮温高于北侧，南侧仓壁附近稻谷黄度值升高。

图4 第一阶段自然储藏结束时的模拟结果（69 d）Fig.4 Simulation results at the end of the first stage of natural storage (69 d)

3.2.2 第二阶段

图5为第二阶段通风3 d结束时的各物理量云图。图5（a）为温度云图，粮仓内部粮温约为5 ℃，粮仓底部粮温约为1.3 ℃，原因是通风采用全孔板垂直通风，通风结束前夜间大气温度偏低，导致通风道附近稻谷温度低于粮堆内部温度。图5（b）为水分云图，稻谷粮堆在进风口处和粮堆上层粮面附近水分较高，在粮面以下的中上层区域水分较低，因为通风空气相对湿度高于粮堆内部平衡湿度，而粮堆具有吸湿特性，所以在进风口处粮堆吸湿后水分升高，而在中上层，失去部分水分后的湿空气随着温度上升，相对湿度会进一步下降，明显低于所接触粮堆的内部平衡湿度，所以中上层粮堆的水分仍然处于失水状态，72 h后，粮堆的整体平均水分为13.45%，较通风前降低0.55%。图5（c）为黄度云图，与第一阶段相比，黄度值基本没变化，分析原因是稻谷黄变是一个不可逆的过程，黄度值不会因稻谷温度下降而下降，且稻谷发生黄变是一个极其漫长的过程，短时间内黄度变化不明显。

图5 第二阶段通风结束时的模拟结果（72 d）Fig.5 Simulation results at the end of the second stage of ventilation (72 d)

3.2.3 第三阶段

该指标常用于度量两个聚类结果的相近程度.设U和V是对N个样本数据的分配情况，两种分布的熵(代表不确定程度)分别为：

图6为第三阶段自然储藏293 d结束时的各物理量云图。图6（a）为温度云图，发现仓顶空气区域温度远高于粮堆内部温度，原因是此时为夏季，外界环境温度高，空气的热惯性较小，装粮面上部空气区域易受外界环境大气温度的影响。对比图6（a）与图3（a）温度云图发现，工况一自然储藏365 d后，整个粮仓处于“热芯”状态，而分阶段自然储藏的工况二，整个粮仓处于“冷芯”状态，只有仓壁附近及装粮面上部为高温区，与工况一形成了鲜明对比，证实增加通风的分阶段自然储藏更有利稻谷长期安全储藏。

图6 第三阶段自然储藏结束时的模拟结果（365 d）Fig.6 Simulation results at the end of the third stage of nat ural storage (365 d)

图6（b）为水分云图，装粮面上部空气区域温度远高于粮食温度，导致粮堆内部空气向下流动，粮堆上层的水分随气流向下输送，使粮堆底部水分升高，此时粮堆的平均水分为13.54%，第三阶段与第二阶段水分相比升幅为 0.09%，粮堆平均水分仍在安全储藏水分之内。对比图6（b）与图3（b）水分云图发现，工况二分阶段自然储藏365 d时的平均水分比工况一低约0.5%，可知增加通风阶段更利于稻谷粮堆安全储藏。

图6（c）为第三阶段自然储藏293 d结束时的黄度云图，仓壁附近稻谷黄度值12.0以上，变化明显，粮仓内部平均黄度值约为11.59，原因是该阶段通风结束时，粮食平均温度不高，低于粮堆初始温度，稻谷平均黄度变化很小；仓壁侧和装粮面附近的稻谷粮食温度高，对应的黄度值变化大。对比图6（c）与图3（c）黄度云图，发现工况一全年自然储藏结束时，稻谷粮仓内部黄度值增加约 1.0，稻谷黄变情况加重；而工况二储藏365 d结束时，南侧仓壁和装粮面附近因受太阳辐射与外界环境温度的影响温度偏高，致使稻谷黄变值增加明显，但粮仓内部稻谷黄度值增加仅为0.08，可以通过在仓壁外侧和仓顶刷涂层，减少太阳辐射，降低仓壁附近稻谷的温度，减轻稻谷黄变情况。对比工况一与工况二可知，通风可以降低稻谷温度，缓解稻谷黄变进程，有利于稻谷安全储藏。

3.3 工况一与工况二结果对比分析

3.3.1 各物理场分析

图7为全年自然储藏与分阶段自然储藏的工况在各个监测点的各物理场变化。图7（a）发现，无通风时，靠近仓壁的监测点温度与外界气温变化趋势一致，且靠近仓壁南侧的1、2和3监测点温度比靠近北墙壁附近的7、8和9监测点温度高，而粮堆内部和仓底中间的监测点温度一直上升，这是因为仓壁直接与外界空气接触，并且南墙壁长期受太阳辐射的影响，所有获得的热量将通过仓壁传递到附近的粮堆，所以南边仓壁的温度比北边仓壁的温度要高，然而浅圆仓直径很大，大气温度很难影响到粮堆内部。图7（d）发现，有通风时，通风时粮堆内各监测点温度骤降，第三阶段自然储藏时，靠近壁面的监测点温度因外界气温回升的影响逐渐上升，粮堆内部的温度因稻谷发生呼吸作用而上升。

图7 各监测点参数变化Fig.7 Parameter changes at each monitoring point

图7（b）为全年自然储藏的各监测点水分变化，发现大部分监测点的变化幅度不是很大，未超过稻谷的安全储藏水分，但装粮面中心处监测点4的含水率一直呈上升趋势，变化幅度超过了1%，超过了安全储藏水分，不利于稻谷的安全储藏。7（e）为全年自然储藏的各监测点水分变化，通风时，位于粮面的监测点1、4、7的水分升高，位于粮堆中部的监测点2、3、5、6、8、9的水分降低，且不管各监测点的水分如何变化，各监测点水分全部处于安全水分之内；第三阶段自然储藏后期，观察到粮堆中心监测点5处的水分呈下降趋势，分析原因是受微气流影响，气流自上往下流动，带走水分至粮仓底部。图7（c）可知，无通风时，监测点5、6的黄度值升高幅度很大，储藏后期，监测点1、2、3、4处的黄度值也有不同程度的上升，不利于稻谷的安全储藏。图7（f）为有通风条件的分阶段自然储藏工况下黄度值随时间的变化曲线，监测点5、6处黄度值增长幅度不超过0.15，储藏后期，监测点1、2、4黄度值有所增加，但幅度不超过 0.4，各监测点的黄度值变化均小于全年自然储藏的工况，相比之下，分阶段自然储藏更利于稻谷的安全储藏。

3.3.2 多物理场综合分析

图8为全年自然储藏与分阶段自然储藏的多物理场综合分析。图8（a）为工况一与工况二的整个稻谷粮堆的平均温度对比，无通风条件时，粮堆的平均温度呈上升趋势，最大增幅约 5 ℃；有通风时，粮堆平均温度骤降，降幅约23 ℃，储藏后期，受外界环境温度和谷物呼吸作用的影响，稻谷平均温度逐渐上升至15 ℃左右，仍处于安全储藏温度。图8（b）为平均水分对比，无通风条件时，粮堆的平均水分变化幅度很小，基本维持在14.0%；有通风时，粮堆水分明显下降，降幅约0.55%，即使第三阶段自然储藏期间，稻谷粮堆小幅度增水，上升至13.5%，仍可保证稻谷的安全储藏。图8（c）为平均黄度对比，无通风条件时，稻谷平均黄度一直在增加，储藏结束时，增幅为0.9；有通风时，稻谷黄度基本不发生变化，维持在11.59，原因是分阶段自然储藏工况平均温度降低，稻谷黄度值受温度影响，变化缓慢，全年自然储藏的工况，平均温度一直在上升，黄度值增加。

图8 工况一与工况二模拟结果对比Fig.8 Comparison of simulation results between operating condition 1 and operating condition 2

4 结论

文章通过对稻谷分别进行大时间跨度、多物理场的全年自然储藏与分阶段自然储藏两种过程的数值模拟，得到以下结论：

无通风条件自然储藏情况下，稻谷粮堆的平均温度受仓外大气温度及自身呼吸作用影响持续缓慢增加，储藏结束时平均粮温高达 34 ℃；稻谷整体的平均水分变化不大，但在装粮面附近水分受微气流运动影响明显升高；稻谷平均黄度受温度升高的影响，一直呈上升趋势，储藏结束时黄度值增至12.4，进一步说明高温区域的稻谷极易黄变，不利于安全储藏。

有通风条件分阶段自然储藏情况下，第一阶段自然储藏69 d时稻谷平均粮温为28.2 ℃，黄度值为 11.57，平均水分为 14%；第二阶段通风储藏3 d时，稻谷平均温度受进风冷空气的影响降至 5.9 ℃，黄度值仍为 11.57，平均水分降至13.45%；第三阶段自然储藏293 d时，南侧仓壁附近的稻谷受太阳辐射的影响，温度与黄度远高于粮堆内部稻谷的平均温度、黄度值，稻谷平均温度受仓外大气温度和自身呼吸作用的影响升高至16.3 ℃，黄度值基本维持在11.59，变化幅度很小，平均水分略有升高，为13.54%。对比全年自然储藏与分阶段自然储藏的模拟结果，发现分阶段自然储藏更有利于安全储藏。

综上，分阶段自然储藏的平均温度比全年自然储藏低约18 ℃，平均水分低约0.5%，平均黄度低约 0.8，低温空气通入高温粮仓既可以降低温度，又可以降低粮食水分，缓解稻谷黄变的进程，分阶段自然储藏更有利于稻谷的安全储藏；仓壁附近的稻谷极易受外界环境温度的影响，可通过在仓壁外侧刷涂层，降低太阳辐射的影响，避免出现局部粮温升高等现象。

备注：本文的彩色图表可从本刊官网（http//lyspkj.ijournal.cn）、中国知网、万方、维普、超星等数据库下载获取。