数值模拟技术预测风机两种摆放方式对冷库堆垛货物的影响

杜子峥 谢 晶 朱进林

（1.上海海洋大学食品学院，上海 201306；2.上海水产品加工与贮藏工程技术研究中心，上海 201306）

果蔬等农产品在进入超市，供给消费者前，需在冷库中进行较长时间贮藏，冷库内部合理的温度分布可显著延长商品货架期，减少商品损耗，提高经济效益。为研究冷库内气流的分布机理以改善货物贮藏条件，国内外研究人员［1－6］对冷库内气流组织分布做了大量的研究，主要通过计算流体力学模拟和实验验证两种方式相互验证。初期的研究主要通过对冷库空库模型进行数值模拟，使学者对冷库内部气流组织的温度、风速、涡旋分布有了初步认识。然而冷库主要功能是贮藏货物，货物不同的摆放方式、堆垛尺寸、潜在的蒸发潜热、呼吸热均会显著影响其内部气流组织，因此仅对空库进行模拟研究远远不够，必须考虑货物对冷库气流组织的影响。刘妍玲等［7］在果蔬摆放方式研究中发现中间有通道式果蔬摆放方式，气流组织分布较合理。Chourasia等［8］利用数值模拟技术，围绕冷库土豆堆垛宽高比、堆垛间距做了大量的研究，发现加大货物垂直间距可显著改善货物温度分布。Delele等［9］在对侧吹风式水果冷库加湿系统CFD模拟研究中，建立了单个堆垛货架模型，并对货架采用多孔介质模型模拟，成功预测了冷库内部空气相对湿度、速度以及货架温度分布，发现货物顶部温度最高，货物最热区域与冷库空气侧存在0.92℃温差。Tanaka等［10］对多种货物装载方式进行数值模拟研究，结合试验测量值，确定了货物最优化堆放方式。但是以上研究并没有在冷库满载情况下进行。

本研究拟以笔者实际参观走访的一家中型冷库为蓝本，建立模型，利用数值模拟技术，在货物满载下研究对吹风式和下吹风式两种风机摆放方式对库内气流分布以及堆垛货物温度分布的影响，旨在寻找适合中型冷库的最佳风机设置方案，为未来冷库建设提供一定的参考依据。

1 模型和方法

1.1 物理模型建立

研究对象为一间24.3m×21.6m×7.2m的中型冷库，库内贮藏果蔬等农产品，配备有货架以方便货物堆垛，冷库实际安装4台吊顶式的Kuba GEA冷风机，采用侧面对吹方式，风机总功率为8.96kW，每台风机有3个直径为0.45m的圆形风口，出风温度为0℃，风速12m/s。库内共有42个货架，货架间距0.5m，离壁面距离1m，货架底部距离地面0.3m，货物堆垛尺寸为4.5m×1.2m×5.5m。冷库中间设有宽度为2m走道，以方便人员进行堆垛和捯垛操作。为研究风机不同摆放方式对库内货物影响，分别建立两种模型：① 吊顶风机对送风模型，该模型以实际冷库作为模板，风机设置在两侧墙面，结构见图1（a）；② 风机上回下送风模型，风机设置在冷库顶部，且出风口下侧处于货架间隙，模型结构见图1（b）。由于两种模型均具有对称结构，不影响计算结果的前提下，为缩短计算时间，只对模型取对称面一侧进行模拟。

图1 两种不同风机摆放冷库模型Figure 1 Two different cold store models with various fan arrangement

为简化计算，对模型进行如下假设：

（1）冷库壁面假设为无滑移恒温壁面，壁面温度恒定为3℃，冷库内部密封良好，不考虑外界渗透作用的影响；

（2）冷库气体为不可压缩的理想气体，符合Boussinesq假设；

（3）货物侧存在呼吸热，且呼吸热为恒定值；

（4）货物侧为多孔结构，且内部气体是层流流动；

（5）不考虑内部因潜热导致的温度变化。

1.2 堆垛货物模型方程

冷库内部堆放货物被视为多孔介质，多孔介质模型的动量方程是在标准动量方程的后面加上动量源项。货物侧流动阻力可以参照Darcy-Forchheimer公式：

其中右侧第一项为考虑粘性损失的达西公式项，右侧第二项为惯性损失项。Darcy-Forchheimer公式通常也可以表示为：

式中：

