相場法模擬強各向異性作用下二元合金枝晶生長

袁訓鋒, 丁雨田

相場法模擬強各向異性作用下二元合金枝晶生長

袁訓鋒, 丁雨田

(蘭州理工大學 甘肅省有色金屬新材料省部共建國家重點實驗室 蘭州 730050)

基于Wheeler模型和Eggleston修正強界面能各向異性的方法，建立耦合溶質場和溫度場的相場模型，模擬強界面能和界面動力學各向異性下Ni-Cu合金的枝晶生長過程。結果表明：兩種各向異性均顯著影響枝晶生長，在強界面動力學各向異性下，固相以類矩形方式沿〈110〉方向生長；在強界面能各向異性及同時存在兩種各向異性下，固相以枝晶方式沿〈100〉方向生長，界面方向不連續，枝晶臂主枝尖端出現棱角。在各向異性強度取值相同情況下，僅有界面能各向異性時，〈100〉方向枝晶尖端溫度梯度大，生長迅速，穩態生長速度比同時存在兩種各向異性的大 32.26%；僅有界面動力學各向異性時，〈100〉方向枝晶尖端溫度梯度小且溶質濃度高，生長緩慢，穩態生長速度比同時存在兩種各向異性的小48.92%。

二元合金；界面能各向異性；界面動力學各向異性；相場法；枝晶生長

過冷熔體中的枝晶生長是凝固科學中的一個基礎問題，同時也是理解其它復雜凝固現象的前提[1]。枝晶的形成過程是一個涉及到熱量、質量和動量傳輸以及界面動力學和毛細作用效應相耦合的自由邊界問題，采用理論分析方法進行研究會遇到巨大的數學困難，而傳統試驗方法又無法揭示枝晶結構形成的有關物理機理。隨著計算科學的發展和凝固理論的完善，采用數值模擬技術研究過冷熔體中的枝晶生長已成為凝固微觀組織研究的熱點。

相場方法在模擬凝固過程中微觀組織形成方面具有巨大的優勢。相場法以Ginzburg- Landau理論為基礎，通過引入序參量場把液固界面由尖銳界面拓展為彌散界面，使得在整個計算區域中所有變量都連續可微，并且可以使用統一的控制方程描述，不必區分相和界面[2]。自 LANGER[3]第一次提出采用擴散界面模型研究凝固現象以來，經過國內外學者不懈努力，相場法模擬凝固微觀組織經歷了從二維到三維[4?6]、從二元合金到多元合金[7]、從單晶粒到多晶粒[8?9]，從自由枝晶到定向凝固[10]、從單相場到多相場[11]、從無流場到包含流場[12?15]的逐步深入的發展歷程，使得模擬結果越來越接近真實凝固過程。在 2001年，EGGLESTON等[16]提出了修正強界面能各向異性的方法，從而解決晶向缺失問題，建立強界面能各向異性下的枝晶生長相場模型。KASAJIMA等[17]和 KIM等[18]采用此相場模型研究了強各向異性強度下過冷純Si枝晶生長特性。在國內，李俊杰等[19]、ZHANG等[20]和CHEN等[21]采用修正的相場模型，對強界面能各向異性下的枝晶生長行為進行了模擬。趙達文等[22?23]采用有限元法求解相場方程研究了弱界面能和界面動力學各向異性下過冷熔體中的枝晶生長行為。到目前為止，很少見到有關強界面能和界面動力學各向異性下枝晶生長的報道，人們對強界面能和界面動力學各向異性作用下的枝晶生長行為仍缺乏了解。

本文作者基于Wheeler模型和Eggleston修正強界面能各向異性的方法，建立耦合溶質場和溫度場的相場模型，模擬強界面能和界面動力學各向異性作用下Ni-Cu合金的枝晶生長過程。研究強界面能各向異性、強界面動力學各向異性及同時存在兩種各向異性下的枝晶生長形貌、溶質和溫度分布及〈100〉方向枝晶尖端的生長行為。

1 相場模型

1.1 相場控制方程

