Materials Science and Engineering A 421 (2006) 208–216

Stochastic simulation of grain growth during continuous casting

A. Ramirez a,∗ , F. Carrillo b , J.L. Gonzalez c , S. Lopez d
a Department of Aerounatical Engineering, S.E.P.I., E.S.I.M.E., IPN, Instituto Politécnico Nacional (Unidad Profesional Ticomán),
Av. Ticomán #600, Col. Ticomán, C.P.07340, México, Mexico
b Department of Processing Materials, CICATA-IPN Unidad Altamira Tamps, México, Mexico
c Department of Metallurgy and Materials Engineering, E.S.I.Q.I.E.-IPN, México, Mexico
d Department of Molecular Engineering of I.M.P., AP 14-805, México, Mexico

Received 7 December 2004; received in revised form 18 January 2006; accepted 19 January 2006

Abstract
The evolution of microstructure is a very important topic in material science engineering because the solidification conditions of steel billets
during continuous casting process affect directly the properties of the final products. In this paper a mathematical model is described in order to
simulate the dendritic growth using data of real casting operations; here a combination of deterministic and stochastic methods was used as a
function of the solidification time of every node in order to create a reconstruction about the morphology of cast structures.
© 2006 Elsevier B.V. All rights reserved.

Keywords: Continuous casting; Dendrite; Mathematical simulation; Chill; Columnar and equiaxed zones

1. Introduction ers ways are combination of these two. In deterministic methods,

the behavior of each grain is determined by established parame-
The final quality and properties of castings products depends ters like size, number of neighboring grains and others, while the
on the microstructure formed during solidification. The simu- stochastic methods were first used by Louat and then modified
lation on grain scale is used to provide a reliable foundation by some other authors [4–7].
for the process control and the improvement of the casting Deterministic approach refers explicitly to the driving force
products so, that many authors have been developing mathe- of grain growth, which is due to the decrease in total grain bound-
matical models, algorithms and simulation systems to predict it ary energy. Nevertheless, this approach is standard but leads to
[1–20]. wrong predictions of the grain size distribution. Therefore it
In the 1960s Oldfield [1] was one of the beginners who tried is necessary to increase the state parameters and obtain better
to simulate the solidification structure. In the 1980s, the Monte results. On the other hand, the stochastic approach is much more
Carlo method was used for modeling the nucleation and growth successful in predicting size distributions because this approach
of the grains. In the early 1990s, cellular automaton model was can reproduce the heterogeneity of materials using random grain
employed in which the physical mechanism of heterogeneous sizes and populations.
nucleation and growth was taken into account and in a similar A very important problem involved in grain growth simula-
way deterministic methods were introduced to deal with the tion is that the number of cells needed to simulate the microstruc-
distribution of the nuclei. In recent years, the phase-field model ture of a casting is tremendous due to the scale of the grain size.
has been developed to simulate the formation of microstructures Consequently, the efficiency of a personal computer will be very
[2,3]. low, the reason why some authors have explored using parallel
In general, there are two basic theoretical approaches to computing techniques [3]. Nevertheless, improving the algo-
describe the grain growth: deterministic and stochastic, the oth- rithms to simulate accurately the phenomenon is other possible
way. It is important to remember that the number of cells is a
function of the scale used to simulate the micro- or macro-scale
∗ Corresponding author. Tel.: +55 57 29 6000x55127. phenomena, and very large array sizes require longer times to
E-mail address: adalop123@mailbanamex.com (A. Ramirez). analyze and create a good approach.

0921-5093/$ – see front matter © 2006 Elsevier B.V. All rights reserved.
doi:10.1016/j.msea.2006.01.077
 A. Ramirez et al. / Materials Science and Engineering A 421 (2006) 208–216 209

