王金锋 李文俊 谢 晶

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



数值模拟在食品冻结过程中的应用

王金锋1,2李文俊1,2谢 晶1,2

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

文章主要通过综述冻结过程的数值模型、求解微分方程、预测冻结时间及分析送风速度、温度和送风方式等，讨论数值模拟技术在食品冻结过程中的应用现状；总结了国内外研究者针对不同冻结对象所采用的数值模拟方法，为今后数值模拟方法在食品传热过程中进一步发挥作用提供理论参考。

冻结过程；速冻食品；数值模拟；食品冷加工

随着速冻产业在中国的飞速发展，中国自主研发并生产的速冻设备和食品加工技术已有了明显提高，但与发达国家相比，还存在一定的差距。

食品的冻结速度直接决定了速冻食品的质量，因此，采用科学的研究方法对冻结过程中食品温度的变化进行准确预测，可以有效地预测冻结时间，并且对指导速冻设备优化设计、提高设备能效、降低人员消耗具有重要意义。然而，大多数食品存在非均匀的物性参数，以及冻结过程存在相变导致食品冻结过程是一个严格非线性变化过程，因此很难找到一个普遍适用的方法对冻结过程进行描述[1]。

通常预测模型采用：理论求解、经验公式和数值求解三种方法[2]。其中理论求解和经验公式法通常只用于求解简单的研究对象，对复杂的研究对象和流场环境，前两种方法并不能进行有效地描述。比较前两种求解方法，数值求解能够结合食品的热物性及流场环境进行分析求解，且不受模拟对象的限制，尤其在模拟复杂的流场环境时优势更为突出。本文主要综述了数值模拟技术在食品冻结过程中的实施方法及其应用情况，以期为今后该技术在食品传热过程研究中进一步发挥作用供理论参考。

1 食品冻结过程分析

食品冻结是通过使用制冷装置将食品中游离态的水分子降温至结晶形成冰的过程[3]。食品的形状、尺寸、热导率、比热容，冷却介质的温度、速度，气流组织的均匀性等都是影响食品冻结速率的重要因素。在食品工业中，冻结过程不仅影响食品的质量，而且与节能息息相关，因此对冻结过程进行模拟分析是研发高质高效冻结装置的重要前提。

冻结过程的热力学分析属于低温传热领域，而相变传热和食品随时间变化的热物性在一般的传热过程中往往被忽略，但在低温传热领域，上述两方面是不可忽略的问题[1]。食品冻结过程中，热量由内而外通过导热的形式传递到食品外表面。因此，预冷过程中食品内部存在较大的温差，且越靠近食物中央，温度梯度越大；随着食品热量散失，食品最外层最先形成一层冻结层，并且随着时间推移，冻结层向食品中心移动，最终完全冻结。由于冻结过程中涉及相变，而相变传热是严格非线性过程[1]，普通的试验测温法难以对冻结过程实现直观的数学描述，因此数值模拟技术已成为食品冻结科研工作者有效的科研工具。

2 数值模拟技术在食品冻结过程中的应用

20世纪60年代以后，随着计算机技术的迅猛发展，速冻工业也迅速发展起来。大多数研究人员以禽肉、牛肉、水产作为研究对象，并根据食品的几何特征建立一维、二维、三维物理模型(表1)。在食品冻结过程中对试验对象的边界条件、导热系数、热质转移进行合理的假设，通过数值模拟对速冻时间、送风速度和温度及送风方式等进行预测和优化。 由于相变传热过程完全非线性，其过程十分复杂，通过分析进行精确求解这一问题是不可能的。因此，为了提高分析结果的准确性，采用数值模拟法来近似求解冻结过程。过去的研究工作中，无限大平板、无限长圆柱和球形的试验对象常常按一维问题处理[2]；然而，实际情况下物体的形状都十分复杂，因此一维模型并不能普遍用来解决实际问题；目前科究人员针对不同的物体形状，采用二维或三维模型对实际问题求解[4-5,9]。

表1 数值模拟在食品冻结过程中的应用

2.1 建立物理模型

