Equiaxed and columnar dendrite growth simulation in Al-7Si- Mg ternary alloys using cellular automaton method

In this paper, a modified cellular automaton (MCA) model allowing for the prediction of dendrite growth of Al-Si-Mg ternary alloys in two and three dimensions is presented. The growth kinetic of S/L interface is calculated based on the solute equilibrium approach. In order to describe the dendrite growth with arbitrarily crystallographic orientations, this model introduces a modified decentered octahedron algorithm for neighborhood tracking to eliminate the effect of mesh dependency on dendrite growth. The thermody namic and kinetic data needed for dendrite growth is obtained through coupling with Pandat software package in combination with thermodynamic/kinetic/equilibrium phase diagram calculation databases. The effect of interactions between various alloying elements on solute diffusion coefficient is considered in the model. This model has first been used to simulate Al-7Si (weight percent) binary dendrite growth followed by a validation using theoretical predictions. For ternary alloy, Al-7Si-0.5Mg dendrite simulation has been carried out and the effects of solute interactions on diffusion matrix as well as the differences of Si and Mg in solute distribution have been analyzed. For actual application, this model has been applied to simulate the equiaxed dendrite growth with various crystallographic orientations of Al-7Si-0.36Mg ternary alloy, and the predicted secondary dendrite arm spacing (SDAS) shows a reasonable agreement with the experimental ones. Furthermore, the columnar dendrite growth in directional solidification has also been simulated and the predicted primary dendrite arm spacing (PDAS) is in good agreement with experiments. The simulated results effectively demonstrate the abilities of the model in prediction of dendritic microstructure of Al-Si-Mg ternary alloy.


Introduction
Due to an excellent combination of castability, high corrosion resistance and comprehensive mechanical properties of Al-7Si-Mg ternary casting alloy, typically A356, A357, ZL114, it has been widely used in the automotive and aerospace industries as structural component [1,2] . Since the characteristic of dendrite has a dominant effect on the mechanical properties of castings, its understanding of the dendritic microstructural formation is of great importance to control the desirable microstructure, and thereby to modify the performance of castings. Numerical simulation has provided a powerful tool to study the dendrite structure during solidification, and cellular automaton (CA) model, as an effective method with high calculation efficiency and clear physical meaning, has been widely applied in the investigation of dendrite growth process [3~13] . However, most of the existing CA models are confined to describe the dendrite of binary alloys with well-defined crystallographic axes, while in practice, most commercial alloys are multicomponent with three or more elements, and the solidification process involves complex polycrystalline growth with various crystallographic orientations. Besides, the work on threedimensional dendrite simulation using CA model is also far from satisfaction and further effort is still needed. In addition, carrying out the simulation combining with the practical industrial production and enhancing the practicability of CA model is one of the further developments [14] .
In this paper, through coupling with Pandat software package in combination with thermodynamic/kinetic/equilibrium phase diagram calculation databases, a modified cellular automaton (MCA) model allowing for the prediction of dendrite growth of ternary alloys in two and three dimensions is presented. This model has been first used to simulate Al-7Si binary dendrite for model validation, and then applied to Al-7Si-0.5Mg ternary alloy simulation for the investigation of the different solute distribution behavior between Si and Mg as well as the effects of solute interactions on dendrite growth. Finally, to show the model ability in practical application, it is applied to simulate the equiaxed dendrite growth and columnar dendrite growth of Al-7Si-0.36Mg alloy.

Solute distribution calculation
To describe mass transfer of Al-7Si-Mg ternary system, the solute fields of Si and Mg need to be calculated separately. In this model, the effect of solute interaction between Si and Mg elements on solute diffusion coefficient is considered, and the solute diffusion is calculated based on the following equation where w is the composition with a superscript φ, denoting the liquid or solid phase, and the subscript i refers to the alloying element Si or Mg. k i is the composition and temperature dependent equilibrium partition coefficient which is obtained through coupling with equilibrium phase diagram database. f S is the solid fraction. J i is the diffusion flux of element i given by Fick-Onsager law as: where D ij φ is the solute diffusivity matrix for phase φ. Due to the fact that Si and Mg atoms are substitutional, D ij φ can be expressed as [15] Al,Si,Mg Al where x ν φ is the mole fraction of component ν, and the parameter M ν φ is the atomic diffusion mobility of element ν in the given phase φ. M ν φ can be determined from the experimental data by means of the Einstein relation D ν * =RTM ν φ (D ν * is the tracer diffusivity of element ν), or can also be obtained by using the kinetics database. δ ji is the Kronecker delta and it equals to 1 when j=i, otherwise δ ji =0. The term ∂μ ν φ /∂x j φ corresponding to the thermodynamic factor can be obtained by thermodynamic database.

