CFD模拟优化三传过程

动量传递（Momentum Transfer）是指在流动着的流体中动量由高速流体层向相邻的低速流体层转移的过程，如流体输送、过滤、沉降、固体流态化等，遵循流体力学基本规律。动量传递过程包含多相湍流的复杂过程，涉及变量的范围较广，因此动量传递过程一直未能被精确地预测。应用CFD对动量传递过程中浓度场、温度场、压力场和速度场进行模拟，能够优化动量传递过程，对于建立动量传递的理论基础具有重大的现实意义。

（一）典型案例概述

（1）Yang运用CFD对清洗食品加工设备时所需中间冲洗水的驱替时间、最小耗水量、最小产生量和最小需要量进行了评价，并提出了一种用于计算低用水量冲洗步骤的创新算法。与传统的冲洗方法相比，清洗剂的回收率可达体积的89.3%，中间冲洗水的节约率可达55%以上。

（2）Getahun研究了用于鲜果运输的T形和平板型冷藏集装箱内的气流分布，建立了气流的计算流体力学模型，并进行了实验验证。结果表明，T形冷藏集装箱与平板型冷藏集装箱相比，空气再循环区明显减小，垂直均匀气流运动增强。

（3）Dantas采用耦合流体流动、传热、扩散和失活/降解动力学的方法，对扩散层流条件下的食品热加工过程进行了数值模拟，创建了一个对设计、分析、控制和优化有用的过程虚拟原型。

（4）Ameur研究了无挡板圆柱容器中复杂非牛顿流体的混合特性，比较了U形、双U形、V形和W形切割叶轮叶片的优劣。结果表明，与普通叶片相比，在叶片中引入切口不仅减小了孔洞尺寸，而且降低了功耗。V形或W形切割叶片的节能和混合效果最佳。

（5）Bo利用Fluent软件对高压微射流均质器内部流场进行了数值模拟，研究了均质压力、油水比和温度。正交实验数值模拟得到的最佳参数组合为20℃、70MPa、1∶12。模型的有效性可以为均质器的开发和使用提供重要的参考。

（6）Rinaldi研究了罐体几何形状和淀粉浓度对食品模型传热和流体流动的影响，通过实验数据验证了数学模型的正确性，并将其用于传热和流体流动的比较。整体传热系数受罐体几何形状的影响，在低淀粉浓度下，喇叭形罐体的传热系数最高，在高淀粉浓度时，高圆筒型罐体传热系数最高。对于自然对流，方形罐体平均速度最高，而喇叭形罐体平均速度最低。

（7）Šcˊepanovicˊ 对面团模型进行了非均匀流动的计算研究，通过小振幅振荡剪切、应力松弛、简单剪切启动和压缩实验，得到了模型参数的值。又使用OpenFOAM软件中的有限体积程序对流动问题进行了数值求解。结果显示，平面拉伸对显微组织的破坏作用最大。

（二）典型案例解析

[例4-1] 搅拌介质磨机湿法制备荷叶粉及其动力学研究（俞建峰，2018）。

采用搅拌介质磨机湿法研磨制备荷叶粉。考察不同搅拌转速下30min内荷叶粉的粒径变化情况，借助研磨过程解析模型对荷叶粉湿法研磨动力学进行研究；并基于CFD-DEM（离散单元模型）耦合方法，分析不同搅拌转速下研磨剪切率分布和研磨介质平均碰撞能量的变化情况。

1.模型构建及网格划分

棒销式搅拌介质磨机研磨腔基本结构如图4-11所示。研磨腔内壁直径95mm，长度为100mm，研磨腔容量为0.5L，搅拌器转轴转速200~3000r/min（可调）。

图4-11 介质磨机研磨腔

应用Fluent 17.0的前处理软件ICEM建立搅拌研磨流体计算域网格模型。为简便运算和节省时间，在三维建模过程中只截取研磨腔一段进行模拟，且省略了圆角、倒角等细节。为提高计算效率和计算精度，仿真计算中采用结构化网格，共划分139544个，流体区域结构网格划分如图4-12所示。计算区域分为两部分，包含搅拌器在内的旋转区域和静止区域。

