镁合金微观偏析相场法定量数值模拟

尧军平 ，李翔光，龙文元，张 磊，王法军，李 文

(1. 南昌航空大学 轻合金加工科学与技术国防重点学科实验室，南昌 330063；2. 俄亥俄州立大学 工程学院，哥伦布，俄亥俄州 43210，美国)

镁合金具有比强度和比刚度高、密度低、阻尼性和切削加工性好等优点，是目前实际应用最轻的金属结构材料，用镁合金生产的铸件具有优良的减振和降噪性能、耐磨性能、低惯性能、电磁屏蔽性、易回收特性等，被称为21世纪的“绿色”工程材料，在电子、汽车和航空工业领域有广泛的应用[1]。由于合金元素的加入，镁合金凝固过程中溶质元素的重新分配会造成合金元素在晶粒内部和晶界的分布不均匀，从而产生微观偏析[2]。微观偏析可以引起非平衡第二相、气孔和裂纹等的形成，这些都严重损害镁合金铸件的力学性能。因此，对镁合金凝固过程中的微观偏析进行定量分析预测就显得十分重要。以前预测合金的微观偏析程度常使用实验方法和分析方法，但是，采用这些方法很难观察和测定小枝晶区域内的溶质偏析。此外，分析模型对实际问题做了很多假设，因此，该方法不能特别准确地对实际现象进行描述。相场法作为目前国内外很受重视的模拟凝固微观组织的有力工具，通过相场与溶质场以及温度场的耦合，可以对镁合金凝固过程中的固相和液相的溶质偏析进行直接模拟预测，并据此调整优化生产工艺，从而得到力学性能更优良的铸件，所以，开展镁合金微观偏析相场法定量模拟分析具有很强的现实与学术意义。

目前，KUNDIN等[3]和袁训锋等[4]所进行的微观组织模拟绝大部分是针对具有FCC晶胞的金属 (例如Ni、Al和Cu)。而镁合金是典型的HCP晶体结构，具有不同于FCC结构晶体的各向异性特征。因此，要对镁合金凝固过程中的微观偏析进行模拟分析就必须建立适合HCP晶体结构金属的数值模型。

在镁合金微观组织相场法模拟方面，EIKEN等[5-7]对镁合金进行了二维、三维以及定向凝固微观组织的模拟。WANG等[8-9]模拟了镁合金三维等轴晶的生长。MONTIEL等[10]对镁合金焊接条件下的微观组织进行了模拟研究。缪家明等[11]研究了模拟参数对镁合金枝晶生长过程的影响。袁训锋等[12]基于WHEELER等[13]提出的模型研究了强迫对流对镁合金凝固过程中单晶粒和多晶粒生长行为的影响。但以上工作主要集中于研究镁合金的枝晶生长形貌，采用相场法定量研究镁合金微观偏析方面的工作尚未见报道。本文作者基于KIM等[14-16]在WMB模型的基础上改进获得的KKS相场模型，发展一个耦合相场、溶质场及温度场且适合于HCP晶系的相场模型，以AZ91D镁合金为例，模拟研究镁合金凝固时的溶质分布情况，定量分析过冷度对枝晶尖端前沿最高溶质浓度、枝晶干轴对称中心溶质浓度、枝晶尖端溶质扩散层厚度以及溶质偏析比的影响。

1 相场模型

1.1 相场控制方程

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

式中：相场变量φ=0代表液态，φ=1代表固态；c为溶质浓度，上标e表示平衡状态，下标S和L分表表示固相和液相；R是气体常数；T是温度；Vm是摩尔体积；h′(φ)为势函数的导数，h′(φ)=30φ2(φ2-2φ+1)；g′(φ)为剩余自由能函数的导数，g′(φ)=2φ(2φ2-3φ+1)；W为相场参数，Mφ是与界面动力学有关的相场参数，ε为与界面能有关的参数。

为了考虑界面能各向异性，可将梯度项系数ε表示为

式中：k为各向异性的模数，表示系统的对称次数，为适应六重对称性镁合金的求解，其值取6；v为各向异性强度；θ表示界面与优先生长方向之间的夹角，θ=arctan(φy/φx)。引入界面能各向异性后，方程(1)变成：

1.2溶质场控制方程

在合金的模拟中，还需要同时求解溶质扩散方程，其控制方程可表示为

式中：D(φ)为体系的溶质扩散系数，D(φ)=DL+h(φ)(DS-DL)，DS和DL分别代表固、液相中的溶质扩散系数。

1.3 温度场控制方程

在研究中，同时耦合计算了温度场。由凝固过程的传热主导方程式可以得到温度场的控制方程为

式中：T为温度；L为潜热；cp为比热容；D=λ/(ρcp)；λ为导热系数。

1.4 扰动

扰动可引发枝晶侧向分枝的生长，为了模拟实际凝固过程中界面处的波动和二次或更高次枝晶臂的生长，本文作者在溶质扩散方程中加入随机扰动，即

式中：χ为取值在-1~+1之间的随机数；ω˙为与时间有关的扰动强度因子。