Nomenclature
a1 chill zone size (mm)
a2 columnar zone size (mm)
a3 equiaxed zone size (mm)
d1 –d4 distances from the pivoting node to billet surfaces
(mm)
nR random number
t simulation time (s)
ta time defining if a node is nucleated or growth
t1 time defining limit of chill zone (s)
t2 time defining limit of the columnar zone (s)
tI,J pivoting node
tliq,I,J time for which liquid state of the pivoting node
exists (s)
tmushy,I,J time for which mushy state of the pivoting node
exists (s)
tsol,I,J final point of solidification of the pivoting node
(s)
Fig. 1. Sketch of ingot grain structure showing chill, columnar and equiaxed
tx,y neighbor of the pivoting node zones [1].
Xsol,I,J solid fraction in pivoting node (%)
x,y positions of the neighbors of the pivoting node
after nucleation, the dendrite was considered to grow along four
directions (branches at 90◦ ) and finally tend to adopt a square
2. Nucleation and growth model shape. The result is a group of randomly oriented grains that
grow until their borders make contact.
The first step to develop a model to describe the grain growth
process is to understand the phenomena physically involved in L(θ) = L0 [1 + (a − 1) cos 4θ] (1)
it.
When a liquid metal is quenched, the nucleation and growth η = 1 + γ cos 4β (2)
process will appear. The atoms in the domain are moving in
Other authors have been developing models to create one
liquid conditions, but they form clusters as their temperature
in combination with some phenomena like recrystallization
decreases below liquidus. When the cluster grows to a certain
[7,10–12], but a very important point is the necessity for pro-
size it becomes a crystalline nucleus, and this process is called
gramming a micro-scale model without losing the macro-scale
nucleation.
viewpoint. In recent years, many models have been developed
When a nucleation point has been established, one or more
[13–20] to simulate grain growth but no one has joined the cal-
of the nearest neighbors may tend to join with it; this process is
culation of a thermal behavior and the solidification times as a
called grain growth.
base to build a probable structure of the cast products using real
The solidification grain structure of a cross-section of an ingot
operating conditions.
shown in Fig. 1. It was described by Flemings [8]. Here three dif-
ferent zones can be identified, one of them is the outer zone also
called chill zone, it comprises of fine grains with random ori- 3. Mathematical model
entation. The second, an intermediate zone is a columnar zone,
with many elongated and oriented grains from the billet surfaces The mathematical model used solves the equations for heat
to the centre. Finally, there is a central equiaxed zone compris- transfer involved during the continuous casting process; firstly,
ing of less randomly oriented grains. Steel solidification during the steel is discretizated using a regular square grid that repre-
continuous casting is dendritic. Each grain contains one dendrite sents volumes of liquid steel. Heat removal is calculated accord-
with a main arm and many secondary arms. ing to the mechanism of the steel position along the cast machine.
Feng et al. [3] and Lan et al. [9] used a cellular automaton For the billet surface different border conditions were defined
model to describe the microstructure of some alloys based on to calculate heat removal.
a very similar equation for drawing the profile of a dendrite. In the mold, heat removal is calculated using Eq. (3)
To predict its evolution in a micro-scale they used Eqs. (1) and √
(2) to describe the geometry adopted by the nucleated points, q = A 0 + B0 t (3)
where θ and β are the polar angles of the dendrite, a = 1.25
a constant and L and η are the average dendrite radii. These where the values for the coefficients A0 and B0 are 2680 and
authors show the evolution of the dendrite shape in a liquid pool; 335, respectively, as obtained by Savage and Pritchard [21].
210 A. Ramirez et al. / Materials Science and Engineering A 421 (2006) 208–216

as is shown in Eq. (6)

P+1 P
1 TI,J − TI,J 1
= (TI−1,J + TI+1,J − 2TI,J )
αI,J t ( x)2
1
+ (TI,J−1 + TI,J+1 − 2TI,J ) (6)
( y)2
where α is the value of thermal diffusivity
The corresponding thermal behavior of steel is shown in
Fig. 2. In the present study, to test the model, we considered
a square billet of 1044 steel (160 mm × 160 mm) cast under the
cooling conditions described in Table 1. The cast speed was
2.0 m/min. After determining the information of the nodes, the
values for changing of liquid to mushy and solid states can be
classified as a function of their solidification times. In Fig. 3a–c
Fig. 2. Temperature and solidification profiles. are shown the corresponding profiles of residence time on the
liquid, mushy and solid states.
When the steel is traveling below the spray zones, conduction In the phases and cellular automaton models, the grain struc-
force [22] of Eq. (4) is used to solve it, ture depends on the energy loss of every node as described by
many authors [13–20]. The simulation process uses energy loss
q = h(TI,J − TW ) (4)
to define the solidification times of every node of the billet dur-
where TW is the water temperature of every spray zone, TI,J the ing thermal behavior calculation; for this reason the procedure
temperature of the node in the billet surface and h is a coefficient employed in the present work can be considered as an indirect
corresponding to the intensity of the water flow rate. method to simulate a solidification process.
Finally, at the end of the running, steel passes on a free zone
where the equation of Stephan–Boltzman [23] is used
q = σε(TI,J
4
− Tamb
4
) (5)
Inside the billet, heat removal is transported only by conduc-
tion and finite difference method is employed to solve the region

