Journal of Materials Science 【-逻*辑*与-】amp; Technology, 2020, 49(0): 15-24 doi: 10.1016/j.jmst.2020.01.047

Research Article

## Numerical simulation for dendrite growth in directional solidification using LBM-CA (cellular automata) coupled method

Wonjoo Leea, Yuhyeong Jeonga, Jae-Wook Leeb, Howon Leeb, Seong-hoon Kangb, Young-Min Kimb, Jonghun Yoon,c,*

a Department of Mechanical Design Engineering, Hanyang University, 222, Wangsimni-ro, Seongdong-gu, Seoul 04763, Republic of Korea

b Materials Deformation Department, Korea Institute of Materials Science, 797 Changwondaero, Seongsan-gu, Changwon-si, Gyeongnam-do 51508,Republic of Korea

c Department of Mechanical Engineering, Hanyang University, 55 Hanyangdaehak-ro, Sangnok-gu, Ansan-si, Gyeonggi-do 15588, Republic of Korea

Corresponding authors: * E-mail address:yooncsmd@gmail.com(J. Yoon).

Received: 2019-10-24   Revised: 2019-12-11   Accepted: 2020-01-15   Online: 2020-07-15

Abstract

To predict the dendrite morphology and microstructure evolution in the solidification of molten metal, numerically, lattice Boltzmann method (LBM) - cellular automata (CA) model has been developed by integrating the LBM to solve the mass transport by diffusion and convection during solidification and the CA to determine the phase transformation with respect to the solid fraction based on the local equilibrium theory. It is successfully validated with analytic solutions such as Lipton-Glicksman-Kurz (LGK) model in static melt, and Oseen-Ivantsov solution under the fluid flow conditions in terms of tip radius and velocity of the dendrite growth. The proposed LBM-CA model does not only describe different types of dendrite formations with respect to various solidification conditions such as temperature gradient and growth rate, but also predict the primary dendrite arm spacing (PDAS) and the secondary dendrite arm spacing (SDAS), quantitatively, in directional solidification (DS) experiment with Ni-based superalloy.

Keywords： Cellular automata (CA) ; Lattice Boltzmann method (LBM) ; Dendritic growth ; Directional solidification

Export EndNote| Ris| Bibtex

Wonjoo Lee, Yuhyeong Jeong, Jae-Wook Lee, Howon Lee, Seong-hoon Kang, Young-Min Kim, Jonghun Yoon. Numerical simulation for dendrite growth in directional solidification using LBM-CA (cellular automata) coupled method. Journal of Materials Science & Technology[J], 2020, 49(0): 15-24 doi:10.1016/j.jmst.2020.01.047

## 1. Introduction

Conventional casting process utilizes the solidification of the molten metal to manufacture a product having extremely complex shape with high productivity. The mechanical properties of the as-cast product are not only substantially influenced by the process parameters such as the cooling rate, temperature, and target shape, but also determined by the microstructure induced by the thermal and solute constraints [[1], [2], [3], [4]]. Especially, the dendrite morphology plays an important role in exhibiting the mechanical properties of the as-cast product since the distance of the dendrite arm is closely related to the ductility and the strength of the as-cast alloy [[5], [6], [7]]. Likewise, it is essential to control and predict the microstructure of the final product in the other processes such as the welding [[8], [9], [10], [11]] and 3D printing [[12], [13], [14]] undergoing the melting and solidification. Unfortunately, the dynamic evolution of microstructure is hardly observed, experimentally, in these processes since it is highly sensitive to the rate of heat extraction and the heat flow in the solid-melt interface during solidification. Therefore, there are extensive numerical approaches using the phenomenological [15,16], analytical [17,18], and probabilistic [[19], [20], [21]] models to take into consideration of the nucleation and grain growth based on the physical mechanism during the melting and solidification processes. The phase field (PF) [[22], [23], [24], [25]] and cellular automata (CA) [[26], [27], [28], [29]] models are the most widely used approaches to demonstrate the evolution of the dendrite formation during the solidification process since they are able to keep track of dynamic motion of interface between the solid and liquid phases based on the solute and thermal distributions. While the interface is expressed with a smooth change in the phase field variable in the PF model, continuously, it consists of a single layer of discrete cells between the solid and the liquid phase in the CA model. Then, the PF model tends to have larger computational cost than the CA model since it requires fine grid to demonstrate the smallest length scale for representing the smooth change in the phase field variable. Zaeem et al. [30] and Choudhury [31] have shown that the CA model requires less computational cost than the PF model even though it exhibits compatible accuracy in the simulation results for the dendrite growth supported by the analytical model (Lipton-Glicksman-Kurz, LGK).

