鎂合金微觀偏析相場法定量數(shù)值模擬

堯軍平 ，李翔光，龍文元，張 磊，王法軍，李 文

(1. 南昌航空大學(xué) 輕合金加工科學(xué)與技術(shù)國防重點學(xué)科實驗室，南昌 330063；2. 俄亥俄州立大學(xué) 工程學(xué)院，哥倫布，俄亥俄州 43210，美國)

鎂合金具有比強度和比剛度高、密度低、阻尼性和切削加工性好等優(yōu)點，是目前實際應(yīng)用最輕的金屬結(jié)構(gòu)材料，用鎂合金生產(chǎn)的鑄件具有優(yōu)良的減振和降噪性能、耐磨性能、低慣性能、電磁屏蔽性、易回收特性等，被稱為21世紀的“綠色”工程材料，在電子、汽車和航空工業(yè)領(lǐng)域有廣泛的應(yīng)用[1]。由于合金元素的加入，鎂合金凝固過程中溶質(zhì)元素的重新分配會造成合金元素在晶粒內(nèi)部和晶界的分布不均勻，從而產(chǎn)生微觀偏析[2]。微觀偏析可以引起非平衡第二相、氣孔和裂紋等的形成，這些都嚴重損害鎂合金鑄件的力學(xué)性能。因此，對鎂合金凝固過程中的微觀偏析進行定量分析預(yù)測就顯得十分重要。以前預(yù)測合金的微觀偏析程度常使用實驗方法和分析方法，但是，采用這些方法很難觀察和測定小枝晶區(qū)域內(nèi)的溶質(zhì)偏析。此外，分析模型對實際問題做了很多假設(shè)，因此，該方法不能特別準確地對實際現(xiàn)象進行描述。相場法作為目前國內(nèi)外很受重視的模擬凝固微觀組織的有力工具，通過相場與溶質(zhì)場以及溫度場的耦合，可以對鎂合金凝固過程中的固相和液相的溶質(zhì)偏析進行直接模擬預(yù)測，并據(jù)此調(diào)整優(yōu)化生產(chǎn)工藝，從而得到力學(xué)性能更優(yōu)良的鑄件，所以，開展鎂合金微觀偏析相場法定量模擬分析具有很強的現(xiàn)實與學(xué)術(shù)意義。

目前，KUNDIN等[3]和袁訓(xùn)鋒等[4]所進行的微觀組織模擬絕大部分是針對具有FCC晶胞的金屬 (例如Ni、Al和Cu)。而鎂合金是典型的HCP晶體結(jié)構(gòu)，具有不同于FCC結(jié)構(gòu)晶體的各向異性特征。因此，要對鎂合金凝固過程中的微觀偏析進行模擬分析就必須建立適合HCP晶體結(jié)構(gòu)金屬的數(shù)值模型。

在鎂合金微觀組織相場法模擬方面，EIKEN等[5-7]對鎂合金進行了二維、三維以及定向凝固微觀組織的模擬。WANG等[8-9]模擬了鎂合金三維等軸晶的生長。MONTIEL等[10]對鎂合金焊接條件下的微觀組織進行了模擬研究?？娂颐鞯萚11]研究了模擬參數(shù)對鎂合金枝晶生長過程的影響。袁訓(xùn)鋒等[12]基于WHEELER等[13]提出的模型研究了強迫對流對鎂合金凝固過程中單晶粒和多晶粒生長行為的影響。但以上工作主要集中于研究鎂合金的枝晶生長形貌，采用相場法定量研究鎂合金微觀偏析方面的工作尚未見報道。本文作者基于KIM等[14-16]在WMB模型的基礎(chǔ)上改進獲得的KKS相場模型，發(fā)展一個耦合相場、溶質(zhì)場及溫度場且適合于HCP晶系的相場模型，以AZ91D鎂合金為例，模擬研究鎂合金凝固時的溶質(zhì)分布情況，定量分析過冷度對枝晶尖端前沿最高溶質(zhì)濃度、枝晶干軸對稱中心溶質(zhì)濃度、枝晶尖端溶質(zhì)擴散層厚度以及溶質(zhì)偏析比的影響。

