高温合金定向凝固过程中枝晶生长与溶质对流数值模拟

张勇佳，周建新，殷亚军，沈旭，计效园，李文

高温合金定向凝固过程中枝晶生长与溶质对流数值模拟

张勇佳，周建新*，殷亚军，沈旭，计效园，李文

（华中科技大学 材料成形与模具技术全国重点实验室，武汉 430074）

针对高温合金叶片在定向凝固过程中容易出现雀斑缺陷，从而导致叶片报废的问题，对定向凝固枝晶生长与溶质对流进行模拟研究，以揭示雀斑缺陷的形成规律。针对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）所示。

采用格子Boltzmann模型对流动场进行求解，选取Bhatnagar-Gross-Krook（BGK）碰撞算子形式的格子Boltzmann方程，如式（5）所示。

式中：为网格所在的空间位置矢量；f为方向的粒子分布函数；feq为方向的平衡态粒子分布函数；Δ为时间步长；为松弛时间；F为外力项；为方向的离散格子速度。平衡态分布函数表达式如式（6）所示。

式中：为流体密度；为流体速度矢量；s为格子声速；w为权系数。这里采用D2Q9格子速度模型，w的取值如式（7）所示。

采用Boussinesq假设来处理溶质差引起的浮力，其计算式如式（8）所示。

采用有限体积法对相场方程和溶质场的对流扩散方程进行求解。采用一阶Euler格式对时间导数进行离散，采用迎风格式离散对流项，以文献[24]中的通量形式对式（4）的右端项进行处理。为减小计算量，分别在粗、细2套网格上求解格子Boltzmann方程和相场方程。参考文献[25]，粗细网格的边长比值为2，采用GPU并行算法进行求解，粗、细网格上的场变量需要进行相应的插值。

2 定向凝固过程中枝晶生长的计算参数

选取CM247LC合金进行定向凝固枝晶生长模拟，合金成分如表1所示。采用伪二元近似方法[26]得到相场模型所需的液相线斜率和平衡分配系数等参数。相场的计算网格尺寸为1 μm，流动场的计算网格尺寸为2 μm，界面宽度0为1.25 μm，时间步长为2.0×10−5s，界面能各向异性系数为0.02，液相溶质扩散系数为3.6×10–9m2/s，液相的运动黏度为1.0× 10–6m2/s，重力加速度常量为−9.81 m/s2，耦合系数=33.27。前期研究表明[25]，该耦合系数的取值能够保证相场模拟结果具备收敛性。当耦合系数采用该取值时，CMSX-4合金定向凝固枝晶生长的模拟结果表明，在溶质羽流发起位置的枝晶生长速度呈周期性振荡，与同步辐射实验[6]观测到的振荡周期与幅值均具有较好的一致性[25]。因此，本文计算模型在模拟枝晶生长与溶质对流相互作用方面具有较高的准确性。

表1 CM247LC合金成分

Tab.1 Composition of CM247LC alloy wt.%

模拟的工艺参数如下：=10 K/mm，p=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中的为无量纲溶质浓度。可以看到，在凝固初期（40 s时），初始条件的平界面发展为树枝晶，且由于竞争生长，一次枝晶臂间距逐渐变大。当凝固时间为100 s时，凝固前沿的溶质在对流作用下富集，并形成溶质羽流，如图2b中左起第2、3个枝晶间。可以看出，该位置的枝晶臂间距较大。当枝晶臂间距较大时，溶质对流的阻力较小，因而更容易形成溶质羽流。随着凝固过程的进行（150 s和200 s），形成的溶质羽流由于具有较低的密度而上浮，在上浮过程中，羽流中心的溶质浓度在对流扩散作用下有所降低，因而羽流上浮的速度有所减缓，开始在环流式对流的作用下向两侧运动，使得溶质羽流呈现为“烟囱”状。模拟结果表明，最大对流速度可达160 μm/s，远大于该工艺参数条件下枝晶生长的稳态速度50 μm/s。

图1 计算区域初始条件与边界条件

图2 晶体取向角度为5 °时的模拟结果

晶体取向角度为15°时和30°时的模拟结果分别如图3和图4所示。可以看到，模拟结果中均出现了溶质羽流。当晶体取向角度为30°时，一次枝晶臂间距大于取向角度为15°时的一次枝晶臂间距，且在凝固200 s时，晶体取向角度为30°时的对流速度最大可达180 μm/s，略高于15°时的150 μm/s。