晶體原子按照特定對稱性排列，導致固液界面能、界面動力學具有不同程度的各向異性。在枝晶生長過程中液固界面各向異性不僅決定枝晶生長方向，還在很大程度上影響枝晶生長行為[24]。在液固界面上，通常界面能與界面動力學各向異性同時存在。在相場模型中與界面能相關的是彌散界面厚度 ε(θ)，與界面動力學效應相關的是原子馳豫時間τ(θ).這里選擇界面各向異性具有四重對稱性，即彌散界面厚度ε(θ)=ε( 1 + ε4cos(4θ))、 原 子 馳 豫 時 間 τ(θ)=τ0(1+εkcos(4θ))2，其中，θ為界面法向與X正方向之間的夾角，ε4和εk分別為界面能和界面動力學各向異性強度。

界面能各向異性對平衡形狀的影響可以用Gibbs-Thomson方程描述：

式中：R(θ)為固液界面的曲率半徑，fL和fS分別為液相和固相的自由能密度。在弱界面能各向異性時(ε4＜1/15)，方程的兩邊都為正，枝晶形貌光滑連續；在強界面能各向異性時(ε4＞1/15)，方程左邊為負，導致出現缺失方向，界面變的不連續，凹下去的部分出現耳子。為了模擬強界面能各向異性下的枝晶生長，必須將耳子消除。采用EGGLESTON等[16]的方法修正缺失方向得到的界面厚度，如下：

式中：i=0～3，θm為平衡形貌出現耳子時的角度，可以通過下式求得。

基于 WHEELER等[15]的相場模型和 EGGLESTON等[16]的修正方法建立了強界面能各向異性作用下枝晶生長的修正相場方程。

式中：Φ=0代表固相，Φ=1代表液相，在固/液界面上Φ在0→1之間連續變化；溫度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/cm3K，σ=3.38×10?5J/cm2， κ=0.27 cm2/s， Ds=1.0×10?9cm2/s ， Dl=1.0×10?5cm2/s ， μ=175.3 cm/Ks，ω=2.1×10?4cm，mL=0.821 1，k0=0.86，α=419，m=0.017，ε=0.005，τ0=1.0，ε4=0.08，ΔT=0.1，θm=0.258，x0=40.38%，ΔX=ΔY=0.005，Δt=2.0×10?5。

2.2 數值計算方法

對于傳輸方程(7)采用交替隱式格式求解，控制方程(4)～(6)和(8)同時采用顯示差分格式求解。由于溫度場方程所采用的交替隱式格式具有任意的穩定性，整個計算過程的穩定性受以下條件約束：

式中：Δt為時間步長，ΔX=ΔY為空間步長，m′=max(Dl，m)，k′為考慮方程中非線性項而設計的修正系數，一般取1～2。

2.3 初始條件和邊界條件

開始計算時假設網格數為R的球形初始晶核位于充滿均勻過冷合金熔體的中心(X0, Y0)，即：

在計算區域的邊界上，相場、溫度場、溶質場采用Neumann條件。

3 結果與分析

3.1 強各向異性下的枝晶形貌

圖1所示為Ni-Cu40.38%合金在不同各向異性作用下ΔT為0.1時的枝晶生長形貌。模擬計算網格數為800×800，圖1顯示為模擬中心400×400區域。從圖1(a)可以看出，在強界面能各向異性作用下，固相以枝晶方式沿〈100〉方向生長，枝晶臂主枝細長，根部明顯頸縮，并且由于某些方向的生長消失，界面方向不連續，枝晶臂主枝尖端出現棱角。在強界面動力學各向異性作用下，固相沿〈110〉方向生長，如圖1(b)所示，主枝不發達，整體為類矩形形貌。在強界面能和界面動力學各向異性共同作用下，固相以枝晶方式沿〈100〉方向生長，如圖 1(c)所示，枝晶臂主枝發達，尖端出現棱角，根部無明顯頸縮現象。

3.2 強各向異性下的枝晶溶質分布

