马铃薯超声解冻过程的数值模拟研究

2019-09-04刘东红吕瑞玲王丹丽王文骏周建伟

中国食品学报 2019年8期

李慢刘东红，2* 吕瑞玲王丹丽丁甜王文骏周建伟

（1 浙江大学生物系统工程与食品科学学院杭州310058 2 浙江大学馥莉食品研究院杭州310058）

我国食品速冻工艺有着较快的发展，相比于其它贮藏方法，速冻能很好地保持其原有的风味、色泽和营养价值[1]。而解冻工艺对食品最终品质也有一定的影响。不当的解冻方式会增大汁液流失，加速品质的劣变等一系列不利的影响。如果能精确掌握解冻食品内部的温度分布和变化趋势以及解冻时间，这将对优化解冻工艺，加快生产效率，提高食品品质有重要的意义。

数值法结合计算机知识对一些现象或者过程进行模拟的技术可简称为数值模拟技术[2]，它在食品冻结和解冻过程中的应用非常广。Birla 等[3]用FEMLAB 软件来模拟水果的射频解冻过程，设计试验验证结了冷凝胶水果模型解冻过程中的温度分布，结果证明模拟结果与试验结果十分吻合，该模型能成功预测解冻过程中的瞬时温度分布。Kiani 等[4]用CFD 软件来模拟超声辅助土豆冷冻的过程，基于有限体积的热焓法来进行数值模拟，其中冷冻时间的预测值和真实值的误差最大值不超过0.05，中心温度的模拟值和试验值也较为一致。Pitchai 等[5]用COMSOL 模拟在微波作用下土豆泥的解冻情况，模拟值和试验值拟合较好，REMS 为7.5 ℃。Uyar 等[6]用COMSOL 来模拟射频解冻牛肉，模拟值和试验结果的REMS 值仅为0.9℃，表明该模型的精确性和可靠性。

早期对食品解冻的模拟大多假设食品热物性在解冻过程中为常数，而实际上食品的热物性在解冻过程中有较大的变化，特别是含水量高的食品。食品解冻发生在一定的温度范围内而非某一特定温度。由于解冻涉及食品中水的相变过程，需要吸收大量的相变潜热，同时食品的密度、比热容和热导率等热物理性质会随着解冻的进行不断发生变化，因此本文采用显热容法来解决潜热问题，将冰晶融化需要的潜热表达在比热容上，采用变热物性方程来反映马铃薯在解冻过程中热物性参数随温度的变化[7]。

现阶段对解冻过程中的传热模拟已有很多研究，大都集中在对射频解冻和微波解冻的模拟[3，5-6]，而对超声解冻这种新型的解冻方法的有限元模拟几乎没有。本文以马铃薯为对象，应用COMSOL Multiphysics 5.2a 模拟超声解冻马铃薯的过程，建立食品的热物理学性质随时间变化的超声－热传导双向耦合模型，并对模拟结果进行试验验证。

1 数学物理模型

1.1 超声解冻装置与模型

解冻装置如图1所示。超声解冻池是半径为51 mm，高为150 mm 的圆柱状水池，其液面高度为120 mm，解冻介质是纯净水。待解冻的土豆是半径为17 mm，高为25 mm 的圆柱体，放在高度为25 mm 的镂空金属网支架上，中心插入data trace温度探头，用来记录中心温度的变化。

建立模型时，选用固体传热和压力声场两个物理场，根据实验装置建立几何模型，相关参数如表1所示。其中解冻介质水的基本参数由COMSOL Multiphysics 软件提供。整个装置均匀对称，为了简化计算，将解冻装置视为二维轴对称的声热场。为了简化研究，建立模型时，有几点假设：

1）马铃薯在解冻过程中与外界环境仅有热交换；

2）马铃薯各组成部分均为各向同性的均匀体；

3）马铃薯在解冻过程中，内部热交换只考虑导热，不考虑对流传热；

图1 超声辅助解冻装置示意图Fig.1 Schematic drawing of the experimental setup of ultrasound-assisted thawing

4）解冻开始时，马铃薯温度和水分分布均匀。

表1 相关参数Table 1 Relevant parameters

1.2 理论基础

超声解冻是多物理场问题，本文研究传热场和声场的耦合。

在二维轴对称坐标体系中波动方程采用霍姆赫兹方程：

式中：P——声压（Pa）；r、z——分别是径向坐标和轴向坐标（mm）；ρ——密度（kg/m3）；c——物料的声速（m/s）；ω——角频率（rad/s）。

传热控制方程为：

式中：T——马铃薯的温度（℃）；t——时间（s）；ρ——马铃薯的密度（kg/m3）；Cp——马铃薯的比热容（J/（kg·K）；k——马铃薯的热导率[W/（m·K）]；Q0——单位体积物料所吸收的超声能量（W/m3），可由式（3）式计算：

其中声压P 和声强I 分别是：

由式（4）、（5）得：

