Comparison of Cellular Automaton and Phase Field Models to Simulate Dendrite Growth in Hexagonal Crystals

J. Mater. Sci. Technol., 2012, 28(2), 137–146.

Comparison of Cellular Automaton and Phase Field Models to Simulate Dendrite Growth in Hexagonal Crystals Mohsen Asle Zaeem1,2)† , Hebi Yin3) and Sergio D. Felicelli1,2) 1) Center for Advanced Vehicular Systems, Mississippi State University, Starkville, MS 39759, USA 2) Mechanical Engineering Department, Mississippi State University, Starkville, MS 39759, USA 3) Materials Processing Group, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA [Manuscript received September 13, 2011, in revised form December 7, 2011]

A cellular automaton (CA)–finite element (FE) model and a phase field (PF)–FE model were used to simulate equiaxed dendritic growth during the solidification of hexagonal metals. In the CA–FE model, the conservation equations of mass and energy were solved in order to calculate the temperature field, solute concentration, and the dendritic growth morphology. CA–FE simulation results showed reasonable agreement with the previously reported experimental data on secondary dendrite arm spacing (SDAS) vs cooling rate. In the PF model, a PF variable was used to distinguish solid and liquid phases similar to the conventional PF models for solidification of pure materials. Another PF variable was considered to determine the evolution of solute concentration. Validation of both models was performed by comparing the simulation results with the analytical model developed by Lipton–Glicksman–Kurz (LGK), showing quantitatively good agreement in the tip growth velocity at a given melt undercooling. Application to magnesium alloy AZ91 (approximated with the binary Mg–8.9 wt% Al) illustrates the difficulty of modeling dendrite growth in hexagonal systems using CA–FE regarding mesh-induced anisotropy and a better performance of PF–FE in modeling multiple arbitrarily-oriented dendrites growth. KEY WORDS: Dendrite growth; Cellular automaton; Phase-field model; Finite element; Magnesium alloy

1. Introduction

tion such as phase field (PF)[3–17] and cellular automaton (CA)[18–30] methods.

Dendrite growth during solidification process forms the microstructure and has significant effects on the material properties of cast alloys[1] . The multi-scale nature of solidification patterns, which depend on temperature distribution, solute concentration, capillary forces, and kinetic length[2] , makes the understanding and prediction of these patterns exceptionally difficult. In addition to the experimental solidification research[1,2] , recent numerical simulations have emerged as a powerful tool in predicting microstructural evolution during solidification. Different mathematical approaches have been employed to investigate the dendritic growth during solidifica-

The originally proposed PF models by Cahn and Hilliard[31] for concentration evolution in binary systems and later by Allen and Cahn[32] for phase boundary motions were adopted by Kobayashi[3] for simulating dendrite grain growth in pure materials by considering a non-conserved PF variable as the marker for solid and liquid phases. Other thermodynamically consistent PF models for dendritic solidification of pure materials somehow recover Kobayashi s model[4] . During the last two decades, different PF models have been developed for solidification of binary[9–13] and ternary[14,15] alloy systems. Warren and Boettinger[9] synthesized the entropy formulation by Wang et al.[4] and the isothermal energy formulation by Wheeler et al.[10] to develop a PF model for dendritic solidification for binary alloys. Beckermann et al.[12] presented

† Corresponding author. Assist. Prof., Ph.D.; Tel.: +1 662 3250126; Fax: +1 662 3255433; E-mail address: [email protected] (M. Asle Zaeem).

138

M. Asle Zaeem et al.: J. Mater. Sci. Technol., 2012, 28(2), 137–146.

