仓储粮堆机械通风时压力场的模拟研究

任广跃 彭 威 张忠杰 吴子丹 王 芳

(河南科技大学食品与生物工程学院1，洛阳 471003)

(国家粮食局科学研究院2，北京 100037)

仓储粮堆机械通风时压力场的模拟研究

任广跃1彭 威1张忠杰2吴子丹2王 芳1

(河南科技大学食品与生物工程学院1，洛阳 471003)

(国家粮食局科学研究院2，北京 100037)

以仓储粮堆机械通风试验为基础，以实仓试验相关数据为依据，建立了仓储粮堆机械通风过程中内部压力场分布的计算流体动力学(CFD)模型，并进行了不同通风风量条件下开环流熏蒸口和开人孔时的数值模拟计算。最后通过实仓试验数据和模拟结果进行验证。结果表明，仓储粮堆机械通风过程CFD模型能够真实反映机械通风过程中粮堆内部压力场的分布情况，将计算流体力学模拟技术应用于仓储粮堆压力场的预测分析是可行的。

仓储粮堆 机械通风 压力场 计算流体动力学

粮食是国民经济重要的战略物质，粮食安全储藏关系到国家安全稳定、国民经济发展以及可持续发展等重大战略问题，粮食生产的季节性与消费的常年性矛盾长期存在，加之粮食是有生命的活体［1］，因此如何实现粮食的安全储存成为世界性难题。我国粮食产量大，储备量高，同时由于储藏时间长，导致储备期间的损耗也较高。据调查，一个储藏周期内，仅因储粮损耗减量就高达1%～2%，如以我国粮食储备动态数量1.5亿吨来计算，储粮损失约为150～300万吨，折合人民币23～45 亿元，造成巨大的经济损失［2－3］。

目前，国内粮库在储粮过程中广泛采用了机械通风技术［4］，作为智能粮情检测、低剂量环流熏蒸、智能机械通风和高效谷物冷却“四合一”储粮新技术之一［5］，其在降温降水、降温干燥和增湿调质等方面发挥越来越重要的作用。机械通风技术主要是利用通风机产生的压力，将外界低温、低湿的空气送入粮堆，促使粮堆内外气体进行湿热交换，降低粮堆内的温度与水分，增加储粮稳定性的一种安全储粮技术［6］。通过大量实践表明，机械通风在粮食储藏中的应用，使粮食整年处于低温储粮状态，粮食的品质明显提高，延缓陈化速度，抑制新陈代谢，降低储粮损耗，有效地控制储粮害虫和微生物的活动、危害，它代表着未来粮食安全储藏的方向。但在机械通风时仓储粮堆内部流场的分布情况比较复杂，很难进行有效地测量和计算，而计算流体动力学(CFD)方法不但能突破一些现实条件的限制，而且通过计算模拟达到减少操作成本的目的，基于此，CFD方法可以有效的应用于粮食储藏中内部流场模拟研究［7］，如张忠杰等［8］利用FLUENT软件对平房仓的机械通风情况进行了数值模拟，获得了仓储粮堆的温度场分布及其变化规律;王远成等［9］采用CFD方法对就仓通风时粮堆内部热湿耦合传递过程的数值预测，建立了通风储粮过程中控制粮堆内部热量和水分传递的数学模型，并对就仓通风时粮堆内部热湿耦合传递过程进行了预测，获得了就仓通风时粮堆内部热量和水分迁移的基本规律。

粮食仓储内部生态系统决定储粮的安全，仓储粮堆内部压力场的CFD模拟研究，目的就是应用CFD模拟，采用对比分析实仓试验的方法，对仓储粮堆内部压力场分布变化情况进行验证，掌控粮食在机械通风过程中其压力场的分布情况，以期为储粮安全提供新的检测方法和手段。

1 材料与方法

1.1 试验材料与设备

试验粮仓位于北京市昌平区国家粮食局科学研究院小汤山中试基地内，该仓为室外180 t浅圆钢板仓，仓高为10 m，其中圆柱形仓高8.5 m，锥形顶部高1.5 m，仓半径为3 m，仓内粮食堆高7.5 m，其储粮品种为小麦;试验中通过U型压力计注水读mm水柱数后转化为Pa进行对比验证。

1.2 试验方法

1.2.1 实仓试验

在仓储粮堆机械通风过程中对其外界及内部状态参数变化进行测量，分析测量数据并最终得到一定结论。本研究所采用的CFD模拟研究就是基于实仓试验基础之上的，因此，实仓试验既是模型模拟的前提，也是对数值计算结果进行对比验证的标准。