Rappaz and Gandin [32,33] have developed a two-dimensional CA model based on the consideration of physical mechanisms on nucleation, growth kinetics of a dendrite tip, and crystallographic orientations. Zhu and Hong [34] proposed a modified CA model by applying the solute redistribution both in liquid and solid during solidification. Sanchez and Stefanescu [35] applied the solute conservation equation to determine the velocity of the solid-liquid (SL) interface while the previous studies [[32], [33], [34]] adopted the analytical solution based on the KGT (Kurz-Giovanola-Trivedi) model. Krane et al. [36] have proposed a numerical method to reduce an artificial anisotropy introduced by a Cartesian CA grid in the growth kinetics. However, the conventional approaches utilizing the CA model have assumed the undercooled static alloy melt, only, although the dendrite growth is substantially influenced by the fluid flow by various external forces. Therefore, there are increasing demands for the accurate prediction of the dendrite formation under the fluid flow condition [37,38].

In the current study, a coupled lattice Boltzmann method (LBM) - cellular automata (CA) model has been constructed to predict the microstructure evolution during the solidification process with taking into consideration of fluid flow of molten metal, which is able to demonstrate the dendrite morphologies and solute redistributions. To validate the constructed numerical model, the analytical solution of LGK model [39,40] is adopted for representing the tip velocity and radius of the dendrite when assuming the static melt without fluid flow condition. Oseen-Ivantsov solution [41] has been derived to examine the relation between growth and flow Peclet numbers under the fluid flow condition during the solidification. With supporting the validated LBM-CA model, as-casted microstructure in a directional solidification (DS) process has been examined with experimental results in terms of the dendrite formation and texture distribution.

## 2. Numerical model for solidification

This chapter mainly consists of two formulations concerning the computational fluid dynamics (CFD) for taking into consideration of fluid flow of the molten metal and solute transport by convection and diffusion, which is represented with the lattice Boltzmann method (LBM) [42]. Then, the CA determines the phase transformation to obtain the dendrites morphology according to the composition field previously calculated by the LBM based on the local equilibrium theory.

### 2.1. Lattice Boltzmann method (LBM)

The lattice Boltzmann method, a class of CFD methods, is adopted to solve the fluid flow and solute transport of the molten metal. Instead of solving Navier-Stokes equations for mass, momentum and energy conservation in traditional CFD, the LBM iteratively solves a set of nonlinear algebraic equations for fluid flow and solute transport by regarding those as collision and streaming processes of pseudo particles. The kinetic equations for fluid flow with the BGK (Bhatnagar-Gross-Krook) approximation [43], widely used for the collision operator, can be represented with Eq. (1) [42],

${{f}_{k}}\left( x+{{c}_{k}}\Delta t,t+\Delta t \right)={{f}_{k}}\left( x,t \right)\left[ 1-1/{{\tau }_{v}} \right]+f_{k}^{eq}\left( x,t \right)/{{\tau }_{v}}+F(x,t)$

where fk is the particle distribution function (PDF) to characterize the probability of a particular molecule at certain location (x) and time (t) along k -direction based on the lattice arrangement as depicted in Fig. 1, and fkeq, τv, F(x,t) indicate the equilibrium distribution function, relaxation time, and the external force, respectively. In the current research, Eq. (1) can be reduced to Eq. (2) by letting F(x,t) be zero since the external forces such as gravity and magnetic force are neglected.

${{f}_{k}}\left( x+{{c}_{k}}\Delta t,t+\Delta t \right)={{f}_{k}}\left( x,t \right)\left[ 1-1/{{\tau }_{v}} \right]+f_{k}^{eq}\left( x,t \right)/{{\tau }_{v}}$

### Fig. 1.

Fig. 1.   Lattice arrangement of D2Q9 model [46].

Similarly, the distribution function, sk, of the composition for a given alloy can be expressed with Eq. (3) in which Sk(x,t) denotes amount of rejected solute at interface.

${{s}_{k}}\left( x+{{c}_{k}}\Delta t,t+\Delta t \right)={{s}_{k}}\left( x,t \right)\left[ 1-1/{{\tau }_{D}} \right]+s_{k}^{eq}\left( x,t \right)/{{\tau }_{D}}+{{S}_{k}}\left( x,t \right)$

where ${{S}_{k}}\left( x,t \right)={{w}_{k}}\Delta {{f}_{s}}{{C}_{l}}(1-{{k}_{0}})$ [44].

Since the LBM equations calculate the macroscopic quantities such as density, temperature and fluid velocity with respect to k -direction in the single lattice, it is not only effective in terms of computing cost, but also possible to represent the distributions of each characteristics along k -direction, separately [45]. The D2Q9 lattice arrangement [46] is adopted in the current study. The discrete velocities, ck and weight coefficients, wk, in D2Q9 are expressed with Eqs. (4) and (5), where c=Δx/Δt indicates the lattice speed.

