仓储粮堆内热湿耦合传递的数值模拟研究

王远成 白忠权 张中涛 亓 伟

（山东建筑大学热能工程学院1，济南 250101）（河南省国防工业设计院有限公司2，郑州 450003）

仓储粮堆内热湿耦合传递的数值模拟研究

王远成1白忠权2张中涛1亓 伟1

（山东建筑大学热能工程学院1，济南 250101）（河南省国防工业设计院有限公司2，郑州 450003）

通过理论分析与数值模拟相结合的方法，以典型吸湿性多孔介质—小麦为研究对象，根据小麦吸湿和解吸湿曲线，建立了吸湿性多孔介质内部热湿耦合传递的数学模型，通过与相关试验数据比较验证了数学模型的合理性。基于有限元的方法模拟分析了外界气温和小麦分别为273 K（0℃）和293 K（20℃），小麦初始水分为14%和18%时直径为10 m高度为10 m的充满小麦的圆柱仓内部热湿迁移过程，重点探究了近似冬季和夏季仓储粮堆内部温度和水分的动态变化规律。研究结果表明，冬季工况下水分从粮堆内部向顶部和右侧壁面迁移，聚集在相对狭窄的壁面边界附近，形成低温高水分区域，最大水分值出现在右上部侧壁面附近。夏季工况下水分从顶部和侧面向内迁移，而内部又朝右上部区域扩散，形成较宽的高温高水分区域（x＜4.5 m，8.0 m＜y＜9.5 m）。

热湿耦合传递 自然对流 谷物储藏 数值模拟 有限元法

入仓储藏是粮食收获以后普遍采用的一种储存方式，粮食收获以后，通常被晾晒或机械干燥至12%－14.5%的水分含量，才能保证安全储藏。然而，粮堆是吸湿性多孔介质，仓外大气温度的季节性变化能使水分迁移带来安全隐患。仓储粮堆的水分迁移直接与粮堆内部的温湿度相关。一方面，由于仓外气温的季节性变化和昼夜变化，粮堆通过粮仓外围护结构与外界环境不断进行热量交换，使得其温度发生动态变化；另一方面，粮堆具有热惰性，使得粮堆内部容易产生温度梯度。在温度梯度的作用下，粮堆内部会形成一个与温度梯度方向相同的水蒸气压力梯度，导致水蒸气在粮堆内部的扩散，并通过粮堆的吸湿和解吸湿作用引起粮堆内部水分的迁移和再分配。同时，由于温度梯度而产生的粮堆内部空气的密度差，形成自然对流运动，它也会促进粮堆内部的水分迁移，如图1所示。由于粮堆内部的热量传递和水分迁移，可能会在粮堆局部形成高温高湿区域，不仅为害虫和霉菌生长提供了有利条件，还加剧粮食的呼吸作用而致使粮食品质劣变；另外，热量传递和水分迁移有时还可能会在粮堆局部形成低温高湿区域，一旦低于露点温度就会在粮仓内部出现结露现象，使得粮堆局部水分迅速增加，造成储粮发热、霉变、发芽的严重后果。

图1 仓储粮堆中热湿迁移的成因及其影响因素

