Cellular Automaton Modeling of Dendritic Growth Using a Multi-grid Method

A two-dimensional cellular automaton model with a multi-grid method was developed to simulate dendritic growth. In the present model, we used a triple-grid system for temperature, solute concentration and solid fraction fields as a new approach of the multi-grid method. In order to evaluate the validity of the present model, we carried out simulations of single dendritic growth, secondary dendrite arm growth, multi-columnar dendritic growth and multi-equiaxed dendritic growth. From the results of the grid dependency from the simulation of single dendritic growth, we confirmed that the larger grid can be used in the simulation and that the computational time can be reduced dramatically. In the simulation of secondary dendrite arm growth, the results from the present model were in good agreement with the experimental data and the simulated results from a phase-field model. Thus, the present model can quantitatively simulate dendritic growth. From the simulated results of multi-columnar and multi-equiaxed dendrites, we confirmed that the present model can perform simulations under practical solidification conditions.


Introduction
Recently, with the advances in numerical simulation techniques and computers, numerical models have been developed for the simulation of dendritic growth during solidification. Models based on a cellular automaton (CA) method can quantitatively reproduce most solidification microstructures observed experimentally within an acceptable computational time. CA models have been developed for models of grain growth [1][2][3] and of dendritic growth [4][5][6][7][8]. However, in the model of dendritic growth, a smaller grid size needs to be used, which is usually less than 1 µm, because the smaller radius of curvature at the dendrite tip has to be simulated. Therefore, it is difficult to perform large scale simulations of dendritic growth due to high calculation cost. In order to overcome this difficulty, methods using adaptive mesh [9] and multi-grid methods [6] have been reported. The adaptive mesh technique is effective for reducing the computational time, but the programing is not simple. On the other hand, the programing for the multi-grid method is comparatively simple. Thus, in the present work, we adopted the multi-grid method as an approach for reducing the computational time and carried out two-dimensional CA simulations of dendritic growth. In the present model, different grid sizes were used for the calculations of temperature, solute concentration, and solid fraction fields; in other word, this is a triple-layered multi-grid method. In order to evaluate the validity of the present model, we performed simulations for comparison with the Lipton-Glicksman-Kurz (LGK) model [10], and experimental secondary dendrite arm spacing (SDAS) [11], and phase-field (PF) simulation [12].
where T is temperature, t is time, ΔH is latent heat, C p is specific heat and f S is the solid fraction in each grid.

Concentration field
The governing equation for solute diffusion in solid and liquid phases is given as follows: where c i is the solute concentration and D i is the diffusion coefficient, and the subscript i indicates solid or liquid, respectively. For the solid/liquid interface, the solute concentration c i , and the diffusion coefficient, D i are defined as c i = c L and D i = D L between interface and liquid, and D i = D S between interface and solid. And the solute distribution between solid and liquid is given by the following relationship: where k is the equilibrium partition coefficient, and c S and c L are the solute concentrations of solid and liquid phases, respectively.

Solid fraction
The driving force for dendritic growth is directly related to the difference between the equilibrium solute concentration of liquid, c L * , and the actual solute concentration of liquid, c L , at the solid/liquid interface. In the local equilibrium condition, the equilibrium solute concentration of liquid is given by: where c 0 is the initial solute concentration, T * is the equilibrium temperature at the interface, T L is the equilibrium liquidus temperature, m L is the slope of liquidus, Γ is the Gibbs-Thomson coefficient, K is the interface curvature, ε is the anisotropy parameter, θ is the growth angle between the normal to the interface and the x-axis and θ 0 is the angle between the preferential growth direction and the x-axis. The interface curvature is calculated by the following equation: The equilibrium solute concentration of liquid is compared with the actual solute concentration of liquid that is calculated by Eq. (2). The increase (or decrease) in solid fraction of an interface cell can be given as follows: When the sum of the solid fraction in an interface cell becomes one, that is f S = 1, this interface cell changes its state to solid. When f S = 0, this interface cell changes its state to liquid. The capture rule of the interface cell used in the present model was based on the CA approach. When an interface cell has been fully solidified (melted), its neighboring liquid (solid) cells are captured as interface cells.

Multi-grid method
Generally, the grid size of the used solid fraction field has been the same as that of the solute concentration field. Thus, the multi-grid method of dual-layer for temperature and concentration fields has been adopted. As a new approach, we used a different size grid between the solute concentration and the solid fraction fields in the present model; in other word, this is a triple-grid type CA model. Figure 1 shows a schematic illustration of the multi-grid system used in the present model. The temperature field is divided into the largest grid (Layer-1), the solute concentration field is divided into the second largest grid (Layer-2), and the solid fraction field is divided into the smallest grid (Layer-3). The grid sizes of Layer-1, Layer-2, and Layer-3 are defined as Δx L1 , Δx L2 , and Δx L3 , respectively. The relationships among Δx L1 , Δx L2 , and Δx L3 are Δx L1 = n 12 Δx L2 and Δx L2 = n 23 Δx L3 , that is, number of cells in Layer-2 in Layer-1 is (n 12 ) 2 and number of cells in Layer-3 in Layer-2 is (n 23 ) 2 . The solid fraction field is directly related to the morphology of the solid/liquid interface because it is used in the calculation of interface curvature. Thus, the small grid size is preferred. On the other hand, the grid size of temperature and solute concentration fields is expected to be large because a larger time step is required for the simulations of the dendritic growth under the practical solidification conditions. Using the triple-grid approach, the smaller radius of curvature at the dendrite tip could be accurately simulated by calculating the grid size of Layer-3, and the larger time step obtained from the grid size of Layer-1 and Layer-2 can be used in the simulation. In order to reduce error, the temperature value of each grid of Layer-2 divided into a grid of Layer-1 is given by a linear approximation of the neighboring grids of Layer-1. The relationship between Layer-2 and Layer-3 is also the same.