${{c}_{k}}=\left\{ \begin{matrix} c\left( 0,0 \right) \\ c\left( 1,0 \right),c\left( 0,1 \right),c\left( -1,0 \right),c\left( 0,-1 \right) \\ c\left( 1,1 \right),c\left( -1,1 \right),c\left( -1,-1 \right),c\left( 1,-1 \right) \\ \end{matrix}\begin{matrix} k=0, \\ k=1-4, \\ k=5-8 \\\end{matrix} \right.$
${{w}_{k}}=\left\{ \begin{matrix} 4/9 \\ 1/9 \\ 1/36 \\ \end{matrix} \right.\begin{matrix} k=0, \\ k=1-4, \\ k=5-8 \\ \end{matrix}$

The macroscopic variables such as fluid density, ρ, velocity, u, and composition, C, are represented with summation of all the distribution functions with respect to k-directions, respectively, as shown in Eq. (6).

$\rho =\underset{k=0}{\overset{8}{\mathop \sum }}\,{{f}_{k}},\rho u=\underset{k=0}{\overset{8}{\mathop \sum }}\,{{f}_{k}}{{c}_{k}},C=\underset{k=0}{\overset{8}{\mathop \sum }}\,{{s}_{k}}$

The equilibrium distribution functions in Eq. (2) and (3) are defined by Eq. (7) and 8) [42].

$f_{k}^{eq}\left( \text{x},\text{t} \right)={{w}_{k}}\rho \left[ 1+3({{c}_{k}}\cdot u)+4.5{{({{c}_{k}}\cdot u)}^{2}}-1.5(u\cdot u) \right]$$f_{k}^{eq}\left( \text{x},\text{t} \right)={{w}_{k}}\rho \left[ 1+3({{c}_{k}}\text{ }\!\!\cdot\!\!\text{ }u)+4.5{{({{c}_{k}}\text{ }\!\!\cdot\!\!\text{ }u)}^{2}}-1.5(u\text{ }\!\!\cdot\!\!\text{ }u) \right]$
$s_{k}^{eq}\left( \text{x},\text{t} \right)={{w}_{k}}C\left[ 1+3({{c}_{k}}\cdot u)+4.5{{({{c}_{k}}\cdot u)}^{2}}-1.5(u\cdot u) \right]$

According to the Chapman-Enskog theory [47], it is possible to obtain the transport parameters such as viscosity and diffusion coefficients based on the relation between macro and meso scales as expressed in Eqs. (9) and (10), where v is kinematic viscosity, D is diffusion coefficient.

$v=\frac{\Delta {{x}^{2}}\left( 2{{\tau }_{v}}-1 \right)}{6\Delta t}$
$D=\frac{\Delta {{x}^{2}}\left( 2{{\tau }_{D}}-1 \right)}{6\Delta t}$

### 2.2. Cellular automata (CA) model

The CA model partitions the simulation domain into discrete numbers of regular grids or cells to evaluate the current phase and the phase transformation, numerically, in order of each cell by taking into account the nucleation and grain growth during the solidification. The state of each cell is classified into the liquid, solid, and interface depending on the solid fraction. When the calculated solid fraction of the current cell is greater than unity, it is regarded as solid, while it is treated as the interface with the solid fraction between zero and unity. A liquid cell without solid fraction starts to transform to the solid state at the nucleation sites by capturing the neighboring cells in the state of the interface. Since there occurs the heterogeneous nucleation such that the mold surface or impurity tends to promote the nucleation formations than the homogeneous nucleation in the casting process, the probabilistic nucleation [32] was adopted to describe the heterogeneous nucleation in the mold surface and molten metal, simultaneously, as shown in Eq. (11).

$\frac{dn}{d\left( \Delta T \right)}=\frac{{{n}_{max}}}{\sqrt{2\pi }\Delta {{T}_{\sigma }}}\exp \left[ -\frac{1}{2}{{\left\{ \frac{\Delta T-\Delta {{T}_{m}}}{\Delta {{T}_{\sigma }}} \right\}}^{2}} \right]$

It represents the increase of the nuclei density, dn, with respect to the increase of the undercooling, d(ΔT) [33], where n is the nuclei density, nmax is the maximum density of nuclei, ΔTσ is the standard deviation of nucleation, and ΔTm is the mean nucleation undercooling for representing the Gaussian distribution of nucleation. Then, the nuclei density is calculated by integrating the distribution of the nucleation sites for the specified undercooling as expressed in Eq. (12).

$\delta n\left( \Delta {T}' \right)=\underset{\Delta {T}'}{\overset{\Delta {T}'+\delta (\Delta {T}')}{\mathop \int }}\,\frac{dn}{d\left( \Delta T \right)}d\left( \Delta T \right)$

If the solid and liquid compositions at their interface are in equilibrium, the solute partition between liquid and solid is determined by Eq. (13),

$C_{s}^{*}={{k}_{0}}C_{l}^{*}$

where Cs* and Cl* are the interface equilibrium composition in solid and liquid, respectively. The equilibrium composition at the interface Cl* is calculated with Eq. (14) [35],

$C_{l}^{\text{*}}={{C}_{0}}+\frac{{{T}^{\text{*}}}-T_{l}^{eq}+\Gamma Kf\left( \varphi ,\theta \right)}{{{m}_{l}}}$

