石玉峰，许庆彦，柳百成

(清华大学 机械工程系 先进成形制造教育部重点实验室，北京 100084)

通过数值模拟手段来模拟枝晶形貌及微观偏析，可以达到控制组织形态并且预测铸件性能的目的。目前在微观组织模拟方面主要的两种方法是元胞自动机(Cellular automaton)和相场法(Phase field)。其中，CA方法通过耦合实际温度场和溶质场，结合枝晶生长动力学理论以及随机形核与长大物理机理，能够描述晶粒度以及微观枝晶形貌，同时能够模拟确定铸造条件下的等轴晶生长、定向凝固柱状晶生长以及柱状晶-等轴晶转变等现象，因此，CA方法在最近的十几年以来得到了较为广泛的应用。很多学者对CA方法进行了不同的改进，出现了很多改进CA模型，已经可以用来模拟多种合金在不同铸造工艺条件下枝晶形貌的演变[1-8]。对于立方晶系合金来说，〈100〉是枝晶的择优生长取向，固液界面能在〈100〉方向上具有最低值，〈100〉方向上生长速率较高，故枝晶一次臂最容易沿着择优生长取向进行生长。实际枝晶的择优生长取向与坐标轴水平方向可能不平行，也即择优生长取向角在0°～90°之间，但是CA模型在模拟立方晶系合金的枝晶生长过程中，由于正交元胞带来的网格各向异性导致CA模型在描述任意择优生长取向角的枝晶生长时会出现计算偏差，从而在最终枝晶一次臂生长过程中逐渐丢失自身固有的择优生长取向，仅能沿着0°或 45°方向角生长。国内外许多学者通过对 CA捕获算法进行改进以消除网格各向异性的影响，GANDIN和RAPPZ[9]首次提出了偏心立方捕获算法，以最大限度消除网格各向异性对枝晶择优生长取向的影响，来实现枝晶一次臂沿枝晶择优生长取向的生长，但是偏心立方捕获算法需要人为对枝晶尖端进行实时修正，算法较为复杂；DONG和LEE[10]在CA-FD模型中采用偏心立方捕获算法[9]并耦合溶质场模拟了 Al-Cu合金定向凝固过程中的CET转变，再现了择优生长取向与坐标轴成多个不同夹角的枝晶生长形貌；LIU等[11]在CA模型中通过坐标旋转法将枝晶择优生长取向角旋转到任意值，模拟了Al-4.5%Cu合金的单个等轴枝晶形貌；于靖等[12]在CA模型中采用一种“父子单元法”的捕获算法，耦合溶质扩散方程计算了不同择优生长取向角的单枝晶生长形态；ZHU和STEFANESCU[13]在CA模型中使用界面前沿追踪方法来跟踪枝晶生长过程中的固液界面推移，模拟了Al-Cu合金多个不同择优生长取向角的枝晶形貌演化。

本文作者在已有CA模型的基础上，在计算界面生长速率方程时考虑了择优生长取向的长程作用，并对界面生长速率方程采用斜中心差分格式进行离散，从而建立了可模拟多个不同择优生长取向角的枝晶生长的MCA模型。选取NH4Cl-H2O系透明合金模拟了单个等轴枝晶和多个不同择优生长取向角的枝晶生长，模拟结果与实验结果符合较好；同时模拟了NH4Cl-H2O系透明合金的择优生长取向与坐标轴不平行时的柱状枝晶生长过程，并进行了实验验证。

1 数学模型和数值算法

1.1 模型假设

1)在本研究微观组织模拟MCA模型中，假设微观区域与外界无溶质交换，计算域中溶质质量守恒；

2) 凝固过程中的凝固速率较小，其溶质Péclet数小于 1，故可假设凝固过程中固液界面处于局部热力学平衡状态；同时由于动力学过冷一般很小，故本模型忽略动力学过冷；

3) 由于液相中溶质扩散系数比固相中的大几个数量级，因此忽略固相溶质扩散，凝固过程中仅考虑液相溶质扩散。

