The effect of natural convection on the dendritic tip stability during directional solidification: insights from phase-field-lattice Boltzmann simulations

INTRODUCTION

Nickel (Ni)-based superalloys have been extensively employed in critical hot-section components such as turbine blades of aero-engines and industrial gas turbines due to their exceptional high-temperature mechanical properties and superior corrosion resistance^[1,2]. The service performance of these components is largely dictated by the microstructure formed during casting, which is strongly influenced by dendritic growth behavior and its coupling with solutal convection during solidification^[3,4]. Under terrestrial gravity, local density variations arising from temperature and solute gradients drive natural convection, which not only alters the morphology and distribution of dendrites but also profoundly affects dendrite arm spacing and elemental segregation^[5,6]. Such thermosolutal coupling-induced disturbances often trigger severe solidification defects, such as freckles, which undermine the mechanical integrity of the material and shorten the service life of components^[7,8]. Therefore, the fundamental understanding of the effect of natural convection on the dendritic growth dynamics mechanisms during directional solidification is of great engineering and theoretical significance for mitigating microstructural defects in turbine blade casting.

During directional solidification, the dendritic growth behavior and solute distribution characteristics of multicomponent high-temperature alloys are significantly influenced by processing parameters such as the thermal gradient and pulling velocity^[5,9]. Previous studies demonstrated that for relatively low thermal gradient and pulling velocity, the solutal-driven natural convection leads to the formation of solutal plumes, and subsequently induces classic solidification defects such as freckles^[10]. With the continued advancement of synchrotron X-ray imaging technologies, in-situ observation of dendrite growth during solidification of alloys has become increasingly feasible^[11-14]. For instance, Shevchenko et al.^[15,16] employed X-ray microscopy to investigate the directional solidification of Ga-In alloys, revealing that natural convection induced by density differences between Ga and In leads to solutal plumes, which in turn generate segregation channels and initiate freckle formation. Ruvalcaba et al.^[17] utilized synchrotron radiation to observe dendrite fragmentation in Al-Cu alloys, which is attributed to localized solute enrichment during directional solidification. Furthermore, Reinhart et al.^[18] investigated the solidification behavior of high-temperature nickel-based single crystal alloy CMSX-4 using synchrotron X-ray techniques. They reported that solute convection markedly influences dendrite growth velocity, with dendritic fragmentation occurring during the final solidification of residual liquid. These experimental findings consistently demonstrate the prevalence of solutal plumes, dendrite remelting, and overgrowth phenomena during alloy solidification. However, it is still quite challenging to deeply reveal the effect of convection on the dendritic growth mechanism in directional solidification of nickel-based superalloys by experimental means.

With the continuous advancement of computational technologies, numerical simulation has gradually emerged as a critical approach for investigating alloy solidification processes. In parallel with physics-based modeling, machine learning techniques have also been increasingly explored to assist in microstructure characterization and process-property correlation in metallic systems, benefiting from their ability to handle complex, high-dimensional datasets^[19,20]. Nevertheless, a mechanistic understanding of microstructural evolution under coupled thermosolutal convection still relies predominantly on physically grounded numerical frameworks. As a widely accepted and accurate numerical method, the phase-field method that recovers the Gibbs-Thomson effect with high accuracy can simulate the evolution of complicated phase change microstructures. In the phase-field method, the interface tracking of the complex dendritic contour can be avoided by introducing the phase-field variable φ that smoothly but sharply varies from -1 (liquid phase) to 1 (solid phase). By treating the melt as an incompressible fluid flow, incorporating the Navier-Stokes (N-S) equations into the phase-field method is a direct and conventional approach for investigating the effect of convection on dendritic growth during solidification^[21-23]. However, N-S solvers, which are typically based on global schemes such as finite element or finite difference methods, often suffer from low computational efficiency and numerical instability, particularly in regions with high solid fraction^[24].