图4-12 流体区域结构化网格划分

2.计算流体力学模型

基于欧拉—拉格朗日方法的CFD-DEM耦合模型对搅拌介质磨机研磨过程进行模拟研究。应用流体力学仿真软件Fluent 17.0进行流场模拟，研磨腔内流体选用水，水是牛顿流体，流体雷诺数Re与流体密度ρ、流体黏度μ有关，计算公式为式（4-21）：

式中 Re——雷诺数；

ρ——流体密度，kg/m³；

u——流速，m/s；

d——当量直径，m；

μ——流体黏度，Pa·s。

当量直径d=0.095m。u取棒销末端线速度，棒销半径r=0.04m，经计算工程中的临界雷诺数为Re_c=2300，当雷诺数Re>Re_c时，流体运动状态属于湍流。不同转速下的流体运动状态如表4-1所示。

表4-1 不同转速时流体流动状态

采用动参考系下滑移网格方法来解决旋转流动问题，滑移网格将计算区域分为两部分，包含搅拌器在内的旋转区域和静止区域。划分网格后，定义静止区域与旋转区域的动静耦合交界面、旋转区域与旋转元件的接触表面均为无相对运动。由于RNG-k-ε在湍流能耗散率计算精度较高，故采用RNG-k-ε湍流计算模型，假设无进出口边界条件，考虑到流体的黏性作用，固壁表面边界条件采用无滑移边界条件。收敛残差设定为0.001。CFD-DEM耦合迭代计算到收敛，Fluent 17.0软件后处理得到流场特性云图，分析流场运动的规律。

3.离散元模型

应用离散单元法模拟软件EDEM 2.7对研磨介质的运动情况进行仿真，离散单元法基于牛顿运动定律来描述每一个颗粒的运动。EDEM和Fluent模拟中采用同一网格模型。搅拌介质磨机研磨腔体及搅拌器材料为钢，研磨介质材料为氧化锆球。表4-2为EDEM中物料属性。

表4-2 颗粒模型的物理属性

EDEM模拟在物理属性中选择颗粒与颗粒、颗粒与几何体的接触模型为赫兹—明德林（Hertz-Mindlin），即无滑动接触模型。设置好全局参数后，利用EDEM软件Simulator模块进行仿真计算。利用EDEM软件后处理部分对研磨介质质量、研磨介质碰撞总次数、研磨介质相对法向平均速度等数据进行提取。

4.流体力学黏性能量耗散率表征

在不可压缩的各向同性湍流能量流动过程中，黏性能量耗散率P是黏度和平均速度梯度的函数。可以用来分析研磨腔中各部分的研磨效果，其定义为式（4-22）：

式中 φ_v——能量耗散函数。

能量耗散函数φ_v的定义为式（4-23）：

式中 u——x方向分速度，m/s；

v——y方向分速度，m/s；

w——z方向分速度，m/s。

在Fluent中不能直接取得黏性能量耗散率P的定义，研究认为可以选用剪切率S来替代表征：

式中 S——剪切率。

由于水为牛顿流体，黏度是一常量，剪切率S与黏性能量耗散率P平方根成正比，可以用来表征搅拌介质磨机研磨腔局部研磨效果。

5.结论

通过研究搅拌转速对剪切率分布（图4-13）、荷叶粉粒径、荷叶粉研磨速率和研磨介质碰撞能量的影响，得出结论：荷叶粉湿法搅拌研磨过程符合一级研磨动力学方程。研磨过程中，随着荷叶粉粒径由大变小，荷叶粉研磨速率也随之减小。研磨过程中存在两种破碎方式：冲击破碎与摩擦破碎。提高搅拌转速，冲击破碎方式更加显著，荷叶粉颗粒获得能量增加，荷叶粉研磨速率也会随之增加。研磨介质平均碰撞能量与不同粒级研磨速率之间存在线性关系。提高搅拌转速，荷叶粉各粒级的研磨速率也会随之增加。模拟仿真结果可以应用于对荷叶粉研磨速率的预测。影响荷叶粉研磨速率的因素还有研磨介质颗粒大小、研磨介质填充率以及荷叶粉初始粒径等，相关内容有待进一步深入研究。