where C0 is the initial solute concentration, T* is the interface equilibrium temperature, Tleq is the liquidus temperature at C0, ml is the liquidus slope, Γ is the Gibbs-Thomson coefficient, K is the curvature of the interface of the solid and liquid, and f(φ,θ) is the surface tension with a function of the angle of the surface to the interface φ and the angle of preferred growth orientation θ [32,33].

To update the state of the current cells with increase of solidification time, the change of solid fraction, Δfs, at each time step is calculated by the equilibrium composition [34,35] as expressed in Eq. (13), where Cl is the local liquid composition calculated by Eq. (15).

$\Delta {{f}_{s}}=\frac{C_{l}^{\text{*}}-{{C}_{l}}}{C_{l}^{\text{*}}\left( 1-{{k}_{0}} \right)}$

The anisotropy of the interface energy depending on the types of crystal structure is represented with Eq. (16) [48],

$f\left( \varphi ,\theta \right)=1-\delta \text{cos}\left[ 4\left( \varphi -\theta \right) \right]$

where δ=15ε is the amplitude of the directional anisotropy. To keep track of the change in the solid fraction, the interface curvature, K, was calculated, with Eq. (17) [21].

$K=\frac{2\frac{\partial {{f}_{s}}}{\partial x}\frac{\partial {{f}_{s}}}{\partial y}\frac{{{\partial }^{2}}{{f}_{s}}}{\partial x\partial y}-{{\left( \frac{\partial {{f}_{s}}}{\partial y} \right)}^{2}}\frac{{{\partial }^{2}}{{f}_{s}}}{\partial x\partial x}-{{\left( \frac{\partial {{f}_{s}}}{\partial x} \right)}^{2}}\frac{{{\partial }^{2}}{{f}_{s}}}{\partial y\partial y}}{{{\left[ {{\left( \frac{\partial {{f}_{s}}}{\partial x} \right)}^{2}}+{{\left( \frac{\partial {{f}_{s}}}{\partial y} \right)}^{2}} \right]}^{\frac{3}{2}}}}$

Likewise, the change of the solid fraction is updated at current time step based on the interface equilibrium condition, which makes it possible to predict the formation of the dendrite morphology, in detail. Overall iterative calculation is summarized in flowchart as depicted in Fig. 2.

### Fig. 2.

Fig. 2.   Overall flowchart for LBM-CA analysis.

## 3. LBM-CA validations

### 3.1. LGK model

To validate the simulation results of the proposed LBM-CA model, the LGK analytical model [39,40] has been adopted in terms of the tip velocity and radius of the dendrites when assuming the static melt without effect of fluid flow, which analytically expresses the tip radius and velocity with a function of the undercooling with assuming the shape of the dendrite tip as a paraboloid. Total undercooling is partitioned into the thermal, solutal, and curvature undercooling, respectively, as expressed in Eq. (18), where L is the latent heat, cp is the specific heat, ${{P}_{t}}\left( =\frac{RV}{2\alpha } \right)$. is the thermal Peclet number, ${{P}_{c}}\left( =\frac{RV}{2{{D}_{l}}} \right)$. is the solutal Peclet number. And Iv(P) is the Ivantsov function [47] in which E1(P) is the exponential integral function, and R is the tip radius obtained by the stability criterion [49,50] and V is the tip velocity of the dendrite.

$\Delta T=\frac{L}{{{\text{c}}_{p}}}{{I}_{v}}\left( {{P}_{t}} \right)+m{{C}_{0}}\left\{ 1-\frac{1}{1-\left( 1-{{k}_{0}} \right){{I}_{v}}\left( {{P}_{c}} \right)} \right\}+\frac{2\Gamma }{R}$
$\text{where}{{I}_{v}}\left( P \right)=P\text{exp}\left( P \right){{E}_{1}}\left( P \right),R=\frac{\Gamma }{{{\sigma }^{*}}}{{\left[ \frac{L}{{{c}_{p}}}{{P}_{T}}-\frac{{{P}_{C}}m{{C}_{0}}\left( 1-{{k}_{0}} \right)}{1-\left( 1-{{k}_{0}} \right){{I}_{v}}\left( {{P}_{C}} \right)} \right]}^{-1}}$

It is possible to solve Eq. (18) by the Newton-Raphson method, incrementally, since Pt and Pc in the Ivantsov function are nonlinearly correlated with each other in terms of R and V for a fixed undercooling and composition [51]. To make a comparison between the LGK analytic and numerical solutions from the LBM-CA model, the material properties of Al-Cu 4 wt% alloy has been adopted as summarized in Table 1. The size of LBM-CA entire domain is 200 × 200 μm, which is divided into 400 × 400 uniform square cells with a unit length of 0.5 μm. It is assumed that the initial nucleation is activated at the center of the domain as shown in Fig. 3(a) in which four boundary cells corresponding to the square solid cell are automatically designated to the interface cells. When the solid fraction of the interface cells becomes unity, they change to solid cells by capturing the interface cells due to the phase transformation as depicted in Fig. 3(b). It is possible to predict the dendrite tip velocity, numerically, since the phase transformation is calculated by dividing the capturing distance of the unit cells with respect to time. The tip radius also can be calculated with parabolic fitting method [35]. Fig. 4 demonstrates comparison results of the dendrite tip velocity and radius between the LGK analytic solution and LBM-CA prediction with increase of the undercooling, which exhibits good agreement with each other. It is concluded that the current CA model is successfully implemented to analyze the microstructure evolution during the solidification process.