1 相場模型

1.1 相場控制方程

本文作者采有KIM等[14-16]提出的KKS相場模型，其相場控制方程表示為

式中：相場變量φ=0代表液態(tài)，φ=1代表固態(tài)；c為溶質(zhì)濃度，上標e表示平衡狀態(tài)，下標S和L分表表示固相和液相；R是氣體常數(shù)；T是溫度；Vm是摩爾體積；h′(φ)為勢函數(shù)的導(dǎo)數(shù)，h′(φ)=30φ2(φ2-2φ+1)；g′(φ)為剩余自由能函數(shù)的導(dǎo)數(shù)，g′(φ)=2φ(2φ2-3φ+1)；W為相場參數(shù)，Mφ是與界面動力學(xué)有關(guān)的相場參數(shù)，ε為與界面能有關(guān)的參數(shù)。

為了考慮界面能各向異性，可將梯度項系數(shù)ε表示為

式中：k為各向異性的模數(shù)，表示系統(tǒng)的對稱次數(shù)，為適應(yīng)六重對稱性鎂合金的求解，其值取6；v為各向異性強度；θ表示界面與優(yōu)先生長方向之間的夾角，θ=arctan(φy/φx)。引入界面能各向異性后，方程(1)變成：

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

在合金的模擬中，還需要同時求解溶質(zhì)擴散方程，其控制方程可表示為

式中：D(φ)為體系的溶質(zhì)擴散系數(shù)，D(φ)=DL+h(φ)(DS-DL)，DS和DL分別代表固、液相中的溶質(zhì)擴散系數(shù)。

1.3 溫度場控制方程

在研究中，同時耦合計算了溫度場。由凝固過程的傳熱主導(dǎo)方程式可以得到溫度場的控制方程為

式中：T為溫度；L為潛熱；cp為比熱容；D=λ/(ρcp)；λ為導(dǎo)熱系數(shù)。

1.4 擾動

擾動可引發(fā)枝晶側(cè)向分枝的生長，為了模擬實際凝固過程中界面處的波動和二次或更高次枝晶臂的生長，本文作者在溶質(zhì)擴散方程中加入隨機擾動，即

式中：χ為取值在-1~+1之間的隨機數(shù)；ω˙為與時間有關(guān)的擾動強度因子。

2 數(shù)值求解

2.1 合金熱物性參數(shù)

模擬所用AZ91D鎂合金物性參數(shù)如表1[14]所列。

表1 AZ91D鎂合金物性參數(shù)[17]Table 1 Physical properties of AZ91D magnesium alloy[17]

2.2 數(shù)值求解方法

流場控制方程采用Simple算法求解，采用顯式差分法求解相場和溶質(zhì)場控制方程，采用ADI方法求解溫度場控制方程，時間步長Δt=Δx2/(5Dl)，模擬區(qū)域網(wǎng)格劃分為1 000×1 000，相場、濃度場的網(wǎng)格尺寸為1.0×10-8m。為了減少計算量，本文作者在數(shù)值求解時采用了雙重網(wǎng)格技術(shù)，溫度場的計算網(wǎng)格數(shù)為100×100，網(wǎng)格尺寸為1.0×10-7m。計算程序是基于VC++6.0平臺在Windows xp系統(tǒng)下開發(fā)，采用tecplot 360實現(xiàn)計算結(jié)果的可視化。

2.3 初始條件和邊界條件

假設(shè)初始晶核半徑為R(R=10Δx)，則

式中：x、y分別代表橫坐標、縱坐標；(a，b)代表晶核中心坐標。在計算區(qū)域的邊界上，溫度場采用第一類溫度邊界條件，相場、溶質(zhì)場均采用Zero-Neumann邊界條件。

3 模擬結(jié)果與分析

3.1 枝晶形貌與溶質(zhì)分布

