界面厚度对果汁等温结晶冰晶生长影响的相场法模拟

田妍基，陈锦权，陈梅英,3*

(1.宁德职业技术学院，福建 福安 355000；2. 福建农林大学食品科学学院，福建 福州 350002； 3. 福建农林大学管理学院旅游学院，福建 福州 350002)



界面厚度对果汁等温结晶冰晶生长影响的相场法模拟

田妍基1，陈锦权2，陈梅英2,3*

(1.宁德职业技术学院，福建 福安 355000；2. 福建农林大学食品科学学院，福建 福州 350002； 3. 福建农林大学管理学院旅游学院，福建 福州 350002)

将果汁体系视为水和溶质组成的二元系统，采用相场法模拟二元系统的相变微观结构，研究了等温结晶过程中，界面厚度对冰晶生长的形貌及溶质分布的影响，确定模型中界面厚度的合理取值。结果表明：界面厚度对冰晶的形貌有正负影响。在扰动干扰的范围内，二次枝晶随着界面厚度的增大而越来越发达，枝晶尖端也尖锐；在扰动干扰范围外，晶粒不产生二次枝晶，同时主轴枝晶变细，部分主轴枝晶消失。界面厚度取值越大，晶粒生长速率越快，而枝晶尖端半径减小。界面厚度取值在3 dx左右能较好地反应枝晶生长形貌及提高计算效率。

相场法；界面厚度；等温结晶；冰晶生长

当前, 相场法已成为模拟凝固组织、界面形貌的重要方法[1]。相场法[2]的理论基础是金兹堡-朗道相变理论，系统的模拟二元系统的相变微观结构，研究了等温结晶过程中界面厚度对冰晶生长的扩散、有序化势以及热力学驱动力的综合作用，这是一种通过该理论的微分方程来反映的一种计算方法。近年来，国内外研究者利用相场法模拟凝固微观组织经历了从二维→三维[3-4]、 二元合金→多元合金[5]、单晶粒→多晶粒[6-8]、自由枝晶→定向凝固[9-10]、单相场→多相场[11]、无流场→包含流场[12-16]逐步深入的发展历程，从而对凝固过程进行真实的模拟[17]。一般情况下，凝固体微观组织的形成取决于溶质扩散、界面曲率与热扩散[1，18]。迄今关于利用相场法模拟在果汁冷冻浓缩过程中的冰晶生长情况的相关研究较少，目前仍处于初步摸索的阶段[19-20]，亟待进行深入研究。

如今，研究果汁冷冻浓缩尚未构建冰晶生长特征的数学物理模型[21-22]，还是主要以实验为基础，冰晶夹带的影响因素主要有溶液浓度、降温速率、结晶时间等[23-25]，另外，还包括处于非平衡态热力学条件下的动力学规律[18,26]。根据目前已有研究结果表明，研究者何立群等[27]提出：异相成核后，研究呈现分形生长的冰晶形态结构和冰晶的演变过程，并基于分形的相关理论，对水分子的结晶能力与降温的速率、溶液的浓度之间的关系进行分析。陈梅英等[19-20,28-29]首次把相场模型引入食品领域，对冷冻浓缩过程冰晶生长机制进行初步探索，理论分析了冷冻浓缩过冷时间、过冷度等对冰晶生长的影响。

按理界面厚度应取其真实值，但在相场法模拟中需要一个能计算的界面厚度。本研究基于Kim S G等提出的KKS相场模型[26]，结合Allen[30]提出的弥散型界面厚度的定义，探寻能反映真实情况且提高计算效率的界面厚度取值，将果汁的多元成分视为水和溶质二元系统，对二元系统冷冻浓缩等温结晶过程进行相场法模拟，分析冰晶生长过程中界面厚度与冰晶形貌和溶质分布之间的关系。

1 建立冰晶生长的相场模型

1.1 相场法的控制方程

冰晶生长相场模型以自由能减小原理为基础，基于kim S G等[26]提出的KKS相场模型，并结合近似稀溶液时的基础方程，从而建立了果汁在冷冻浓缩过程中，冰晶生长的相场模拟模型。相场和溶质场方程采用自由能密度形式的相场法分别推导出。自由能密度的定义是：利用固相与液相的自由能密度，分别乘以各自的分数，之后再加上剩余自由能的和，可将其表示为：

f(c，φ)=h(φ)fS(cS)+[1-h(φ)]fL(cL)+Wg(φ)

(1)