Table 1   Material properties of Al-4 wt% Cu alloy [35].

Property and symbolValueUnits
Density, ρ2.475×103kg m-3
Solute diffusivity in liquid, Dl3×10-9m2 s-1
Gibbs-Thomson coefficient, Γ2.4×10-7mK
Partition coefficient, k00.17
Liquidus slope, mL-2.6K wt%-1
Liquidus temperature, Teq921.15K
Initial composition, C04wt%
Anisotropy coefficient, ε0.0267

### Fig. 3.

Fig. 3.   Schematic representation of each phases based on cell structure in CA analysis: (a) single solid cell surrounded by interface cells, (b) growth of solid and interface cells.

### Fig. 4.

Fig. 4.   Comparison between LGK analytic solution and LBM-CA prediction in terms of dendrite tip velocity and radius: (a) dendrite tip velocity, (b) dendrite tip radius.

### 3.2. Oseen-Ivantsov solution

Bouissou and Pelce [41], theoretically, examined the effect of a forced flow on the dendritic growth velocity under the solvability condition in which they derived a relation between Peclet number and amount of undercooling as expressed in Eq. (20) with assumptions for isothermal and parabolic dendrite tip with low Reynolds number,

$\Omega ={{P}_{c}}\text{exp}\left( {{P}_{c}}-{{P}_{f}} \right)\underset{1}{\overset{\infty }{\mathop \int }}\,\text{exp}\left\{ -{{P}_{c}}\eta +{{P}_{f}}\left( 2+\underset{1}{\overset{\eta }{\mathop \int }}\,\frac{g\left( \zeta \right)}{\sqrt{\zeta }}d\zeta -\eta \right) \right\}\frac{d\eta }{\sqrt{\eta }}$

where $\Omega \left( =\frac{C_{l}^{\text{*}}-{{C}_{\infty }}}{C_{l}^{\text{*}}\left( 1-{{k}_{0}} \right)} \right)$ is the dimensionless supersaturation, ${{P}_{c}}\left( =\frac{RV}{2D} \right)$ is the growth Peclet number, ${{P}_{f}}\left( =\frac{UR}{2D} \right)$ is the flow Peclet number, and a function, g(ζ) can be defined [41] by Eq.(21), in which $Re\left( =\frac{UR}{\nu } \right)$ is the Reynolds number.

$g\left( \zeta \right)=\frac{\sqrt{\zeta }erfc\left( \sqrt{\frac{Re\zeta }{2}} \right)+\sqrt{\frac{2}{\pi Re}}\left[ \text{exp}\left( -\frac{Re}{2} \right)-\text{exp}\left( -\frac{Re\zeta }{2} \right) \right]}{erfc\left( \sqrt{\frac{Re}{2}} \right)}$

Based on the Eqs. (20) and (21), the relation between flow Peclet number and growth Peclet number can be represented with specific dimensionless supersaturation and Reynolds number. To examine the relation between flow Peclet number and growth Peclet number of the proposed LBM-CA model under the fluid flow condition with Fe-0.82 wt% C alloy [52], it is assumed that a nucleus with crystallographic orientation of zero degree is set in center of the domain and forced flow of 0.1-150 mm/sec is applied to the horizontal direction as shown in Fig. 5(a). The dimensionless supersaturation value, Ω=0.226, equilibrium composition, Cl*/C0=1.175, are applied, and material properties of Fe-0.82 wt% C alloy are summarized in Table 2. Fig. 5(b) demonstrates asymmetric dendrite morphology with convectional melt since rejected solutes at the SL interface are washed away from the upstream to the downstream direction by the flowing melt [44] in which the tip radius is calculated based on the fitting with the 2nd order polynomial function as depicted in Fig. 5(b). By manipulating the analysis results in terms of Eq. (20) and (21), the relation between flow Peclet number and growth Peclet number of the proposed LBM-CA model has been plotted as shown in Fig. 6, which remarkably exhibits good agreement with the Oseen-Ivantsov solution. Therefore, it can conclude that the proposed LBM-CA model is successfully implemented in both of the static and convectional melt conditions to simulate the solidification procedure of the molten metals.

Table 2   Material properties of Fe-0.82 wt% C alloy [52].

Property and symbolValueUnits
Density, ρ7.020×103kg m-3
Solute diffusivity in liquid, Dl6.36×10-9m2 s-1
Gibbs-Thomson coefficient, Γ1.9×10-7mK
Partition coefficient, k00.34
Liquidus slope, mL-78.0K wt%-1
Liquidus temperature, Teq1809K
Initial composition, C00.82wt%
Anisotropy coefficient, ε0.04
Viscosity, μ5.5×10-3kg m-1 s-1