In recent years, the lattice Boltzmann method (LBM) has emerged as an alternative computational tool to predict the fluid flow in the melt by streaming and collision of a collection of pseudoparticles^[25]. The first phase-field lattice-Boltzmann (PF-LBM) scheme was adopted by Miller et al.^[26,27]. Subsequently, numerous attempts that combine the phase-field method and the LBM have been carried out to effectively simulate the dendritic growth during solidification of alloys under forced and natural convection^[28,29]. For instance, Sun et al.^[30,31] developed a PF-LBM model that couples dendrite growth with the flow field, enabling in-depth investigation of thermosolutal convection effects during alloy solidification. Zhang et al.^[32,33] integrated the PF-LBM with parallel adaptive mesh refinement to investigate the effects of natural and forced convection on dendritic evolution in Al-Cu alloys. The results demonstrate that convection significantly modifies solute transport dynamics, thereby influencing the development of dendritic arms and the competitive growth among multiple dendrites. Qiu et al.^[34] employed a phase-field-lattice Boltzmann coupling method combined with parallel adaptive mesh refinement to simulate dendritic growth in Al-Cu alloys under fully coupled thermal-solute-convection conditions, thereby elucidating the coupled influence of latent heat release and convection on dendrite morphology and growth kinetics. Zhang et al.^[35] adopted a phase-field-lattice Boltzmann coupled framework to simulate dendritic growth during the directional solidification of nickel-based high-temperature alloys, revealing the mechanisms of oscillatory dendrite growth and freckle defect formation induced by solutal convection, and systematically analyzing the effects of temperature gradients and pulling rates. Yang et al.^[36] analyzed the dendritic morphology and fluid velocity within the mushy zone during directional solidification of nickel-based superalloys. The study revealed that freckles tend to form when the average fluid velocity in the mushy zone exceeds the withdrawal velocity. Similarly, Takaki et al.^[37-39] successfully employed a PF-LBM approach to simulate the influence of natural convection on dendritic morphologies during directional solidification of binary alloys. These studies have not only validated the effectiveness of the PF-LBM method in simulating complex solidification dynamics involving fluid flow but also provided new theoretical perspectives and computational pathways for understanding multiphysical interactions in alloy solidification. Using the PF-LBM method, Ma et al.^[40] have dynamically reproduced dendritic evolution during solidification and quantitatively analyzed the nonlinear coupling between solute segregation and melt convection. These simulation results provide essential theoretical support for understanding the formation mechanisms of freckle defects and lay a solid foundation for the optimization of directional solidification processes.

Despite significant advancements in experimental and numerical investigations, the dynamic mechanism through which natural convection affects the stability of the dendritic tip during directional solidification is unclear. Herein, a quantitative multi-order parameter phase-field model is coupled with the LBM to simulate dendritic growth under natural convection during directional solidification of a multicomponent alloy. A pseudo-binary approximation is employed, with a specific focus on the morphological instability of dendritic tips coupled with solute transport behavior under varying gravitational accelerations. Section "MATHEMATICAL MODEL" provides a concise mathematical description of the pseudo-binary approximation, the phase-field model and the LBM. We present results and discussion in Section "RESULTS AND DISCUSSION". In Section "CONCLUSIONS", conclusions from numerical simulations are presented.

MATHEMATICAL MODEL

In this paper, we adopt the quantitative multi-order parameter phase-field model proposed by Yang et al.^[41] to investigate dendritic growth during directional solidification. Here, we focus on the directional solidification of the Ni-based superalloy CMSX-4. A pseudo-binary approximation is employed to simplify the multi-component superalloy. The CMSX-4 alloy can be treated as eight Ni-X binary alloys (Cr 6.5 wt.%, Mo 0.6 wt.%, Ti 1 wt.%, W 6.4 wt.%, Co 9.6 wt.%, Al 5.6 wt.%, Ta 6.5 wt.%^[42]), and further reduced to an equivalent Ni-X binary alloy. Hence, the corresponding equivalent mass fractions c₀, the liquidus slope m, partition coefficient k, and solutal partition coefficient β_c are determined using^[43]

where i is the physical properties in the i-th binary system, and N is the number of elements in the alloys. Table 1 lists chemical composition, slopes of liquidus, partition coefficients, and solutal expansion coefficients of the elements in CMSX-4. Hence, in this paper, c₀ = 38.2 wt.%, m = 1.7023 K/wt.%, k = 0.688, β_c = 2.986 × 10^-2 wt.%^-1.

In this model, the phase-field variable ϕ_i is introduced to distinguish the region within the i-th grain (ϕ_i = 1) from the other regions including the liquid and other grains (ϕ_i = -1). The evolution equations for the ϕ_i and solute concentration are expressed as

with