圖1所示為對AZ91D鎂合金過冷熔體的凝固過程進行模擬而獲得的溶質(zhì)分布圖。由圖1可見，在{0001}基面上，整個枝晶沿著〈12 10〉晶向以對稱形式生長，一次枝晶主干間互成60°夾角，一次枝晶干與該主干上的二次枝晶臂間也約成60°夾角。貧Al區(qū)主要集中于一次枝晶干和二次枝晶臂的軸對稱中心位置，這主要是由于凝固過程枝晶尖端曲率效應(yīng)引起過冷，使固相線向下移動，而固相中溶質(zhì)的擴散速度又遠遠落后于枝晶的生長速度。在液固界面前沿出現(xiàn)了溶質(zhì)Al的富集，其中以一次枝晶干根部液固界面前沿的溶質(zhì)富集程度最高，這主要是由于鎂合金的凝固屬于非平衡結(jié)晶過程，合金在凝固過程中需要對溶質(zhì)進行再分配，這導(dǎo)致固相中溶質(zhì)的濃度低于初始濃度，而液相中溶質(zhì)的擴散速度也要小于枝晶生長速度，凝固析出的溶質(zhì)不能及時充分地擴散到液相中，從而在枝晶前沿富集。一次枝晶干根部處的溶質(zhì)擴散通道由于二次枝晶臂的阻擋而變得狹長，嚴重阻礙了該處溶質(zhì)的擴散，從而使該處溶質(zhì)富集程度最高。另外，由于枝晶尖端生長速度最大，溶質(zhì)來不及擴散，溶質(zhì)擴散層最薄，因此，該處的濃度梯度最大。

圖1 單個枝晶溶質(zhì)分布模擬圖Fig. 1 Simulated solute distribution scheme of equiaxed dendrite (ΔT=103 K, t=6 000Δt)

圖2所示為砂型鑄件AZ91D在光學(xué)顯微鏡偏振光下的枝晶結(jié)構(gòu)金相照片[18]。對比圖1和2的可以發(fā)現(xiàn)，模擬結(jié)果與文獻[18]中報道的實驗結(jié)果在枝晶形貌上十分吻合，這也在一方面驗證了文中模型的可行性與正確性。

圖2 砂型AZ91D 鑄態(tài)的微觀組織[18]Fig. 2 As-casting microstructure of sand mould AZ91D(under polarized light)[18]

3.2 過冷度對枝晶尖端前沿最高濃度的影響

圖3 不同過冷度下一次枝晶干軸對稱中心沿〈1210〉晶向的濃度分布Fig. 3 Concentration distribution under different supercooling(t=6 000Δt)

圖3顯示的為不同過冷度下凝固時間t=6 000Δt時一次枝晶干軸對稱中心沿〈1210〉晶向的濃度分布情況。圖中、、和分別對應(yīng)過冷度為108、103、98和93 K時枝晶尖端前沿固液界面處液相一側(cè)的最高溶質(zhì)濃度，cS*1、cS*2、cS*3以及cS*4則分別對應(yīng)過冷度為108、103、98和93 K時沿〈1210〉晶向一次枝晶干軸對稱中心溶質(zhì)Al的濃度，而δ1、δ2、δ3和δ4則分別對應(yīng)過冷度為108、103、98和93 K時枝晶干尖端溶質(zhì)擴散層的厚度。從圖中可以看出，過冷度 ΔT越大，枝晶干尖端前沿固液界面處液相一側(cè)的最高溶質(zhì)濃度越高，即c*L1＞c*L2＞c*L3＞cL*4。這主要是由于鎂合金凝固過程中通過溶質(zhì)再分配向液相中析出溶質(zhì)Al，而溶質(zhì)Al在液相中的擴散速度小于枝晶生長速度，析出的Al來不及充分擴散，而富集在液固界面前沿，而ΔT越大，枝晶的生長速度越快(圖4)，溶質(zhì)Al擴散得也越不充分，則在液固界面前沿的富集也越嚴重。

圖4 不同過冷度下的枝晶尖端生長速度與凝固時間的關(guān)系Fig. 4 Relationship between tip growth velocity and solidification time of dendrite at different supercoolings