a PF model for solidification in both pure and binary systems involving convection in the liquid phase. Kobayashi et al.[14] developed a multi-component PF model for dendrite solidification in a ternary alloy, which later was used by Wang and Yang[15] to study the effects of solid fraction on the isothermal dendritic coarsening behaviors and mechanisms in a Ni–Al–Nb ternary system. Simultaneously, different numerical methods have been developed and used to solve the governing evolution equations of dendrite solidification in PF models, which are mostly based on finite difference (FD) and FE methods[3,7,9–17] . Although the PF method seems to be a reliable computational approach in the area of dendritic solidification, it requires significant computer resources and for large scale domains, using supercomputers and parallel computing methods is inevitable. The CA method is another computational approach to simulate the dendritic growth during solidification. Combined with FE or FD, this method produces results similar to those of the PF method by obtaining the temperature and solute fields and then determining the solid/liquid (S/L) interface. Rappaz and Gandin[18,19] developed a two dimensional probabilistic model to simulate the dendritic grain formation during solidification based on the CA technique. They coupled the CA algorithm and FE method to obtain the thermal field and the microstructure of an Al–Si alloy, but this work did not provide the details of the growing process of dendritic grain. Sanchez and Stefanescu[20] developed a dendrite growth model, which proposed a solution for the artificial anisotropy induced by the CA grid, but the dendrite could grow only aligned with the mesh or in a 45 deg. orientation. Then an improved model[21,22] was proposed by introducing a new virtual front tracking method, which was able to simulate dendrites growing in any preferential orientation. Zhu and Hong[23] developed a CA model to simulate the solidification microstructure of both eutectic and hypoeutectic Al–Si alloy and provided good insight into the eutectic nucleation and growth behaviors. Later, a CA model, which considers the influence of fluid flow on the dendrite growth was also developed[24,25] . In addition to these two-dimensional simulations, three-dimensional models[26–28] were also reported, which combined the CA and FE methods to simulate the dendritic growth in binary alloys controlled by solute diffusion. Difficulties associated with CA method are a result of the artificial anisotropy introduced by the CA mesh. Although different alternatives have been adopted to solve this problem for cubic crystal structures, no practical solutions have been found for hexagonal materials. Most of the numerical studies in the area of dendritic solidification focus on growth of cubic crystal materials (FCC/BCC) during solidification[9,11–28] . Bottger et al.[33] used the PF method to simulate the microstructure during the equiaxed solidification of

hexagonal crystal materials, particularly the magnesium alloy AZ31, but their results did not show the single dendritic morphology with six-fold symmetry. Liu et al.[34] simulated the two-dimensional dendrite profiles by employing a mathematical construction method. But this mathematical description method of the grain contour did not consider the tertiary and above arm branching. In this paper, two models, one based on CA and the other based on PF were developed to simulate the dendritic growth of hexagonal metals during the solidification process. In particular, for the case of dendrite growth in a binary Mg–8.9 wt% Al alloy, a simulated microstructure with perfect six-fold symmetry was obtained in both CA and PF models. Although the CA model is capable of predicting the six-fold geometry of the HCP crystal structure, in the present version the dendrite growth is limited to the grid orientation and arbitrary orientations are not currently possible. On the other hand, the PF model does not have these limitations and is capable of simulating multiple arbitrarily-oriented dendrites. 2. Description of Models The dendrite growth for Mg–8.9 wt% Al alloy, which is similar in composition to alloy AZ91, was simulated. The properties of this alloy in the simulations are listed in Table 1[35–37] . 2.1 Cellular automaton model A single nucleus with the preferential orientation of 0 deg. is placed at the center of the calculation domain (100 μm×100 μm) discretized with a 200×200 CA grid and a 20×20 regular FE mesh of quadrilateral bilinear elements. The calculated morphology and composition field for a single equiaxed dendrite is shown in Fig. 1. It is observed that there are six primary arms, but the angle between the primary arms is not 60 deg. Some arms grow aligned with the axis of the mesh, and others grow at 45 deg. They should grow with perfect six-fold symmetry, without hindrance from other grains or walls. Fu et al. also experienced a similar problem, as indicated in literature[38] , and they attributed this problem to the methodology of defining neighboring cells. As evidenced in the previous example, the CA method used to simulate dendrite evolution has a disadvantage: grid dependent anisotropy[39,40] , which is that the grain growth process is very sensitive to the mesh shape and size. For the simulation of dendrite growth with cubic crystal structure materials, meshes of rectangular elements have been widely used, and hexagonal meshes have also been introduced by some articles using front-tracking techniques[41,42] . Based on the fact that Mg has an HCP crystal structure with six-fold symmetry and previous reports that the use of hexagonal elements seems to reduce mesh-induced

M. Asle Zaeem et al.: J. Mater. Sci. Technol., 2012, 28(2), 137–146.

Table 1 Properties of Mg alloy used in calculations Property Thermal expansion coefficient (βT ) Density of liquid (ρl ) Density of solid (ρs ) Diffusivity of alloy elements in liquid (Dl ) Diffusivity of alloy elements in solid (Ds ) Thermal conductivity in liquid (λl ) Thermal conductivity in solid (λs ) Average specific heat of liquid (cl ) Average specific heat of solid (cs ) Latent heat of fusion (L) Liquidus temperature (TL ) Solidus temperature (TS ) Gibbs–Thomson coefficient Equilibrium liquidus slop (ml )

139

[35,36,44]

