虾仁热物性的计算及冻结时间的数值模拟

唐 婉

谢 晶1,2,3,4

王金锋1,2,3,4

周洪剑5

(1.上海海洋大学，上海 201306；2.上海水产品加工与贮藏工程技术研究中心，上海 201306；3.上海冷链装备性能与节能评价专业技术服务平台，上海 201306；4.食品科学与工程国家级实验教学示范中心﹝上海海洋大学﹞，上海 201306；5.上海宝丰机械制造有限公司，上海 201306)

速冻食品因其具有降低酶的活性、抑制微生物的生长以及延长食品的货架期等优点[1]，产量逐年递增，而食品的冻结时间以及冻结过程中食品内部的温度分布是直接影响食品品质以及冻结设备能耗的重要因素之一[2-3]。因此准确预测食品冻结时间对冻结设备的优化设计以及食品品质的提升均具有重要的意义。为获取食品冻结时间的信息，传统的方法是理论计算法、经验公式法，但只适合于规则的、简单的几何模型，局限性大并且获得信息量少。而数值模拟法适用范围广，对于一些不规则、复杂的冷冻过程，数值模拟法是首要选择，并且可以实时监测冻品任何位置的温度，获取大量的信息[4-6]。刘庆等[7]对鲢鱼肉冷冻过程进行模拟计算，预测冻结时间与实验冻结时间的相对误差为-1.64%。Cleland D J 等[8]研究表明在进行食品冻结时间或解冻时间的模拟计算时，有限元法准确性较好。李晓宇等[9]利用COMSOL软件模拟不同风速条件下，马铃薯冻结过程中相变界面的移动规律，模拟结果与试验值吻合较好。由于虾仁的形状复杂，模拟时网格划分较难等，前人关于虾仁冻结的模拟研究较少，且都对模型大大简化。杨小梅等[10]研究了南美白对虾在冷冻柜中冻结过程的模拟，选取对虾的横截面为研究对象，建立了虾仁的二维模型，进行模拟验证。李杰等[11]模拟了虾仁在鼓风冻结设备中冻结的过程，简化了虾仁的模型，以长方形代替虾仁的纵截面，并建立了二维模型，但模拟结果仅局限于某个面，无法直观获取虾仁整体、任意截面的温度云图，以及实时获取虾仁内部各个位置的温度，存在模型的局限性。

本试验针对冰箱冻结虾仁的冻结过程进行了数值模拟和实验验证研究。利用CFD数值模拟技术，建立虾仁的三维有限元非稳态模型，测得虾仁的含水率，按经验公式计算得出虾仁的导热系数、比热容等参数，调入到模拟环境，模拟计算虾仁在冰箱不同位置处的冻结时间，并通过实验验证模型及参数设置的准确性，以期为冰箱冻结虾仁以及提高冻结食品的品质相关模拟研究提供参考。

1 虾仁各参数的测定与计算

1.1 虾仁含水率的测量

采购新鲜的活虾10只，擦干表面水分，放置于干净无水塑料盒子中，去掉虾头、虾尾、虾线、虾皮等，迅速切碎制成虾糜备用，按照GB 5009.3—2011直接干燥法测量各个虾仁的含水率，取平均值，计算得出虾仁含水率为78.4%。

1.2 虾仁比热容的计算

食品冻结过程中比热容实时变化，一般很难准确测量，本试验参照文献[12]39-42的计算公式进行计算，具体如下：

冻结前：C1=0.837+3.34ω，

(1)

冻结后：C2=0.837+1.256ω，

(2)

式中：

ω——含水量，%；

C1——冻结前的比热容，kJ/(kg·K)；

C2——冻结后的比热容，kJ/(kg·K)。

经计算：冻结前虾仁的导热系数为3.45 kJ/(kg·K)，冻结后虾仁的导热系数为1.82 kJ/(kg·K)。

食品的相变过程是在一个温度范围内发生，而非某一个温度值。从食品刚到达初始冻结温度时，一部分自由水冻结成冰，食品内部的浓度增加，冻结点温度下降，直至食品完全冻结。食品冻结过程中水变成冰相变逐步进行，而相变过程中产生大量的潜热。