图3 晶体取向角度为15°时的模拟结果

图4 晶体取向角度为30 °时的模拟结果

图5 不同晶体取向角度时的平均对流速度vtavg（a）及其时间平均值vavg（b）

统计结果表明，当晶体取向角度从0°变化到30°时，在凝固初期（凝固时间小于50 s左右时），平均对流速度均保持在较小值（低于10 μm/s）。随着凝固过程的进行，平均对流速度迅速增大至峰值，速度峰值超过70 μm/s。之后，平均对流速度开始衰减，并保持周期性振荡。除了晶体取向角度为15°时平均对流速度增大至峰值较为滞后（在120 s左右速度出现明显的增大）外，当晶体取向角度为其他取值时，平均对流速度开始迅速增大的时间均较早（在60～80 s之间）。结果表明，当晶体取向角度不同时，平均对流速度时间平均值均大于50 μm/s，即大于抽拉速度，即对应枝晶生长速度的稳态值。平均对流速度振荡幅值较大（见图5a），其时间平均值的分散度较大，这表明，平均对流速度时间平均值与晶体取向角度之间没有明显的依赖性。

溶质羽流的发起与一次枝晶臂间距密切相关。为此，统计了凝固终态时不同晶体取向角度对应的一次枝晶臂间距的最大值、最小值及均值，结果如图6所示。

图6 不同晶体取向角度时的一次枝晶臂间距的最大值、最小值及均值

晶体取向差为5°时的溶质分布与速度分布如图7所示。结果表明，当晶体取向差为5°时，在右侧发散型晶界处，对流速度较大。当凝固时间为74 s时，在发散型晶界两侧出现明显的溶质浓度起伏（见图7a）。随着凝固过程的进行，在79 s时，发散型晶界右侧率先形成溶质羽流（见图7b）。随后，溶质羽流逐渐向上运动（见图7c），之后，溶质羽流从两侧向下方回流（见图7d），溶质羽流中心的浓度减小。在溶质羽流发起的位置，倾斜枝晶生长出的二次枝晶及其分枝的生长被抑制，形成的枝晶间隙较大，进而有利于溶质羽流的发展。取向差为10°和15°时的模拟结果也呈现出类似的规律。

图7 不同凝固时间下晶体取向差为5°时的模拟结果

统计了不同晶体取向差时的平均对流速度，如图8所示。可以看到，当晶体取向差不同时，平均对流速度均在增大至峰值后，在一定范围内呈周期性振荡，但达到峰值的时间不同。随着晶体取向差变大，平均对流速度达到峰值的时间逐渐减小。从图8可以看出，当晶体取向差为5°、10°、15°时，达到峰值的时间（如图8中竖直黑色虚线所示）分别为98、84、63 s。平均对流速度达到峰值的时间越小，表明对流过程中溶质羽流出现得越早。因此，模拟结果表明，晶体取向差越大，溶质羽流越容易发起。

图8 不同晶体取向差时的平均对流速度

4 结论

CM247LC合金在温度梯度为10 K/mm、抽拉速度为50 μm/s条件下的相场模拟结果表明，当晶体取向角度不同时，在枝晶生长过程中，液相区域的平均对流速度均表现出周期性变化。在溶质对流条件下，当晶体取向角度较大时，定向凝固的一次枝晶臂间距随晶体取向角度的变大而变大，但与Gandin等[27]提出的幂函数关系存在一定偏离。晶体取向差越大，溶质羽流越容易发起。

[1] POLLOCK T M, MURPHY W H. The Breakdown of Single-Crystal Solidification in High Refractory Nickel-base Alloys[J]. Metallurgical and Materials Transactions A, 1996, 27(4): 1081-1094.

[2] TIN S, POLLOCK T M, MURPHY W. Stabilization of Thermosolutal Convective Instabilities in Ni-based Single-Crystal Superalloys: Carbon Additions and Freckle Formation[J]. Metallurgical and Materials Transactions A, 2001, 32(7): 1743-1753.

[3] TIN S, POLLOCK T M, LEE P D, et al. Predicting Freckle Formation in Single Crystal Ni-base Superalloys[J]. Journal of Materials Science, 2004, 39(24): 7199-7205.