式(1)中：f为自由能密度，J·m-2；fS表示固相，J·m-2；fL表示液相，J·m-2；且其固相分数为h(φ)，h(φ)=φ3(6φ2-15φ+10)，液相分数为1-h(φ)；c为溶质浓度；φ为相场变量，φ=-1或0时为液相，φ=1时为固相，在液-固的界面上，φ值在-1～1或0～1的范围内发生连续的变化；Wg(φ)表示剩余自由能，其中，W是相场参数，g(φ)是剩余自由能函数，g(φ)=φ2(1-φ)2。可将相场方程表示为：

(2)

式(2)中：M为相场迁移率参数[式(9)]； f(φ)是自由能密度对相场φ的一阶导数；ε为与界面能有关的相场参数。当利用稀溶液来近似时，则表达式为：

(3)

式(3)中： R表示气体常数， R的值为8.31J·kmol-1；T表示温度,K；Vm表示摩尔体积,L·mol-1；ce表示平衡浓度；下标L和S分别表示液相和固相； h′(φ)为h(φ)的导函数， h′(φ)=30φ2(1-φ)2； g′(φ)为 g(φ)的导函数，g′(φ)=2φ(1-φ)(1-2φ)。

1.2 溶质场扩散方程

采用自由能密度的形式，可以将溶质场的扩散方程，改写表示为：

(4)

式(4)中：D(φ) 表示溶质扩散率；fc表示自由能对浓度的一阶导数；fcc表示自由能对浓度的二阶导数。界面区域的溶质浓度c，是指固、液相的摩尔分数之和，在固相和液相达到平衡时，在界面区域中，任意一点的固相和液相的化学势相等，即：

c=h(φ)cS[1-h(φ)]cL

(5)

μS[cS(x,t)]=μL[cL(x,t)]

(6)

式(6)中μL和μS分别为液相和固相的化学势。

1.3 加入扰动

计算时，在界面处施加微小的扰动也会影响真实的枝晶形貌，在相场方程中加入人为的随机扰动：

(7)

式中χ在-1到1之间随机取值， ϖ是与时间相关的相扰动强度因子，16g(φ)用以控制在固液界面中出现扰动，最大扰动可能出现在φ=0.5处，若远离界面则扰动迅速减小。

2 确定相场参数及模拟计算

2.1 参数的确定

相场参数ε和W与界面厚度及界面能都有关系，相场迁移率参数M则与界面动力学系数有关，它们可以表示为：

(8)

(9)

式(8)和(9)中：σ表示固液两相的界面能，J·m-2；R表示气体常数；T表示温度；Vm表示摩尔体积；ke表示平衡常数；μk表示动力系数；me表示液相线斜率；ε表示与界面能有关的参数, λ表示界面厚度，μm。

2.2 物性参数的取值

从已有的研究[31-34]得知，冰晶结构在一般情况下，均表现为密排的六方结构。模拟冰晶结构各向异性模数k的取值为6。当控制其他相关的物性参数取值不发生改变时，只改变界面厚度λ，观察及分析界面厚度对冰晶形貌的影响。选用质量分数为9.6%的糖水溶液视为水和溶质二元系统，将其作为研究对象。在进行计算时，所使用的物性参数，具体如表1所示。

表1 蔗糖水溶液的物性参数

注：表中数据参照文献[35]计算得到。

2.3 数值计算

(1)边界条件与初始条件。当计算界面区域的边界时，相场与溶质场的边界条件都取绝热边界条件。r表示初始晶核半径，则

(10)

式(10)中： x是模拟中的横坐标；y是模拟中的纵坐标；T表示有量纲温度，Tm表示模拟的初始温度，ΔT表示过冷度。

(2)数值计算方法。相场方程利用显示有限差分求解，时间步长Δt受到溶质场计算的限制，即Δt<(Δx)2/4DL，DL表示液相中的溶质扩散系数，取值为：Δt=1.5×10-9s。晶体轴生长主轴与直角坐标系的x轴和y轴相对应。模拟计算网格数为1400×1400，空间步长为5×10-8m，初始晶核为半径r=10个网格数的球。

3 模拟结果与分析

在相变过程中, 固体与液体之间存在明显分界面。将固液界面厚度扩大至可计算的尺度，在不同的界面厚度下利用相场法进行模拟，研究界面厚度这个因素对冰晶形貌与其冰晶溶质分布的影响情况，这样就能获得稳定的界面推进速度。

3.1 界面厚度对冰晶形貌的影响

