食品速凍保鮮過程冰晶生長的相場模擬

楊 燕 袁訓鋒 喬希民 范 娜

(1. 商洛學院數(shù)學與計算機應用學院，陜西 商洛 726000；2. 商洛市分布式新能源應用技術研究中心，陜西 商洛 726000；3. 商洛學院健康管理學院，陜西 商洛 726000；4. 西北大學生命科學學院，陜西 西安 710029)

速凍保鮮方法能最大限度保持食品的原有風味、固有品質和營養(yǎng)價值，是目前食品保鮮領域最安全、經(jīng)濟、實用的手段[1-4]。在速凍保鮮過程中，熱量傳遞過程比水分滲透過程快，食品細胞內外的水被過冷凝固成冰晶分布其中，大冰晶易破壞細胞結構，降低食品的品質。因此，深入研究食品速凍保鮮過程中的冰晶生長行為，有助于改進速凍保鮮工藝、提升速凍食品品質。

冰晶形成過程中，無序狀態(tài)的液態(tài)水分子按特定對稱性向有序排列轉變，形成不同類型的冰晶形態(tài)，嚴重影響冰晶生長形態(tài)的這種特定對稱性稱為界面能各向異性。Liu等[5]曾采用試驗方法揭示冷表面冰晶生長模式，認為水滴遇冷凝結成對稱的雪花狀六邊形冰晶，界面能各向異性具有六重對稱性；當冰沿著表面生長時形成不同的形狀，界面能各向異性具有不同類型對稱性。

準確描述晶體生長必須考慮固液界面處的界面能各向異性。在弱界面能各向異性條件下，晶體形貌光滑，所有界面方向都可以生長；隨著界面能各向異性增加, 晶體形貌出現(xiàn)晶向缺失。近年來，針對凝固過程的強界面能各向異性相場模型得到深入發(fā)展。Eggleston等[6]建立強界面能各向異性相場模型，對小平面枝晶生長行為進行研究。隨后，肖榮振等[7-9]研究了強制流動條件下面心立方(FCC)結構Ni-Cu合金小平面枝晶生長行為。陳志等[10-12]深入研究了各向異性強度、流動速度對純物質和合金小平面枝晶生長的影響。袁訓鋒等[13-14]采用強界面能各向異性相場模型模擬密排六方(HCP)結構材料的枝晶生長行為，在界面能各向異性值通過臨界值(1/35)時主枝尖端穩(wěn)態(tài)生長速率降低，通過設置初始橢圓晶核長短半軸比值<1能夠消除主枝尖端分裂現(xiàn)象。針對冰晶生長的研究主要集中在四重對稱和六重對稱弱界面能各向異性方面。李方方等[15]嘗試引入相場變量描述系統(tǒng)各點是處于固態(tài)、液態(tài)或固液界面，采用各向同性相場方程再現(xiàn)冰晶的二維生長形貌。陳梅英等[16-19]建立冷凍濃縮過程冰晶生長相場模型，確定液態(tài)食品冷凍濃縮工藝中各向異性系數(shù)取值范圍為0.010～0.025，同時適當控制結晶時間能夠降低溶質夾帶損失。袁訓鋒等[20]以Wheeler模型為基礎，建立速凍保鮮過程冰晶生長相場模型，模擬結果顯示在高界面能各向異性強度下，易形成大冰晶破壞果蔬的細胞結構。這些研究表明，界面能各向異性嚴重影響冰晶生長形貌及穩(wěn)定性。王明華[21]指出英德國際研究小組發(fā)現(xiàn)冰在納米尺度上的平面結構為五邊形，這項發(fā)現(xiàn)將啟發(fā)人們不再局限六面體冰晶的認知，為了解冰晶的形成過程提供新的思路。然而，至今未見五重對稱強界面能各向異性下的冰晶生長行為的相關研究報道。

本研究擬在構建五重對稱強界面能各向異性相場模型的基礎上，將食品細胞內外溶液視為水和溶質二元系統(tǒng)，模擬再現(xiàn)冰晶生長演化過程，探討冷凍時間對冰晶生長形貌、溶質濃度、無量綱溫度的影響。

1 強界面能各向異性相場模型

在冰晶生長中，采用能量梯度系數(shù)反映界面能各向異性的影響，其表達形式為[13]：

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

(1)

式中：

