楊 燕 袁訓(xùn)鋒 喬希民 范 娜

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

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

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

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

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

1 強(qiáng)界面能各向異性相場模型

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

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

(1)

式中：

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

w0——常數(shù)；

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

εk——界面能各向異性強(qiáng)度；

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

1.1 相場控制方程

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

(2)

式中：

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

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

ω——參考長度，cm；

A——常數(shù)。

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

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

(3)

出現(xiàn)晶向缺失現(xiàn)象的臨界晶向角θm滿足：

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

(4)

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

圖1 界面能各向異性強(qiáng)度取不同值時(shí)的平衡晶體形貌Figure 1 Equilibrium crystal shape at various interfacial energy anisotropies

基于Wheeler模型建立五重對稱結(jié)構(gòu)的強(qiáng)界面能各向異性相場方程，表達(dá)式為：

當(dāng)θm≤|θ|≤π/5時(shí)

(5)

當(dāng)|θ|<θm時(shí)

(6)

式中：

φ——相場變量，φ=0代表固相，φ=1代表液相，固/液界面上φ在0→1之間連續(xù)變化；

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

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

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

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

T——溫度，K；

Tm——熔點(diǎn)，K；

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

ΔT——熱過冷度，K；

δ——界面層厚度，cm；

t′——時(shí)間，s；

κ——熱擴(kuò)散率, cm2/s；

μ——界面遷移率；

mL——液相線斜率；

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

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

σ——界面能，J/cm2；

x0——過冷熔體初始濃度(摩爾分?jǐn)?shù))；

x——熔體實(shí)際濃度(摩爾分?jǐn)?shù))。

1.2 溫度場控制方程

溫度滿足的控制方程為：

(7)

式中：

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

1.3 溶質(zhì)場控制方程

溶質(zhì)濃度滿足的控制方程為：

(8)

式中：

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

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

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

2 數(shù)值計(jì)算方法

2.1 初始條件和邊界條件

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

(9)

式中：

X——[100]方向；

Y——[010]方向。

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

2.2 模擬方法

(10)

式中：

Δt——時(shí)間步長；

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

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

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

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

3 結(jié)果與討論

3.1 強(qiáng)界面能各向異性下冰晶生長形貌

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

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

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

1. 橫坐標(biāo)350 2. 橫坐標(biāo)600 3. 橫坐標(biāo)850

3.2 強(qiáng)界面能各向異性下冰晶生長溶質(zhì)濃度和溫度分布

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

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

3.3 冷凍時(shí)間對冰晶生長的影響

將速凍保鮮食品細(xì)胞內(nèi)外溶液體系視為水和溶質(zhì)二元系統(tǒng)，選用摩爾分?jǐn)?shù)0.56%的糖水溶液作為研究對象，在過冷度T為0.9(71 K)時(shí)，取界面各向異性模數(shù)為5，各向異性強(qiáng)度為0.1，模擬計(jì)算獲得不同冷凍時(shí)間下相場形貌、溶質(zhì)濃度和無量綱溫度分布，如圖6～7所示。

圖4 強(qiáng)界面能各向異性冰晶生長溶質(zhì)濃度和無量綱溫度分布Figure 4 Solute concentration and dimensionless temperature distribution at strong interface energy anisotropy for ice crystal growth

1. 橫坐標(biāo)350 2. 橫坐標(biāo)600 3. 橫坐標(biāo)850圖5 不同橫坐標(biāo)情況下溶質(zhì)濃度和無量綱溫度隨縱坐標(biāo)網(wǎng)格數(shù)的變化關(guān)系Figure 5 Value of solute concentration and dimensionless temperature vs. vertical coordinate grids number for ice crystal at various horizontal grids

圖6 冰晶生長形貌隨冷凍時(shí)間變化過程Figure 6 Morphology of ice crystal growth for various freezing times

圖7 冰晶生長溶質(zhì)濃度、溫度分布隨冷凍時(shí)間變化過程Figure 7 Solute concentration and dimensionless temperature of ice crystal growth for various freezing times

從圖7可以看出，在冰晶生長初期，溶質(zhì)和溫度擴(kuò)散層緊緊包裹五邊形冰晶，在五邊形頂點(diǎn)處溶質(zhì)和溫度擴(kuò)散層薄，溶質(zhì)梯度和溫度梯度大，冰晶生長迅速；在五邊形各邊中部溶質(zhì)和溫度擴(kuò)散層厚，溶質(zhì)梯度和溫度梯度小，冰晶生長緩慢。隨著冷凍時(shí)間的推移，冰晶主枝生長排擠的溶質(zhì)和熱量在主枝之間的根部區(qū)域固液界面前沿富集，此區(qū)域的冰晶生長受到抑制，冰晶縮頸現(xiàn)象加劇。繼續(xù)延長冷凍時(shí)間，主枝之間區(qū)域的溶質(zhì)和熱量富集嚴(yán)重，形成半封閉的高溶質(zhì)濃度和高溫區(qū)域；在主枝和側(cè)枝的部分固液混合區(qū)，擴(kuò)散通道暢通，冰晶凝固排擠的溶質(zhì)擴(kuò)散充分，形成低溶質(zhì)濃度固液混合區(qū)。

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

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

4 結(jié)論

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

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