潜热值与含水率有关，潜热量Q=ωγ=0.784×334.9=262.6 kJ/kg，本试验采用等效比热容法，将潜热换算到比热容中，再带入模型中进行换算。

1.3 虾仁热导率的计算

食品的导热系数是一个复杂的参数，不仅与食品的组成成分有关，还与食品的结构特点(如是否多孔介质、是否属于均相物质等)有关。虾肉可视为均相物质，含水率高，而水的导热系数远远大于其内部其他组成成分的导热系数，因此其内部的导热系数主要与含水率有关，可用如下公式计算[12]49-50：

冻结前：λ1=0.26+0.34ω，

(3)

冻结后：λ2=2ω+0.22(1-ω)，

(4)

式中：

ω——含水量，%；

λ1——冻结前的导热系数，W/(m·K)；

λ2——冻结后的导热系数，W/(m·K)。

计算得：冻结前虾仁的导热系数为0.53 W/(m·K)，冻结后虾仁的导热系数为1.62 W/(m·K)。

2 数值模拟

2.1 物理模型

本研究采用西门子KK25F55TI型冰箱，虾仁置于冷冻室的底层(尺寸为400 mm×200 mm×260 mm)，该冷冻室有上下两层排管式蒸发器进行降温，虾仁的质量为2.542 5 g，不考虑冷冻室内的管道设计对流场的影响。具体模型图与食品摆放位置见图1、2。

图1 冰箱冷冻室以及食品摆放示意图Figure 1 The freezer and food

图2 食品摆放位置的主视图Figure 2 Food location in main view

2.2 控制方程及离散

本试验是研究虾仁悬浮于冰箱冷冻室内的冻结过程，虾仁四周与冷空气直接接触，虾仁表面通过与空气自然对流传热，虾仁内部通过导热进行传热。由于虾仁是一种多组分物质，其内部材料各向异性，且在冻结发生相变时，液相和固相成分会不断变化，使计算比较困难，为简化模型作出以下假设：① 虾仁的初始温度均匀一致，且冻结环境温度一定；② 虾仁内部各向同性，相变温度一定；③ 虾仁的内部只靠导热进行热传递，且虾仁表面对流换热系数保持不变；④ 整个冻结过程，虾仁的密度保持不变，为1 065 kg/m3[9]；⑤ 整个冻结过程中，边界条件保持不变；⑥ 虾仁内部非凝固区不考虑传质与流动。

基于上述假设，食品冷冻过程的热传导方程为：

(5)

式中：

T——虾仁各几何点的温度，K；

t——时间，s；

ρ——虾仁的密度，kg/m3；

C——虾仁的比热容，kJ/(kg·K)；

λ——虾仁的导热系数，W/(m·K)；

q——内热源，W/m3。

边界条件：

(1) 食品表面第三类边界条件：

q=h(Tw-Tf)，

(6)

式中：

Tw——虾仁表面的温度， K；

Tf——冰箱冷冻室内空气的温度(由试验测量获得)，K。

表面对流换热通过系统耦合处理[13-14]。

(2) 在冷冻室上下蒸发盘管板：Tban=244 K，按试验测量赋值。

(3) 在冷冻室侧壁绝热处理：q=0。

(4) 在移动相界面上满足质量守恒和能量守恒：

TS[S(t),t]=Tl[S(t),t]=Tp，

(7)

(8)

式中：

TS——相界面上固相的温度，K；

Tl——相界面上液相的温度，K；

Tp——相界面上相变的温度，K；

λS——相界面上固相的导热率，W/(m·K)；

λl——相界面上液相的导热率，W/(m·K)；

h——比焓，kJ/kg；

ρ——密度，kg/m3。

初始条件：

Tfood(x,y,z)=T0=290 K，食品温度按试验测量赋值。

Tf(x,y,z)=246 K，冷冻室气流温度按试验测量赋值。

2.3 虾仁冻结过程的数值模拟

