基于COMSOL的平房仓冷却通风过程中粮堆热湿耦合传递研究

陈桂香 张宏伟 王海涛 刘超赛 尹 君

(河南工业大学土木建筑学院1, 郑州 450001)(国家粮食局科学研究院2, 北京 100037)

我国每年粮食总产量高达5亿t，仓内害虫、粮堆发热和谷物霉变是引起储藏过程中粮食损失的主要原因[1-2]。粮堆内的温度和水分是影响谷物霉变、粮堆发热和仓内害虫的重要因素[3]。粮仓冷却干燥通风可以把粮堆温度降低至安全温度以下，把粮堆水分减小到安全水分以下，是一种有效减少粮食在储存过程中损失的技术措施。研究冷却干燥通风过程中粮堆内热湿耦合传递规律，对指导粮仓的通风冷却和改善粮食储存环境具有重要意义。

在冷却干燥通风过程中，粮堆内部热湿迁移是一个复杂的多场耦合过程[4]，粮堆内部热湿迁移过程的影响因素很多，包括粮仓所在地气象条件、粮堆物性参数和粮食生物特性等。冷却干燥通风过程中粮堆热湿耦合迁移规律研究涉及湍流、湿空气传热、多孔介质传热、多孔介质传质和环境大气等多个物理场。

数值模拟作为一种经济有效的研究手段，在探索粮堆内部的热湿传递规律过程中日益受到学者的重视和关注[5]。Chang等[6-7]提出了一个预测小麦通风状态下温度和水分变化的数学模型，并据此对粮堆热湿传递规律进行了数值模拟研究。尹君等[8-9]利用MATLAB模拟软件重现了不同仓型的小麦粮堆温度场和水气分压场的变化过程，并研究了不同季节粮堆的“冷芯”“热芯”和结露现象。陈桂香等[10]采用CFD软件数值模拟研究了通风过程中粮堆内热湿传递及霉变预测。李祥利等[11]通过数值模拟方法研究了通风状态下“圭”字形风道的高大平房仓粮堆内部温度以及水分的变化规律。王雪等[12]模拟研究了粮食温度场在仓内部热源作用下的分布及变化情况。张前等[13]采用回归模型方法预测了高大平房仓内部粮堆的温度变化。Daniela de Carvalho Lopes等[14]采用数值模拟软件AERO研究了粮堆内发热点对粮堆一维传热的影响。Carrera-Rodríguez等[15]采用Fortran软件研究了周围环境温度对高粱粮堆空隙内热环境的影响。

本研究建立了冷却干燥通风过程中高大平房仓COMSOL三维物理模型，以实测的送风空气温度和湿度为入口边界条件，通过修改COMSOL内置材料方程、耦合渗流、能量守恒、动量守恒和水分迁移控制方程，数值模拟粮堆温度场、湿空气速度场、湿空气压力场和水分变化过程，得出了冷却干燥通风过程中粮堆内热湿耦合传递规律。

1 粮堆内热湿耦合传递的COMSOL数值模拟

1.1 高大平房仓三维物理模型

以位于华东地区长宽高为42 m×24 m×11 m的高大平房仓为研究对象。该粮仓采用六个“一机两道”的压入式通风系统，采用“U”字型地笼通风道。粮仓装备了一套无线传感器的粮情监控系统。粮仓布置了83根测温电缆，共安装了四层332个温度传感器。高大平房仓进行冷却干燥通风时，粮仓内小麦堆粮高度为8 m。粮仓冷却干燥通风时间为70 h。在COMSOL软件中，建立了如图1所示的高大平房仓三维物理模型。整个粮仓物理模型采用非结构网格形式共划分为350多万个网格。

图1 高大平房仓的COMSOL三维物理模型

1.2 COMSOL软件开发接口与参数设置