2.1.1 创建几何模型 对研究对象进行几何建模，是将对象进行数学描述的过程。常规几何模型有以下三种：无限大平板、无限长圆柱和球体[2]；当研究对象的几何形状和三种基本模型类似时，求解过程中加入合适的几何形状因子，近似地按标准形状处理，这样不仅可以节约大量计算时间，且收敛精度较高。然而，在实际问题中，多数研究对象几何外形复杂，采用近似处理的方法误差较大。近年来，一些科研工作者开始把目光投向图形处理，如李宝方等[10]通过数码设备对研究对象进行图像采集，使用MATLAB对图像进行分析处理，将图像划分为二维或三维离散区域，再导入网格生成软件划分网格、设置边界条件。Goni等[11-12]使用图像处理软件对食品表面进行色区分割，然后通过后处理软件对分割模型进行形态调整，并构建食品几何模型，最终结果表明：通过线性差值构建的食品几何模型比真实样品表现出良好的一致性，该建模方法可以有效地减小数值模拟过程中采用简化几何模型所带来的误差。

2.1.2 创建数值模型 数值模型建立的准确性，将直接影响试验结果和模拟值之间的接近程度。在模拟计算过程中，首先需要通过网格生成软件将所建立的食品几何模型离散化，为下一步求解微分方程提供前提；通常商业软件所采用的离散方法有：有限差分法、有限元法和有限体积法[2](表2)。

有限差分法是数值模拟计算中最早采用的网格划分方法，其基本思想在于，通过对模拟区域进行网格划分使连续区域离散化，并通过网格上的节点对连续区域上的连续变量函数使用离散后网格上定义的离散变量函数来近似替代。考虑时间因子在计算过程中的作用，可将有限差分法简单细化为：显式、隐式两种格式。采用显式格式对研究对象进行数值计算时，其计算工作量较小，但对时间、空间步长有一定的限制，否则将出现震荡解；隐式格式对步长没有限制，不会出现震荡现象，但计算工作量大[13]。有限差分法在数值计算过程中是逐层推进的，误差也会逐层累积，因此有限差分法使用于较为复杂的模拟区域中时，最终计算结果误差偏大，同时对计算过程中收敛精度要求较高，导致计算时间花费较大，不利于对实际问题的求解；因此，通常不采用有限差分法建立数值模型。但具体在冷冻传热过程中，Botheju等[6]采用有限差分法建立一维数学模型，不仅极大地缩短了计算时间，而且将标准误差控制在0.63%～1.23%，使用有限差分法对特殊研究对象创建简单的物理模型，计算结果令人满意。

表2 数值离散法分类原理及优劣

Table 2 The classification principle of numerical discretization method and its advantages and disadvantages

数值离散方法原理及优劣有限差分法将连续区域离散化后,通过网格上的节点对连续区域上的连续变量函数使用离散后网格上定义的离散变量函数来近似替代。求解速度快、精度较低有限元法 将计算域转换为一系列的子域,通过每个子域内的近似函数来描述所求解的未知场函数。求解速度慢、精度高有限体积法计算区域内所包含的每个节点的控制体都分布于网格上,通过对控制体进行积分来求解微分方程。求解速度快、精度较高

有限元法其基本思想是将计算域转换为一系列的子域，通过每个子域内的近似函数来描述所求解的未知场函数。有限元法不受空间的限制，可以实现对任意模型进行预测，尤其在食品以导热的形式进行热量传递时，有限元法表现出优良的精度[1]。因此，有限元法被Ansys公司广泛地应用在数值计算模块。Sabliov等[4]对一个轴对称非稳态传热的有限元在低温冷却中开发了壳蛋低温冷却模型，并通过对预测模型输出温度与试验测量温度进行比对；结果表明，有限元法可以有效地预测相变过程热物性随时间的变化规律，从而精确地预测冷却时间。Pham等[14]以牛肉为研究对象建立有限元模型，将基于二维网格上的热传导方程和一维网格上的水分扩散方程总结成一组回归方程，对食品的各项物性参数进行预测，最终通过实验得到了很好的验证；但由于有限元法求解过程复杂，并且对计算机的内存要求高，因此并未得到推广，而仅适用于一些研究精度较高的工程领域。