图4-13 不同转速下剪切率分布

热量传递（Heat Transfer）是指在物体内部或者物体之间，热量因温差的存在而自发地由高温处向低温处传递的过程，如加热、冷却、蒸发、冷凝等，遵循热量传递基本规律。应用CFD对热对流、热传导和热辐射进行模拟，能清晰地了解热量传递的过程，对计算食品的传热系数、优化杀菌过程、设计冷冻柜等都具有重要的现实意义。

（一）典型案例概述

（1）Andreas建立CFD模型探究高压均质器中液滴的温度分布，并研究在均质器中产生的摩擦热是否导致对温度敏感的分子（如营养物质）的热降解。结果表明，由于停留时间短，无论在何种工况下，高压均质器内都不会发生热降解。

（2）Zhang将评价指标与CFD技术相结合，研究不同尺寸、不同位置的开舱门对冷藏车内部温度的影响。这不仅对打开冷藏车舱门时车厢内的流量分布和温度变化有了更详细地了解，而且为冷藏车的优化设计提供了可靠的理论依据。

（3）Woo将传热模拟模型与基于模糊C均值聚类方法的分级系统成功地结合起来，优化了包装整粒玉米的杀菌过程。该方法根据玉米穗的大小和重量调整并获得最佳杀菌工艺，改善了玉米的质地和感官特性。

（4）Mukama通过实验和模拟的双重方法，研究了石榴强制风冷过程中，容器设计、塑料衬套和堆置方向对气流、冷却速度和能耗的影响，为减少预冷过程中的用电量提供了一种较为经济可行的方法。

（5）Gruyters利用基于X射线计算机断层成像图像的几何模型生成器，建立了可变的三维苹果模型，用CFD模拟了苹果的强制空气冷却过程，并与用等效球表示的结果进行了比较。结果表明，改进的计算机辅助设计方法有助于模拟更精确的对流冷却过程。在下一步，这些模拟将用于多目标优化包装方面的冷链效率和冷却均匀性。

（6）Pasban提出了一种模拟苹果切片对流干燥过程中三维传热、传质耦合过程的数值方法，利用Fluent软件对外流场和温度场进行了数值模拟。结果显示，模拟结果与实验结果吻合较好，验证了该方法的鲁棒性、计算效率和精度。

（7）Schaer提出了一种利用流变学测量和CFD方法确定换热器中酸乳结构损失的耦合方法。结果表明，该模型可以诊断机械应力，评估酸乳样品的纹理损失，并提供生产商无法获得的信息，以改进他们的工艺，优化酸乳等乳制品的生产线。

（二）典型案例解析

[例4-2] 基于CFD数值模拟的豆腐干软罐头杀菌工艺优化（王磊，2017）。

运用CFD对豆腐干软罐头恒温热杀菌工艺进行优化，能够解决杀菌过程中存在的因杀菌过度而品质不佳问题。在传热学和品质动力学的基础上建立优化计算，构建豆腐干杀菌非稳态固体热传导模型，进行CFD数值模拟，获取豆腐干各节点的温度历史，再寻找杀菌值F不低于安全限值而蒸煮值C最小的升温方式，然后通过实验验证可行性，最终结合生产实际情况确定一种最为合理的升温方式。

1.非稳态固体热传导模型

豆腐干软罐头传热的完整数学描述由固体内部传热控制方程、初始温度分布条件和对流传热边界条件构成。

（1）控制方程 豆腐干软罐头内部导热过程适用简化的三维非稳态导热微分方程：

