仓储稻谷黄变的数值模拟研究

余海，王远成，*，陈朝晖，尹君，杨开敏

(1.山东建筑大学热能工程学院，山东济南 250101；2.山东省公共资源交易中心，山东济南 250014；3.国家粮食和物资储备局科学研究院，北京 100037)

0 引言

稻谷的颜色是衡量其品质和价格的重要标准之一，但稻谷储藏一段时间后会产生黄变，其影响因素主要为温度、湿度和储藏时间。 新收获的稻谷水分含量＞20%，该水分含量下稻谷易生长霉菌，因此对新收获的稻谷采用干燥处理，将水分含量降至＜15%，方可保证稻谷安全储藏。 稻谷储藏可分为自然储藏阶段(即非人工干预阶段)和人工干预阶段，自然储藏时圆筒仓是封闭的，与外界不可渗透，又称为密闭储藏。

GRAS 等[1]和 BASON 等[2]分别通过实验研究了温度、相对湿度和空气成分对两个品种白米和稻谷颜色的影响，发现温度和相对湿度对白米和稻谷黄变速率有影响，氧气对黄变速率有轻微影响，而二氧化碳对黄变速率无影响。 PHILLIPS 等[3]研究了稻谷在干燥前的黄变现象，指出在一定时间内，散装高水分的稻谷比袋装的低水分稻谷黄变速度更快。戚禹康等[4]、俞晓静等[5]、王小萌等[6-7]和鲁子枫等[8]采用计算流体动力学(Computational Fluid Dynamics，CFD)的方法模拟了粮堆温度和水分的变化，发现两种参数互有影响。 欧阳毅等[9]研究了粮堆温度和水分的变化与真菌生长的关系，得到了粮堆储藏温度和水分与真菌生长的预测关系曲线。 刘进吉等[10]探索了高大浅圆仓储粮温度的变化，发现在壁面附近温度梯度较大。 唐为民[11]分析了黄粒米的成因并提出了一些预防方法。 万立昊等[12]在探究稻谷储藏质量影响因素的过程中发现，黄粒米随温度和水分含量升高而增多。

虽然学者们对稻谷黄变进行过大量的实验研究[1-3，12]，结果较为准确且合乎实际，但是稻谷的黄变周期很长，且在高温高湿的情况下才容易黄变，随着时间的积累，黄变将越来越严重。 国内对稻谷黄变的研究与国外相比相对较少，对安全储粮的研究主要停留于温度和水分的变化规律，少数研究稻谷黄变还只是通过实验进行的，但是通过实验研究将耗费大量的时间和人力物力。 近年来，随着电脑技术的高速发展，涌现出一种用电脑模拟流体流动、传热传质的新技术，即CFD 的方法，为粮食储藏问题提供了一个良好的工具[13-14]。对比大量的研究现状发现，数值模拟时主要为设定常壁温边界条件，然而大气温度是变化的，因此采用此条件模拟的结果与真实粮情有偏差。

文章建立了仓储稻谷自然储藏过程中粮堆内部流动和热湿耦合传递的数学模型，并结合已有的稻谷黄变动力学模型，通过多物理场仿真软件COMSOL Multiphysics 数值模拟了稻谷黄变过程，温度边界条件设置为动态变化的大气环境温度，模拟环境更加接近真实粮情，研究分析了局地气候条件下仓内粮堆局部温度、水分的分布规律以及稻谷黄变规律，以期为安全储粮提供精确的理论指导。

1 稻谷黄变模型的建立与条件设置

1.1 稻谷黄变物理模型

以中试圆筒仓为研究对象，由于仓壁厚度不大，通过仓壁的传热影响较小，为了简化模型未考虑圆筒仓壁厚和其对粮食的热传导。 圆筒仓的高度为6 m、筒仓半径为1.5 m。 粮食在储藏过程中，内部粮情参数是动态变化的，因此在粮堆内部设定监测点来记录参数变化情况，其中监测点1 位于粮堆左下角(-1.25 m，0.25 m)处，监测点2 位于粮堆中心(0 m，3 m)处，而监测点3 位于粮堆右上角(1.25 m，5.75 m)处，位置示意图如图1 所示。 由于圆筒仓为轴对称图形，沿直径各截面都一样，因此选取沿筒仓直径的截面作为研究对象。 数值计算区域的网格划分如图2 所示，并且对粮仓壁面进行网格加密处理。