4) 二维计算域被划分成M×N个正方形元胞，元胞尺寸Δx=Δy；元胞状态的标识为元胞的固相分数fS：液态(fS=0)、固态(fS=1)和界面状态(0＜fS＜1)。

1.2 溶质传输控制方程

二元合金枝晶生长过程中，其液相和固相溶质传输方程形式分别为

式中：CL和DL分别是液相溶质成分和液相溶质扩散系数；k0为溶质平衡分配系数；fS是固相分数。式(1)右端第二项表示界面元胞由于固相分数的增加所导致的界面溶质再分配。

式中：CS和DS分别是固相溶质成分和固相溶质扩散系数。固液界面上的溶质再分配满足如下方程

1.3 界面热力学平衡

在二元合金枝晶生长过程中，由于界面元胞满足局部热力学平衡条件，因此界面元胞的平衡温度满足

式中：T为界面平衡温度，也即界面元胞的实际温度场；TL是初始溶质成分C0所对应的液相线温度；ΔTR和ΔTC分别表示曲率过冷和成分过冷。

1) 曲率过冷

在二维条件下，考虑Gibbs-Thomson效应的二元合金凝固界面曲率过冷表达式为[14]

式中：Г是Gibbs-Thomson系数；f(ψ)是立方晶系界面能各向异性函数；ψ是固液界面元胞单位法向矢量n所对应的平面角；K是界面曲率。其中二维条件下f(ψ)计算公式为

式中：ε为界面能各向异性系数；ψ0是枝晶的择优生长取向角；通过式(5)～(6)可以得到界面曲率过冷计算公式为

单位法向矢量n与界面元胞的固相分数fS有关，其表达式为

式中：nx和ny为法向矢量n在x和y方向的分量。

CA界面元胞曲率K受到元胞自身以及周围元胞的固相分数控制，二维条件下界面元胞曲率的计算方法见参考文献[15]。

2) 成分过冷

对于二元合金枝晶生长过程，由于界面前沿液相溶质扩散系数是有限值，因此在界面前沿溶质再分配造成的溶质富集形成了界面元胞的成分过冷，其表达式如下

式中：mL是液相线斜率；CL*是界面平衡溶质成分。在二元合金CA模型中假设液相线斜率为常数，不随温度和溶质成分变化。

3) 界面平衡溶质成分

通过联立式(5)和(10)同时根据式(3)可以得到界面平衡液相溶质成分CL*和固相溶质成分CS*的计算公式为

1.4 界面元胞固相分数增量

1) 界面位置矢量与择优生长取向之间的夹角β

枝晶生长过程中的择优生长取向与坐标轴x不一定重合，即 0°≤ψ0＜90°。在CA 枝晶生长模拟过程中，界面元胞的生长仅受到周围最近邻或者次近邻的4个元胞的影响，故随着枝晶生长，会逐渐丢失自身择优生长取向而趋向于沿着坐标轴或与坐标轴成 45°方向生长。为了消除这种网格各向异性带来的副作用，需要考虑界面元胞位置矢量与择优生长取向间的夹角β，β越接近0°或者90°，则表明界面元胞越接近择优生长取向，其生长越快。由此，当界面元胞与形核核心的位置越来越远时，可以避免枝晶逐渐丢失其自身固有的择优生长取向。图1所示为择优生长取向与界面元胞位置矢量夹角示意图。择优生长取向px=(cosψ0, sinψ0)，py=(-sinψ0, cosψ0)；计算域中枝晶形核核心为O(i0,j0)，界面元胞位置坐标为(i,j)，定义界面元胞与x轴方向夹角为α，其单位方向矢量p0=(cosα,sinα)，计算方法如下

式中：i0和j0为形核核心O的x和y方向坐标值，i和y是界面元胞P的坐标值。

图1 择优生长取向与界面元胞位置矢量夹角示意图Fig.1 Schematic illustration of intersection angle between preferred growth orientation and interface cell position vector