θ——X軸與冰晶固液界面法向的夾角，rad；

w0——常數(shù)；

k——各向異性模數(shù)，對于五重對稱結構值取5；

εk——界面能各向異性強度；

f(θ)——能量梯度系數(shù)。

1.1 相場控制方程

將界面能各向異性函數(shù)f(θ)代入平衡晶體形態(tài)的參數(shù)方程[22]，整理可得：

(2)

式中：

X、Y——無量綱的距離，且X=X′/ω，Y=Y′/ω；

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

ω——參考長度，cm；

A——常數(shù)。

界面能各向異性強度取不同值時的平衡晶體形貌，如圖1 所示。從圖1中可以看出，當界面能各向異性強度ε5≤1/24(0.04)時，平衡晶體形態(tài)處處光滑連續(xù)，整體形貌為類五邊形。當界面能各向異性強度ε5>1/24時，五邊形頂點附近晶向方向在平衡晶體形貌中消失，優(yōu)先生長方向的界面不連續(xù)，出現(xiàn)失真的“耳子”。經(jīng)過修正后，這些“耳子”退化為一個點。

在晶向角周期|θ|≤π/5進行修正，當θm≤|θ|≤π/5時，晶體正常生長；當|θ|<θm時，晶向發(fā)生缺失。為了獲得準確的平衡晶體形貌，根據(jù)Eggleston等[6]提出修正界面能各向異性函數(shù)的方法，則：

(3)

出現(xiàn)晶向缺失現(xiàn)象的臨界晶向角θ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之間連續(xù)變化；

u——無量綱溫度，u=(T-Tm)/(Tm-T0)；

t——無量綱時間，t=t′/(ω2/κ)；

Ω——無量綱過冷度，Ω=cp[ΔT+mL(x-x0)]/L；

m——界面動力學系數(shù)，m=μσTm/(κL)；

T——溫度，K；

Tm——熔點，K；

T0——系統(tǒng)初始溫度，K；

ΔT——熱過冷度，K；

δ——界面層厚度，cm；

t′——時間，s；

κ——熱擴散率, cm2/s；

μ——界面遷移率；

mL——液相線斜率；

cp——定壓比熱容，J/(cm3·K)；

L——單位體積的結晶潛熱，J/cm3；

σ——界面能，J/cm2；

x0——過冷熔體初始濃度(摩爾分數(shù))；

x——熔體實際濃度(摩爾分數(shù))。

1.2 溫度場控制方程

溫度滿足的控制方程為：

(7)

式中：

p′(φ)——勢函數(shù)p(φ)=φ2(10-15φ+6φ2)對φ的導數(shù)。

1.3 溶質場控制方程

溶質濃度滿足的控制方程為：

(8)

式中：

Ds——固相溶質擴散系數(shù)，cm2/s；

Dl——液相溶質擴散系數(shù)，cm2/s；

k0——溶質平衡分配系數(shù)。

2 數(shù)值計算方法

2.1 初始條件和邊界條件

在二維空間選擇正方形模擬計算區(qū)域，網(wǎng)格節(jié)點數(shù)為1 200×1 200，設置初始晶核半徑為R0且位于正方形模擬區(qū)域中心。即：

(9)

式中：

X——[100]方向；

Y——[010]方向。

相場和溫度場變量均采用黎曼邊界條件。

2.2 模擬方法

(10)

式中：

Δt——時間步長；

ΔX、ΔY——空間步長；

m′——穩(wěn)定性系數(shù)，且m′=max(Dl，m)；

k′——修正系數(shù)，一般取1～2。

2.3 物性參數(shù)和模擬參數(shù)

3 結果與討論

3.1 強界面能各向異性下冰晶生長形貌