图1 监测点位置示意图

图2 计算区域网格图

1.2 稻谷黄变数学模型

1.2.1 连续性方程

密闭储藏过程中，由于浅圆仓是封闭的，与外界不可渗透，并假定空气不可压缩，根据质量守恒，连续性方程[15]由式(1)表示为

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

1.2.2 动量方程

考虑浮力效应，空气密度随温度的变化近似于布森内斯克方程(Boussinesq equation)[16]，根据动量守恒，其方程由式(2)表示为

式中P 为多孔介质中的压力，Pa；μ 为流体的动力黏度，Pa·s；g 为重力加速度，m/s2；β 为空气的体积膨胀系数，β ＝ 1/T0，K-1；T0和 T 分别为粮堆初始温度和实时温度，K；d 为稻谷颗粒当量直径，mm；K 为粮堆的渗透率；ui(i ＝ 1，2) 为空气的流动速度，m/s；δij为粮粒空隙间距，m；xi为坐标方向，x1＝x 时，为水平方向，x2＝ y 时，为垂直方向。

1.2.3 能量方程

密闭储藏过程中，由于浅圆仓是封闭的，与外界不可渗透，根据能量守恒定律，能量方程由式(3)表示为

式中 ρs为粮堆的密度，kg/m3；Cs和 Ca分别为稻谷和空气的比热容，J/(kg·K)；keff为粮堆的有效导热系数，W/(m·K)；Wg为稻谷粮堆的干基水分，为吸湿或解吸湿热，J。

1.2.4 水分守恒方程

自然储藏时粮堆内的水分迁移满足质量守恒，由式(4)表示为

式中ω 为稻谷颗粒间空气的绝对含湿量，kg/kg；Deff吸湿产生的干基水分质量，kg。

1.2.5 稻谷黄变模型

稻谷在仓内储藏时，在温、湿度的共同影响下，稻谷将慢慢变黄，黄度由b 值计量[17]，由式(5)表示为

其中，

式中φ 为相对湿度，%。

1.3 数值模拟参数设置

模型采用有限元法研究稻谷，粮堆的孔隙率ε 为0.48、密度ρs为580 kg/m3、渗透率K ＝ 7.27 × 10-9m3、比热容Cs为(1269 ＋ 34.89M) W/m·k-1，其中M 为粮堆的湿基水分，M ＝Wg/ (1 ＋Wg) × 100。 文章对不同温度和相对湿度的9 种工况进行模拟研究，具体参数见表1，其中恒温恒湿工况1、2、3 的模拟为对文献[17] 实验结果的验证(该工况数据源自SOPONRONNARIT 的实验数据)，其他工况为秋冬季(2016 年 10 月～2017 年 2 月)与春夏季(2017 年4 月～2017 年 8 月)各 5 个月的工况模拟，秋冬季与春夏季温度为大气温度(简称气温)(数据取自浙江省某粮库记录的温度)如图3 所示。

图3 秋冬季与春夏季温度图

表1 9 种模拟工况参数表

2 模拟结果与分析

2.1 工况 1、2 和 3 模拟结果分析

SOPONRONNARIT 等[17]采用实验的方法测量了稻谷黄变的黄度，步骤为(1) 将饱和盐溶液装在一个玻璃材质的器皿中，依次在器皿上方铺上一层塑料网和初始湿基水分为19.35%的稻谷；(2) 将该玻璃器皿放置在塑料箱内，将塑料箱进行密封后放于烤箱中。 实验时，将烤箱温度、饱和盐溶液相对湿度分别调为 35 ℃与 80%、60 ℃与 80%、60 ℃与95.2%，实验总时长为6 d，用型号为JP7100 的色差仪测量稻谷的黄度，初始黄度为11.50。

文章基于建立的数学模型，采用 COMSOL Multiphy 软件模拟上述3 种实验工况，将模拟结果与实验测量的结果对比如图4 所示，两者的变化趋势基本相同，最大误差＜3%，说明通过数值模拟的方法研究稻谷黄变是可行的。