Table 1
Operating conditions of the spray zones
Case 1

Zone 1
Water flow rate (l/min) 175
Spray number 1
θ 2
Ω 60
dbs (mm) 160
Zone 2
Water flow rate (l/min) 225
Spray number 9
θ 12.6
Ω 60
dbs (mm) 160
Zone 3
Water flow rate (l/min) 215
Spray number 9
θ 16
Ω 60
dbs (mm) 160
Zone 4
Water flow rate (l/min) 195
Spray number 9
θ 16
Ω 60
dbs (mm) 160
Fig. 3. (a–c) Profiles of the residence time of billet solidification.
 A. Ramirez et al. / Materials Science and Engineering A 421 (2006) 208–216 211

4. Identification of the dendritic zone

In the present study grain growth was treated as a case of

drunken walk problem. The flow chart for the solidification pro-
cess is shown in Fig. 4, where the time loop t will command the
simulation; when the simulation time is identical to the solidi-
fication time of a node (t = tSol,I,J ), this node is considered as a
solidified one, The time criterions for defining the dimensions
of chill, columnar and equiaxed zones (a1 , a2 and a3 ) and their
morphologies (t1 and t2 ) are shown in Fig. 5.
The corresponding zone of every node is identified according
to the shortest time difference between the pivoting node and its
neighbors (tshort,I,J ) using Eq. (7). For the analysis the nearest

Fig. 5. Distances of the chill columnar & equiaxed zones.

Fig. 6. Pivoting node (tI,J ) and its neighbors.

eight neighbors to the pivoting node are considered as is shown

in Fig. 6. The neighbors’ positions on the array are defined with
the subscript (x, y)
tshort,I,J = tmushy,I,J − tmushy,x,y (7)
Fig. 7a–c shows a classification for the liquid, mushy and
solid time differences. The time criterions t1 and t2 used for
command the process are 2 and 12 s, respectively. Here it can be
seen that, near the billet surfaces, the residence times are very
short. Therefore this area is short too, but inside the billet the
residence time of steel in the mushy zone is longer, and the grain
sizes will be bigger too. It can be appreciated here that the best fit
for describing the distances of the zones can be obtained using
the values of the mushy zone (tmushy,I,J ).

5. Algorithms for dendritic zones

After this, the next step is to create three different algorithms

Fig. 4. Flow chart for describing the solidification process and defining dendrit- to simulate the grain structure of every zone which must be based
ical zones. on the nucleation and grain growth process.
212 A. Ramirez et al. / Materials Science and Engineering A 421 (2006) 208–216

Fig. 8. Logical flow chart to select a nucleation or growth mechanism.

The flow chart of Fig. 9 shows the procedure to find the