[4] MA D X, WU Q, BÜHRIG-POLACZEK A. Some New Observations on Freckle Formation in Directionally Solidified Superalloy Components[J]. Metallurgical and Materials Transactions B, 2012, 43(2): 344-353.

[5] SHEVCHENKO N, BODEN S, GERBETH G, et al. Chimney Formation in Solidifying Ga-25Wt Pct in Alloys under the Influence of Thermosolutal Melt Convection[J]. Metallurgical and Materials Transactions A, 2013, 44(8): 3797-3808.

[6] REINHART G, GRANGE D, ABOU-KHALIL L, et al. Impact of Solute Flow During Directional Solidification of a Ni-Based Alloy: In-Situ and Real-Time X-Radio­graphy[J]. Acta Materialia, 2020, 194: 68-79.

[7] SCHNEIDER M C, GU J P, BECKERMANN C, et al. Modeling of Micro- and Macrosegregation and Freckle Formation in Single-Crystal Nickel-base Superalloy Directional Solidification[J]. Metallurgical and Materials Transactions A, 1997, 28(7): 1517-1531.

[8] FELICELLI S D, HEINRICH J C, POIRIER D R. Three-Dimensional Simulations of Freckles in Binary Alloys[J]. Journal of Crystal Growth, 1998, 191(4): 879- 888.

[9] YUAN L, LEE P D. A New Mechanism for Freckle Initiation Based on Microstructural Level Simulation[J]. Acta Materialia, 2012, 60(12): 4917-4926.

[10] KARAGADDE S, YUAN L, SHEVCHENKO N, et al. 3-D Microstructural Model of Freckle Formation Validated Using in Situ Experiments[J]. Acta Materialia, 2014, 79: 168-180.

[11] KAO A, SHEVCHENKO N, ALEXANDRAKIS M, et al. Thermal Dependence of Large-Scale Freckle Defect Formation[J]. Philosophical Transactions A, 2019, 377(2143): 20180206.

[12] SAKANE S, TAKAKI T, ROJAS R, et al. Multi-GPUs Parallel Computation of Dendrite Growth in Forced Convection Using the Phase-Field-Lattice Boltzmann Model[J]. Journal of Crystal Growth, 2017, 474: 154-159.

[13] GUO Z P, MI J W, GRANT P S. An Implicit Parallel Multigrid Computing Scheme to Solve Coupled Thermal-Solute Phase-Field Equations for Dendrite Evolution[J]. Journal of Computational Physics, 2012, 231(4): 1781-1796.

[14] GUO Z P, MI J W, XIONG S M, et al. Phase Field Study of the Tip Operating State of a Freely Growing Dendrite against Convection Using a Novel Parallel Multigrid Approach[J]. Journal of Computational Physics, 2014, 257: 278-297.

[15] YANG C, XU Q Y, LIU B C. Study of Dendrite Growth with Natural Convection in Superalloy Directional Solidification via a Multiphase-Field-Lattice Boltzmann Model[J]. Computational Materials Science, 2019, 158: 130-139.

[16] TAKAKI T, SAKANE S, OHNO M, et al. Competitive Growth during Directional Solidification of a Binary Alloy with Natural Convection: Two-Dimensional Phase-Field Study[J]. Modelling and Simulation in Materials Science and Engineering, 2019, 27(5): 54001.

[17] TAKAKI T, SAKANE S, OHNO M, et al. Large-Scale Phase-Field Lattice Boltzmann Study on the Effects of Natural Convection on Dendrite Morphology Formed During Directional Solidification of a Binary Alloy[J]. Computational Materials Science, 2020, 171: 109209.

[18] SAKANE S, TAKAKI T, OHNO M, et al. Acceleration of Phase-Field Lattice Boltzmann Simulation of Dendrite Growth with Thermosolutal Convection by the Multi-GPUs Parallel Computation with Multiple Mesh and Time Step Method[J]. Modelling and Simulation in Materials Science and Engineering, 2019, 27(5): 54004.

[19] WATANABE S, AOKI T, TAKAKI T. A Domain Partitioning Method Using a Multi-Phase-Field Model for Block-based AMR Applications[J]. Parallel Computing, 2020, 97: 102647.

