基于COMSOL Multiphysics对含顶隙液态罐头食品热杀菌的模拟研究

陈舒航，周建伟，刘东红，4*，刘 懿

（1 浙江大学生物系统工程与食品科学学院 杭州310058 2 浙大宁波理工学院 浙江宁波315100 3 江苏省食品先进制造装备技术重点实验室 江苏无锡214122 4 浙江大学馥莉食品研究院 杭州310058）

罐藏是一种常见的食品贮藏方法，广泛运用于食品工业中。热力杀菌因其具有高效、方便的特点，是食品加工中最常用的杀菌手段[1]。罐藏食品的热杀菌涉及流体传热，用传统方法难以得到罐头内的任意时刻每个位置的温度分布情况，对罐头内冷点和最低杀菌强度判断的准确度不足[2]。借助计算机对热杀菌过程进行仿真模拟，能有效预测不同类型食品在该过程中温度、速度、压力、最慢加热区（Slowing heating zone，SHZ）及杀菌强度的变化[3-8]。将仿真结果应用于调整几何尺寸、热处理温度、处理时间、旋转速度等参数[9-12]，在保证杀菌强度的同时，减少杀菌时间，提高产品风味和营养水平，节约能耗。精准的仿真模型可提高试验的精准性和有效性，为设备研发及确定最优控制方案提供依据。

罐头顶隙是指罐藏食品的表面或液面与罐头容器上盖间所留空隙，可维持一定的真空度，平衡罐头内外的压力差，避免在热杀菌过程中出现凸角、瘪听等外观变形的情况[13]。目前对罐头热杀菌过程的模拟大多忽略了顶隙对传热的影响，且未考虑顶隙处存在的液体蒸发和冷凝现象，因而无法对杀菌过程中的压力变化进行准确预测[14-15]。忽略顶隙也会导致对温度的预测产生偏差，使得最小杀菌值的计算不够准确。

羧甲基纤维素钠（Carboxymethylcellulose，CMC）具有与一些液态食品相似的流变学特性，常被用于液态罐头食品热杀菌的模拟过程[16]。本文选取1.0% CMC 溶液，借助COMSOL Multiphysics软件，构建了包含及不包含顶隙的2 种模型模拟罐头食品的热杀菌过程。通过热穿透试验，由模拟值与真实值之间的平均相对误差验证模型可靠性，以期获得一种较准确的罐头热杀菌过程的预测模型，更精准地分析热杀菌过程中速度场、温度场、致死率场及压力的变化情况。在含顶隙模型的基础上，探究不同顶隙高度及保温温度对1.0%CMC 罐头内相对压力的影响，为工厂生产罐头时所预留的顶隙高度及不同工作温度下杀菌釜的压力控制提供参考。

1 材料与方法

1.1 材料试剂

金属罐（罐号7113），山东环球印铁制罐有限公司；纯净水，娃哈哈公司；羧甲基纤维素（化学纯级），国药集团化学试剂有限公司。

1.2 仪器设备

Data Trace MP Ⅲ无线温度传感器及工作站，美国Mesa Laboratories 公司；YJ-C20 型手动封罐机，张家港易捷自动化设备；50 L 安全智能型反压高温蒸煮锅，济南兰光机电技术有限公司。

1.3 试验方法

采用水和1% CMC 溶液分别对数值模型进行试验验证，具体试验步骤如下：

1）水的模型验证：用自制支架将温度传感器固定在罐头中心，探针高度位于距底部15 mm 处。在室温下，将纯净水灌入马口铁罐头中，使液面高度为95 mm，准备3 罐，用封罐机将马口铁罐封口后放入反压高温蒸煮锅中，同时在杀菌釜内放置1 个温度传感器。随后杀菌釜开始进行保温温度121 ℃-保温时间30 min 的自动杀菌。冷却完成后，取下4 个温度传感器，读取杀菌锅温度及罐内液体预估冷点温度。将温度探针固定于距底部100 mm 的位置，进行相同试验，获取罐头顶隙处温度变化曲线。

图1 热渗透试验装置示意图Fig.1 Schematic diagram of thermal penetration experimental device

2）1.0% CMC 溶液的模型验证：配制质量分数为1.0%的CMC 溶液，用同样方法测定1.0%CMC 罐头液体及顶隙处温度随时间变化的曲线。

计算试验所得实际温度与数值模型的模拟温度的平均相对误差，以评价模型的准确性。

2 数值模型的搭建

在COMSOL Multiphysics 5.3a 软件中构建了2 种二维轴对称模型，对液体罐头食品的热杀菌过程进行模拟。模型一忽略了顶隙的影响，认为罐内被液体充满，而模型二液面高度为95 mm，顶部为湿空气，同时考虑湿空气与液体表面水分的蒸发冷凝。利用层流、湿空气中的传热、空气中的水分输送3 个物理场模拟传热传质过程，并添加了热湿、水分输送、非等温流动节点来实现传热与蒸汽输送之间的多物理场耦合。

2.1 模型假设