图1所示为不同界面厚度下的冰晶生长形貌。当λ=1.5dx时，晶粒形貌呈现等轴枝晶形貌，且二次分支不明显；随着界面厚度的增加，当λ=3dx时随着界面厚度的增加，晶粒开始出现分支，二次分支更加发达；当λ=4dx时枝晶的主干变细，二次枝晶生长发达，枝晶整体形貌呈六边形生长。当λ=4.5dx及λ=5dx时，界面厚度进一步增大，主轴枝晶变得更细，部分主轴枝晶消失，枝晶整体形貌呈球形。这是因为扰动对枝晶的影响存在一定的范围。在扰动的影响范围内，当界面厚度较小的时候，受到扰动的影响较小，二次枝晶不明显；随着界面厚度的增大，扰动的影响随之增大，晶体的界面前沿不稳定，二次枝晶开始产生，并随着界面厚度的增加愈加发达。当界面厚度增大到一定程度时，能有效地限制扰动的对晶粒的影响减弱，二次枝晶减少最终晶粒没有二次枝晶的产生，如图1-d、e。这种现象分别与于艳梅[36]、王智平[37]等模拟的结果相吻合。

3.2 界面厚度对冰晶溶质分布的影响

图2 显示了不同界面厚度条件下的溶质分布情况,图1和图2对比可知相场和溶质场相吻合。从述两幅图中可以发现，枝晶的固相中夹带着部分的溶质，凝固界面的溶质浓度比周围环境低，且枝晶前沿富集了大量的溶质。这是溶质扩散的过程。一方面水分子和溶质分子都在不断运动中，水分子之间容易形成氢键而结晶，同时溶质分子与水分子相互作用，造成了夹带；另外一方面，在结晶过程中，水结成冰是一个膨胀的过程，因为冰晶生长的速率大于溶质扩散的速率，所以在溶质来不及充分扩散之前，冰晶内的浓缩溶质会被迅速挤向尖端，从而在枝晶前端富集，形成一股流行尖端的溶质流。由图2可以发现，随着界面厚度的逐渐增大，枝晶的生长速率逐渐加快，而枝晶尖端半径则逐渐减小。同时，从方程(9)可以看出，随着界面厚度取值的增大，界面的迁移速率减小，溶质的扩散较为充分，溶质梯度较少，枝晶生长凝固时间越短，所以生长速度越快。所以，方程(9)的理论结果与相场模拟结果吻合。

由图1和图2可知，当界面厚度取值过大时，枝晶生长形貌变异。所以采用相场法模拟冰晶的生长时，界面厚度取值应在3dx左右，这样不仅可以得到较可靠的冰晶生长形貌，又可以提高计算机的模拟效率。

4 讨论与结论

4.1 将果汁视作水和溶质的二元系统，利用相场法模拟二元系统冷冻浓缩过程中的冰晶生长情况，探讨界面厚度对冰晶形貌的影响。界面厚度对晶粒的形貌有正负影响，当界面厚度尺度达到一定取值时，其数值误差随之被放大，从而在界面处形成误差噪音，误差噪音在界面上形成明显的扰动，使计算得到的冰晶形态发生变化。在扰动干扰的范围内，随着界面厚度的增大,二次枝晶越来越发达，枝晶尖端也尖锐；在扰动干扰范围外，晶粒不产生二次枝晶，同时主轴枝晶变细，部分主轴枝晶消失。由于误差噪音的性质难以描述，故利用相场法模拟冰晶生长时，界面厚度大小应有适当的取值，从而达到有效控制误差噪音。

4.2 相场和溶质场相吻合，溶质分布情况随着界面厚度发生变化。随着界面厚度取值越大，界面的迁移速率越小，导致溶质的扩散速度与冰晶的生长速度相差较小，溶质扩散有相对足够的时间，避免大量溶质聚集在固-液界面前沿，削弱了界面前沿出现溶质截留现象。因此溶质的扩散较为充分，溶度梯度较少，枝晶生长凝固时间越短，所以生长速度越快。

4.3 晶粒生长速率和枝晶尖端半径受界面厚度取值的影响，界面厚度取值越大，晶粒生长速率越快，而枝晶尖端半径减小。刚开始，晶粒生长速度呈指数上升，且随着界面厚度取值的增大而迅速上升。这是由于开始结晶的时候伴随着能量起伏波动，从而影响初始速度，界面厚度较小的时候，扰动对其的影响较少，在扰动干扰的范围内，随着界面厚度的增大，二次枝晶越来越发达，枝晶尖端半径减小。因此界面厚度的取值应在3 dx左右，较能真实反映晶粒的生长形貌；界面厚度取值不恰当，将使冰晶形貌发生变异，也会受到误差噪音的影响。

4.4 在果汁冷冻浓缩过程中，忽略冰晶结晶时释放的潜热，本研究采用相场法对冰晶生长进行等温模拟，但模拟的结果很大程度上受到参数取值的影响。今后将继续在探讨热扩散率、初始晶核半径等其他他相关参数取值方面努力，以取得更优的模拟结果。

