高溫合金定向凝固過程中枝晶生長與溶質對流數值模擬

張勇佳，周建新，殷亞軍，沈旭，計效園，李文

高溫合金定向凝固過程中枝晶生長與溶質對流數值模擬

張勇佳，周建新*，殷亞軍，沈旭，計效園，李文

（華中科技大學 材料成形與模具技術全國重點實驗室，武漢 430074）

針對高溫合金葉片在定向凝固過程中容易出現雀斑缺陷，從而導致葉片報廢的問題，對定向凝固枝晶生長與溶質對流進行模擬研究，以揭示雀斑缺陷的形成規(guī)律。針對CM247LC合金定向凝固過程，采用相場模型模擬凝固過程枝晶生長，采用格子Boltzmann模型模擬溶質濃度差引起的自然對流。采用基于雙重網格的GPU并行算法對相場-格子Boltzmann模型進行數值求解。研究在不同晶體取向角度與取向差條件下的枝晶形貌、對流速度及溶質羽流的演變規(guī)律。當晶體取向角度不同時，在枝晶生長過程中，液相區(qū)域的平均對流速度均表現為周期性變化。當晶體取向角度較大時，隨著晶體取向角度的變大，一次枝晶臂間距變大。當枝晶間存在晶體取向差時，溶質羽流傾向于在發(fā)散型晶界附近發(fā)起；隨著晶體取向差的增大，溶質羽流發(fā)起時間提前。溶質羽流的形成阻礙了枝晶尖端及附近枝晶側臂的生長。晶體取向角度對溶質羽流形成的影響較小，較大的晶體取向差對溶質羽流的形成有促進作用。

高溫合金；定向凝固；枝晶生長；溶質對流；相場模擬

高溫合金定向凝固技術是制造航空發(fā)動機和燃氣輪機渦輪葉片的主要成形工藝。高溫合金葉片包括等軸晶葉片、定向凝固柱狀晶葉片及單晶葉片。其中，定向凝固柱狀晶葉片與單晶葉片消除了橫向晶界，具有優(yōu)異的高溫力學性能。由于高溫合金葉片的結構較為復雜，其凝固過程中的溫度場難以穩(wěn)定控制，所以在定向凝固過程中容易出現雀斑缺陷。定向凝固雀斑缺陷的形成與凝固過程中枝晶組織與溶質對流之間的相互作用有關。

目前，主要采用實驗與數值模擬的方法對雀斑缺陷的形成機制與規(guī)律進行研究。Pollock等[1]研究了高溫合金定向凝固形成的雀斑缺陷，研究發(fā)現，當冷卻速率較低時，一次枝晶臂間距較大，雀斑缺陷更容易形成，且不同高溫合金之間的雀斑缺陷形成傾向存在顯著差異。Tin等[2]研究發(fā)現，將碳的質量分數增大至0.1%能夠顯著降低雀斑缺陷數量，這得益于富Ta的MC型碳化物析出。Tin等[3]采用回歸分析的方法得到了高溫合金試樣中雀斑缺陷數量與糊狀區(qū)密度差及一次枝晶臂間距的關系式，該關系式可用于高溫合金雀斑缺陷形成傾向的評估。Ma等[4]研究了試樣幾何形狀對雀斑缺陷的影響，發(fā)現在試樣截面擴張與收縮位置容易形成雀斑缺陷。Shevchenko等[5]采用同步輻射技術觀察了Ga-In合金定向凝固過程中溶質羽流的形成，結果表明，富集的溶質促進了偏析通道的形成。Reinhart等[6]采用同步輻射技術觀察了CMSX-4合金凝固過程中枝晶生長與溶質羽流的相互作用。