图4 模拟结果与实验结果的比较图

3 种工况下稻谷黄度的分布结果如图5 所示。可以看出，恒温恒湿条件下，仓内稻谷黄变分布大致均匀；在相对湿度恒定的情况下，温度越高，黄度值越大；在恒温情况下，相对湿度越大，稻谷黄度值也越大，结果与图4 一致，说明稻谷在高温高湿的环境中黄变越严重，稻谷的黄变与时间呈线性关系，且随着时间的增加黄变越来越严重。

图 5 6 d 时工况 1、2、3 的黄度分布图

2.2 秋冬季与春夏季各工况模拟结果分析

2.2.1 自然对流速度场分布

秋冬季和春夏季工况下90 d 时的速度场分布结果如图6 所示。 自然对流流动速度大小的数量级大概为10-5～10-4m/s，秋冬季工况下，在对称轴(x ＝0)右侧，微气流为顺时针流动，对称轴左侧，微气流为逆时针流动，且壁面附近与粮堆内部流动相对更强；春夏季工况下，气体流动方向与秋冬季工况相反。

图6 秋冬季与春夏季工况90 d 时速度场图

2.2.2 粮堆温度分布

秋冬季1 与春夏季1 工况下温度场及温度变化图如图7 所示。

图7 秋冬季1 与春夏季1 工况下温度场及温度变化图

由图7(a)和(b)可以看出，整个储藏周期内，秋冬季1 工况下粮堆平均温度由20 ℃下降至14.5 ℃，其降幅为5.5 ℃；春夏季1 工况下粮堆平均温度由 0 ℃上升至 28.2 ℃，其升幅为 28.2 ℃。150 d时，秋冬季工况粮堆处于“热芯粮”状态，春夏季工况时粮堆为“冷芯粮”状态，这是由于秋冬季时外界环境温度较低，春夏季时其温度较高，而仓壁附近的粮食与仓顶粮面受外界环境的影响较大，因此秋冬季时仓壁与仓顶附近的粮温较低，粮堆内部温度相对较高；春夏季时仓壁与仓顶附近的粮温较高，粮堆内部温度相对较低。 由于外界温度呈波动状，温度的传递有一定的延滞性，且粮食本身导热系数比较小，所以外界温度对粮堆影响区域较小。 由图7(c)和 (d)可以看出，秋冬季1 工况下3 个监测点温度均呈下降趋势，春夏季1 则呈上升趋势，由于监测点1 和3 靠近仓壁，而仓壁附近受外界环境影响很大，因此监测点1 和3 的温度变化趋势与大气温度的变化趋势大致一致。 同时还可看出，无论是秋冬季还是春夏季，监测点3 处的温度变化最大，在春夏季时，监测点3 处温度最高。

2.2.3 粮堆湿基水分分布

秋冬季1 与春夏季 1 粮堆初始湿基水分为14%，150 d 湿基水分场及水分变化图如图8 所示。由图8(a)和(b)可知，秋冬季时，仓顶粮面及仓壁附近湿基水分较高，粮堆内部和仓底湿基水分较低；在春夏季，仓顶粮面及仓壁附近、粮堆底角处湿基水分较低，粮堆内上半部水湿基水分较高；由图8(c)和 (d)可知，秋冬季1 工况下，监测点1 和2 处湿基水分均呈降低趋势，其中监测点1 处湿基水分下降最快，监测点3 处湿基水分略微升高；春夏季1 工况时，监测点1 处湿基水分降低，监测点2 和3 处湿基水分略微升高，这是因为在自然对流的作用下，秋冬季时湿基水分从粮堆内部和底部向仓壁和仓顶粮面迁移，其中仓底附近湿基水分最低；春夏季时湿基水分从仓壁和仓顶粮面向粮堆上半部迁移，其中粮堆左下角与右下角湿基水分最低。 同时还可知，无论是秋冬季还是春夏季，监测点3 处都是高水分区。

图8 秋冬季1 与春夏季1 工况下湿基水分场及水分变化图

