Three-Dimensional Lattice Boltzmann Modeling of Dendritic Solidification under Forced and Natural Convection

Abstract

A three-dimensional (3D) lattice Boltzmann (LB) model is developed to simulate the dendritic growth during solidification of Al-Cu alloys under forced and natural convection. The LB method is used to solve for solute diffusion and fluid flow. It is assumed that the dendritic growth is driven by the difference between the local actual and local equilibrium composition of the liquid in the interface. A cellular automaton (CA) scheme is adopted to capture new interface cells. The LB models for solute transport and fluid flow are first validated against two benchmark problems. The dendrite growth model is also validated with available analytical solutions. The evolution of a 3D dendrite affected by melt convection is investigated. Also, density inversion caused by solute concentration gradient is studied. It is shown that convection can change the kinetics of growth by affecting the solute distribution around the dendrite. In addition, the growth features of two-dimensional (2D) and 3D dendrites are briefly compared. The results show that decreasing undercooling and increasing solute concentration decelerates the growth in all branches of the dendrite. While increasing fluid velocity does not significantly influence upstream and transverse arms, it decreases the growth rate in the downstream direction considerably. The size ratio of the upstream arm to the downstream arm rises by increasing inlet velocity and solute content, and decreasing undercooling. Similarly, in the case of natural convection, redistribution of solute due to buoyancy-induced flow suppresses the growth of the upward arm and accelerates the growth of the downward arm. Considering the advantages offered by the LB method, the present model can be used as a new tool for simulating 3D dendritic solidification under convection.

1. Introduction

Dendrites are common type of microstructure observed in most crystalline materials, especially metallic alloys. Understanding the solidification process and insight into the kinetics of dendritic growth is extremely valuable, as the microstructural characteristics strongly influence the properties of the solidified alloy materials.

Numerical studies were undertaken by several researchers to simulate the microstructural evolution during alloy solidification. The variety of numerical models of dendrite growth at the microscale falls into one of the three dominant categories: those based on the Phase-Field (PF) method [1,2,3,4,5]; models based on the Level Set (LS) method [6,7,8,9], and models that perform a Direct Interface Tracking (DIT) [10,11,12,13]. The PF is likely the strongest of these methods, because its field variables can directly represent any complex morphology, thus eliminating the need to identify and track the interface. However, it is very computationally taxing even when combined with adaptive meshing methods. The LS method also interchanges the interface with a field variable (the level set), but requires knowledge of the direction in which the solid front is advancing, its velocity, and the calculation of the vector normal to the interface; which makes it less effective than the PF method in complex three-dimensional geometries. DIT is the simplest and computationally the most efficient of the methods, but it is also less powerful because it requires the calculation of the temperature gradients at the interface in addition to the normal velocity and curvature of the interface. Moreover, the complexity for handling interfaces in all possible solidification conditions limits the applicability of DIT methods.

Despite the extensive work that has been done in two dimensions, fewer three-dimensional calculations of dendritic solidification at the microscale have been reported [14,15,16,17,18,19,20,21,22,23]. Some of these works include convection effects, but only the recent publication [14,19,20,23] consider binary alloys. Regardless of the method of choice, the computational effort is significantly larger than in two-dimensions, this is due in part to the fact that the growth velocities are higher in three-dimensions, and also to the presence of the capillary length scale, that cannot be fully resolved.

For the reasons explained above, it is apparent that none of the PF, LS, or DIT methods is very well suited for three-dimensional calculations of dendritic growth. On the other hand, cellular automaton (CA) based methods have been shown to offer an alternative [24,25,26,27,28,29,30,31,32] allowing to achieve reasonable accuracy for calculations in bigger domains. These techniques do not capture the same level of detail at the solid-liquid interface as methods based on PF or DIT. In particular, initial CA models could not account for the effect of undercooling on the growth velocity except ahead of the dendrite tips. However, in recent years they have been combined with finite element and/or finite difference/volume solutions of the energy and solute transport equations and developed to a point where solute concentration, undercooling, and latent heat dissipation can be properly modeled [28,29]; crystal anisotropy effects have also been incorporated.

Cellular automaton (CA) models are attractive due to being simple in their construction and yet capable to represent complex phenomena. Therefore, CA models have been deployed in numerical simulations to characterize the development of microstructure throughout the process of solidification.

One of the most elaborate early CA models for calculation of grain growth, combining stochastic nucleation, diffusional growth and macroscopic heat conduction was proposed by Gandin and Rappaz [24,25] as an extension to the stand-alone CA model. Sanchez and Stefanescu [27] and Zhu and Stefanescu [31] proposed a model based on the cellular automaton technique for the simulation of dendritic growth controlled by solutal effects in the low Peclet number regime. One of the innovative aspects of this model is that it does not use an analytical solution to determine the velocity of the solid-liquid interface as is common in other models.

