食品速冻保鲜过程冰晶生长的相场模拟

杨 燕 袁训锋 乔希民 范 娜

(1. 商洛学院数学与计算机应用学院，陕西 商洛 726000；2. 商洛市分布式新能源应用技术研究中心，陕西 商洛 726000；3. 商洛学院健康管理学院，陕西 商洛 726000；4. 西北大学生命科学学院，陕西 西安 710029)

速冻保鲜方法能最大限度保持食品的原有风味、固有品质和营养价值，是目前食品保鲜领域最安全、经济、实用的手段[1-4]。在速冻保鲜过程中，热量传递过程比水分渗透过程快，食品细胞内外的水被过冷凝固成冰晶分布其中，大冰晶易破坏细胞结构，降低食品的品质。因此，深入研究食品速冻保鲜过程中的冰晶生长行为，有助于改进速冻保鲜工艺、提升速冻食品品质。

冰晶形成过程中，无序状态的液态水分子按特定对称性向有序排列转变，形成不同类型的冰晶形态，严重影响冰晶生长形态的这种特定对称性称为界面能各向异性。Liu等[5]曾采用试验方法揭示冷表面冰晶生长模式，认为水滴遇冷凝结成对称的雪花状六边形冰晶，界面能各向异性具有六重对称性；当冰沿着表面生长时形成不同的形状，界面能各向异性具有不同类型对称性。

准确描述晶体生长必须考虑固液界面处的界面能各向异性。在弱界面能各向异性条件下，晶体形貌光滑，所有界面方向都可以生长；随着界面能各向异性增加, 晶体形貌出现晶向缺失。近年来，针对凝固过程的强界面能各向异性相场模型得到深入发展。Eggleston等[6]建立强界面能各向异性相场模型，对小平面枝晶生长行为进行研究。随后，肖荣振等[7-9]研究了强制流动条件下面心立方(FCC)结构Ni-Cu合金小平面枝晶生长行为。陈志等[10-12]深入研究了各向异性强度、流动速度对纯物质和合金小平面枝晶生长的影响。袁训锋等[13-14]采用强界面能各向异性相场模型模拟密排六方(HCP)结构材料的枝晶生长行为，在界面能各向异性值通过临界值(1/35)时主枝尖端稳态生长速率降低，通过设置初始椭圆晶核长短半轴比值<1能够消除主枝尖端分裂现象。针对冰晶生长的研究主要集中在四重对称和六重对称弱界面能各向异性方面。李方方等[15]尝试引入相场变量描述系统各点是处于固态、液态或固液界面，采用各向同性相场方程再现冰晶的二维生长形貌。陈梅英等[16-19]建立冷冻浓缩过程冰晶生长相场模型，确定液态食品冷冻浓缩工艺中各向异性系数取值范围为0.010～0.025，同时适当控制结晶时间能够降低溶质夹带损失。袁训锋等[20]以Wheeler模型为基础，建立速冻保鲜过程冰晶生长相场模型，模拟结果显示在高界面能各向异性强度下，易形成大冰晶破坏果蔬的细胞结构。这些研究表明，界面能各向异性严重影响冰晶生长形貌及稳定性。王明华[21]指出英德国际研究小组发现冰在纳米尺度上的平面结构为五边形，这项发现将启发人们不再局限六面体冰晶的认知，为了解冰晶的形成过程提供新的思路。然而，至今未见五重对称强界面能各向异性下的冰晶生长行为的相关研究报道。

本研究拟在构建五重对称强界面能各向异性相场模型的基础上，将食品细胞内外溶液视为水和溶质二元系统，模拟再现冰晶生长演化过程，探讨冷冻时间对冰晶生长形貌、溶质浓度、无量纲温度的影响。

1 强界面能各向异性相场模型

在冰晶生长中，采用能量梯度系数反映界面能各向异性的影响，其表达形式为[13]：

w(θ)=w0f(θ)=w0[1+εk(kθ)]，

(1)

式中：

θ——X轴与冰晶固液界面法向的夹角，rad；

w0——常数；

k——各向异性模数，对于五重对称结构值取5；

εk——界面能各向异性强度；

f(θ)——能量梯度系数。

