## Services on Demand

## Journal

## Article

## Indicators

- Cited by SciELO
- Access statistics

## Related links

- Cited by Google
- Similars in SciELO
- Similars in Google

## Share

## Materials Research

##
*Print version* ISSN 1516-1439

### Mat. Res. vol.13 no.4 São Carlos Oct./Dec. 2010

#### http://dx.doi.org/10.1590/S1516-14392010000400010

**REGULAR ARTICLES**

**Micro-macroscopic coupling in the cellular automaton model of solidification**

**Vinicius Bertolazzi Biscuola; Marcelo Aquino Martorano ^{*}**

Department of Metallurgical and Materials Engineering, University of São Paulo - USP, Av. Prof. Mello Moraes, 2463, CEP 05508-900, São Paulo, SP, Brazil

**ABSTRACT**

A cellular automaton (CA) model to predict the formation of grain macrostructure during solidification has been implemented and the coupling between the microscopic and the macroscopic submodels has been investigated. The microscopic submodel simulates the nucleation and growth of grains, whereas the macroscopic solves the heat conduction equation. The directional solidification of an Al-7 wt. (%) Si alloy was simulated, enabling the calculation of the temperature and solid fraction profiles. The calculated temperature was used to obtain the solid fraction profile by an application of Scheil equation. This solid fraction disagrees with that calculated in the micro-macro coupling of the model, although this coupling is completely based on Scheil equation. Careful examination of the discrepancies shows that it is a result of the undercoolings for nucleation and growth of grains and also of the interpolations of enthalpy change and temperature from the finite volume mesh to the CA cell mesh.

**Keywords: **mathematical modeling, cellular automaton, solidification, coupling

**1. Introduction**

The cellular automaton (CA) technique is an algorithm constructed to simulate the spatial/temporal evolution of a system by applying simple transformation rules to the sites of a lattice^{1}, which is a set of fixed points. The lattice is used to map some important parameter of the system, such as atomic particles, grains of a metallic alloy, etc. Each point of the lattice must have a state, selected among a set of finite states that represent important information about the phenomena under study. In the case of modeling the grain structure of a given material, the points of the lattice might have a number that indicates the orientation of the grains. Therefore, lattice points within the same grain would have the same number and, an approximate image of all the grain structure could be constructed from the states and positions of all lattice points.

In the CA algorithm, provided the state values of all lattice points are known at the beginning of the simulation, the time evolution of these values is calculated after each time step using transformation rules based on the state of a lattice and its neighbors and, sometimes, on parameters such as temperature and concentration, depending on the type of model. The CA technique has been used in different areas of materials science. Some important applications are the mathematical modeling of the microstructure evolution during recrystallization^{2}, sintering^{3}, and solidification^{4}. Further applications are given by Raabe^{1}. One important result of these mathematical models based on the CA is an image of the microstructure during its evolution and, obviously, the final microstructure.

Gandin and Rappaz^{4} proposed one of the first models of solidification based on the CA technique. This model evolved into the popular CAFE (cellular automaton - finite element) model^{5}, which was applied to predict grain macrostructures of alloys in the as-cast state^{6}. Although modified CA models have been proposed to simulate detailed microstructures^{7}, the CAFE model is still frequently used, because it is more computationally efficient when only the grain structure (without the details of dendrite arms) is required.

In most solidification models, including CAFE, the calculation of the solid fraction field depends on the coupling between microscopic and macroscopic submodels. The time evolution of the solid fraction field during solidification is of paramount importance to: the latent heat release, the solute rejection by the solid to the liquid, and the fluid flow in the mushy zone. Nevertheless, in the cellular automaton model the calculation of this field has never been examined, because the main interest has always been to obtain an image of the grain macrostructure. The objective of this work is to examine the micro-macroscopic coupling and the resulting solid fraction field calculated with the cellular automaton model.

**2. Cellular Automaton Model**

A stochastic model based on CAFE^{5} was developed in the present work by implementing and coupling a microscopic and a macroscopic submodel. A brief description of the model is presented here, emphasizing the important features of the coupling between the submodels (micro-macro coupling), summarized in Figure 1.