在计算任意位置界面元胞固相分数增加量的过程中，首先计算界面元胞位置矢量p0与择优生长取向px的夹角β，其计算公式如下：

2) 界面生长速率

二元合金枝晶固液界面上存在溶质质量守恒，为了消除网格各向异性对界面元胞法向生长速率的影响，需要考虑择优生长取向的作用，故界面元胞法向生长速度w的表达式为

式中：ψ是界面元胞的法向矢量和x轴的夹角。

3) 界面元胞固相分数增量计算

界面元胞固相分数增量 ΔfS与固液界面法向生长速率w之间的关系式如下[2]

式中：L是考虑择优生长取向作用下的界面元胞空间步长；Δt是时间步长；fSn和fSn+1分别表示当前时刻和下一时刻界面元胞的固相分数。当界面元胞固相分数fS达到1时，开始对周围液相元胞进行捕获，具体捕获规则见参考文献[12]。

1.5 数值求解方法

在MCA模型求解过程中，对界面生长速率方程(16)中的一阶导数项采用斜中心差分格式进行离散，如图2所示为考虑择优生长取向角ψ0后的斜中心差分格式示意图，其中d是元胞P的空间步长，其大小与Δx相等。

图2 斜中心差分格式示意图Fig.2 Schematic diagrams of central difference scheme:(a) 0°≤ψ0≤45°; (b) 45°＜ψ0＜90°

根据图2，可得到方程(16)中一阶导数项的差分离散形式为式中：CL*是元胞P的界面平衡液相溶质成分；CL和CS是液相和固相溶质成分；N、S、W、E分别为元胞P上、下、左、右4个方向的邻居元胞；d0是斜中心差分格式下的元胞步长，当0°≤ψ0≤45°时，d0=dsecψ0；当45°＜ψ0＜90°时，d0=dcscψ0。

当0°≤ψ0≤45°时，斜中心差分格式的E、W、N、S 4个元胞的液相溶质成分和固相分数的计算公式如下

对于斜中心差分格式的E、W、N、S 4个元胞的固相溶质成分计算方法同式(20)。

当45°＜ψ0＜90°时，斜中心差分格式的 E、W、N、S 4个元胞的液相溶质成分和固相分数的计算方法与式(20)和(21)有所不同，离散形式为

对于斜中心差分格式的E、W、N、S 4个元胞的固相溶质成分计算方法同式(22)。

MCA模型的求解流程如下所示：

1) 计算域划分成M×N个正方形元胞，元胞状态包含元胞大小、初始溶质成分和初始温度；在计算域中形成若干具有一定择优生长取向角ψ0的固相晶核，晶核捕获周围的液相邻居元胞，使之成为界面元胞；

2) 确定计算域时间步长，Δt=Δx2/(5DL)；

3) 假设微观计算域初始温度均匀，其值为T0，且计算域以一定的冷却速率CR下降，可以计算元胞任意时刻的界面平衡温度T*。根据式(4)～(12)计算得到界面平衡液相溶质成分CL*和固相溶质成分CS*，其中方程(9)的一阶导数项差分格式为

离散后的E、W、N、S 4个元胞的固相分数计算方法见式(21)和(23)。

4) 通过斜中心差分格式求解方程(16)，得到界面法向生长速率w；

5) 通过方程(17)计算界面元胞固相分数增量ΔfS；

6) 求解方程(1)～(3)，获得液相、界面和固相元胞中的溶质成分分布；

7) 判断界面元胞固相分数fS，若fS小于1，则重复步骤(3)～(6)；若fS等于 1，则转变成固相元胞，并对周围液相元胞进行捕获。

2 模拟结果与实验结果

2.1 氯化铵水溶液单个等轴晶生长过程模拟

实际合金凝固过程的直接实时观察较为困难且实验成本很高，而可在低温下凝固的透明合金具有和实际合金类似的凝固潜热、枝晶形貌特性和动力学行为，可代替实际合金来研究枝晶生长形貌演化机理。为了验证本研究的MCA模型，选取NH4Cl-70%H2O(质量分数)溶液进行枝晶生长实验。实验过程中的冷却速率为2.0 K/s，当温度低于液相线温度TL以下1 K时，在显微镜视野中产生晶核并开始长大。