圖2所示為圖1中不同各向異性作用下枝晶相場形貌對應的溶質分布。為了方便比較分析，對強界面能各向異性、強界面動力學各向異性及同時存在強界面能和界面動力學各向異性作用下的溶質分布采用統一的標尺，即相同顏色對應同一溶質濃度值。由圖 2可以看出，溶質分布情況與枝晶形貌是相吻合的。從圖2可以看出，由于在凝固過程中存在溶質再分配，先凝固的枝晶主枝中心對應的溶質濃度最低，而液相中溶質的擴散速度小于枝晶生長速度，凝固析出的溶質不能充分擴散到液相中，從而大量溶質富集在枝晶臂前沿，其溶質濃度高。枝晶臂根部存在微觀偏析，各枝晶臂的溶質擴散層分布相同。

圖3所示為圖2中的枝晶縱向主軸方向溶質分布情況。從圖3可以看出，由于Ni-Cu合金的固液界面溶質平衡分配系數小于 1，枝晶生長速度愈大，實際的固液界面溶質分配系數也越大，界面能各向異性為0.08時的枝晶生長速度比界面能和界面動力學各向異性均為0.08時的大，故在僅有界面能各向異性時的固相及固液界面前沿溶質濃度比界面能和界面動力學各向異性均存在時的要高。界面動力學各向異性為0.08時的固相及固液界面前沿溶質濃度最高，這是由于只有界面動力學各向異性時，在〈110〉方向原子的馳豫時間最小，固相沿〈110〉方向生長快，而在縱向主軸〈100〉方向原子馳豫時間長，其生長緩慢，溶質富集嚴重，溶質濃度高?？傊?，〈100〉方向的固相及固液界面前沿溶質濃度在僅有界面動力學各向異性時的最高，僅有界面能各向異性時次之，界面能和界面動力學各向異性均存在時最低。

3.3 強各向異性下的枝晶溫度分布

圖3 不同各向異性下枝晶縱向主軸方向的溶質分布曲線Fig.3 Solute distribution of dendrite along longitudinal principal axis at different anisotropy

圖4 Ni-40.38%Cu合金在不同各向異性下ΔT=0.1時的枝晶溫度分布Fig.4 Temperature distribution of dendrite of Ni-40.38%Cu alloy with different anisotropy at melt undercooling of ΔT=0.1:(a) ε4=0.08; (b) εk=0.08; (c) ε4=εk=0.08

圖4 所示為圖1中不同各向異性作用下枝晶相場形貌對應的溫度分布。為了方便比較分析，對3種情況下的溫度分布采用統一的標尺，即相同顏色對應同一溫度值。從圖4中可以看出，枝晶周圍由等溫線包圍，枝晶臂尖端等溫線緊實，溫度梯度大，其生長速度快；枝晶臂之間的區域等溫線稀疏，溫度梯度小，其生長緩慢。在界面能各向異性為0.08時，枝晶的4個尖端生長迅速，枝晶臂發達，等溫線為類星形；在界面動力學各向異性為0.08時，4個優先生長方向的生長優勢不明顯，各向異性程度較弱，等溫線為類球形；在界面能和界面動力學各向異性均為0.08時，4個優先生長方向上枝晶臂優先生長，等溫線為類矩形。圖5所示為圖4中的枝晶縱向主軸方向溫度分布情況。從圖5中可以看出，曲線中心平坦部分對應已凝固的枝晶固相，溫度較高，溫度梯度??；固相兩側為固液界面及液相，在固液界面上，曲線斜率最大，溫度梯度最大，隨后溫度梯度減小。在僅有界面能各向異性時，〈100〉方向枝晶尖端固液界面溫度梯度最大，枝晶生長最快；在僅有界面動力學各向異性時，〈100〉方向枝晶尖端固液界面溫度梯度最小，枝晶生長最慢；在界面能和界面動力學各向異性均存在時，〈100〉方向枝晶尖端固液界面溫度梯度和枝晶生長速度位于二者之間。

圖5 不同各向異性下枝晶縱向主軸方向溫度分布曲線Fig.5 Temperature distribution of dendrite along longitudinal principal axis at different anisotropy

