强界面能各向异性下二元Ni-Cu合金枝晶生长过程的相场法模拟

袁训锋, 丁雨田

强界面能各向异性下二元Ni-Cu合金枝晶生长过程的相场法模拟

袁训锋, 丁雨田

(兰州理工大学 甘肃省有色金属新材料省部共建国家重点实验室，兰州 730050)

基于Wheeler模型和Eggleston修正强界面能各向异性的方法，建立耦合溶质场和温度场的相场模型，模拟强界面能各向异性下Ni-Cu合金枝晶生长过程。结果表明：在强界面能各向异性作用下，界面方向枝晶生长不连续且枝晶出现棱角；由于枝晶尖端温度梯度和溶质梯度较大，枝晶生长迅速。当界面能各向异性强度低于临界值时，枝晶尖端生长速度随界面能各向异性强度的增加而增大；当界面能各向异性强度等于临界值时，枝晶尖端生长速度下降4.34%；当界面能各向异性强度大于临界值时，枝晶尖端生长速度随界面能各向异性强度的增加先增大到极大值后逐渐减小。当无量纲热过冷度较小时，晶体平衡形貌为类矩形；随着无量纲热过冷度的增加，晶体平衡形貌向枝晶转变，枝晶尖端生长速度先呈幂指数增加，然后呈线性增加，枝晶生长由扩散控制转变为动力学控制。

Ni-Cu合金；界面能；各向异性强度；过冷度；相场法；枝晶生长

枝晶的形成过程是一个涉及热量、质量和动量传输以及界面动力学和毛细作用效应相耦合的自由边界问题。理论分析方法会遇到巨大的数学困难，传统试验方法也无法揭示枝晶结构形成的机理。随着计算科学的发展和凝固理论的完善，数值模拟技术克服了理论分析和实验研究的困难，已成为研究凝固过程中枝晶生长的重要手段。

相场方法将界面曲率、各向异性和动力学隐含在相场方程中，只需要求解相场方程就可以得到固−液界面的状态、曲率及界面的移动，从而避免跟踪固−液界面[1]。LANGER[2]提出采用扩散界面模型研究凝固现象。KOBAYASHI[3]利用含有各向异性的相场模型，实现了对纯金属过冷熔体中枝晶生长过程的模拟。WARREN 和 BOETTINGER[4]及 SUZUKI等[5]通过对界面区域自由能密度函数的不同定义，分别建立两种不同描述合金凝固中枝晶生长的相场模型，即 WBM和KKS模型，从此开始了二元合金枝晶生长的相场法研究。近20年来，国内外研究者利用相场法模拟凝固微观组织经历了从二维到三维[6−7]、 从二元合金到多元合金[8]、从单晶粒到多晶粒[9−11]、从自由枝晶到定向凝固[12−13]、从单相场到多相场[14]、从无流场到包含流场[15−19]逐步深入的发展历程，使得模拟结果越来越接近真实凝固过程，并已取得了显著成果。晶体原子按照特定对称性排列，导致固−液界面能和界面动力学具有不同程度的各向异性。在枝晶生长过程中液−固界面各向异性不仅决定枝晶生长方向，还在很大程度上影响枝晶的生长行为[20]。在相场模型中，一般通过引入与界面能有关的因子 η(θ)=1+γcos(4θ)来描述界面能各向异性。其中，θ为界面法向与X正方向之间的夹角；γ为界面能各向异性强度。界面能各向异性对晶体平衡形状的影响可以用 Gibbs-Thomson方程描述：η(θ)+ η′(θ)=fL− fS=1−15γcos (4θ)。其中，fL和 fS分别为液相和固相的自由能密度。在弱界面能各向异性时(γ＜1/15)，方程两边都为正，枝晶形貌光滑连续；在强界面能各向异性时(γ＞1/15)，方程左边为负，导致出现缺失方向，界面变得不连续，凹下去的部分出现耳子。2001年，EGGLESTON等[21]提出了修正强界面能各向异性的方法，从而解决了晶向缺失问题，建立了强界面能各向异性下的枝晶生长相场模型。KIM等[22]采用此相场模型研究了强界面能各向异性下过冷纯Si枝晶生长特性。李俊杰等[23]、ZHANG等[24]、CHEN等[25−26]采用修正的相场模型，对强界面能各向异性下的枝晶生长行为进行了研究。但是，这些研究都是在Karma模型基础上进行的，此模型中部分相场参数是强制假定的，与实际情况相差较大。而Wheeler模型中的相场参数与实际物理参数相关联，更接近真实物理过程。