where λ is the coupling parameter between the phase-field equation and the solute concentration equation, D_l is the solute diffusivity in the liquid phase, W₀ = λd₀/α₁ and τ₀ = α₂λW₀²/D_l are the interface width and the relaxation time scale, respectively, α₁ = 0.8839 and α₂ = 0.6267 are two constants in the model, d₀ =Γ/[|m|c₀(1/k -1)] is the capillary length, Γ is the Gibbs-Thomson coefficient, T_M is the melting temperature of pure nickel, T is the temperature field with T = T₀ + Gy - R_pt in the y-direction, T₀ is the initial temperature (t = 0) at the bottom of the temperature field, f_l = (1 - Σ_iϕ_i)/2 is the liquid fraction, G is the temperature gradient along the y-axis, R_p is cooling rate, τ(n) = τ₀a_s(n)∇ϕ_i², a_s(n) = 1 + ε_scos4(θ - θ_i) is the anisotropic function, ε_s is the anisotropic strength of the surface tension, θ is the angle between the normal to the solid-liquid interface and a fixed crystalline axis, θ_i is the crystallographic orientation of i-th grain.

The LBM is employed to simulate fluid flow driven by temperature and solute concentration gradients. By using the Bhatnagar-Gross-Krook collision operator, the fluid flow can be expressed by the following lattice Boltzmann equation^[51]

where f_i(r,t) is the particle distribution, i denotes the i-th lattice direction, r is position vector, τ_f = 3ν/c²δt + 0.5 is the relaxation time with respect to the kinematic viscosity ν, and f_i^eq(r,t) is equilibrium distribution function. The imposed external force, i.e., the force term F_i(r,t) in Equation (8) drives the fluid flow in the melt. The equilibrium particle distribution function is given by^[51]

where the particle density $$\rho_{L B}= {\textstyle \sum_{𝑖=0}^{𝑄−1}}{f_i}$$, c_i is the discrete velocity in the i-direction, ν_LB = Σ_if_ic_i/ρ + δtF_i/(2ρ) is the dimensional flow velocity, and c = δx/δt is the lattice speed. In the D2Q9 model^[52], ω₀ = 4/9, ω_1-4 = 1/9, ω_5-8 = 1/36, and the discrete velocity is given as

The discrete volumetric force F_i is given by^[53]

with

where F_b and F_d are the dissipative drag force and the buoyancy force. The drag force F_b is introduced to account for momentum dissipation in regions with high solid fraction. This term effectively suppresses unphysical penetration of melt flow into solid or near-solid regions and ensures a smooth transition from free flow in liquid channels to full flow blockage in the solid phase, thereby recovering the correct no-slip boundary condition at the sharp-interface limit^[54]. Here, h = 2.757 is a dimensionless constant^[53], g = [0, g_y] is gravitational acceleration, and hence g_y > 0 for the gravitational acceleration along y⁺, and g_y < 0 for the gravitational acceleration along y^-, β_T and β_C are the expansion coefficients for temperature and solute, respectively. The coupled equations including Equations (5), (6) and (8) are implemented for parallel computing on multi-graphics processing units (Multi-GPU) by using in-house code written in modular Compute Unified Device Architecture (CUDA). In our numerical simulations, the grid size and the time step are Δx/W₀ = 0.7 and Δt/τ₀ = 0.001, respectively.

RESULTS AND DISCUSSION

In this paper, we assume that the physical parameters of CMSX-4 are constant, and physical properties of CMSX-4 and numerical conditions are presented in Table 2. The supersaturated liquid has u₀ = -0.45, and the temperature on the bottom boundary is T = 1,662.88 K initially. The thermal gradient is applied along the vertical direction (y-axis) with a magnitude of 0.005 K/μm, and the cooling rate is set to 5 K/s. Hence, the growth velocity of the columnar dendrite is V_p = R/G = 1,000 μm/s. Figure 1 shows the computational domains and initial conditions employed in this study. Two typical types of simulations were conducted: single crystalline configuration and bi-crystalline configuration. No-flux boundary conditions are applied on the top and bottom boundaries of the numerical domain, while periodic boundary conditions are applied on the left and right boundaries. For the single-crystalline configuration, the computational grid size is 200 × 1,000, and a semicircular solid seed with an initial radius r₀ = 15W₀ and θ_i = 0° was placed on the center of the bottom boundary. For the bi-crystalline configuration, the computational grid size was 2,400 × 1,000, and twelve semicircular solid seeds are uniformly distributed along the bottom boundary of the computational domain. Six adjacent seeds share the same crystallographic orientation angle θ₁, forming one grain, while the remaining six adjacent seeds have an orientation angle θ₂, forming the neighboring grain. The centers of all seeds are uniformly spaced along the horizontal direction.