研究者采用數值模擬的方法研究了定向凝固枝晶生長與溶質對流現象，以揭示雀斑缺陷的形成機制。Schneider等[7]基于高溫合金熱力學數據庫對定向凝固過程中溶質羽流的發(fā)起進行了模擬，結果表明，隨著冷卻速率的減小，凝固前沿的溶質分布轉變?yōu)椴环€(wěn)定狀態(tài)，進而出現了通道偏析。Felicelli等[8]對Pb-10%（質量分數）Sn合金定向凝固過程中的溶質對流現象進行了三維模擬，模擬結果表明，在壁面位置出現了間隔一定距離的偏析通道。Yuan等[9]采用元胞自動機模型研究了Pb-Sn合金定向凝固過程中的溶質對流現象，結果表明，由對流引起的枝晶間富集的溶質使枝晶干和二次枝晶重熔，進而形成了通道偏析。Karagadde等[10]采用元胞自動機模型對Ga-25%（質量分數）In合金定向凝固過程進行了模擬，發(fā)現溶質羽流的形成與枝晶取向及一次枝晶臂間距密切相關，且形成溶質羽流的臨界瑞利數為150～170。Kao等[11]研究了固液界面形狀對定向凝固枝晶生長與溶質對流的影響，模擬結果表明，當固液界面為凸界面時，溶質羽流傾向于在中間位置發(fā)起，而當固液界面為凹界面時，則在兩側壁面位置發(fā)起。由于溶質對流過程中的流動場求解計算量較大，開展大規(guī)模模擬需要極大的計算量。為提高模擬速度，Sakane等[12]提出了基于多個GPU的并行求解算法，采用格子Boltzmann模型對溶質差引起的浮力驅動流動過程進行了快速求解。Guo等[13-14]采用自適應網格算法對熱溶質對流條件下的枝晶生長過程進行了并行求解。Yang等[15]提出了基于耦合熱力學數據庫的GPU并行算法，并對CMSX-4合金凝固過程中的枝晶生長與溶質對流現象進行了模擬。Takaki等[16]采用二維相場模型研究了自然對流條件下Al-Cu合金枝晶生長過程中的競爭生長規(guī)律。Takaki等[17]對Al-Cu合金自然對流條件下的枝晶生長過程進行了大規(guī)模相場模擬，得到了不同重力加速度條件下一次枝晶臂隨時間的變化規(guī)律。Sakane等[18]采用多套網格系統(tǒng)對枝晶生長過程中的相場、流動場與溫度場模型進行求解，有效提高了模擬計算速度。Watanabe等[19]提出了一種用于多相場模型的自適應結構化網格區(qū)域分解算法，以提高計算速度。此外，也有模型考慮了凝固過程中的枝晶運動[20]、碰撞[21]及枝晶與氣孔的相互作用[22-23]等。

目前，晶體取向角度與取向差對高溫合金定向凝固枝晶生長與溶質對流的影響規(guī)律尚不清晰，本文采用相場-格子Boltzmann模型對高溫合金定向凝固過程進行了枝晶尺度模擬，得到了不同晶體取向角度與取向差時溶質羽流的形成傾向以及溶質對流速度的變化規(guī)律。

1 定向凝固過程中枝晶生長數學模型

在定向凝固過程中，枝晶生長的數學模型包括相場模型與格子Boltzmann模型。其中，相場模型用于求解凝固過程中相界面的演化與溶質濃度分布，格子Boltzmann模型用于求解溶質差驅動的流動場。相場模型的控制方程如式（1）所示。

采用格子Boltzmann模型對流動場進行求解，選取Bhatnagar-Gross-Krook（BGK）碰撞算子形式的格子Boltzmann方程，如式（5）所示。

式中：為網格所在的空間位置矢量；f為方向的粒子分布函數；feq為方向的平衡態(tài)粒子分布函數；Δ為時間步長；為松弛時間；F為外力項；為方向的離散格子速度。平衡態(tài)分布函數表達式如式（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合金定向凝固枝晶生長的模擬結果表明，在溶質羽流發(fā)起位置的枝晶生長速度呈周期性振蕩，與同步輻射實驗[6]觀測到的振蕩周期與幅值均具有較好的一致性[25]。因此，本文計算模型在模擬枝晶生長與溶質對流相互作用方面具有較高的準確性。

表1 CM247LC合金成分

Tab.1 Composition of CM247LC alloy wt.%

模擬的工藝參數如下：=10 K/mm，p=50 μm/s。前期的計算結果表明[25]，在該工藝條件下，CM247LC合金中溶質羽流的發(fā)起處于臨界穩(wěn)定條件。在該工藝條件下，有利于溶質羽流的形成，從而研究溶質對流條件下晶體取向角度與取向差的影響規(guī)律。在不同晶體取向角度的模擬過程中，在整個計算區(qū)域底部設置一層連續(xù)的固相，取向角度為0°～30°，取值間隔為5°。計算區(qū)域大小為2 048 μm× 4 096 μm。在不同晶體取向差的模擬過程中，計算區(qū)域初始條件與邊界條件設置如圖1所示。在底部區(qū)域設置2種不同取向的晶粒，計算區(qū)域的左右邊界均為周期性邊界條件，底部為等間距放置的初始晶粒，晶粒之間的間距0為穩(wěn)態(tài)時的一次枝晶間距，具體數值由前期的模擬結果得到[25]。在底部左右各設置4個初始晶粒，取向角度為0°，在中間設置8個初始晶粒，取向角度分別為5°、10°和15°，計算區(qū)域大小為4 096 μm×4 096 μm，采用移動計算區(qū)域的方法以減小計算量。當抽拉距離為10 mm時，終止計算。模擬計算在單個NVIDIA V100 GPU上進行。

3 結果與分析