有限体积法是商用CFD软件最常使用的方法，其与有限差分法和有限元法的本质区别在于，采用有限体积法离散的计算区域内所包含的每个节点的控制体都分布于网格上，通过对控制体进行积分来求解微分方程。Moraga等[5]通过以三种不同的肉柱为研究对象，在研究冻结过程中比较了有限差分法和有限体积法对二维定常条件下时间及热物性的预测，结果表明有限体积法更能精确地模拟冻结过程以及估计该产品冻结时间。

冻结过程中由于食品存在热物性，将导致冻结结束后食品形态发生剧烈变化；比如，相变过程伴随着膨胀、收缩，现有的网格技术并不能解决物性变化对几何模型所产生的影响。因此，近几年一些学者如：Belytschko等[15-17]开始将目光投向无网格技术，并指出无网格方法和空间区域离散法的本质区别在于，无网格法将模拟对象近似成点的集合，此种方法将有效处理冻结过程中由于食品物性的突变所产生的大变形以及裂纹等问题，在模拟过程中不需要对网格进行重构和划分，避免了网格设置等技术难题，不仅提高了计算结果的准确性，而且降低了计算难度。目前，通过各国科学家的努力已提出了十余种无网格法，如：伽辽金法[18-20]、最小二乘法[21-22]、配点法[21]等。Amantea等[23]开发了一种基于三次径向基函数的无网格数值方法求解非线性耦合对流占优的偏微分方程组；将此方程运用于干燥过程中，模拟结果与前人的研究相比表现出优良的精度。目前，无网格技术在食品干燥领域已初步涉及，而在冷冻领域还没有可查阅的相关文献资料[24]。

2.2 求解微分方程

使用数值模拟计算是在质量、能量守恒方程满足的前提下进行的，求解过程中应充分考虑导热和扩散对能量的影响。Welti等[25]认为由于食品冻结过程涉及到相变，其物理性质(密度、比热容、热导率等)将随温度的变化发生突变，因此求解过程中需要对模拟对象建立间断微分方程，通过将物理模型连续区域离散化，求解微分方程得到仿真解，仿真解和试验结果进行对比即可评价建立模型的优劣。当计算区域较复杂时，有限差分法性能较差，长期研究发现，选用合适的时间步长可以改善这种情况。食品在冻结过程中由于一些物理性质发生突变，其非线性使得数值计算难度加大，且准确性较低，基于这一情况，基尔霍夫转换将所有的非线性因素都统筹为一个影响因子，使食品中不同温度节点取值得以解决。Scheerlinck等[26]发现在数值模拟计算过程中采用基尔霍夫转换时可以节省大量的迭代时间。此外，还有表观比热容法、热质传递耦合以及准焓法；其中表观比热容法是商业模拟软件常采用的方法，此方法特征在于：将食品冻结过程中相变所释放的潜热以显热的形式在传热方程进行表达，但由于潜热的叠加使表观比热与温度呈现非线性关系，导致食品冻结点附近的比热容难以确定[2]。Santos等[27]使用DSC测得食品比热，并将其作为模型参量输入数值计算模型，解决了冻结点比热容难以确定的问题。

2.3 预测冻结时间

20世纪60年代以来，数值模拟技术在冻结领域不断发展，其中最重要的任务就是预测冻结时间，因为冻结时间直接关系到食品品质和能耗。Norton等[28-30]提出了对食品冻结时间预测的方法，并指出冻结时间预测公式的验证需要高质量的试验数据，然而研究表明，常用的试验方法有不足之处，特别是在维持表面换热系数的均匀性、试验过程中引入的误差等，都将影响预测结果的准确性。李杰等[9]对鼓风式速冻装置中的食品冻结过程以及冷冻时间进行了非稳态模拟，结果表明，食品在冻结过程中预测值与试验值之间绝对温度误差最大值为1.5 ℃，食品冻结时间相对误差为3.8%，从试验角度验证了数值模拟计算的可行性。针对不同的冻结阶段，李杰等[31]以冰箱为试验平台，采用各项物性较为均匀的土豆为试验对象，建立三维数值模型，但由于冻结模型只适用于液体，因此数值模型设置时需要将土豆四周设置成固定壁面的封闭体，且内部设置成和土豆物性参数(热导率、密度、比热容)相同的液体，同时不考虑非凝固区域的传质与流动，通过模拟土豆在冷冻过程中温度场随时间的变化规律，较为准确地预测了土豆的冻结时间。