3.3 過冷度對枝晶干軸對稱中心溶質(zhì)濃度的影響

同時，從圖3中還可以看出，過冷度ΔT越大，一次枝晶干軸對稱中心溶質(zhì)Al的濃度也越高，即＞＞＞。對于這一現(xiàn)象，可做如下解釋：1) 過冷度ΔT增大，凝固點在相圖上的液相線與固相線右移，凝固時平衡濃度增大，從而導(dǎo)致凝固過程固液界面溶質(zhì)再分配時固相中的濃度增大，過冷度 ΔT越大，固相中溶質(zhì)濃度越高；2) 雖然合金的凝固是非平衡結(jié)晶過程，但液固界面處的局部平衡狀態(tài)是一直存在的，不管過冷度是大還是小，液固界面處固相側(cè)溶質(zhì)濃度和液相側(cè)溶質(zhì)濃度的比值都滿足溶質(zhì)平衡擴散系數(shù)k0，即k0=cS*/c*L，也就是說=，結(jié)合上面的＞＞＞，即可得＞＞＞；3) BOWER等[19]通過假設(shè)液相均勻混合固相有限擴散，推導(dǎo)出了液固界面處固相側(cè)溶質(zhì)濃度分布情況，即BF模型，其具體表達式為[20]

式中：為固相溶質(zhì)濃度；為液相溶質(zhì)濃度；fS為固相分數(shù)；α為固相反擴散因子，其值決定固相中溶質(zhì)擴散狀況；DS為固相溶質(zhì)擴散系數(shù)；tf為局部凝固時間，λ為枝晶間距的一半；c0為初始溶質(zhì)濃度，k0為溶質(zhì)平衡分配系數(shù)。

由式(9)和(10)可知，局部凝固時間tf與α成正比，α參數(shù)決定固相中溶質(zhì)的擴散狀況，也就是說，局部時間tf越長，擴散時間越長，固相中溶質(zhì)元素擴散越充分，溶質(zhì)濃度越低。當過冷度 ΔT增大時，枝晶的生長速率也增大(見圖4)，這使得局部的凝固時間tf減小，溶質(zhì)擴散時間縮短，從而使得固相中溶質(zhì)元素擴散越不充分，溶質(zhì)濃度也就越高，從而有＞＞＞。

3.4 過冷度對枝晶尖端溶質(zhì)擴散層厚度的影響

由圖3可以看出，過冷度會對枝晶尖端前沿的溶質(zhì)擴散層厚度δ產(chǎn)生影響。為更直觀地顯示出各枝晶尖端溶質(zhì)擴散層厚度的關(guān)系，截取相同凝固時間(t=6 000Δt)不同過冷度下的枝晶尖端并進行同等倍數(shù)放大后再進行比較，如圖5所示。由圖5可見，過冷度ΔT越大，枝晶尖端前沿溶質(zhì)擴散層越薄，即δ1＜δ2＜δ3＜δ4。這是由于鎂合金等軸晶生長時，從固相中析出的溶質(zhì)Al需通過溶質(zhì)擴散層釋放到液相中去，溶質(zhì)擴散的長度λs=DL/v[21]。由于枝晶的生長速度v隨過冷度ΔT的減小而減小(見圖4)，故熔體中溶質(zhì)擴散的長度λs隨過冷度的減小而增大，這導(dǎo)致包圍枝晶尖端的溶質(zhì)擴散層厚度δ更大，從而有δ1＜δ2＜δ3＜δ4。

圖5 不同過冷度下枝晶尖端溶質(zhì)Al擴散層的厚度Fig. 5 Thickness of solute Al diffusion layer of dendrite tip under different supercoolings and t=6 000Δt: (a) ΔT=108 K; (b)ΔT=103 K; (c) ΔT=98 K; (d) ΔT=93 K

3.5 過冷度對溶質(zhì)偏析比的影響