The macro-model (macroscopic submodel) consists of the numerical solution by the explicit finite volume method^{8} of the following one-dimensional energy conservation equation

where *H* is the volumetric enthalpy; *T* is the temperature; *t* is time; *y* is the spatial coordinate in the direction of heat flow; and κ is the thermal conductivity, defined as an average of the conductivities of the liquid (κ_{l}) and solid (κ_{s}) phases, weighted by their corresponding local volume fractions. Equation 1 was discretized by a one-dimensional mesh of equal rectangular finite volumes with a node at the centre of each volume and aligned in the *y* direction. From the discretized equation, a volumetric enthalpy change Δ*H*^{t + }^{δ}^{t} was calculated for each finite volume after a time step δ*t*.

The micro-model (microscopic submodel) was developed mainly to simulate the nucleation and growth of dendritic grains in two dimensions during solidification. A two-dimensional microscopic mesh of CA square cells with a site at each cell centre was adopted. Each CA cell had two possible states, namely, liquid and solid. At the beginning of simulation, all cells were liquid and particles representing substrates for heterogeneous nucleation of grains were randomly located in some of these cells. Nucleation of a new grain occurred at a particle when the liquid cell reached a predetermined local undercooling Δ*T*_{N} = *T*_{L} - *T*_{N}, where *T*_{L} and *T*_{N} are the liquidus and nucleation temperatures, respectively. In the present work, the same Δ*T*_{N} value was adopted for all particles, which is equivalent to assume instantaneous nucleation of grains^{9}. At nucleation, the cell state is changed from liquid to solid and a square with a random orientation and centered at the cell site is created to represent a new grain envelope.

The growth of a newly created grain is simulated in the micro-model by increasing the diagonals of the square at a velocity *V* = *A* Δ*T*^{m}, where Δ*T* is the undercooling at the CA site in relation to *T*_{L}, and *A* and *m* are constants dependent on the alloy (Figure 2a). After sufficient growth, some part of the square eventually reached the position of one of the four nearest neighbor sites of the two-dimensional mesh (Figure 2b). Then, the state of the corresponding neighbor cell was changed to solid and a new square was created with an orientation and size defined as in the CAFE model^{5} to represent another part of the same growing envelope (Figure 2c). This whole process is shown schematically in Figure 2.

The strong coupling between the micro-model and macro-model is established in the solid fraction calculation. Immediately after nucleation at a liquid CA cell of the micro-model, its solid fraction ε_{s} is calculated by Scheil equation^{10}, given below:

where *T*_{f} is the melting point of the pure metal and *k* is the solute partition coefficient. The temperature *T* is obtained by interpolation, at the position of the cell site, of the *T* values at the four nearest nodes of the finite volume mesh (macro-model). During grain growth, after each time step (δ*t*) the change in the solid fraction (δε_{s}^{t}^{ + }^{δ}^{t}) is calculated at every CA cell by an equation derived from a combination of Scheil equation and the change in volumetric enthalpy at the cell site position, δ*H*^{t}^{ + }^{δ}^{t}, as follows^{4}

where ρ is the alloy density; *c*_{p} is the specific heat; and *L*_{f} is the latent heat of fusion. The value of δ*H*^{t}^{ + }^{δ}^{t} at the cell site position is also obtained by a linear interpolation of Δ*H*^{t}^{ + }^{δ}^{t}, calculated at the nodes of the finite volume mesh (macro-model). During solidification, the interpolated cell temperature eventually reaches the eutectic temperature, *T*_{E}. At this moment, the temperature is kept constant at *T*_{E} and the change in solid fraction at the cell is calculated as follows^{4}

At a finite volume of the macroscopic model, the change in solid fraction, De_{s}^{t}^{ + }^{d}^{t}, is obtained as in CAFE^{5}, i.e., by averaging the values of Δε_{s}^{t}^{ + }^{δ}^{t} at all cells within this finite volume, weighted by the same interpolation coefficients used to interpolate *T* and δ*H*^{t + }^{δ}^{t}. Finally, the temperature at a node of the finite volume mesh (macro-model) is calculated by the following equation, completing the micro-macro coupling cycle (Figure 1).