1.1 相场控制方程

将界面能各向异性函数f(θ)代入平衡晶体形态的参数方程[22]，整理可得：

(2)

式中：

X、Y——无量纲的距离，且X=X′/ω，Y=Y′/ω；

X′、Y′——距离，cm；

ω——参考长度，cm；

A——常数。

界面能各向异性强度取不同值时的平衡晶体形貌，如图1 所示。从图1中可以看出，当界面能各向异性强度ε5≤1/24(0.04)时，平衡晶体形态处处光滑连续，整体形貌为类五边形。当界面能各向异性强度ε5>1/24时，五边形顶点附近晶向方向在平衡晶体形貌中消失，优先生长方向的界面不连续，出现失真的“耳子”。经过修正后，这些“耳子”退化为一个点。

在晶向角周期|θ|≤π/5进行修正，当θm≤|θ|≤π/5时，晶体正常生长；当|θ|<θm时，晶向发生缺失。为了获得准确的平衡晶体形貌，根据Eggleston等[6]提出修正界面能各向异性函数的方法，则：

(3)

出现晶向缺失现象的临界晶向角θm满足：

f(θm)sin(θm)+fθ(θm)cos(θm)=0。

(4)

通过求解方程可以得到θm。

图1 界面能各向异性强度取不同值时的平衡晶体形貌Figure 1 Equilibrium crystal shape at various interfacial energy anisotropies

基于Wheeler模型建立五重对称结构的强界面能各向异性相场方程，表达式为：

当θm≤|θ|≤π/5时

(5)

当|θ|<θm时

(6)

式中：

φ——相场变量，φ=0代表固相，φ=1代表液相，固/液界面上φ在0→1之间连续变化；

u——无量纲温度，u=(T-Tm)/(Tm-T0)；

t——无量纲时间，t=t′/(ω2/κ)；

Ω——无量纲过冷度，Ω=cp[ΔT+mL(x-x0)]/L；

m——界面动力学系数，m=μσTm/(κL)；

T——温度，K；

Tm——熔点，K；

T0——系统初始温度，K；

ΔT——热过冷度，K；

δ——界面层厚度，cm；

t′——时间，s；

κ——热扩散率, cm2/s；

μ——界面迁移率；

mL——液相线斜率；

cp——定压比热容，J/(cm3·K)；

L——单位体积的结晶潜热，J/cm3；

σ——界面能，J/cm2；

x0——过冷熔体初始浓度(摩尔分数)；

x——熔体实际浓度(摩尔分数)。

1.2 温度场控制方程

温度满足的控制方程为：

(7)

式中：

p′(φ)——势函数p(φ)=φ2(10-15φ+6φ2)对φ的导数。

1.3 溶质场控制方程

溶质浓度满足的控制方程为：

(8)

式中：

Ds——固相溶质扩散系数，cm2/s；

Dl——液相溶质扩散系数，cm2/s；

k0——溶质平衡分配系数。

2 数值计算方法

2.1 初始条件和边界条件

在二维空间选择正方形模拟计算区域，网格节点数为1 200×1 200，设置初始晶核半径为R0且位于正方形模拟区域中心。即：

(9)

式中：

X——[100]方向；

Y——[010]方向。

相场和温度场变量均采用黎曼边界条件。

2.2 模拟方法

(10)

式中：

Δt——时间步长；

ΔX、ΔY——空间步长；

m′——稳定性系数，且m′=max(Dl，m)；

k′——修正系数，一般取1～2。

2.3 物性参数和模拟参数

3 结果与讨论

3.1 强界面能各向异性下冰晶生长形貌