实仓试验的主要内容为测量仓储粮堆在全底板通风形式下不同通入风量时内部静压分布变化情况。对于压力的测量采用的是特制测压管，测压管长有3种规格，分别为35、55、85 cm，样图如图1所示。

图1 特制测压管示意图

在实仓试验中，3种规格的特制测压管横向入仓深度分别为5、20、50 cm，分别将6根管布置和编号如下:1(西面短管)，2(西面中管)，3(西面长管)，4(南面长管)，5(东面短管)，6(东面中管)，其中3号管位于粮高5.16 m处，2号和 6号管位于粮高2.64 m处，1 号，4 号和5 号管位于粮高0.48 m 处，具体布置如图2所示。

图2 测压管布置示意图

1.2.2 数学模型建立及求解

1.2.2.1 模型的建立及选择

因研究对象为粮食储藏在通风状态下粮堆内压力分布的问题，故模型模拟的区域为包括堆积散粮的整个储仓，仓房内的压力场分布是三维问题，所以对仓房的整体结构建立笛卡坐标。

通风过程中，整个散粮堆积区域内部粮层温度场的分布状态均随时间而变化，整个空间内的相关参数大多也随时间而变，因此粮食仓储的通风过程必须作为一个非定常问题来进行计算。建立通风条件下的粮食仓储的物理模型，须对整体仓房作相应的合理简化，对于粮仓的空气入口及出口采用相应的边界条件处理，忽略通风地槽的导热影响，底面热传导可以忽略不计，内部环境的变化主要由通风地槽的气流输入和仓房壁面与外界的热传导作用所导致，同时将粮食堆积区域看作与粮食物理特性等效的多孔介质区域，通风过程中气体穿过该区域的运动采用FLUENT中的多孔介质模型进行计算［10－13］。

仓储粮堆机械通风过程中，气流由通风地槽进入并穿过整个散粮堆积区域，从通风管道到仓房内，空气流动充分发展，风速较大具有典型的湍流特性，因此必须考虑湍流的问题。

在实际数值模拟中，要根据具体问题的特点来决定所选用的模型。本研究在模拟计算中选用目前应用最为广泛且计算量相对较小、精度合适的标准k－ε模型。k－ε模型主要是基于湍流动能和扩散率，其中湍流动能方程是精确方程，湍流扩散方程是由经验公式导出的方程，该模型对较小压力梯度下的自由剪切流具有较好的拟合度［14］。

1.2.2.2 边界条件及计算区域

入口及出口:流体入口条件按照速度入口边界条件处理，需确定风速大小、入口方向以及湍流参数情况。由于模拟主要涉及压力场的分布，则空气入口温度的设定对模拟的结果没有影响，可将其设为与外界温度一致。流体出口按照压力出口边界条件处理，需确定出口表压及湍流参数情况，采用湍流强度及湍流长度尺度(I－l)来表示湍流情况，其具体相关参数的计算公式如下所示:

湍流强度I的计算公式为:

式中:ReDH为按水力直径DH计算得到的雷诺数，对于圆管，水力直径DH等于圆管直径。

湍流长度尺度l的计算公式为:

式中:L为关联尺寸，对于完全发展的湍流流动，选择在水力学直径流场中指定L=。

仓房壁面边界条件:粮仓壁面为不锈钢结构钢板材料，固体壁面为无滑移边界，靠近壁面处采用壁面函数法进行处理。由于通风时间相对较短，壁面换热边界对整个仓储粮堆机械通风过程中温度场变化的影响很小，因此在模拟过程中将壁面设置固定温度且与外界温度一致。

多孔介质区域:在散粮堆积区域内，谷物颗粒作为多孔介质的骨架在谷物彼此之间构成了一定的空隙，同时固体骨架遍及整个多孔介质所占的体积空间，空隙空间相互连通，在机械通风过程中，气流充满整个介质区域并在谷物空隙之间流动。试验中将散粮堆积区域作为多孔介质区域进行数值模拟。粮食颗粒作为多孔介质的骨架，空气在多孔介质孔隙中流动时要受到粮食颗粒的阻力，包括黏性阻力和惯性阻力。

多孔介质的模拟是通过在流体流动标准运动方程中添加一个运动源项来实现的。黏性损失项(Darcy)和惯性损失项组成运动源项［16］:

式中:Si是i方向(x，y，z)动量源项，D 和C 是指定的系数矩阵。在多孔介质单元中，动量损失和压力梯度相关联，压降与流体速度(或速度方阵)成比例，此运动损失造成了多孔介质内的压力梯度，从而产生了与流体速度(或速度的平方)成比例的压力降。