3.4 強各向異性作用下枝晶的尖端行為

為了定量分析強各向異性對枝晶生長行為的影響，本文作者計算了在僅有界面能各向異性和僅有界面動力學各向異性情況下，〈100〉方向枝晶尖端生長速度，尖端溶質摩爾分數和尖端溫度隨時間的變化關系，并與界面能和界面動力學各向異性均存在情況下進行了對比，結果見圖6。

從圖6可以看出，在凝固初始階段，驅動枝晶生長的過冷度為初始過冷度，在枝晶尖端附近的初始溫度梯度大，結晶潛熱向過冷熔體中擴散速度大，枝晶尖端以較快的速度生長。隨著生長過程的進行，熱量在凝固前沿不斷積累，熱擴散邊界層厚度逐漸增加，溫度梯度不斷減小，使得熱擴散速度降低，枝晶尖端生長速度逐漸降低[22]。同時，由于枝晶尖端生長速度降低，實際的固液界面溶質分配系數減小，導致僅有界面能各向異性及界面能和界面動力學各向異性均存在的〈100〉方向枝晶尖端溶質濃度降低。而僅有界面動力學各向異性時〈100〉方向枝晶尖端溶質濃度不斷增加，這是由于〈100〉方向原子馳豫時間長，該方向生長緩慢，隨著枝晶生長，溶質的富集越來越嚴重。經過一段時間后，溶質和潛熱在界面前沿的釋放與通過擴散的遷移基本達到動態平衡，枝晶尖端生長速度、溶質濃度和溫度逐漸趨于穩定。在僅有界面能各向異性時的〈100〉方向枝晶尖端穩態生長速度和溶質摩爾分數分別比界面能和界面動力學各向異性均存在時的大32.26 %和0.01%，穩態溫度比界面能和界面動力學各向異性均存在時的低13.35%；在僅有界面動力學各向異性時的穩態生長速度比界面能和界面動力學各向異性均存在時的小48.92%，穩態溶質摩爾分數和穩態溫度分別比界面能和界面動力學各向異性均存在時的高0.1%和35.17%。

圖6 枝晶尖端生長速度、尖端溫度和尖端溶質濃度與凝固時間的關系Fig.6 Relationship among tip velocity (a), tip solute concentration (b) and tip temperature (c) and time for dendrite growth

4 結論

1)基于Wheeler模型和Eggleston修正強界面能各向異性的方法，建立耦合溶質場和溫度場的相場模型，再現Ni-Cu合金強界面能和界面動力學各向異性下的枝晶生長過程。

2)在強界面動力學各向異性作用下，固相沿〈110〉方向生長，整體為類矩形形貌；在強界面能各向異性及同時存在強界面能和界面動力學各向異性作用下，固相以枝晶方式沿〈100〉方向生長，界面方向不連續，枝晶臂主枝尖端出現棱角。

3) 在各向異性強度取值相同的情況下，僅有界面動力學各向異性時的〈100〉方向枝晶尖端溶質濃度最高；僅有界面能各向異性時的〈100〉方向枝晶尖端溶質濃度次之；界面能和界面動力學各向異性均存在時的〈100〉方向枝晶尖端溶質濃度最低。

4) 在界面能各向異性強度為 0.08時，〈100〉方向枝晶尖端枝晶溫度梯度大，生長迅速，穩態生長速度比界面能和界面動力學各向異性強度均為 0.08的大32.26%；在界面動力學各向異性強度為0.08時，〈100〉方向枝晶尖端溫度梯度小，生長緩慢，穩態生長速度比界面能和界面動力學各向異性強度均為 0.08的小48.92%。

REFERENCES

[1] KURZ W, FISHER D J. Fundamentals of solidification [M].Switzerland: Trans Tech, 1998: 45.

[2] CAGINALP G, FIFE P. Higher-order phase field models and detailed anisotropy[J]. Phys Rev B, 1986, 34(7): 4940?4943.

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

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