Value −2.6×10−5 K−1 1650 kg·m−3 1750 kg·m−3 5.0 × 10−9 m2 ·s−1 5.0 × 10−13 m2 ·s−1 80 J·K−1 ·m−1 ·s−1 105 J·K−1 ·m−1 ·s−1 1350 J·kg−1 ·K−1 1200 J·kg−1 ·K−1 3.7 × 105 J·kg−1 868 K 743 K 2.0 ×10−7 K·m −5.69 K/wt%

outward direction normal to the boundary. The FE method is employed to solve Eqs. (1) and (2) and then to obtain the temperature field. 2.1.2 Solute distribution calculation The solute conservation in solid and liquid phases is obtained by solving the governing equations for each phase separately as illustrated below: ∂Ci ∂fs = Di · ∇2 Ci + Ci · (1 − k) ∂t ∂t

Fig. 1 Equiaxed dendrite of AZ91 Mg alloy simulated with CA–FE method with finite element mesh of quadrilateral bilinear elements, showing problem of mesh-induced anisotropy [40,43]

anisotropy , an FE mesh consisting of hexagonal elements was developed to address this problem with the CA technique[44] . 2.1.1 Temperature field calculation The two-dimensional transient differential equation governing the heat transfer within the calculation domain is given by: ∂T L ∂fs = α · ∇2 T + ∂t ρCP ∂t

(1)

where T is temperature, t is time, α is thermal diffusivity, L is latent heat of solidification, ρ is density, CP is specific heat, and fs is volume solid fraction. A forced boundary condition, heat flux q, is prescribed on the four walls (equiaxed grain growth simulation) or single wall (columnar grain growth simulation) as shown below: ∂T (2) q = −λ ∂n where λ is the thermal conductivity, and n is the

(3)

where C is concentration of solute, D is solute diffusivity, the subscript i indicates solid or liquid, and k is the equilibrium partition coefficient (k=0.2 in this work). At the solid/liquid interface, the solute partition between liquid and solid is calculated by Cs∗ = k · Cl∗

(4)

where Cs∗ and Cl∗ are interface solute concentrations in solid and liquid phases, respectively. The FE method is employed to solve Eq. (3), using the same CA grid as an FE mesh of bilinear elements. A zeroflux boundary condition is applied to the cells located at the boundaries of the calculation domain. 2.1.3 Kinetics parameters for the CA model The solute distribution ahead of the interface is used as the driving force to simulate the dendrite growth in the CA model. The process of dendrite growth is predominantly controlled by the difference between the local interface equilibrium solute concentration and the local actual liquid solute concentration. The changing rate of solid fraction determines the velocity and morphology of dendrite growth. Based on the calculation of the local actual liquid concentration Cl from Eq. (3) and the interface equilibrium composition Cl∗ , the increase of solid fraction, Δfs , at the interface cells can be obtained as: Δfs = (Cl∗ − Cl )/(Cl∗ · (1 − k)) (5)

M. Asle Zaeem et al.: J. Mater. Sci. Technol., 2012, 28(2), 137–146.

140