式中 ，拉普拉斯算子，x、y、z代表笛卡尔坐标系；

k——固体的导热系数，W/（m·℃）；

T——固体微元的温度，℃；

ρ——固体密度，kg/m³；

C_p——固体的比热容，J/（kg·℃）；

t——传热时间，s。

（2）初始条件和边界条件

①初始条件：初始温度为T₁=15.8℃。

②边界条件：流体—颗粒对流加热过程中，其边界控制方程为式（4-26）：

式中 ，哈密顿算子，i、j、k分别代表x、y、z坐标轴上的单位矢量；

h——对流传热系数，W/（m²·℃）；

T_s——固体表面温度，℃；

T_f——对流温度，℃。

（3）几何模型的构建及网格划分 豆腐干建模如图4-14（1）所示；豆腐干模型网格划分如图4-14（2）所示，其中单元数量为99284个。

图4-14 豆腐干建模及网格划分

2.限制条件

为确保豆腐干软罐头的安全性，优化方法以微生物热致死模型作为限制条件。使用Z值模型，F_z值表示参考温度在121.1℃的微生物等效致死时间，见式（4-27）：

式中 t——杀菌时间，s；

T_sh——样品冷点温度，℃；

T_ref——参考温度，取121.1℃；(www.xing528.com)

z——微生物对热的敏感性，其值为D值变化一个对数值所需温度，通常取10℃。

3.目标函数

在满足限制条件的基础上，为进一步筛选能够有效降低品质破坏的升温方式，采用蒸煮值模型作为品质指标，即优化的目标函数。

（1）蒸煮值 蒸煮值（Cooking Value）简称C值，表征品质破坏的程度。C值越大，则品质破坏越严重，反之说明品质保存率高，见式（4-28）：

式中 T_ref——100℃；

z_q——品质对热的敏感性，整体品质通常取33℃。

（2）中心C值 中心C值，即C_c，表征中心品质破坏的程度，见式（4-29）：

式中 T_c——样品中心温度，℃。

（3）表面C值 表面C值，即C_s，表征表面品质破坏的程度，见式（4-30）：

式中 T_s——样品表面温度，℃。

（4）体积平均C值 由于数值计算可以获得模型几何空间中任意节点的温度时间历史，因此可对食品进行整体的动力学分析。C_avg能够较为全面地评价样品整体的品质保持率，见式（4-31）：

式中 V——样品体积，m³。

4.温度数据采集及传热学、动力学参数数值计算

（1）杀菌全程温度数据采集

①采集杀菌釜对流温度：杀菌完成后，根据采集的数据拟合时间—温度函数，并取最慢加热区的温度函数。

②采集样品表面及杀菌冷点温度：表面、中心温度采集装置如图4-15（1）所示。温度采集位置如图4-15（2）所示。灭菌完成后，表面和中心的温度均取最慢加热区温度历史，且该样品中心为灭菌冷点。

图4-15 温度采集示意图

（2）计算对流传热系数h根据采集的对流温度函数，假设系列h值，应用CFD软件输入式（4-25）、式（4-26）及物理参数进行杀菌全程的模拟，将模拟的与试验采集的冷点温度数据按照最小温度目标总体差平方和法（Least Summation of the Squared Temperature Difference for Overall Target，LSTD）计算并搜索升温、保温和降温的对流传热系数h，见式（4-32）：

式中 LSTD——温度差平方和，℃²；

T_sn、T_cn——在共为m个的第n个时间点分别由数值模拟和采集获得的温度，℃。

（3）计算F值、C_s值和C_avg值 根据采集的冷点和表面温度历史分别计算F值和C_s值。应用CFD软件，输入式（4-25）、式（4-26），构建三维非稳态固体传热方程，并将前两步获得的对流温度和对流传热系数输入边界条件，求解计算，后处理器导出所有网格单元的温度历史计算C_avg值。

5.优化搜索模型