对2 种模型均做以下假设以简化求解过程：罐体内液体、气体均匀对称；忽略罐壁、探针及固定装置等对传热的影响；液体在罐内壁无滑移；加热过程中罐体不发生形变；罐体外壁的温度分布均匀且等于杀菌釜的温度。

2.2 控制方程

用Navier-Stokes 方程模拟空气域中的流体流动并用动量守恒方程、能量守恒方程模拟相应的对流传热[17]。由于涉及水分蒸发，需将蒸发的水蒸气质量作为潜热源添加到输运方程中[18]。蒸发的水蒸气质量按式（1）计算：

式中，mvap——蒸发的水蒸气质量，kg；K——蒸发率，s-1；aw——水分活度；csat——饱和水分浓度，mol/L；c——水分浓度，mol/L。

式中，psat——饱和压力，Pa；R——理想气体常数，8.314 J/（mol·K）；T——温度，K。

将蒸发热作为潜热源插入到传热方程中：

式中，Q——蒸发热，J；Hvap——蒸发潜热，44 172 J/mol。

2.3 计算条件

初始温度为室温，初始绝对压力（-40 kPa）作为罐内的真空度。纯水的热物理特性参数由COMSOL Multiphysics 内置材料库给出，1.0%CMC 溶液的参数如表1[19]。

表1 1.0% CMC 溶液的热物理参数特性Table 1 Thermal physical properties of 1.0% CMC solution

2.4 网格剖分及求解

2 种模型均利用细化的物理场控制的自由三角形网格进行剖分，模型一与模型二分别将几何体剖分为5 429 及5 971 个单元进行求解。

在模型验证试验中，将实际测量得到的保温温度121 ℃，保温时间为30 min 的杀菌釜内温度作为内插函数，导入软件中作为罐头外壁的温度函数，以表征罐体在杀菌釜内的受热情况。对液面高度为95 mm 的水及1.0% CMC 罐头进行仿真求解及模型验证后，通过改变温度函数以适应不同的杀菌工艺。利用参数化扫描，将液面高度由75 mm 以5 mm 为间隔增加至100 mm，保温温度由85 ℃以5 ℃为间隔递增至130 ℃进行求解。

表2 模拟值与试验值的平均相对误差Table 2 Average relative error between simulated and experimental values

3 结果与讨论

3.1 模型验证

图2比较了水和1.0% CMC 溶液在罐头中心高度为15 mm（内容物）及高度为100 mm 处（顶隙）的试验温度变化及不同模型在相同位置的温度模拟值，表2显示了模拟值与试验值之间的平均相对误差，其误差均在5%以下。由于实际生产过程中，罐头内容物主要为黏性液体，故本研究主要对1.0% CMC 的模拟结果进行分析。

由图2及表2可看出，模型二对罐头顶隙处的温度变化尤其是冷却过程描述得更为准确，说明考虑顶隙的模型能更真实地反映热杀菌过程中罐头内的传热情况，精准地分析罐头内物理量的变化，准确地预测罐内冷点位置及能达到的杀菌程度。

图2 实际测量及不同模型仿真得到的罐内冷点及罐头顶隙处温度变化曲线Fig.2 Temperature variation curve of the cold spot and the head gap of the can obtained by different model simulations and experimental measurements

3.2 物理场分析

3.2.1 速度场 选取保温温度131 ℃，升温时间10 min，保温时间15 min，冷却时间10 min 的理想杀菌曲线，对加热条件下热杀菌过程中的物理场进行分析。图3～4 为2 种模型在杀菌全过程中不同阶段的速度场模拟结果，用处处都与速度矢量方向相切的流线及面上箭头来表征速度大小，流线分布越密，箭头越长表明该处速度越快。

模型一仅考虑液体的非等温流动引起的速度场变化。在升温阶段，靠外壁处的液体升温较快，密度降低而上浮，罐内温度相对较低的液体下沉。由于液体和罐头底面的碰撞产生挤压，在靠近底部中心位置出现速度向上的二级流场。冷却过程的流场分布与升温过程类似，二级流场出现在顶部靠近中心轴线处。保温阶段的液体流动速度较低。

模型二模拟结果的液体部分流线分布情况与模型一的模拟结果类似。升温阶段，顶隙处的湿空气与底部液体大致呈逆时针流动，降温阶段与之相反，且存在二级流场。与底部液体的流动场相比，顶隙处流线较密，箭头较长，可能是由于气体密度低，动力黏度小，流动速度远大于底部液体。同时，受水分蒸发冷凝的影响，液面处流线较密，速度较大，且存在垂直液面的速度分量。

图3 模型一模拟的速度场分布图Fig.3 Velocity field distribution simulated by model 1

图4 模型二模拟的速度场分布图Fig.4 Velocity field distribution simulated by model 2