基于以上研究，本次模拟冷冻室内不同高度虾仁的冻结过程，3只虾分别距冷冻室底部15，50，100 mm进行冷冻，为研究虾仁在冰箱冷冻室内冻结过程温度分布情况，首先用CAD建立了冰箱冷冻室以及虾仁的几何模型，并在ANSYS 15.0中自带的ICEM进行网格划分并细化。划分好网格后，调入到模拟环境，选择标准的k-ε计算模型，将计算虾仁的物性参数添加到材料，激活能量项，能量方程的收敛精度取10-6。其次建立虾仁在冰箱冻结的三维隐式非稳态模型，并采用有限元法离散控制方程。最后设置初始条件，分别给予虾仁的初始温度以及冰箱冷冻室内环境温度赋值，初始化流场后进行计算，时间步长取0.5 s，为了获取虾仁的冻结时间，分别对3组虾仁的最高温度进行监测，实时导出数据，直至最高温度低于255.15 K，模拟结束。为清晰观察虾仁内部的温度分布情况，截取不同高度虾仁的3个对称横截面，分别为15，50，100 mm。

由图3～5可知，纵向比较同一位置虾仁冻结过程中，虾仁尾部温度下降最快，虾仁顶部和侧面次之。由图5可知，在虾仁冻结330 s时，横向比较不同位置虾仁的冻结过程中，虾仁热中心区域位于虾仁上半部分。这是因为虾尾和虾仁侧表面与冷冻室内冷空气直接接触进行流固耦合，传热较快，而虾尾的截面尺寸小，降温速度较快。虾仁内部温度变化是通过导热作用从虾仁表面向内部传热，因而降温相对较慢。横向比较不同高度虾仁的冻结过程，可以看出距冷冻室底部15 mm处的虾仁降温最快，50 mm次之，100 mm降温最慢。这是因为15 mm处的虾仁距离冷源蒸发盘管距离近，蒸发盘管温度低于冷冻室环境温度，位于下部的虾仁与环境温度发生换热，高温气体会往上流动，致使虾仁降温速度随着距离蒸发盘管距离的增加而逐渐下降。

图3 冷冻时间为60 s时虾仁内部的温度分布云图Figure 3 Temperature distribution in shrimp with freezing time of 60 s

图4 冷冻时间为330 s时虾仁内部的温度分布云图Figure 4 Temperature distribution in shrimp with freezing time of 330 s

图5 冷冻时间为1 358 s时虾仁内部的温度分布云图Figure 5 Temperature distribution in shrimp with freezing time of 1 358 s

图6 不同时刻虾仁剖面温度分布云图Figure 6 Section temperature distribution in shrimp with different freezing time

3 实验验证

取距冰箱冷冻层底部15 mm的虾仁为对象，实验前冰箱空载运行1 d，保证冰箱冻结室内温度恒定(分别用标定好的热电偶，测量冰箱冷冻室的温度，达到稳定状态后取平均值)。基于以上内容先将Fluke 2640A网络型多点温度采集仪与计算机连接好，设定温度采集时间间隔为10 s，其次将标定好的美国OMEGA四氟测温线TT-T-36型热电偶温度探头固定距虾仁顶部10 mm处(经预试验后确定)，经测量虾仁的最初温度290.65 K。然后将固定好热电偶的虾仁放到冰箱指定位置，迅速关好冰箱门开始冻结。由于虾仁形状较小，实验测点不易固定，因此重复上述实验，分别测出不同组虾仁在该点的温度变化曲线，结果表明虾仁在该点的温度曲线基本一致，最大误差在30 s之内。