（1）预设一系列不同的梯度升温模式，随后应用CFD软件模拟原恒温杀菌过程和预设模式，后处理器导出冷点、表面及所有网格单元的温度历史。动力学计算模拟得到F值，在F值达到标准后计算该升温方式的C_s值，相比于原有的恒温杀菌方式（116℃保温30min），当C_s值有效减少后（>5%），则认定该升温方式能够满足安全标准并能有效减小品质损失。

（2）将选出的梯度升温杀菌方式进行试验验证，根据实际采集的数据计算实际F值、C_s值；根据模拟结果计算C_avg值。

（3）根据计算结果，结合实际生产的要求，最终确定一种梯度升温模式作为豆腐干新杀菌工艺。

6.结论

研究最终得到一种既能满足安全指标，品质劣化程度又显著降低的豆腐干软罐头梯度升温杀菌工艺。该杀菌工艺温度峰值能够满足工业化生产需求。研究结果为固体软罐头食品的变温杀菌工艺优化提供参考。

质量传递（MassTransfer）是指物质在介质中因化学势差（浓度差）的作用发生由化学势高的部位向化学势低的部位迁移的过程，如蒸馏、吸收、萃取、干燥等，遵循质量传递基本规律。直接测定质量传递中物质的流量、流速是相当困难的。应用CFD可以解决这些困难，实现整个质量传递过程的数值模拟，更好地优化质量传递过程及其设备。

（一）典型案例概述

（1）Olenskyj采用CFD方法分析了玉米醇溶蛋白纳米粒凝聚过程中等值面的对流通量、扩散通量、速度大小、压力和剪切应力，发现只有压力始终与纳米粒直径的减小和多分散指数的增加相关。这一发现可能会影响微流体抗溶剂沉淀工艺的未来发展。此外，CFD模拟可用于预测随着工艺参数的变化而变化的纳米粒性能，从而减少物理试验。

（2）Ramachandran采用商用CFD软件包，将湿筒干燥模型与过热蒸汽外流模型相结合，对酒糟颗粒干燥过程中的传热、传质现象进行了数值研究。结果表明，采用这种耦合方法进行模拟，可以在不确定质量边界层的情况下，得到球团内部和界面处质量通量的完整分布。该模型可作为大型过热蒸汽干燥系统设计和优化的初步工具。

（3）Orona采用耦合的CFD-DPM（离散相模型）-DEM模拟方法，研究了一种具有运动食物球和多个射流的流体化系统，研究了操作变量（流量、温度）和球孔数对动量、传热和传质的影响。结果表明，CFD-DPM-DEM是对食品加工系统进行相对较低成本仿真的强大工具，可以更真实地描述该系统。

（4）Khampakool研究了红外辅助冷冻干燥生产香蕉零食的效果，利用指数模型、佩奇（Page）模型和扩散模型对干燥动力学进行了评估。结果表明，红外辅助冷冻干燥具有生产高品质香蕉零食的潜力，并在干燥过程中节省了大量的时间成本、人力成本。

（5）Woo研究了大豆在25℃复水过程中水分含量的变化。采用法勒（Peleg）模型估计大豆达到目标含水率（33.33%）所需的浸泡时间为14.59min。采用大本（Omoto）模型和仿真模型，通过两个计算步骤确定了传质模拟的传质系数k。确定最适合模拟的k值为6.0×10^-7m²/s，高于Omoto模型得到的表观k值。

（二）典型案例解析

[例4-3] 基于Fluent的菊花热风干燥流场特性仿真分析（李赫，2018）。

以ZDG230型负压式电加热干燥机为试验样机，运用FLUENT软件分析菊花烘干过程中干燥机内部的流场分布特性，为今后优化菊花干燥工艺、提高菊花干燥品质提供参考。

1.几何建模和网格划分

ZDG230型负压式干燥机干燥室的内尺寸为125cm（长）×93cm（宽）×127cm（高），建模时默认选用干燥机的中心为坐标系的原点。干燥室的菊花层数设计为10层，则菊花堆积区域（多孔介质区域）在z坐标轴上的对应区间为（46.5cm，78.2cm），干燥机的物理模型如图4-16所示。

