随着相场方法在单相凝固体系中的运用，其在多相凝固系统中的应用也越来越广泛，通常借助多相场模型解决多相凝固问题，如共晶凝固[1]、包晶凝固[2]等。 一般地，使用多个相场变量来描述系统的多个组成相，故建立了多相场模型。 多相场模型将序参量视为各相的体积分数进而定义一组多组元序参量，并以约束形式限制这组序参量关系[3]。 对于相场模型来说，多相的化学自由能输入是计算体系热力学驱动力的关键，其在模型中的表达方式可有多种形式，如伪自由能多项式[4]、依据热力学数据进行多项式拟合[5]或者直接将热力学数据库中数据导入到相场模型[6]。 相对于稀溶液模型，开发非稀溶液模型也将为创建多相化学自由能模型带来一定困难；另外，多相场模型内多相界面的处理也是关键环节之一，当前主流处理方法如Wheeler-Boettinger-Mc-Fadden(WBM)模型[7]、Kim-Kim-Suzuki(KKS)模型[8]及Steinbach 模型[9]等。 考虑到相场界面宽度增加所带来的溶质截流现象，在成分场控制方程中需添加反溶质截流项[10]；此外，模型还应考察冷速、合金过冷度等主要控制参数对凝固过程的影响规律[11]。 若综合考虑上述多种关键模型元素，将为创建某合金体系下的多相凝固模型带来极大挑战。 本文以钛铝(TiAl)合金为研究体系，鉴于其具有比强度高、高温性能优异等特点，且密度仅为镍基高温合金的1/2，可替代高温合金制作航空航天转动或往复运动结构件,实现推重比及燃油效率大幅度提高[12]。 考虑到包晶相变使钛铝合金在铸造过程中容易发生成分偏析、组织粗大等现象，将导致塑性较差且各向异性明显进而显著影响其使役性能；且包晶转变温度高且反应速度快，通过实验手段很难进行观察与表征， 采用相场模拟无疑是探索包晶转变过程的有效手段。
目前， 国内外研究包晶相变开展的相场模拟较少，如Wang 等[13]利用相场法研究了Fe-Ni 合金内包晶反应机制；Pan 等[14]与Cai 等[2]模拟了Fe-C 体系内铁素体(δ)与液相(L)发生包晶反应生成奥氏体(γ)的过程，即δ+L→γ 转变；Luo 等[15]模拟了Fe-C 合金中包晶相γ 在δ/L 界面处形核， 然后围绕δ 相生长并增厚。 近期，Yang 等[16]建立了多相场模型并实现对低碳钢中包晶转变的模拟；Shen 等[17]基于相场法研究了高锰钢中的包晶转变；Parida 等[18]在模拟Fe-C 系包晶相变时还添加了δ/γ 之间弹性作用的影响；在钛铝合金包晶反应方面，即β+L→α；Eiken 等[19]借助MICRESS 软件半定量地模拟了Ti-45Al 合金中包晶α 相的形成过程；Viarin 等[20]模拟了发生包晶相变后的重熔过程， 提出了包晶反应逆过程可诱发断裂的新机制。总之，国内外已经开展的包晶相变相场模拟工作主要集中在Fe-C 系包晶合金的研究，而其它合金体系，如钛铝合金[21]，相关包晶相变模拟研究极为有限[22-23]，有待开展系统研究。
本文拟构建可模拟钛铝合金中包晶相变过程的多相场模型， 旨在明确钛铝合金铸造过程中包晶相变的限制性环节， 为有益调控包晶相变过程提供模型基础及数据支撑。
1 模拟方法
1.1 多相场模型
本文基于Steinbach 等[9,24]提出的多相场模型，Ginzburg-Landau 型自由能泛函采用Ohno 等[25]提出的改进形式，即：
式中，ε 为梯度能系数；ω({准i})为双阱势的高度；gi({准i})为插值函数；fi(ci,T)为i 相的吉布斯自由能。
此外，相场控制方程为：
式中，准i(i=1,2,3)分别表示体系任一位置上α、β 和L 相的体积分数；si 为阶跃函数， 表示i 相在系统中某点的存在状态(当i 相存在时，si 被赋值为1，否则赋值为0)；N
为系统中某一点存在相的总数；n 为总相数， 在包晶相变的相场模拟时，n=3；Mij 为相场动力学系数。
成分场的扩散满足守恒条件， 为了消除界面厚度增加所带来的溶质截留效应， 在成分场的演化方程中引入反溶质截留项，具体的控制方程如下：
式中，Di 为扩散系数；aij 为反溶质截留项系数；ni 为i相等值线外法线的单位向量，其定义为ni=-塄准i/|塄准i|。
1.2 自由能输入
对于Ti-45Al 合金， 各相的吉布斯自由能计算不能采用传统的稀溶液模型， 通常采用正规溶液模型或亚正规溶液模型。 通过Pandat 热力学计算软件[26]及Ti-Al 合金热力学数据库[27]获取所需的定量化热力学数据，计算得到的部分相图如图1 所示。图中的实线为平衡成分，短划线为平衡成分的扩展。从计算结果可知，包晶相变温度为Tp=1 763.9 K，此时对应的L、β 和α 相的平衡成分分别为0.496、0.446 及0.467。
图1 钛铝二元合金部分相图
Fig.1 Partial phase diagram of the Ti-Al binary alloy
1.3 模拟参数设置
相场动力学系数Mij 表达式如下[25]：
式中，a2 为常数，这里取0.627 6[28]；Tm,ij 为纯组元发生i-j 相转变时的温度， 可以根据图1 中某一温度下平衡成分和该点的斜率求得；wij 为i-j 相的双阱势的高度，其值与单相场模型相同；Vm 为摩尔体积；kij 为平衡分配系数
为平衡成分。
对式(2)和(3)采用标准的二阶有限差分进行离散，且相场和成分场均采用零通量边界条件，主要的模拟参数如表1 所示。
表1 多相场模型使用的模拟参数
Tab.1 Simulation parameters used in the multiphase-field model
DefinitionunitValue Grid size dxμm0.01 Time step dts1.0×10-9 Interface energy σβL[19,29]J/m20.10 Interface energy σαL[19]J/m20.10 Interface energy σαβ[19]J/m20.10 Diffusion coefficient in the α phase Dα[30]m2/s6.6×10-3×exp[(-328 971.020 4/(RT)]Diffusion coefficient in the β phase Dβ[30]m2/s1.94×10-4×exp[(-32 4147.398 4/(RT)]×exp[(122 134 109.04/(R×T2)]Diffusion coefficient in the L phase DL[29,31]m2/s1.0×103×max(Dα,Dβ)
2 实验结果及讨论
2.1 自由能曲线的二阶抛物线近似及其验证
Ti-Al 合金包晶凝固的多相场模拟， 在忽略了α/β 界面的弹性相互作用， 最关键的问题是界面区域的溶质分布和相变驱动力的计算。KKS 模型根据准平衡假设鄣fL/鄣cL=鄣fβ/鄣cβ=鄣fα/鄣cα=μc 和混合规则c=准LcL+准βcβ+准αcα，可计算虚拟相组分，进而求出相变驱动力[8]。
对于正规溶液模型或亚正规溶液模型，KKS 相场模型在求解虚拟相成分过程中需要求解复杂的非线性方程组。 对于高阶非线性方程组一般采用牛顿-拉夫逊迭代法，这极大的限制了数值模拟效率。 为了提高计算效率， 在求解虚拟相组分时对吉布斯自由能曲线进行二阶抛物线近似， 将求解非线性方程组转换为求解线性方程组，极大地提高了计算效率。本文采用Pandat 热力学计算软件获取不同温度条件下L、β 及α 相的自由能曲线，并对自由能曲线进行二阶抛物线拟合，其结果如图2 所示。从图中可以看出，拟合曲线和原始数据点基本完全重合，且所有曲线的拟合精度均为R2≥0.99。
图2 Ti-Al 合金不同温度下各相吉布斯自由能曲线：(a)1 760 K;(b)1 750 K
Fig.2 Gibbs free energy curves of phases in Ti-Al alloys at different temperatures:(a)1 760 K;(b)1 750 K
在模拟过程中， 为了确保采用二阶抛物线近似算法的可行性， 比较了传统的牛顿-拉夫逊迭代法和二阶抛物线近似算法求得的虚拟相组分和相变驱动力，其结果如图3 所示。 图中的实线为牛顿-拉夫逊迭代法计算得到的虚拟相组分和相变驱动力，短划线为二阶抛物线近似算法所求得的虚拟相组分和相变驱动力。对比结果表明，这两种方法计算得到的虚拟相组分和相变驱动力在扩散势合理变化区间内高度一致。 因此采用二阶抛物线近似算法不仅可以保证计算的准确性，还可以提高计算效率，节省计算时间。
图3 牛顿-拉夫逊迭代法和二阶抛物线近似算法对比：(a)相组分；(b)相变驱动力
Fig.3 Comparison between the Newton-Raphson iteration method and second-order parabola approximation method:(a)phase constituent;(b)phase transformation driving force
2.2 过冷度对一维包晶转变平界面运动的影响
在一维等温包晶转变模拟过程中， 采用的过冷度分别为ΔT=3,5 及10 K。 为了避免系统尺寸对界面运动的影响，使用了相对较大的系统尺寸lx=10 μm。设置初始计算域左侧为β 相，右侧为L 相，中间为α 相(初始厚度为0.2 μm)。 β、L 和α 相初始浓度分别对应不同温度条件下的平衡成分cLβ,e、cβL,e 和cLα,e。
图4 给出了1 760.9 K(ΔT=3 K)条件下一维包晶转变的模拟结果，其中图4a 为α 相的相场变量，图4b 为成分分布曲线。 从图中可以看出，α 相同时向L 相和β 相中生长，且α/L 界面移动速率远高于α/β 界面的移动速率， 这与Fe-Ni 合金包晶转变的模拟相一致，与Fe-C 合金的模拟结果相反[14]。 在模拟时间达到10 ms 时，α/L 界面几乎停止迁移，而α/β 界面仍以相对稳定的速度进行迁移。 对相应的成分场进行分析， 由于液相中的扩散系数远大于α相和β 相中的扩散系数，在10 ms 以后α/L 界面成分近乎达到两相的平衡状态， 因此α/L 界面几乎停止演化。
图4 T=1 760.9 K 下一维包晶转变模拟结果：(a)相场；(b)成分场
Fig.4 Simulation results of one-dimensional peritectic transformation at 1 760.9 K:(a)phase field;(b)composition field
对包晶转变过程平界面的运动进行定量分析，其结果如图5 所示。从图5a 中可以看出，i-j 界面的移动距离xij 与模拟时间t 遵循关系，其中Aij 为抛物线速率常数，可见包晶转变过程是受扩散控制的。在不同过冷度条件下， 计算得到α/L 界面和α/β 界面的抛物线速率常数，结果如图5b 所示。可以看出，随着过冷度的增加，α/L 界面和α/β 界面的抛物线速率常数逐渐增大， 且α/L 界面抛物线速率常数的增大速率远大于α/β 界面。
图5 包晶转变过程平界面运动的定量分析：(a)T=1 760.9 K 下各界面移动距离与时间的关系；(b)各界面抛物线速率常数与过冷度的关系
Fig.5 Quantitative analysis of planar interface motion during peritectic transformation:(a)relationship between the moving distance and time of each interface at 1 760.9 K;(b)relationship between the parabolic rate constant of each interface and the degree of supercooling
2.3 过冷度对包晶反应的影响
借助多相场模型， 在4 μm×10 μm 的二维计算区域内对Ti-45Al 合金单个片层进行了不同过冷度条件下的等温包晶相变模拟。 模拟开始时， 设置β相和L 相共存，初始成分分别为cLβ,e 和cβL,e。 将1 个半径(R0)为0.2 μm 的α 相晶核置于计算区域底部的β/L 界面上，初始成分为cLα,e。图6 给出了ΔT=3 K条件下α 相的生长形貌。 从图中可以看出，α 相沿β/L 界面生长，逐渐将L 相与β 相分开。 与此同时，α 相沿x 方向也不断增厚， 但由于α 相向L 相和β相生长速度不同，表现出了非对称性生长。
图6 T=1 760.9 K 下α 片层的生长形貌
Fig.6 Growth morphology of the α plate at 1 760.9 K
图7 给出了与图6 相对应的α 相生长过程中Al 原子的溶质分布情况。 从图中可以看出，溶质场形貌与相场相匹配。 在α 相生长过程中，Al 原子从α 相再分配到液相中， 所以溶质在α/L 界面L 相一侧富集。然而α/β 界面侧则相反，β 相向α 相转变过程中，需要消耗β 相中的Al，从而导致α/β 界面β相侧溶质含量降低。但由于溶质Al 原子在α 相和β相中扩散较慢， 所以α/β 界面侧溶质含量变化并不明显。
图7 T=1 760.9 K 下α 片层生长时的溶质分布
Fig.7 Solute distribution during the growth of the α plate at 1 760.9 K
图8 给出了不同过冷度条件下包晶反应三相区附近界面形貌。 从图中可以看出α 片层相对于初始的β/L 界面(x=0 处红色点划线)呈现处明显的非对称特征。 且在过冷度ΔT=3 K 时，三相区前沿可以观察到β/L 界面向β 相侧发生偏移，这种偏移是由于α 相生长过程中排出溶质原子使该区域溶质浓度升高，导致α 片层前方的β 相重熔造成的。 然而，随着过冷度的增加，这种偏移现象逐渐消失，但三相区形状基本维持不变。 这与Ohno 等[32]的研究结果相一致，三相交界处的局部形状主要由界面能(也即界面张力)之间的平衡决定，与界面移动速度的差异无关。
图8 不同过冷度条件下包晶反应三相区附近界面形貌：(a)ΔT=3 K;(b)ΔT=5 K;(c)ΔT=10 K
Fig.8 Interface morphology near the three-phase region of the peritectic reaction under different undercooling conditions:(a)ΔT=3 K;(b)ΔT=5 K;(c)ΔT=10 K
3 结论
(1)在求解相场界面处虚拟相成分时，对吉布斯自由能曲线进行二阶抛物线近似， 将求解非线性方程组转换为求解线性方程组， 极大地提高了计算效率。 该策略将为后续大尺寸体系的包晶相变模拟提供有益思路。
(2)发生包晶转变时，α/L、α/β 相界面的移动距离xij 与时间t 遵循xij=Aijt1/2 关系，表明包晶转变过程受扩散控制。 随着过冷度的增大，α/L 界面和α/β 界面的迁移速率均增大， 特别是α/L 界面的迁移速率增幅更为显著。
(3)在包晶反应过程中，α 相沿着β/L 界面快速增长。由于α 相向L 相和β 相的生长速度存在差异，α 相片层表现出了非对称性生长；较低过冷度下，α 片层前端出现了β 相的重熔。 随着过冷度增大，重熔现象逐渐消失， 但是三相区的形状基本保持不变。 可见， 三相交界处的形状主要由界面能(也即界面张力)之间的平衡决定，与界面迁移速率的差异无关。
[1] KIM S G,KIM W T,SUZUKI T,ODE M.Phase-field modeling of eutectic solidification[J].Journal of Crystal Growth,2004,261(1):135-158.
[2] CAI Y H,WANG F,ZHANG Z H,NESTLER B.Phase-field investigation on the peritectic transition in Fe-C system[J].Acta Materialia,2021,219:117223.
[3] YANG C,WANG X T,WANG J S,HUANG H B.Multiphase-field approach with parabolic approximation scheme[J].Computational Materials Science,2020,172:109322.
[4] CHEN L Q.A novel computer simulation technique for modeling grain growth[J].Scripta Metallurgica et Materialia, 1995, 32(1):115-120.
[5] ZHANG J H,TENG C Y,MENG Z C,XU H S,YANG L,XU D S,YANG R.Selection and mechanical evaluation of γ/γ boundary in γ-TiAl alloy[J].Intermetallics,2020,126:106946.
[6] CHEN Q,MA NN,WU K S,WANG Y Z.Quantitative phase field modeling of diffusion-controlled precipitate growth and dissolution in Ti-Al-V[J].Scripta Materialia,2004,50(4):471-476.
[7] WHEELER A A, BOETTINGER W J, MCFADDEN G B.Phasefield model for isothermal phase transitions in binary alloys[J].Physical Review A,1992,45:7424-7439.
[8] KIM S G,KIM W T,SUZUKI T.Phase-field model for binary alloys[J].Physical Review E,1999,60:7186-7197.
[9] STEINBACH I, PEZZOLLA F, NESTLER B, SEE覻ELBERG M,PRIELER R,SCHMITZ G J,REZENDE J L L.A phase field concept for multiphase systems[J].Physica D:Nonlinear Phennomena,1996,94(3):135-147.
[10] KAVOUSI S, ZAEEM M A.Quantitative phase-field modeling of solute trapping in rapid solidification[J].Acta Materialia, 2021,205:116562.
[11] XU X J,LIN J P,XIANG C,LIANG Y F.Element distribution and its induced peritectic reaction during solidification of Ti-Al-Nb alloys[J].Metals,2021,11(9):1386.
[12] 杨锐.钛铝金属间化合物的进展与挑战[J].金属学报，2015,51(2):129-147.YANG R.Advances and challenges of TiAl base alloys[J].Acta Metallurgica Sinica,2015,51(2):129-147.
[13] WANG J H,PAN S Y,LI Y S,SHEN X P,ZHANG Q Y,JIA D.A phase-field study of the peritectic reaction mechanisms in Fe-Ni alloys[J].Computational Materials Science,2023,230:112491.
[14] PAN S Y, ZHU M F, RETTENMAYR M.A phase-field study on the peritectic phase transition in Fe-C alloys[J].Acta Materialia,2017,132:565-575.
[15] LUO S,LIU G G,WANG P,WANG X H,WANG W L,ZHU M Y.In situ observation and phase-field modeling of peritectic solidification of low-carbon steel[J].Metallurgical and Materials Transactions A,2020,51:767-777.
[16] YANG Y M, LUO S, WANG P, WANG W L, ZHU M Y.Multiphase field modeling of dendritic solidification of low-carbon steel with peritectic phase transition[J].Metallurgical and Materials Transactions B,2021,52:3708-3719.
[17] SHEN Y Z,LIU J H,XU H.Study of peritectic phase transition in high-Mn steel using phase-field method [J].Metallurgical and Materials Transactions B,2022,53:121-135.
[18] PARIDA H,KUNNDINN J,ALVES C L M.Study of the peritectic phase transformation kinetics with elastic effect in the Fe-C system by quantitative phase-field modeling[J].Computational Materials Science,2023,224:112160.
[19] EIKEN J,APEL M,WITUSIEWICZ V T,ZOLLINGER J,HECHE U.Interplay between α(Ti)nucleation and growth during peritectic solidification investigated by phase-field simulations[J].Journal of Physics:Condensed Matter,2009,21:464104.
[20] VIARDIN A,BOUSSINOT G,ZOLLINGER J.Phase field modeling of partial remelting during reheating of a multiphase peritectic solidification microstructure[J].Materialia,2022,26:101590.
[21] 傅恒志，李新中，刘畅，苏彦庆，李双明，郭景杰.Ti-Al 包晶合金定向凝固及组织选择[J].中国有色金属学报，2005(4):495-505.FU H Z,LI X Z,LIU C,SU Y Q,LI S M,GUO J J.Directional solidification and microstructure selection for Ti-Al peritectic alloy[J].The Chinese Journal of Nonferrous Metals,2005(4):495-505.
[22] 李新中，苏彦庆，郭景杰，吴士平，傅恒志.定向凝固包晶相变微观组织演化的相场方法研究I.三相交节点的延伸[J].金属学报，2006,42:599-605.LI X Z,SU Y Q,GUO J J,WU S P,FU H Z.Phase-field research of microstructure evolution for directionally solidified peritectic transition I.Extension of Trijunction[J].Acta Metallurgica Sinica,2006,42:599-605.
[23] 苏彦庆，李新中，郭景杰，吴士平，傅恒志.定向凝固包晶相变微观组织演化的相场方法研究II.形核控制的微观组织模拟[J].金属学报，2006,42:606-610.SU Y Q,LI X Z,GUO J J,WU S P,FU H Z.Phase-field research of microstructure evolution for directionally solidified peritectic transition II.Simulation of Nucleation-Controlled Microstructure[J].Acta Metallurgica Sinica,2006,42:606-610.
[24] STEINBACH I, PEZZOLLA F.A generalized field method for multiphase transformations using interface fields[J].Physica D:Nonlinear Phenomena,1999,134(4):385-393.
[25] OHNO M,MATSUURA K.Quantitative phase-field modeling for two-phase solidification process involving diffusion in the solid[J].Acta Materialia,2010,58(17):5749-5758.
[26] CHEN S L,ZHANG F,XIE F Y,DANIEL S,YAN X Y,CHANG Y A, SCHMID-FETZER R, OATES W A.Calculating phase diagrams using PANDAT and panengine[J].JOM,2003,55:48-51.
[27] WITUSIEWICZ V T, BONDAR A A, HECHT U, REX S, VELIKANOVA T Y.The Al-B-Nb-Ti system: III.Thermodynamic re-evaluation of the constituent binary system Al-Ti[J].Journal of Alloys and Compounds,2008,465(1-2):64-77.
[28] OHNO M,MATSUURA K.Quantitative phase-field modeling for dilute alloy solidification involving diffusion in the solid[J].Physical Review E,2009,79:031603.
[29] VIARDINN A, SOUHAR Y, FERN魣NDEZ M C, ALEL M, ZALO譕NIK M.Mesoscopic modeling of equiaxed and columnar solidification microstructures under forced flow and buoyancy-driven flow in hypergravity:Envelope versus phase-field model[J].Acta Materialia,2020,199:680-694.
[30] MISHIN Y, HERZIG C.Diffusion in the Ti-Al system[J].Acta Materialia,2000,48(3):589-623.
[31] 介万奇，傅恒志，周尧和.多元多相合金凝固理论模型的研究进展[J].中国材料进展，2010,29(6):1-11.JIE W Q,FU H Z,ZHOU Y H.Progress of solidification researches on the multi-component multi-phase alloys[J].Materials China,2010,29(6):1-11.
[32] OHNO M,MATSUURA K.Motion and morphology of triple junction in peritectic reaction analyzed by quantitative phase-field model[J].ISIJ International,2010,50(12):1879-1885.