2 数值求解

2.1 合金热物性参数

模拟所用AZ91D镁合金物性参数如表1[14]所列。

表1 AZ91D镁合金物性参数[17]Table 1 Physical properties of AZ91D magnesium alloy[17]

2.2 数值求解方法

流场控制方程采用Simple算法求解，采用显式差分法求解相场和溶质场控制方程，采用ADI方法求解温度场控制方程，时间步长Δt=Δx2/(5Dl)，模拟区域网格划分为1 000×1 000，相场、浓度场的网格尺寸为1.0×10-8m。为了减少计算量，本文作者在数值求解时采用了双重网格技术，温度场的计算网格数为100×100，网格尺寸为1.0×10-7m。计算程序是基于VC++6.0平台在Windows xp系统下开发，采用tecplot 360实现计算结果的可视化。

2.3 初始条件和边界条件

假设初始晶核半径为R(R=10Δx)，则

式中：x、y分别代表横坐标、纵坐标；(a，b)代表晶核中心坐标。在计算区域的边界上，温度场采用第一类温度边界条件，相场、溶质场均采用Zero-Neumann边界条件。

3 模拟结果与分析

3.1 枝晶形貌与溶质分布

图1所示为对AZ91D镁合金过冷熔体的凝固过程进行模拟而获得的溶质分布图。由图1可见，在{0001}基面上，整个枝晶沿着〈12 10〉晶向以对称形式生长，一次枝晶主干间互成60°夹角，一次枝晶干与该主干上的二次枝晶臂间也约成60°夹角。贫Al区主要集中于一次枝晶干和二次枝晶臂的轴对称中心位置，这主要是由于凝固过程枝晶尖端曲率效应引起过冷，使固相线向下移动，而固相中溶质的扩散速度又远远落后于枝晶的生长速度。在液固界面前沿出现了溶质Al的富集，其中以一次枝晶干根部液固界面前沿的溶质富集程度最高，这主要是由于镁合金的凝固属于非平衡结晶过程，合金在凝固过程中需要对溶质进行再分配，这导致固相中溶质的浓度低于初始浓度，而液相中溶质的扩散速度也要小于枝晶生长速度，凝固析出的溶质不能及时充分地扩散到液相中，从而在枝晶前沿富集。一次枝晶干根部处的溶质扩散通道由于二次枝晶臂的阻挡而变得狭长，严重阻碍了该处溶质的扩散，从而使该处溶质富集程度最高。另外，由于枝晶尖端生长速度最大，溶质来不及扩散，溶质扩散层最薄，因此，该处的浓度梯度最大。

图1 单个枝晶溶质分布模拟图Fig. 1 Simulated solute distribution scheme of equiaxed dendrite (ΔT=103 K, t=6 000Δt)

图2所示为砂型铸件AZ91D在光学显微镜偏振光下的枝晶结构金相照片[18]。对比图1和2的可以发现，模拟结果与文献[18]中报道的实验结果在枝晶形貌上十分吻合，这也在一方面验证了文中模型的可行性与正确性。

图2 砂型AZ91D 铸态的微观组织[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〉晶向的浓度分布情况。图中、、和分别对应过冷度为108、103、98和93 K时枝晶尖端前沿固液界面处液相一侧的最高溶质浓度，cS*1、cS*2、cS*3以及cS*4则分别对应过冷度为108、103、98和93 K时沿〈1210〉晶向一次枝晶干轴对称中心溶质Al的浓度，而δ1、δ2、δ3和δ4则分别对应过冷度为108、103、98和93 K时枝晶干尖端溶质扩散层的厚度。从图中可以看出，过冷度 ΔT越大，枝晶干尖端前沿固液界面处液相一侧的最高溶质浓度越高，即c*L1＞c*L2＞c*L3＞cL*4。这主要是由于镁合金凝固过程中通过溶质再分配向液相中析出溶质Al，而溶质Al在液相中的扩散速度小于枝晶生长速度，析出的Al来不及充分扩散，而富集在液固界面前沿，而ΔT越大，枝晶的生长速度越快(图4)，溶质Al扩散得也越不充分，则在液固界面前沿的富集也越严重。

图4 不同过冷度下的枝晶尖端生长速度与凝固时间的关系Fig. 4 Relationship between tip growth velocity and solidification time of dendrite at different supercoolings

3.3 过冷度对枝晶干轴对称中心溶质浓度的影响

同时，从图3中还可以看出，过冷度ΔT越大，一次枝晶干轴对称中心溶质Al的浓度也越高，即＞＞＞。对于这一现象，可做如下解释：1) 过冷度ΔT增大，凝固点在相图上的液相线与固相线右移，凝固时平衡浓度增大，从而导致凝固过程固液界面溶质再分配时固相中的浓度增大，过冷度 ΔT越大，固相中溶质浓度越高；2) 虽然合金的凝固是非平衡结晶过程，但液固界面处的局部平衡状态是一直存在的，不管过冷度是大还是小，液固界面处固相侧溶质浓度和液相侧溶质浓度的比值都满足溶质平衡扩散系数k0，即k0=cS*/c*L，也就是说=，结合上面的＞＞＞，即可得＞＞＞；3) BOWER等[19]通过假设液相均匀混合固相有限扩散，推导出了液固界面处固相侧溶质浓度分布情况，即BF模型，其具体表达式为[20]