shortest time difference between the pivoting nodes and their
neighbors, here the solidification time of the neighbors are tx,y ,
while tB is the time difference between the pivoting node and its
neighbors, tA is the shortest time and nn is the number of neigh-
bors with less time difference. This procedure is very important
for the next zones to establish a first approach about the direction
of the dendrites as a function of the heat flow direction.
In the columnar zone, dendrites grow from the end of the chill
zone to the billet centre; and their orientation is perpendicular
to the billet faces. For this reason, it is necessary to improve
the coalescence and orientation mechanisms of dendrites. This
is done by calculating the distances from the pivoting node to
the billet surfaces (d1 , d2 , d3 , and d4 ) shown in Fig. 10. After
this, the shortest distance (d3 for the example of the figure) will
Fig. 7. (a–c) Classification of pivoting nodes as a function of their difference receive a preferential probability to grow with the pivoting node.
with their neighbors.
The logical process of this operation is shown in Fig. 11. It is
important to mention that in the columnar zone the probability
of growth is bigger in contrast with the chill zone due to the
When the billet is quenched quickly in the cast operations, longer residence time in the mushy zone.
chill is the first zone to appear. In this zone, there are many nucle- The equiaxed zone has many randomly oriented grains but
ation points, because the steel remains in the mushy zone for a they are bigger than those in the chill zone. This zone is devel-
very short time. The change from liquid to solid state is almost oped when the chill and columnar zones will enclose a liquid
instant, and the differences between the pivoting nodes and their pool in the billet centre and the heat removal is less intense;
neighbors are very short too; for this reason many nucleation this phenomena was considered to be possible to appear only
points will appear in this zone and every node will tend to adopt if the nodes remain for a very long time in the mushy zone.
the nucleation process. The structure obtained here will have Then a combination of mechanisms is necessary to simulate
many fine grains randomly oriented; this kind of structure is grain growth. Here nucleation and growth mechanism occurs
found at very high solidification speeds. Chill zone is formed simultaneously due to the very slow speed of solidification. The
when the billet is in the mold and it must be thick enough to algorithm to simulate it is described below.
contain the pressure of the liquid steel inside.
The process to select a solidification mechanism is shown in 5.1. Nucleation process
Fig. 8, where ta is the time criterion to determine if the pivoting
nodes will be nuclei or growth. This process is employed in the In this algorithm, when the solid fraction in a node Xsol,I,J
columnar and equiaxed zones but the value of the criterion is exceeds 50%, the probability to form a nucleation point is
different. assigned at every iteration and it is compared with a random
 A. Ramirez et al. / Materials Science and Engineering A 421 (2006) 208–216 213

Fig. 11. Logical flow chart to adopt a preferential probability for a neighbor.

number of a stochastic process (Monte Carlo method). The result

is the appearance of many random nucleation points at every
time step in the simulation. The solid fraction in each node is
obtained for every time step using Eq. (8):
t t
Xsol,I,J = = (8)
(tsol,I,J − tliq,I,J ) tmushy,I,J

Fig. 9. Flow chart for finding the shortest time difference.

Fig. 12. Calculation of the solid and liquid fractions.

Fig. 10. Distances from pivoting node to the billet surfaces. Fig. 13. Flow chart to know if one of the mechanisms occurs.
214 A. Ramirez et al. / Materials Science and Engineering A 421 (2006) 208–216

Fig. 14. Evolution of the solidification profile (time in seconds).

A. Ramirez et al. / Materials Science and Engineering A 421 (2006) 208–216 215

Fig. 15. Evolution of the grain structure in the billet (time in seconds).
216 A. Ramirez et al. / Materials Science and Engineering A 421 (2006) 208–216