从虾仁的冻结曲线可以看出。虾仁冻结过程分为3个阶段，初始阶段，虾仁温度降速较快直至冻结点附近，测点冻结速率斜率有所下降，可能是虾仁表面与环境发生热交换，使周围环境温度略有上升造成的。当冻结过程到达相变温度时，相变过程释放大量的潜热，从而大大降低了冻结速率，温度趋于平缓。最后随着冻结过程进一步深入，虾仁内部大量的水结成冰，而冰的比热低于水，但是由于虾仁内部还残留少部分水分结冰释放潜热，使得降温速度低于初始阶段[15-16]。从图7模拟值与实验值的对比可以看出，冻结初期模拟值与实验值吻和度较高，随着冻结进入相变期，模拟值冻结温度略高于实验值，可能是模拟时设置的相变温度恒定，而实验过程中随着相变过程的发生，虾仁溶液浓度增加，相变温度逐渐下降产生的。冻结末期模拟值与实验值差距相对较大，可能是冻结后虾仁经计算得出的比热容比实际虾仁的比热容高造成的。经计算模拟冻结时间与实验冻结时间曲线相关系数为0.996，虾仁切面中心测点的最大误差为1.85 K，误差平均值为1.105 K，由此可知模拟技术能较为准确地预测食品的冻结时间。

图7 模拟值与实验值的对比Figure 7 Comparison of simulation value and experiment value

4 结论

本试验对冰箱冻结虾仁进行数值模拟，采用等效热容法处理相变潜热，建立描述三维食品冻结过程的微分方程，预测了虾仁冻结过程中温度分布以及冻结时间，并考虑了虾仁在冷冻室的位置对冷冻过程的影响，所得的温度云图能直观地观察到随着距冰箱冷冻室底部蒸发盘管距离的增加，冻结时间逐渐增加；对模拟结果进行实验验证表明数值模拟能很好地反映虾仁内部温度场的变化，模拟冻结时间与实验冻结时间曲线相关系数为0.996，虾仁切面中心测点的最大误差值为1.85 K，由此可知数值模拟可以有效预测食品冻结过程中温度的分布，对虾仁的冷冻加工、品质控制以及速冻设备的设计优化具有重要意义。

[1] AWAD T S, MOHARRAM H A, SHALTOUT O E, et al.Applications of ultrasound in analysis, processing and quality control of food: A review[J].Food Research International, 2012, 48(2): 410-427.

[2] JI Min, NIN Xin-chao, TANG Bo-ning, et al.Numerical simulation and experimental study on thawing time of cylindrical frozen food[J].Advanced Materials Research, 2014, 989-994: 3 513-3 517.

[3] 王金锋, 李文俊, 谢晶.数值模拟在食品冻结过程中的应用[J].食品与机械, 2016, 32(10): 200-204.

[4] NORTON T, SUN Da-wei.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.

[6] 唐婉, 王金锋, 李文俊, 等.牛肉冻结过程中模拟及热值传递理论分析[J].食品与机械, 2017, 33(7): 117-121.

[7] 刘庆, 熊善柏, 鲁长新, 等.鲢鱼肉冷冻过程中传热预测模型的建立与验证[J].现代食品科技, 2017(3): 110-115.

[8] CLELAND D J, CLELAND A C, EARLE R L, et al.Prediction of freezing and thawing times for multi-dimensional shapes by numerical methods[J].International Journal of Refrigeration, 1987, 10(1): 32-39.

[9] 李晓宇, 刘斌, 邸倩倩, 等.马铃薯冻结过程中不同风速对相变界面的影响[J].食品科技, 2015(4): 76-81.

[10] 杨小梅.对虾冷冻过程的温度分析及其数值模拟[D].杭州: 浙江大学, 2014: 50-61.

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

[12] 华泽钊, 李云飞, 刘宝林.食品冷冻冷藏原理与设备[M].北京: 机械工业出版社, 1999.

[13] 王贵强, 邹平华, 刘明生, 等.包装对食品冻结过程的影响研究[J].农业机械学报, 2014, 45(7): 171-176.

[14] MORGAN V T.The overall convective heat transfer from smooth circular cylinders[J].Advances in Heat Transfer, 1975, 11: 199-264.

[15] XANTHAKIS E, HAVET M, CHEVALLIER S, et al.Effect of static electric field on ice crystal size reduction during freezing of pork meat[J].Innovative Food Science & Emerging Technologies, 2013, 20(4): 115-120.

[16] PHAM Q T.Freezing time formulas for foods with low moisture content, low freezing point and for cryogenic freezing[J].Journal of Food Engineering, 2014, 127(4): 85-92.