Multi-scale coupling simulation in directional solidification of superalloy based on cellular automaton-finite difference method
Multi-scale coupling simulation in directional solidification of superalloy based on cellular automaton-finite difference method
Document code: A
0 Jian-xin Zhou, Ya-jun Yin, Dong-qiao Zhang, Xiao-yuan Ji and Xu Shen State Key Laboratory of Materials Processing and Die & Mould Technology, Huazhong University of Science and Technology , Wuhan 430074
Casting microstructure evolution is difficult to describe quantitatively by only a separate simulation of dendrite scale or grain scale, and the numerical simulation of these two scales is difficult to render compatible. A three-dimensional cellular automaton model couplling both dendritic scale and grain scale is developed to simulate the microstructure evolution of the nickel-based single crystal superalloy DD406. Besides, a macro-mesoscopic/ microscopic coupling solution algorithm is proposed to improve computational efficiency. The simulation results of dendrite growth and grain growth of the alloy are obtained and compared with the results given in previous reports. The results show that the primary dendritic arm spacing and secondary dendritic arm spacing of the dendritic growth are consistent with the theoretical and experimental results. The mesoscopic grain simulation can be used to obtain results similar to those of microscopic dendrites simulation. It is indicated that the developed model is feasible and effective.
multi-scale coupling; dendritic growth; grain growth; directional solidification; cellular automata
irectional solidification technology is one of
in solidification microstructure simulation at the
mesothe key technologies for manufacturing high
grain and micro-dendritic scale.
performance and complex parts such as aero engines and
gas turbines blades [
]. In the manufacturing process of
In multi-scale solidification microstructure simulation,
the macroscopic temperature and the microstructure
blades, microstructure defects such as misoriented grains
evolution of solidification are generally calculated by
and freckles are easily formed, which seriously reduce
using the coupled method of the CA method and the
the overall performance of the blades. The solidification
finite element (FE) method or the finite difference (FD)
simulation method can help to analyze the mechanism
method. For instance, based on CA-FE method, the 3D
of the microstructure evolution and optimize the
decentered square method, proposed by Rappaz and
solidification process parameters during casting.
Gandin et al.[
], was used to simulate the directional
At present, the most widely used simulation methods
solidification process of the blade. Then, based on the
for solidification microstructure include the
Volumefront tracking method, the CA-FE method was modified
averaged, the Phase Field, the Monte Carlo and the
to solve the heat flow problem in the solidification mushy
Cellular Automaton (CA) methods [
]. Among them, the
zone of casting. Similarly, based on the commercial
CA method has the advantages of clearly describing the
casting simulation software ProCAST, the numerical
physical meaning and high computational efficiency of
simulation of the directional solidification process of the
solidification phase transition, and has been widely used
and Meng et al.[
]. In addition, the CA-FD method was
also developed by Xu et al. [
] to simulate the directional
solidification process of the superalloy turbine blades.
It is important to note that the CA method mentioned
above is only applied to the mesoscopic grain scale.
The CA model of mesoscopic grain scale is based on
where σ is the Stefan-Boltzmann radiation constant, T1, T2 are
the temperature of any two finite gray bodies, and ε1, ε2 are the
emissivity, and A1, A2 are the occupied boundary surface area,
their values equal to the area of the grid unit, and X1,2 is the
radiation angle factor, the value is taken as constant 1.
1.2 Meso-microscopic model of dendritic and grain evolution
In the solidification microstructure simulation, the continuous
nucleation model proposed by Rappaz and Gandin et al. [
used to describe the random nucleation process of dendrites/
grains, as shown in Eq. (3):
where ΔT is the undercooling, n(ΔT) is the nucleation number
density at ΔT, nmax is the initial nucleation base number, ΔTσ is
the undercooling variance in nucleation, ΔTN is the undercooling
mean in nucleation, respectively.
The growth rate of the dendritic interface and the non-dendritic
interface are calculated by using a unified and simplified analytical
model proposed by Nastac and Stefanescu [
]. At the same time,
considering the influence of solute component supercooling and
curvature on the growth rate of dendrites in the simulation, the
phase equilibrium equation and the solute diffusion equation at the
interface are integrated to calculate. The main governing equations
of the model are given by
the theory of solidification interface stability. It is the solution of
the dendritic growth velocity, the tip radius of the dendrites and
the supercooling degree of the interface, ignoring the dendritic
morphology of the grain itself; the numerical calculation is
very fast. For the CA model of the microscopic dendritic scale,
solving the solidification interface growth rate [
] is based on
the solute or energy conservation relationship at the front of the
dendritic growth interface, taking full account of the influence of
the solute composition and interface curvature on the dendritic
], and the numerical solution is very accurate.
However, the former is based on the analytical model, and the
latter is based on conservation model. It is difficult to couple the
two in solidification simulation.
Accordingly, a typical nickel-based single crystal superalloy
DD406 was studied. Based on the advantages of the analysis
model and the conservation model, a three-dimensional CA model
coupling dendrite scale and grain size is developed. First of all,
considering the solute undercooling, curvature undercooling
and other factors, the analytical model of the dendrite tip growth
velocity in front of the interface is used to calculate the solid
fraction. Secondly, the modified decentered square method was
used to calculate the growth rate of arbitrary dendritic/grain
interface, and the solid fraction in the interface was coupled with
the modified decentered square method to track the interface.
Finally, the 3D CA model proposed in this study is compared
with the theoretical and experimental results given in previous
studies. Combined with the actual macro scale flow, the effect of
different directional solidification parameters on the microstructure
evolution of the superalloy DD406 was studied.
1 Model and governing equations
1.1 Macroscopic model of heat transfer
The governing equation for macroscopic temperature field is
where ρ is density, cp is the specific heat capacity, T is the
temperature, t is the solidification time, ki is the thermal
conductivity, L is the solidification latent heat, FS is the solid
fraction, and QR is the radiant heat source during the directional
In addition, for the calculation of the radiation heat source,
the method of High Rate Solidification (HRS) in the directional
solidification is taken as an example to simplify the radiation
heat transfer between the heating zone, the cooling zone and
the simulated calculation area in the furnace. The heat transfer
between any two finite gray bodies is calculated by
StefanBoltzman's law. The governing equation is given by
where V is the interface growth rate, Γ is the Gibbs-Thomson
coefficient, m is the liquidus slope, k0 is the partition coefficient,
C0 is the initial composition at the temperature T0eq, CL* is the
interface equilibrium composition at the temperature TLeq, DL is
the solute diffusion coefficient in liquid, fS is the solid fraction,
K is the curvature of the interface calculated by an approximate
2 Numerical algorithm and solution
2.1 Macro–mesoscopic/microscopic coupling solution
The CA-FD method is used to simulate the directional
solidification process of superalloy. Two sets of uniform
numerical grid ΔX, Δx (Δx=ΔX/10) are adopted, with the same
time step Δt. For the microscopic scale simulation, the cell
states are defined as three types: liquid (fs = 0), solid (fs = 1) or
interface (0<fs<1). And the cell neighbours are characterized
by the first-order neighbours (6 neighbourhoods), second-order
neighbours (12 neighbourhoods) and third-order neighbours (8
The solid phase fraction of mesoscopic/microscopic cell [
calculated as Eq. (7):
In order to ensure the convergence of the numerical
calculation, the solid phase fraction of the cell is limited to 1 in a
time step Δt, that is to say, the solidification interface is no more
than one unit mesh.
The solid phase fraction of macroscopic cell is calculated as
In the first step, the macroscopic temperature field is
calculated by the finite difference method. Secondly, the
macroscopic temperature field is mapped to the temperature field
of the mesoscopic/microscopic cell by the linear interpolation
method. Then, the dendritic growth rate and solid fraction are
calculated by the supercooling functions. Finally, the solid
fraction of the mesoscopic/microscopic cells are mapped back
to the solid fraction of the macroscopic cell, accompanied
by the release of latent heat, which affects the distribution of
the macroscopic temperature field. The macro-mesoscopic/
microscopic coupling solution is shown in Fig. 1.
2.2 Cell transformation or capture rules in cellular automaton
The modified decentered square method is adopted. The
general algorithm is used to track the six vertices of the
octahedron. When one of six vertices touches the liquid cell
in the surrounding neighbor, it is considered that the neighbor
liquid cell is captured and transformed into the interface cell.
At the same time, the growth center of the new interface cell
is calculated from the decentered square growth of the initial
interface cell growth center. In time step Δt, the six vertices of
the octahedron may touch one of the surrounding 26 neighbors’
liquid cells, and even touch two neighbor liquid cells, resulting
in the complexity of decentered square algorithm. In order to
avoid the complex calculation of the newly captured interface
cell, it is modified as follows:
(1) When the vertices of the octahedron reach only the
firstorder nearest neighbor liquid cell within time step Δt, the
vertex at this time is used as the growth center of the newly
captured interface cell, as shown in Fig. 2.
(2) In a time step Δt, while the first-order nearest neighbor
liquid cell and the second-order neighbor liquid cell are
touched at the same time, the vertex at this time is used as
the growth center of the two newly captured interface cells.
However, the growth rate is still calculated according to the
undercooling of their respective cells.
Fig. 2: Schematic diagram of modified decentered
3 Results and discussion
The material properties of DD406
(Ni-3.92Cr-9.34Co1.85Mo-8.88W-5.92Al-7.49Ta-2.0RE, mass fraction wt.%) is
calculated by using commercial JMatPro software.The material
properties and model parameters used in the simulations are
given in Table 1.
3.1 Dendritic morphology and dendritic spacing
To verify the correctness of the simulation, the dendrite
evolution of directional solidification under the boundary
conditions of 6 crystal seeds (n) with a temperature gradient (G)
of 12100 K•m-1 and a cooling rate Cr of 8.4 K•s-1 was simulated
and analyzed, as shown in Fig. 3.
Figure 3 shows that the initial spacing of seed is equivalent
to the primary dendritic arm spacing at the initial time. But with
the increase of the undercooling, the secondary dendrites are
formed on the trunk of dendrite, even the tertiary dendrite arms
are also formed on the secondary dendritic arms. The tertiary
dendrite arms lead to a decrease in the primary dendritic arm
spacing, as shown in Fig. 3 (a) and (b). Besides, the competition
growth between dendrites can lead to the increase of local
primary dendrite spacing, as shown in Fig. 3 (b) and (c). This
is caused by the fluctuation of the influencing factors of the
solidification process, such as solute composition, temperature,
etc. Finally, at t = 2.5 s, the height of the primary dendritic arms
is approximately the same. It is considered that the dendritic
growth tends to be stable and the size of the primary dendritic
spacing of the subsequent dendritic growth is approximately
constant. It can be concluded that the growth behavior of
dendrites is not only dependent on the parameters of initial
solidification boundary conditions, but also related to the results
of competitive growth history, which is consistent with the
relevant theoretical and experimental conclusions [
In order to quantitatively study the dendritic growth, the average
spacing of primary dendrites and secondary dendrites were
simulated and calculated under the boundary conditions such as
initial spacing, temperature gradient and cooling rate. According to
the classical primary dendritic arm spacing models of Hunt,
KurzFisher and Trivedi [
] and the secondary dendritic arm spacing of
Feurer-Wunderlin, the dendritic growth is analyzed and compared.
The comparison between the simulation results and the theoretical
calculation results are shown in Fig. 4 and Fig. 5.
From Fig. 4 and 5, it is known that the simulation results of
primary dendritic arm spacing and secondary dendritic arm
spacing are consistent with those of classical theory, which show
that the proposed and developed model and algorithm are feasible.
(a) t=1.0 s (b) t=1.5 s (c) t=2.0 s (d) t=2.5 s
Fig. 3: Dendritic morphology evolution during directional
solidification at n = 6, G = 12,100 K•m-1 and
Cr = 8.4 K•s-1
3.2 Evaluation of grain growth simulation
The same calculation model of dendritic tip growth rate is
used in the microscopic dendrite simulation and mesoscopic
grain simulation. Based on the same parameters and boundary
conditions, the evolution of solidification microstructure between
grain and dendrite simulation are compared and analyzed. The
mesoscopic grain growth and the microscopic dendrite growth at
temperature gradient of 12,100 K•m-1 and cooling rate of 14 K•s-1
are shown in Fig. 6 and 7, respectively.
Figure 6 and 7 show that, at the beginning of grain or
dendritic growth, the growth rate of grain or dendrites is
fast, and it is decreases gradually. This is because the grain
or dendrite growth behavior is impeded by the boundary of
computation domain. Finally, the two grains or dendrites
with orientation a1 and a2 are obtained. It appears that the
mesoscopic grain simulation can be used to obtain results
similar to those of microscopic dendrites simulation. But the
biggest difference between the two is that the growth time
of mesoscopic grain does not match the time of microscopic
dendrite growth, that is to say, the growth rate is faster in
the mesoscopic grain simulation. The reason is that the
mesoscopic grain growth rate is obtained by the KGT model
fitting under the supercooling function. But the model is only
suitable for describing the stable growth and evolution of the
equiaxed crystal of the ideal unrestricted region. Besides,
the solute concentration has an important influence on the
evolution of solidification microstructure. In the present work,
the computation domain is limited. The tip growth rate of each
grain is unstable, which leads to the difference of the results.
Therefore, in order to accurately simulate the growth behavior
of the mesoscopic grain, the growth rate of the dendritic tip of
the mesoscopic grain must be calculated more accurately.
3.3 Simulation of directional solidification
process of single crystal superalloy
In the directional solidification process of superalloy, the drawing
speed is the key process parameter. The directional solidification
process of single crystal superalloy under different drawing
speeds of 6, 9 and 12 mm•min-1 are simulated by the mesoscopic
grain simulation method. In order to accurately simulate the
mesoscopic grain growth, the mesoscopic grain growth rate is
obtained indirectly by the fitting function of the average growth
rate in the microscopic dendritic growth simulation.
The distribution of the temperature field, the solid-liquid
interface and the mesoscopic grain structure during the
directional solidification are shown in Fig. 8.
Figure 8 shows that the temperature gradient of the spiral
selector is approximately along the Z-axis direction at the initial
time, when the solid-liquid interface enters the spiral selector,
and its temperature gradient is similar to the helix centerline
of the spiral selector. But at time t=590 s, it can be seen that
the solid-liquid interface is no longer along the helix centerline
of the spiral selector, but forms the mushy zone near the
(a) t=250 s
(b) t=430 s
(c) t=510 s (d) t=590 s
Fig.8: Evolution of temperature field, solid-liquid interface and mesoscopic grain structure
during directional solidification at 6 mm•min-1
central wall, which may affect the changes in the solidification
microstructure. Finally, due to the efficient selection of the
spiral selector, only one grain in the solidification microstructure
grows along the spiral selector and forms a single crystal.
In order to further quantitatively characterize the growth and the
selection efficiency of the directional solidification microstructure
along the Z-axis at different drawing speeds, the grains number
along the Z-axis direction is analyzed. Because of the randomness
of nucleation, the grain orientation is characterized by a random
integer. The maximum integer is taken as the maximum grain
number. The growth distance along the Z axis of different grains
at different drawing speeds is shown in Fig. 9.
(a) Geometry model
(b) 6 mm•min-1
(d) 9 mm•min-1 (d) 12 mm•min-1
Fig. 9: Growth distance along Z axis of different grains at different drawing speeds
From Fig. 9, most of the grains can be eliminated in the starter
block. At drawing speeds of 6, 9, and 12 mm•min-1, there are 5,
8 and 7 grains with different crystallographic orientations that
can be selected out respectively. That is to say, reducing the
drawing speed can cut down the number of grains selected from
the starter block, thus improving the crystal selection efficiency.
In order to quickly and accurately predict the evolution of
solidification microstructure, a three-dimensional cellular
automaton model coupled dendritic scale and grain scale is
developed, and a macro-mesoscopic/microscopic coupling
solution algorithm is proposed. The nickel-based single crystal
superalloy DD406 was researched and the conclusions can be
drawn as follows:
(1) The results of the dendritic scale simulation are consistent
with the theoretical and experimental results given in previous
studies. It can be considered that the developed
threedimensional cellular automaton model is feasible.
(2) Even though the growth time of mesoscopic grain does not
match the time of microscopic dendrite growth, similar results
can be obtained from the two kinds of simulations.
(3) The parameters of drawing speed and geometry size of
starter block and grain selector all can affect the solidification
microstructure. Under the same nucleation and random
orientation, the final growth distance along the Z axis of
different grains at different drawing speeds is similar.
 Li Y , Liu L , Huang T , et al. Simulation of stray grain formation in Ni-base single crystal turbine blades fabricated by HRS and LMC techniques . China Foundry , 2017 , 14 ( 2 ): 75 - 79 .
 Meng X , Li J , Zhu S , et al. Method of stray grain inhibition in the platforms with different dimensions during directional solidification of a ni-base superalloy . Metallurgical and Materials Transactions A: Physical Metallurgy and Materials Science , 2014 , 45 ( 3 ): 1230 - 1237 .
 Cao L , Liao D , Lu Y , et al. Heat transfer model of directional solidification by lmc process for superalloy casting based on finite element method . Metallurgical and Materials Transactions A: Physical Metallurgy and Materials Science , 2016 , 47 ( 9 ): 4640 - 4647 .
 Zhang H , Xu Q. Multi-scale simulation of directional dendrites growth in superalloys . Journal of Materials Processing Technology , 2016 , 238 : 132 - 141 .
 Szeliga D , Kubiak K , Motyka M , et al. Directional solidification of Ni-based superalloy castings: Thermal analysis . Vacuum , 2016 , 131 : 327 - 342 .
 Zhi X , Liu J , Xing J , et al. Effect of cerium modification on microstructure and properties of hypereutectic high chromium cast iron . Materials Science and Engineering: A , 2014 , 603 : 98 - 103 .
 Stefanescu D M. Science and Engineering of Casting Solidification. Switzerland: Springer, 2015 .
 Gandin C , Desbiolles J , Rappaz M , et al. A three-dimensional cellular automaton-finite element model for the prediction of solidification grain structures . Metallurgical and Materials Transactions A , 1999 , 30 ( 12 ): 3153 - 3165 .
 Gandin C , Rappaz M. A 3D Cellular Automaton algorithm for the prediction of dendritic grain growth . Acta Materialia , 1997 , 45 ( 5 ): 2187 - 2195 .
 Carozzani T , Digonnet H , Bellet M , et al. 3D CAFE simulation o f a m a c r o s e g r e g a t i o n b e n c h m a r k e x p e r i m e n t . I O P Conference Series: Materials Science and Engineering, 2012 , 33 ( 1 ): 12087 - 12096 .
 Guillemot G , Gandin C , Combeau H , et al. A new cellular a u to ma to n -fi n i te e l e me n t co u p l i n g sch e me fo r a l l o y solidification . Modelling and Simulation in Materials Science and Engineering , 2004 , 12 ( 3 ): 545 - 556 .
 Carter P , Cox D C , Gandin C A , et al. Process modelling of grain selection during the solidification of single crystal superalloy castings . Materials Science and Engineering A , 2000 , 280 ( 2 ): 233 - 246 .
 Meng X , Lu Q , Li J , et al. Modes of grain selection in spiral selector during directional solidification of nickel-base superalloys . Journal of Materials Science & Technology , 2012 , 28 ( 3 ): 214 - 220 .
 Xu Q , Zhang H , Qi X , et al. Multiscale Modeling and Simulation of Directional Solidification Process of Turbine Blade Casting with MCA Method. Metallurgical and Materials Transactions B , 2014 , 45 ( 2 ): 555 - 561 .
 Nastac L. Numerical modeling of solidification morphologies and segregation patterns in cast dendritic alloys . Acta Materialia , 1999 , 47 ( 17 ): 4253 - 4262 .
 Zhu M , Stefanescu D. Virtual front tracking model for the quantitative modeling of dendritic growth in solidification of alloys . Acta Materialia , 2007 , 55 ( 5 ): 1741 - 1755 .
 Pan S , Zhu M. A three-dimensional sharp interface model for the quantitative simulation of solutal dendritic growth . Acta Materialia , 2010 , 58 ( 1 ): 340 - 352 .
 Zhang X , Zhao J , Jiang H , et al. A three-dimensional cellular automaton model for dendritic growth in multi-component alloys . Acta Materialia , 2012 , 60 ( 5 ): 2249 - 2257 .
 Yuan L , Lee P D. A new mechanism for freckle initiation based on microstructural level simulation . Acta Materialia , 2012 , 60 ( 12 ): 4917 - 4926 .
 R a p p a z M , G a n d i n C A . P r o b a b i l i s t i c m o d e l l i n g o f microstructure formation in solidification processes . Acta Metallurgica et Materialia , 1993 , 41 ( 2 ): 345 - 360 .
 Nastac L , Stefanescu D M. An analytical model for solute redistribution during solidification of planar, columnar, or equiaxed morphology . Metallurgical Transactions A , 1993 , 24 ( 9 ): 2107 - 2118 .
 Wang W , Lee P D , Mclean M. A model of solidification microstructures in nickel-based superalloys: predicting primary dendrite spacing selection . Acta Materialia , 2003 , 51 ( 10 ): 2971 - 2987 .
 Kurz W , Fisher D J. Dendrite growth at the limit of stability: tip radius and spacing . Acta Metallurgica , 1981 , 29 ( 1 ): 11 - 20 .
 Trivedi R. Interdendritic Spacing: Part II. A Comparison of Theory and Experiment . Metallurgical Transactions A , 1984 , 15 ( 6 ): 977 - 982 .