本文作者基于Wheeler模型和Eggleston修正强界面能各向异性的方法，建立耦合溶质场和温度场的相场模型，模拟强界面能各向异性作用下Ni-Cu合金枝晶的生长过程，研究界面能各向异性强度和过冷度对枝晶形貌、尖端稳态生长速度、温度和溶质浓度的影响。

1 相场模型

相场法通过引入一个相场变量φ 来描述系统中各点的物理状态。在研究中，φ=0代表固相，φ =1代表液相，在固/液界面上φ在0～1之间连续变化。

1.1 相场控制方程

采用 Eggleston方法[21]修正缺失方向得到的界面能各向异性因子如下：

式中：i=0～3；θm为晶体平衡形貌出现耳子时的角度，可以通过方程(2)求得：

基于Wheeler相场模型[18]和Eggleston修正方法[21]建立强界面能各向异性作用下枝晶生长的修正相场方程：

式中：温度u、时间t及距离X和Y均为无量纲变量，u=(T−Tm)/(Tm−T0)； t=t′/(ω2/κ)； X=X′/ω ， Y=Y′/ω ；Ω=cp[Ωt+mL(x−x0)]/L为无量纲过冷度；α=( 2 ωL2)/(12cpσTm)为系统的物理参数；m=μσTm/(κL)为界面动力学系数；ε=δ/ω为与界面层厚度有关的参数；T为温度；Tm为熔点；T0为系统初始温度；Ωt为无量纲热过冷度；X′和Y′分别为距离；δ为界面层厚度；ω为参考长度；t′为时间；κ为热扩散率；μ为界面迁移率；ml为液相线斜率；cp为定压比热容；L为单位体积的结晶潜热；σ为界面能；x0为过冷熔体初始浓度(摩尔分数)；x为熔体实际浓度(摩尔分数)。

1.2 温度场控制方程

温度场控制方程为

式中：p′(φ )为势函数 p′(φ )=φ2(10−15 φ +6 φ2)对φ的导数。

1.3 溶质场控制方程

溶质场控制方程为

式中：D′为有效溶质扩散系数。若 Ds和 Dl分别表示固相和液相的溶质扩散系数，k0为溶质平衡分配系数，则D′可表示为

2 计算方法

2.1 相场参数及材料物性参数

选择Ni-Cu合金作为模拟对象，结合Ni-Cu合金的物性参数，使用相场法对过冷合金熔体中强界面能各向异性作用下的枝晶生长行为进行模拟。Ni-Cu合金的热物性参数及模拟参数如下：Tm=1 594.5 K，L=2 100 J/cm3，cp= 4.83 J/(cm3·K)，σ=3.38×10−5J/cm2，κ=0.27 cm2/s，Ds=1.0×10−9cm2/s，Dl=1.0×10−5cm2/s，μ=175.3 cm/(K·s)，ω=0.21 µm，ml=0.821 1，k0=0.86，α=419，m=0.017，ε=0.005，γ=0.1，Ωt=0.1，θm=0.378，x0=40.38%，空间步长 ΔX=ΔY=0.005，时间步长Δt=2.0×10−5。

2.2 数值计算方法

对传输方程(6)采用交替隐式格式求解，控制方程(3)～(5)和(7)，同时采用显示差分格式求解。由于温度场方程所采用的交替隐式格式具有任意的稳定性，整个计算过程的稳定性受以下条件约束：

式中：m′=max(Dl，m)；k′为考虑方程中非线性项而设计的修正系数，一般取1～2。

2.3 初始条件和边界条件

