強界面能各向異性下二元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

(編輯 陳衛萍)