[1]HUA HOU，YUHONG ZHAO，YUHUI ZHAO．Simulation of the precipitation process of ordered inter metallic compounds in binary and ternary Ni- Al based alloys by the phase-field model[J]．Materials Science and Engineering，2009，499：204- 207.

[2]王洁玉，陈长乐，陈志，等．相场法模拟 Ni-Cu 合金枝晶生长中过冷度的影响[J]．铸造技术，2007，28(4)：538-540.

[3]KARMA A, RAPPEL W J. Quantitative phase-field modeling of dendritic growth in two and three dimensions[J]．Phys Rev E,1998, 57(4)： 4323-4349.

[4]朱昌盛, 冯力, 王智平, 等．三维枝晶生长的相场法数值模拟研究[J]．物理学报，2009，58(11)：8055-8061.

[5]ZHANG R J，JING T，JIE W Q，et al．Phase-field simulation of solidification in multi-component alloys coupled with thermodynamic and diffusion mobility databases[J]．Acta Mater，2006，54(8)：2235-2239.

[6]LI M E，XIAO Z Y，YANG G C，et al．Anisotropic growth of multigrain in equiaxial solidification simulated with the phase field method[J]．Chin Phys Soc，2006，15(1)：219-223.

[7]SUN Q，ZHANG Y T，CUI H X，et al．Phase field modeling of multiple dendrite growth of Al-Si binary alloy underisothermal solidification[J]．China Foundry，2008, 5(4)：265-267.

[8]FENG L，WANG Z P，ZHU C S，et al．Phase-field model of isothermal solidification with multiple grain growth[J]．Chin Phys B，2009，18(5)：1985-1990.

[9]TAKAKI T，FUKUOKA T，TOMITA Y．Phase-field simulation during directional solidification of a binary alloy using adaptive finite element method[J]．J Cryst Growth，2005，283(1/2)：263-278.

[10]WANG Z J，WANG J C，YANG G C．Phase field investigation on the initial planar instability with surface tension anisotropy during directional solidification of binary alloys[J]．Chin Phys B，2010，19(1)：492-496.

[11]TIADEN J，NESTLER B，DIEPERS H J，et al．The multiphase-field model with an integrated concept for modeling solute diffusion[J]．Physica D，1998，115(1/2)：73-86.

[12]TONG X，BECKERMANN C，KARMA A．Velocity and shape selection of dendritic crystals in a forced flow[J]．Phys Rev E，2000，61(1)：49-52.

[13]CHEN Z，CHEN C L，HAO L M．Numerical simulation of succinonitrite dendritic growth in a forced flow[J]．Acta Metall Sin, 2008，21(6)：444-450.

[14]龙文元，吕冬兰，夏春，等．强迫对流影响二元合金非等温凝固枝晶生长的相场法模拟[J]．物理学报．2009.58(11)：7802-7808.

[15]袁训锋，丁雨田，郭廷彪，等．对流作用下枝晶生长行为的相场法[J]．中国有色金属学报，2010，20(4)：681-687.

[16]袁训锋，丁雨田，郭廷彪，等．强制对流作用下镁合金枝晶生长的相场法数值模拟[J]．中国有色金属学报，2010，20(8)：1474-1480.

[17]余德洋，刘宝林，王伯春．超声波对蔗糖溶液中树枝冰晶生长速度的影响[J]．食品与发酵业，2011，37(8)：92-94.

[18]李刚，刘新田．相场模型及其在凝固组织模拟中的研究进展[J]．铸造技术，2005，26(10)：974-976.

[19]陈梅英，陈锦权，陈永雪．荔枝汁冷冻浓缩冰晶生长的动力学分析和相场模拟[J]．中国农业大学学报：自然科学版，2011，16(2)：153-157.

[20]陈梅英，陈永雪，王文成，等． 相场法模拟[J]．福建农林大学学报：自然科学版，2010，39(5)：548-551.

[21]KIM S G，KIM W T，SUZUKI T．Phase-field model for bianry alloys[J]．Phys Rev E，1999，60(6)：71-86.

[22]吕福扣，刘宝林，李维杰．HA纳米微粒对PEG-600低温保护剂反玻璃化结晶的影响[J]．低温物理学报，2012，34 (4)：315-320.

[23]LEVENT B．Mathematical analysis of freeze concentration of apple Juice[J]．Journal of Food Engineering，1993(19)：95-107.

[24]肖旭霖，李慧．苹果汁冷冻浓缩工艺的研究[J]．农业工程学报，2006，22(1)：192-194.

