摘要：目的 针对高温合金叶片在定向凝固过程中容易出现雀斑缺陷，从而导致叶片报废的问题，对定向凝固枝晶生长与溶质对流进行模拟研究，以揭示雀斑缺陷的形成规律。方法 针对CM247LC合金定向凝固过程，采用相场模型模拟凝固过程枝晶生长，采用格子Boltzmann模型模拟溶质浓度差引起的自然对流。采用基于双重网格的GPU并行算法对相场-格子Boltzmann模型进行数值求解。研究在不同晶体取向角度与取向差条件下的枝晶形貌、对流速度及溶质羽流的演变规律。结果 当晶体取向角度不同时，在枝晶生长过程中，液相区域的平均对流速度均表现为周期性变化。当晶体取向角度较大时，随着晶体取向角度的变大，一次枝晶臂间距变大。当枝晶间存在晶体取向差时，溶质羽流倾向于在发散型晶界附近发起；随着晶体取向差的增大，溶质羽流发起时间提前。溶质羽流的形成阻碍了枝晶尖端及附近枝晶侧臂的生长。结论 晶体取向角度对溶质羽流形成的影响较小，较大的晶体取向差对溶质羽流的形成有促进作用。
关键词：高温合金；定向凝固；枝晶生长；溶质对流；相场模拟
高温合金定向凝固技术是制造航空发动机和燃气轮机涡轮叶片的主要成形工艺。高温合金叶片包括等轴晶叶片、定向凝固柱状晶叶片及单晶叶片。其中，定向凝固柱状晶叶片与单晶叶片消除了横向晶界，具有优异的高温力学性能。由于高温合金叶片的结构较为复杂，其凝固过程中的温度场难以稳定控制，所以在定向凝固过程中容易出现雀斑缺陷。定向凝固雀斑缺陷的形成与凝固过程中枝晶组织与溶质对流之间的相互作用有关。
目前，主要采用实验与数值模拟的方法对雀斑缺陷的形成机制与规律进行研究。Pollock等[1]研究了高温合金定向凝固形成的雀斑缺陷，研究发现，当冷却速率较低时，一次枝晶臂间距较大，雀斑缺陷更容易形成，且不同高温合金之间的雀斑缺陷形成倾向存在显著差异。Tin等[2]研究发现，将碳的质量分数增大至0.1%能够显著降低雀斑缺陷数量，这得益于富Ta的MC型碳化物析出。Tin等[3]采用回归分析的方法得到了高温合金试样中雀斑缺陷数量与糊状区密度差及一次枝晶臂间距的关系式，该关系式可用于高温合金雀斑缺陷形成倾向的评估。Ma等[4]研究了试样几何形状对雀斑缺陷的影响，发现在试样截面扩张与收缩位置容易形成雀斑缺陷。Shevchenko等[5]采用同步辐射技术观察了Ga-In合金定向凝固过程中溶质羽流的形成，结果表明，富集的溶质促进了偏析通道的形成。Reinhart等[6]采用同步辐射技术观察了CMSX-4合金凝固过程中枝晶生长与溶质羽流的相互作用。
研究者采用数值模拟的方法研究了定向凝固枝晶生长与溶质对流现象，以揭示雀斑缺陷的形成机制。Schneider等[7]基于高温合金热力学数据库对定向凝固过程中溶质羽流的发起进行了模拟，结果表明，随着冷却速率的减小，凝固前沿的溶质分布转变为不稳定状态，进而出现了通道偏析。Felicelli等[8]对Pb-10%（质量分数）Sn合金定向凝固过程中的溶质对流现象进行了三维模拟，模拟结果表明，在壁面位置出现了间隔一定距离的偏析通道。Yuan等[9]采用元胞自动机模型研究了Pb-Sn合金定向凝固过程中的溶质对流现象，结果表明，由对流引起的枝晶间富集的溶质使枝晶干和二次枝晶重熔，进而形成了通道偏析。Karagadde等[10]采用元胞自动机模型对Ga-25%（质量分数）In合金定向凝固过程进行了模拟，发现溶质羽流的形成与枝晶取向及一次枝晶臂间距密切相关，且形成溶质羽流的临界瑞利数为150～170。Kao等[11]研究了固液界面形状对定向凝固枝晶生长与溶质对流的影响，模拟结果表明，当固液界面为凸界面时，溶质羽流倾向于在中间位置发起，而当固液界面为凹界面时，则在两侧壁面位置发起。由于溶质对流过程中的流动场求解计算量较大，开展大规模模拟需要极大的计算量。为提高模拟速度，Sakane等[12]提出了基于多个GPU的并行求解算法，采用格子Boltzmann模型对溶质差引起的浮力驱动流动过程进行了快速求解。Guo等[13-14]采用自适应网格算法对热溶质对流条件下的枝晶生长过程进行了并行求解。Yang等[15]提出了基于耦合热力学数据库的GPU并行算法，并对CMSX-4合金凝固过程中的枝晶生长与溶质对流现象进行了模拟。Takaki等[16]采用二维相场模型研究了自然对流条件下Al-Cu合金枝晶生长过程中的竞争生长规律。Takaki等[17]对Al-Cu合金自然对流条件下的枝晶生长过程进行了大规模相场模拟，得到了不同重力加速度条件下一次枝晶臂随时间的变化规律。Sakane等[18]采用多套网格系统对枝晶生长过程中的相场、流动场与温度场模型进行求解，有效提高了模拟计算速度。Watanabe等[19]提出了一种用于多相场模型的自适应结构化网格区域分解算法，以提高计算速度。此外，也有模型考虑了凝固过程中的枝晶运动[20]、碰撞[21]及枝晶与气孔的相互作用[22-23]等。
目前，晶体取向角度与取向差对高温合金定向凝固枝晶生长与溶质对流的影响规律尚不清晰，本文采用相场-格子Boltzmann模型对高温合金定向凝固过程进行了枝晶尺度模拟，得到了不同晶体取向角度与取向差时溶质羽流的形成倾向以及溶质对流速度的变化规律。
1 定向凝固过程中枝晶生长数学模型
在定向凝固过程中，枝晶生长的数学模型包括相场模型与格子Boltzmann模型。其中，相场模型用于求解凝固过程中相界面的演化与溶质浓度分布，格子Boltzmann模型用于求解溶质差驱动的流动场。相场模型的控制方程如式（1）所示。
式中：为序参量，在固相中，
=+1，在液相中，
=−1；W0为界面宽度；τ0为松弛时间；k为溶质分配系数；t为时间；
，为各向异性函数的表达式，其中，ε为固液界面能各向异性系数，
和
分别为序参量在x和y方向的偏导；λ为耦合系数；U为过饱和度；θ为无量纲温度。U和θ的表达式分别如式（2）和式（3）所示。在定向凝固条件下，存在沿着抽拉方向的温度梯度。为此，本文采用“冻结温度”方法设置整个计算域的温度分布，即T(x,y,t)=T0+G(y–vpt)，其中G为抽拉方向的温度梯度，vp为定向凝固抽拉速度，x和y分别为二维区域的横、纵坐标，T0为计算区域底部（对应y=0位置）的初始温度。
(2)
式中：TM为溶剂元素熔点；m为液相线斜率；c为溶质浓度；为平衡时固液界面液相一侧的溶质浓度。模型考虑了溶质对流，溶质浓度的控制方程如式（4）所示。
(4)
式中：为速度矢量；D为溶质在液相中的扩散系数；
为插值函数，
；
为反溶质截留项的系数，
=
。
采用格子Boltzmann模型对流动场进行求解，选取Bhatnagar-Gross-Krook（BGK）碰撞算子形式的格子Boltzmann方程，如式（5）所示。
式中：x为网格所在的空间位置矢量；fi为i方向的粒子分布函数；fieq为i方向的平衡态粒子分布函数；Δt为时间步长；τ为松弛时间；Fi为外力项；ci为i方向的离散格子速度。平衡态分布函数表达式如式（6）所示。
(6)
式中：ρ为流体密度；u为流体速度矢量；cs为格子声速；wi为权系数。这里采用D2Q9格子速度模型，wi的取值如式（7）所示。
采用Boussinesq假设来处理溶质差引起的浮力F，其计算式如式（8）所示。
(8)
式中：g为重力加速度矢量；fs为单元的固相分数，满足；βc为溶质膨胀系数。由此可得外力Fi的表达式如式（9）所示。
采用有限体积法对相场方程和溶质场的对流扩散方程进行求解。采用一阶Euler格式对时间导数进行离散，采用迎风格式离散对流项，以文献[24]中的通量形式对式（4）的右端项进行处理。为减小计算量，分别在粗、细2套网格上求解格子Boltzmann方程和相场方程。参考文献[25]，粗细网格的边长比值为2，采用GPU并行算法进行求解，粗、细网格上的场变量需要进行相应的插值。
2 定向凝固过程中枝晶生长的计算参数
选取CM247LC合金进行定向凝固枝晶生长模拟，合金成分如表1所示。采用伪二元近似方法[26]得到相场模型所需的液相线斜率和平衡分配系数等参数。相场的计算网格尺寸为1 μm，流动场的计算网格尺寸为2 μm，界面宽度W0为1.25 μm，时间步长为2.0×10−5 s，界面能各向异性系数为0.02，液相溶质扩散系数为3.6×10–9 m2/s，液相的运动黏度为1.0× 10–6 m2/s，重力加速度常量为−9.81 m/s2，耦合系数λ=33.27。前期研究表明[25]，该耦合系数的取值能够保证相场模拟结果具备收敛性。当耦合系数采用该取值时，CMSX-4合金定向凝固枝晶生长的模拟结果表明，在溶质羽流发起位置的枝晶生长速度呈周期性振荡，与同步辐射实验[6]观测到的振荡周期与幅值均具有较好的一致性[25]。因此，本文计算模型在模拟枝晶生长与溶质对流相互作用方面具有较高的准确性。
表1 CM247LC合金成分
Tab.1 Composition of CM247LC alloy wt.%
模拟的工艺参数如下：G=10 K/mm，vp=50 μm/s。前期的计算结果表明[25]，在该工艺条件下，CM247LC合金中溶质羽流的发起处于临界稳定条件。在该工艺条件下，有利于溶质羽流的形成，从而研究溶质对流条件下晶体取向角度与取向差的影响规律。在不同晶体取向角度的模拟过程中，在整个计算区域底部设置一层连续的固相，取向角度为0°～30°，取值间隔为5°。计算区域大小为2 048 μm× 4 096 μm。在不同晶体取向差的模拟过程中，计算区域初始条件与边界条件设置如图1所示。在底部区域设置2种不同取向的晶粒，计算区域的左右边界均为周期性边界条件，底部为等间距放置的初始晶粒，晶粒之间的间距Λ0为稳态时的一次枝晶间距，具体数值由前期的模拟结果得到[25]。在底部左右各设置4个初始晶粒，取向角度为0°，在中间设置8个初始晶粒，取向角度Θ分别为5°、10°和15°，计算区域大小为4 096 μm×4 096 μm，采用移动计算区域的方法以减小计算量。当抽拉距离为10 mm时，终止计算。模拟计算在单个NVIDIA V100 GPU上进行。
3 结果与分析
当晶体取向角度为5°时，溶质浓度分布与对流速度分布模拟结果如图2所示，图2中的C为无量纲溶质浓度。可以看到，在凝固初期（40 s时），初始条件的平界面发展为树枝晶，且由于竞争生长，一次枝晶臂间距逐渐变大。当凝固时间为100 s时，凝固前沿的溶质在对流作用下富集，并形成溶质羽流，如图2b中左起第2、3个枝晶间。可以看出，该位置的枝晶臂间距较大。当枝晶臂间距较大时，溶质对流的阻力较小，因而更容易形成溶质羽流。随着凝固过程的进行（150 s和200 s），形成的溶质羽流由于具有较低的密度而上浮，在上浮过程中，羽流中心的溶质浓度在对流扩散作用下有所降低，因而羽流上浮的速度有所减缓，开始在环流式对流的作用下向两侧运动，使得溶质羽流呈现为“烟囱”状。模拟结果表明，最大对流速度可达160 μm/s，远大于该工艺参数条件下枝晶生长的稳态速度50 μm/s。
图1 计算区域初始条件与边界条件
Fig.1 Initial conditions and boundary conditions of the computational domain
图2 晶体取向角度为5 °时的模拟结果
Fig.2 Simulation results of the case with crystal orientation angle of 5°
晶体取向角度为15°时和30°时的模拟结果分别如图3和图4所示。可以看到，模拟结果中均出现了溶质羽流。当晶体取向角度为30°时，一次枝晶臂间距大于取向角度为15°时的一次枝晶臂间距，且在凝固200 s时，晶体取向角度为30°时的对流速度最大可达180 μm/s，略高于15°时的150 μm/s。
统计了不同晶体取向角度时的平均对流速度及其时间平均值，结果如图5所示。平均对流速度定义为，其中，
为(x,y)位置在t时刻的流体速度矢量，dS为空间微元，S为液相区域，S0为液相区域的面积。vtavg可以反映不同时刻溶质对流的强度。相应地，对平均对流速度在不同时刻求平均值，将其定义为
，其中，t1、t2分别为计算时间平均值的开始和终止时刻。
图3 晶体取向角度为15°时的模拟结果
Fig.3 Simulation results of the case with a crystal orientation angle of 15°
图4 晶体取向角度为30 °时的模拟结果
Fig.4 Simulation results of the case with a crystal orientation angle of 30°
图5 不同晶体取向角度时的平均对流速度vtavg（a）及其时间平均值vavg（b）
Fig.5 (a) Average convection velocity vtavg and (b) corresponding time-averaged value vavg under different crystal orientation angles
统计结果表明，当晶体取向角度从0°变化到30°时，在凝固初期（凝固时间小于50 s左右时），平均对流速度均保持在较小值（低于10 μm/s）。随着凝固过程的进行，平均对流速度迅速增大至峰值，速度峰值超过70 μm/s。之后，平均对流速度开始衰减，并保持周期性振荡。除了晶体取向角度为15°时平均对流速度增大至峰值较为滞后（在120 s左右速度出现明显的增大）外，当晶体取向角度为其他取值时，平均对流速度开始迅速增大的时间均较早（在60～80 s之间）。结果表明，当晶体取向角度不同时，平均对流速度时间平均值均大于50 μm/s，即大于抽拉速度，即对应枝晶生长速度的稳态值。平均对流速度振荡幅值较大（见图5a），其时间平均值的分散度较大，这表明，平均对流速度时间平均值与晶体取向角度之间没有明显的依赖性。
溶质羽流的发起与一次枝晶臂间距密切相关。为此，统计了凝固终态时不同晶体取向角度对应的一次枝晶臂间距的最大值、最小值及均值，结果如图6所示。
图6 不同晶体取向角度时的一次枝晶臂间距的最大值、最小值及均值
Fig.6 Maximum, minimum and average values of primary dendrite arm spacing under different crystal orientation angles
当晶体取向角度为0°～20°时，一次枝晶臂间距均值Λavg基本保持在200 μm，且当晶体取向角度为0 °时，一次枝晶臂间距均值Λavg最小。当晶体取向角度为20°～30°时，随着晶体取向角度的变大，一次枝晶臂间距均值Λavg变大。Gandin等[27]给出了一次枝晶臂间距Λ与晶体取向角度的幂函数关系：，其中，ΔT0为合金的凝固温度范围，Θ为晶体取向角度，a、b、c、d和e为拟合系数，且拟合系数d、e均为正数。由拟合可得，d=0.08，e=12。基于该关系的预测结果表明，一次枝晶臂间距均值Λavg随着晶体取向角度的变大而变大。而本文的模拟结果表明，当晶体取向角度较大时，一次枝晶臂间距随晶体取向角度的变化符合上述关系；而当晶体取向角度较小（低于20°）时，一次枝晶臂间距随晶体取向角度的变化不显著，与上述关系存在一定的偏差。由于Gandin等[27]是在不考虑溶质对流条件下得到的一次枝晶臂间距预测模型。因此，这种偏差表明，当枝晶间存在溶质对流时，强烈的对流会改变凝固前沿的溶质分布，使不同位置的枝晶生长速度存在显著差异，进而影响一次枝晶臂间距的发展。
晶体取向差为5°时的溶质分布与速度分布如图7所示。结果表明，当晶体取向差为5°时，在右侧发散型晶界处，对流速度较大。当凝固时间为74 s时，在发散型晶界两侧出现明显的溶质浓度起伏（见图7a）。随着凝固过程的进行，在79 s时，发散型晶界右侧率先形成溶质羽流（见图7b）。随后，溶质羽流逐渐向上运动（见图7c），之后，溶质羽流从两侧向下方回流（见图7d），溶质羽流中心的浓度减小。在溶质羽流发起的位置，倾斜枝晶生长出的二次枝晶及其分枝的生长被抑制，形成的枝晶间隙较大，进而有利于溶质羽流的发展。取向差为10°和15°时的模拟结果也呈现出类似的规律。
图7 不同凝固时间下晶体取向差为5°时的模拟结果
Fig.7 Simulation results of the case with a crystal orientation difference of 5° for different solidification time
统计了不同晶体取向差时的平均对流速度，如图8所示。可以看到，当晶体取向差不同时，平均对流速度均在增大至峰值后，在一定范围内呈周期性振荡，但达到峰值的时间不同。随着晶体取向差变大，平均对流速度达到峰值的时间逐渐减小。从图8可以看出，当晶体取向差为5°、10°、15°时，达到峰值的时间（如图8中竖直黑色虚线所示）分别为98、84、63 s。平均对流速度达到峰值的时间越小，表明对流过程中溶质羽流出现得越早。因此，模拟结果表明，晶体取向差越大，溶质羽流越容易发起。
图8 不同晶体取向差时的平均对流速度
Fig.8 Average convection velocities under various orientation differences
4 结论
CM247LC合金在温度梯度为10 K/mm、抽拉速度为50 μm/s条件下的相场模拟结果表明，当晶体取向角度不同时，在枝晶生长过程中，液相区域的平均对流速度均表现出周期性变化。在溶质对流条件下，当晶体取向角度较大时，定向凝固的一次枝晶臂间距随晶体取向角度的变大而变大，但与Gandin等[27]提出的幂函数关系存在一定偏离。晶体取向差越大，溶质羽流越容易发起。
文章引用：张勇佳, 周建新, 殷亚军, 等. 高温合金定向凝固过程中枝晶生长与溶质对流数值模拟[J]. 精密成形工程, 2023, 15(10): 13-20.