模拟采用与实验相同的冷却速率，在选定的计算域中心形核，CA 元胞步长为 3 μm，计算域大小为1 500 μm×1 500 μm，枝晶择优生长取向角ψ0为 0°。其余参数见表1。

如图3(a)～(c)所示为模拟的NH4Cl-70%H2O透明合金单个等轴晶的形貌演化过程，图3(d)所示为光学显微镜照片。从模拟结果可以看到，等轴晶形貌具有四重对称性特征，符合立方晶系枝晶形貌特征，4个一次枝晶臂之间夹角互为 90°，且枝晶的择优生长取向与x轴夹角为0°，即枝晶沿着初始确定的择优生长取向生长。在凝固过程中，由于枝晶间溶质富集随着过冷度增大而增加，导致一次枝晶臂上逐渐长出二次枝晶臂，与一次枝晶臂之间的夹角保持 90°。模拟得到的二次枝晶臂间距为51 μm，实验得到的二次枝晶臂间距为 44.5 μm，模拟结果与实验结果基本一致，从而说明MCA模型是可靠的。

表1 NH4Cl-70%H2O溶液的热物性参数[16-18]Table1 Thermal physical properties of NH4Cl-70%H2O solution[16-18]

图3 NH4Cl-70%H2O透明合金单枝晶生长模拟结果和实验结果Fig.3 Simulated dendrite growth and experimental results of NH4Cl-70%H2O transparent alloy: (a) fS =0.5%; (b) fS=2.5%;(c) fS=9.0%; (d) Experimental result

2.2 不同择优生长取向角枝晶生长模拟及实验结果

采用本研究的MCA模型，在计算界面元胞生长过程中考虑枝晶择优生长取向的作用，可以最大限度地消除网格各向异性的影响，实现立方晶系合金枝晶生长过程中保证自身的择优生长取向不变。本节针对NH4Cl-74%H2O透明合金，进行枝晶生长过程的数值模拟，凝固过程的热物性参数见文献[15]。选取计算域大小为1.8 mm×2.4 mm，CA 元胞尺寸为5 μm，计算域冷却速率CR为2 K/s，当温度降低到液相线温度TL以下5 K时，计算区域形核。

图4(a)～(c)所示为模拟具有不同择优生长取向角的枝晶形貌演化，其择优生长取向角ψ0分别为50°、0°、75°和60°(自左至右)。如图4(a)所示，生长初期二次枝晶臂并不发达，随着生长的继续，枝晶界面前沿和枝晶间溶质逐渐富集，如图4(b)和4(c)所示，导致各处成分过冷出现差别，从而二次枝晶臂逐渐出现。从图4(a)～(c)可以看出，各个枝晶在生长过程中自身的择优生长取向并没有因网格各向异性的影响而改变，一次枝晶臂始终沿着择优生长取向生长。如图4(d)所示为实验获得的NH4Cl-74%H2O透明合金枝晶结果，经过测量可得到实验获得的四枝晶的择优生长取向角ψ0分别为 48°、0°、78°和63°，两者吻合较好。模拟过程中4个枝晶的一次枝晶臂上出现的二次枝晶臂与一次臂主干的夹角均保持 90°，符合 NH4Cl-74%H2O透明合金枝晶特征，与实验结果一致，说明本研究的MCA模型可用来模拟不同择优生长取向角的多个枝晶的生长形貌演化。

2.3 定向凝固柱状晶生长过程数值模拟与实验验证

本节对NH4Cl-74%H2O透明合金柱状枝晶定向凝固过程进行模拟，凝固过程中选取与实验过程一致的温度梯度G=1 K/mm，方向沿着竖直方向。模拟采用CA元胞数量为720×480，元胞尺寸为5 μm，热物性参数取值见文献[15]。当计算域温度降低到液相线温度TL以下5 K时，计算域底部形核，共产生9个形核核心，具有相同的择优生长取向角ψ0=60°。在抽拉速度vp=8 μm/s的作用下，柱状枝晶开始生长。