Figure 1. The initial condition and boundary condition of the computational domain: (A) Single crystalline configuration; (B) bi-crystalline configuration.

In this study, the gravitational acceleration g is used as a control parameter to systematically regulate both the intensity and direction of solutal buoyancy-driven natural convection. The gravitational acceleration in the simulation is defined as g = [0, g_y], and g_y = -g₀, -g₀/2, -g₀/5, g₀/5, g₀/2, and g₀. When all other parameters are fixed, the thermal Rayleigh number Ra_T = gβ_TΔcL³/να and the solutal Rayleigh number Ra_c = gβ_cΔcL³/νD_l, which characterize buoyancy-driven flow induced by temperature and solute concentration gradients, respectively, are linearly proportional to g, thereby representing the strength of the buoyancy-driven forces, where ΔT and Δc denote the temperature difference and solute concentration difference, L is a characteristic length, andα is the thermal diffusivity. Consequently, both positive and negative values of the gravitational acceleration are selected, corresponding to different convection intensities and directions, in order to elucidate the underlying mechanisms by which natural convection influences solute transport, local undercooling, and dendritic tip stability.

Firstly, we focus on the effect of natural convection on the dendritic tip instability in the single crystalline configuration. In this study, the dendrite with θ₁ = 0° grows from the initial seed in the stationary state until a steady dendritic morphology forms (t′ = 3.56 × 10⁶); then, the natural convection was introduced. Note that this numerical processing differs from practical directional solidification, where buoyancy-driven flow develops continuously from the onset of solidification. However, this strategy enables a clear understanding of the instability mechanism of the dendritic tip under natural convection during directional solidification. Figure 2 shows the snapshots of a single dendrite in the downward buoyancy ((A) g_y = -g₀/5, (B) g_y = -g₀/2 and (C) g_y = -g₀) and the upward buoyancy ((D) g_y = g₀/5, (E) g_y = g₀/2 and (F) g_y = g₀) after the natural convection is introduced, where the grayscale indicates the solutal concentration field, and the flow velocity is marked by the colored arrows. Meanwhile, Figure 3 shows the evolution of the solute concentration profiles extracted along the dendrite symmetry axis for the simulation cases presented in Figure 2. Obviously, the strength of natural convection increases with the magnitude of the imposed gravitational acceleration g_y. Because the thermal gradient used in our simulations is relatively small (0.005 K/μm), natural convection is mainly caused by the uneven distribution of solute. The solute expansion coefficient β_c = 2.986 × 10^-2 wt.%^-1 implies that solute atoms are lighter, and hence the rejected solute during dendritic growth will move in the opposite direction of gravity.

Figure 2. Snapshots of a single dendrite in the downward buoyancy ((A) g_y = -g₀/5, (B) g_y = -g₀/2 and (C) g_y = -g₀) and the upward buoyancy ((D) g_y = g₀/5, (E) g_y = g₀/2 and (F) g_y = g₀) at t′ = 3.92 × 10⁶ to 6.67 × 10⁶, where the grayscale indicates the solutal concentration field, the flow velocity is marked by the colored arrows.

Figure 3. Evolution of the solute distribution along the crystal symmetry axis: (A) g_y = -g₀/5, (B) g_y = -g₀/2 and (C) g_y = -g₀, (D) g_y = g₀/5, (E) g_y = g₀/2 and (F) g_y = g₀.

As shown in Figure 2A-C, for g_y < 0, i.e., the gravitational acceleration along y^-, the presence of the natural convection transports the solute atoms from the dendritic root towards the dendritic tip. The enrichment of the solute concentration near the dendritic tip lowers the local undercooling and hinders the growth of the dendritic tip. For g_y = -g₀/5 (see Figure 2A), the effect of natural convection on the dendritic shape is not significant as the convection intensity is relatively weak. As shown in Figure 3A, although the maximum of the solute concentration increases and the solute boundary layer thickens in front of the dendritic tip with the time, the increase is not significant. However, the convection intensity increases with the increase of g_y. For g_y = -g₀/2 (see Figure 2B), more solute atoms are transported from the root towards the dendritic tip. As shown in Figure 3B, the maximum of the solute concentration in front of the dendritic tip significantly increases with the time. Therefore, the enrichment of the solute concentration can significantly slow down the growth of the dendritic tip according to the corresponding equivalent liquidus slope and solute partition coefficient of the equivalent binary alloy. When g_y = -g₀ (see Figure 2C), the enrichment of the solute concentration can result in the occurrence of dendritic tip splitting. Moreover, as solute-rich liquid flows up, the so-called “plume” forms in front of the dendritic tip, as shown in Figure 3C.