开始计算时，假设网格数为R的球形初始晶核位于充满均匀过冷合金熔体的中心(X0, Y0)，在球形区域内节点的相场变量初始值为0，其余节点为1，即

在计算区域边界时，相场、温度场和溶质场的采用Neumann条件。

3 结果与分析

3.1 强界面能各向异性作用下枝晶生长

图1(a)、(b)和(c)所示分别为强界面能各向异性下γ=0.1及Ωt=0.1时枝晶生长相场、温度场和溶质场的模拟结果。模拟计算网格数为800×800，图1所示为模拟中心400×400区域的模拟结果。可以看出，在强界面能各向异性下，晶体以枝晶方式生长，4个一次枝晶沿〈100〉方向生长，互成90°夹角，枝晶主干细长，各向异性明显。由于某些方向的生长消失，界面方向不连续，枝晶主枝尖端出现棱角，一次枝晶根部出现明显颈缩现象。

从图1(b)和(c)可以看出，温度和溶质的分布情况与枝晶生长相吻合。枝晶尖端深入熔体，热量和溶质扩散较快，温度梯度和溶质浓度梯度大，枝晶生长迅速；在枝晶根部，由于溶质再分配析出的溶质堆积和凝固释放潜热的富集引起过冷度减小，枝晶生长缓慢。此外，先凝固的枝晶主干中心位置溶质浓度最低，后凝固的位置溶质浓度较高。

图1 强界面能各向异性下γ=0.1、Ωt =0.1时的枝晶生长模拟结果Fig.1 Simulation results of dendrite growth under anisotropy of strong interface energy at γ=0.1 and Ωt =0.1: (a) Phase field;(b) Temperature field; (c) Solute field

3.2 界面能各向异性强度对枝晶生长的影响

界面能各向异性强度表示合金熔体凝固过程中界面的张力、厚度及动力学各向异性的程度。 准确理解晶体生长界面能各向异性强度是理解和明晰枝晶生长理论的关键。但到目前为止，界面能各向异性强度与实际物性参数的关系无法用实验方法准确获得。为此，本文作者将研究界面能各向异性强度与凝固过程中枝晶生长行为的关系，为确定界面能各向异性强度与实际物性参数的关系提供理论依据。

图2所示为在Ωt=0.1及t=30 000Δt时不同界面能各向异性强度下的枝晶生长形貌。从图2可以看出，当界面能各向异性强度γ为0.060时，界面所有生长方向稳定，晶体界面方向连续，枝晶光滑；随着界面能各向异性强度的增加，枝晶主干变得细长，枝晶尖端生长速度增加，曲率半径减小；当界面能各向异性强度增加到0.068时，枝晶某些方向的生长开始出现缺失，界面方向不连续，枝晶尖端出现棱角。

为了定量分析界面能各向异性强度对枝晶生长行为的影响，计算了枝晶尖端稳态生长速度、尖端温度和溶质浓度与界面能各向异性强度的关系，结果如图3所示。

从图 3可以看出，当界面能各向异性强度小于0.065时， 随着界面能各向异性强度的增加，枝晶尖端生长速度迅速增加，因为枝晶尖端的温度和溶质扩散层变薄，使热扩散和溶质扩散变得容易，尖端温度和溶质浓度降低；当界面能各向异性强度接近临界值0.067时，枝晶尖端生长速度继续增加达到最大值，由于凝固析出的大量溶质富集在固−液界面前沿，溶质扩散速度小于枝晶生长速度，出现溶质截留现象，枝晶尖端溶质浓度增加，而枝晶尖端温度逐渐趋于稳定；当界面能各向异性强度等于临界值0.067时，枝晶尖端部分生长方向出现缺失，枝晶尖端生长速度和溶质浓度分别下降4.34%和0.002 7%，枝晶尖端温度变化不大；当界面能各向异性强度大于临界值时，随着界面能各向异性强度的增加，枝晶尖端生长速度先增大到极大值后逐渐减小，尖端温度和溶质浓度缓慢降低。

3.3 无量纲热过冷度对枝晶生长的影响