式中：为固相溶质浓度；为液相溶质浓度；fS为固相分数；α为固相反扩散因子，其值决定固相中溶质扩散状况；DS为固相溶质扩散系数；tf为局部凝固时间，λ为枝晶间距的一半；c0为初始溶质浓度，k0为溶质平衡分配系数。

由式(9)和(10)可知，局部凝固时间tf与α成正比，α参数决定固相中溶质的扩散状况，也就是说，局部时间tf越长，扩散时间越长，固相中溶质元素扩散越充分，溶质浓度越低。当过冷度 ΔT增大时，枝晶的生长速率也增大(见图4)，这使得局部的凝固时间tf减小，溶质扩散时间缩短，从而使得固相中溶质元素扩散越不充分，溶质浓度也就越高，从而有＞＞＞。

3.4 过冷度对枝晶尖端溶质扩散层厚度的影响

由图3可以看出，过冷度会对枝晶尖端前沿的溶质扩散层厚度δ产生影响。为更直观地显示出各枝晶尖端溶质扩散层厚度的关系，截取相同凝固时间(t=6 000Δt)不同过冷度下的枝晶尖端并进行同等倍数放大后再进行比较，如图5所示。由图5可见，过冷度ΔT越大，枝晶尖端前沿溶质扩散层越薄，即δ1＜δ2＜δ3＜δ4。这是由于镁合金等轴晶生长时，从固相中析出的溶质Al需通过溶质扩散层释放到液相中去，溶质扩散的长度λs=DL/v[21]。由于枝晶的生长速度v随过冷度ΔT的减小而减小(见图4)，故熔体中溶质扩散的长度λs随过冷度的减小而增大，这导致包围枝晶尖端的溶质扩散层厚度δ更大，从而有δ1＜δ2＜δ3＜δ4。

图5 不同过冷度下枝晶尖端溶质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 过冷度对溶质偏析比的影响

在通常的铸件生产中常采用偏析比SR来表示微观偏析的大小，偏析比在数值上等于枝晶间最大溶质浓度与枝晶干最小溶质浓度的比值。为更直观地显示过冷度对微观偏析程度的影响，计算了不同过冷度下溶质的偏析比曲线(见图6)。从图6中可明显看出，溶质的偏析比SR1＞SR2＞SR3＞SR4，即过冷度ΔT越大，微观偏析越严重。这是由于镁合金凝固过程中会通过溶质再分配向液相中析出溶质Al，而溶质Al在液相中的扩散速度小于枝晶生长速度，析出的Al来不及充分扩散从而富集在液固界面前沿。过冷度 ΔT越大，枝晶生长速度越快(见图4)，单位时间内析出到液相中的溶质就越多，二次枝晶也越发达，二次枝晶臂对富集溶质扩散的阻碍作用也越大，富集的溶质越难扩散，从而使得最高溶质浓度越大(见图7)。而固相中溶质扩散系数比液相中的要小得多，过冷度 ΔT对枝晶干最低溶质浓度的影响要小得多，因此，这最终导致溶质的偏析比越高，微观偏析越严重。

由以上分析可知，模拟结果与基于凝固理论分析得到的结论十分吻合，这也从另一方面验证了模型的可行性与正确性。

微观偏析的存在会降低合金的塑性、韧性和拉伸性能，因此，必须对微观偏析予以消除或者减轻。由模拟结果可知，减小过冷度可以降低合金微观偏析的程度，但这会使枝晶间距增大，材料力学性能下降，均匀化退火时间延长。合理的方法应该是降低初始温度，使过冷熔体快速凝固得到细小枝晶，然后再进行长时间的均匀化退火。

图6 不同过冷度下元素Al的偏析比曲线Fig. 6 Curves of segregation ratio of Al under different supercoolings

图7 不同过冷度下的溶质分布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 结论

1) 镁合金贫Al区主要集中于一次枝晶干和二次枝晶臂的中心，溶质富集于液固界面前沿，枝晶尖端浓度梯度最大。

4) 模拟结果和实验结果的对比以及基于凝固理论的分析验证了模型的可行性与正确性。

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] 袁训锋, 丁雨田. 强界面能各向异性下二元Ni-Cu合金枝晶生长过程的相场法模拟[J]. 中国有色金属学报, 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]. 金属学报, 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] 袁训锋, 丁雨田, 郭廷彪, 胡 勇. 强制对流作用下镁合金枝晶生长的相场法数值模拟[J]. 中国有色金属学报, 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]. 中国有色金属学报, 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] 朱昌盛, 王智平, 荆 涛, 肖荣振. 二元合金微观偏析的相场法数值模拟[J]. 物理学报, 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.