where t is the simulation time and tsol,I,J and tliq,I,J are the time a function of the time during the simulation process. Here it is
for changing the state of every node; for example in Fig. 12, possible to see that during the initial seconds of the solidification
the points (a–c) have the following fractions of solid and liquid process the chill zone is formed around the billet surface, because
concentrations: the solidification speed is very high. After this, the columnar
zone with oriented dendrites begins to grow, and finally after
a is 75% liquid and 25% solid; 30 s many stationary nuclei points begin to appear and begin to
b is 50% liquid and 50% solid; grow to create the equiaxed zone.
c is 25% liquid and 75% solid.
6. Conclusions
When the simulation time (t) is equal to the time at point ‘a’,
the node is not considered for the solidification process, when t A first comparison of the simulation profile obtained using the
is equal to the time in point ‘b’ the node begins to be consider algorithm with a real macro-etching ingot shows a very similar
for solidification process and this node can probably become representation of the continuous casting structure.
a nuclei. When t is equal to the time in point ‘c’, the node is
involved directly in the solidification algorithm and it can be Acknowledgments
considered as a nuclei or a probable point for coarsening with
any other neighbor. The authors wish to thank Consejo Nacional de Cien-
cia y Tecnologı́a (CONACyT), Instituto Politécnico Nacional
5.2. Growth process (IPN-ESIME-ESIQIE) and Instituto Mexicano del Petróleo
(IMP).
When the time of simulation and the solidification time of
the pivoting node are the same, growth occurs. Here the process References
to find the shortest time difference of a neighbor is repeated as
in the columnar zone; in addition, the first solidified neighbor [1] W. Oldfield, ASM Trans. 59 (1966) 945.
must be found; but it is not necessary to assign a preferential [2] N. Provatas, N. Golden feld, A. Danzing, in: B.G. Thomas (Ed.), Mod-
possibility to one of them, so that all the solidified neighbors eling of Casting Welding & Advanced Solidification Process, vol. VIII,
will have the same chance to coalesce with the pivoting node; as TMS, Warrendale, PA, 1998, p. 533.
[3] W. Feng, Q. Xu, B. Liu, ISIJ 42 (7) (2002) 702–707.
a result of this algorithm, the morphology, orientation and size [4] N.P. Luoat, Acta Metall. 22 (1974) 721.
of the dendrites will be missed (anisotropic). [5] N.P. Luoat, Mater. Sci. Forum 94–96 (1992) 67–76.
The process described in Fig. 8 is used again in the equiaxed [6] O.M. Ivasishin, S.V. Shevchenko, N.L. Vasiliev, Acta Mater. 51 (2003)
zone; if the pivoting node (tI,J ) has been solidified, this node is 1019–1034.
[7] A.M. Gusak, K.N. Tu, Acta Mater. 51 (2003) 3895–3904.
ignored; but if the node remains in the mushy region, one of the
[8] M.C. Flemings, Solidification Processing, New York Ed., McGraw Hill
mechanism can be present. The flow chart is shown in Fig. 13; Book Co., 1974.
here a random number, nR , is generated as a function of tS2,I,J ; [9] C.W. Lan, Y.C. Chang, C.J. Shih, Acta Mater. 51 (2003) 1857–1869.
then nR is compared with the parameter tS3,I,J to take a decision. [10] E.A. Holm, M.A. Miodownik, A.D. Rollett, Acta Mater. 51 (2003)
tS1,I,J is the time when solidification of the node pivoting has 2701–2716.
[11] S. Mishra, T. DebRoy, Acta Mater. 52 (2004) 1183–1192.
begun and it is equal to 50% of the solid fraction. It is calculated
[12] C.E. Krill III, L.Q. Chen, Acta Mater. 50 (2002) 3057–3073.
using Eq. (9), while the parameters tS2,I,J and tS3,I,J are calculated [13] A. Carpinteri, B. Chiaia, P. cornetti, Mater. Sci. Eng. A 365 (2004)
using Eqs. (10) and (11): 235–240.

 [14] H. Yoshioka, Y. Tada, Y. Hayashi, Acta Mater. 52 (2004) 1515–1523.
tsol,I,J − tliq,I,J
tS1,I,J = tliq,I,J + (9) [15] Y.J. Lan, D.Z. li, Y.Y. Li, Acta Mater. 52 (2004) 1721–1729.
2 [16] K. Yee, C.P. Hong, ISIJ 37 (1) (1997) 38–46.
[17] Y.H. Shin, C.P. Hong, ISIJ 42 (4) (2002) 359–367.
tS2,I,J = tsol,I,J − tS1,I,J (10) [18] M.F. Zhu, C.P. Hong, ISIJ 41 (5) (2001) 436–445.
[19] M.F. Zhu, J.M. Kim, C.P. Hong, ISIJ 41 (9) (2001) 992–998.
tS3,I,J = t − tS1,I,J (11)
[20] M. Qan, Z.X. Guo, Mater. Sci. Eng. A 365 (2004).
In order to understand the meaning and working of the algo- [21] J. Savaje, W.H. Pritchard, Iron Steel Inst., vol. 179, pp. 269–277.
[22] Heat Transfer Problem Solvers, Research & Education Association, Pis-
rithms described here, the solidification process is shown in cataway, New Jersey, 1993. pp. 17–81, 385–410, 424–461, 505–520.
Fig. 14, the liquid, mushy and solid regions are shown. Fig. 15 [23] C.E. Wicks, R.E. Wilson, J.R. Welty, Fundamentals of Momentum, Heat
shows some sketches of the evolution of the grain growth as and Mass Transfer, John Wiley & Sons, USA, 1984, pp. 269–380.