As for g_y > 0, as shown in Figure 2D, the presence of natural convection can transport rejected solute atoms from the dendritic tip region towards the dendritic root, resulting in the enrichment of the solute concentration in the interdendritic region. As shown in Figure 3A-C, the solute concentration has an extremum at the dendritic tip, and all the solute concentration versus distance curves overlap after the natural convection is introduced for the cases of g_y. This means that natural convection has negligible influence on the solute distribution along the crystal symmetry axis for g_y > 0, and the morphology of dendritic tips remains almost unchanged with time although the sidebranches develop just after the convection is introduced.

In order to investigate the influence of natural convection on the dendritic tip instability, the evolution of the dendritic tip radius with time is analyzed after the natural convection is introduced. In this study, the dendritic tip radius is calculated by fitting the contour of ϕ_i = 0 to a fourth-order interpolation function y = a₂x² + a₄x⁴ with R_tip = 1/2a₂. The evolution of the dendritic tip radius over time for g_y > 0 and g_y < 0 is presented in Figure 4A and B, respectively. Here, the evolution of the dendritic tip radius over the time in the static state, i.e., g_y = 0, is also presented for comparative analysis. When g_y = -g₀/5, the dendritic tip radius continuously but slightly decreases. For g_y = -g₀/2 and g_y = -g₀, it can be clearly observed that the dendritic tip radius firstly decreases and then increases rapidly until the tip splitting occurs for g_y < 0. However, for g_y > 0, the dendritic tip radius first increases and then remains constant. This indicates that the natural convection with g_y > 0 can stabilize the dendritic tip for lighter solute atoms while the natural convection with g_y < 0 can result in the dendritic tip instability.

Figure 4. Evolution of the dendritic tip radius over the time: (A) for g_y < 0 and (B) for g_y > 0.

Then, let us turn to the competitive growth of bi-crystalline grains regarding the evolution of the converging grain boundary under natural convection. Similar to the single configuration, the dendrites grow from the initial seeds in the static state until t′ = 3.56 × 10⁶, and then the natural convection was introduced. Figure 5 shows the snapshots of the bi-crystalline grains with converging grain boundary (θ₁, θ₂) = (-15°, 15°) in the downward buoyancy ((A) g_y = -g₀/5, (B) g_y = -g₀/2 and (C) g_y = -g₀) and the upward buoyancy ((D) g_y = g₀/5, (E) g_y = g₀/2 and (F) g_y = g₀) after the natural convection was introduced. The primary spacing of the dendritic arrays is history-dependent; i.e., it is determined by the spacing of the initial solid seeds. It can be seen clearly that the natural convection has negligible effect on the dendritic shape and the positions of the dendritic tips in the downward buoyancy except for the obvious sidebranches near the left and right boundaries for lower convection intensity (g_y = -g₀/5). With the increase of the convection intensity (g_y = -g₀/2 and g_y = -g₀), solute-rich liquid flows up and plumes form in front of the dendritic arrays, resulting in the instability of the dendritic tips and the occurrence of the continuous tip splitting. In upward buoyancy, regardless of gravity strength, the morphology of dendrites and the position of the grain boundary remain barely affected by natural convection.

Figure 5. Snapshots of the bi-crystalline grains with the converging grain boundary (θ₁, θ₂) = (-15°, 15°) in the downward buoyancy ((A) g_y = -g₀/5, (B) g_y = -g₀/2 and (C) g_y = -g₀) and the upward buoyancy ((D) g_y = g₀/5, (E) g_y = g₀/2 and (F) g_y = g₀) after natural convection is introduced.