### Fig. 5.

Fig. 5.   Analysis results with proposed LBM-CA model under convectional melt: (a) geometrical boundary condition, (b) convectional melt.

### Fig. 6.

Fig. 6.   Comparison of Ossen-Ivantsov solution and proposed LBM-CA result.

## 4. Application to directional solidification (DS) process

In the previous chapter, the effects of diffusion and convection on dendrite growth in LBM-CA simulation were confirmed by comparison with the analytical solutions. In the current chapter, performance of LBM-CA model in prediction of microstructure is validated by experimental observations in directional solidification (DS) with nickel-based superalloy in terms of morphologies of equiaxed and columnar dendritic structure with respect to boundary conditions and effect of cooling rate during directional solidification. The DS process is widely used to manufacture a turbine blade with a single crystal by reducing thermal stress parallel to the axis of the thermal gradient to produce a sharp <001> crystallographic texture, which enhances the creep resistance under high operating temperature [[58], [59], [60], [61]] as shown in Fig. 7.

Fig. 7.   Directional solidification process: (a) schematic diagram, (b) cylindrical DS specimen.

Table 3 shows the chemical composition of Ni-based superalloy. The crucible was preheated by induction coil and maintained 1500 °C to melt the ingot as shown in Fig. 7(a) which is poured into the mold. After having a holding time of 30 s for stabilizing in the heat chamber, water-cooled chill plate with temperature of 20 °C is lowered with velocity of 200 μm/s to generate the cylindrical specimen consisting of directional grains as shown in Fig. 7(b). The cylindrical specimen has a diameter and length of 30 mm and 60 mm, respectively, in which a selected single grain grows to produce a blade through pigtail grain selector. Fig. 8 demonstrates the microstructures of the as-casted with optical microscopy (OM) along axial and radial directions of the cylindrical specimen, which were etched with etchant (3 g CuCl2, 30 ml HCl, 70 ml ethanol). The primary dendrite arm spacing (PDAS), λ1, and secondary dendrite arm spacing (SDAS), λ2, are measured as 428 μm and 103 μm, respectively, based on the OM observation.

Table 3   Chemical composition of the Ni-based superalloy.

ElementsCrAlTiTaWMoZrHfCoNi
wt.%7.945.220.673.079.070.600.021.239.21Bal.

Fig. 8.   As-casted microstructures in DS process along: (a) radial, (b) axial direction.

The commercial simulation package, ProCAST [62], for the DS has been utilized to reduce computational cost as well as validate the performance of the proposed LBM-CA in terms of the dendrite morphology prediction. Since ProCAST solve the fluid flow based on the full Navier-Stokes equation coupled with thermal analysis, it is effectively able to simulate the overall casting process including solidification defect and part distortion. However, ProCAST tends to demonstrate relatively low accuracy in predicting the dendrite growth, where the continuity condition of fluid velocity is not guaranteed due to solidification with multi-phase and micro-scale flow [44]. Then, temperature field in DS is calculated from ProCAST as shown in Fig. 9, which is applied to the LBM-CA analysis as a boundary condition as depicted in Fig. 10. With taking into consideration of the temperature and cooling rates with respect to solidification time during DS, the LBM-CA analysis has been carried out with the nucleation and material parameters presented in Table 3, Table 4. It is noted that computational cost of the current simulation substantially decreases to 3.3 % of stand-alone LBM-CA by integrating LBM-CA with ProCAST. The simulation domain is confined with 1000 × 1000 cells with a unit cell size of 4 μm, which represents 4 × 4 mm specimen. Fig. 11, Fig. 12 demonstrate the comparison of microstructure morphologies between experiment and numerical results including ProCAST and LBM-CA analyses in which colors represent the growth orientation, respectively. It is successively possible to simulate the dendrite morphologies with respect to growth directions with the LBM-CA code since the primary dendrites are promoted along the axial direction, dominantly, with exhibiting the secondary dendrites along the transverse direction as shown in Fig. 12(c) while it tends to grow along the radial direction as shown in Fig. 11(c). It is noteworthy that LBM-CA precisely describes the dendrite formation compared with ProCAST in terms of the PDAS and SDAS. The quantitative values of PDAS and SDAS from the LBM-CA are evaluated with 382 μm and 98 μm about 8.4 % and 4.8 % error, respectively, in experiment (Table 5).

Fig. 9.   Numerical results with ProCAST [62] on distribution of temperature and solid fraction for directional casing.

Fig. 10.   Temperature and cooling rate with respect to solidification time in ProCAST analysis.

Fig. 11.   Microstructure prediction along radial direction between experiment and numerical results: (a) experiment, (b) ProCAST, (c) LBM-CA.

Fig. 12.   Microstructure prediction along axial direction between experiment and numerical results: (a) experiment, (b) ProCAST, (c) LBM-CA.