Grid dependency and comparison with LGK model
In order to evaluate the grid dependency of the present model, we carried out the simulations for a single dendrite of Fe-0.3wt.%C alloy growing at constant melt undercooling, ΔT, with various grid sizes. Table 1 shows the material properties of the Fe-C alloy used in the simulations. The calculation domain was 0.5 x 1 mm, and the initial solid was set on the center of the left wall of the calculation domain. In this calculation, the released latent heat was not considered.   Figure 2 shows the change in the tip growth velocity for various grid sizes of Layer-2 with values of n 23 ranging from one to five at ΔT = 2 K. In the case of n 23 = 1 (Δx L2 = Δx L3 ), the tip growth velocity rapidly decreases when Δx L2 becomes larger than 1.0 µm. On the other hand, for n 23 = 2 and n 23 = 3 (Δx L2 > Δx L3 ), the tip growth velocity nearly maintains the converged velocity (52 µm/s) until the value of Δx L2 reached approximately 2.0 µm, and gradually decreases after that. For n 23 = 5, the tip growth velocity obtained for all of Δx L2 is significantly smaller than the converged velocity. Figure 3 shows the simulated dendrite morphologies for the various Δx L2 and n 23 at a growth time of 2.0 s. In case of n 23 = 1, although the necking of the root of the primary arm was observed when Δx L2 was small, it was not observed as Δx L2 increased. This difference of growth morphology is related to the change in the tip growth velocity. In case of n 23 = 2, 3, and 4, necking was observed except for the case where Δx L2 = 4.0 µm and n 23 = 4. The dendrite morphologies were in good agreement with the dendrite morphology for the smallest Δx L2 . In case of n 23 = 5, the size of dendrites was smaller than that for the other values of n 23 .
Next, we compared the present model and the LGK model. In the LGK model, the thermal diffusion due to the latent heat release has been considered. But, in the present model, there is no thermal diffusion term in Eq. (1). In the condition of a lower melt undercooling, the melt undercooling far away from the dendrite tip is kept almost constant [6]; the thermal gradient is very small. Thus, we assumed that the thermal diffusion could be negligible. Figure 4 shows the comparison of tip growth velocity and tip radius of curvature predicted by the LGK model and the present model for various amounts of melt undercooling. The simulated tip growth velocities were lower than those predicted by the LGK model. These changes in tip growth velocity over undercoolings predicted by the two models show the same trend. The agreement of the simulated tip radius of curvature appears to be very good. Therefore, we confirmed that the simulation of dendritic growth can be quantitatively performed by the present model. It is especially effective to choose the values of parameter n 23 from 2 to 4.    In order to evaluate the reduction of computational cost, we examined the computational time required to perform the simulations for n 23 = 1, 2, and 3 under a condition of Δx L3 = 0.5 µm. Because the calculation domain of 0.5 x 1.0 mm is used for the whole of conditions, in the case of n 23 = 1, the number of total grids is the largest and the time step is the smallest. As a result, the computational time of the simulation for n 23 = 2 and 3 were 1/11 and 1/50, respectively, of that for n 23 = 1. Therefore, the triple-grid method developed in the present model is very effective in reducing the computational cost.