p——压强，Pa；

u——速度矢量，m/s；

K——通过多孔区域的达西渗透率，D；

C——惯性阻力系数。

当局部雷诺数大于1时，必须考虑惯性损失项。局部雷诺数公式：

式中：

u——速度矢量，m/s；

ρ——流体密度，kg/m3；

dpe——微粒的有效当量直径，m。

当量直径可由式（4）进行计算：

式中：

V——颗粒的体积，m3。

为简化模型，将货物看做均质的体积相同的个体。在多孔介质模型的使用上，模拟多孔板或者管束系统，可以忽略渗透项，只使用惯性损失项。但是对堆垛货物忽略渗透项会导致模拟精度下降，无法实现货物内部气流的准确预测，需考虑堆垛货物内部渗透项系数，使用ergun方程［11］的填充床模型对货物侧进行计算，ergun方程表示为：

式中：

L——填充床厚度，m。

渗透阻力系数和惯性阻力系数可表示为：

式中：

ε——多孔介质区域的孔隙率，%；

λ——形状系数。

孔隙率及形状系数公式：

式中：

A——颗粒表面积，m2；

ρb——货物侧区域密度，kg/m3；

ρp——货物密度，kg/m3。

1.3 湍流模型选择

工程应用中湍流的数值模拟主要分为三类：直接数值模拟（DNS）、大涡模拟（LES）和基于雷诺平均N—S方程组（RANS）模型。由于计算机条件约束，前两种模拟方式仍在探索阶段，RANS模型在工程应用中使用最广泛。冷库湍流模型主要使用k—ε湍流模型、SSTk—ω湍流模型以及RSM湍流模型。早期研究一般采用工程通用的k—ε模型，但其预测精度仍有待检验，Delete等［12］对多种两方程模型（standardk—ε、RNGk—ε、realizablek—ε、standardk—ε和SSTk—ω湍流模型）预测精度进行检测，上述模型对平均风速预测的相对误差分别为24.3%，22.4%，23.5%，18.2%，SSTk—ω湍流模型预测较准确，因此，本次模拟采用SSTk—ω湍流模型。

1.4 多孔介质区域参数设定

将货物侧视为多孔介质区域，考虑货物内部渗透性较差以及其几何尺度比湍流涡的尺度要大，因此必须抑制多孔介质区域湍流影响，使用laminar zone模型对湍流生成进行抑制，多孔介质区域视为各项同性。货物侧堆垛为有效直径80 mm的苹果，其呼吸热可依据公式Sef＝1.91Tf＋4.64（Anno）［13］计算。苹果物性参数［14］：密度为837kg/m3，热导率为0.558W/（m·℃），比热为3 665.6J/（kg3·℃）。苹果表面积A可根据公式A＝8.849－0.835 9dpe＋0.032 5［15］进行计算，货物侧渗透系数、惯性系数以及其他参数见表1。

表1 多孔介质区域具体参数设定Table 1 Specific parameters of the porous medium region

1.5 边界条件设定及计算方法

使用非结构网格对模型进行网格划分，计算域网格尺寸为1.4cm，下吹风式冷库生成3.3×106个网格，两种模型分别生成3.25×106和3.26×106个网格。风机出口采用Velocity inlet速度出口，回风边界设为outflow，紊流强度和水力直径定义湍流，湍流强度根据公式I＝0.16Re－1/8求得为5%，水利直径为出风口特征尺寸0.45m。所有固体表面设为无滑移壁面，冷库外墙和地面为现场聚氨酯发泡，其热导率为0.022W/（m·K），墙体厚度0.15m，外墙和地面采用定热流密度边界条件，热流密度可根据式（10）进行计算。

式中：

λw——墙体热导率，W/（m·℃）；

tw——墙体厚度，m；

Tw、T——分别为室外温度（10℃）和冷库内部温度（0℃），℃。

经计算其热流密度为1.47W/m2。

考虑因浮力驱动的自然对流的影响，引入Boussinesq假设。利用fluent对两种模型进行稳态模拟以获得冷库内部气流相对稳定时的气流分布，连续性残差设定为10－3，其余各值设为10－4。

2 稳态模拟结果分析

通过fluent计算之后，得到两种风机摆放方式下冷库内部速度及温度数据，分别对它们主流截面、货架底部截面及货物侧进行速度及温度的对比分析。