當晶體取向角度為5°時，溶質濃度分布與對流速度分布模擬結果如圖2所示，圖2中的為無量綱溶質濃度?？梢钥吹?，在凝固初期（40 s時），初始條件的平界面發(fā)展為樹枝晶，且由于競爭生長，一次枝晶臂間距逐漸變大。當凝固時間為100 s時，凝固前沿的溶質在對流作用下富集，并形成溶質羽流，如圖2b中左起第2、3個枝晶間?？梢钥闯觯撐恢玫闹П坶g距較大。當枝晶臂間距較大時，溶質對流的阻力較小，因而更容易形成溶質羽流。隨著凝固過程的進行（150 s和200 s），形成的溶質羽流由于具有較低的密度而上浮，在上浮過程中，羽流中心的溶質濃度在對流擴散作用下有所降低，因而羽流上浮的速度有所減緩，開始在環(huán)流式對流的作用下向兩側運動，使得溶質羽流呈現為“煙囪”狀。模擬結果表明，最大對流速度可達160 μm/s，遠大于該工藝參數條件下枝晶生長的穩(wěn)態(tài)速度50 μm/s。

圖1 計算區(qū)域初始條件與邊界條件

圖2 晶體取向角度為5 °時的模擬結果

晶體取向角度為15°時和30°時的模擬結果分別如圖3和圖4所示?？梢钥吹剑M結果中均出現了溶質羽流。當晶體取向角度為30°時，一次枝晶臂間距大于取向角度為15°時的一次枝晶臂間距，且在凝固200 s時，晶體取向角度為30°時的對流速度最大可達180 μm/s，略高于15°時的150 μm/s。

圖3 晶體取向角度為15°時的模擬結果

圖4 晶體取向角度為30 °時的模擬結果

圖5 不同晶體取向角度時的平均對流速度vtavg（a）及其時間平均值vavg（b）

統(tǒng)計結果表明，當晶體取向角度從0°變化到30°時，在凝固初期（凝固時間小于50 s左右時），平均對流速度均保持在較小值（低于10 μm/s）。隨著凝固過程的進行，平均對流速度迅速增大至峰值，速度峰值超過70 μm/s。之后，平均對流速度開始衰減，并保持周期性振蕩。除了晶體取向角度為15°時平均對流速度增大至峰值較為滯后（在120 s左右速度出現明顯的增大）外，當晶體取向角度為其他取值時，平均對流速度開始迅速增大的時間均較早（在60～80 s之間）。結果表明，當晶體取向角度不同時，平均對流速度時間平均值均大于50 μm/s，即大于抽拉速度，即對應枝晶生長速度的穩(wěn)態(tài)值。平均對流速度振蕩幅值較大（見圖5a），其時間平均值的分散度較大，這表明，平均對流速度時間平均值與晶體取向角度之間沒有明顯的依賴性。

溶質羽流的發(fā)起與一次枝晶臂間距密切相關。為此，統(tǒng)計了凝固終態(tài)時不同晶體取向角度對應的一次枝晶臂間距的最大值、最小值及均值，結果如圖6所示。

圖6 不同晶體取向角度時的一次枝晶臂間距的最大值、最小值及均值

晶體取向差為5°時的溶質分布與速度分布如圖7所示。結果表明，當晶體取向差為5°時，在右側發(fā)散型晶界處，對流速度較大。當凝固時間為74 s時，在發(fā)散型晶界兩側出現明顯的溶質濃度起伏（見圖7a）。隨著凝固過程的進行，在79 s時，發(fā)散型晶界右側率先形成溶質羽流（見圖7b）。隨后，溶質羽流逐漸向上運動（見圖7c），之后，溶質羽流從兩側向下方回流（見圖7d），溶質羽流中心的濃度減小。在溶質羽流發(fā)起的位置，傾斜枝晶生長出的二次枝晶及其分枝的生長被抑制，形成的枝晶間隙較大，進而有利于溶質羽流的發(fā)展。取向差為10°和15°時的模擬結果也呈現出類似的規(guī)律。

圖7 不同凝固時間下晶體取向差為5°時的模擬結果

統(tǒng)計了不同晶體取向差時的平均對流速度，如圖8所示。可以看到，當晶體取向差不同時，平均對流速度均在增大至峰值后，在一定范圍內呈周期性振蕩，但達到峰值的時間不同。隨著晶體取向差變大，平均對流速度達到峰值的時間逐漸減小。從圖8可以看出，當晶體取向差為5°、10°、15°時，達到峰值的時間（如圖8中豎直黑色虛線所示）分別為98、84、63 s。平均對流速度達到峰值的時間越小，表明對流過程中溶質羽流出現得越早。因此，模擬結果表明，晶體取向差越大，溶質羽流越容易發(fā)起。

圖8 不同晶體取向差時的平均對流速度

4 結論

CM247LC合金在溫度梯度為10 K/mm、抽拉速度為50 μm/s條件下的相場模擬結果表明，當晶體取向角度不同時，在枝晶生長過程中，液相區(qū)域的平均對流速度均表現出周期性變化。在溶質對流條件下，當晶體取向角度較大時，定向凝固的一次枝晶臂間距隨晶體取向角度的變大而變大，但與Gandin等[27]提出的冪函數關系存在一定偏離。晶體取向差越大，溶質羽流越容易發(fā)起。

[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

國家重點研發(fā)計劃（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.

責任編輯：蔣紅晨