2.4 分析送风速度、温度及送风方式

冻结装置的送风速度和温度直接影响食品的冻结时间，对于最佳送风速度的研究，过去已有一些学者[32-34]进行了相关的试验。Dominguez等[34]以馒头为试验对象，对流场内低温空气的温度和速度进行了研究，结果表明：食品的冻结速度与空气冷却器进口风速成正比。Song等[7]以马铃薯为研究对象，使用COMSOL仿真软件对马铃薯建立二维瞬态导热固体模型，并分别以1，2，4，8 m/s的风速进行模拟计算，将仿真结果和试验结果对比分析，结果表明，当冷却空气为-40 ℃时，随着风速的增加，相变表面移动速度降低，当无冷却风时，背风面相变界面移动最快，当风速大于4 m/s时，迎风面的相变表面移动速度大于背风面。张珍等[35-37]通过对送风方式和流场区域进行优化，将传统的单侧吹风方式变为上下两侧吹风，并且研究了不同网带孔隙率对食品的冻结速度的影响，结果表明：上下两侧协同送风的送风方式能更快地加速食品的冻结。随着研究的深入，李杰等[37]以虾为试验对象建立了二维数值模型，试验结果与模拟结果存在较高的一致性，在此基础上还重点讨论了影响虾冻结时间的多个物理参数，结果表明：上下两侧吹风与单侧吹风相比能够有效地缩短食品冻结时间，并且当上吹风风速一定时，提高下吹风风速能够明显提高冻结速率，同时使冻结过程中食品内部温度均匀降低，有效地提高了食品冻结质量。Sarkar等[38-39]对鼓风式速冻设备进行了一系列的设计研究，模拟了不同喷嘴直径与喷嘴长度比值α，不同的冲击方式对食品冻结过程效率的影响。研究结果表明：当α在6～8时效率最高，双向喷射比单向喷射效率高。

3 结论与展望

使用数值模拟技术预测食品冻结过程中的各项物理参数是目前普遍使用的科研工具，但由于食品存在复杂的几何形状以及食品冻结过程中多变的热物理性质，这些因素将直接影响数值模拟结果的准确性。针对食品种类与外形的多样性，准确地选择研究工具有助于降低软件本身带来的误差。总结前人的研究可以发现：当食品通过导热的形式进行冻结，有限元分析法可以有效地和试验结果吻合；当食品以对流传热的形式传热，有限体积法预测结果令人满意；当食品几何形状复杂时，使用图像处理分析法建立几何模型可以直观地对模拟对象进行准确地描述。

将来对冻结食品采用数值模拟方法的研究工作可以从以下几个方面展开：① 运用图像处理技术对研究对象建立更准确的物理模型；② 根据送风条件，对冻结过程进行优化设计；③ 开发更为精确的数值模型方法(软件)，从技术角度降低预测结果与试验之间误差；④ 总结前人对同类研究对象采用的不同模拟方法，分析预测结果筛选误差较小的模拟方法作为解决此类对象的最优解法，为后续研究工作提供参考依据。⑤ 借鉴无网格技术在干燥领域的应用，将无网格技术应用于冻结过程的分析计算，开拓一种研究冻结过程的新工具。

[1] HUAN Zhong-jie, HE Shao-shu, MA Yi-tai. Numerical simulation and analysis for quick-frozen food processing[J]. Journal of Food Engineering, 2004, 60(3): 267-273.

[2] 成芳, 杨小梅, 由昭红, 等. 食品冻结过程的数值模拟技术[J]. 农业机械学报, 2014(7): 162-170.

[3] 赵金红, 胡锐, 刘冰, 等. 几种冷冻新技术对食品冻结过程中冰晶形成的影响[J]. 食品与机械, 2012, 28(6): 241-245.