Results obtained with the computer code implemented in the present work showed excellent agreement with those presented by Rappaz and Gandin^{4,11} for the CAFE model.

**3. Analysis of the Micro-macro Coupling**

The unidirectional solidification of an Al-7 wt. (%) Si alloy was simulated in a rectangular two-dimensional domain of size 0.15 m in the heat flow direction (*y*) and 0.05 m perpendicular to it (*x*). An adiabatic boundary was assumed at the top of the domain (*y* = 0.15 m), while through the bottom (*y* = 0), the heat flux was given by *q* = *h*(*T* - *T*_{∞}), where *h* = 250 W.m^{-2}.K^{-1} is the heat transfer coefficient and *T*_{∞} = 298 K is a reference temperature. A uniform temperature of 991 K was assumed as the initial condition. The number density of substrate particles for the nucleation of grains was defined as *n*_{T,bulk} = 5 × 10^{5 }m^{-3} for the bulk liquid and *n*_{T,layer} = 3 × 10^{6 }m^{-2} for the liquid layer adjacent to the lower boundary (*x* = 0). Since the mesh of CA cells is two-dimensional, these values were converted by^{11} and The alloy properties were: ρ = 2452 kg.m^{-3}; *c*_{p} = 1126 J.kg^{-1}.K^{-1}; κ_{l} = 60.5 W.m^{-1}.K^{-1}; κ_{s} = 137.5 W.m^{-1}.K^{-1}; *L*_{f} = 387400 J.kg^{-1}; *k* = 0.13; *T*_{L} = 891 K; *T*_{E} = 850 K; *T*_{f} _{ }= 933 K; *A* = 3 × 10^{-6} m.s^{-1}.K^{-}^{m}; m = 2.7; and Δ*T*_{N} = 3K. Convergence of the numerical solution was achieved with a finite volume mesh (macro-model) of 30 volumes and a CA mesh (micro-model) of 300 (*y*) × 100 (*x*) cells.

The calculated grain structure and the profiles of ε_{s} and *T* at the finite volume nodes (indicated by dots) are shown in Figure 3. The solid fraction given by Scheil model, ε_{s}^{Scheil}, was obtained with Equation 2, in which *T* is the temperature at the finite volume nodes. It is possible to see that ε_{s}^{Scheil} does not agree with the ε_{s} profile calculated with the CA model equations (Equations 2, 3 and 4), which were based on the Scheil model.

To confirm the discrepancy between ε_{s}^{Scheil} and ε_{s}, the implemented CA code was adjusted to simulate the unidirectional solidification problem studied by Guillemot et al.^{12}, in which the growth of only one grain, positioned at the CA cell adjacent to the lower boundary, is simulated,. Therefore, after nucleation the grain envelope grows throughout the domain from the lower to the upper boundary. The calculated *T* and ε_{s} as a function of time are compared with Guillemot et al.^{12} results in Figure 4. The solid fraction ε_{s}^{Scheil} obtained with Scheil equation (Equation 2) using the calculated temperature at 80 mm from the lower boundary is significantly different from ε_{s} obtained with both the implemented model and Guillemot et al.^{12} model, confirming the type of discrepancy observed in Figure 3.

**4. Effect of Undercooling and Interpolation on the Micro-macro Coupling**

Three types of discrepancy are observed in Figure 3 and Figure 4b: (1) ε_{s}^{Scheil} > ε_{s} = 0 at the undercooled liquid ahead of the solid-liquid interface; (2) ε_{s}^{Scheil} < ε_{s} in the mushy zone, down to *T*_{E}; and (3) ε_{s} < 1 below *T*_{E}, although solidification should end at *T*_{E}. The first discrepancy exists because the Scheil equation (Equation 2), used to obtain ε_{s}^{Scheil}, neglects the undercooling for nucleation and growth of grains^{10}. Therefore, it predicts ε_{s}^{Scheil} > 0 immediately below *T*_{L}. In the stochastic model, however, the nucleation undercooling Δ*T*_{N} (= 3K) is partly responsible for the undercooled liquid ahead of the solid-liquid interface. The other part is the undercooling necessary for grains to grow up at a velocity *V*, approximately equal to the *T*_{L} isotherm velocity.