Table 4   Nucleation parameters for numerical simulations.

Property and symbolValueUnits
Maximum density, nmax5.0×107m-3
Standard deviation, ΔTσ0.1K
Mean undercooling, ΔTm10K

Table 5   Material properties of Ni-based superalloy [53,54].

Property and symbolValueUnits
Density, ρ8.780×103kg m-3
Solute diffusivity in liquid, Dl3.6×10-9m2 s-1
Gibbs-Thomson coefficient, Γ3.65×10-7mK
Partition coefficient, k00.788
Liquidus slope, mL-3.95K wt%-1
Liquidus temperature, Teq1672K
Initial composition, C037.0wt%
Anisotropy coefficient, ε0.02
Viscosity, μ8.5×10-3kg m-1 s-1

In addition, to examine the effects of temperature gradient (G) and growth rate (V) of dendrite tip on the dendrites morphology in DS, they are varied from G = 50,000 K/m to G = 5,000 K/m with cooling rate of 1.5 K/s as depicted in Fig. 13(a-c) and 3.0 K/s as depicted in Fig. 13(d-f). The growth rates (V) of dendrite tips are calculated with $\text{V}=\text{\dot{T}}/\text{G}$, briefly, where $\dot{T}$ is the cooling rate. There are three types of dendrite morphologies in the axial direction, i.e., planar, cellular, and dendritic with respect to the process conditions. The planar type morphology tends to appear when temperature gradient is sufficiently high with low growth rate while the dendritic type morphology is dominant when temperature gradient is low with high growth rate, which has been supported by the other literatures based on the conventional solidification map [63,64]. Under these close examinations, it is concluded that the current LBM-CA code does not only predict the various types of dendrite morphologies, in detail, but also enable representation of the PDAS and SDAS with respect to the process conditions, precisely.

Fig. 13.   Various types of dendrite formation with respect to the processing conditions (cell size = 4 μm, domain size = 4 mm): (a-c) $\dot{T}$ = 1.5 K/s, (d-f) $\dot{T}$ =3.0 K/s.

## 5. Conclusions

This paper mainly concerns novel approach to simulate the solidification of the molten metal, numerically, with supported by LBM-CA integration in terms of dendrite morphology and microstructure evolution. The LBM is adopted to solve the mass transport by diffusion and convection during solidification, and the CA model is utilized to determine the phase transformation with respect to the solid fraction based on the local equilibrium theory. The proposed LBM-CA model has been validated by analytic solutions such as LGK model in static melt, and Oseen-Ivantsov solution under the fluid flow conditions, which shows good agreement in terms of dendrite growth and morphology. In order to evaluate the performance of the proposed LBM-CA model, DS experiment with Ni-based superalloy has been conducted to make a comparison with numerical prediction by LBM-CA concerning the PDAS and SDAS, etc. Under these circumstances, final conclusion can be drawn:

(1) Implementation of LBM does not only enable precise representation of dendrite morphology by regarding fluid flow and solute transport as collision and streaming processes of pseudo particles in meso-scale, but also possible to represent the distributions of each properties along k -direction, separately, in a single lattice.

(2)LBM-CA model makes it possible to describe different types of dendrite formations with respect to various solidification conditions, which were confirmed by precise predictions in PDAS and SDAS with those in DS experiment about 8.4 % and 4.8 % error.

(3)With supporting numerical prediction in dendrite formation, the effects of process conditions such as temperature gradient and growth rate are able to be evaluated and calibrated, which enables tremendous cost reduction induced by repetitive experiments and trials.

## Declaration of Competing Interest

The authors have no conflict of interests to declare.

## Acknowledgments

This research was financially supported by the Ministry of Trade, Industry, and Energy (MOTIE), Korea, under the “Digital manufacturing platform (DigiMaP)” (reference number N0002598) supervised by the Korea Institute for Advancement of Technology (KIAT). This work was also supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (2019R1A2C4070160).

## Reference By original order By published year By cited within times By Impact factor

S.G. Shabestari, F. Shahri, J. Mater. Sci. 39 (2004) 2023-2032.

S.C. Huang, M.E. Glicksman, Acta Metall. 29 (1981) 717-734.

S. Terzi, L. Salvo, M. Suery, A.K. Dahle, E. Boller, Acta Metall. 58 (2010) 20-30.

M. El-Bealy, B.G. Thomas, Metall. Mater. Trans. B 27B (1996) 689-693.

J.M.V. Quaresma, C.A. Santos, A. Garcia, Metall. Mater. Trans. A 31A (2000) 3167-3177.

C.C. Hays, C.P. Kim, W.L. Johnson, Mater. Sci. A-Struct. 304 (2001) 650-655.

Y.S. Oh, C.P. Kim, S. Lee, N.J. Kim, Metall. Mater. Trans. A 43 (2012) 1911-1920.