采用强界面能各向异性相场模型计算获得冰晶充分生长的相场形貌和冰晶沿着固体表面的生长实际形貌，如图2所示。从冰晶生长相场形貌可以看出，当晶粒生长限制在二维平面内，晶粒在5个方向上优先生长形成冰晶主枝，各主枝之间呈72°夹角，冰晶形貌具有典型五重对称性；在主枝中部具有发达的侧枝，侧枝与主枝呈72°夹角，侧枝之间竞争生长，相邻主枝的侧枝甚至彼此接触呈现融合趋势；在主枝和侧枝的尖端，由于存在晶向缺失其界面方向不连续，固液界面变得不稳定且在相应位置出现棱角；主枝根部向内凹陷其生长缓慢形成“缩颈”，主枝和侧枝的根部表面光滑包含全部的晶向方向。从冰晶实际形貌可以看出，主枝上具有发达的侧向分支，主枝根部向内凹陷生长缓慢形成“缩颈”，主枝和侧枝尖端轮廓不是抛物线且具有棱角，冰晶形貌不具备六重对称。这些冰晶沿着固体表面生长的实际形态特征与模拟得到冰晶生长相场形貌一致。

图2 冰晶生长的模拟结果与实际形貌对比Figure 2 Comparison between simulated result and actual growth morphology of ice crystal

横坐标网格数分别为350，600，850的位置，绘制相场变量值随纵坐标网格数的变化关系，如图3所示。从图3中可以看出，相场变量值曲线整体呈现左右对称，表明冰晶生长呈现X轴对称，X轴为五重对称冰晶生长的对称轴之一；横坐标网格数350和850位置的相场变量值曲线差异较大，表明冰晶生长不具备Y轴对称。相场变量值曲线波谷对应计算区域固相，曲线波峰对应冰晶生长的液相。当横坐标网格数为350时，从一段平坦曲线液相区开始，穿过曲线波谷冰晶主枝，进入曲线波峰半封闭液相区，随后再次经过曲线波谷冰晶主枝，在平坦曲线液相区结束；当横坐标网格数为600时，仅有1个固相区域的曲线波谷，在固液界面区域相场变量值存在微小波动；当横坐标网格数为850时，从一段平坦曲线液相区开始，穿过5个曲线波谷冰晶侧枝，进入曲线波峰半封闭液相区，随后经过曲线波谷冰晶主枝，再次穿过5个曲线波谷冰晶侧枝，在平坦曲线液相区结束。

1. 横坐标350 2. 横坐标600 3. 横坐标850

3.2 强界面能各向异性下冰晶生长溶质浓度和温度分布

采用强界面能各向异性相场模型计算获得冰晶充分生长的溶质浓度和无量纲温度分布，如图4所示。从冰晶溶质浓度分布可以看出，由于速冻保鲜过程降温速率快，冰晶形成和生长速率快，冰晶凝固排挤的溶质在固液界面前沿扩散不充分，导致固液相溶质浓度相差不大，固相中心溶质浓度偏高；在主枝和侧枝之间扩散通道不畅通，冰晶凝固排挤在固液界面前沿的溶质得不到及时迁移，导致形成高溶质浓度的液相区，冰晶生长受到抑制；在主侧枝部分固液混合区，扩散通道畅通，冰晶凝固排挤的溶质扩散充分，形成低溶质浓度的固液混合区。从冰晶无量纲温度分布可以看出，主枝之间的区域，由于凝固释放的热量得不到及时导走，使其温度较高呈现深红色，远离固相的液相区域温度较低为蓝色；在固液界面前沿具有一定厚度的温度扩散层，主枝尖端温度扩散层薄其温度梯度大，冰晶生长迅速，主枝之间的区域温度扩散层厚其温度梯度小，冰晶生长缓慢。

横坐标网格数分别为350，600，850的位置，冰晶生长溶质浓度和温度随纵坐标网格数的变化关系，如图5所示。在溶质浓度曲线中，主侧枝之间半封闭液相区溶质浓度达到1.76% 为初始溶质浓度的3倍以上，主侧枝的固液界面前沿存在溶质富集其溶质浓度较高，冰晶固相和液相溶质浓度在初始溶质浓度附近。在无量纲温度曲线中，冰晶固相温度低，冰晶主枝之间的热量传导不及时，区域温度高。横坐标网格数为850曲线靠近冰晶主枝根部的固液界面区域，曲线斜率小，温度变化平缓；横坐标网格数为350，600曲线靠近冰晶主枝尖端固液界面区域，曲线斜率大，温度变化急剧。

3.3 冷冻时间对冰晶生长的影响