图4-16 干燥机物理模型示意图

1—热交换管 2—进气口 3—排湿口 4—风机叶片 5—风机电机 6—风机变频器 7—自动控制器 8—多空风道

9—多空盘 10—门锁夹 11—多空突出盘 12—托盘架 13—托盘 14—地轮 15—门合页

16—托盘支架 17—温度传感器 18—门密封材料 19—电加热器固定台 20—电加热器

通过前处理软件GAMBIT进行网格划分，为兼顾计算精度和运算效率，采用六面体网格单元分块划分网格的方案。在干燥机的菊花堆积区域（多孔介质区域）采用较密的网格划分格式，而在其他区域则采用相对稀疏的网格划分格式，共产生2.45×10⁶个网格单元。

2.基本假设

根据负压式电加热干燥机两侧送风的特点，设定湿菊花的热风干燥过程为典型的多孔介质热质传递过程，对传热过程做以下假设：

（1）干燥机的干燥室内分布有装料门，从简化数学模型的角度出发，认为其物性参数和聚氨酯的物性参数是一致的。

（2）干燥机内空气为不可压缩的理想气体。

（3）菊花与菊花之间的热传导忽略不计。

（4）由于热空气温度不高，所以仅考虑对流传热，不计辐射传热。

（5）干燥机的通风排湿顺畅。

（6）干燥机的聚氨酯箱壁为绝热体，其热容量忽略不计。

3.计算模型建立

研究在菊花烘干过程中干燥机内部流场分布及变化规律。式（4-33）为质量守恒方程的表达式：

引入哈密顿算子，式（4-33）可写为：

式中 ρ——密度，kg/m³；

t——时间，s；

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

u、v、w——速度矢量μ在x、y、z方向的分量，m/s。

标准k-ε两方程模型是由湍流动能k方程和湍流动耗散率ε的方程共同构成的：

式中

G_b——由于浮力引起的湍流动能k的产生项；

G_k——由于平均速度梯度引起的湍流动能k的产生项；

Y_M——代表可压湍流中脉动扩张的贡献，Y_M=0；

σ_k、σ_ε——与湍流动能k和耗散率ε对应的普朗特数，分别为1.0、1.3；

S_ε、S_k——用户自定义的源项；

C_1ε、C_2ε、C_3ε——经验常数，分别为1.44、1.92、0.00。

4.边界条件

选用湍流参数为湍流强度I和水力直径D_H来定义流场边界上的湍流。根据流体力学的经典理论，湍流强度I和水力直径D_H，可由式（4-37）、式（4-38）得到：

式中 Re_DH——以水力直径D_H为特征长度求出的雷诺数；

A——截面面积，mm²；

P_w——湿周（即过流断面上流体与固体壁面接触的周界线），mm。

选用压力出口边界条件，出口气流的压力为50.1kPa，出口气流的温度为实测温度40℃干燥机壁面为绝热壁面。

5.结论

应用Fluent软件对模型的求解计算实现了菊花热风干燥过程中流场分布的数值模拟，得到了菊花干燥过程中的速度场（图4-17）、压力场（图4-18）及温度场（图4-19）的分布规律。在菊花放置区域干燥机内流场速度随烤箱高度而增大，风速大致在0.47~8.04m/s时呈梯度变化；干燥室内压力沿干燥机z轴方向不断减小，不同高度平面上差别明显，在同一高度上各托盘内物料压力场分布较为均匀；在干燥初始阶段菊花放置区域初期温度场分布最大温差接近10℃，且其温度场分布在竖直轴心方向温度较低，随着干燥时间的推移，其温度场分布逐渐趋于均匀。

图4-17 干燥机内部速度场分布图

图4-18 干燥机内部压力场分布图

图4-19 干燥机内部温度场分布图