A. Farzadi, M. Do-Quang, S. Serajzadeh, A.H. Kokabi, G. Amberg, Model Simul. Mater. Sci. Eng. 16 (2008), 0650052.

Y.S. sato, M. Urata, H. Kokawa, Metall. Mater. Trans. A 33 (2002) 625-635.

K.V. Jata, K.K. Sankaran, J.J. Ruschau, Metall. Mater. Trans. A 31 (2000) 2181-2192.

P. Cavaliere, R. Nobile, F.W. Panella, A. Squillace, Int. J. Mach. Tool. Manu. 46 (2006) 588-594.

S. Sahoo, K. Chou, Acta Mater. 77 (2014) 85-95.

A multi-scale model that combines the finite element method and stochastic analysis is developed to simulate the evolution of the microstructure of an Nb-bearing nickel-based superalloy during laser additive manufacturing solidification. Through the use of this model, the nucleation and growth of dendrites, the segregation of niobium (Nb) and the formation of Laves phase particles during the solidification are investigated to provide the relationship between the solidification conditions and the resultant microstructure, especially in the morphology of Laves phase particles. The study shows that small equiaxed dendrite arm spacing under a high cooling rate and low temperature gradient to growth rate (G/R) ratio is beneficial for forming discrete Laves phase particles. In contrast, large columnar dendrite arm spacing under a low cooling rate and high G/R ratio tends to produce continuously distributed coarse Laves phase particles, which are known to be detrimental to mechanical properties. In addition, the improvement of hot cracking resistance by controlling the morphology of Laves phase particles is discussed by analyzing the cracking pattern and microstructure in the laser deposited material. This work provides valuable understanding of solidification microstructure development in Nb-bearing nickel-based superalloys, like IN 718, during laser additive manufacturing and constitutes a fundamental basis for controlling the microstructure to minimize the formation of deleterious Laves phase particles. (C) 2014 Acta Materialia Inc. Published by Elsevier Ltd.

A. Zinoviev, O. Zinovieva, V. Ploshikhin, V. Romanova, R. Balokhonov, Mater. Des. 106 (2016) 321-329.

O. Zinovieva, A. Zinoviev, V. Ploshikhin, Comput. Mater. Sci. 141 (2018) 207-220.

S. Majaniemi, N. Provatas, Phys. Rev. E 7 (2009), 011607.

M. Chiumenti, M. Cervera, E. Salsi, A. Zonato, J. Heat Transfer 140 (2018), 082301.

I.L. Ferreira, C.A. Santos, V.R. Voller, A. Garcia, Metall. Mater. Trans. B 35 (2004) 285-297.

The present work focuses on the influence of alloy solute content, melt superheat, and metal/mold heat transfer on inverse segregation during upward solidification of Al-Cu alloys. The experimental segregation profiles of Al 4.5 wt pct Cu, 6.2 wt pct Cu, and 8.1 wt pct Cu alloys are compared with theoretical predictions furnished by analytical and numerical models, with transient h i profiles being determined in each experiment. The analytical model is based on an analytical heat-transfer model coupled with the classical local solute redistribution equation proposed by Flemings and Nereo. The numerical model is that proposed by Voller, with some changes introduced to take into account different thermophysical properties for the liquid and solid phases, time variable metal/mold interface heat-transfer coefficient, and a variable space grid to assure the accuracy of results without raising the number of nodes. It was observed that the numerical predictions generally conform with the experimental segregation measurements and that the predicted analytical segregation, despite its simplicity, also compares favorably with the experimental scatter except for high melt superheat.]]>

W.U. Mirihanage, D.J. Browne, Comput. Mater. Sci. 46 (2009) 777-784.

L. Nastac, S. Sundarraj, K. Yu, Y. Pang, JOM 50 (1998) 30-35.

C. Charbon, M. Rappaz, Modell. Simul. Mater. Sci. Eng. 1 (1993) 455-466.

L. Nastac, D.M. Stefanescu, Modell. Simul. Mater. Sci. Eng. 5 (1997) 391-420.

W.J. Boettinger, J.A. Warren, C. Beckermann, A. Karma, Annu. Rev. Mater. Res. 32 (2002) 163-194.

S.L. Wang, R.F. Sekerka, A.A. Wheeler, B.T. Murray, S.R. Coriell, R.J. Braun, G.B. McFadden, Phys. D 69 (1993) 189-200.

Y. Zhao, B. Zhang, H. Hou, W. Chen, M. Wang, J. Mater. Sci. Technol. 35 (2019) 1044-1052.

H. Pan, Z. Han, B. Liu, J. Mater. Sci. Technol. 32 (2016) 68-75.

Y. Zhao, R. Qin, D. Chen, X. Wan, Y. Li, M. Ma, Steel Res. Int. 86 (2015) 1490-1497.

ISSN: 1005-0302
CN: 21-1315/TG
Editorial Office: Journal of Materials Science & Technology , 72 Wenhua Rd.,
Shenyang 110016, China
Tel: +86-24-83978208
E-mail:JMST@imr.ac.cn