Eckert等［1］通过试验证明，谷物具有一定的吸湿和解吸湿能力，在其内部的水分迁移过程中不会像沙堆那样出现零水分区域。Thorpe［2］采用一维数学模型描述了因温度梯度引起的水蒸气扩散，利用该模型的计算结果与Griffiths［3］的试验结果吻合很好。Tanaka等［4］假设谷物颗粒为球形，且谷粒表面水蒸气扩散与谷物水分和平衡水分之差成正比的关系，以水蒸气质量分率作为因变量对粮堆内部水分进行模拟分析。Nguyen［5］发展了一个二维的瞬态模型来描述由于粮堆内部的自然对流而引起的水分迁移过程，但忽视了扩散引起的水分迁移。Casada等［6］发展了一种有限差分模型预测带壳花生的热湿传递。借助薄层对流干燥速率方程模拟了花生中水分的变化率，但这种方法不能保持空气和花生之间的水分平衡。Smith等［7］建立了一个自然对流传热模型，通过平衡水分方程计算了水分为15%时粮堆内部空气的湿度，然后根据干燥速率方程间接得到粮堆的水分值。由于干燥速率只与局部温度有关，因此，其水分的计算违背基于基本守恒原理的传输机理。Khankari等［8］和 Jimenez－Islas等［9］对仓储粮堆内部由于温度梯度引起的温度和水分变化规律进行了研究，但他们忽略了由温度梯度引起的自然对流对储粮中热湿迁移的影响。Abe等［10］通过小型实验仓试验观测和数值模拟的方法考察了粮仓内谷物温度和水分在自然冷却时的变化规律，虽然模拟结果与试验观测结果基本相符，但在计算仓内谷物水分分布时，仅仅考虑了谷物的吸湿作用，忽略了热湿耦合关系对粮堆内部水分变化的影响。Iguaz等［11］和Marcelino等［12］根据多孔介质热质传递理论建立了扩散型的热湿耦合传递模型，模拟分析了仓储小麦和稻谷的温度和水分随外界气温变化而变化规律，模拟结果与试验测量结果基本相符。Ali等［13］和Alabadan等［14］采用二维柱坐标的导热微分方程模拟了圆筒仓内谷物温度随仓外气温变化而变化的规律，但研究没有涉及谷物水分的变化。

本研究是基于谷物吸湿／解吸湿和多孔介质传热传质学的理论基础，依据局部热质平衡原理和多孔介质宏观体积平均法，建立具有吸湿性作用的深层仓储粮堆内部热湿耦合传递过程的数学模型，并对建立的数学模型进行验证。利用验证后的模型对仓储粮堆内部空气自然对流流动、热量传递和水分迁移进行模拟分析，考察近似冬夏两个季节仓储粮堆中温度场和水分场的动态变化规律，从而为制定预防和消除储粮结露、发热霉变事故的措施提供理论依据。

1 数学模型及其验证

1.1 数学模型的建立

假设多孔介质内部的空气是不可压缩的，且在温度梯度作用下而产生一定的浮升力。同时，假设粮堆中谷物颗粒是各向同性的均质多孔介质，并且认为谷物颗粒与周围空气之间是局部热质平衡的。

式中：uj（j＝1，2）为xj方向空气的流动速度，u1＝u，u2＝v。xj（j＝1，2）为系统坐标，在笛卡尔坐标系统，x1＝x，x2＝y，x为水平方向，y为垂直方向。

动量方程：粮堆内部空气流动采用扩展的达西定律，即Brinkman－Darcy方程，假设浮升力随温度呈线性变化，满足Boussinesq近似。

式中：ρa为空气的密度，ρa取 1.225 kg·m－3，ui（i＝1，2）为空气的流动速度，t为时间，p为空气的压强，ε和K为多孔介质的孔隙率和渗透率，对于粮堆（小麦）来说，ε＝0.4，K＝5.96×10－9m2［8］。T为绝对温度，ρ0为温度为T0时空气的密度，T0为参考温度。g是重力加速度，β为空气的体积膨胀系数，μ为空气的动力黏度，μ＝1.79×10－5Pa·s。

能量方程：由于谷物颗粒与周围空气之间是局部热平衡的，谷物颗粒和周围空气中的化学势皆相等［8］，于是：

式中：ca空气的比热容，ρb、cb和kb为粮堆的密度、比热容和有效导热系数，对于粮堆（小麦）来说，ρb＝650 kg·m－3。W为粮堆干基水分，hs为吸湿性多孔介质的吸湿或解吸湿热，为了计算方便可以近似取为2 476 550（J／kg）［15］。对于小麦，cb＝（1．398＋0．0409×M）J·kg－1·K，kb＝（0．1170＋0．0013×M）W·m－1·K－1。M是粮堆的湿基水分含量，M＝W／（1＋W）×100［16］。

水分迁移方程：根据组分质量守恒原理，即谷粒含湿量（干基水分）以及谷粒间空气中的绝对含湿量在对流扩散过程中满足守恒定律，即

式中：τ为粮堆中孔隙的迂曲度，τ取1.53［2］。w为谷物颗粒间空气中绝对含湿量。Dv为谷物颗粒间空气中水蒸气的扩散系数，对方程（4）进行数量级分析知，方程（4）左边第一项中的可以忽略不计。另假设空气－水蒸气混合物为理想气体，则 ρaw＝pv／RvT，pv为谷粒间空气的水蒸气分压，Rv为水蒸气的气体常数，大小为461.52J／（kg·K）。令上述方程（4）变成：