采用強界面能各向異性相場模型計算獲得冰晶充分生長的相場形貌和冰晶沿著固體表面的生長實際形貌，如圖2所示。從冰晶生長相場形貌可以看出，當晶粒生長限制在二維平面內，晶粒在5個方向上優(yōu)先生長形成冰晶主枝，各主枝之間呈72°夾角，冰晶形貌具有典型五重對稱性；在主枝中部具有發(fā)達的側枝，側枝與主枝呈72°夾角，側枝之間競爭生長，相鄰主枝的側枝甚至彼此接觸呈現(xiàn)融合趨勢；在主枝和側枝的尖端，由于存在晶向缺失其界面方向不連續(xù)，固液界面變得不穩(wěn)定且在相應位置出現(xiàn)棱角；主枝根部向內凹陷其生長緩慢形成“縮頸”，主枝和側枝的根部表面光滑包含全部的晶向方向。從冰晶實際形貌可以看出，主枝上具有發(fā)達的側向分支，主枝根部向內凹陷生長緩慢形成“縮頸”，主枝和側枝尖端輪廓不是拋物線且具有棱角，冰晶形貌不具備六重對稱。這些冰晶沿著固體表面生長的實際形態(tài)特征與模擬得到冰晶生長相場形貌一致。

圖2 冰晶生長的模擬結果與實際形貌對比Figure 2 Comparison between simulated result and actual growth morphology of ice crystal

橫坐標網(wǎng)格數(shù)分別為350，600，850的位置，繪制相場變量值隨縱坐標網(wǎng)格數(shù)的變化關系，如圖3所示。從圖3中可以看出，相場變量值曲線整體呈現(xiàn)左右對稱，表明冰晶生長呈現(xiàn)X軸對稱，X軸為五重對稱冰晶生長的對稱軸之一；橫坐標網(wǎng)格數(shù)350和850位置的相場變量值曲線差異較大，表明冰晶生長不具備Y軸對稱。相場變量值曲線波谷對應計算區(qū)域固相，曲線波峰對應冰晶生長的液相。當橫坐標網(wǎng)格數(shù)為350時，從一段平坦曲線液相區(qū)開始，穿過曲線波谷冰晶主枝，進入曲線波峰半封閉液相區(qū)，隨后再次經(jīng)過曲線波谷冰晶主枝，在平坦曲線液相區(qū)結束；當橫坐標網(wǎng)格數(shù)為600時，僅有1個固相區(qū)域的曲線波谷，在固液界面區(qū)域相場變量值存在微小波動；當橫坐標網(wǎng)格數(shù)為850時，從一段平坦曲線液相區(qū)開始，穿過5個曲線波谷冰晶側枝，進入曲線波峰半封閉液相區(qū)，隨后經(jīng)過曲線波谷冰晶主枝，再次穿過5個曲線波谷冰晶側枝，在平坦曲線液相區(qū)結束。

1. 橫坐標350 2. 橫坐標600 3. 橫坐標850

3.2 強界面能各向異性下冰晶生長溶質濃度和溫度分布

采用強界面能各向異性相場模型計算獲得冰晶充分生長的溶質濃度和無量綱溫度分布，如圖4所示。從冰晶溶質濃度分布可以看出，由于速凍保鮮過程降溫速率快，冰晶形成和生長速率快，冰晶凝固排擠的溶質在固液界面前沿擴散不充分，導致固液相溶質濃度相差不大，固相中心溶質濃度偏高；在主枝和側枝之間擴散通道不暢通，冰晶凝固排擠在固液界面前沿的溶質得不到及時遷移，導致形成高溶質濃度的液相區(qū)，冰晶生長受到抑制；在主側枝部分固液混合區(qū)，擴散通道暢通，冰晶凝固排擠的溶質擴散充分，形成低溶質濃度的固液混合區(qū)。從冰晶無量綱溫度分布可以看出，主枝之間的區(qū)域，由于凝固釋放的熱量得不到及時導走，使其溫度較高呈現(xiàn)深紅色，遠離固相的液相區(qū)域溫度較低為藍色；在固液界面前沿具有一定厚度的溫度擴散層，主枝尖端溫度擴散層薄其溫度梯度大，冰晶生長迅速，主枝之間的區(qū)域溫度擴散層厚其溫度梯度小，冰晶生長緩慢。

橫坐標網(wǎng)格數(shù)分別為350，600，850的位置，冰晶生長溶質濃度和溫度隨縱坐標網(wǎng)格數(shù)的變化關系，如圖5所示。在溶質濃度曲線中，主側枝之間半封閉液相區(qū)溶質濃度達到1.76% 為初始溶質濃度的3倍以上，主側枝的固液界面前沿存在溶質富集其溶質濃度較高，冰晶固相和液相溶質濃度在初始溶質濃度附近。在無量綱溫度曲線中，冰晶固相溫度低，冰晶主枝之間的熱量傳導不及時，區(qū)域溫度高。橫坐標網(wǎng)格數(shù)為850曲線靠近冰晶主枝根部的固液界面區(qū)域，曲線斜率小，溫度變化平緩；橫坐標網(wǎng)格數(shù)為350，600曲線靠近冰晶主枝尖端固液界面區(qū)域，曲線斜率大，溫度變化急劇。