To confirm these effects, the nucleation undercooling was removed (Δ*T*_{N} = 0) and the growth undercooling was eliminated by inserting a substrate particle within every CA cell of the mesh. In these conditions, the profile of ε_{s} calculated with the CA model shows very good agreement with ε_{s}^{Scheil} near *T*_{L} (Figure 5), but some discrepancy still persists within the mushy zone and below *T*_{E}.

This discrepancy is related to the interpolations of Δ*H*^{t}^{ + }^{d}^{t} and δε_{s}^{t}^{ + }^{δ}^{t} to calculate δ*H*^{t}^{ + }^{δ}^{t} and Δε_{s}^{t + }^{δ}^{t}, respectively (Figure 1). The Scheil model is implicitly given by Equation 3; consequently, δε_{s}^{t}^{ + }^{δ}^{t} calculated at a CA cell is exactly related to δ*H*^{t}^{ + }^{δ}^{t} by the Scheil model. Nevertheless, the node temperature is obtained from Equation 5, using Δ*H*^{t}^{ + }^{δ}^{t} and Δε_{s}^{t}^{ + }^{δ}^{t}, which are not related by Equation 3. Actually, the relation between Δ*H*^{t}^{ + }^{δ}^{t} and Δε_{s}^{t}^{ + }^{δ}^{t} is complex, because Δε_{s}^{t}^{ + }^{δ}^{t} is obtained from the weighted average of δε_{s}^{t}^{ + }^{δ}^{t} values at all cells within the finite volume, which are individually calculated by Scheil equation (Equation 3) using an interpolation of Δ*H*^{t}^{ + }^{δ}^{t} at each cell site. Therefore, Δε_{s}^{t}^{ + }^{δ}^{t}, Δ*H*^{t}^{ + }^{δ}^{t}, and the temperature at the finite volume node are not related by Scheil model, explaining the difference between ε_{s} and ε_{s}^{Scheil}.

To eliminate the interpolation effect described above, the number of finite volumes in the *y* direction was increased and made equal to the number of CA cells in the same direction, resulting in δε_{s}^{t}^{ + }^{δ}^{t} = Δε_{s}^{t}^{ + }^{δ}^{t} and δ*H*^{t}^{ + }^{δ}^{t} = Δ*H*^{t}^{ + }^{δ}^{t}. Now, Δε_{s}^{t}^{ + }^{δ}^{t}, Δ*H*^{t + }^{δ}^{t}, and the temperature at the finite volume node are all linked by the Scheil model. The nucleation and growth undercoolings were also removed as in the previous simulation. Figure 6 shows that the resulting ε_{s} profile closely agrees with that from the Scheil model, ε_{s}^{Scheil}, confirming the interpolation effect.

Examining the equations used to calculate the solid fractions ε_{s}^{Scheil} (Equation 2) and ε_{s} (Equation 3) might give futher insight into the discrepancies discussed in the previous paragraphs. Equation 3 is a differential form of Scheil equation, in which a relationship between δ*T* and δ*H* (based on thermodynamics) was substituted. Equation 2, which is the traditional Scheil equation, can actually be obtained after numerical integration of Equation 3, using as initial condition that solidification begins at the liquidus temperature of the alloy, i.e., ε_{s} = 0 for *T* = *T*_{L}. When this Equation 3 is used in the CA microscopic submodel, however, this initial condition is not used, because when *T* = *T*_{L} at a CA cell, nucleation or grain growtn might not have activated the site of this cell. Therefore, solidification does not necessarily begin locally at the liquidus temperature, a behavior that differs from the hypothesis adopted in the tradicional Scheil equation. Furhtermore, another reason for the discrepancy is that ε_{s}^{Scheil} is calculated at the finite volume node, directly from its temperature, whereas ε_{s} is obtained from an average of solid fraction changes calculated in all cells located within the finite volume, as explained in section 2.