由于局部热平衡原理，谷粒中水蒸气分压力和周围空气的水蒸气分压力是相同的，根据修正的Henderson方程［17］，谷物颗粒吸附等温线关系有如下形式：

式中：KH＝1.2299×10－5，C＝64.346，N＝2.5558（关于小麦的常数）［8］，ps是饱和蒸汽压，ps＝

根据方程（6），水蒸气分压力梯度能表示成谷物颗粒水分和温度的梯度形式，关系式如下：

设定，方程（5）改写为：

式中：DM＝Deffσ，DT＝Deffω，DM为水分梯度引起的水分扩散系数，DT为温度梯度引起的水分扩散系数，σ为一定温度时由于谷物水分引起的水蒸气分压，ω为一定水分时由于谷物温度引起的水蒸气分压。方程（7）的派生值σ与ω可以根据方程（6）可以得到：

由于变量DM和DT是温度和湿度的因变量，方程（8）就变成了非线性的，并与能量方程（3）相互耦合。

1.2 数学模型验证

本文通过比较研究对前面建立的数学进行验证。Anderson等［18］通过试验方法研究了装满了初始湿基水分为14.6%的450 kg小麦的长方体形的试验仓内热量传递以及水分迁移规律，试验仓体积为0.57m3。试验仓顶部和底部为绝热，试验仓左右两侧壁面温度不同，从而形成横向方向的温度梯度，且仓壁是不渗透的。Khankari等［8］使用扩散模型对Anderson等［18］实工况进行了数值模拟比较。图2是本研究对Anderson试验工况的数值模拟结果与Anderson等［18］以及 Khankari等［8］研究结果的比较，从图2可以看出三者非常相符。

图2 模拟计算的水分比较

2 模型的模拟结果与分析

模型为一密闭的装满小麦圆筒仓，圆筒仓半径R为5 m，高度L为10 m，如图3。由于圆筒仓为轴对称空间结构，所以选择圆柱仓轴向横截面的一半为数值模拟几何对象，其网格划分如图4，径向为x轴方向，垂直向上的高度方向y轴方向。

图3 物理模型

图4 数值模拟计算区域

表1 模拟计算的4个工况的粮堆与大气的初温和初始湿基

2.1 粮堆内部自然对流流动

近似冬夏季工况下仓内空气自然对流流动如图5所示。从图5中可以看出，近似夏季工况时由于粮仓外部气温较高，粮堆初始温度较低，在粮堆内部形成温度梯度，由于温度梯度而产生的粮堆内部空气的密度差，导致粮堆内部形成自然对流运动，在粮堆中形成逆时针方向的空气自然对流流动；近似冬季工况时由于粮仓外部气温较低，粮堆初始温度较高，在粮堆中形成顺时针方向的空气自然对流流动。从图5还可以看出，粮堆右侧边界附近和右下角的自然对流相比较其他区域要强，粮堆内部自然对流流动的速度大小在10－5～10－4m／s数量级范围内，这也说明了粮堆内部空气的自然对流流动为层流。

图5 水分为14%和18%冬夏季工况150 d时粮堆内部自然对流流场

2.2 粮堆内部热量传递和温度分布

近似冬夏季工况下仓内温度分布如图6所示。夏季工况时，由于粮堆初始温度小于外界大气温度，而且粮堆具有热惰性，使得粮堆内部产生温度梯度。在温度梯度的作用下，造成粮堆温度从外到内逐渐升高，达到150 d时，粮堆右上半部的温度都超过了278 K（5℃）。近似冬季工况时，由于粮堆初始温度大于外界大气温度，使得粮堆内部也会产生温度梯度。在温度梯度的作用下，造成粮堆温度从外到内逐渐降低，达到150 d时，粮堆右上半部的温度降低到了276 K（3℃）。比较冬夏两季工况，可以发现由于冬季是粮堆内部形成了顺时针方向的空气自然对流，促进了热量的传递，导致顶部、壁面和右下角粮堆的温度下降得较快；而在夏季由于粮堆内部形成了逆时针方向的空气自然对流，使得热量不断向粮堆右上部传递，导致右上部粮堆的温度升高得较快。同时，可以看到无论是冬季还是夏季工况，在粮堆温度既有横向变化又有纵向变化。其原因是粮堆底部边界为绝热，顶部和右侧面为定常温度的缘故。