2.1 速度分布及对比

由于风机设置形式不同，出风主流截面风速分布差异显著（见图2）。下吹式风机在出风口左右两侧形成对称的回风主流，风速为0.6m/s；对吹式风机回风主流出现在风机下侧，风速为0.4m/s。风机下吹冷库在货物区域平均流速高于风机对吹型冷库，风机对吹冷库货物侧流速主要处于0.1m/s的较低值，原因有两方面：① 由于侧吹风机在重力和货物呼吸热产生的冷热气流密度差影响下，流速衰减严重，导致风机远端存在大面积送风死角；② 由于货物对侧吹气流产生阻碍作用。这一结果可能导致出风主流远端货物区由于冷量不够影响换热，以致部分货物区域温度较高。

图2 两种不同风机摆放下吹风主流速度分布Figure 2 Comparisons of velocity distribution at the mainstream with two different fan arrangement

图3为出风主流处涡旋分布截面图。风机下吹库在垂直方向近地面处发现对称的回流涡旋（图3（a））。风机对吹库在水平方向靠近风机左右两侧存在对称的回流涡旋（图3（b）），风机对吹库形成的回风循环较短，大部分冷空气只通过货物中部便加入回风循环，并未实现均匀送风。堆垛贮藏中，一般设立垫仓板使货物与地面间隔一段距离，以实现货物底部充分回流换热并方便堆垛货物。图4为两种风机摆放在冷库底部的风速分布，风机下吹方式在垂直方向未受到货物阻碍，且货架的间隔形成天然风道，对冷空气进行分流，平均流速高于风机对吹库，风速呈现放射状衰减，不存在回流死角，风机对吹库风速呈现阶梯状衰减，在货物底部形成大面积回流死角（见图4（b）），风速仅为0.01～0.02m/s。在同样风量条件下，风机下吹库回流更加均匀，有利于冷库内气流通畅流动。

2.2 温度分布及对比

图3 两种不同风机摆放下吹风主流处涡旋分布Figure 3 Comparisons of vorticity distribution at the mainstream with two different fan arrangement

图4 两种不同风机摆放下冷库底部流速分布Figure 4 Comparisons of velocity distribution at the bottom of cold store

图5 两种不同风机摆放下货物侧温度分布三维图Figure 5 Temperature distributions at the cargo side

图5为冷库货物温度分布三维视图。由图5可知，下吹风方式在货物区温度均匀程度优于对吹风方式，风机下吹风摆放方式在货物侧气流分布不但均匀，且平均流速高于侧吹风冷库，使得冷空气能够较好地与货物进行换热，下送风冷库货物最高温度在近风机两侧，温度为274.8K，对吹式冷库风机远端货物温度较高且高温区面积较大，最高温度可达276K，而该区域恰恰存在上述提到的低流速回流死角，可见在相同风量条件下，冷库内部充分均匀回流对温度场合理分布起关键作用，这同时也印证了Delete等［12］在研究中的发现。

图6为两种风机摆放方式下冷库底部温度分布对比，下吹风式风机在底部温度较低并且分布均匀（图6（a）），冷库气流在货物间隔底部形成两个低温带，多方向主流回风循环。对吹式冷库在底部存在温度梯度，风机射流远端温度较高，中间区域温度较低（图6（b））；在射流末端货物区域形成了多个涡旋（图7），涡旋对冷库气流产生动能耗散作用，使射流末端货物换热进一步恶化，导致末端货物温度较高。

图6 两种不同风机摆放下货物底部温度分布截面图Figure 6 Comparisons of temperature distribution at the bottom of cold store

图7 对吹风冷库涡旋速度分布截面图Figure 7 Vorticity distributions of the airflow in cold store with back blowing fans

3 讨论

在本研究中，与风机对吹方式比较，下吹方式可以显著改善冷库及货物内部的温度分布。对吹风冷库在风机吹风主流下部及射流末端货物区温度较高，这是由于对吹风主流在受到热压的影响下流速衰减，在射流末端产生多个涡旋，加上气流受货物阻碍形成较短回流循环，导致流经风机远端形成大面积送风死角，影响货物换热。单台下吹式风机在库内可形成多方向主流回风循环，冷库内回流方向增加无疑也对气流分布起到改善作用。此外，下吹式风机直接对货物进行送风降温，冷量的利用率高于侧吹式风机，货物区域温度更低也更均匀。在本研究中，冷库采用风机下吹方式冷却效果优于对吹方式。目前，中国大部分冷库采用侧吹风型风机，下吹风型风机虽在送风均匀性上有优势，但在实际使用中却有众多限制，比如出风主流易受货物阻碍，货物区域流速过高干耗严重等等。

