考虑压力和含水率影响的小麦粮堆导热特性数值研究

蒋敏敏 曹会敏 王清山，3 柴威丽

(1.河南省粮油仓储建筑与安全重点实验室，河南 郑州 450001；2.河南工业大学土木工程学院，河南 郑州 450001；3.中国建筑第七工程局有限公司工程研究院，河南 郑州 450004)

0 引言

2022 年，我国粮食全年产量已突破6.8 亿t［1］，粮食储备安全与国家安全、国计民生密切相关。在粮食储存过程中，受环境温度、仓内压力、孔隙率等因素的影响，易出现粮食霉变、生虫等问题，导致我国粮食储藏的平均损耗率达8%，其中，国家粮食储备损失约占全国粮食储备总量的0.2%［2］，农民粮食储备损失占农户粮食储备总量的6%～9%。粮堆是由粮食颗粒和颗粒间的孔隙组成的二相体，由于粮食颗粒与空气的导热系数相差数十倍，且在压力等因素的作用下，粮堆的孔隙率会发生较大变化，导致由粮食颗粒与空气组成的粮堆等效导热系数发生较大变化。因此，对不同条件下粮堆的导热特性进行研究具有重要的学术与应用价值。

目前，国内外学者对颗粒材料的导热系数进行了大量研究。Kayisoglu 等［3］通过试验验证了含水率对大米导热系数的影响；程龙等［4］使用球形测量技术来研究小麦、水稻、玉米、大豆的导热系数随水分迁移、稳态时间、籽粒加载均匀性的变化情况；周祖锷等［5］在不同含水率和堆积密度下，对小麦、玉米、水稻等8 种散体物料的导热系数、导温系数和比热进行测定，总结出粮食的比热、导温系数随含水率的变化规律；杨洲等［6］、鞠卫平等［7］、金文等［8］对稻谷导热系数、热扩散系数、比热与含水率的关系进行了研究，即导热特性参数与稻谷含水率之间呈线性关系；Huetter 等［9］通过热线法试验得到了不同压力和孔隙率条件下玻璃球堆积体的等效导热系数；Sun［10］研究了砂岩的热导率特性，砂岩导热系数受矿物组成、微观结构、孔隙度、饱和度、压力和温度的影响，即导热系数随孔隙水饱和度和压力的增加而增加，随孔隙率和温度的增加而减少；Suleiman 等［11］研究了奇亚籽、小麦和金岩的导热系数与水分含量之间的关系。

目前，已有大量学者对各类散体材料导热特性进行研究［12-14］，在粮堆导热特性研究方面，主要是在无上覆压力条件下进行的，不考虑粮堆压缩变形后对导热系数的影响，而部分考虑压力对导热特性的研究主要针对玻璃球、砂岩等材料。在现代大型粮仓内部，粮堆受到的压力较大，整体压缩变形和孔隙率变化也较大，导致由粮食颗粒和孔隙组成的粮堆整体等效导热特性发生较大变化，而目前相关研究仍不足。本研究通过压缩试验获得不同压力和含水率下小麦粮堆的孔隙率，利用数值法来建立粮堆等效导热特性的数值模型，从而分析压力、含水率等因素影响下的小麦粮堆导热特性，以期为粮堆通风冷却问题、降温时间控制、粮食干燥计算及干燥过程提供重要依据。本研究对提升我国储粮安全水平、实现绿色储粮等具有重要意义。

1 压力、含水率对粮堆孔隙率影响试验

1.1 试验材料与试验方案

选用河北产小麦作为试验材料，其初始含水率为11.67%，使用鼓风干燥器对其烘干。材料的形状用长、宽、高来度量，试验材料物理参数见表1［15］。将烘干后的小麦分别配制含水率为8.50%、11.21%、14.45%、16.97%、20.72%等5种状态，在压力等级为0 kPa、50 kPa、100 kPa、200 kPa、300 kPa等作用下进行压缩试验。粮堆的初始孔隙比见式（1）。

表1 压缩试验小麦试样的物理参数

孔隙率与孔隙比的关系见式（2）。

式中：Gs为散体颗粒的比重；ρw为水的密度，g∕cm3；ρ0为试样初始密度，g∕cm3；w0为试样的初始含水率。

在不同压力下，试样被不断压缩变形，压缩试验后试样孔隙率式（3）。

式中：ni为试样在i级压力下的孔隙率；h0为试样的初始高度，mm；n0为试样初始孔隙率；Δhi为粮堆高度变化，mm。

1.2 压缩试验结果分析

根据小麦粮堆压缩试验结果，得到不同压力和含水率下小麦粮堆的孔隙率，如图1 所示。绘制孔隙率—压力—含水率三维图，如图2 所示。当小麦颗粒含水率为8.50%～20.72%时，最松散状态时的孔隙率为0.41～0.484 7。随着压力的增大，粮食籽粒重新排布，粮食籽粒间的孔隙空间也会重排，从而导致孔隙减小。随着水分含量的增加，粮食籽粒膨胀，导致籽粒自身体积增大，从而导致粮堆系统的整体孔隙率增大。综上所述，压力和含水率增大会导致粮堆系统整体孔隙率增大。根据试验结果进行拟合，得出不同压力及含水率下的孔隙率的拟合公式，见式（4）。

图1 不同含水率和压力下小麦粮堆孔隙率

图2 不同压力及含水率下的孔隙率

2 粮堆温度压力耦合数值模型

2.1 数值模型的建立