式中：α——超声在物料中的衰减系数（1/m）；I——声强（W/m2）；A——超声探头的振幅（m）。

1.3 边界条件及初始条件

马铃薯底面在解冻过程中一直与镂空金属支撑网接触，因其导热作用小，故忽略不计。马铃薯与解冻装置中的水发生对流换热，视为第三类边界条件[8]：

式中：h——对流传热系数，本模拟300 取W/（m2·K）[9]；T——马铃薯的表面温度；Twater——介质水的温度，17 ℃。

冻结的马铃薯内部各项均匀，温度一致，马铃薯的初始条件方程：

式中：T0——马铃薯的初始温度，试验测得为-18 ℃。

1.4 马铃薯的热物性参数

马铃薯在解冻过程中会发生相变，大部分冰晶融化成水，其热物性参数会发生较大变化。为准确模拟和预测食品在冻结过程中的温度场及解冻时间，采用分段多项式数学模型来反映马铃薯的热物性变化规律[10]。确定马铃薯的主要成分是合理计算马铃薯热物理参数的基础，其主要成分为：

X’wate=0.8；Xcarb=0.17；Xpro=0.02；Xash=0.01

其中，Xash：马铃薯中灰分的质量分数；Xcarb：马铃薯中碳水化合物的质量分数；X’wate：马铃薯中总的水分质量分数；Xpro：马铃薯中蛋白质的质量分数；Xw：马铃薯在解冻过程中已融化的水的分质量分数；Xice：马铃薯中在解冻过程中冰的质量分数。各组分的热物理性质参考食品工程原理[11]。

由Miles 等[12]可知，含冰率：

式中，Tf——初始冻结温度。根据冻结曲线可知，马铃薯的初始冻结温度为-2℃。

已解冻的水的质量分数：

在解冻过程中，冰融化为水，会释放大量的潜热。本验模拟采用固定网格法中的显热容法对潜热释放进行处理，即把相变潜热看作是在相变区域有一个很大的显热容，这样就能分区描述的控制方程及界面能量守恒条件转化成在整个区域适用的单一非线性导热方程。

马铃薯的比热容的计算式如下[12]：

式中：xi——各组分的质量分数；Cpi——各组分对应的比热容；L——水的潜热。

马铃薯的热导率的计算公式[13]：

式中：ki——各组分对应的热导率。

马铃薯的密度计算公式：

式中：ρi——各组分对应的密度。

马铃薯的比热容、热导率以及密度随温度变化的趋势如图2所示。

图2 马铃薯随温度变化的热物性参数Fig.2 Thermal properties of potato as functions of temperature

2 试验验证

2.1 材料与设备

马铃薯，杭州沃尔玛超市；1 000 mL 玻璃瓶；纯净水，哇哈哈集团有限公司；Data Trace MPⅢ无线温度传感器及工作站，Mesa Laboratories 公司；超声波细胞粉碎机，宁波新芝生物科技股份有限公司；计算机，因特尔处理器，2.93 GHz，24 564 M RAM。

2.2 试验方法

选用新鲜马铃薯，将其切割成半径为17 mm，高为25 mm 的圆柱体。在土豆中心位置插入温度探头，记录土豆在解冻过程中中心温度的变化。在-18 ℃下冻藏24 h，取出，置于超声解冻装置中超声解冻。对照组是在相同条件下仅去掉超声作用的水解冻。

2.3 试验验证

为了验证有限元模拟模型的精度，采用均方根误差（RMSE）对模拟及实测结果间的差异进行统计检验。

式中：Tm——模拟温度（℃）；T——马铃薯的实测中心温度（℃）；i——数据节点编号；n——数据点总数。

2.4 模型预测

验证后的模型用于预测马铃薯的半径（12.5，15.0，17.5，20.0，22.5，25.0 mm）和高度（17.5，20.0，22.5，25.0，27.5，30.0，32.5，35.0 mm），以及在声场中的位置（上表面距离探头45，50，55，60，65，70，75，80 mm）对马铃薯解冻效果的影响，为进一步优化解冻工艺提供理论依据。

3 结果与分析

3.1 马铃薯在解冻过程中的温度分布

马铃薯解冻所需热量以超声衰减产热和与介质的热交换为主，取30，80，270，422s 时马铃薯中心切面的温度分布来观察。

超声解冻和对照组解冻时马铃薯内部的温度变化如图3所示。由模拟结果可知，解冻30 s 时马铃薯与介质水发生剧烈的热交换，其边缘部分迅速升温，而内部温度依旧较低；超声解冻的马铃薯内部整体温度高于对照组，且低温区的面积小于对照组的低温面积。解冻80 s 时，马铃薯表面的温度依旧上升较快，超声解冻的马铃薯整体进入相变区，最低温度为-5 ℃，低温区逐渐减少并逐渐由几何中心向下移；对照组解冻的马铃薯边缘升温迅速，而内部的低温区面积仍然较大，且仍位于几何中心，热量由外向里逐渐传导。当解冻270 s 时，超声解冻的马铃薯上侧的升温速率明显高于下侧，低温区显著下移，马铃薯由外向内逐渐升温，仅一小片区域还处于相变区；对照组解冻的马铃薯解冻区域从边缘到中心逐渐增加，以较对称的形式解冻，中心的一大片区域仍处于相变阶段。解冻422 s 时，超声解冻的马铃薯最低温度达-1 ℃，解冻完成；对照组解冻的马铃薯中心部分仍处于相变期，根据计算可知，当解冻时间为822 s 时，对照组完成解冻。由此可知，超声解冻过程中，马铃薯内部组织吸收超声衰减产生的能量，外部和介质进行热交换，因而解冻速率显著大于对照组。