Interface growth kinetics and solid fraction increment
Assuming that local thermodynamic equilibrium exists at the interface due to the influence of constitutional undercooling and curvature undercooling, the interface equilibrium temperature is characterized by the following equation neglecting the kinetic undercooling is the liquidus temperature obtained from the multi-component phase diagram at a given set of concentration (w Si L ,w Mg L ). ∆T C and ∆T R are the constitutional and curvature undercooling, respectively. In order to consider the influence of different alloying elements on the constitutional undercooling, the commonly used method is to superpose the effects of each alloying element. Therefore, the constitutional undercooling can be expressed as: in which the term ∂T L liq (w Si L ,w Mg L )/∂w i L (i denotes Si or Mg) is the slope of liquidus surface with respect to the solute element i (i.e., m i L ) at the composition ( w Si L , w Mg L ). w i L* is the local equilibrium composition at the interface. The curvature undercooling can be expressed approximately by assuming that the two principal curvatures are equal (6) in which ∆S F is the melting entropy, and κ is the local principal curvature of S/L surface. θ i is the standard spherical angle of the interface normal n, and γ(n) is the anisotropic surface energy function. The function γ(n) describing the dependence of the surface energy on normal direction n in the case of cubic symmetry is given by where γ 0 is the isotropic interfacial energy with the average Gibbs-Thomson coefficient given by . For two dimension, the detailed calculation of local curvature undercooling ∆T R can be found in Ref. [11] . For the interface curvature of a cell with solid fraction f S (0<f S <1), the counting method is used with the expression given by ( ) where Δx is the mesh size and f s i is the solid fraction of the neighboring cells, and N is total number of the first layer neighboring cells which is equal 26 in three dimension and 8 in two dimension.
Assuming that on the length scale of an interfacial cell, the concentrations are uniform, and the solid and liquid mixed well with a local equilibrium, therefore the lever rule can be adopted to calculate the increment of solid fraction. The solute conservation of each solute element during each time step for interfacial cells yields the following defining equation for solid fraction variation The analytical solution expressions for interfacial equilibrium concentrations (w Si L* , w Mg L* ) and Δf S can be directly obtained by a combination of Eq.(4)~Eq.(10).

Coupling Pandat software package
In the present model, the approach for microstructural simulation of Al-7Si-Mg ternary alloys is coupled with Pandat software package in combination with thermodynamic/kinetic/equilibrium phase diagram calculation databases to obtain the required data for simulations [16] . The atomic diffusion mobility and thermodynamic factor needed for the calculation diffusion matrix are obtained through thermodynamic and kinetic databases respectively, and equilibrium phase diagram data including liquidus temperature T L liq (w Si L ,w Mg L ), the partition coefficients (k Si ,k Mg ) and the slope of liquidus surface is generated by PanEngine [16] .