[4] SABLIOV C M, FARKAS B E, KEENER K M, et al. Cooling of shell eggs with cryogenic carbon dioxide: a finite element analysis of heat transfer[J]. LWT-Food Science and Technology, 2002, 35(7): 568-574.

[5] MORAGA N O, VEGA-GALVEZ A, LEMUS-MONDACA R. Numerical simulation of experimental freezing process of ground meat cylinders[J]. International Journal of Food Engineering, 2012, 7(6): 639-646.

[6] BOTHEJU W S, AMARATHUNGE K S P, ABEYSINGHE I S B, et al. Modeling trough withering system to predict the moisture content of tea leaves at real time using one dimensional heat and mass transfer finite difference model[J]. Journal of Tea Science, 2010, 75(1): 27-41.

[7] SONG Xiao-yan, ZOU Tong-hua, HONG Qiao-di, et al. Effect of different positions on quick-frozen results of bean buns and steamed bread in quick-freezer[J]. Storage and Process, 2014, 14(6): 40-42.

[8] CLELAND D J, CLELAND A C, JONES R S. Collection of accurate experimental data for testing the performance of simple methods for food freezing time prediction[J]. Journal of Food Process Engineering, 1994, 17(1): 93-119.

[9] 李杰, 谢晶. 上下冲击式速冻装置内静压箱流场的均匀性研究[C]// 第六届全国食品冷藏链大会论文集. 上海: 中国制冷学会、全国商业冷藏科技情报站, 2008: 6.

[10] 李宝方, 邢玉清. 基于MATLAB GUI的数字图像处理仿真系统设计[J]. 计算机与数字工程, 2014(11): 2 177-2 180.

[11] GONI S M, PURLIS E, SALVADORI V O. Three-dimensional reconstruction of irregular foodstuffs[J]. Journal of Food Engineering, 2007, 82(4): 536-547.

[12] GONI S M, PURLIS E, SALVADORI V O. Geometry modelling of food materials from magnetic resonance imaging[J]. Journal of Food Engineering, 2008, 88(4): 561-567.

[13] 杨世铭, 陶文铨. 传热学[M]. 4版. 北京: 高等教育出版社, 2006: 176-179.

[14] PHAM Q T, TRUJILLO F J, MCPHAIL N. Finite element model for beef chilling using CFD-generated heat transfer coefficients[J]. International Journal of Refrigeration, 2009, 32(1): 102-113.

[15] BELYTSCHKO T, KRONGAUZ Y, ORGNA D, et al. Meshless methods: an overview and recnet developments[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1):3-47.

[16] 张雄, 刘岩. 无网格法[M]. 北京: 清华大学出版社, 2004: 1-13.

[17] 张雄, 宋康祖, 陆明万. 无网格法研究进展及其应用[J]. 计算力学学报, 2003, 20(6): 730-742.

[18] ATLURI S N, SHEN S P. The meshless local Petorv-Galerkin method[M]. [S. l.]: Tech Science Press, 2002: 93-124.

[19] ATLURI S N. The meshless local Petorv-Galerkin(MLPG) method for domain & boundary discretizations[M]. [S. l.]: Tech Science Press, 2004: 105-145.

[20] 张俊贤, 朱风风, 王金田. 题解有限元法和无网格伽辽金法[J].山西建筑, 2010, 36(1): 241-245.

[21] 陶文铨, 吴学红, 戴艳俊. 无网格数值求解方法[J]. 中国电机工程学报, 2010, 30(5): 1-10.

[22] 张雄, 胡炜, 潘小飞. 加权最小二乘无网格法[J]. 力学学报, 2003, 35(4): 425-431.

[23] AMANTEA R P, FORTES M, MARTINS J H, et al. Numerical simulation techniques for optimizing thermodynamic efficiencies of cereal grain dryers[J]. Drying Technology, 2013, 31(6): 672-683.

[24] KARUNASENAA H C P, BROWNA R J, GUA Y T, et al. Application of meshfree methods to numerically simulate microscale deformations of different plant food materials during drying[J]. Journal of Food Engineering, 2015, 146: 209-226.

[25] WELTI-CHANES J, VERGARA-BALDERAS F, BERMUDEZ-AGUIRRE D. Transport phenomena in food engineering: basic concepts and advances[J]. Journal of Food Engineering, 2005, 67(1/2): 113-128.