[20] ROJAS R, TAKAKI T, OHNO M. A Phase-Field-Latt­ice Boltzmann Method for Modeling Motion and Growth of a Dendrite for Binary Alloy Solidification in the Presence of Melt Convection[J]. Journal of Computational Physics, 2015, 298: 29-40.

[21] TAKAKI T, SATO R, ROJAS R, et al. Phase-Field Lattice Boltzmann Simulations of Multiple Dendrite Growth with Motion, Collision, and Coalescence and Subsequent Grain Growth[J]. Computational Materials Science, 2018, 147: 124-131.

[22] ZHANG Q, SUN D, PAN S, et al. Microporosity Formation and Dendrite Growth During Solidification of Aluminum Alloys: Modeling and Experiment[J]. International Journal of Heat and Mass Transfer, 2020, 146: 118838.

[23] ZHANG Q, SUN D, ZHANG S, et al. Modeling of Microporosity Formation and Hydrogen Concentration Evolution during Solidification of an Al-Si Alloy[J]. Chinese Physics B, 2020, 29(7): 78104.

[24] PROVATAS N, ELDER K. Phase-Field Methods in Materials Science and Engineering[M]. Weinheim: Wiley-VCH, 2010.

[25] ZHANG Y, ZHOU J, YIN Y, et al. Study on the Solutal Convection during Dendrite Growth of Superalloy under Directional Solidification Condition[J]. Journal of Materials Research and Technology, 2023, 23: 3916-3927.

[26] RAGHAVAN S, SINGH G, SONDHI S, et al. Construction of a Pseudo-Binary Phase Diagram for Multi- Component Ni-base Superalloys[J]. Calphad, 2012, 38: 85-91.

[27] GANDIN C A, ESHELMAN M, TRIVEDI R. Orientation Dependence of Primary Dendrite Spacing[J]. Metallurgical and Materials Transactions A, 1996, 27(9): 2727-2739.

Numerical Simulation of Dendrite Growth and Solute Convection during Directional Solidification of Superalloy

ZHANG Yong-jia, ZHOU Jian-xin*, YIN Ya-jun, SHEN Xu, JI Xiao-yuan, LI Wen

(State Key Laboratory of Materials Processing and Die & Mould Technology, Huazhong University of Science and Technology, Wuhan 430074, China)

Superalloy blades are prone to freckle defects during directional solidification, which leads to scrap. Therefore, the work aims to simulate the dendrite growth and solute convection during directional solidification to reveal the formation of freckle defects. The dendrite growth of CM247LC alloy was simulated by phase field model and natural convection caused by solute concentration difference was simulated by the lattice Boltzmann model. The phase field-lattice Boltzmann model was solved by GPU-based parallel algorithm based on two different meshes. The dendrite morphology, convection velocity and evolution of solute plumes under different crystal orientation angles and orientation differences were studied. The average fluid velocity in the liquid phase during dendrite growth varied periodically with different crystal orientation angles. When the crystal orientation angle was large, the distance between the primary dendrite arm spacings increased with the increase of the crystal orientation angle. When there was a difference in crystal orientation between dendrites, solute plumes tended to start near divergent grain boundaries. With the increase of crystal orientation difference, the onset time of solute plume advanced. The formation of solute plumes inhibited the growth of dendrite tips and adjacent dendrite side arms. The angle of crystal orientation has slight effect on the formation of solute plumes, and the larger difference of crystal orientation can promote the formation of solute plumes.

superalloy; directional solidification; dendrite growth; solute convection; phase field simulation

10.3969/j.issn.1674-6457.2023.010.002

TG244.3

A

1674-6457(2023)010-0013-08

2023-08-04

2023-08-04

国家重点研发计划（2020YFB1710100）

The National Key R&D Program of China (2020YFB1710100)

张勇佳, 周建新, 殷亚军, 等. 高温合金定向凝固过程中枝晶生长与溶质对流数值模拟[J]. 精密成形工程, 2023, 15(10): 13-20.

ZHANG Yong-jia, ZHOU Jian-xin, YIN Ya-jun, et al. Numerical Simulation of Dendrite Growth and Solute Convection during Directional Solidification of Superalloy[J]. Journal of Netshape Forming Engineering, 2023, 15(10): 13-20.

责任编辑：蒋红晨