图5所示为模拟得到的NH4Cl-74%H2O透明合金定向凝固柱状枝晶溶质场分布的模拟结果和实验结果。其中柱状枝晶的择优生长取向角ψ0为60°，在温度梯度G和抽拉速度vp的作用下，枝晶界面前沿过冷度不断增大；如图5(a)和5(c)所示，枝晶间和枝晶界面前沿随着过冷度的增大出现了溶质富集，由于成分过冷和曲率过冷的作用，柱状枝晶一次枝晶臂上出现二次枝晶臂甚至三次枝晶臂，且三次枝晶臂与二次枝晶臂之间以及二次枝晶臂与一次枝晶臂之间的夹角均保持 90°，符合立方晶系合金特征。柱状枝晶在生长过程中，一次枝晶臂沿着择优生长取向方向生长速度最快，图5(b)和5(d)的实验结果证明了定向凝固过程中柱状枝晶一次枝晶臂始终沿着择优生长取向的方向生长，与水平方向夹角保持 60°，与模拟结果吻合较好。同时，实验过程中测量的平均一次枝晶臂间距为32 μm，与模拟得到的平均一次枝晶臂间距36 μm相差较小，说明MCA模型可以预测定向凝固过程中择优生长取向与坐标轴不平行时的柱状枝晶的形貌演化。

图4 NH4Cl-74%H2O透明合金不同择优生长取向角的多个枝晶生长的模拟结果和实验结果Fig.4 Simulated and experimental results of multiple dendrite growth of NH4Cl-74%H2O transparent alloy with random preferred growth orientation angles: (a) t=1.2 s, fS=1.5%; (b) t=2.4 s, fS=5.5%; (c) t=3.4 s, fS=12.0%; (d) OM photograph

图5 NH4Cl-74%H2O透明合金定向凝固柱状枝晶形貌的模拟和实验结果Fig.5 Simulated and experimental columnar dendrite morphologies of NH4Cl-74%H2O transparent alloy directionally solidified with constant preferred growth orientation angles: (a) Simulated result, t=200 s, G=1 K/mm, vp=8 μm/s; (b) Experimental result;(c) Simulated result, t=270 s, G=1 K/mm, vp=8 μm/s; (d) Experimental result

3 结论

1) 考虑枝晶择优生长取向角的随机性对枝晶界面法向生长速率的影响，建立了MCA模型以模拟不同择优生长取向角的枝晶形貌演化。

2) 使用MCA模型对NH4Cl-70%H2O溶液单个等轴枝晶的生长过程进行模拟，枝晶的择优生长取向与坐标轴平行，模拟结果很好地再现了一次枝晶臂的生长以及二次枝晶臂的产生与发展。同时，模拟和实验得到的二次枝晶臂间距吻合较好，枝晶形貌也基本一致。

3) 选择 NH4Cl-74%H2O溶液进行多个等轴枝晶自由生长的实验与模拟。模拟得到的不同择优生长取向的等轴枝晶形貌与实验结果吻合较好，说明 MCA模型可以模拟不同择优生长取向角的立方晶系合金的枝晶生长过程。

4) 采用MCA模型对NH4Cl-74%H2O溶液定向凝固柱状枝晶形貌演化进行了模拟，模拟结果与实验结果基本一致，说明MCA模型可以模拟一定择优生长取向角的柱状枝晶的形貌演化。

[1]NASTAC L.Numerical modeling of solidification morphologies and segregation patterns in cast dendritic alloys[J].Acta Materialia, 1999, 47: 4253-4262.

[2]BELTRAN-SANCHEZ L, STEFANESCU D M.A quantitative dendrite growth model and analysis of stability concepts[J].Metall Mater Trans A, 2004, 35: 2471-2485.