Figure 6 shows the snapshots of the bi-crystalline grains with diverging grain boundary (θ₁, θ₂) = (-15°, 15°) in the downward buoyancy ((A) g_y = -g₀/5, (B) g_y = -g₀/2 and (C) g_y = -g₀) and the upward buoyancy ((D) g_y = g₀/5, (E) g_y = g₀/2 and (F) g_y = g₀) after the natural convection is introduced. In the case of the diverging grain boundary, a large liquid zone exists between the two adjacent grains. In the downward buoyancy, as shown in Figure 6A-C, the high solute concentration liquid flows from the grain boundary region to the dendritic tip region and the large liquid zone between the two grains via the two elliptical vortex rings. This results in the development of sidebranching and the generation of primary arms at the diverging grain boundary. With the increase of the convection intensity (g_y = -g₀/2 and g_y = -g₀), the complex natural convection leads to the instability of dendritic tips. As for the upward buoyancy, as shown in Figure 6D-F, the high solute concentration liquid flows from the dendritic tip region into the grain boundary region. Although sidebranches are developed at the diverging grain boundary, the high solute concentration prevents the generation of primary arms at the diverging grain boundary.

Figure 6. Snapshots of the bi-crystalline grains with diverging grain boundary (θ₁, θ₂) = (-15°, 15°) in the downward buoyancy ((A) g_y = -g₀/5, (B) g_y = -g₀/2 and (C) g_y = -g₀) and the upward buoyancy ((D) g_y = g₀/5, (E) g_y = g₀/2 and (F) g_y = g₀) after the natural convection is introduced.

Then, the natural convection is introduced at the beginning of solidification (t′ = 0). The converging and diverging grain boundaries of the bi-crystalline grains are considered. Figure 7 shows the snapshots of the bi-crystalline grains with the converging grain boundary (θ₁, θ₂) = (-15°, 15°) in the downward buoyancy ((A) g_y = -g₀/5, (B) g_y = -g₀/2 and (C) g_y = -g₀) and the upward buoyancy ((D) g_y = g₀/5, (E) g_y = g₀/2 and (F) g_y = g₀). In the downward buoyancy, pronounced vortex structures form in front of the columnar arrays. These vortices act as efficient transport channels, drastically altering the local solute field. The flow sweeps solute from the interdendritic regions and advects it towards the dendrite tips, creating localized solute enrichment or depletion that directly influences the local undercooling. The significant solute-rich plumes form and the number of vortices increases with the convection intensity, enhancing this advective transport. The resulting modification of the solute boundary layer around the tips destabilizes the diffusion-controlled growth, resulting in the splitting of the dendritic tip and the decrease of primary spacing. Moreover, for larger convection intensity (g_y = -g₀), the columnar front becomes undulating and the symmetry of the columnar dendrites is broken. Similar to the predicted results presented in Figure 5, except the development of the sidebranches, the convection intensity has a negligible effect on the primary spacing of the dendritic arrays in the upward buoyancy. Although the natural convection can modify the growth pattern, the grain boundary remains almost unaffected. The red dashed line in Figure 7 marks the position five grids above the dendrite leading tip along the y-direction, corresponding to the grid point N_y = 750. Figure 8 shows the liquid flow variations along the red dashed line in Figure 7, specifically the distribution of the velocity components U and V. In Figure 8, in front of dendrites A and B, downward and inward flows occur, which remove the solute from the tips of dendrites A and B, thereby locally decreasing the solute concentration and increasing the local undercooling, which promotes their preferential growth. On the other hand, in front of dendrites C and D, very fast upward and outward flows are present. These flows transport the solute-rich liquid toward the tips of dendrites C and D. This advective supply of solute creates a localized solute-rich layer at the tips, raising the local liquidus temperature and effectively suppressing the driving force for solidification. As a result, the tip positions of dendrites C and D are lower than those of the other dendrites.

Figure 7. Snapshots of the bi-crystalline grains with the converging grain boundary (θ₁, θ₂) = (-15°, 15°) in the downward buoyancy ((A) g_y = -g₀/5, (B) g_y = -g₀/2 and (C) g_y = -g₀) and the upward buoyancy ((D) g_y = g₀/5, (E) g_y = g₀/2 and (F) g_y = g₀) when the natural convection is introduced at the beginning of the solidification (t′ = 0).

Figure 8. Liquid flow variations along the red dashed line in Figure 7, where the distribution of the velocity components U in the x-direction and V in the y-direction, and arrows labeled with capital letters indicate the positions of the dendrite tips shown in Figure 7.