对于简单的各向同性的均匀多孔介质［15］:

式中:α是渗透性，C2时内部阻力因子，简单指定D和C分别为对角阵和C(其他项为零)的矩2阵。

定义渗透性和内部损失系数的计算公式为［16］:

式中:Dp为谷物的平均直径，φ为空隙率。

1.2.3 CFD 模拟

1.2.3.1 建模及网格划分

选取整体仓房和全底板通风系统为模拟计算区域，采用CFD模拟专用前处理软件GAMBIT进行模拟区域的三维模型建立和网格划分，由于通风为全底板通风模式，故将仓房区域和机械通风区域一起处理，分别建立开环流熏蒸口和开人孔的三维模型并进行网格的划分，具体通风模型和网格划分参见图3至图6。

图6 开人孔网格划分

1.2.3.2 边界条件及参数的确定

整个计算区域的边界条件主要涉及通风过程中进口及出口的边界参数，同时包括在全部流动区域内的流体—空气的状态参数以及壁面热边界相关参数，粮堆小麦本身物性参数及多孔介质相关参数也要进行设置。

入口边界:在通风计算过程中，假设空气为不可压缩流体，计算时采用速度入口边界条件，在FLUENT的边界条件设定中需要确定空气入口风速、方向、温度、密度、黏度以及相应的湍流情况。

机械通风过程中，假设空气温度与外界温度一致为27℃，空气密度 ρ=1.293 kg/m3，动力黏度μ =0.000 018 Pa·s。

根据不同通入风量计算得到入口风速，方向垂直于底板，再由通风入风口平面尺寸计算入口湍流情况，对于入口处的湍流情况采用标准k－ε模型的湍流强度及湍流长度尺度(I－l)来描述，其值可由公式(1)和公式(2)计算得出，具体风速及湍流相关参数的设置见表1。

表1 风速及湍流相关参数

出口边界:流体出口按照压力出口边界条件处理，需确定出口表压及湍流参数情况。因出口边界与外界大气相通，故设置出口表压值为0，相关湍流参数设置亦参考表1。

墙体壁面边界:粮仓不锈钢钢板厚度为5 mm;热传导系数为48 W/m·K;热容为480 J/kg·K及钢板密度为7 800 kg/m3。

多孔介质参数:谷物堆积区域采用FLUENT中的多孔介质模型，需对涉及多孔介质的一些特性参数进行设定［17］。设仓储粮堆为小麦，小麦的热传导系数为 λg=0.13 W/m·K，比热容为 Cg=1 780 J/kg·K，容积密度为ρg=750 kg/m3。同时在FLUENT中需要定义黏性和内部阻力系数，即需要输入渗透性α和惯性阻力因素C2。小麦是具有吸渗特性的颗粒，采用密度法可测得其孔隙度φ为0.468 6，同时可测得小麦平均颗粒直径Dp，代入公式(3)和公式(4)可得到模拟小麦物料物性参数 C2=

2 结果与讨论

2.1 实仓试验

按“1.2.1实仓试验”方法，开环流熏蒸口时各测压管所测定的粮堆内压力如表2所示，开人孔时各测压管所测定的粮堆内压力如表3所示。

表2 开环流熏蒸口时测得粮堆内压力数值

表3 开人孔时测得粮堆内压力数值

2.2 CFD 模拟结果

根据已知边界条件及相关参数，采用定常法SIMPLE算法进行了仓储机械通风的计算，计算时各项参数的收敛度设为10－4，在开环流熏蒸口和人孔时分别计算不同风量条件下仓储粮堆内部压力场的分布。其中开环流熏蒸口时的压力场分布见图7，开人孔时压力场的分布见图8。

2.3 结果分析

依据仓房整体结构建立的笛卡坐标，可以得到6个测压管对应的三维坐标值分别为1(0，0.15，0.48)，2(0，0.2，2.64)，3(0，0.5，5.16)，4( －2.5，0，0.48)，5(0，－2.95，0.48)，6(0，－2.8，2.64)，通过Fluent后处理软件CFD－Post分析模拟结果，分别得到各点的模拟数值，结合实仓试验结果(表2及表3)分别对开环流熏蒸口和人孔时的数据进行验证分析，具体对比分析结果如表4和表5所示。

表4 开环流熏蒸口试验值与模拟值对比

表5 开人孔试验值与模拟值的对比

从表4和表5可以看出，尽管由于受到测压人工读数的误差影响及所建模型部分简化等因素的限制，模拟结果与试验结果还存在一定的误差，但整体看来，模拟结果与试验结果有较高的拟合度，试验结果能够真实的体现机械通风过程中粮堆内部压力场的分布情况。

