固/液界面稳定性是晶体生长和金属凝固过程中的重要问题[1]，直接决定了实际材料和铸件的凝固组织特征和性能，是物理学家和材料学家关注的焦点之一[2]。20 世纪60 年代中期，Mullins 和Sekerka[3]将扰动振幅与时间的关系引入到平界面稳定性的判断中，借鉴流体运动稳定性分析方法，提出了界面稳定性理论(M-S 理论)，为研究固/液界面形态特征提供了重要依据。然而，M-S 理论是在假设界面前沿溶质扩散场稳态的基础上建立，而在实际凝固过程中，其界面前沿的溶质扩散不能迅速达到稳态，对于平界面失稳后界面形态的进一步演化问题还需要建立平界面失稳的非线性动力学理论[4]。 基于此，Wollkin 等[5]给出了在一维方向扰动振幅很小的界面失稳的非线性微分方程式。 随后Langer 和Turski[6]发现非平面界面与流体力学中对流模式具有相似性，进一步建立了非线性不稳定的对称理论模型，但模型中无非线性项的具体控制方程，只能定性分析非线性项对固/液界面的影响情况。 基于此，国内学者[7]构建了一种耗散结构的控制模型，并在该模型中引入了非线性项，成功地定量分析了扰动作用下固/液界面形态的演化过程。 然而，这些研究结果均是假设合金相图线性化，即在合金相图特征参数为常数下所得。 Garandet 等[8]和刘少军等[9]分析了Sn-Bi 合金和Al-Zn 二元合金平界面稳定性，结果表明，采用相图线性化参数导致模型预测结果与实验结果之间存在较大误差，甚至有些误差高达10%。
随着热力学理论和计算机技术的发展，使用热力学理论来获取合金相图及其相关参数的方法，如CALPHAD 等逐渐发展起来，并应用于材料的设计和制备，使得理论模型与实验结果更加接近[10]。Chen[11]和Reuther[12]等耦合热力学理论计算分析了SCN-acetone 和Si-Sn 二元合金的平面失稳参数， 发现与实验吻合较好。 刘少军等[9]对比了线性化参数与热力学参数下Al-0.083Zn 二元合金的界面失稳情况，发现热力学参数下模型结果与实验吻合更好。 然而，在实际应用中，商业合金通常为多元合金，溶质之间存在强烈的相互作用， 极大地影响了合金的凝固行为。许多学者[13-15]尝试将凝固界面稳定性理论扩展到多元合金体系中， 但大部分模型仅是将二元合金的参数直接应用于多元体系， 这导致了理论预测与实验结果之间的较大偏差[16]。 基于此，本文利用耦合热力学计算技术，以Al-0.3%Si-0.1%Cu(质量分数)三元合金为例，确定出热力学平衡条件下的相图特征相关参数，讨论多元合金凝固过程中的凝固界面线性稳定性、非线性稳定性和平胞转变临界分叉等问题， 分析线性化参数和热力学计算所得三元参数的适用范围。
1 固/液界面稳定性模型
在合金凝固过程中，凝固组织形貌随生长速度增大，固/液界面将经历平面晶→胞晶→枝晶→细胞晶→准平面。胞晶、枝晶、细胞晶、准平面晶都是平界面失稳后形成的新的有序结构，分别对应着一种非平衡定态，由非平衡非线性动力学方程控制[4]：
式中，δ 为固/液界面前沿的扰动振幅；t 为时间；a0 为线性稳定性系数；a1，a2，…，an 为非线性项系数，其中a1 为Landau 常数。
1.1 线性稳定性模型
假设式(1)中非线性项系数均为0，仅分析线性系数a0 对凝固固/液界面稳定性的影响作用[3]。 联合凝固过程中温度场、溶质场、固/液界面凝固速率，可获得固/液界面前沿存在正弦波动扰动条件下，固/液界面扰动振幅随时间的变化率：
式中
为正弦波振幅对时间的导数；ω 为振动频率；ν 为固/液界面凝固速率；Tm 为合金熔点；DL,j 为溶质i 在液相中的扩散系数；κs、κL 分别为固相、液相的热导率；GL、Gs 分别为固相、 液相中的温度梯度；mi 为溶质i 的液相线斜率
为液相中溶质i 的浓度梯度。
据此， 得到判断多元合金平界面稳定性的界面稳定性动力学判别式[7]：
当S(ω)<0 时(即δ˙/δ<0 或δ˙<0)，平面凝固界面是稳定的；反之，平面凝固界面将失稳。
1.2 固/液界面非线性稳定分析
忽略式(1)高次非线性项，仅分析线性系数a0 及非线性系数a1 的综合作用对固/液界面稳定性的影响。 考虑到非线性项a1 的存在，综合分析固/液界面前沿浓度场、温度场及曲率的非线性项，按照推导线性稳定性方程中δ˙<0 的方法，可得固/液界面非线性稳定性的控制方程为[6]：
当δ2ω2＞1 时，
当δ2ω2＜1时，
当dδ/dt>0 时，凝固固/液界面不稳定；反之，凝固固/液界面稳定。 将式(4)和(5)与式(3)对比，可以发现，固/液界面非线性稳定方程是线性稳定方程的扩展，固/液界面的稳定性将由式(4)和(5)中的线性项和非线性项共同决定。
2 三元合金热力学模型
根据Redlich-Kister-Muggianu 模型， 三元合金的液相和FCC 固溶体相Gibbs 自由能可以表示为以下形式[17]：
式中
为φ 相的摩尔吉布斯自由能；ci 为组元i 的摩尔分数
为纯组元i 在参考态下的摩尔自由能
为二元系i-j 的第(n+1)个Redlich-Kiste 系数
为三元系的三元相互作用系数。
根据热力学原理， 体系在恒温恒压条件下达到平衡的一般条件为：
联立式(6～11)，即可获取三元合金的相图特征参数，如三元系的溶质分凝系数ki 及液相面斜率mi 等。
3 界面稳定性计算结果与讨论
3.1 Al-Si-Cu 相图特征参数
Al-Si-Cu 系合金，如A380、ADC2、A390 等铝合金，具有良好的铸造性能，几乎不产生铸造裂纹，热胀系数小，耐磨性好，是当前重要的高强度铸造合金之一，也是材料学者研究的热点合金之一，被广泛应用于航空航天和汽车工业。 同时，Al-Si-Cu 合金富Al 角相图比较简单，有比较准确的物性参数，在许多理论研究中广泛采用。 该合金优化后的热力学模型参数见表1[18]。
表1 Al-Si-Cu 合金热力学参数[18]
Tab.1 Thermodynamic parameters of the Al-Si-Cu ternary alloy[18]
?
图1a 是使用表1 热力学模型参数计算得到的Al-0.3%Si-0.1%Cu 合金溶质分凝系数。 可以看出，Al-0.3%Si-0.1%Cu 合金的溶质Si 分凝系数随着温度的降低而升高，在913～931 ℃范围内变化约3%，溶质Cu 的分凝系数则随着温度的降低而降低，变化1.7%， 显然在该合金中溶质Si 的偏析将更严重。 图1b 是使用表1 热力学模型参数计算所得Al-0.3%Si-0.1%Cu 合金的液相线斜率，液相线斜率随着温度的增大而增大。 表2 给出了计算过程中Al-Si-Cu 合金的其余物理参数， 同时给出了线性化溶质分凝系数及液相面(线)斜率参数[19-20]。 研究发现[21]，相图特征参数，如溶质分凝系数、液相线斜率使用常数来定量描述。 在低浓度时是合理的， 在高浓度和多元合金中，由于这些参数不仅与温度，还和其他组元的种类和浓度密切相关，且和熔体结构、点阵位置和溶质相互作用有关，是变量，使用热力学计算相图特征参数具有更好的物理意义。
图1 Al-0.3%Si-0.1%Cu 合金Si、Cu 溶质分凝系数和液相线斜率随温度的变化曲线：(a)溶质分凝系数；(b)液相线斜率
Fig.1 Partition coefficients and slopes of the liquidus lines of Si and Cu during solidification of Al-0.3 wt.%Si-0.1 wt.%Cu：(a)partition coefficient;(b)slope of the liquidus line
表2 Al-Si-Cu 合金计算过程中采用的物理参数[19-20]
Tab.2 Physical parameters of Al-Si-Cu used in the numerical calculation[19-20]
?
3.2 线性稳定性计算结果
图2a 是温度梯度G=1 000 K/m 条件下，Al-0.3%Si-0.1%Cu 合金在不同生长速度ν 和扰动频率ω 作用下， 分别采用线性化参数和热力学计算所得固/液界面稳定区域和失稳区域图。在ν 和ω 所形成的曲线内部固/液界面是不稳定的，而且不稳定区域存在于2 个临界速度νc 和νa 之间， 其中νc 是平面晶向胞状晶转变的临界生长速度，νa 是平界面绝对稳定性生长速度，超出这个范围的所有ω 值所对应的S(ω)值均小于0，平界面将保持稳定生长。 图2b是温度梯度G=1 000 K/m， 不同固/液界面生长速度条件下，分别使用热力学计算技术计算参数和线性化参数， 得到Al-0.3%Si-0.1%Cu 合金固/液界面线性稳定性判据S(ω)和扰动频率ω 的关系曲线。 实线为采用热力学计算所得，虚线为采用线性化参数所得。S(ω)>0 的区域为固/液界面失稳区域，S(ω)<0 的区域为固/液界面稳定区域。 可以看出，两种参数条件下计算所得固/液界面稳定区域有很大差别。 当固/液界面生长速度较小时，如ν=10 μm/s 条件下，在整个扰动频率范围内，两种参数计算所得均显示凝固固/液界面处于稳定状态； 而当固/液界面生长速度较大时，如ν=30 μm/s 条件下，两种参数所得固/液界面随ω 的增加都经历固/液界面稳定到固/液界面失稳再到固/液界面稳定的变化规律。 当ν=20 μm/s 时，两种参数计算所得固/液界面线性稳定性结果却截然不同：采用线性化参数计算结果表明，在整个扰动频率范围内凝固固/液界面均处于稳定状态；而采用热力学计算参数，计算所得凝固固/液界面将随着扰动频率的增加经历固/液界面稳定、固/液界面失稳、固/液界面再次稳定3 个阶段。 可见，这与相图特征参数的选取有极大的关系。
图2 Al-0.3%Si-0.1%Cu 合金固/液界面线性稳定性区域图(G=1 000 K/m)：(a)ν 与ω 形成的稳定区域图；(b)界面稳定参数S(ω)随ω 的变化曲线
Fig.2 Stability map of the planar front predicted for the Al-0.3 wt.%Si-0.1 wt.%Cu alloy(G=1 000 K/m)：(a)stable region formed by ν and ω;(b)variation in S(ω)with ω at various growth velocities
图2b 中同时标出了S(ω)=0 时，固/液界面线性稳定性临界扰动频率ω 值ω-和ω+。 可以看出，采用热力学计算所得固/液界面稳定区域明显比采用线性化参数计算所得区域小。 同时，当ω 值很大或很小时，两种参数计算所得线性稳定性判据S(ω)值相近，且均处于固/液界面稳定性区域。 究其原因在于当ω 很大时，式(3)中Гω2 项值很大，即固/液界面能约束固/液界面起伏，维护固/液界面的稳定。 而当ω很小， 固/液界面及溶质梯度成为微不足道的因素，此时固/液界面稳定性取决于固/液界面处温度梯度[4,9]。在这两种条件下， 无论使用线性化参数还是热力学计算参数， 对凝固固/液界面稳定性的研究影响不大，固/液界面均处于稳定状态。
可以看出， 无论是对临界生长速度还是对扰动频率对界面稳定的影响，采用热力学计算所得固/液界面稳定区域明显比采用线性化参数所得区域小，对界面的失稳的判定更为严苛。
3.3 非线性稳定性计算结果
图3为Al-0.3%Si-0.1%Cu 合金在温度梯度G=1 000 K/m，ν=30 μm/s，δ2ω2>1 条件下，采用非线性稳定性判据式(4)计算所得凝固固/液界面的稳定区域和失稳区域图。实线为采用热力学计算所得，虚线为采用线性化参数所得， 垂直线为非线性判据中线性项为零时所对应的扰动频率值ω0。 可以看出，当扰动ω 在图中0～ω0 区域时， 即非线性判据中线性项小于零， 此时非线性稳定性判据dδ/dt 始终小于零，固/液界面扰动将随着时间的增大而逐渐减小，固/液界面将始终保持稳定。 当扰动ω 大于ω0 时，非线性判据中线性项将大于零， 由于非线性项始终小于零， 此时必然会在ω0 的右侧产生一个固/液界面稳定区域， 即图3 局部放大图的三角区域及扰动频率ω 达到ωCritical 时图中所示稳定区域。 在这两个区域，固/液界面首先发生失稳，并在一定时间内形成新的有序结构，如胞状结构或枝晶结构，随后根据耗散理论[4]，当胞晶或枝晶结构被消耗后，固/液界面又重新恢复到原来的平界面，该区域被称为亚临界分叉区[7]。
图3 Al-0.3%Si-0.1%Cu 合金在G=1 000 K/m，ν=30 μm/s 下固/液界面非线性稳定曲线(δ2ω2>1)
Fig.3 Non-linear instability curves of Al-0.3 wt.%Si-0.1 wt.%Cu at G=1 000 K/m-1 and ν=30 μm/s when δ2ω2>1
图4为Al-0.3%Si-0.1%Cu 合金在温度梯度G=1 000 K/m，ν=30 μm/s，δ2ω2<1 条件下， 采用非线性稳定性判据式(5)计算所得凝固固/液界面的稳定区域。实线为采用热力学计算所得，虚线为采用线性化参数所得， 垂直线为非线性判据中线性项为零时所对应的扰动频率值。 不同于δ2ω2>1 的情况，δ2ω2<1 条件下其非线性判据中线性项增加了曲率对固/液界面稳定性的作用， 从而导致固/液界面的稳定区域因扰动频率不同而不同，具体结果如图4a 和b 所示。 当线性项大于零时，非线性项恒为正，则固/液界面将失稳，形成新的有序结构。而当线性项小于零时，由于非线性项恒为正，则固/液界面很有可能出现失稳区域，如图4a 和b 中三角区域，这些区域内固/液界面会失稳，形成新的有序结构，这些区域也被称为亚临界分叉区[7]。
图4 Al-0.3%Si-0.1%Cu 在G=1 000 K/m、ν=30 μm/s下固/液界面非线性稳定曲线(δ2ω2<1)：(a)ω<8×103 m-1;(b)ω>3×104 m-1
Fig.4 Nonlinear instability curves of Al-0.3 wt.%Si-0.1 wt.%Cu at G=1 000 K/m and ν=30 μm/s when δ2ω2<1：(a)ω<8×103 m-1;(b)ω>3×104 m-1
对比图4 和图2b 可以发现， 图4 中扰动频率ω-、ω+ 值与图2b 中ν=30 μm/s 条件下固/液界面稳定性判据S(ω)=0 时，所对应的固/液界面线性稳定性临界扰动频率ω 值ω-和ω+相等。不同点在于，图2 中ω-左侧、ω+右侧的区域， 线性稳定性判据计算所得结果表明固/液界面均处于稳定状态。 而在图4中，ω-左侧、ω+右侧区域，却显示出不同大小程度的固/液界面失稳区域，即亚临界分叉区。 主要原因在于ω-左侧扰动频率越来越小，ω+右侧扰动频率越来越大时，非线性项逐渐发挥了主导作用，非平衡方程就有对应的定态解，有些定态解是稳定的，有些是不稳定的，从而导致最终固/液界面稳定或是失稳[4]。可以认为， 非线性稳定性判据是对线性稳定性临界扰动频率ω-和ω+值两侧稳定区域的细化。
4 结论
(1)基于固/液界面线性稳定性和非线性稳定判据，获得三元合金固/液界面稳定区域和失稳区域及亚临界分叉区域范围。
(2)采用热力学计算所得稳定区域较线性化参数所得稳定区域小，判断固/液界面稳定性更严苛。
(3)扰动频率很大或者很小时，热力学计算所得参数和线性化参数可以互相替代，对于最终固/液界面稳定性判断影响不大。
[1] TOURRETD,KLEMM-TOOLEJ,CASTELLANOSAE,RODGERS B,BECKER G,SAVILLE A,ELLYSON B,JOHNSON C,MILLIGAN B, COPLEY J, OCHOA R, POLONSKY A, PUSCH K,HAINES M, FEZZAA K, SUN T, CLARKE K, BABU S, POLLOCK T, KARMA A, CLARKE A. Morphological stability of solid-liquid interfaces under additive manufacturing conditions[J].Acta Materialia,2023,250：118858.
[2] 孙思杰， 蒋晗. 各向异性界面动力学对深胞晶生长形态稳定性的影响[J].物理学报，2024,73(11)：331-341.SUN S J, JIANG H. Effects of anisotropic interfacial kinetics on morphology stability of deep cellular crystal growth[J].Acta Physica Sinica,2024,73(11)：331-341.
[3] MULLINS W W, SEKERKA R F. Stability of a planar interface during solidification of a dilute binary alloy[M]//PELCÉ P.Dynamics of Curved Fronts.Cambridge：Academic Press,1988：345-352.
[4] 徐鉴君. 凝固过程动力学与交界面稳定性理论导引[M]. 北京： 科学出版社，2006.XU J J.Introduction to kinetics of solidification and stability theory of the interface[M].Beijing：Academic Press,2006.
[5] WOLLKIND D, RAISSI S. A nonlinear stability analysis of the freezing of a dilute binary alloy[J].Philosophical Transactions of the Royal Society of London.Series A：Mathematical and Physical Sciences,1970,268：351-380.
[6] LANGER J S, TURSKI L A. Studies in the theory of interfacial stability—I. Stationary symmetric model[J]. Acta Metallurgica,1977,25(10)：1113-1119.
[7] 胡汉起.金属凝固原理[M].北京：机械工业出版社，2021.HU H Q.The Solidification theory of metals[M]. Beijing： China Machine Press,2021.
[8] GARANDET J P, BOUTET G, LEHMANN P, DREVET B,CAMEL D,ROUZAUD A,FAVIER J J,FAIVRE G,CORIELL S,ALEXANDER J I D, BILLIA B. Morphological stability of a solid-liquid interface and cellular growth： Insights from thermoelectric measurements in microgravity experiments[J].Journal of Crystal Growth,2005,279(1-2)：195-205.
[9] 刘少军，杨光昱，张胜利，介万奇. 二元稀溶液合金凝固界面稳定性分析[J].特种铸造及有色合金，2011,31(3)：278-281,195.LIU S J,YANG G Y,ZHANG S L,JIE W Q.A interface stability analysis of a dilute binary alloy[J].Special Casting&Nonferrous Alloys,2011,31(3)：278-281,195.
[10] GUILLEMOTG,SENNINGERO,HARELANDCA,VOORHEES P W,GANDIN C A.Thermodynamic coupling in the computation of dendrite growth kinetics for multicomponent alloys[J]. Calphad,2022,77：102429.
[11] CHEN M W, WANG Z D. The evolution and morphological stability of a particle in a binary alloy melt[J]. Journal of Crystal Growth,2023,607：127113.
[12] REUTHER K,GALENKO P K,RETTENMAYR M. Dynamic instability of the steady state of a planar front during non-equilibrium solidification of binary alloys[J].Journal of Crystal Growth,2019,506：55-60.
[13] LUDWIG A,PUSTAL B,HERLACH D M.General concept for a stability analysis of a planar in terface under rapid solidification conditions in multi-component alloy systems[J].Materials Science and Engineering：A,2001,304-306：277-280.
[14] ALTIERI A L, DAVIS S T. Instabilities in solidification of multi-component alloys[J]. Journal of Crystal Growth, 2017, 467：162-171.
[15] ALTIERI A L, DAVIS S T. Instabilities in rapid solidification of multi-component alloys[J].Journal of Crystal Growth,2017,476：78-89.
[16] 张瑞杰，王晓颖，陈忠伟，何志，介万奇. 耦合热力学计算研究多元合金的凝固界面稳定性[J].铸造技术，2004(2)：117-119.ZHANG R J,WANG X Y,CHEN Z W,HE Z,JIE W Q.Research on the solid-liquid interface stability of multicomponent alloys with CALPHAD method[J].Foundry Technology,2004(2)：117-119.
[17] SUNDMAN B,JANSSON B,ANDERSSON J O.The thermo-calc databank system[J].Calphad,1985,9(2)：153-190.
[18] HE C Y,DU Y,CHEN H L, XU H H. Experimental investigation and thermodynamic modeling of the Al-Cu-Si system[J].Calphad,2009,33(1)：200-210.
[19] REBOW M,BROWNE D J.On the dendritic tip stability parameter for aluminium alloy solidification[J].Scripta Materialia,2007,56(6)：481-484.
[20] KURZ W, FISHER D J, RAPPAZ M. Fundamentals of solidification[M].Switzerland：Trans Tech Publications Ltd.2023.
[21] GUI L T,LONG M J,CHEN D F,ZHAO J J,WANG Q Z,DUAN H M.Thermodynamic study on the solute partition coefficients on L/δ and L/δ+γ phase interfaces for 1215 high-sulfur steel solidification by orthogonal design[J].Journal of Materials Research and Technology,2020,9(1)：89-103.