4 结论

选择将货物视为多孔介质，即将货物视为有大量密积微小空隙构成的物质，通过计算货物内部粘性阻力、惯性阻力、孔隙率及其他参数以实现货物内部气流组织准确预测，这种方式的缺点是使建立的冷库堆垛货物模型复杂化，模拟计算时间明显延长，收敛速度显著变慢，但更接近实际情况。通过对模型模拟计算结果对比，得到以下主要结论：

（1）建立了两种不同风机摆放方式下冷库堆垛货物模型，并利用多孔介质模型对货物侧温度和速度分布进行模拟，获得冷库内部气流组织速度和温度分布，可为优化冷库货物摆放和风机设置提供理论参考。

（2）货物贮藏过程中，气流组织分布均匀性是影响货物温度分布的关键因素，对吹风式冷库由于受热压和重力影响，风机下侧货物和风机远端易形成送风死角，温度较高，下吹式冷库在货物区域的流速和温度分布较对吹式冷库均匀。

（3）本研究对比分析了两种风机摆放方式下冷库内部温度场和速度场分布情况，未来研究中，可对比研究两种方式下冷库相对湿度的变化。

1 荆有印，王长海，丁桂艳，等.旋转气流及货物对冷库流场影响的研究［J］.冷藏技术，2008（4）：34～37，41.

2 缪晨，谢晶.冷库空气幕流场的非稳态数值模拟及验证［J］.农业工程学报，2013（7）：246～253.

3 Xie Jing，Qu Xiao-hua，Shi Jun-ye，et al.Effects of design parameters on flow and temperature fields of a cold store by CFD simulation［J］.Journal of Food Engineering，2006，77（2）：355～363.

4 汤毅，谢晶，王金锋，等.三维计算流体力学技术用于预测小型冷库内气流分布［C］//中国制冷学会.第八届全国食品冷藏链大会论文集.北京：中国制冷学会，2012.

5 缪晨，谢晶.空气幕的研究进展［J］.食品与机械，2012，28（4）：237～262.

6 汤毅，谢晶，王金锋，等.计算流体力学在冷库优化中的应用研究进展［J］.食品与机械，2011，27（5）：186～189.

7 刘妍玲，张岩，王世清，等.果蔬摆放形式对冷库内气流场分布影响的研究［J］.青岛农业大学学报（自然科学版），2008，25（1）：24～27.

8 Chourasia M K，Goswami T K.Simulation of effect of stack dimensions and stacking arrangement on cool-down characteristics of potato in a cold store by computational fluid dynamics［J］.Biosystems Engineering，2007，96（4）：503～515.

9 Delele M A，Schenk A，Tijskens E，et al.Optimization of the humidification of cold stores by pressurized water atomizers based on a multiscale CFD model［J］.Journal of Food Engineering，2009，91（2）：228～239.

10 Tanaka F，Konishi Y，Kuroki Y，et al.The use of CFD to improve the performance of a partially loaded cold store［J］.Journal of Food Process Engineering，2012，35（6）：874～880.

11 Ergun S，Orning A A.Fluid flow through randomly packed columns and fluidized beds［J］.Industrial ＆Engineering Chemistry，1949，41（6）：1 179～1 184.

12 Delele M A，Schenk A，Tijskens E，et al.Optimization of the humidification of cold stores by pressurized water atomizers based on a multiscale CFD model［J］.Journal of Food Engineering，2009，91（2）：228～239.

13 Delele M A，Vorstermans B，Creemers P，et al.CFD model development and validation of a thermonebulisation fungicide fogging system for postharvest storage of fruit［J］.Journal of Food Engineering，2012，108（1）：59～68.

14 Lisowa H，Wujec M，Lis T.Influence of temperature and variety on the thermal properties of apples［J］.International Agrophysics，2002，16（1）：43～52.

15 Schotsmans W，Verlinden B E，Lammertyn J，et al.Factors affecting skin resistance measurements in pipfruit［J］.Postharvest Biology and Technology，2002，25（2）：169～179.