过冷度对枝晶生长具有决定性的影响。晶核生长过程需要克服能垒，只有在过冷条件下才能实现枝晶生长过程，且过冷度的大小直接决定晶核的生长方式。在此，对无量纲热过冷度为0.01～0.125时的Ni-Cu合金枝晶生长形貌及尖端稳态行为进行研究。

图4所示为γ=0.1及t=30 000Δt时，不同无量纲热过冷度下枝晶生长形貌。从图4可以看出，当无量纲热过冷度较小时，晶体优先生长方向上没有形成突出的尖端，优先生长方向之间的界面近似为平面，晶体平衡形貌为类矩形；随着无量纲热过冷度的增大，优先生长方向上的生长速度增加，优先生长方向之间的界面下凹形成曲面，尖端形成枝晶形貌；从类矩形形状向枝晶形貌转变时的无量纲热过冷度为 0.05(即 22 K)。

图2 不同界面能各向异性强度下的枝晶生长形貌Fig.2 Dendrite morphologies at different anisotropy strengths of interface energy: (a) γ=0.060; (b) γ=0.065; (c) γ=0.068; (d)γ=0.100

图3 界面能各向异性强度对枝晶尖端生长行为的影响Fig.3 Effect of interface energy anisotropy strength γ on dendrite growth: (a) Dendrite tip velocity; (b) Dendrite tip temperature and solute concentration

图4 γ=0.1及t=30 000 Δt时不同无量纲热过冷度下的枝晶生长形貌Fig.4 Dendrite morphologies at γ=0.1, t=30 000 Δt and different dimensionless thermal supercooling degrees: (a) Ωt =0.01; (b) Ωt=0.03; (c) Ωt =0.05; (d) Ωt =0.1

为了定量分析无量纲热过冷度对枝晶生长行为的影响，考察枝晶尖端稳态生长速度、尖端温度和溶质浓度与无量纲热过冷度的关系，结果如图5所示。

从图5可以看出，当无量纲热过冷度小于0.03时，随着无量纲热过冷度的增加，枝晶尖端生长速度呈幂指数增加，由于枝晶尖端的温度和溶质扩散层较厚，枝晶生长释放的潜热和溶质再分配析出的溶质不易扩散，使其在界面前沿富集，尖端温度和溶质浓度增加；随着无量纲热过冷度的继续增大，快速生长的枝晶尖端深入过冷熔体中，热扩散加快，枝晶尖端温度扩散层变薄，枝晶尖端的温度降低；当无量纲热过冷度增加到0.05时，优先生长的尖端形成枝晶形貌，枝晶尖端生长速度由呈幂指数增加转变为呈线性增加，溶质扩散速度小于枝晶生长速度，出现溶质截留现象，枝晶尖端溶质浓度增加、温度降低，枝晶生长由扩散控制转变到动力学控制。

模拟结果显示：当无量纲热过冷度低于 0.05(即22 K)时，晶粒以类矩形生长，生长速度与过冷度之间呈幂函数关系；而当无量纲热过冷度高于0.05时，晶粒以枝晶生长，生长速度与过冷度呈线性关系，这种转变规律与 WILLNECKER等[27]的实验结果基本符合。

图5 无量纲热过冷度Ωt对枝晶尖端生长行为的影响Fig.5 Effect of dimensionless thermal supercooling degree on dendrite tip growth: (a) Dendrite tip velocity; (b) Dendrite tip temperature and solute concentration

4 结论

1) 基于Wheeler模型和Eggleston修正强界面能各向异性的方法，建立耦合溶质场和温度场的相场模型，研究二元合金强界面能各向异性下的枝晶生长过程。

2) 在强界面能各向异性作用下，界面方向变化不连续且枝晶出现棱角，枝晶尖端温度梯度和溶质梯度较大，枝晶生长迅速，枝晶主干中心溶质浓度较低。