The lattice-Boltzmann method (LBM) is a numerical tool suitable for the simulation of multiphase flows subject to complex time-evolving boundary conditions. Contrary to classic computational fluid dynamics (CFD) methods that numerically represent the continuum-based Navier-Stokes equations, the LBM builds on a statistical model that incorporates micromechanical characteristics into mesoscale. The flow representation assumes the fluid to consist of fictitious particles on a lattice. These particles displace to the neighboring lattice points, and collide with particles arriving from other directions. The collisions are characterized by the rate of relaxation to equilibrium flow configuration. When compared with traditional CFD methods, the LBM offers several advantages, such as simple physical interpretation, strictly local calculations suitable for large-scale parallelization, straightforward dealing with complex geometries and time evolving boundary conditions, reasonable accuracy, and an explicit algorithm with regular grid convenient for coupling across multiple processes with varying spatial and temporal scales.

In addition to modeling of basic fluid flow, the LBM has been recently applied to model multi-component and multiphase flows, non-Newtonian flow with shear thickening, heat flow with phase transformation, chemical reactions, and other engineering applications. Jiaung et al. [33] introduced an extended lattice Boltzmann equation governed by the heat conduction equation in conjunction with enthalpy method. A new variant of LBM for heat conduction problem with phase change was presented by Eshraghi and Felicelli [34]. Unlike earlier explicit algorithms, the latent heat term was represented implicitly in the energy equation, improving the stability and efficiency of the solution. Raj et al. [35] combined LBM with the discrete transfer method to study solidification with radiation in a semitransparent material. Medvedev and Kassner [36] and Miller et al. [37] used LBM in combination with the phase-field method to simulate dendritic growth and alloy solidification in the presence of fluid flow. In one of the few known applications combining LBM with a classical cellular automaton, Chopard and Dupuis [38] simulated the erosion and deposition process in water channels.

Cellular automaton (CA) based methods, utilizing explicit local algorithms on regular grids, can be easily incorporated into LBM framework, as it was done in Sun’s et al. [39] 2D LBM-CA simulations of dendrite evolution. In their work, CA was deployed to capture the solid-liquid interface while LBM solved the mass and solute transport. Yin et al. [40] compared their LB-CA solidification model with an FE-CA model. Nevertheless, 2D models assume an in-depth homogeneity of the system, a condition that prohibits modeling of fully-realistic 3D microstructural features. For example, the melt flow can be completely blocked by the dendrite arms in two-dimensional simulations, while it can bypass the 3D arms, producing different distribution of solute, growth kinetics, and dendrite morphology. The dendritic growth in 2D can significantly differ from the 3D growth, as demonstrated in studies [14,15,16,19,23]. A three-dimensional LB-CA model for simulating dendrite growth under convection was developed by Eshraghi et al. [41]. However, the model did not include fluid flow calculations. Their recent work [42], which was obtained by the mentioned 3D LB-CA model, presents large-scale simulation results for columnar dendrite growth in a 1 mm³ region. The model showed an excellent speed-up and scale-up performance on supercomputers.

This paper presents a three-dimensional LB-CA model for simulation of dendrite growth under forced and natural convection. The model is verified against available analytical solutions for diffusion-advection, fluid flow and dendrite growth. The results of 3D and 2D simulations are compared. In addition, the change in growth kinetics and morphology of Al-Cu dendrites is studied by altering melt undercooling, alloy composition, and inlet flow velocity.

2. Model Description

In the LB model developed in this work, the time domain is divided into equal time steps while the spatial domain is discretized with a uniform 3D cubic grid. The present work deploys commonly used D3Q15 lattice, which is a three-dimensional lattice with fifteen discrete velocity vectors; six of them point towards the faces of the cell under consideration, eight towards the cell corners, and one with zero magnitude resides in the cell center. The coordinates of the D3Q15 lattice velocities are:

$$e_i=c \text{×} \{\begin{array}{c}\left(0, 0, 0\right) \mathrm{for} i=0\\ (\pm 1, 0, 0), (0, \pm 1, 0), (0, 0, \pm 1) \mathrm{for} i=1\text{–}6\\ (\pm 1, \pm 1, \pm 1) \mathrm{for} i=7\text{–}14\end{array}$$

(1)

where $c=∆x/∆t$ is the lattice speed, $∆x$ is lattice spacing, and $∆t$ is time step size.

It is assumed that the temperature is constant and the undercooling does not change throughout the simulation. Therefore, the heat transport is absent—only the fluid flow and solute diffusion are considered. Also, solute transport in the solid phase is neglected.

The motion of an incompressible Newtonian fluid is described by the Navier-Stokes equations (NSE):

$$\rho \left(\frac{\partial u}{\partial t}+u.\nabla u\right)=-\nabla p+\nabla .\left(\mu \nabla u\right)+F$$

(2)

where ρ and u are density and velocity, respectively. $F$ is the source term for the natural buoyancy force that can be defined as:

$$F=-{\rho }_{0}g{\mathsf{\beta }}_{\mathrm{C}}\left(C_l-C_{\mathrm{ref}}\right)$$

(3)

where $g$ is the gravitational vector, ${\rho }_0$ is the fluid density at a reference temperature and composition $C_{\mathrm{ref}}$. ${\mathsf{\beta }}_{\mathrm{C}}$ is the expansion coefficients for composition, and $C_{\mathrm{l}}$ is the solute concentration. The direction of the gravitational vector is towards the bottom of the domain. The strength of the buoyancy force is evaluated by Rayleigh (Ra) number, which depends on the concentration gradient and can be defined as:

$$R{\mathrm{a}}_C=\frac{g{\mathsf{\beta }}_C\left(C_{\mathrm{l}}-C_{\mathrm{ref}}\right)L^3}{\upsilon D_{\mathrm{l}}}$$

(4)

where υ is kinematic viscosity and L is the characteristic length in the direction parallel to the gravitational vector. $D_{\mathrm{l}}$ is the diffusivity of solute in the melt. The characteristic length, L, represents the depth of the mushy zone and is set to be equal to the height of the domain. Ra number represents the maximum force due to the difference between the local maximum and the reference solute concentration.

According to the BGK (Bhatnagar-Gross-Krook) approximation [43], the LB equation for the evolution of the fluid density distribution functions can be written as:

$$f_i\left(x+e_i∆t, t+∆t\right)-f_i\left(x, t\right)=-\frac{1}{{\mathsf{\tau }}_u}\left(f_i\left(x,t\right)-f_i^{\mathrm{eq}}\left(x,t\right)\right)+∆t{F}_i\left(x,t\right)$$

(5)

where $e_i$, and ${\mathsf{\tau }}_u$ are the discrete velocities, and the non-dimensional relaxation time, respectively. $f_i$ is the distribution function that is defined as the probability of finding a particle at position x moving in direction i and $F_i\left(x,t\right)$ is the force term produced by natural buoyancy:

$$F_i=-3w_i\rho e_{i}F/c^2$$

(6)

$f_i^{\mathrm{eq}}$ is the equilibrium distribution function that is defined as:

$$f_i^{\mathrm{eq}}=w_i\rho \left(1+3\frac{e_i.u}{c^2}+\frac{9{\left(e_i.u\right)}^2}{2c^4}-\frac{3u^2}{2c^2}\right)$$

(7)

$w_i$ is the weight coefficient that depends on the type of lattice and direction. For the D3Q15 lattice used in this work, the weight coefficients are given as:

$$w_i=\{\begin{array}{c}16/72 \mathrm{for} i=0,\\ 8/72 \mathrm{for} i=1–6,\\ 1/72 \mathrm{for} i=7–14\end{array}$$

(8)

Using the Chapman-Enskog expansion, it can be shown that the LB model can approximate NSE if the following equation is satisfied.

$$\upsilon =\frac{c^2∆t}{6}\left(2{\mathsf{\tau }}_u-1\right)$$

(9)

Then, the macroscopic density and velocity of fluid can be obtained from the distribution functions:

$$\rho ={\displaystyle \sum }_{i=0}^{14}{f}_i$$

(10)

$$u=\frac{1}{\rho }{\displaystyle \sum }_{i=0}^{14}{f}_i{e}_i$$

(11)

The solute within the liquid phase can be transported by diffusion and by advection. The advection-diffusion equation governs the solute transport in presence of diffusion and melt flow:

$$\frac{\partial C_{\mathrm{l}}}{\partial t}=\nabla .\left(D_{\mathrm{l}}\nabla C_{\mathrm{l}}\right)-u\nabla C_{\mathrm{l}}$$

(12)

The equivalent LB equation for solute transport can be written as:

$$g_i\left(x+e_i∆t, t+∆t\right)-g_i\left(x, t\right)=-\frac{1}{{\mathsf{\tau }}_c}\left(g_i\left(x,t\right)-g_i^{\mathrm{eq}}\left(x,t\right)\right)$$

(13)

where $g_i$ is the distribution function for solute concentration and ${\mathsf{\tau }}_c$ is the relaxation time for the solute transport model. The equilibrium distribution function for solute concentration, $g_i^{\mathrm{eq}}$, is given by:

$$g_i^{\mathrm{eq}}=w_i{C}_{\mathrm{l}}\left(1+3\frac{e_i.u}{c^2}\right)$$

(14)

Then, the macroscopic solute concentration, C_l (x, t), can be obtained as:

$$C_{\mathrm{l}}={\displaystyle \sum }_{i=0}^{14}{g}_i$$

(15)

The LB Equation (13) approximates the advection-diffusion Equation (12) if:

$$D_{\mathrm{l}}=\frac{c^2∆t}{6}\left(2{\tau }_{\mathrm{c}}-1\right)$$

(16)

In the present model, the difference between local equilibrium and actual solute concentrations in the liquid determine the kinetics of dendrite growth. The change in the fraction of solid, $∆f_s$ in an interface cell is calculated according to the equilibrium condition at the interface [31]:

$$∆f_s=\left(C_{\mathrm{l}}^{\mathrm{eq}}-C_{\mathrm{l}}\right)/\left(C_{\mathrm{l}}^{\mathrm{eq}}\left(1-k\right)\right)$$

(17)

where k is the solute partition coefficient from the phase diagram, and $C_{\mathrm{l}}$ is the concentration of solute in the liquid phase given by Equation (15). Moreover, $C_{\mathrm{l}}^{\mathrm{eq}}$ is the equilibrium concentration at the interface according to [44]:

$$C_{\mathrm{l}}^{\mathrm{eq}}=C_0+\frac{T_{\mathrm{l}}^{\ast }-T_{\mathrm{l}}^{\mathrm{eq}}+\mathsf{\Gamma } \mathrm{wmc}}{m}$$

(18)

where $T_{\mathrm{l}}^{\ast }$ is the temperature at the interface, $T_{\mathrm{l}}^{\mathrm{eq}}$ is the equilibrium liquidus temperature at the initial solute concentration ($C_0$), 𝑚 is the slope of the liquidus, and Γ is the Gibbs-Thomson coefficient. The weighted mean curvature (wmc) incorporates the effect of anisotropic surface energy [45] as follows:

$$\begin{gathered} \mathrm{wmc}=\left(3\mathsf{\epsilon }-1\right)\left({\partial }_x{n}_x+{\partial }_y{n}_y+{\partial }_z{n}_z\right)-48\mathsf{\epsilon }\left(n_x^2{\partial }_x{n}_x+n_y^2{\partial }_y{n}_y+n_z^2{\partial }_z{n}_z\right)+ \\ 12\mathsf{\epsilon }Q\left({\partial }_x{n}_x+{\partial }_y{n}_y+{\partial }_z{n}_z\right)+12\mathsf{\epsilon }\left(n_x{\partial }_{x}Q+n_y{\partial }_{y}Q+n_z{\partial }_{z}Q\right) \end{gathered}$$

(19)

where $\mathsf{\epsilon }$ is the anisotropy of the surface energy. Considering $\widehat{n}$ as the unit vector normal to the interface, its components can be evaluated as $n_x={\partial }_x{f}_s/\left|\nabla f_s\right|$, $n_y={\partial }_y{f}_s/\left|\nabla f_s\right|$, and $n_z={\partial }_z{f}_s/\left|\nabla f_s\right|$, where $\left|\nabla f_s\right|=\sqrt{{\left({\partial }_x{f}_s\right)}^2+{\left({\partial }_y{f}_s\right)}^2+{\left({\partial }_z{f}_s\right)}^2}$. The parameter $Q$ is defined as $Q=n_x^4+n_y^4+n_z^4$ [44].

A CA algorithm is used to capture new interface cells. The CA mesh is identical to the LB one. Each CA cell is marked as solid, liquid, or as interface, and equipped with the attributes of the solid fraction and crystallographic orientation. The state of each cell at each time step is updated according to the state of itself and its neighbors. Following the Moore’s neighborhood rule, the first nearest neighbors of a completely solidified cell are marked as new interface cells. The final solid/liquid interface is then determined by changing the fraction of solid according to $∆f_s$ from Equation (17) as follows:

$$f_s^n=f_s^{n-1}+∆f_s$$

(20)

During solidification, the solute is redistributed between solid and liquid phases. The solute amount determined as:

$$∆C_{\mathrm{l}}=C_{\mathrm{l}}\left(1-k\right)∆f_s$$

(21)

is rejected from the solidifying cells and then evenly redistributed among the surrounding liquid and interface cells.

The physical properties for the alloy material considered in the present simulations are shown in

Table 1. Physical properties of Al-Cu alloy used in the simulations.

Density, ρ (kg m−3)	Diffusion Coefficient, D (m2 s−1)	Viscosity, μ (N s m−2)	Liquidus Slope, m (°C wt %−1)	Partition Coefficient, k	Gibbs-Thomson Coefficient, Г (m °C)	Degree of Anisotropy, ε
2475.0	3.0 × 10−9	0.0024	−2.6	0.17	2.4 × 10−7	0.04

. The same parameters are employed for all simulations, unless otherwise mentioned.

The dendrite growth is initialized by placing nucleation sites with a chosen crystallographic orientation into the undercooled melt. The solute is then rejected from the solid to the interface cells and the local liquid composition increases to reach the equilibrium composition.

For the case of forced convection, a uniform flow with velocity in the x-direction enters the domain from the west boundary, as shown in Figure 1. The non-slip bounce-back boundary condition is applied at the solid/liquid interface for both fluid flow and solute diffusion calculations. In bounce back, the incoming distribution functions at the solid particles are reflected back to the fluid. Using this simple method, the interaction of fluid with complex boundaries of the dendrite can be efficiently modeled. The side walls are assumed to be periodic and a zeroth extrapolation is used to simulate the open wall on the east side. Besides, all boundary walls are insulated against the solute transport. The dendrite remains stationary in the center of the simulation domain.

For the case of natural convection, the solid seed is located at the same place as the forced convection case. Also, the boundary conditions for concentration profile are similar to the forced convection case. For the flow field, periodic boundary condition is applied for all walls. The buoyancy force induced by solute concentration gradient is added to the momentum conservation equation by Rayleigh number (Ra).

The non-dimensional relaxation parameter is chosen to be 1 for both fluid flow and solute transport models and the mesh spacing is 0.3 μm for all dendrite growth simulations.

3. Results and Discussion

3.1. Validation

Two benchmark problems are implemented to validate the LBM fluid flow and solute transport model by comparison with analytical solutions. The dendrite growth model is then validated by comparing the LBM results with the Lipton-Glicksman-Kurz (LGK) [46,47] dendrite growth characterization.

3.1.1. Fluid Flow

A benchmark problem of a 3D laminar flow around an off-centered cylinder inside a rectangular channel is implemented within the present LBM framework. This benchmark has been presented by Schäfer and Turek [48] together with results using a variety of numerical methods. The geometry, initial conditions, and boundary conditions employed here are identical to those in the reference.

Figure 2 shows the configuration and boundary conditions for flow around a circular cylinder. The inflow condition is $U\left(0,y,z\right)=16U_m yz\left(H-y\right)\left(H-z\right)/H^4 , V=W=0$ with U_m = 0.45 m/s. For no slip boundary condition on the walls and the cylinder, the bounce-back rule is applied. The inflow condition is imposed by specifying the given velocity profile at the entrance, and a zeroth-order extrapolation for the distribution function is imposed at the exit. Density and viscosity are assumed to be 1 kg/m³, and 0.001 m²/s, respectively. The method suggested by Mei et al. [49] is used to calculate drag and lift forces, and consequently drag and lift coefficients, as well as pressure drop.

Table 2. Fluid characteristic quantities for steady state laminar flow over a cylinder. Simulation results are listed for three different grid sizes. Lower and upper bounds present the estimated interval for the exact results [48].

Grid Spacing (m)	Cd	Cl	ΔP
0.005	6.6315	0.0225	0.1730
0.0025	6.3503	0.0148	0.1702
0.00166	6.2550	0.0082	0.1681
Lower Bound	6.0500	0.0080	0.1650
Upper Bound	6.2500	0.0100	0.1750

shows the LBM results for three different grid spacings. The last two rows present the intervals of the exact solution as reported by Schäfer and Turek [48]. It is observed that the error in the drag and lift coefficients is larger for coarser grids. The reason is that the discretization of the nodes using a rectangular grid cannot precisely reproduce the circular outline of the cylinder. However, the pressure drop across the cylinder is within the bounds for all three grid sizes.

3.1.2. Solute Transport

For the case of an instantaneous volume source of mass over the range −∞ < x < 0 with concentration C_i, and for a steady velocity field u = (U₀, 0, 0), the solution of the diffusion-advection equation is given by [50]:

$$\frac{C\left(x,t\right)}{C_i}=\frac{1}{2}\left(1-\mathrm{erf}\left(\frac{x-U_{0}t}{\sqrt{4Dt}}\right)\right)$$

(22)

This case is examined for a 3D channel with a square cross section. A one-dimensional flow of U₀ = 10 mm/s velocity enters the domain from the negative x-direction and exits through the positive x-boundary. Periodic boundary conditions are applied on the y- and z-boundaries. The concentration of the source is set to C_i = 4 wt %. Figure 3 shows the concentration profile in the x-direction after 2 ms. The channel length and the analysis time are chosen so that the concentration remains unchanged in the regions close to right and left boundaries and the infinity assumption is satisfied. A good agreement is found between LB and analytical solutions.

3.1.3. Dendrite Growth

The present dendrite growth model is validated by comparison with the LGK analytical approximation [46,47]. The present model is able to predict the steady state dendrite’s tip velocity, radius, growth Peclet number, and interface concentration at the tip as a function of melt undercooling, diffusion coefficient, initial concentration, liquidus slope, partition coefficient, Gibbs-Thomson coefficient and a selection parameter, σ*, that depends on the anisotropy of surface energy. According to 3D linearized solvability theory [51], the σ* value of 0.085 for ε = 0.03 is calculated and applied for verifying the dendrite growth model. A solid seed is initially placed at the center of a cubic domain. There is no fluid flow or thermal transport in the melt and dendrite growth is driven solely by undercooling with solute redistribution. The tip properties are measured after the initial transient stage, but before the side branches start to develop, as the LGK model assumes a parabolic needle tip.

The interface position is determined by the location of the interface cell and the fraction of solid in that cell. Then, the interface velocity is calculated using a simple difference formula.

Figure 4a compares the tip velocities simulated by the LB-CA model and the tip velocities from the LGK approximation. The tip velocities presented in this paper are obtained as an average of the interface velocities of the tips of six primary dendrite arms over a time interval when growth is in the steady state. The interface velocity increases with undercooling, as expected. Although small differences are present, the LB-CA values are in a good agreement with the LGK approximation.

The tip radius of the dendrites is measured from 3D fraction-of-solid contour plots at f_s = 0.5. These 3D dendrite tips were sectioned by two planes passing through the center of the dendrite to obtain 2D fin shapes. The fin shapes were fitted by fourth-order polynomials, and the curvature at the dendrite’s tip was obtained as K = d²y/dx²(1 + (dy/dx)²)^−3/2 [52]. Then, the tip radius was evaluated as R = 1/K. The same measurement was taken for six different branches at different simulation times after the dendrite growth reached the steady state. The plot in Figure 4b shows the tip radii from present simulations along the LGK estimates. A larger difference between LGK and simulated values can be noticed in the tip radius especially at low undercooling. Nonetheless, the simulation results are correct in the same order of magnitude and exhibit the same trend as LGK approximation. The primary cause of the discrepancy in the tip radius is the dendrite tip discretization in the CA method, which results in an inaccurate representation of the tip and introduces the error in calculated curvature and radius. As undercooling decreases, the tip radius increases and the tip curvature becomes smaller and more erroneous. Also, at small undercoolings, the dendritic shape deviates from the parabolic shape assumed in the LGK model. A similar issue has been reported in [44,53].

The solutal Peclet number, P_c, is evaluated as P_c = VR/(2D_l), where V and R are the tip velocity and radius, respectively. The simulated Peclet numbers along with LGK predictions for a range of undercoolings are plotted in Figure 4c, demonstrating a good match.

Furthermore, the interface equilibrium liquid composition, $C_{\mathrm{l}}^{\mathrm{eq}}$, is compared with the LGK approximation in Figure 4d. The plotted values are obtained by averaging the composition of six different tips over several time steps. Since the temperature is constant in the domain, according to Equation (15), the equilibrium composition is only affected by the tip’s weighted mean curvature, wmc. A reasonable match between the analytical and simulation results was found.

3.2. Dendrite Growth under Melt Convection

In this section, the simulation results for dendritic growth under forced and natural melt convection are discussed. The growth characteristics are compared for 2D and 3D simulations, and the effects of altering melt undercooling, inlet flow velocity, and alloy composition on the kinetics of 3D dendrite growth is explored.

3.2.1. Kinetics of Growth under Forced Convection

Figure 5 depicts the growth of a single dendrite in the presence of the melt flow. The snapshots taken after 2 ms, 4 ms, 6 ms, and 8 ms of simulated time. The melt enters the domain from the left boundary with the velocity U₀ = 7 mm/s. The melt flow induces the solute convection that affects solute distribution around the dendrites and therefore influences the kinetics of dendritic growth. The cube-shaped simulation domain contains 2883 cubic cells with the grid spacing ∆x = 0.3 μm. The domain temperature is set to 4.5 °C below the melting point (the undercooling ∆T = 4.5 °C). The streamlines show how the melt flows around the developing dendrite. In addition, the 2D sections of the dendrite morphology perpendicular to z-direction passing through the dendrite center are shown in Figure 6. The contours, from inside to outside, show the morphologies after 0.75 ms, 2.00 ms, 3.25 ms, 4.50 ms, 5.75 ms, 7.00 ms, 8.25 ms, and 9.00 ms. At the early stages of solidification, the primary arms grow along their preferred crystallographic. As solidification progresses, the primary arms grow and coarsen, and then, the secondary arms appear perpendicular to the primary arms. Melt convection washes the solute from upstream primary and secondary arms and carries it downstream. This leads to a lower concentration in the upstream area and a higher concentration downstream and, consequently, a higher growth rate upstream and a lower growth rate downstream. This matches the findings of previous studies [14,18,19,23,39,40]. The transverse arms (the primary arms in the y- and z-directions) are not significantly affected by convection, only the secondary branches grow slightly faster on the upstream side of the transverse arms.

Figure 7 shows the velocity of the melt flow around the dendrite. The velocity is higher in the areas away from the dendrite and lower in the dendrite’s vicinity. The bottom portion of Figure 7 displays the solute concentration in two perpendicular sections through the center of the dendrite. As the solubility of solute in solid is lower than its solubility in liquid, an increasing amount of solute is rejected to the melt as the solidification progresses. At later times, when solidification is nearly complete, there are liquid regions enclosed by the solid phase, holding a high solute concentration. These regions may end up causing micro-segregation or forming eutectic phases.

3.2.2. Comparison of 2D and 3D Simulations

As mentioned earlier, studies have shown that the kinetics of dendrite growth differs considerably in two and three dimensions. Even without convection, the diffusion in three dimensions is more effective than the diffusion in two dimensions. Therefore, solute concentration and its gradient are higher close to 2D dendrites, slowing down the 2D growth. Figure 8 compares the solute distributions in 2D and 3D domains. The results are captured after the simulated time of 7.5 ms for both cases. 2D simulations are performed with a similar methodology as 3D, but using the D2Q9 lattice. Interested readers are referred to Reference [40] for more details on the 2D model. Melt flow around the 3D dendrite produces a thinner solute layer than the flow around 2D dendrite. Therefore, the tip growth velocity is higher in 3D, the observation that was also confirmed by other studies [23]. This effect restrains the side branches from forming in 2D, but promotes their growth in 3D, especially the upstream growth. Comparison of 2D and 3D sections reveals that the center of the 3D dendrite in presence of convection is noticeably displaced in the upstream direction, something that is not detected in 2D, as 2D flow cannot sweep the solute downstream as effectively as 3D flow.

Figure 9 depicts the x-component of velocity around 2D and 3D dendrites. Note that the color ranges differ for 2D and 3D cases. While the melt can pass between the side arms of the 3D dendrite, the flow cross-section is mostly blocked by the side arms in 2D and the melt cannot pass around the side arms. Due to smaller passage for the flow in 2D, the maximum velocity is much higher in the 2D case than in the 3D case. It can be seen that the downstream side of the 2D dendrite faces much less convection than the 3D one, resulting in more solute accumulation and less growth. There is no apparent recirculation downstream of the 3D dendrite because liquid can flow more freely while recirculation past the dendrite is very likely in 2D. These effects altogether amplify the size difference between the arms growing in upstream and downstream directions in a 2D domain. The length ratio of the upstream arm to the downstream arm in Figure 8 is 2.75 for the 2D case, while it is 1.43 for the 3D case. The results of the simulations indicate that the length ratio increases with time for both cases.

3.2.3. Effect of Melt Flow and Undercooling Strength

The undercooling effects are examined for Al-3 wt %Cu alloy with the inlet melt velocity of U₀ = 7 mm/s. Figure 10 shows the variation of tip velocity with undercooling. The growth rate and tip velocity of all arms increase with stronger undercooling. Tip growth velocity of the upstream arm is the highest whereas the tip velocity of the downstream arm is lowest. The tip velocity of the transverse arms is in between the velocities of the upstream and downstream arms. The melt flow washes the solute from the upstream tip, reducing the local solute concentration in front of the tip. Therefore the difference between local and equilibrium solute concentrations increases and leads to a faster growth. On the other hand, melt convection causes an increase in the solute concentration in front of the downstream tip and consequently decreases the growth rate. Fluid flow has no substantial impact on the solute concentration at the tip of the transverse arms because the fluid flow does not produce a solutal gradient in the transverse tip growth direction as it does along the upstream arm tip. It can be seen that the average tip velocity of the transverse arms is almost same as the tip velocity for the case without melt flow.

Figure 11 shows 2D sections of the 3D dendrite growing in melts with different undercooling strength. The bottom and top halves present morphologies of the dendrite growing in the melt with 3 °C and 4 °C of undercooling, respectively. Different colors represent different simulated times. As expected, the dendrite in the melt with higher undercooling, ∆T = 4 °C, grows faster than the dendrite in the melt with lower undercooling, ∆T = 3 °C. However, the impact of melt convection is more significant for the dendrite with lower undercooling, because convection has more time to accumulate solute downstream. The difference between the lengths of upstream and downstream branches is larger for the lower undercooling. The difference increases with time, as the downstream arm grows very slow for ∆T = 3 °C. As a result of faster growth, secondary branches also form earlier in the melt with larger undercooling.

3.2.4. Effect of Inlet Velocity

Variation of tip growth velocities with the velocity of melt flow is shown in Figure 12. The solute concentration and undercooling are set to 3 wt % and 4.5 °C, respectively. Increasing flow velocity slightly increases the growth velocity of the upstream arm. As mentioned earlier, fluid flow does not seem to have a significant effect on the transverse arms. Even for high inlet velocities, the transverse arms don’t seem to be much affected. Increasing flow velocity shows a more significant effect on the tip velocity of the downstream arm. By increasing the flow velocity, more solute is accumulated around the downstream arm, increasing the local liquid composition. According to Equation (17), the change in the fraction of solid is reduced as the local liquid composition increases, which results in a lower growth rate. A tip splitting phenomenon is observed for inlet velocities higher than 10 mm/s, leading to inaccurate and unstable measurement of the tip velocities.

Figure 13 illustrates the morphologies of the dendrites growing under different melt flow velocities. As the melt velocity increases, the dendrite grows faster in the direction opposite to the flow direction. Increasing flow velocity does not show a significant impact on the length of the upstream arm. Interestingly, the flow washing away the solute layer at the interface causes enhanced solute growth on all the upstream faces of the dendrite, producing a net displacement of the dendrite center in the upstream direction. This displacement increases with the increasing magnitude of the flow velocity. It can be seen that the distance from the center of the simulation domain to the dendrite’s tip is slightly longer, but there is not a large increase in the distance from the center of dendrite to the upstream tip between the dendrites growing under different flow velocities. As mentioned earlier, increasing the flow velocity has a more significant effect on the downstream arm. Flow velocity slightly decreases the tip velocity of the transverse arms, but the effect is not pronounced. However, as inlet velocity increases, growth of secondary arms on the transverse arms is promoted in the upstream direction.

3.2.5. Effect of Alloy Composition

Effect of alloy composition on the growth kinetics is displayed in Figure 14. The undercooling and melt velocity at the inlet are set to 4.5 °C and 7 mm/s, respectively. In general, the growth rate decelerates with increasing concentration of the solute. As the alloy’s solute content increases, all upstream, downstream, and transverse tip growth velocities decrease at a comparable rate. This is an expected trend due to the constant partition coefficient and liquidus slope used in the simulations.

Figure 15 depicts 2D wireframes of the 3D dendrites growing in melts with different compositions. The bottom and top sections represent Al-3 wt % Cu and Al-9 wt % Cu alloys, respectively and the colors indicate different simulated times. Note that since the growth rates are significantly different, showing the outlines at identical times for both compositions may not illustrate the detail of dendritic morphologies. The only time at which the morphology is displayed for both compositions is t = 3.5 ms, which is shown in red color in Figure 15. As can be seen, the growth speed is considerably faster in the alloy with lower solute concentration. In addition, the size ratio of the upstream arm to the downstream arm is significantly larger in the alloy with higher solute concentration. Because the dendrite grows very slowly in Al-9 wt % Cu alloy, convection can effectively wash away the solute from upstream branch tip and accumulate it downstream. Hence, the local liquid concentration stays always high in front of the downstream arm, hindering the arm’s advancement and increasing the length ratio between upstream and downstream arms.

3.2.6. Effect of Natural Convection

The evolution of an equiaxed dendrite under the effect of natural convection is examined. The calculation domain consisted of a 200 × 200 × 200 grid which is equivalent to a 60 μm × 60 μm × 60 μm domain. Figure 16 and Figure 17 present the simulated morphology of the Al-3 wt % Cu dendrite in an undercooled melt with ∆T = 4.5 Κ under the effect of natural convection. Simulation time is 6.4 ms for the case presented in the figures. Figure 16 shows the liquid motion induced by buoyancy force in the vicinity of the dendritic structures, while Figure 17 depicts the morphological changes of the 3D dendrite growing under natural convection. The Rayleigh number is taken as $R{\mathrm{a}}_C$ = 85 × 10³. Solute gradient due to natural convection forces leads to the onset of vertical flows in the upward direction, as can be seen in Figure 16. The buoyancy-induced flow washes the solute from the interface of the downward arm increasing its growth rate. On the other hand, the vertical motion of flow accumulates the solute around upward arm and retard its growth based on the same mechanism (Figure 17).

4. Conclusions

A three-dimensional lattice Boltzmann model for simulating dendrite growth under forced and natural convection is introduced. The fluid flow, solute diffusion, and dendrite growth models are validated against analytical solutions. The results show that growth kinetics and dendritic morphology are significantly affected by the presence of melt flow. In the presence of melt flow, the primary and secondary dendrite arms grow faster in the upstream direction, and significantly slower in the downstream direction. The melt flow does not have a noticeable impact on the dendrite growth in the direction transverse to the flow; the solidification rate in the transverse directions is similar to the case without melt flow. A comparison with 2D results revealed that 3D dendrites grow faster in general, but the difference between the growth velocities of the upstream and downstream arms is more significant in 2D simulations. Moreover, 3D dendrites are found to be more likely to form secondary arms. It was also revealed that under the influence of melt convection the center of the 3D dendrite will slightly displace in the upstream direction. The displacement grows with increasing inlet melt flow velocity. Increasing the degree of undercooling accelerates the growth rate in all directions. The convection effects are intensified by increasing the magnitude of flow velocity. The flow impact on the downstream arm seems more significant than on the upstream arm. The growth rate slows down in all branches when the solute concentration in the alloy is higher. The size ratio of the upstream arm to the downstream arm grows by increasing inlet velocity and solute content, and by decreasing undercooling. When natural convection is present, the buoyancy-induced flow suppresses the growth of the upward arm and accelerates the evolution of the other arms, especially the downward arm.

Considering the special capabilities of the lattice Boltzmann method, e.g., simple implementation, computational efficiency, local calculations, and inherent parallel structure, the model presented here provides a great tool for simulating the solidification of 3D dendritic structures under convection.