Algorithmic description for 3-D polycrystals with various crystallographic orientations
In CA models, due to the influence of mesh induced anisotropy, the traditional capture process will make the dendrite tend to grow along the axis independent of the original crystallographic orientation.
In practice, however, the growth orientations of the dendrites in the undercooled melt vary from each other, as shown in Fig.1a. In order to deal with different crystallographic orientations for polycrystalline growth in three dimension, a space coordinate conversion algorithm is adopted. As demonstrated in Fig.1a, any individual dendrite is assigned to a local coordinate system (i.e., x 0 , y 0 , z 0 axis in Fig.1), which is chosen to be parallel to the [100] directions of the dendrite arm. Therefore, the relative coordinate conversion from (x 0 , y 0 , z 0 ) to (x, y, z) can be uniquely determined by a rotation matrix (define as Rotate 1 , Rotate 2 , Rotate 3 ) which is expressed in the terms of three Euler angles (α, β, γ) as shown in Eq.11, Eq.12 and Fig.1b.
For two dimension condition (i.e., α=0, β=0), the relationship of Eq.12 can be expressed as The underlying capture approach is similar to the 3-D decentered octahedron growth or 2-D decentered square growth technique. In the paper, the 3-D dendrite capture process was analyzed in detail and the 2-D capture rule can be found in Ref. [11] . A misoriented grain with the preferential angle (α, β, γ) nucleates at the center and the cell turns into interface and begins to grow into an octahedron shape. At each time step, the evolution of distance L(t) from the center to the vertex of the octahedron is related to the increment of solid fraction Δf S , and can be calculated based on the following equation (14) Once the value L(t) is calculated, the coordinates of the six vertexes of the octahedron in the local coordinate system is known, thus their corresponding coordinates in the (x, y, z) coordinate system can be easily obtained through coordinate conversion. Once the octahedron grows large enough and the vertexes touch any of the neighboring liquid cells, the touched cell will be changed to interface state, meanwhile is assigned with the same preferential orientation as the nuclei's. A new octahedron is generated for the new interface cell and can keep growing along the same direction as the nuclei'.

Verification of the CA model
Model validation was carried out in this section by comparing the simulated results with classic analytical model for Al-7Si binary alloy. The physical parameters used in simulations were referred to the literature [17] . This simulation was carried out in a domain of 0.3mm×0.3mm×0.3mm with a cell size Δx=2μm. A single solid seed with initial composition k Si w Si 0 and the preferential crystallographic orientation (0 o , 0 o , 0 o ) was set at the center of the domain. Various constant undercooling ΔT (defined as T L liq (w Si 0 )-T * in which T * is the imposed temperature) in the range of 2K to 8K was used to simulate the dendrite growth. Fig.2 shows the simulated dendrite morphologies at the undercooling of 2K and 6K. It can be seen that the dendrite under 2K exhibits a branchless needle shape while the dendrite under 6K appears well-developed side branching, being consistent with the fact that higher undercooling tends to make the dendrite appearing side branching [18] . The secondary dendrite arms grow perpendicular to the primary arms and the initial secondary dendrite arm spacing λ 2 ，which develops immediately behind the primary dendrite tip, has a random spacing. For the dendrite tip morphology characteristics, Langer and Müller-Krumbhaarhas proposed the scaling law (Eq.15) based on the detailed numerical analysis of the wavelength of instabilities along the sides of the dendrite tip [19] λ 2 /R=2.1±0.03 (15) where R is the steady-state dendrite tip radius. In order to measure the tip radius of the simulated dendrite, the parabola fitting approaches is used, and the 3-D isosurface plot is depicted for f s =0.5. For ΔT=6K, the parabolic line is evaluated as y = 0.1434x 2 + 0.032x + 5.973 as shown in Fig3, and the tip radius is determined as R=3.49µm. The secondary dendrite arm spacing λ 2 , near the dendrite tip was approximately 8.0µm. Therefore, the value of λ 2 /R equals to 2.29, which is close to the theoretical value in Eq.15.
In order to make a comparison with the predictions of the Lipton-Glicksman-Kurz (LGK) analytical model, the steady-state tip growth velocity under different undercooling was measured. Since the selection parameter σ * varies with the anisotropy parameter ε predicted by microsolvability theory, according to the 3-D linearized solvability theory, the value of σ * for 3-D dendrite is 0.085 corresponding to ε=0.03 as used in the simulations. Fig.4 compares the simulated dendrite tip velocity with the values calculated by LGK model as a function of undercooling. As expected, the tip velocity increases with the undercooling increasing, and both of the model predictions follow the same trend with a reasonable agreement.