2.2.4 粮堆内部黄度变化

不同初始湿基水分时，各工况不同时间点下计算得到的粮堆平均黄度值，见表2。

表2 秋冬季与春夏季各工况不同时间点的模拟黄度表

由表2 可知，当稻谷初始湿基水分为14%时，秋冬季1 工况和春夏季1 工况的稻谷平均黄度值均变化较小，但春夏季和秋冬季相比，春夏季环境下稻谷黄变略微严重。 当稻谷初始湿基水分为16%时，在秋冬季2 工况下，稻谷的平均黄度值由11.50 上升至11.54，在春夏季 2 工况下则由 11.50 上升至11.62。 当稻谷初始湿基水分为19%时，秋冬季3 工况下的稻谷平均黄度值由11.50 上升至11.64，和秋冬季1、2 工况相比，稻谷黄变程度相对更加严重；春夏季3 工况下，稻谷的平均黄度值由11.50 上升至11.87。 表明初始湿基水分和气温越高，则稻谷黄变程度越严重。 另外，从表2 还可知，相同初始湿基水分时，前期时间段内秋冬季工况下比春夏季工况下的稻谷平均黄度值大，这是因为前期时间段内，秋冬季工况下稻谷的平均温度高于春夏季，但到了后期，春夏季工况下稻谷的平均温度高于秋冬季，因此到150 d 时，春夏季工况下的稻谷平均黄度值更大。无论是秋冬季工况还是春夏季工况，稻谷的初始湿基水分越高，稻谷的黄度速度越快，与初始湿基水分为14%和16%相比，当稻谷的初始湿基水分为19%时，稻谷的黄变速度明显更快。

秋冬季与春夏季工况下3 个监测点处的黄度变化如图9～11 所示，由图9～11 的(a)可知，秋冬季工况下，3 个监测点的黄度值最终相差较小。 由图9 ～11 的(b)可知，春夏季工况下，监测点1 和2 的黄度值变化也较小，监测点3 处的黄度变化较大，在春夏季1 工况(初始湿基水分为14%)，监测点3 处的黄度由11.50 上升至11.67，升幅为0.17，在春夏季2工况(初始湿基水分为16%)，监测点3 处的黄度由11.50 上升至 11.95，升幅为 0.45，在春夏季 3 工况(初始湿基水分为19%)，监测点3 处的黄度由11.50上升至12.8，升幅为1.3，说明在春夏季工况下，监测点3 处即靠近仓顶和仓壁交界附近区域黄度受影响较大，结合图7(c)和 (d)、图8(c)和 (d)也可知，监测点3 处温度和水分变化比较大，而黄变与温湿度密切相关，在高温度和高水分的双重影响下，该区域的稻谷黄变更加严重，因此，靠近仓顶和仓壁交界附近区域为危险区域。

图9 初始湿基水分为14%时各监测点黄度变化情况图

图10 初始湿基水分为16%时各监测点黄度变化情况图

图11 初始湿基水分为19%时各监测点黄度变化情况图

3 结论

文章采用有限元法数值模拟了圆筒仓内稻谷自然储藏过程中粮堆温度、水分和黄变的变化规律，分析了3 种恒温恒湿工况(温度35 ℃和相对湿度80%、温度60 ℃和相对湿度80%、温度60 ℃和相对湿度95.2%)，以及春夏季与秋冬季大气温度环境6种工况，得出以下结论:

(1) 稻谷的黄变情况受储藏环境的影响较大，恒温恒湿工况下，稻谷的黄变与时间呈线性关系，随着时间的增加黄变越来越严重；低温环境下稻谷黄变缓慢，高温环境中稻谷黄变较快，高温高湿环境下稻谷极易黄变。

(2) 在秋冬季，由于粮温较低，稻谷黄变缓慢；在春夏季，随着粮温的升高，稻谷黄变较快；在相同的储藏环境中，和稻谷的初始水分为14%和16%相比，当稻谷的初始水分为19%时，稻谷的黄变速度明显更快。

(3) 稻谷储藏过程中，靠近仓顶和仓壁交界附近区域为危险区域，该区域温度和水分变化较大且黄变严重。