COMSOL 为用户提供了多种开发接口方式，其主要开发接口方式有：利用MATLAB 进行二次开发、利用应用程序接口API 进行二次开发、利用物理模型创建器进行二次开发和利用修改COMSOL 内置方程进行二次开发。用户可以利用MATLAB 的数据处理功能、利用COMSOL 提供的API 函数、利用COMSOL 提供的创建器开发出自定义物理场模型或修改COMSOL 内置方程，根据研究需要进行COMSOL 软件二次开发。在本研究中，COMSOL 数值模拟涉及的粮食物性参数和边界条件可为任意的常数、变量或函数表达式等多种形式，因此本文主要采用两种COMSOL 软件开发接口形式：通过修改COMSOL 内置方程进行二次开发和利用应用程序接口API 进行二次开发。

以实测的送风湿空气的温湿度拟合函数作为模拟研究的入口边界条件。参数的正确设置对粮堆内热湿耦合传递COMSOL 数值模拟十分重要。表1给出了COMSOL数值模拟时的相关条件设置和参数设置。

表1 COMSOL数值模拟的相关条件和参数设置

1.3 粮堆内热湿耦合传递COMSOL数值模拟控制方程

冷却干燥通风过程中粮堆内热湿耦合传递COMSOL数值模拟需要湿空气和多孔介质两种材料，利用COMSOL内置材料方程将小麦粮堆物理参数赋值给多孔介质材料。COMSOL数值模拟涉及湍流、湿空气传热、湿空气传质、多孔介质传热和多孔介质传质等多个物理场，还需要热湿耦合、温度耦合和流动耦合多个耦合物理场。多孔介质热湿耦合传递数值模拟需要描述谷粒吸附解析水分时的热量交换和水分迁移，该过程的实现需将COMSOL 内置方程修改为粮堆内热湿耦合传递控制方程。基于非稳态的多孔介质热湿传递偏微分方程数值解，以含湿量、温度和压力为驱动势，建立了冷却干燥通风过程中粮堆内热湿耦合传递动态数学模型。

1.3.1 粮堆内水分迁移控制方程

COMSOL数值模拟将粮堆视为各项同性的连续多孔介质。根据单元体质量守恒，式(1)给出了粮堆内水分迁移控制方程，用于描述粮堆内热湿耦合传递的水分迁移情况。

(1)

式中：ω为体积含湿量/kg/m3；t为时间/s；jv为水蒸气传递速率/kg/(m2·s)；jl为液态水传递速率/kg/(m2·s)；Sw为水分源项[16]，Sw利用Thorpe[17]介绍的方法计算确定。

根据菲克定律和达西定律：

jv=-δpPv+jaxa

(2)

jl=KlPk

(3)

式中：δp为水蒸气渗透率/kg/(m·s·Pa)；Pv为水蒸气分压力/Pa；ja为空气流动速率/kg/(m2·s)；xa为粮食颗粒间湿空气含湿量/kg/kg；Kl为液态水渗透率/kg/(m·s·Pa);Pk为毛细水压力/Pa。

将式(2)和式(3)代入式(1) 得：

(4)

假设温度对粮堆平衡含湿量的影响忽略不计，则有：

(5)

式中：ζ是等温吸放湿曲线的斜率/kg/m3；φ是粮食颗粒间湿空气相对湿度，φ的开尔文关系式可以表示为：

(6)

粮食颗粒间湿空气按理想气体处理，则有：

Pv=φPs

(7)

xa=6.2×10-6Pv

(8)

式中：ρl为水的密度/kg/m3；RD为水蒸气气体常数/J/(kg·K)；T为开尔文温度/K；Ps为饱和水蒸气压力(Pa)。

将式(5)～式(8)代入式(4)，可以得到：

(9)

式中：Dφ和DT可以分别利用式(10)和式(11)计算确定。

(10)

(11)

1.3.2 粮堆内热量平衡控制方程

COMSOL数值模拟假设粮堆内的水分只有汽、液两相[18]，根据单元体能量守恒，式(12)给出了对流传热过程热量守恒偏微分方程。

(12)