将速冻保鲜食品细胞内外溶液体系视为水和溶质二元系统，选用摩尔分数0.56%的糖水溶液作为研究对象，在过冷度T为0.9(71 K)时，取界面各向异性模数为5，各向异性强度为0.1，模拟计算获得不同冷冻时间下相场形貌、溶质浓度和无量纲温度分布，如图6～7所示。

图4 强界面能各向异性冰晶生长溶质浓度和无量纲温度分布Figure 4 Solute concentration and dimensionless temperature distribution at strong interface energy anisotropy for ice crystal growth

1. 横坐标350 2. 横坐标600 3. 横坐标850图5 不同横坐标情况下溶质浓度和无量纲温度随纵坐标网格数的变化关系Figure 5 Value of solute concentration and dimensionless temperature vs. vertical coordinate grids number for ice crystal at various horizontal grids

图6 冰晶生长形貌随冷冻时间变化过程Figure 6 Morphology of ice crystal growth for various freezing times

图7 冰晶生长溶质浓度、温度分布随冷冻时间变化过程Figure 7 Solute concentration and dimensionless temperature of ice crystal growth for various freezing times

从图7可以看出，在冰晶生长初期，溶质和温度扩散层紧紧包裹五边形冰晶，在五边形顶点处溶质和温度扩散层薄，溶质梯度和温度梯度大，冰晶生长迅速；在五边形各边中部溶质和温度扩散层厚，溶质梯度和温度梯度小，冰晶生长缓慢。随着冷冻时间的推移，冰晶主枝生长排挤的溶质和热量在主枝之间的根部区域固液界面前沿富集，此区域的冰晶生长受到抑制，冰晶缩颈现象加剧。继续延长冷冻时间，主枝之间区域的溶质和热量富集严重，形成半封闭的高溶质浓度和高温区域；在主枝和侧枝的部分固液混合区，扩散通道畅通，冰晶凝固排挤的溶质扩散充分，形成低溶质浓度固液混合区。

为了定量分析食品冷冻保鲜过程中冰晶生长尖端行为，绘制冰晶生长尖端无量纲温度、溶质浓度、无量纲生长速率和曲率半径随冷冻时间的变化关系，如图8所示。从图8中可以看出，随着冷冻时间的延长，无量纲温度和溶质浓度先快速增加后逐渐趋于稳定，无量纲温度稳态值为-0.12，溶质浓度稳态值为0.63%；尖端生长速率和曲率半径先急剧减小后逐渐趋于稳定，无量纲生长速率稳态值为19.45，曲率半径稳态值为0.015。

图8 冰晶生长尖端行为随冷冻时间的变化关系Figure 8 Tip behavior vs freezing time in ice crystal growth

4 结论

强界面能各向异性相场模型能够成功再现速冻保鲜过程冰晶生长形貌，模拟结果与冰晶实际形貌对比发现，两者获得的冰晶主枝中部具有发达的侧向分支，主枝根部向内凹陷生长缓慢形成“缩颈”，主枝和侧枝的尖端具有棱角，整体形貌不具备六重对称。在主枝和侧枝之间扩散通道不畅通，冰晶凝固排挤在固液界面前沿的溶质和释放的热量得不到及时迁移，半封闭液相区溶质浓度达到1.76%，冰晶生长受到抑制。冰晶主枝根部的固液界面区域无量纲温度变化平缓，冰晶主枝尖端固液界面区域无量纲温度变化急剧。随着冷冻时间延长，无量纲温度和溶质浓度先快速增加后逐渐趋于稳定，尖端生长速率和曲率半径先急剧减小后逐渐趋于稳定，冰晶从初期近似五边形经过5主枝结构向五重对称多分支小平面结构转变。速冻保鲜过程中食品细胞内外多分支小平面冰晶，容易破坏细胞结构，导致食品的品质下降。要避免细胞结构破坏，需把冷冻时间控制在合理的范围。

在模拟速冻保鲜过程冰晶生长时，忽略了冰晶生长的其他模式，仅进行五重对称界面能各向异性冰晶生长模拟，并且视食品细胞内外溶液体系为水和溶质二元系统。今后需要考虑细胞溶液的多元结构以及冰晶生长的其他模式，进行多元结构多模式冰晶生长模拟。