The interface equilibrium composition is calculated by: T ∗ − Tleq + ΓKf (ϕ, θ0 (6) Cl∗ = C0 + ml where C0 is the initial solute concentration, Tleq is the equilibrium liquidus temperature at the initial solute concentration, ml is the liquidus slope, Γ is the GibbsThomson coefficient, and K is the curvature of the solid/liquid interface. The function accounting for the anisotropy of the surface tension is denoted by f (ϕ, θ0 ), where ϕ is the growth angle between the normal to the interface and the x-axis, and θ0 is the angle of the preferential growth direction with respect to the x-axis. The interface equilibrium temperature calculated by Eq. (1) is denoted by T ∗ . For hexagonal crystals, the function f (ϕ, θ0 ) exhibits a six-fold anisotropy[33,45] : f (ϕ, θ0 ) = 1 − δ cos[6 · (ϕ − θ0 )]

(7)

where δ is the anisotropy coefficient (δ=0.06 in this work), and the growth angle can be calculated from the following equation: ⎧   ∂fs s /∂x ⎨ ≥0 cos−1 ((∂fs /∂x)∂f 2 +(∂f /∂y)2 )1/2 s   ∂y ϕ= ∂fs /∂x ∂fs ⎩ 2π − cos−1 ∂y < 0 ((∂fs /∂x)2 +(∂fs /∂y)2 )1/2 (8) The interface curvature of a cell with solid fraction fs can be obtained by counting the nearest and second nearest neighboring cells[46] :

K=

1 1−2 a

fs +

N  j=1

fs (j)

N +1

(9)

where a is the length of the CA cell side, N is the number of the nearest and second nearest neighboring cells (N =19), and fs (j) is the solid fraction of neighboring cells. Different time steps, ΔtT and ΔtC , are used for the calculation of heat transfer and mass transfer, reρC (m·a)2 a2 . spectively: ΔtT = p4.5λ , and ΔtC = 4.5D 1 Since the heat diffusivity is much faster than the solute diffusivity, Eq. (1) is solved Nt (Nt = ΔtC /ΔtT ) times with time step ΔtC , per each solution of Eq. (3) with time step ΔtC in order to obtain converged temperature and solute concentration fields. 2.1.4 Rules of capturing interface cells Because an HCP crystal material grows with a sixfold symmetry, a hexagon-shape solid seed is initially placed at the center of a cell and the seed is allowed to grow along its diagonals. Once the corners of the hexagon seed reach any of the six neighboring cells, the neighbor cell will be changed to an interface cell. A new hexagon seed, having the same preferential crystallographic orientation as the original cell, is generated and placed at the center of the new interface

cell. The new hexagon seed starts to grow according to the change of the solid fraction in the new interface cell. After the original hexagon has changed all the neighboring cells into interface cells, the original cell continues to grow until its solid fraction becomes unity, after which the state of the original cell becomes solid and changes any surrounding liquid cells into interface cells. The length of the diagonal (d) of the seed is calculated based on the solid fraction fs d = kd · fs · a/ cos(θ0 )

(10)

In where kd is the diagonal length coefficient. this work, kd is equal to 0.962 for all the CA–FE simulations[44] . 2.1.5 Numerical implementation procedures The simulation of the grain growth process can be summarized in the following steps: (1) Solve Eq. (1) to obtain the temperature field in the domain using time step ΔtT for Nt iterations with constant heat flux as boundary condition; (2) Interpolate in each finite element to obtain the temperature for all cells inside the element and then solve Eq. (3) with time step ΔtC to obtain the solute field in the whole domain; (3) Calculate the increase of solid fraction at the interface cells; (4) Update the thermal field and solute fields based on the release of latent heat and the rejection of solute during solidification; (5) Use the rules indicated in Section 2.1.4 to capture the new interface cells. A detailed description of the CA–FE model is presented in literature[44] . 2.2 Phase field model A PF variable, φ(x, t), is considered in the spatial (x) and time (t) domains similar to the Kobayashi s model[3] . φ(x, t)=0 represents the liquid phase and φ(x, t) = 1 represents the solid phase. For a binary alloy, the Gibbs–Thomson equation for an isotropic surface energy can be written as[12] : 1 ∂φ = Tm − T + ml Cl − Γ∇ · n μ|∇φ| ∂t

(11)

where μ is a linear kinetic coefficient, T is the temperature, Tm is the equilibrium melting temperature, ml is the slope of the liquids line from an equilibrium phase diagram, Cl is the equilibrium concentration in the liquid, Γ is the Gibbs–Thomson coefficient, and n is the unit normal vector exterior to the solid: ∇φ (12) |∇φ| A double-well potential is considered for the Gibbs free energy[12] : n 1 ) (13) φ = (1 − tanh 2 2ω n=−

M. Asle Zaeem et al.: J. Mater. Sci. Technol., 2012, 28(2), 137–146.

where n is the coordinate normal to the interface and 6ω is the diffuse-interface thickness (interface range) where φ varies between 0.05 to 0.95. PF method depends on a diffuse-interface separating the neighboring phases, and it can be shown that this method reduces to the standard sharp-interface approach in the limit of vanishing diffusive-interface thickness[47] . It is imperative to understand how the quality of the solution varies based on the diffuse-interface thickness, because in numerical simulations the grid spacing needs to be of the order of or smaller than this thickness to generate valid results[48] . It was shown by different studies that the interface thickness must be smaller than the capillary length for the solution to converge to the sharp-interface limit[4,10] . Also the PF interface thickness only needs to be small compared to the “mesoscale” of the heat and/or solute diffusion field, and a finite PF interface thickness can satisfy the standard interface conditions[7] . Similarly, we used a PF approach where the diffuse-interface thickness only needs to be small compared to the mesoscale of the diffusion and temperature fields[12] . Using Eqs. (11), (12) and (13), and after some manipulations, the evolution equation for φ becomes[12] :  ∂φ φ(1 − φ)(1 − 2φ)  = μΓ ∇2 φ − + ∂t ω2 φ(1 − φ) (14) ω To include anisotropy in the system, we used a modified version of Eq. (14) similar to the Kobayashi s model[3] :  ∂φ ∂φ ∂   ∂φ ∂ = μΓ − σσ + (σσ  )+ ∂t ∂x ∂y ∂y ∂x μ(Tm − T + ml Cl )

∇ · (σ 2 ∇φ) −

φ(1 − φ)(1 − 2φ)  + ω2

φ(1 − φ) μ(Tm − T + ml Cl ) ω

141

The difference in density of liquid (ρl ) and solid (ρs ), specific heat of liquid (Cpl ) and solid (Cps ), and thermal conductivity in liquid (λl ) and solid (λs ) are included similar to Eq. (18): ρ˜ = ρs + (ρl − ρs )

1−φ 1 − φ + kφ

Cp = Cps + (Cpl − Cps ) λ = λs + (λl − λs )

1−φ 1 − φ + kφ

1−φ 1 − φ + kφ

(20)

(21)

(22)

In the current PF model, we have not used an antitrapping term that can overcome the possible excess solute trapping; although this may affect the calculated solute maps, it may not have a noticeable effect on the tip growth velocity. Certainly to have more accurate solute maps, especially in cases with higher undercoolings or rapid solidification (solidification time less than 10−4 s), using an anti-trapping current term in the phase-field model is inevitable[49,50] . In previous works, we developed a mixed order FE model for the coupled PF equation and elasticity equations for solid state phase transformations[48,51–53] . Although in the present work all the equations in the PF model are 2nd order partial differential equations, a similar algorithm is adopted. A uniform mesh of square elements with quadratic interpolation functions was used to discretize the domain for Eqs. (15), (17) and (19). To have relatively reliable results, it was shown that at least 5 FE elements have to be considered along the PF interface thickness in the PF model[48] . The mesh size is determined to be: (6ω)/5=0.2 μm. COMSOL Multiphysics software is used for the computations[54] . 3. Simulation Results

(15)

3.1 Model validation

where: σ(θ) = 1 − δ cos[6(θ − θ0 )];

θ = tan−1

 ∂φ/∂y

∂φ/∂x (16) Another PF variable, C(x, t), is introduced, which represents the concentration of species. The evolution equation for the concentration is[12] :   ∂C ˜ ∇C + (1 − k)C ∇φ =∇·D ∂t 1 − φ + kφ where ˜ = Ds + (Dl − Ds ) D

1−φ 1 − φ + kφ

(17)

(18)

The evolution of the temperature field is similar to Eq. (1): ∂T L ∂φ = α · ∇2 T + (19) ∂t ρ˜Cp ∂t

The relationship between the primary and secondary dendrite arm spacing (SDAS) with the cooling rate has been amply documented and empirical relations have been proposed for several alloy systems[55,56] . Using our developed CA model, we simulated the growth of an individual dendrite under several intensities of cooling rate. The SDAS was then recorded for different values of cooling rate and plotted in Fig. 2. The experimental data on SDAS vs cooling rate previously reported by Caceres et al.[57] , Labrecque et al.[58] , Dube et al.[59] and Sequeira et al.[60] are also included in Fig. 2, showing a reasonable agreement with our calculated values. Lipton et al.[61] developed an analytical model (the LGK model), which describes free dendrite growth at a given melt undercooling. The tip growth velocity with various temperature undercoolings calculated by

142

M. Asle Zaeem et al.: J. Mater. Sci. Technol., 2012, 28(2), 137–146.

Fig. 2 SDAS vs cooling rate for alloy AZ91 calculated by the present CA model and comparison with the data from literature[57–60] . The square calculation domain has a 400×400 CA grid with side lengths of 0.5 μm

Fig. 3 Tip growth velocity vs temperature undercooling calculated by the present CA and PF models and comparison with the LGK theory[61] . For CA simulations, the square domain has a 400×400 CA grid with side lengths of 0.5 μm, and for PF simulations mesh size of 0.2 μm are used to discretize the domain

Fig. 4 Solute map at different simulation time using CA–FE model with preferred growth orientations of zero degree with respect to the horizontal direction: (a) 0.0212 s, (b) 0.0424 s, (c) 0.0636 s. The square domain has a 400×400 CA grid with side lengths of 0.5 μm, and the heat flux imposed at the four walls was 10 kW·m−2 . Colorbar shows the solute concentration in wt%

M. Asle Zaeem et al.: J. Mater. Sci. Technol., 2012, 28(2), 137–146.

143

Fig. 5 Solute map at different simulation time using PF–FE model with preferred growth orientations of zero degree with respect to the horizontal direction: (a) 0.02 s, (b) 0.04 s, (c) 0.06 s. The calculation domain is 200 μm×200 μm with mesh size of 0.2 μm, and the heat flux imposed at the four walls was 10 kW·m−2 . Colorbar shows the solute concentration in wt%

Fig. 6 Dendrite growth from a single nucleus with preferred growth orientations of 45 deg. with respect to the horizontal direction presented using the PF–FE model. Times: (a) 0.02 s, (b) 0.04 s, (c) 0.07 s. The calculation domain is 400 μm×400 μm with mesh size of 0.2 μm, and the heat flux imposed at the four walls was 15 kW·m−2 . Colorbar shows the solute concentration in wt%

144

M. Asle Zaeem et al.: J. Mater. Sci. Technol., 2012, 28(2), 137–146.

Fig. 7 Phase-field, temperature, and solute maps of multiple arbitrarily-oriented seeds at different simulation time using the PF–FE model. Heat flux imposed at the four walls was 10 kW·m−2 . The calculation domain is 400 μm×400 μm with mesh size of 0.2 μm

the LGK theory is shown in Fig. 3, alongside with values calculated by our CA and PF models. To measure the tip growth velocity, we counted the number of cells captured by the grain at different time, and since the cell size is known, the speed could be obtained. It is observed that our simulation results follow rather well the LGK predictions. We also observed that increasing the cooling rate enhances the branching of the dendrite arms, which is in agreement with the calculations reported by Bottger et al.[33] using the PF method. 3.2 CA vs PF A single nucleus was set at the calculation domain center to simulate the grain growth process during solidification. With constant heat flux (10 kW/m2 ), the calculation domain has uniform initial temperature and composition. The nucleus has an initial composition kC0 and preferred growth orientations

of zero degree with respect to the horizontal direction. In the CA–FE model, the square domain has a 400×400 CA grid with side lengths of 0.5 μm. A uniform mesh of hexagonal elements was used for the FE discretization. Fig. 4(a), (b), and (c) show the simulated evolution of equiaxed dendrite growth at different simulation time of 0.0212, 0.0424, and 0.0636 s, respectively. It can be seen that in the early stages of solidification, dendrites develop primary arms, which follow the crystallographic orientations as shown in Fig. 4(a). As solidification proceeds, the primary arms become larger and the secondary arms begin to occur (Fig. 4(b)). With further solidification, some tertiary dendritic arms form on the secondary arms (Fig. 4(c)). Equiaxed dendrite growth from a single nucleus similar to Fig. 4 is simulated using the PF–FE model in Fig. 5. Uniform square elements with quadratic interpolation functions and with mesh size of 0.2 μm

M. Asle Zaeem et al.: J. Mater. Sci. Technol., 2012, 28(2), 137–146.

are used to discretize the domain of 200 μm×200 μm. The results at different simulation times of 0.02, 0.04, and 0.06 s are shown in Fig. 5(a), (b), and (c), respectively. The results of the dendrite solidification of a single nucleus using CA–FE and PF–FE models are very similar in terms of the tip growth velocity with various undercoolings as presented in Fig. 3. From Figs. 4 and 5, it can be observed that in both the CA and PF simulations of dendrite growth from a single nucleus with preferred growth orientations of zero degree with respect to the horizontal direction, the angles between the primary arms are exactly 60 deg. The results, however, show some differences in the degree of branching, with the CA–FE model producing more ramifications and thinner secondary and tertiary arms than the PF–FE model. Also slight differences can be observed in solute maps simulated by CA–FE and PF–FE models. The simulated structure predicts qualitatively well the measured microstructure of an AZ91 dendrite[62] . In Fig. 6, the dendrite growth from a single nucleus with preferred growth orientations of 45 deg. with respect to the horizontal direction is presented using the PF–FE model. The CA–FE model was unable to obtain meaningful results for orientations other than zero degree due to mesh inducedanisotropy. In Fig. 7, the PF–FE simulation results of the solidification of multiple arbitrarily-oriented dendrites are plotted, which shows phase-field, temperature, and concentration maps at different simulation time. The results show that the PF–FE model can successfully simulate the exact 60 deg. angle between the six-primary arms of the HCP structure with arbitrary preferred growth directions. 4. Conclusion A CA–FE model and a PF–FE model were used to simulate the solute controlled dendrite growth in hexagonal crystal metals. Both the CA–FE and PF– FE models performed reasonably well in simulations of the dendrite growth from a single nucleus with preferred growth orientations of zero degree with respect to the horizontal direction. Both models were able to predict the six-fold symmetry of dendritic growth in magnesium alloy AZ91 (approximated with the binary Mg–8.9 wt% Al). However, the CA–FE model requires an FE mesh of hexagonal elements, which makes the method rather impractical. Even with the hexagonal mesh, the CA–FE model was not able to avoid the mesh-induced anisotropy problems when simulating arbitrarily oriented dendrites. On the other hand, the PF–FE model produced good results for single and multiple dendrites growing in any orientation, unaffected by mesh anisotropy issues. The sensitivity to mesh-induced anisotropy exhibited by CA–FE models in simulating HCP dendrite structures has been reported in previously published articles[29,30] . The PF–FE method seems to be a ro-

145

bust alternative tool for simulating multiple arbitraryoriented dendrites with six-fold symmetry.

Acknowledgements This work was supported by the National Science Foundation (USA) through Grant No. CBET-0931801 and the Department of Energy (USA) through cooperative agreement No. DE-FC-26-06NT42755. The authors also appreciate the sponsorship of the Center for Advanced Vehicular Systems (CAVS) of Mississippi State University. REFERENCES [1 ] D. Askeland: The Science of Materials and Engineering, Thompson, Brooks/Cole, Belmont, CA, 2003. [2 ] R. Trivedi and W. Kurz: Int. Mater. Rev., 1994, 39(2), 49. [3 ] R. Kobayashi: Physica D, 1993, 63, 410. [4 ] S.L. Wang, R.F. Sekerka, A.A. Wheeler, B.T. Murray, S.R. Coriell, R.J. Braun and G.B. McFadden: Physica D, 1993, 69, 189. [5 ] J.A. Warren, W.J. Boettinger, C. Beckermann and A. Karma: Annu. Rev. Mater. Res., 2002, 32, 163. [6 ] X. Tong, C. Beckermann, A. Karma and Q. Li: Phys. Rev. E, 2001, 63, 061601. [7 ] A. Karma and W.J. Rappel: Phys. Rev. E, 1998, 57, 4323. [8 ] I. Steinbach: Model. Simul. Mater. Sci. Eng., 2009, 17, 073001. [9 ] J.A. Warren and W.J. Boettinger: Acta Metall. Mater., 1995, 43, 689. [10] A.A. Wheeler, W.J. Boettinger and G.B. McFadden: Phys. Rev. A, 1992, 45, 7424. [11] N. Provatas and M. Greenwood: Int. J. Modern Phys. B, 2005, 19, 4525. [12] C. Beckermann, H.J. Diepers, I. Steinbach, A. Karma and X. Tong: J. Comput. Phys., 1999, 154, 468. [13] N. Ofori-Opoku and N. Provatas: Acta Mater., 2010, 58, 2155. [14] H. Kobayashi, M. Ode, S.G. Kim, W.T. Kim and T. Suzuki: Scripta Mater., 2003, 48, 689. [15] J. Wang and G. Yang: Acta Mater., 2008, 56, 4585. [16] Y. Xu, J.M. McDonough and K.A. Tagavi: J. Comput. Phys., 2006, 218, 770. [17] J.C. Ramirez, C. Beckermann, A. Karma and H.J. Diepers: Phys. Rev. E, 2004, 69, 051607. [18] M. Rappaz and C. Gandin: Acta Metall. Mater., 1993, 41, 345. [19] C.A. Gandin and M. Rappaz: Acta Metall. Mater., 1994, 42, 2233. [20] L.B. Sanchez and D.M. Stefanescu: Metall. Mater. Trans. A, 2003, 34, 367. [21] L.B. Sanchez and D.M. Stefanescu: Metall. Mater. Trans. A, 2004, 35, 2471. [22] M.F. Zhu and D.M. Stefanescu: Acta Mater., 2007, 55, 1741. [23] M.F. Zhu and C.P. Hong: Metall. Mater. Trans. A, 2004, 35, 1555. [24] M.F. Zhu, S.Y. Lee and C.P. Hong: Phys. Rev. E, 2004, 69, 061610-1-12. [25] Y.H. Shin and C.P. Hong: ISIJ Int., 2002, 42, 359.

146

M. Asle Zaeem et al.: J. Mater. Sci. Technol., 2012, 28(2), 137–146.

[26] W. Wang, P.D. Lee and M. McLean: Acta Mater., 2003, 51, 2971. [27] C.A. Gandin, J.L. Desbiolles, M. Rappaz and P. Thevoz: Metall. Mater. Trans. A, 1999, 30, 3153. [28] S.G.R. Brown and N.B. Bruce: J. Mater. Sci., 1995, 30, 1144. [29] L. Huo, Z.Q. Han and B.C. Liu: Acta Metall. Sin., 2009, 45, 1414. (in Chinese) [30] L. Huo, Z. Han and B. Liu: Mater. Sci. Forum, 2010, 638–642, 1562. [31] J.W. Cahn and J.E. Hilliard: J. Chem. Phys., 1958, 28, 258. [32] S.M. Allen and J.W. Cahn: Acta Metall., 1979, 27, 1085. [33] B. Bottger, J. Eiken and I. Steinbach: Acta Mater., 2006, 54, 2697. [34] Z.Y. Liu, Q.Y. Xu and B.C. Liu: Int. J. Cast Metals Res., 2007, 20, 109. [35] N. Ricketts: Properties of Cast Magnesium Alloys, http://members.tripod.com/Mg/asm prop.htm. [36] H. Yan, W. Zhuang, Y. Hu, Q. Zhang and H. Jin: J. Mater. Process. Technol., 2007, 187–188, 349. [37] H. Hu and A. Yu: Model. Simul. Mater. Sci. Eng., 2002, 10, 1. [38] Z.A. Fu, Y. Qing and S.M. Xong: Mater. Sci. Forum, 2007, 546–549, 133. [39] M.J.M. Krane, D.R. Johnson and S. Raghavan: Appl. Math. Model., 2009, 33, 2234. [40] G.S. Grest, D.J. Srolovitz and M.P. Anderson: Acta Metall., 1985, 33, 509. [41] A. Jacot and M. Rappaz: Acta Mater., 2002, 50, 1909. [42] Q. Du and A. Jacot: Acta Mater., 2005, 53, 3479. [43] W. Yang, L. Chen and G.L. Messing: Mater. Sci. Eng. A, 1995, 195, 179. [44] H. Yin and S.D. Felicelli: Model. Simul. Mater. Sci. Eng., 2009, 17, 075011. [45] T. Borzsonyi and A. Buka: Phys. Rev. E, 1998, 58, 6236.

[46] L. Nastac: Acta Mater., 1999, 47, 4253. [47] K.R. Elder, M. Grant, N. Provatas and J.M. Kosterlitz: Phys. Rev. E, 2001, 64, 021604. [48] M. Asle Zaeem and S.D. Mesarovic: J. Comp. Phys., 2010, 229, 9135. [49] A.M. Mullis, J. Rosam and P.K. Jimak: J. Cryst. Growth, 2010, 312, 1891. [50] D. Danilov and B. Nestler: Acta Mater., 2006, 54, 4659. [51] M. Asle Zaeem and S.D. Mesarovic: Solid State Phenom., 2009, 150, 29. [52] M. Asle Zaeem and S.D. Mesarovic: ASME Conf. Proc., 2008, 12, 267, doi:10.1115/IMECE2008-66767. [53] M. Asle Zaeem and S.D. Mesarovic: Comput. Mater. Sci., 2011, 50, 1030. [54] COMSOL Multiphysics Users’s Guide, COMSOL Inc., 2010. [55] M.C. Flemings: Solidification Processing, McGrawHill, NewYork, 1974, 148–154. [56] J.D. Hunt and S.Z. Lu: Metall. Maters. Trans. A, 1996, 27, 611. [57] C.H. Caceres, C.J. Davidson, J.R. Griffiths and C.L. Newton: Mater. Sci. Eng. A, 2002, 325, 344. [58] C. Labrecque, R. Angers, R. Tremblay and D. Dube: Can. Metall. Q., 1997, 36, 169. [59] D. Dube, A. Couture, Y. Carbonneau, M. Fiset, R. Angers and R. Tremblay: Int. J. Cast Metals Res., 1998, 11(3), 139. [60] W.P. Sequeira, M.T. Murray, G.L. Dunlop and D.H. StJohn: in Proc. TMS Symposium on Automotive Alloys, Feburary 9–13, Orlando, FA, The Minerals Metals and Materials Society (TMS), Warrendale, PA, 1997, 169. [61] J. Lipton, M.E. Glicksman and W. Kurz: Mater. Sci. Eng., 1984, 65, 57. [62] B. Bottger, J. Eiken, M. Ohno, G. Klaus, M. Fehlbier, R. Schmid-Fetzer, I. Steinbach and A. Buhrig-Polaczek: Adv. Eng. Mater., 2006, 8, 241.