由于粮堆内部孔隙结构的复杂性，各类孔隙模型均存在局限性［16-17］，本研究将具有复杂结构的粮堆等效为表征单元体（Representative Elementary Volume，REV），用来分析粮堆的导热特性［18-20］，即将粮堆中的颗粒等效为均一介质，并保持整个REV单元的孔隙率不变。

由于小麦粮堆中颗粒和孔隙均匀分布，因此在数值研究中，可将粮堆的三维导热特性问题简化为二维热传导模型［21］。然而将三维空间问题转变为二维平面问题后，其二维孔隙率必然与三维孔隙率不同，本研究采用何咏睿等［22］基于离散元颗粒流程序建立的二维、三维孔隙比关系，见式（5）［23］。

式中：e3d为三维孔隙比；e2d为二维孔隙比。利用孔隙率与压力、含水率关系及二维、三维孔隙比关系，可模拟出不同压力、含水率下二维热传导模型中等效REV单元的孔隙率。

在采用上述等效的孔隙率模拟法时，可运用COMSOL 有限元程序来模拟粮堆导热特性。在建模中模拟实际粮堆边长为L的单元体进行分析［24］，固相小麦颗粒则根据等效REV 方法，模拟半径为R的圆，则数值模型中粮堆孔隙率见式（6）。

式中：R为等效REV 单元中固相半径；L为整个REV 单元体（包括固相和孔隙）的边长。数值模拟的孔隙率与试验中的孔隙率保持一致，建立的粮堆二维导热特性等效REV 单元体模型如图3所示。

图3 粮堆二维导热特性等效REV单元体模型

2.2 材料参数和边界条件

根据相关研究［25-28］，数值模拟中小麦颗粒及空气的热传导参数见表2。为了便于计算出粮堆二维导热特性等效REV 单元体模型的等效导热系数［21］，将模型上、下表面设置为热绝缘边界，并在左右两个表面对REV 单元体施加温度边界条件T0和T1（T0＞T1），T0为40 ℃，T1为30 ℃。由于T0＞T1，故模型内的热流量会由左向右流动，模拟恒温边界条件下粮堆REV单元体的二维导热特性，如图4所示。

图4 粮堆二维导热特性边界条件

表2 数值模拟中小麦颗粒及空气的热传导参数

2.3 数值模拟方案

在COMSOL 有限元软件中建立粮堆二维导热特性等效REV 模型［29］。不同含水率和压力下粮堆孔隙率的变化情况可通过孔隙率与压力、含水率关系及二维、三维孔隙比来表示［30］，即含水率和压力的变化体现为孔隙率的变化［31］。研究含水率分别为8.50%、11.21%、14.45%、16.97%、20.72%，压力分别为0 kPa、500 kPa、1000 kPa、2000 kPa 和300 kPa下的粮堆导热特性［32］。

3 压力和含水率对导热特性的影响

3.1 数值计算结果

根据建立的粮堆等效REV 导热数值模型，得到粮堆的REV 单元温度分布如图5 所示。由图5可知，在粮堆REV 单元两侧设置恒温边界后，由于中心粮堆等效单元体的导热系数较大，而边缘处空气单元体的导热系数小，因此，更多的热量会通过粮堆进行等效单元体传递。根据粮堆的传热原理，粮堆REV单元的等效导热系数见式（7）。

图5 粮堆REV单元温度分布

3.2 压力对导热特性的影响

由粮堆REV 单元数值计算结果得到不同压力下小麦粮堆等效导热系数，如图6 所示。由图6 可知，粮堆等效导热系数随上覆压力的增大而增大，即粮堆下部（上覆压力大的区域）的导热系数比粮堆表层（上覆压力为零）的要大，含水率为8.5%～20.72%、上覆压力为300 kPa 时，导热系数增大17.5%～27.5%。通过线性关系来拟合粮堆REV 单元数值计算结果［35］，可得到粮堆等效导热系数与上覆压力的关系，见式（8）。

图6 压力对小麦粮堆等效导热系数的影响

式中：c为拟合线的斜率，表示导热系数随着压力增大的增长率，W∕(m·K·kPa)；keff0为拟合线的截距，表示压力为零的导热系数。

参数c随含水率变化而变化的曲线如图7 所示。由图7 可知，参数c取值为0.663×10-4～1.051×10-4W∕(m·K·kPa)，含水率小于14.45%，参数c随含水率的增大而增大；含水率大于14.45%，参数c随含水率的增大而降低。参数keff0随含水率的变化曲线如图8 所示，参数keff0取值为0.107 9～0.137 8 W∕(m·K)，参数keff0随含水率的增大而略有降低。

图7 参数c的变化

图8 参数keff0的变化

4 结论

本研究先通过压缩试验得到不同压力及含水率下小麦粮堆的孔隙率，再通过离散元与有限元相结合的方法，建立小麦粮堆导热特性计算的REV单元模型，分析了不同压力和含水率影响下的小麦粮堆导热特性，得出以下结论。

①通过压缩试验得到小麦粮堆压力、含水率与孔隙率之间的关系，可得出粮堆孔隙率与压力、含水率之间的拟合公式。

②基于三维、二维孔隙率之间的转化关系及粮堆孔隙率与压力、含水率之间关系，运用离散元与有限元相结合的方法，建立小麦粮堆等效REV 单元数值模型。

③粮堆等效导热系数随着上覆压力的增大而增大，当含水率为8.5%～20.72%、上覆压力为300 kPa时，导热系数增大17.5%～27.5%。