**5. Summary and Conclusions**

A cellular automaton (CA) model to predict the grain macrostructure formation during solidification has been implemented and the micro-macro coupling of this model investigated. In the coupling cycle, after each time step of the numerical method, a change in solid fraction is calculated at the CA cells by a Scheil model. The solid fraction (ε_{s}) in the finite volume is calculated from an averaged change in all cells located within this volume. A second type of solid fraction (ε_{s}^{Scheil}) was also calculated in the finite volume using its temperature; this fraction was compared with that obtained from the micro-macro coupling (ε_{s}). Three types of important discrepancies were observed: (1) ε_{s}^{Scheil} > ε_{s} = 0 at the undercooled liquid ahead of the solid-liquid interface; (2) ε_{s}^{Scheil} < ε_{s} in the mushy zone, down to the eutectic temperature (*T*_{E}); and (3) ε_{s} < 1 below *T*_{E}, although solidification should end at *T*_{E}. These discrepancies were examined and their cause were seen to be the undercoolings for nucleation and growth of grains and the interpolations of enthalpy change and temperature from the finite volume mesh to the CA cell mesh, required in the micro-macroscopic coupling.

**Acknowledgements**

The authors thank FAPESP (grant 03/08576-7) and CNPq (grant 475451/04-0).

**References**

1. Raabe, D. *Computational materials science*: the simulation of materials microstructures and properties. New York: Wiley-VCH, Weinheim; 1998. [ Links ]

2. Hesselbarth, HW and Gobel IR. Simulation of recrystallization by cellular automata. *Acta Metallurgica et Materialia*. 1991; 39(9):2135-2143. [ Links ]

3. Pimienta PJP, Garboczi EJ and Carter WC. Cellular automaton algorithm for surface mass transport due to curvature gradients. *Computational Materials Science*. 1992; 1(1):63-77. [ Links ]

4. Gandin CA and Rappaz M. A coupled finite-element cellular-automaton model for the prediction of dendritic grain structures in solidification processes. *Acta Metallurgica et Materialia*. 1994; 42(7):2233-2246. [ Links ]

5. Gandin CA and Rappaz M. A 3D cellular automaton algorithm for the prediction of dendritic grain growth. *Acta Materialia*. 1997; 45(5):2187-2195. [ Links ]

6. Rappaz M, Gandin CA, Desbiolles JL. and Thevoz P. Prediction of grain structures in various solidification processes. *Metallurgical and Materials Transactions a-Physical Metallurgy and Materials Science*. 1996; 27(3):695-705. [ Links ]

7. Zhu MF and Hong CP. A modified cellular automaton model for the simulation of dendritic growth in solidification of alloys. *Isij International*. 2001; 41(5):436-445. [ Links ]

8. Patankar SV. *Numerical heat transfer and fluid flow*. New York: Hemisphere Pub. Corp.; 1980. [ Links ]

9. Stefanescu DM. Methodologies for modeling of solidification microstructure and their capabilities. *Isij International*. 1995; 35(6):637-650. [ Links ]

10. Kurz W and Fisher DJ. *Fundamentals of solidification*. Aedermannsdorf: Trans Tech Publications; 1989. [ Links ]

11. Rappaz M and Gandin CA. Probabilistic modeling of microstructure formation in solidification processes. *Acta Metallurgica et Materialia*. 1993; 41(2):345-360. [ Links ]

12. Guillemot G, Gandin CA, Combeau H and Heringer R. A new cellular automaton - finite element coupling scheme for alloy solidification. *Modelling and Simulation in Materials Science and Engineering*. 2004; 12(3):545-556. [ Links ]

Received: July 9, 2010; Revised: September 29, 2010

* e-mail: martoran@usp.br