图6 水分为14%和18%冬夏工况150 d时粮堆内部温度场

2.3 粮堆内部水分迁移和水分分布

近似冬夏两季工况下仓内粮堆水分分布如图7所示。近似冬季工况时，水分从底部和右下角向粮堆的顶部和右侧面迁移，使得粮堆上部和右侧水分升高，粮堆下部水分降低。其中右下角水分最低，顶部水分最高，形成所谓“结顶”现象。近似夏季工况时，水分从粮堆的顶部、右侧面和右下角向粮堆里部迁移，使得粮堆上半部（x＜4.5 m，8.0 m＜y＜9.5 m）水分升高，顶部、右侧和右下角粮堆水分降低，其中，右下角粮堆水分最低，右上半部（x＜4.5 m，8.0 m＜y＜9.5 m）水分最高。这与 Hellevang等［19］现场调查的水分结果基本一致。比较冬夏两季的水分分布图可以看出，在相同的初始水分以及相同的温差情况下，夏季达到150 d后粮堆内部水分变化幅度比冬季工况要大，但冬季升高幅度（相对于初始水分）要大于夏季。

图7 初始干基水分为0.163（14%）和0.22（18%）时近似冬夏工况时粮堆内部水分场

3 结论

本研究基于有限元的方法模拟分析了外界气温和小麦分别为273 K和293 K、293 K和293 K，小麦初始水分为14%和18%时直径为10 m高度为10 m的充满小麦的圆柱仓内部的自然对流流动和热量传递及水分迁移过程，重点探究了近似冬季和夏季仓储粮堆内部温度和水分的动态变化规律。

3.1 冬夏两季都会在的粮堆内部形成空气自然对流运动，夏季时在粮堆中形成逆时针方向的自然对流流动，冬季时形成顺时针方向的自然对流流动。冬夏两季粮堆右侧边界和右下角的自然对流相比较其他区域要强，粮堆内部自然对流速度在10－5～10－4m／s数量级范围内。

3.2 冬季工况下水分从粮堆内部向顶部和右侧壁面迁移，粮堆顶部和右侧水分升高，粮堆右下角水分降低，顶部水分最高，形成所谓“结顶”现象。夏季工况时，水分从粮堆的顶部、右侧面和右下角向右上半部区域（x＜4.5 m，8.0 m＜y＜9.5 m）迁移，使得顶部、右侧附近和右下角粮堆水分降低。比较冬夏两季的水分分布图可以看出，在相同的初始水分以及相同的温差情况下，夏季达到150 d后粮堆内部水分变化幅度比冬季工况要大，但冬季升高幅度（相对于初始水分）要大于夏季。

［1］Eckert ER G，FaghriM A.A general analysis ofmoisturemigration caused by temperature differences in an unsaturated porousmedium［J］.International Journal of Heat and Mass Transfer，1980，23：1613－1623

［2］Thorpe G R.Moisture diffusion through bulk grain subjected to temperature gradient［J］.Journal of Stored Products Research，1982，18：9－12

［3］Griffiths H J.Bulk storage of grain：A summary of factors governing control of deterioration［G］.Report ED8 of CISRO Division of Mechanical Engineering，Melbourne，Australia，1964

［4］Tanaka H，Yoshida K.Heat and mass transfermechanisms in a grain storage silo［M］／／McKenna B M.Engineering Sciences in the Food Industry：Engineering and Food.London，UK：Elseviser.Science Pubs.，1984：89－98

［5］Nguyen T V.Modeling temperature and moisture changes resulting from natural convection in grain storages［M］／／Champ B R and Highely E H.Preserving Grain Quality by Aeration and In－Store Drying.Canberra，AUS：ACIAR Proc.，1986：81－87

［6］Casada M E，Young JE.Heat and mass transfer during shipment of peanuts［M］／／ASAE.Michigan，USA：Joseph St.，1989：89－6105

［7］Smith E A，Sokhansanj S.Moisture transport caused by natural convection in grain stores［J］.Journal of Agricultural Engineering Research，1990，47：23－34