Solute distribution behavior and interactions between Si and Mg elements
In this section, simulations were performed for Al-7Si-0.5Mg alloy in a domain containing 200×200×200 cells with a cell size Δx=2μm, and isothermal solidification is adopted with a constant temperature T * =881K. The effect of interactions between alloying elements on diffusion coefficients was considered. Fig.5 depicts the dendrite morphologies and solute distribution both for Si and Mg elements. Since the diffusivity of solute is much smaller than the dendrite growth velocity, the solute discharged from the newly solidified solid phase can not sufficiently diffuse to the bulk liquid. Therefore, it can be clearly seen from the Fig.5 that both Si and Mg elements are piled up near the dendrite interface as a result of solute redistribution. It is found that the maximum concentration of Si in the region between the dendrite arms is approximately 7.9wt.% while the value for Mg is 0.54wt.%. If defining (w i L -w i 0 )/w i 0 as the degree of solute enrichment for i element, the value for Si is 12.8% which is greater than that of Mg (8%). This can be accounted by the reason that the liquid diffusion coefficient of Mg is approximately two times of the diffusion coefficient of Si, i.e., D Mg L ≈ 2D Si L , which means that the accumulated Mg element ahead dendrite interface is easier to diffuse out. Furthermore, the difference in equilibrium partition coefficient (k Si ≈0.13, k Mg ≈0.24) also contributes to this phenomenon. The larger equilibrium partition coefficient results in the dissolution of more solutes by the solid phase, and less solutes discharged from the newly solidified solid phase into the liquid phase.   6 shows the diffusion matrix in 2-D slice both for diagonal (D sisi and D mgmg ) and off-diagonal (D simg and D mgsi ). It is found that in the liquid region, the diffusion coefficients are concentration dependent, and the off-diagonal terms are negative which means the attractive interaction between the alloying solutes. It can be seen that the values of off-diagonal terms (D simg and D mgsi ) and diagonal term (D mgmg ) at the interface is smaller than that far from the interface, while the diagonal term (D sisi ) at the interface is larger than that at other positions. Fig.7 represents the diagonal diffusion coefficients of Al-7Si-0.5Mg alloy varying with the concentration of Si and Mg at the isothermal temperature of 881K. As is shown, the diffusion coefficient D mgmg decreases with the increasing of Si and Mg. For diffusion coefficient D sisi , it increases approximately from 2.385×10 -9 m 2 /s to 2.405×10 -9 m 2 /s as the Si varies from 7wt.% to 8wt.%, while it decreases from 2.385×10 -9 m 2 /s to 2.380×10 -9 m 2 /s with the Mg increasing from0.5wt.% to 0.8wt.%. Therefore, it can be seen that the diffusion coefficient D sisi is more sensitive to the change of Si composition and that is why the value of D sisi near the dendrite interface is larger than that far from the interface in Fig.6a. In Fig.6, we can also clearly see that the off-diagonal diffusion coefficients are smaller than that of diagonal term by an order of magnitude. To illustrate the influence of solute interaction on dendrite growth, Fig.8 shows the comparison of the solute concentrations at the dendrite tip in the liquid as a function of dendrite growth time by taking or without taking the solute interaction into consideration. It is interesting to note that the concentration of Si is a little lower in the case with solute interaction than that without solute interaction, while for Mg element, the concentration at dendrite tip is higher with consideration of solute interaction than the case without solute interaction. A brief analysis of this phenomenon can be described as: attractive interactions tend to drive both of the Si and Mg elements toward the solid-liquid interface by the concentration gradient of the other solute, slowing down the diffusion behavior from the interface to the bulk liquid, which tends to make solute concentration at the dendrite tip higher when considering solute interaction. The higher solute concentration at the dendrite tip when considering solute interaction leads to a decrease in the undercooling and a further decrease in the tip growth velocity. Through measuring the average growth velocity during the solidification time of 0.2s~0.7s, the value in the condition of solute interaction is 183µm/s and in the condition of no solute interaction, it is 191µm/s. Hnece the lower growth rate when the effects of solute interaction is taken into consideration makes it more time available for the solute ahead of the dendrite to diffuse out. Besides, due to the composition of Si is much higher than Mg in the Al-7Si-0.5Mg alloy, the contribution of cross diffusion term (D simg ) to Si solute distribution is much smaller compared with the contribution of D mgsi to the distribution of Mg solute. Both of the reasons result in the difference of the solute distribution behavior of Si and Mg observed in Fig.8.  Since most practical industrial casting processes of Al-7Si-Mg casting alloy mainly involve equiaxed polycrystalline growth in random crystallographic orientations, this case is simulated using this model for Al-7Si-0.36Mg (wt.%) alloy. The temperature profiles for simulation were directly obtained from the measured cooling curves as shown in Fig.9, and the temperature was assumed to be uniform in the whole calculation domain. The corresponding experimental procedure as well as results including cooling curves and metallographical microstructure will report in another paper. Two dimensional simulations were performed in a domain with the size of 7.2×7.2mm 2 . Fig.10 shows the simulated final dendrite microstructure of sample 2, sample 4, sample 6, and sample 7. It can be seen that with the increasing of cooling rate, the nucleation density increases and the secondary dendrite arm spacing (SDAS) decreases, both of which result in the refinement of the microstructure.  seen that the numerical results from the model agree well with the experimental ones, and the fitted formula, i.e., SDAS=42.8R c -0.302 , is quite in agreement with the experimental result as reported in Ref. [20] that SDAS in A356/A357 alloys fits well the empirical equation SDAS=39.4R c -0.317 , proving the reliability of present model in predicting SDAS. Figure 11. Comparison between the predicted secondary dendrite arm spacing and the experimental ones solidified in different cooling conditions.