3) 当界面能各向异性强度低于临界值时，随着界面能各向异性强度的增加，枝晶尖端生长速度增大，且在接近临界值处达到最大值；当界面能各向异性强度等于临界值时，枝晶尖端生长速度下降 4.34%；当界面能各向异性强度大于临界值时，随着界面能各向异性强度的增加，枝晶尖端生长速度先增大到极大值后逐渐减小。

4) 当无量纲热过冷度较小时，晶体平衡形貌为类矩形；随着无量纲热过冷度的增加，晶体平衡形貌向枝晶转变，枝晶尖端溶质浓度单调增加，尖端生长速度先呈幂指数增加后呈线性增加，枝晶生长由扩散控制转变为动力学控制。

REFERENCES

[1] TAN L J, ZABARAS N. A level set simulation of dendritic solidification with combined features of front-tracking and fixed-domain methods[J]. J Comput Phys, 2006, 211(1): 36−63.

[2] LANGER J S. Direction in condensed matter physics[M].Singapore: World Science, 1986: 164.

[3] KOBAYASHI R. Modeling and numerical simulations of dendritic crystal growth[J]. Physica D, 1993, 63(3/4): 410−423.

[4] WARREN J A, BOETTINGER W J. Prediction of dendritic growth and micro-segregation patterns in a binary alloy using the phase-field method[J]. Acta Metall Mater, 1995, 43(2): 689−703.[5] SUZUKI T, ODE M, KIM S G. Phase-field model of dendritic growth[J]. J Cryst Growth, 2002, 237(1): 125−131.

[6] 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.

[7] 朱昌盛, 冯 力, 王智平, 肖荣振. 三维枝晶生长的相场法数值模拟研究[J]. 物理学报, 2009, 58(11): 8055−8061.ZHU Chang-sheng, FENG Li, WANG Zhi-ping, XIAO Rong-zhen. Numerical simulation of three-dimensional dendritic growth using phase-field method[J]. Acta Phys Sin, 2009, 58(11):8055−8061.

[8] ZHANG R J, JING T, JIE W Q, LIU B C. Phase-field simulation of solidification in multi-component alloys coupled with thermodynamic and diffusion mobility databases[J]. Acta Mater,2006, 54(8): 2235−2239.

[9] LI M E, XIAO Z Y, YANG G C, ZHOU Y H. Anisotropic growth of multigrain in equiaxial solidification simulated with the phase field method[J]. Chin Phys Soc, 2006, 15(1): 219−223.

[10] SUN Q, ZHANG Y T, CUI H X, WANG C Z. Phase field modeling of multiple dendrite growth of Al-Si binary alloy under isothermal solidification[J]. China Foundry, 2008, 5(4):265−267.

[11] FENG L, WANG Z P, ZHU C S, LU Y. Phase-field model of isothermal solidification with multiple grain growth[J]. Chin Phys B, 2009, 18(5): 1985−1990.

[12] 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.

[13] 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): 017305.

[14] TIADEN J, NESTLER B, DIEPERS H J, STEINBACH I. The multiphase-field model with an integrated concept for modeling solute diffusion[J]. Physica D, 1998, 115(1/2): 73−86.

[15] 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.

[16] 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.

[17] 龙文元, 吕冬兰, 夏 春, 潘美满, 蔡启舟, 陈立亮. 强迫对流影响二元合金非等温凝固枝晶生长的相场法模拟[J]. 物理学报, 2009, 58(11): 7802−7808.LONG Wen-yuan, LÜ Dong-lan, XIA Chun, PAN Mei-man, CAI Qi-zhou, CHEN Li-liang. Phase field simulation of nonisothermal solidification dendrite growth of binary alloy under the force flow[J]. Acta Phys Sin, 2009, 58(11): 7802−7808.

[18] 袁训锋, 丁雨田, 郭廷彪, 胡 勇, 唐向前. 对流作用下枝晶生长行为的相场法[J]. 中国有色金属学报, 2010, 20(4):681−687.YUAN Xun-feng, DING Yu-tian, GUO Ting-biao, HU Yong,TANG Xiang-qian. Phase-field method of dendritic growth under convection[J]. The Chinese Journal of Nonferrous Metals,2010, 20(4): 681−687.