图3 不同解冻时间的中心切面温度场Fig.3 Center section temperature fields of different thawing time

3.2 试验验证

为了判断COMSOL 模拟的准确性，本文设计试验进行验证。如图4所示，超声辅助解冻时，COMSOL 计算得到的中心温度到达-1 ℃的时间为350 s，实测时间为326 s，较为接近。在整个解冻阶段，模拟值与实测值拟合较好，RMSE 为0.393 ℃。图5是对照组中心温度实测值和模拟值的对比，RMSE 为1.162 ℃，表明COMSOL 模拟能较好地反映超声解冻过程中温度分布和变化趋势，模型具有一定的精确性。

图4 超声辅助解冻时马铃薯中心温度模拟值和实测值对比Fig.4 Comparison of center temperatures of potato between measurement and simulation for ultrasound-assisted thawing

图5 对照组解冻时马铃薯中心温度模拟值和实测值对比Fig.5 Comparison of center temperatures of potato between measurement and simulation for control

3.3 设计参数对马铃薯解冻效果的影响

在成功模拟马铃薯解冻的基础上，改变该装置的设计参数，如马铃薯的尺寸、放置位置等，探究其对其解冻效果的影响。

3.3.1 马铃薯外形尺寸对解冻过程的影响为探究马铃薯外形尺寸对解冻过程的影响，分别改变半径以及高度来模拟计算。

如图6所示，完全解冻的时间随着马铃薯半径的变化而变化，半径越大，解冻时间越长。这是因为马铃薯解冻所需热量一部分来源于外界介质环境，通过对流换热的形式将能量从介质传递到马铃薯上，再以热传导的形式从外向里传递热量。当半径越大时，导热热阻越大，解冻时间越长[15-16]；然而当半径较小时，虽然会显著加快解冻速率，但会造成马铃薯对超声能量的不均匀吸收，导致马铃薯内部温度分布不均，温差较大。最大温差是指马铃薯完成解冻时，最高温度和最低温度（-1℃）的差值。解冻完成时的最大温差与马铃薯的半径没有线性关系，这可能是因为不同半径的马铃薯内部分布的声场不同，不同的声压场造成马铃薯内部对超声衰减能量不同程度的吸收。由图6可知，当半径为22.5 mm 时，马铃薯内部的温差最大，分布不均的声场导致能量不均匀地吸收，进而导致局部过热现象；当半径为15mm 时，解冻时间较短且内部温差较小。通过优化出待解冻食品的最佳半径，达到较好的解冻效果。

马铃薯的高度也对解冻时间有着显著的影响。如图7所示，随着马铃薯高度的增加，解冻时间延长；最大温差随着高度的增加依次降低。这是因为马铃薯高度增加，其热阻随之增大，因而解冻时间长；同时高度的增加也改变了内部声压场的分布，对超声能量的吸收程度也发生变化。由图7可知，随着高度的增加，声场在马铃薯内部的分布更加均匀，超声衰减的能量也更均匀，有利于减缓局部过热现象的发生。

3.3.2 马铃薯摆放位置对解冻效果的影响为了探究马铃薯在超声解冻装置中摆放位置对解冻效果的影响，将马铃薯分别放在其上表面距离探头45，50，55，60，65，70，75，80 mm 的地方，用COMSOL 来模拟计算。

如图8所示，解冻时间随着马铃薯放置位置的变化而变化，无线性关系。当距离为65 mm 时，解冻时间最少。这与装置中的声场分布有关，马铃薯的放置位置不仅影响装置中整体的声场分布，也影响其内部声场的分布。当所放位置的声压越大时，马铃薯吸收的能量也越多，因而解冻越迅速。当距离越近时，马铃薯内部温差越大，局部过热的现场也越明显。超声解冻物体时，需要根据物料的性质，优化出一个最佳的放置位置，在较快解冻的同时，避免局部过热导致的品质劣变。

图6 超声辅助解冻时不同半径的马铃薯所需解冻时间和内部最大温差Fig.6 Ultrasound-assisted thawing time and the biggest temperature difference at different radius of potato

图7 不同高度的马铃薯所需解冻时间和内部最大温差Fig.7 Ultrasound-assisted thawing time and the biggest temperature difference at different height of potato（final center temperature at -1 ℃）