3.3 冷凍時間對冰晶生長的影響

將速凍保鮮食品細胞內外溶液體系視為水和溶質二元系統(tǒng)，選用摩爾分數(shù)0.56%的糖水溶液作為研究對象，在過冷度T為0.9(71 K)時，取界面各向異性模數(shù)為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 不同橫坐標情況下溶質濃度和無量綱溫度隨縱坐標網(wǎng)格數(shù)的變化關系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可以看出，在冰晶生長初期，溶質和溫度擴散層緊緊包裹五邊形冰晶，在五邊形頂點處溶質和溫度擴散層薄，溶質梯度和溫度梯度大，冰晶生長迅速；在五邊形各邊中部溶質和溫度擴散層厚，溶質梯度和溫度梯度小，冰晶生長緩慢。隨著冷凍時間的推移，冰晶主枝生長排擠的溶質和熱量在主枝之間的根部區(qū)域固液界面前沿富集，此區(qū)域的冰晶生長受到抑制，冰晶縮頸現(xiàn)象加劇。繼續(xù)延長冷凍時間，主枝之間區(qū)域的溶質和熱量富集嚴重，形成半封閉的高溶質濃度和高溫區(qū)域；在主枝和側枝的部分固液混合區(qū)，擴散通道暢通，冰晶凝固排擠的溶質擴散充分，形成低溶質濃度固液混合區(qū)。

為了定量分析食品冷凍保鮮過程中冰晶生長尖端行為，繪制冰晶生長尖端無量綱溫度、溶質濃度、無量綱生長速率和曲率半徑隨冷凍時間的變化關系，如圖8所示。從圖8中可以看出，隨著冷凍時間的延長，無量綱溫度和溶質濃度先快速增加后逐漸趨于穩(wěn)定，無量綱溫度穩(wěn)態(tài)值為-0.12，溶質濃度穩(wěn)態(tài)值為0.63%；尖端生長速率和曲率半徑先急劇減小后逐漸趨于穩(wěn)定，無量綱生長速率穩(wěn)態(tài)值為19.45，曲率半徑穩(wěn)態(tài)值為0.015。

圖8 冰晶生長尖端行為隨冷凍時間的變化關系Figure 8 Tip behavior vs freezing time in ice crystal growth

4 結論

強界面能各向異性相場模型能夠成功再現(xiàn)速凍保鮮過程冰晶生長形貌，模擬結果與冰晶實際形貌對比發(fā)現(xiàn)，兩者獲得的冰晶主枝中部具有發(fā)達的側向分支，主枝根部向內凹陷生長緩慢形成“縮頸”，主枝和側枝的尖端具有棱角，整體形貌不具備六重對稱。在主枝和側枝之間擴散通道不暢通，冰晶凝固排擠在固液界面前沿的溶質和釋放的熱量得不到及時遷移，半封閉液相區(qū)溶質濃度達到1.76%，冰晶生長受到抑制。冰晶主枝根部的固液界面區(qū)域無量綱溫度變化平緩，冰晶主枝尖端固液界面區(qū)域無量綱溫度變化急劇。隨著冷凍時間延長，無量綱溫度和溶質濃度先快速增加后逐漸趨于穩(wěn)定，尖端生長速率和曲率半徑先急劇減小后逐漸趨于穩(wěn)定，冰晶從初期近似五邊形經(jīng)過5主枝結構向五重對稱多分支小平面結構轉變。速凍保鮮過程中食品細胞內外多分支小平面冰晶，容易破壞細胞結構，導致食品的品質下降。要避免細胞結構破壞，需把冷凍時間控制在合理的范圍。

在模擬速凍保鮮過程冰晶生長時，忽略了冰晶生長的其他模式，僅進行五重對稱界面能各向異性冰晶生長模擬，并且視食品細胞內外溶液體系為水和溶質二元系統(tǒng)。今后需要考慮細胞溶液的多元結構以及冰晶生長的其他模式，進行多元結構多模式冰晶生長模擬。