3 结论

为了验证分析CFD方法模拟机械通风过程中仓储粮堆压力场分布的可行性和准确性，试验以室外浅圆仓为模拟对象，利用CFD方法模拟了打开不同出口及不同通入风量条件下的粮堆压力场，并将模拟结果与实仓试验数据比较分析，得出以下结论:

模拟结果与试验结果存在一定的误差，但两者有较高的符合度，能够真实的体现机械通风过程中粮堆内部压力场的分布情况，证明了利用CFD模拟技术研究仓储粮对机械通风过程中压力场分布是可行的。

由于试验条件和粮仓结构的限制，试验只设定6个测压管进行仓内压力值的监测，获取的数据较少。后期可以选择合适仓房进行大批量的粮堆内部压力的测定，以期得到较多数据，为CFD模拟和后续仓储粮堆粮层阻力的研究提供实例支撑。

［1］胡志超，胡良龙，高刚华，等.我国粮食仓储技术的概况与发展［J］.农机化研究，2007(4):171－173

［2］彭威，张忠杰，任广跃，等.仓储粮堆温度场CFD模拟应用研究［J］.粮油食品科技，2011，19(6):5 －8

［3］任广跃，张忠杰，朱文学，等.粮食干燥技术的应用及发展趋势［J］.中国粮油学报，2011，26(2):124－128

［4］吕新.现代粮食通风技术与粮食储藏生态系统［J］.粮油仓储科技通讯，2003(1):11－13

［5］吴子丹.绿色生态低碳储粮新技术［M］.北京:中国科学技术出版社，2011

［6］隋明波.机械通风对粮食安全储藏的作用及影响［J］.中国科技信息，2010(4):81－82

［7］李琼，汪喜波，杨德勇.CFD方法在仓储粮堆温度场研究中的应用探索［J］.粮食储藏技术，2008，37(3):21 －24

［8］张忠杰，李琼，杨德勇，等.准静态仓储粮堆温度场的CFD模拟［J］.中国粮油学报，2010，25(4):46 －50

［9］王远成，段海峰，武传欣.就仓通风时粮堆内部热湿耦合传递过程的数值预测［C］.中国粮油储藏学会第六届年会，2009(4):167－171

［10］Tsimpanogiannis I N，Yortsos Y C.Scaling theory of drying in porous media［J］.Journals of the American Physical Society，1999，59(4):4353 －4365

［11］Marcus Liew Kai Hoa，Meihua Lu，Yong Zhang.Preparation of porous materials with ordered hole structure［J］.Advances in Colloid and Interface Science，2006，121(1 －3):9 －23

［12］Prat M.Recent advances in pore－scale models for drying of porous media［J］.Chemical Engineering Journal，2002(86):153－164

［13］F Plourde，M Prat.Pore network simulations of drying of capillary porous media［J］.Influence of thermal gradients，International Journal of Heat and Mass Transfer，2003(46):1293－1307

［14］王福军.计算流体动力学分析－CFD软件原理与应用［M］.北京:清华大学出版社，2004

［15］王瑞金，张凯，王刚.Fluent技术基础与应用实例［M］.北京:清华大学出版社，2007

［16］Fluent Incorporated.Fluent5 user's guide［M］.Pittsburgh:FLUENT Incorporated，1998，7:9 －2，6 －17

［17］江帆，黄鹏.Fluent高级应用和实例分析［M］.北京:清华大学出版社，2008.

Study on Simulation for Pressure Field Under the Condition of Mechanical Ventilation in Grain Storage

Ren Guangyue1Peng Wei1Zhang Zhongjie2Wu Zidan2Wang Fang1
(College of Food and Bioengineering，Henan University of Science and Technology1，Luoyang 471003)
(Academy of State Administration of Grain2，Beijing 100037)

This paper is based on the experiments under the condition of mechanical ventilation in grain storage，according to the results of mechanical ventilation in real warehouse experiment，and then established the computational fluid dynamics(CFD)model of pressure field in grain heap.Pressure field distributions are simulated under the air outlet for recirculation fumigation and entrance for workers at the different blast volume respectively.The results showed that the CFD model of pressure field in grain heap could express well with the experimental results.And it prove that the CFD method which is applied to research and analyze of pressure field in grain storage.

grain storage，mechanical ventilation，pressure field，CFD

S226.6

A

1003－0174(2012)09－0090－06

国家科技支撑计划(2009BADA0B04)

2012－01－13

任广跃，男，1971年出生，博士，副教授，农产品干燥技术