[25]方婷，陈锦权，唐凌，等．橙汁冷冻浓缩动力学模型的研究[J]．农业工程学报，2008，24(12)：243-248.

[26]GUILLERMO PETZOLD，JOSÉ M．AGUILERA．Ice Morphology：Fundamentals and technological applications in foods[J]．Food，2009，4(4)：378-396.

[27]何立群，张永锋，罗大为，等．生命材料低温保护剂溶液二维降温结晶过程中的分形特征[J]．自然科学进展，2002，12(11)：1167-1171.

[28]陈梅英，陈永雪，陈威，等．过冷度对冷冻浓缩过程冰晶生长的宏微观研究[J]．西南大学学报：自然科学版，2010，32(5)：140-145.

[29]陈梅英，王文成，陈锦权．相场法模拟冷冻浓缩过程冰晶生长的可行性探讨[J]．江西农业学报，2010，22(3):137-139.

[30]ALLEN S M，CAHN J W．A microscopic theory for anti-phase boundary motion and its application to anti-phase domain coarsening[J]．Acta Metall，1979，27(6)：1085-1092.

[31]陶乐仁，华泽钊．低温保存时的气泡形成现象试验研究[J]．食品科学，2000，21(6)：12-15.

[32]陶乐仁，华泽钊．低温保护剂溶液结晶过程的显微实验研究[J]．工程热物理学报，2001，22(4)：481-484.

[33]刘占杰，华泽钊，陶乐仁．脂质体悬浮液结晶对其冻干品质影响的研究[J]．青岛海洋大学学报，2001，31(7)：612-618.

[34]陶乐仁．水溶液冻结过程的显微实验研究[D]．上海：上海理工大学，2000.

[35]李建国．凝固原理：第4版[M]．北京：高等教育出版社，2013：171.

[36]于艳梅，杨根仓，赵达文，等．相场法模拟过冷熔体枝晶生长的界面厚度参数的取值[J]．自然科学进展，2001，11(11)：1192-1197.

[37]王智平，张殿喜，石可伟，等．多元合金等温凝固相场法模拟[J]．兰州理工大学学报，2008，34(6)：1-4.

(责任编辑：林海清)

Phase-field Simulation of Ice Crystal Growth by Interfacial Thickness Influence on Juice Isothermal Crystallization Process

TIAN Yan-ji1，CHEN Jin-quan2，CHEN Mei-ying2,3*

(1.NingdeVocationalandTechnicalCollege,Fu′an,Fujian355000,China； 2.CollegeofFoodScience,FujianAgricultureandForestryUniversity,Fuzhou,Fujian350002,China；3.ManagementCollege&Tourismcollege,FujianAgricultureandForestryUniversity,Fuzhou,Fujian350002,China)

The phase-field model which coupled the concentration field and phase field is applied to simulate microstructure evolution during crystal growth of juice binary in 2-D, exploring the process of fruit juice crystallization. In this paper, the influence of interfacial thickness on ice crystal growth morphology and solute distribution was investigated, and a reasonable value of interfacial thickness was estimated. The simulated results showed that the interfacial thickness had positive and negative effects on ice crystal morphology.In the range of disturbance, the secondary dendrite became more and more developed and the dendrite tip was also more sharper with the increase of the interfacial thickness. Outside the range of disturbance, the secondary dendrite didn't grow, and the principal axis of the dendrite became thin and even some disappeared. The growth speed of ice crystal would increase and the tip radius would decreased with the increase of interfacial thickness.When the interfacial thickness was around 3dx, the best simulation result of the dendrite morphology and high calculation efficiency could be obtained.

phase-field；interfacial thickness；isothermal crystallization；crystal growth

2016-05-12初稿；2016-06-05修改稿

田妍基(1981-)，女，硕士，讲师，主要从事食品加工与营养研究(E-mail：tyj6285268@126.com)

*通讯作者：陈梅英(1972-)，女，副教授，博士，主要从事食品工程与旅游食品安全研究(E-mail：cmy2816@126.com)

国家自然科学基金项目(31101327)；高等学校博士学科点专项科研基金博导类资助项目(20123515110016)；食品生物技术校级示范专业和校级优秀教学团队建设项目(2015)

TS 275.5

A

1008-0384(2016)09-975-06

田妍基，陈锦权，陈梅英.界面厚度对果汁等温结晶冰晶生长影响的相场法模拟[J].福建农业学报，2016，31(9)：975-980.

TIAN Y-J，CHEN J-Q，CHEN M-Y.Phase-field Simulation of Ice Crystal Growth by Interfacial Thickness Influence on Juice Isothermal Crystallization Process[J].FujianJournalofAgriculturalSciences，2016，31(9)：975-980.