[5] 趙紅兆, 荊 濤, 柳百成. 鋁合金三維枝晶生長相場模擬[J].金屬學報, 2005, 41(5): 491?495.ZHAO Hong-zhao, JING Tao, LIU Bai-cheng. Simulation of 3D dendritic growth of aluminum alloy by phase field model[J].Acta Metall Sin, 2005, 41(5): 491?495.

[6] 朱昌盛, 馮 力, 王智平, 肖榮振. 三維枝晶生長的相場法數值模擬研究[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.

[7] ODE M, LEE J S, KIM W T. Phase-field model for solidification of ternary alloys[J]. ISIJ International, 2000, 40(9): 870?876.

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

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

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

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

[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 (Engl Lett), 2008, 21(6): 444?450.

[14] 龍文元, 呂冬蘭, 夏 春, 潘美滿, 蔡啟舟, 陳立亮. 強迫對流影響二元合金非等溫凝固枝晶生長的相場法模擬[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 non-isothermal solidification dendrite growth of binary alloy under the force flow[J]. Acta Phys Sin, 2009, 58(11):7802?7808.

[15] WHEELER A A, MURRAY B T, SCHAEFER R J. Computation of dendrites using a phase field model [J]. Physica D, 1993,66(1/2): 243?262.

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

[17] KASAJIMA H, NAGANO E, SUZUKI T, KIM S G, KIM W T.Phase-field modeling for facet dendrite growth of silicon[J].Science and Technology of Advanced Materials, 2003, 4(6):553?557.

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

[19] 李俊杰, 王錦程, 楊根倉. 相場法模擬界面能各向異性對枝晶生長行為的影響[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.

[20] 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(Z2): 307?313.

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

[22] 趙達文, 李金富. 相場模型模擬液固界面各向異性作用下自由枝晶生長[J]. 物理學報, 2009, 58(10): 7094?7100.ZHAO Da-wen, LI Jin-fu. Phase-field modeling of the effect of liquid-solid interface anisotropies on free dendritic growth[J].Acta Phys Sin, 2009, 58(10): 7094?7100.

[23] 趙達文, 李金富. 相場法模擬動力學各向異性對過冷熔體中晶體生長的影響[J]. 金屬學報, 2009, 45(10): 1237?1241.ZHAO Da-wen, LI Jin-fu. Phase-field simulation of the effect of kinetic anisotropy on crystal growth in undercooled melts[J].Acta Metall Sin, 2009, 45(10): 1237?1241.

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

Phase-field simulation of dendrite growth for binary alloy with strong anisotropy

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

Based on the Wheeler model and the Eggleston regularization technique of strong anisotropy of interface energy, the phase-field model was built by coupling with the concentration field and temperature field. The dendrite growth process of Ni-Cu alloy with strong surface energy and kinetic anisotropy were simulated. The results show that the dendrite growth depends on the two kinds of anisotropies. Under the strong surface kinetic anisotropy condition, the melt solidifies grow along the 〈110〉 orientation and the crystals grow into a square-like. Under the strong surface energy anisotropy or having two kinds of anisotropies condition, the melt solidifies in a dendrite pattern grow along the 〈100〉orientation and the variation of interface orientation discontinuity can lead to the corners form on the tip of dendrite. In the case of anisotropy strength with the same values, under the strong surface energy anisotropy condition, the thermal gradient along the 〈100〉 orientation is large, and makes the dendrite growth become fast, the tip velocity at steady state increases by about 32.26% compared with the case that having two kinds of anisotropies. Under the strong surface kinetic anisotropy condition, the thermal gradient along the 〈100〉 orientation is small and the concentration of solute is large, and makes the dendrite growth become slow, the tip velocity at steady state decreases by about 48.92% compared with the case that having two kinds of anisotropies.

binary alloy; surface energy anisotropy; surface kinetic anisotropy; phase-field; dendrite growth

TG146

A

1004-0609(2011)09-2216-07

蘭州市科技局資助項目(2009-1-9)；蘭州理工大學博士基金資助項目(SB01200606)

2010-08-04；

2011-01-15

丁雨田，教授，博士；電話：0931-2757285；E-mail: dingyt @lut.cn

(編輯 何學鋒)