[26] SCHEERLINCK N, VERBOVEN P, FIKIIN K A, et al. Finite element computation of unsteady phase change heat transfer during freezing or thawing of food using a combined enthalpy and Kirchhoff transform method[J]. Transactions of the Asae American Society of Agricultural Engineers, 2001, 44(2): 429-438.

[27] SANTOS M V, LESPINARD A R. Numerical simulation of mushrooms during freezing using the FEM and an enthalpy: Kirchh off formulation[J]. Heat and Mass Transfer, 2011, 47(12): 1 671-1 683.

[28] TOMAS N, SUN Da-wen. Computational fluid dynamics (CFD) - an effective and efficient design and analysis tool for the food industry: a review[J]. Trends in Food Science & Technology, 2006, 17(11): 600-620.

[29] CLELAND D J, CLELAND A C, JONES R S. Collection of accurate experimental data for testing the performance of simple methods for food freezing time prediction[J]. Journal of Food Process Engineering, 1994, 17(1): 93-119.

[31] 李杰, 谢晶, 陆方娟. 食品冻结过程温度场及冻结时间的数值模拟与实验研究[J]. 食品工业科技, 2009(2): 123-125.

[32] KOTCHETKOV N D, MIKHAILIN N V, AVERIN G D. The effect of air temperature and velocity at various periods upon the technological characteristics during meat freezing[J]. Proceedings of the International Congress of Refrigeration (13th Washington), 1973(3): 205-213.

[33] MOLEERATANOND W. Heat exchange of boxed ground beef as influenced by packaging, refrigeration temperatures and air velocities[J]. Dissertation Abstracts International, 1982, 42(8): 148-178.

[34] DOMINGUEZ M, ELVIRA C, FUSTER C. Influence of air velocity and temperature on the two-stage cooling of perishable large-sized products[J]. Bulletin de l'Institut International du Froid, 1975(4): 83-90.

[35] 张珍, 谢晶. 带有上下均风孔板的速冻装置中流场及温度场的数值模拟[J]. 制冷学报, 2009, 30(5): 36-40.

[36] 张珍, 谢晶. 上下冲击式高效鼓风冻结装置速度场的数值模拟与验证[J]. 低温工程, 2008(6): 45-50.

[37] 李杰, 谢晶. 鼓风冻结虾仁时间的数值模拟及实验验证[J]. 农业工程学报, 2009, 25(4): 248-252.

[38] SARKAR A, SINGH R P. Air impingement technology for food processing: visualization studies[J]. Lebensmittel-Wissenschaft und-Technologie, 2004, 37(8): 873-879.

[39] HU Ze-hua, SUN Da-wen. Modelling of an experimental air-blast freezer using CFD code, in: advance in refrigeration system, food technologies and cold chain[J]. International Institute of Refrigeration, 1998, 6: 395-400.

Application of numerical simulation in process of food freezing

WANGJin-feng1,2LIWen-jun1,2XIEJing1,2

(1.ShanghaiEngineeringResearchCenterofAquaticProductProcessing&Preservation,Shanghai201306,China;2.FoodCollegeofShanghaiOceanUniversity,Shanghai201306,China)

The freezing process of numerical model, solving differential equations, predicting freezing time and analyzing of air supply velocity, temperature and methods were reviewed. The application of numerical simulation technology in food freezing process was discussed. The numerical simulation methods for different freezing objects were also summarized, which would be useful for the future application of numerical simulation methods in the process of heat transfer of foods.

freezing process; quick-frozen food; numerical simulation; food cold processing

2014年国家农业成果转化资金项目(编号：2014GB2C000081)；2016年上海市科技兴农重点攻关项目[编号：沪农科攻字(2016)第1-1号]；上海市科委平台能力提升项目(编号：16DZ2280300)

王金锋，女，上海海洋大学副教授，博士。

谢晶(1968—)，女，上海海洋大学教授，博士，博士生导师。E-mail: jxie@shou.edu.cn

2016-05-24

10.13652/j.issn.1003-5788.2016.10.044