Primary dendrite spacing prediction of columnar dendrite
Columnar dendrite is the most frequently observed microstructure in directional solidified products, and the primary dendrite spacing, as an important characteristic microstructural length scale, determining the segregation of alloying elements as well as the size and distribution of eutectic particles and solidification defects, shows a significant effect on the mechanical properties of solidified products. Therefore, it is of particular importance to investigate the columnar dendrite evolution and to give a precise prediction of primary dendrite spacing. For a given solidification history, the dendrite spacing mainly controlled by the temperature gradient G and dendrite growth velocity v. Hence, both experimental and simulation investigation on dendrite spacing evolution for Al-7Si-0.36Mg alloy during directional solidification is carried out. The experimental results will published in another paper. The simulations were performed in a domain of 6.25mm×15mm with a mesh size of 5μm. At the beginning of the simulation, a certain amount of crystal seeds with <100>directions aligned with the coordinate axes were nucleated at the bottom of the domain. Fig.12 presents the simulated steady state columnar dendrite arrays for v=200μm/s, v=150μm/s, v=100μm/s, v=50μm/s at a given temperature gradient G=15K/mm. It can be seen that during initial growth stage, developed tertiary arms emitted from the secondary arms are observed as a result of too large initial dendrite spacing λ ini , i.e., λ ini exceeds the upper limit λ max . However, most of the initial tertiary arm spacing between adjacent dendrite arms is lower than the lower limit λ min , and competitive growth occurs by blocking some tertiary arms, and some become the stable primary dendrite arms. After a period of competitive growth, the dendrite array reaches a stable state without dendrite spacing varying. Fig.13 shows the comparison between the CA prediction, analytical model predictions [21~23] , and the experimental results for the average primary dendrite spacing λ ave . It is observed that CA model prediction show quite good agreement with the experimental results, which is much better than the agreement between analytical models' predictions and experiments. The comparisons reveal that the present CA model has a high accuracy in predicting primary dendrite spacing in directional solidification.

Conclusions
A cellular automaton model for the prediction of the ternary dendrite growth in two and three dimensions during solidification processes is developed. The model was coupled with thermodynamic/kinetic/equilibrium phase diagram calculation databases incorporated in the Pandat software package and the solute interactions between alloying elements was considered. To eliminate the influence of mesh-induced anisotropy on dendrite growth, a modified decentered octahedron for neighborhood tracking is introduced. First, detailed model validations have been performed with the simulated steady-state tip parameters being in good agreement with the LGK predictions and the wavelength of instabilities theory. The model is then applied to Al-7Si-0.5Mg dendrite simulation for investigation of Si and Mg solute distribution behavior and the effects of solute interactions on dendrite growth. It shows that Si element is easier to segregated due to smaller diffusion coefficient and smaller equilibrium partition coefficient compared with Mg element. The diffusion coefficient is concentration dependent and interactions between Si and Mg elements can affect the dendrite growth behavior. Finally, this model is applied to simulation of equiaxed dendrite growth and columnar dendrite growth for Al-7Si-0.36Mg alloy in practical solidification process. The predictions of secondary dendrite arm spacing and primary dendrite spacing shows a reasonable agreement with the experimental results. The simulated results effectively demonstrated the abilities of this model in the prediction of dendritic microstructure in ternary alloys.