式中：ρm是小麦颗粒的干基密度/kg/m3；ρm=790.3 kg/m3。Cp,m是粮堆的恒压热容/J/(kg·K)，Cp,m=1 934.2 (J/(kg·K))；qcond是导热热流密度/W/m2；qconv是对流热流密度/W/m2；Sh是热量源项[19]，Sh利用Thorpe[17]介绍的方法计算确定。

qcond=-kT

(13)

qconv=Cp,ajaT+jvhlv

(14)

式中：k是粮堆的导热系数/W/(m·K)；Cp,a是干空气的恒压热容/J/(kg·K)；hlv是水蒸气汽化潜热/J/kg。

将式(13)和式(14)代入式(12)得：

(15)

1.3.3 粮堆内动量守恒控制方程

COMSOL数值模拟根据泊肃叶定律[20]，将粮堆内的空气平均流动速度表示为压力梯度与速度的关系式。

ja=-kaPa

(16)

式中：ka是粮堆内湿空气的渗透率/m2；ka是指在空气流动方向上空气流动速率与压力梯度的比值，其物理意义是指在一定压差下粮堆允许湿空气通过的能力；Pa是湿空气压力/Pa。

根据连续性方程，式(16)可以变形为

(17)

式中：ρa是湿空气密度/kg/m3；ε是粮堆孔隙率，ε=0.4。

在粮堆内热湿耦合传递过程的COMSOL数值模拟研究中，粮堆内的湿空气流速低，湿空气流动的压力低，湿空气流动过程中温度变化小，因此粮堆内湿空气可以作为不可压缩气体处理，则式(17)可简化为：

(18)

1.4 计算机软件和硬件配置

本研究采用AutoCAD软件绘制高大平房仓三维几何模型，然后把高大平房仓三维几何模型导入COMSOL Multiphysics 5.3软件，构建形成联合体的高大平房仓COMSOL三维物理模型。采用COMSOL Multiphysics 5.3软件模拟计算冷却干燥通风过程中粮堆内热湿耦合传递过程，COMSOL Multiphysics 5.3和Origin8.5软件作为数值模拟的绘图及数据分析软件。

2 模拟结果与讨论

数值模拟了冷却干燥通风38 h粮堆内热湿耦合传递过程，分析了粮堆温度场、湿空气速度场和湿空气压力场的变化情况，通过对比预测值与实测值，验证了COMSOL数值模拟结果的正确性，进而得出了冷却干燥通风过程中粮堆内热湿耦合传递的规律。

2.1 数值模拟结果的准确性验证

图2为冷却干燥通风过程中的粮堆平均温度模拟预测值和实测值的变化关系。分析图2可知，冷却干燥通风过程中粮堆平均温度模拟预测值与实测值吻合较好，两者的最大温度偏差为0.211 ℃，表明COMSOL数值模拟的温度预测精度较高。

图2 粮堆实测和模拟预测温度平均值

图3 粮堆4 m高度处实测温度平均值与预测温度平均值

利用粮堆4 m高度处温度传感器测量值验证COMSOL数值模拟结果的准确性。图3为粮堆4 m高度处实测温度与预测温度平均值。由图3可知，冷却干燥通风开始阶段粮堆4 m高度处实测温度平均值与预测温度平均值之间有较大的温差，随着冷却通风时间的增加实测温度平均值与预测温度平均值之间温差逐渐减小，冷却通风结束时粮堆4 m高度处预测温度平均值已十分接近实测温度平均值。这是因为数值模拟假定粮堆初始温度为16.7 ℃的等温体，导致冷却通风开始时粮堆4 m高度处预测温度平均值与实测温度平均值之间有较大的温差。随着冷却通风时间的增加，粮堆4 m高度处预测温度平均值逐渐接近实测温度平均值。这进一步证明COMSOL数值模拟具有较高的预测精度。

2.2 数值模拟结果分析与讨论

利用建立的高大平房仓COMSOL三维物理模型，通过改变初始条件，模拟了粮堆高度相同的新粮粮堆冷却干燥通风过程。新粮粮堆的初始温度为32 ℃，粮堆的初始干基水分为16.7%。冷却干燥通风的数值模拟时间为70 h。数值模拟研究冷却干燥通风过程中粮堆热湿耦合传递过程，给出粮堆温度场、湿空气压力场和湿空气速度场变化情况。