在通常的鑄件生產(chǎn)中常采用偏析比SR來表示微觀偏析的大小，偏析比在數(shù)值上等于枝晶間最大溶質(zhì)濃度與枝晶干最小溶質(zhì)濃度的比值。為更直觀地顯示過冷度對微觀偏析程度的影響，計算了不同過冷度下溶質(zhì)的偏析比曲線(見圖6)。從圖6中可明顯看出，溶質(zhì)的偏析比SR1＞SR2＞SR3＞SR4，即過冷度ΔT越大，微觀偏析越嚴重。這是由于鎂合金凝固過程中會通過溶質(zhì)再分配向液相中析出溶質(zhì)Al，而溶質(zhì)Al在液相中的擴散速度小于枝晶生長速度，析出的Al來不及充分擴散從而富集在液固界面前沿。過冷度 ΔT越大，枝晶生長速度越快(見圖4)，單位時間內(nèi)析出到液相中的溶質(zhì)就越多，二次枝晶也越發(fā)達，二次枝晶臂對富集溶質(zhì)擴散的阻礙作用也越大，富集的溶質(zhì)越難擴散，從而使得最高溶質(zhì)濃度越大(見圖7)。而固相中溶質(zhì)擴散系數(shù)比液相中的要小得多，過冷度 ΔT對枝晶干最低溶質(zhì)濃度的影響要小得多，因此，這最終導(dǎo)致溶質(zhì)的偏析比越高，微觀偏析越嚴重。

由以上分析可知，模擬結(jié)果與基于凝固理論分析得到的結(jié)論十分吻合，這也從另一方面驗證了模型的可行性與正確性。

微觀偏析的存在會降低合金的塑性、韌性和拉伸性能，因此，必須對微觀偏析予以消除或者減輕。由模擬結(jié)果可知，減小過冷度可以降低合金微觀偏析的程度，但這會使枝晶間距增大，材料力學(xué)性能下降，均勻化退火時間延長。合理的方法應(yīng)該是降低初始溫度，使過冷熔體快速凝固得到細小枝晶，然后再進行長時間的均勻化退火。

圖6 不同過冷度下元素Al的偏析比曲線Fig. 6 Curves of segregation ratio of Al under different supercoolings

圖7 不同過冷度下的溶質(zhì)分布Fig. 7 Solute distribution under different supercoolings: (a) ΔT=108 K; (b) ΔT=103 K; (c) ΔT=98 K; (d) ΔT=93 K; t=6 000Δt

4 結(jié)論

1) 鎂合金貧Al區(qū)主要集中于一次枝晶干和二次枝晶臂的中心，溶質(zhì)富集于液固界面前沿，枝晶尖端濃度梯度最大。

4) 模擬結(jié)果和實驗結(jié)果的對比以及基于凝固理論的分析驗證了模型的可行性與正確性。

5) 消除或者減輕微觀偏析的合理方法是降低初始溫度，然后再進行長時間的均勻化退火。

[1] SCHUMANN S, FRIEDRICH H. Current and future use of magnesium in the automotive industry[J]. Mater Sci Forum,2003, 419/422(1): 51-56.

[2] ESKIN D G. Physical metallurgy of direct chill casting of aluminum alloys[M]. London: CRC Press, 2008: 12.

[3] KUNDIN J, SIQUIERI R, EMMERICH H. A quantitative multi-phase-field modeling of the microstructure evolution in a peritectic Al-Ni alloy[J]. Physica D, 2013, 243(1): 116-127.

[4] 袁訓(xùn)鋒, 丁雨田. 強界面能各向異性下二元Ni-Cu合金枝晶生長過程的相場法模擬[J]. 中國有色金屬學(xué)報, 2011, 21(7):1656-1663.YUAN Xun-feng, DING Yu-tian. Phase-field simulation of dendrite growth process for binary Ni-Cu alloy with anisotropy of strong interface energy[J]. The Chinese Journal of Nonferrous Metals, 2011, 21(7): 1656-1663.

[5] EIKEN J, B?TTGER B, STEINBACH I. Simulations of microstructure evolution during solidification of magnesium-based alloys[J]. Trans Indian Inst Met, 2007, 60(2/3):179-184.