［8］Khankari K K，Morey R V，Patankar S.V.Mathematicalmodel formoisture diffusion in stored grain due to temperature gradients［J］.Transactions of the ASAE，1994，37（5）：1591－1604

［9］Jiménez－Islas H，Navarrete－Bolaños J，Betello－Alvarez E.Numerical study natural convection heat and mass 2D cylindrical grain stored in silos［J］.Agrociencia，2004，38：325－342

［10］Abe T，Basunia M A.Simulation of Temperature and Moisture Changes during Storage of Rough Rice in Cylindrical Bins Owing to Weather Variability［J］.Journal of Agricultural Engineering Research，1996，65：223－233

［11］Iguaz A，Arroqui C，Esnoz A，Virseda P.Modeling and Validation of Heat Transfer in Stored Rough Rice without Aeration［J］.Biosystems Engineering，2004，88（4）：429－439

［12］Marcelino Carrera－Rodríguez，Gloria M，Martínez－González，et al.Transient numerical study of the effect of ambient temperature on 2－D cereal grain storage in cylindrical silos［J］.Journal of Stored Products Research，2011，47：106－122

［13］AliM SAl－Amri，Sirelkhatim K Abbouda.Application of a Mass Transfer Model for Simulation and Prediction of Moisture Distribution in Stored Corn Grains［J］.Scientific Journal of King Faisal University（Basic and Applied Sciences），2004（5）：197－213

［14］Alabadan B A.Temperature Changes in Bulk Stored Maize［J］.Assumption University Journal of Technology，2006，9（3）：187－192

［15］Brooker D E，Arkema FW Bakker，Hall CW.Drying Cereal Grains［M］／／Westport，Conn.：The AVI Publishing Co.Inc.，1974

［16］ASAE Standards，37th Ed.Thermal properties of grain and grain products［S］.Michigan：Joseph St.，1990

［17］Chen C C，Morey R V.Comparison of four EMC／ERH equations［J］.Transactions of the ASAE，1989，32（3）：983－990

［18］Anderson JA，Babbitt JD，Meredith W O S.The effect of temperature differential on the moisture content of stored wheat［J］.Canadian Journal of Research，1943，21（C）：297－306

［19］Hellevang K J，Himing H J.Moisture movement in stored grain during summer［M］／／ASAE.Michigan，USA：Joseph St.，1988：88－6052.

Numerical Study on Natural Convection Heat and Moisture Transfer in Stored Grain in A Silo

Wang Yuancheng1Bai Zhongquan2Zhang Zhongtao1QiWei1
（College of Thermal Energy Engineering，Shandong Jianzhu University1，Jinan 250101）（Henan Design＆Research Institute CO.，Ltd of National Defence2，Zhengzhou 450003）

In the paper，a mathematical model for the simultaneously natural convective heat and moisture transfer had been developed by utilizing an incorporating a sorption isotherm relationship to predictmoisturemigration due to natural convection and diffusion within stored grain.Themodelwas solved numerically by the Finite Element Method（FEM）to predict the natural flow patterns，temperature and grain moisture distributions on conditions approximating winter and summer storage periods in a silo with 10 m in diameter a height of16 m.The conditions approximating winter and summer extremeswere independently simulated for a period of150 d in four caseswith two initialmoisture contents as 14%and 18% （w.b.）respectively.Two sets of constant wall and grain temperatures，0℃and 20℃，20℃and 0℃were studied respectively.Natural convection currents had been developed in silos due to temperature gradients in approximating winter storage periods；while natural convection were counterclockwise in approximating summer storage periods.On approximate winter conditions，themoisture was concentrated in relatively narrow bands，near the lateral edge；a small zone near the surface of bulk wheat，and maximum moisture contents occurred in upper left－hand corner in bulk wheat.On summer conditions，moisture was concentrated in the upper zones（x＜4.5 m，8.0 m ＜y＜9.5 m）of the bulk wheat，which waswider and deeper than that created on the approximate winter conditions.

coupled heat andmoisture transfer，natural convection，grain storage，numerical simulation，FEM

TK124

A

1003－0174（2015）11－0097－06

国家自然科学基金（51276102），国家粮食公益专项（201313001），山东省自然科学基金（ZR2011EEM 011）

2014－05－05

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