3.2.2 温度场 液态罐头食品的传热方式以对流传热为主，流体速度场会对温度分布产生影响。图5为2 种模型模拟生成的131 ℃-15 min 的杀菌过程中不同阶段的温度场分布云图。受自然对流传热影响，罐头顶部升温较快，保温后期SHZ 位于离底部约10%罐高，离罐壁10%～15%半径处而非罐头中心位置。对比模型一与模型二得到的温度场分布情况，发现由于顶部湿空气的比热与密度较小，传热速度更快，使得模型二的顶隙部分升温更快，并由此导致靠液体更快升温。同时水分蒸发吸热也会对液面附近的温度分布造成影响，在水蒸气未饱和时靠近液面处的空气温度较低。

图5 不同模型模拟的温度场分布图Fig.5 Temperature distribution simulated by different model

3.2.3 致死率场 工业生产中通常用累计致死率评价特定杀菌条件下的总体杀菌效果，致死率按式（4）计算：

以肉毒杆菌为参考，式中，F——致死率，min；T——食品冷点的温度，℃；Tref——杀菌标准温度，121℃；Z——微生物耐热性特征值，10 ℃。

2 种模型温度场的差别最终会对致死率场产生影响，图6为131℃-15 min 的杀菌过程完成后，罐内的微生物致死率场分布情况，表3为罐内内容物的最大、最小致死率值及最小致死率所在的坐标。结果表明，2 种模型模拟得到的致死率最小值不同，且出现的位置有所差别。模型二模拟的致死率值总体较模型一高，分布更为均匀。随着保温温度的升高及罐头几何尺寸的增大，不同模型对致死率分布的模拟差别更为显著，说明建立模型时考虑顶隙与否会影响罐头致死率场的分布。为更准确地寻找罐内冷点位置，更精准地确定最佳保温温度及最短杀菌时间，需应用含顶隙的模型。

3.2.4 压力 为保证罐头内外压力差小于罐头可承受的临界压力差，避免出现凸罐或瘪听等罐头变形甚至破损等问题，实际热杀菌过程中需对压力进行控制，同时根据罐内外压力差选择合适的罐壁厚度以确保罐头的耐压值大于杀菌过程的最大压力差。由于模型一不考虑顶隙，无法模拟热杀菌过程中罐内的压力变化，只能通过理想气体方程，利用温度变化推算罐内的压力变化情况。图7为121 ℃-30 min 灭菌过程中，通过模型二模拟和通过罐内温度推算的罐内相对压力及杀菌釜压力表读数，结果表明由温度计算罐内的相对压力不准确，而模型二模拟得到的罐内绝对压力与实际情况较为接近。对该杀菌过程而言，罐内外压差最大值出现在刚开始进行补压时，差值为0.025 MPa，需保证罐头的耐压值大于该值。

图6 不同模型模拟杀菌结束时微生物致死率场的分布图Fig.6 Microbial lethality field distribution at the end of sterilization simulated by different models

工业生产中以纯水的饱和蒸气压为参考对杀菌釜进行压力控制，然而由于分压效果及水分活度的影响，完全用饱和蒸气压控制压力并不够准确。利用模型二模拟得到不同保温温度及不同顶隙高度下保温阶段罐头内的绝对压力及相应温度下的饱和水蒸气压。由表4可知，随着温度的升高，顶隙高度的减小，罐内的相对压力逐渐增大，这与理想气体方程相吻合。顶隙过小，罐内压力大，容易产生凸角、假胀等现象，而顶隙过大则可能在冷却过程中出现瘪罐。生产厂家可通过模型模拟辅助设定不同顶隙高度下杀菌工艺中各阶段的压力，以避免出现胀罐、瘪听等罐头外观变形问题。

表3 两种模型模拟的致死率最大值及最小致死点坐标Table 3 Maximum and minimumt lethality and minimum lethal coordinates of the two models

图7 罐内相对压力模拟值及杀菌釜相对压力变化曲线Fig.7 Simulated relative pressure curve in cans and relative pressure curve of sterilizer

4 结论

1）利用COMSOL Multiphysics 软件构建的2种二维轴对称模型，都能较为准确地模拟水、1.0%CMC 溶液的热杀菌过程。其中模型二对该过程模拟的温度变化过程与实测值更为接近，是一种能精准分析物理场变化的仿真模型。

表4 不同保温温度及不同顶隙高度下罐内压力（MPa）与相应温度下饱和水蒸气压Table 4 Pressure （MPa）in the cans at different holding temperatures and different headspace heights and saturated water vapour pressure at corresponding temperature

2）为更准确地预测罐头内速度场、温度场及微生物致死率场的分布情况，不可忽略顶隙的存在。本文建立了更准确地分析杀菌过程中各物理场的变化的模型，适用于不同罐形大小及不同内容物的液态罐头食品，有助于在热渗透测试进行前准确寻找冷点位置，最终实现杀菌工艺的优化。

3）含顶隙模型可对罐内压力变化情况进行分析，本文给出不同顶隙高度及不同保温温度下的罐内相对压力，顶隙高度越小、温度越高，罐内压力越大。该模型可用于多种罐形及不同热杀菌工艺液态罐头食品物理场的预测，以便生产厂家对杀菌釜压力变化情况及罐头最大耐压值进行控制，避免出现胀罐或瘪罐等质量问题。