[3]WANG W, LEE P D, MCLEAN M.A model of solidification microstructures in nickel-based superalloys: Predicting primary dendrite spacing selection[J].Acta Materialia, 2003, 51:2971-2987.

[4]LIU B C, XU Q Y, JING T, SHEN H F, HAN Z Q.Advances in multi-scale modeling of solidification and casting processes[J].JOM, 2011, 63(4): 19-25.

[5]LIU X B, XU Q Y, JING T, LIU B.Microstructure of aluminum twin-roll casting based on cellular automaton[J].Transactions of Nonferrous Metals Society of China, 2008, 18(4): 944-948.

[6]ZHU M F, HONG C P, STEFANESCU D M, CHANG Y A.Computational modeling of microstructure evolution in solidification of aluminum alloys[J].Metall Mater Trans B, 2007,38: 517-524.

[7]LIU B C, XIONG S M, XU Q Y.Study on macro- and micromodeling of the solidi fi cation process of aluminum shape casting[J].Metall Mater Trans B, 2007, 38: 525-532.

[8]李 强, 李殿中, 钱百年.应用连续性方法模拟枝晶生长[J].金属学报, 2004, 40(6): 634-638.LI Qiang, LI Dian-zhong, QIAN Bai-nian.Modelling of the dendrite growth by using continuous method[J].Acta Metall Sin,2004, 40(6): 634-638.

[9]GANDIN C A, RAPPAZ M.A 3D cellular automaton algorithm for the prediction of dendritic grain growth[J].Acta Materialia,1997, 45: 2187-2195.

[10]DONG H B, LEE P D.Simulation of the columnar-to-equiaxed transition in directionally solidi fi ed Al-Cu alloys[J].Acta Materialia, 2005, 53: 659-668.

[11]LIU Y, XU Q Y, LIU B C.A modified cellular automaton method for the modeling of the dendritic morphology of binary alloys[J].Tsinghua Science and Technology, 2006, 11(5): 495-500.

[12]于 靖, 许庆彦, 崔 锴, 柳百成.基于一种改进 CA模型的微观组织模拟[J].金属学报，2007, 43(7): 731-738.YU Jing, XU Qing-yan, CUI Kai, LIU Bai-cheng.Numerical simulation of microstructure evolution based on a modified CA method[J].Acta Metall Sin, 2007, 43(7): 731-738.

[13]ZHU M F, STEFANESCU D M.Virtual front tracking model for the quantitative modeling of dendritic growth in solidification of alloys[J].Acta Materialia, 2007, 55: 1741-1755.

[14]MICHELIC S C, THUSWALDNER J M, BERNHARD C.Polydimensional modelling of dendritic growth and microsegregation in multicomponent alloys[J].Acta Materialia,2010, 58: 2738-2751.

[15]石玉峰, 许庆彦, 龚 铭, 柳百成.定向凝固过程中NH4Cl-H2O枝晶生长的数值模拟[J].金属学报, 2011, 47(5):620-627.SHI Yu-feng, XU Qing-yan, GONG Ming, LIU Bai-cheng.Simulation of NH4Cl-H2O dendritic growth in directional solidification[J].Acta Metall Sin, 2011, 47(5): 620-627.

[16]HANSEN G, LIU S, LU S Z, HELLAWELL A.Dendritic array growth in the systems NH4Cl-H2O and [CH2CN]2-H2O: Steady state measurements and analysis[J].Journal of Crystal Growth,2002, 234: 731-739.

[17]BENNON W D, INCROPERA F P.The evolution of macrosegregation in statically cast binary ingots[J].Metall Mater Trans B, 1987, 18: 611-616.

[18]冯妍会, 聂 红, 张欣欣.定向凝固宏微观热质传递及双扩散流作用[J].工程热物理学报, 2008, 29(2): 301-305.FENG Yan-hui, NIE Hong, ZHANG Xin-xin.Coupled macro-micro transports with double diffusive flow during directional solidification[J].Journal of Engineering Thermophysics,2008, 29(2): 301-305.