[6] EIKEN J. Dendritic growth texture evolution in magnesium-based alloys investigated by phase-field simulation[J]. International Journal of Cast Metals Research,2009, 22(1/4): 86-89.

[7] B?TTGER B, EIKEN J, STEINBACH I. Phase field simulation of equiaxed solidification in technical alloys[J]. Acta Mater,2006, 54: 2697-2704.

[8] WANG M Y, WILLIAMS J J, JIANG L, CARLO F DE, JING T,CHAWLA N. Dendritic morphology of α-Mg during the solidification of Mg-based alloys: 3D experimental characterization by X-ray synchrotron tomography and phase-field simulations[J]. Scripta Materialia, 2011, 65:855-858.

[9] WANG M Y, JING T, LIU B C. Phase-field simulations of dendrite morphologies and selected evolution of primary α-Mg phases during the solidification of Mg-rich Mg-Al-based alloys[J]. Scripta Materialia, 2009, 61: 777-780.

[10] MONTIEL D, LIU L, XIAO L, ZHOU Y, PRIVATAS N.Microstructure analysis of AZ31 magnesium alloy welds using phase-field models[J]. Acta Mater, 2012, 60: 5925-5932.

[11] 繆家明, 荊 濤, 柳百成. 鎂合金枝晶形貌的相場方法模擬[J]. 金屬學(xué)報, 2008, 44(4): 483-488.MIAO Jia-ming, JING Tao, LIU Bai-chang. Numerical simulation of dendritic morphology of magnesium alloys using phase field method[J]. Acta Metallurgica Sinica, 2008, 44(4):483-488.

[12] 袁訓(xùn)鋒, 丁雨田, 郭廷彪, 胡 勇. 強制對流作用下鎂合金枝晶生長的相場法數(shù)值模擬[J]. 中國有色金屬學(xué)報, 2010, 20(8):1474-1480.YUAN Xun-feng, DING Yu-tian, GUO Tin-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.

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

[14] SUZUKIA T, ODE M, KIM S G, KIM W T. Phase-field model of dendritic growth[J]. J Cryst Growth, 2002, 26(237/239):125-131.

[15] KIM S G, KIM W T, SUZUKIA T. Phase-field model for binary alloys[J]. Phys Rew E, 1999, 60(6): 7186-7197.

[16] KIM S G, KIM W T, SUZUKIA T. Interfacial compositions of solid and liquid in a phase-field model with finite interface thickness for isothermal solidification in binary alloys[J]. Phys Rew E, 1998, 58(3): 3316-332.

[17] 付振南, 許慶彥, 熊守美. 基于概率捕獲模型的元胞自動機方法模擬鎂合金枝晶生長過程[J]. 中國有色金屬學(xué)報, 2007,17(10): 1567-1573.FU Zhen-nan, XU Qing-yan, XIONG Shou-mei. Numerical simulation on dendrite growth process of Mg alloy using cellular automaton method based on probability capturing model[J]. The Chinese Journal of Nonferrous Metals, 2007, 17(10):1567-1573.

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

[19] BOWER T F, BRODY H D, FLEMINGS M C. Measurements of solute redistribution in dendritic solidification[J]. Transaction of the Metallurgical Society of AIME, 1966, 236: 624-633.

[20] CHEN Q, SUNDMAN B. Computation of partial equilibrium solidification with complete interstitial and negligible substitutional solute back diffusion[J]. Materials Transactions,2002, 43(3): 551-559.

[21] 朱昌盛, 王智平, 荊 濤, 肖榮振. 二元合金微觀偏析的相場法數(shù)值模擬[J]. 物理學(xué)報, 2006, 55(3): 1502-1507.ZHU Chang-sheng, WANG Zhi-ping, JING Tao, XIAO Rong-zhen. Numerical simulation of solute segregation patterns for a binary alloy using phase-field approach[J]. Acta Physica Sinica, 2006, 55(3): 1502-1507.