[19] 袁训锋, 丁雨田, 郭廷彪, 胡 勇. 强制对流作用下镁合金枝晶生长的相场法数值模拟[J]. 中国有色金属学报, 2010, 20(8):1474−1480.YUAN Xun-feng, DING Yu-tian, GUO Ting-biao, HU Yong.Numerical simulation of dendritic growth of magnesium alloys using phase field method under forced flow[J]. The Chinese Journal of Nonferrous Metals, 2010, 20(8): 1474−1480.

[20] LANGER J S. Existence of needle crystals in local models of solidification[J]. Phys Rev A, 1986, 33(1): 435−441.

[21] EGGLESTON J J, MCFADDEN G B, VOORHEES P W. A phase-field model for highly anisotropic interfacial energy[J].Physica D, 2001, 150(1/2): 91−103.

[22] KIM S G, KIM W T. Phase field modeling of dendrite growth with high anisotropy[J]. J Cryst Growth, 2005, 275(1/2):355−360.

[23] 李俊杰, 王锦程, 杨根仓. 相场法模拟界面能各向异性对枝晶生长行为的影响[J]. 自然科学进展, 2005, 15(11):1312−1317.LI Jun-jie, WANG Jin-cheng, YANG Gen-cang. Phase field modeling of dendritic growth with high interfacial energy anisotropy[J]. Progress in Natural Science, 2005, 15(11):1312−1317.

[24] ZHANG G W, HOU H , CHENG J. Phase field model for strong anisotropy of kinetic and highly anisotropic interfacial energy[J].Transactions of Nonferrous Metals Society of China, 2006,16(s2): 307−313.

[25] CHEN Z, CHEN C L, HAO L M. Numerical simulation of facet dendrite growth[J]. Transactions of Nonferrous Metals Society of China, 2008, 18(4): 938−943.

[26] CHEN Z, CHEN C L, HAO L M. Numerical simulation of facet dendritic growth in a forced flow[J]. Can J Phys, 2009, 87(2):117−123.

[27] WILLNECKER R, HERLACH D M, FEUERBACHER B.Evidence of noneqilibrium processes in rapid solidification of undercooled melts[J]. Phy Rev Lett, 1989, 62(23): 2707−2710.

Phase-field simulation of dendrite growth process for binary Ni-Cu alloy with anisotropy of strong interface energy

YUAN Xun-feng, DING Yu-tian
(State Key Laboratory of Gansu Advanced Nonferrous Metal Materials,Lanzhou University of Technology, Lanzhou 730050, China)

Based on the Wheeler model and the Eggleston regularization technique of the anisotropy of strong interface energy, the phase-field model was built by coupling the solute field with temperature field, and the dendrite growth process of the anisotropy of strong interface energy of Ni-Cu alloy was simulated. The results show that the variation of interface orientation discontinuous and the corners form on the dendrite with the anisotropy of strong interface energy,the dendrites grow fast due to the large temperature gradient and solute gradient. With an increase in anisotropy strength,the growth velocity of dendrite increases when the anisotropy strength is lower than the critical value; when the anisotropy strength equals the critical value, the growth velocity drops down by about 4.34%; with an increase again in anisotropy strength (larger than the critical value), the growth velocity reaches the maximum value and then tends to decrease. The equilibrium morphology of the crystal grows into a square-like one in the situation of low dimensionless supercooling degree; with an increase in dimensionless thermal supercooling degree, the equilibrium morphology of the crystal changes from square-like to dendrite, the growth velocity of dendrite increases exponentially and then increases linearly, and the dendrite growth is under control from diffusion to kinetics.

Ni-Cu alloy; anisotropy strength; interface energy; supercooling degree; phase-field method; dendrite growth

TG146

A

1004-0609(2011)07-1656-08

兰州市科技局资助项目(2009-1-9); 兰州理工大学博士基金资助项目(SB01200606)

2010-07-30；

2010-12-18

丁雨田，教授，博士；电话：0931-2757285；E-mail: yuanyang2011@163.com

(编辑 陈卫萍)