图4 通风70 h粮堆截面y=6 m的湿空气压力场云图

图5 通风70 h粮堆截面y=6 m的速度场云图

图4是冷却干燥通风过程中粮堆湿空气压力场云图。由图4 可知, 通风地笼内的湿空气压力最大，湿空气压力从粮堆底部往上逐渐减小。粮层阻力是造成这一现象的主要原因。湿空气压力分布存在明显的分层现象，而且基本上是对称分布的。

图5 是通风过程中粮堆湿空气速度场云图。由图5 可知，通风地笼内的湿空气速度最大，粮堆中湿空气速度最小，粮堆上部湿空气速度介于两者之间，这一现象出现的原因也主要是粮层阻力的存在。在冷却干燥通风过程中，粮堆温度场变化会受湿空气速度的影响，而湿空气速度分布又取决于其压力的分布。

图6是通风过程中粮堆温度场云图。分析图6可知，粮堆温度场分布存在明显的分层现象，而且基本上是对称分布的。与粮堆的初始温度相比，温度下降较为明显。粮堆平均温度由初始平均温度32 ℃降为14.74 ℃。粮堆底部的送风道周围区域温度最低，自粮堆底部向上粮堆温度逐渐升高，未被冷却的粮层温度最高。相邻两个风道之间粮堆温度下降速度相对较慢，粮堆底部两个风道之间的中心区域比其周围区域的温度高，主要原因是该温度较高区域处于通风死角，通过该区域的冷却干燥通风湿空气量相对较少，该区域的体积大小跟粮层厚度、送风风速、送风道间距和送风道类型等因素有关，可以通过粮仓通风系统优化来减小通风死角区域的体积。

图6 通风70 h粮堆截面y=6 m的温度场云图

图7 冷却干燥通风过程中粮堆干基水分的预测值

图7给出了冷却干燥通风过程中粮堆干基水分变化图。分析图7可知，粮堆干基水分由初始水分16.7%降为14.5%，表明冷却干燥通风可以干燥小麦粮堆。在冷却干燥通风前期阶段水分下降速度较快，随着通风时间的增加粮堆干基水分下降速度不断减小，分析水分迁移控制方程可知，小麦颗粒的水分与颗粒空隙间湿空气的相对湿度和温度之间存在对应关系。随着通风时间的增加，粮堆干基水分和温度不断降低，而粮食颗粒内水蒸气分压力与湿空气水蒸气分压力之间的压力差不断减少，因此粮堆干基水分的下降速度不断减小，直至小麦粮堆水分与湿空气水分达到近似平衡水分状态。小麦粮堆达到平衡水分状态后，继续进行冷却干燥通风，粮堆干基水分将不再发生明显变化。

3 结论

3.1 通过多场耦合有限元软件COMSOL 的二次开发接口，可以实现冷却干燥通风过程中粮堆内热湿耦合传递的数值模拟。通过修改COMSOL 材料内置方程，使多孔介质材料具有了小麦粮堆的物性参数。通过将COMSOL 内置方程修改为粮堆内水分迁移、热量守恒和动量守恒控制方程，可以实现冷却干燥通风过程中粮堆内热湿耦合传递过程的数值模拟。粮堆孔隙率、渗透系数和密度等物理参数会影响数值模拟的预测精度，数值模拟时应该根据实际情况合理设置这些参数。

3.2 在粮仓冷却干燥通风过程中，粮堆降温和降水是同时存在的。粮堆干基水分的变化方向是逐渐趋于粮食的平衡水分状态，粮堆干基水分的变化方向和变化速率主要跟粮堆干基水分值和通风的湿空气参数有关。通风正上方区域粮堆降温较为明显，相邻通风道之间的区域粮堆降温效果较差，冷却干燥通风过程中存在通风死角区域，通过粮仓通风系统优化可以消除或减小通风死角区域的体积，可以更好地保障储粮安全。