The diverging case is presented in Figure 9. In the case of the downward buoyancy, two elliptical vortex rings form in the liquid zone between the two grains. These large-scale vortices dominate the solute redistribution within the liquid zone. They continuously flush the interdendritic solute from the grain boundary region and carry it towards the dendrite tips on the flanks, thereby creating strong lateral solute gradients. With the increase of the convection intensity, the liquid zone becomes shortened. When g_y = -g₀, the liquid zone between grains almost disappears, and the solidification front at the grain boundary lags behind other positions. In the case of the upward buoyancy, the liquid zone between the two grains remains almost unchanged, and the increase of the convection intensity promotes the coarsening of secondary sidebranches. Figure 10 shows the flow distribution along the red dashed line in Figure 9, with the velocity component U in the x-direction and the velocity component V in the y-direction. In Figure 9, the arrows with capital letters indicate the positions of the dendrite tips. In front of dendrites E and H, very fast upward and outward flows occur. These flows transport the solute-rich liquid toward the tips of dendrites E and H, enriching the solute boundary layer and stabilizing it against the development of new sidearms, while simultaneously retarding the primary tip advancement, which results in the solute accumulation. Thus, the tip positions of dendrites E and H are lower than those of the other dendrites. In contrast, in front of dendrites F and G, very fast downward and inward flows occur. These flows remove the solute around the tips of dendrites F and G, effectively thinning the solute boundary layer, increasing the local undercooling, and creating conditions favorable for the initiation and sustained growth of sidearms that allow them to grow preferentially. The asymmetric flow pattern thus establishes a competitive growth environment where solute redistribution via convection directly dictates which dendrite arms are favored or suppressed, leading to the observed sidearm development and tip instability.

Figure 9. Snapshots of the bi-crystalline grains with diverging grain boundary (θ₁, θ₂) = (-15°, 15°) in the downward buoyancy ((A) g_y = -g₀/5, (B) g_y = -g₀/2 and (C) g_y = -g₀ ) and the upward buoyancy ((D) g_y = g₀/5, (E) g_y = g₀/2 and (F) g_y = g₀) when the natural convection is introduced at the beginning of the solidification (t′ = 0).

Figure 10. Liquid flow variations along the red dashed line in Figure 9, where the distribution of the velocity components U in the x-direction and V in the y-direction, and arrows labeled with capital letters indicate the positions of the dendrite tips shown in Figure 9.

These flow characteristics revealed by the velocity distributions in Figures 8 and 10 provide a foundation for a more general mechanistic understanding of how grain-boundary geometry modulates convection-driven solute transport and dendritic growth behavior. For converging grain boundaries, the narrowing liquid channel enhances the interaction between adjacent solute boundary layers and promotes the formation of coherent upward plumes under downward buoyancy. This geometric confinement amplifies solute accumulation and flow instability in front of the dendritic array, making converging boundaries particularly sensitive to convection-induced tip destabilization. In contrast, diverging grain boundaries provide an expanded liquid region that favors the formation of closed recirculating vortices rather than vertically extended plumes. These vortical flows redistribute solute laterally along the grain boundary, promoting sidebranch development while suppressing the sustained solute accumulation required for primary arm formation. Overall, the present results demonstrate that the influence of natural convection on dendritic growth is governed not only by its intensity and direction, but also by the geometrical characteristics of the grain boundary region. Convection-driven solute redistribution near dendrite tips and grain boundaries, together with the associated changes in local undercooling, provides a consistent mechanistic basis for understanding the observed variations in tip stability, sidebranching, and grain-boundary evolution under different gravitational conditions.

CONCLUSIONS

In this paper, the coupled PF-LBM method is employed to reliably simulate complex dendritic growth phenomena including solutal plumes, vortex structures, and tip splitting under natural convection during directional solidification of alloys. We implement the coupled equations by using in-house code written in CUDA. Compared with previous studies, this work elucidates the physical mechanisms by which natural convection regulates dendrite tip stability. It is found that in the single crystalline case the impact of natural convection on the dendritic tip stability depends strongly on gravitational direction and intensity. Downward buoyancy promotes solute accumulation at tips, enhancing tip instability, while upward buoyancy moves solute away from tips to interdendritic regions, thereby stabilizing the dendritic tips. In the bi-crystalline case, the solute transport and dendritic growth depend strongly on the grain boundary type (converging or diverging). Converging boundaries are sensitive to the downward buoyancy-induced instability, while diverging boundaries form vortex-driven solute flow that promotes sidebranching and primary arms in the downward buoyancy or inhibits primary arms in upward buoyancy. The proposed PF-LBM framework provides a robust and efficient platform for quantitatively exploring the interplay between fluid dynamics and microstructural evolution, offering valuable theoretical insights for optimizing solidification parameters and minimizing defect formation in industrial casting processes.