Comparison with experimental SDAS and PF simulation
For the quantitative evaluation of the present model, we predicted the SDAS for Fe-0.1wt.%C, Fe-0.3wt.%C and Fe-0.4wt.%C alloys using the present model, and compared our results with the experimental data by Okamoto et al. [11] and with the PF simulation by Ode et al. [12]. The calculation domain was 2.01 x 0.255 mm, and grid sizes of Δx L1 = 15 µm, Δx L2 = 1.5 µm, and Δx L3 = 0.5 µm were used. Initial planar solid was set on the bottom of the calculation domain, which corresponds to the primary dendrite arm. In this simulation, the released latent heat was also considered. Constant cooling rates from 2 K/s to 50 K/s were used. The values of the SDAS were estimated at the end of solidification before the peritectic temperature (1766 K). Figure 5(a) shows the relationship between the SDAS and the cooling rate in the simulations. For comparison, the experimental data of Fe-0.41wt.%C-6.71wt.%Cr and Fe-0.24wt.%C-14.2wt.%Cr alloys and the results of PF simulation of Fe-0.1wt.%C alloy are also shown in Fig. 5(a). We compared our results with the experimental data for Fe-C-Cr ternary alloys because it was difficult to obtain the experimental data for Fe-C binary alloys. First, we compared it with the experimental data of Fe-0.41wt.%C-6.71 wt.%Cr alloy; the simulated results for Fe-0.4wt.%C alloy were in good agreement with the experimental data. The carbon content of these alloys is very similar. Next, we compared it with the experimental data of Fe-0.24wt.%C-14.2 wt.%Cr alloy; the experimental values of the SDAS were smaller than those of every simulation. This discrepancy is presumably due to the fact that the Cr content added to an alloy has affected the SDAS of the experimental alloy, and this result may be expected that the present model can quantitatively predict the SDAS for a multi-component alloy. We will investigate the validity of a model for multi-component alloys with the present multi-grid approach in near future. In the comparison with the PF simulation for the same carbon content, that is, Fe-0.1wt.%C, the simulated results obtained from the present model were in fairly good agreement with those predicted by the PF model. Figure 5(b) shows the relationship between the SDAS and the local solidification time obtained from the simulations of the present model and the PF simulation for Fe-C binary alloys. The simulated results obtained from both models show a good agreement. From the simulated results obtained from the present model, we estimated the exponent of the relationship between the SDAS and the local solidification time. The value of the exponent was 0.45, which corresponds to the slope of line shown in Fig. 5(b). This value of 0.45 is larger than the value of 1/3 predicted analytically [13]. However, it has been reported that the experimental values of the exponent of some Fe-base alloys or steels can be larger than 1/3 [14,15].

Large-scale simulations of columnar and equiaxed multi-dendritic growth
The present model was used to simulate the formation of both multi-columnar and multi-equiaxed dendrites in Fe-0.15wt%C and Fe-0.3wt%C alloys under practical solidification conditions. The simulations of multi-columnar dendrites were performed within a calculation domain of 3.0 x 2.01 mm and the grid sizes were Δx L1 = 15 µm, Δx L2 = 1.5 µm, and Δx L3 = 0.5 µm. The five initial solids were placed on the left wall of the calculation domain at equal intervals. The temperature gradient, G, was 5 K/mm and the cooling rates were 1 K/s and 2 K/s. The temperature of a grid of Layer-1 reached the peritectic temperature and was maintained at that temperature because the peritectic reaction is not considered in the present model. The simulations were stopped when the temperature of the whole domain reached the peritectic temperature. Figure 6 shows the simulated solute concentration distribution of columnar dendrites for Fe -0.15wt.%C and Fe -0.3wt.%C alloys. Comparing the simulated dendrites growing at cooling rate of 1 K/s (Fig. 6(a1) to (a3)) with those at 2 K/s (Fig. 6(b1) to (b3)) for Fe-0.15wt.%C alloy, the growth of a large number of secondary arms can be observed in the columnar dendrites growing at a higher cooling rate of 2K/s. Comparing the growth of columnar dendrites of Fe-0.15wt.%C alloy with those of Fe-0.3wt.%C alloy at a cooling rate of 1 K/s (ex. Fig.  6(a3) and (c3)), it is obvious that the simulated dendrite morphologies reflect the influence of the initial composition.  The simulations of multi-equiaxed dendrites were performed in a calculation domain of 2.01 x 2.01 mm and the grid sizes were Δx L1 = 15 µm, Δx L2 = 1.5 µm, and Δx L3 = 0.5 µm. The fifty initial solids having random preferred growth orientations were randomly and simultaneously distributed in the calculation domain. The preferred growth orientations were given in the range from -45 to 45 degrees with respect to the horizontal direction. The domain that was assumed to have uniform temperature was initialized at the liquidus temperature and alloy composition, and was cooled at a constant cooling rate of 10 K/s. The simulations were stopped when the temperature of the entire domain reached the peritectic temperature. Figure 7 shows the simulated solute concentration distribution of equiaxed dendrites for Fe-0.15wt.%C and Fe-0.3wt.%C alloys. In the early stage of solidification, the primary dendrite arms developed along their crystallographic orientations. As solidification proceeds, the growth and coarsening of secondary dendrite arms occur. The influence of the initial composition is seen as the simulation of multi-columnar dendrites. These simulations were carried out to investigate the computational performance of the present model; the simulations of multi-columnar dendrites and multi-equiaxed dendrites required a computational time of about one week and two days on a PC Core i7 with 1-CPU 3.4GHz of 6 cores using Open MP parallelization. However, if an approach like a multi-grid method is not used, huge amounts of computational time would be required to perform these large-scale simulations under practical solidification conditions.

Conclusions
A 2-D CA model with a multi-grid method was developed to simulate a dendritic growth. In the present model, we used a triple-grid for temperature, solute concentration and solid fraction fields as a new approach of the multi-grid method. In order to evaluate the present model, we carried out simulations of single dendritic growth, secondary arm growth, multi-columnar dendritic growth and multi-equiaxed dendritic growth. From the simulation results of the present model using the triple-grid, we confirmed that the present model can perform a quantitative simulation of dendritic growth at dramatically reduced computational cost. These simulated results have shown that a practical simulation of dendritic growth in a short computational time can be performed by the introduction of the multi-grid method using triple-grid.