APCOM’07 in conjunction with EPMESC XI, December 3-6, 2007, Kyoto, JAPAN

Static Recrystallization Phase-Field Simulation Based on The
Predicted Subgrain Structures
T. Takaki1*, A. Yamanaka2, Y. Tomita3
1

Graduate School of Science and Technology, Kyoto Institute of Technology, Matsugasaki, Sakyo, Kyoto,606-8585, Japan
2
Graduate School of Science and Technology, Kobe University, 1-1, Rokkodai, Nada, Kobe, 657-8501,
Japan
3
Graduate School of Engineering, Kobe University, 1-1, Rokkodai, Nada, Kobe, 657-8501, Japane-mail: takaki@mech.kit.ac.jp, yamanaka@solid.mech.kobe-u.ac.jp, tomita@mech.kobe-u.ac.jp
Abstract In this study, a coupled model that estimates recrystallization kinetics and structure is developedby employing a phase-field method and crystal plasticity theory. In this model, the deformation subgrain
structure which is the initial structure of recrystallization simulation is predicted from the stored energy
and crystal orientation calculated from crystal plasticity finite element simulation. Then, based on thepredicted subgrain structures, static recrystallization are simulated by multi-phase-field method taking
account of the orientation dependence of boundary energy and mobility. The model developed in thisstudy enable to consider the deformation structures and the spontaneous nucleation through the abnormal
subgrain growth, unlike the conventional model driven by the stored energy and using the nucleation
criteria.
Key words: Static Recrystallization, Multi-Phase-Field Method, Subgrain Structure, Abnormal Subgrain
Growth, Crystal Plasticity, Coupled Simulation
INTRODUCTION
The development of recrystallization model is essential to predict the recrystallization texture,
microstructure and kinetics during annealing. There are some computational models for recrystallization[1], such as vertex or network model, Monte Carlo Potts model [2], cellular automata model [3] andphase-field model [4, 5]. In all these models, the introduction of the realistic deformation microstructure
which is the initial structure of recrystallization simulation is indispensable for the accurate estimations ofrecrystallization process. Therefore, the deformation microstructures are obtained experimentally [6] ornumerically [7].
There are two types of recrystallization model: One is the grain growth model driven by the stored energy[4, 6, 7], in which the deformation microstructure is characterized by the crystal orientation and stored
energy. Since this model needs nucleation criteria, the calculated recrytallization kinetics and structurechange depending on the predefined nucleation conditions. The other model describes grain growth by thebalance of grain boundary energies, and considers the subgrain structure in which there are no dislocations
inside individual subgrain [5]. Although this model requires smaller computational grids and muchcomputational cost comparing to the first model in order to express subgrain microstructures, byintroducing the misorientation dependence of boundary energies and mobility to the model, the natural
nucleation of recrystallized grains is enabled through the abnormal grain growth [8].
In this study, we develop a coupled model for static recrystallization of a phase-field method that describes
second model mentioned above and crystal plasticity finite element method. Since the normal crystal
plasticity finite element simulation never give subgrain microstructures, the subgrain structures are
predicted from crystal orientation and stored energy calculated by the deformation simulation. The
multi-phase-field model of Steinbach et al. [9] is employed as a grain growth model and the misorientation
dependency of grain boundary energies and mobility are introduced into the phase-field model.

This paper consists of five chapters: Following this introduction, a multi-phase-field model for grain
growth and algorithm for the efficient computation are described. Next, we explain the procedure ofrecrystallization model developed in this study. In chapter of numerical results, first, the effects ofmisorientation and inhomogeneous deformation on subgrain growth are confirmed by using simple
polycrystalline model, and then the results for every step following the developed recrystallization model
are shown in detail. And, we discuss the capability of the developed model against the predictability of realprocesses. Finally we conclude this study.
PHASE-FIELD MODELIn grain growth simulation during recrystallization, the multi-phase-field model of Steinbach et al. [9] isemployed. This model has an advantage that the phase field parameters, such as gradient coefficient andphase field mobility, can be perfectly related to the grain boundary energy and mobility withmisorientation dependency. Furthermore, this model removes some difficulties in the treatment of triple
points or higher order interactions. The computational difficulties caused by using many phase field
variables are overcome by employing an efficient algorism of Kim et al. [10].
1. Multi-phase-field model Considering a system containing N different grains φ1, φ2,.... φN, the freeenergy functional can be taken as⎡ N N ⎛ a2⎞⎤
F = ∫ ⎢∑ ∑ ⎜⎜ − αβ ∇φα ⋅ ∇φβ + Wαβ φαφβ ⎟⎟⎥dV ,V⎠⎦⎥⎣⎢ α =1 β =α +1 ⎝ 2(1)where, aαβ and Wαβ are the gradient coefficient and the height of double well potential between α-th andβ-th grains. The phase field φα takes a value of 1 inside the α-th grain, 0 inside other grains, and 0 < φα < 1
at grain boundary. φα is not an independent variable and must satisfy the following condition:N

φα = 1 .
∑
α

(2)

=1

Here, we define a step function σα with σα = 1 if 0 < φα < 1 and σα = 0 elsewhere. Using this step function,
the number of locally present phases n is expressed asNn = ∑ σ α ( x, t ) .

(3)

α =1

Using Eq. 3, N in Eqs. 1 and 2 can be replaced by n.
Following Ref. 9, the time evolution equation of phase field is obtained asn Mφn
M ijφ ⎛ δF δF ⎞⎜⎟ = −∑ ij ⎡⎢∑ ⎧⎨(Wik − W jk )φk + 1 aik2 − a 2jk ∇ 2φk ⎫⎬⎤⎥ ,−⎜δφ j ⎟⎠
2⎭⎦
j =1 n ⎝ δφi
j =1 n ⎣ k =1 ⎩(nφ&i = −∑)

(4)

where i, j and k are the local grain number from 1 to n at a numerical grid point.
2. Phase field parameters The phase-field parameters are related to the material parameters as follows
[11]:
4γ
2
π2
2δγ ij and M ijφ =Wij = ij , aij =M ij ,(5)πδ
4δ
where δ is the grain boundary thickness and γij and Mij are the grain boundary energy and mobility between
i-th and j-th grains. To introduce the misorientation dependency of grain boundary energy and mobility,
the following equations are used [8].

γ ij (Δθ ij ) = γ m

Δθ ij ⎛Δθ ⎞⎜⎜1 − ln ij ⎟⎟ .
Δθ m ⎝
Δθ m ⎠(6)⎡⎧⎪ ⎛ Δθ ij ⎞ 4 ⎫⎪⎤⎟⎟ ⎬⎥ .
M ij (Δθ ij ) = M m ⎢1 − exp ⎨− 5⎜⎜θΔ⎢⎪⎩ ⎝ m ⎠ ⎪⎭⎥⎦⎣

(7)

Here, Δθij is the misorientation between i-th and j-th grains, γm and Mm are respectively the boundary
energy and mobility at Δθij = Δθm which is the misorientation when the boundary becomes a high angle
boundary. The variations in grain boundary energy and mobility with misorientation expressed by Eqs. 6
and 7 are illustrated in Fig. 1.

Fig. 1 Misorientation dependency of grain boundary energy and mobility expressed by Eqs. 6 and 7
3. Algorithm for efficient computation As we can see from Eq. 4, we do not have to solve Eq. 4 on thegrid points with n = 1 and it is sufficient to save the values of not N but n phase fields. Furtheremore, byintroducing the following algorithm [10], we can solve Eq. 4 very efficiently.
(1) Solve Eq. 4 for n phase fields.
(2) During (1), if φi at previous time step is zero and the calculated increment Δφi has a minus value, i-thphase field is put away from the phase field group on the grid point. And, steps (1) and (2) are repeateduntil all phase fields satisfied above condition are removed.
(3) The phase fields are rearranged in the order of largest to smallest. Phase fields from the largest until
nφ -th are recorded. Here, nφ is a predefined maximum number of the recorded phase field.

(4) The phase fields are replaced by following φi* = φi

nφ

nφ

j =1

i =1

∑φ j so as to satisfy ∑ φi* = 1 .

(5) If a phase field is not saved at a grid point of (l, m) and its value at grid points of (l±1, m±1), or thenearest four neighbors, is not zero, the phase field is added to the phase field group on the grid (l, m).
Therefore, the number of phase fields solved in the step (1) sometime becomes more than nφ.
Although this algorithm basically follows Ref. [10], the iteration in steps (1) and (2) is added to achieve
higher accuracy.
RECRYSTALLIZATION COUPLED MODEL

When we obtain the stored energy Estore from subgrain structures measured experimentally, such as EBSD
method, we often use the equation
Estore =

Kγ s
,
Ds

(8)

where γs is the grain boundary energy, Ds is the subgrain diameter and K is a constant depending on thespace dimension [12]. Since, in crystal plasticity finite element simulation, the stored energy can becalculated by

Estore =1
ρμb 2 ,
2(9)

where ρ is the dislocation density, μ is the shear modulus and b is the magnitude of Burgers vector, thesubgrain diameter can be calculated as
Ds =

Δθ ⎛
Δθ ⎞
K
K
⎜⎜1 − ln⎟,γs =γm
Δθ m ⎝
Δθ m ⎟⎠
EstoreEstore(10)where Eq. 6 is substituted. However, usually it is impossible to obtain the misorientation Δθ directly fromcrystal plasticity finite element simulation. Furthermore, in the case using triangular element in finite
element simulation, the stored energy and crystal orientation are constant inside individual element andabruptly change between neighboring two triangular elements as a step function. Therefore, the smoothing
operation is carried out using the equationqθ (l , m ) = ∑ θ p w p
p =1

q

∑w ,
p =1

(11)

p

where θ(l, m) is the smoothed orientation at a grid point (l, m) shown in Fig. 2, q is the number of grid points
used in the smoothing operation (or the number of grid points located inside dashed circle with radius Rs in
Fig. 2), θp is the crystal orientation at local point p, and wp is a weight function expressed by 1/ rp , where rpis the distance between points (l, m) and p. The same operation is performed for stored energy. From thesmoothed orientation, the amount of local orientation gradients ∇θ are calculated at all grid points.
Assuming Δθ = ∇θ Ds as shown in Fig. 3, from Eq. 10, the misorientation Δθ can be calculated using the
local orientation gradient ∇θ and smoothed stored energy Estore as⎡⎧
K γ m ⎫⎤
Δθ = Δθ m exp ⎢1 − 1 ⎨ ∇θ⎬⎥ .
Estoer Δθ m ⎭⎦
⎩
⎣

(12)

Here, we summarize the recrystallization model developed in this study:(1) Crystal orientation and stored energy of deformed polycrystalline material are calculated throughcrystal plasticity finite element simulation.
(2) The crystal orientation and stored energy are smoothed on the regular grids for phase-field
simulation following Eq. 11.
(3) The amount of local orientation gradient ∇θ =(∂θ ∂x )2 + (∂θ ∂y )2 is calculated on all grid

points.
(4) The misorientation Δθ is calculated by Eq. 12 using the smoothed stored energy and local orientation
gradient.
(5) Subgrain size at all grid points is determined using Eq. 10.
(6) Nuclei to generate the subgrain structure are seeded on the regular grids for phase-field simulation
taking into account the calculated subgrain size.
(7) Subgrain structures are prepared by performing a normal grain growth simulation from the nuclei.
(8) Recrystallization phase-field simulation is carried out from the initial subgrain structure.

Fig. 2 Schematical illustration of relationship between triangle finite element and regular grids

Fig. 3 Schematical illustration of relationship between smoothed orientation and predicted subgrain

NUMERICAL RESULTS
1. Abnormal subgrain growth The effects of misorientation and inhomogeneous deformation on thesubgrain growth are investigated using simple polycrystalline model. Figure 4 show the initial subgrain
structures used here. There are four hexagonal grains that are classified by coloring and of whichboundaries are high-angle boundaries. Here, periodic boundary conditions are applied to the all sides ofcomputational domain. Model A has 718 subgrains with almost same diameter of Dave = 0.94 μm. Models
B and C have two kinds of subgrain size. In model B with 657 subgrains, only central red grain has small
subgrains of Dave = 0.80 μm and the other three grains have large subgrains Dave = 1.07 μm. On the other
hand, in model C with 846 subgrains, central red grain has large subgrains of Dave = 1.07 μm and the othergrains have small subgrains of Dave = 0.80 μm. The computational domain size is 500Δx ×394Δx where
grid size is Δx = 0.05 μm. The other conditions are as follows: the interface thickness δ = 7Δx, interface
energy γ = 1 J/m2, time increment Δt = 0.0075 s and grain boundary mobility M = M0exp(-Q/kT) where
pre-exponential factor M0 = 6.2x10-6 m2/Ns, activation energy Q = 2.08x10-19 J, Boltzmann constant k =1.38x10-23 J/K and temperature T = 800 K [13].
Figures 5 (a), (b) and (c) show the time slices of subgrain microstructual evolution for Model A withdifferent misorientations Δθ. In Fig. 5 (c), normal grain growth is observed, because all boundaries are
high-angle boundary. On the other hand, in Fig. 5 (a), only grain boundaries with high-angle
misorientation migrate, or the abnormal grain growth is observed. In Fig. 5 (b), intermediate behaviors of(a) and (c) can be observed: although the migrations of initial high-angle grain boundaries are prominent,
the coarsening of subgrains inside initial grain is also occurred.
Figures 6 (a) and (b) demonstrate the time evolution of subgrain structures for Model B and C,
respectively. In Models B and C, the deformation inhomogeneity is simply modeled by changing thesubgrain size. Although , in Fig. 5 (a), the high-angle grain boundaries migrate into both grains with same
subgrain structures, in Figs. 6 (a) and (b), those migrate only into the grain with smaller subgrains or larger
deformation. Therefore, we can confirm that the present model can reproduce the experimental fact thatrecrystallization occurs much faster in highly deformed grain than lower one.

(a) Model A(b) Model B(c) Model C
Fig. 4 Initial subgrain structures with four hexagonal grains

(a) Δθ = 4 degree(b) Δθ = 7 degree(c) Δθ > 15 degree
Fig. 5 Time evolution of subgrain structures for Model A. t = 7.5, 75, 150, 225, 300 and 375 s from left.

(a) Model B(b) Model C
Fig. 6 Time evolution of subgrain structures for models B and C in the case of Δθ = 4 degree. t = 7.5, 75,
150, 225, 300 and 375 s from left.
2. Recrystallization simulation Following the procedure of recrystallization model described in previous
chapter, recrystallization simulation taking into account deformation microstructures is performed and theresults in each computational step are shown in detail.

Fig. 7 Crystal orientation distributions at 50% compression by crystal plasticity finite element
simulation [14]Figure 7 shows the crystal orientation distributions at 50% compression calculated by crystal plasticity
finite element simulation. The used crystal plasticity theory and numerical conditions are shown in Ref.
[14]. In Fig. 7, the solid lines indicate the grain boundaries before deformation. Here, we use the squareregion of 80 μm x 80 μm in Fig. 7 as a computational domain for recrystallization simulation. This region
has 626 triangle elements and is divided into 800 x 800 finite difference grids. The close-up views of squareregion in Fig. 7 are illustrated in Fig. 8. Since, in crystal plasticity finite element simulation, crossed
triangle elements are employed, (a) crystal orientation and (b) stored energy are colored element byelement. The misorientation shown in Fig. 8 (c) indicates the orientation difference of neighboring two
triangle elements.

(a)
(b)
(c)
Fig. 8 Close-up views of square region in Fig. 7: (a) crystal orientation θ, (b) stored energy Estore and(c) misorientation Δθ
Figures 9 (a) and (b) are the smoothed orientation and stored energy, respectively. The smoothing is carriedout using Eq. 11 with Rs = 3Δx. In this case, since the Rs is relatively small, the smoothed regions arelimited around sides of triangular elements. Next, the magnitude of orientation gradient at grid point (l,m)
1(θ(l +1,m ) − θ(l −1,m ) )2 + (θ(l ,m+1) − θ(l ,m−1) )2 using smoothed orientation. Then,is calculated by ∇θ (l ,m ) =2Δxusing smoothed stored energy shown in Fig. 9 (b) and calculated gradient ∇θ , the misorientations on all
grid points can be obtained from Eq. 12. The calculated misorientation distributions are shown in Fig. 9 (c).
Here, the minimum misorientation is set to be 4 degrees, because we need very small numerical grid inorder to reproduce small subgrain structures caused by small misorientation angle, as can be seen from Eq.
10. Comparing Fig. 9 (c) to Fig. 8 (c) , the area with high-angle misorientation increases in Fig. 9 (c). As thereason, it is considered that the highly localized orientation gradient at around sides of triangular elements
causes such high misorientation through Eq. 12. Finally, the subgrain size on all grid points are determinedby Eq. 10 using smoothed stored energy and calculated misorientation. The subgrain diameters calculated
at all grid points are the range 0.97 to 3.46 μm.

(a)
(b)
(c)
Fig. 9 Smoothed results of (a) crystal orientation θ and (b) stored energy Estore and (c) misorientationΔθ calculated by Eq. 12
The nuclei are sowed so as that the circular subgrains with diameter Ds of grid points selected randomly fill
the computational domain. The distributes of created 2229 nuclei are shown in Fig. 10 (a). From these
nuclei, the normal grain growth is simulated to obtain the subgrain structures. Figures 9 (a), (b) and (c) are
the time slices during the simulation. Figures 9 (c) and (d) are the predicted subgrain structures, where (d)indicates grain boundaries and crystal orientations. From Fig. 9 (d), the crystal orientations are almost samewith Fig. 9 (a). The small subgrains are observed at the area where the stored energy is high. Therefore, it isconfirmed that the plausible subgrain structures can be predicted by the proposed model.

(a)
(b)
(c)
(d)
Fig. 10 (a) Nuclei sowed randomly, (b) during grain growth simulation to prepare subgrain structure,(c) predicted subgrain structures and (d) its crystal orientations
The static recrystallization simulation is carried out from the predicted subgrain structures shown in Fig. 9(d). The computational domain size is 80 μm x 80 μm (800 x 800 grids) with grid size Δx = 0.1 μm. Zero
Neumann boundary conditions are applied for all sides of computational domain. The other conditions are
as follows: the interface thickness δ = 7Δx, interface energy γ = 0.6 J/m2, time increment Δt = 0.05 s andgrain boundary mobility M = M0exp(-Q/kT) where pre-exponential factor M0 = 6.2x10-6 m2/Ns, activationenergy Q = 2.08x10-19 J, Boltzmann constant k = 1.38x10-23 J/K and temperature T = 800 K [13].
Figure 11 shows the time evolution of subgrain structures and the growth of recrystallized grains during
annealing. At the beginning of simulation, the high-angle grain boundaries become smooth owing to theircurvature, as can be seen in Fig. 11 (a) with comparing to Fig. 10 (d). After then, some grains withhigh-angle boundaries grow abnormally and become recrystallized grains. In particular, large
recrystallized grains grow into upper and lower grains, because their two grains are deformed intensely andhave small subgrain structures. The large deformation area on the diagonal line of computational domain isalso preferred nucleation site. However, since the deformation is localized around the diagonal line, thegrowth rate of recrystallized grains is slower than those of upper and lower regions. In the large initial
subgrain regions in Fig. 10 (d), no recrystallization occurs in this simulation period. As a result, thespontaneous nucleation by abnormal grain growth and the recrystallized grain growth depending on the
deformation structures are observed in the present model developed here.

(a)
(b)
(c)
(d)
(e)
Fig. 11 Time evolution during static recrystallization simulation at (a) 300, (b) 1500, (c) 3000, (d) 4500
and (e) 6000 s
CONCLUSIONS

We developed a coupled model for static recrystallization of the phase-field method and crystal plasticity
finite element method. In this model, the subgrain structures are predicted from crystal orientation andstored energy calculated by crystal plasticity finite element simulation. The multi-phase-field model
proposed is employed as a grain growth model and the misorientation dependency of grain boundary
energy and mobility are introduced into the phase-field model to reproduce abnormal grain growth whichenable spontaneous nucleation. By performing a series of procedure, it is confirmed that the proposedmodel can predict the desired subgrain structures depending on the stored energy and crystal orientation
distribution. Furthermore, during static recrystallization simulation, spontaneous nucleation by abnormal
grain growth and the growth of recrystallized grain depending on the deformation microstructures are
confirmed.

Acknowledgements Financial support from the Ministry of Education, Culture, Sports, Science andTechnology of Japan through a Grant-in-Aid for Scientific Research is highly appreciated.
REFERENCES[1] M. A. Miodownik, A review of microstructural computer models used to simulate grain growth andrecrystallisation in aluminium alloys, J. Light Metals Vol. 2-3 (2002), 125-135.
[2] D. J. Srolovitz, G. S. Grest, M. P. Anderson, Computer simulation of recrystallization—I.
Homogeneous nucleation and growth, Acta metal., 34, (1986), 1833-1845.
[3] H. W. Hesselbarth, I. R. Göbel, Simulation of recrystallization by cellular automata, Acta Metall.
Mater., 39, (1991), 2135-2143.
[4] T. Takaki, Y. Yamanaka, Y. Higa, Y. Tomita, Phase-Field Model during Static Recrystallization
based on Crystal-Plasticity Theory, Scientific Modeling and Simulation, (2008), accepted.
[5] Y. Suwa, Y. Saito, H. Onodera, Phase field simulation of stored energy driven interface migration ata recrystallization front, Mater. Sci. Eng. A, 457, (2007), 132-138.
[6] D. Raabe, L. Hantcherli, 2D cellular automaton simulation of the recrystallization texture of an IF
sheet steel under consideration of Zener pinning, Comp. Mater. Sci., 34, (2005), 299-313.
[7] D. Raabe, R. C. Becker, Coupling of a Crystal Plasticity Finite Element Model with a Probabilistic
Cellular Automaton for Simulating Primary Static Recrystallization in Aluminum, Model. Simul.
Mater. Sci. Eng., 8, (2000), 445-462.
[8] F. J. Humphreys, A unified theory of recovery, recrystallization and grain growth, based on thestability and growth of cellular microstructures—I. The basic model, Acta Mater., 45, (1997),
4231-4240.
[9] I. Steinbach, F. Pezzolla, A generalized field method for multiphase transformations using interface
fields, Physica D, 134, (1999), 385-393.
[10] S. G. Kim, D. I. Kim, W. T. Kim, Y. B. Park, Computer simulations of two-dimensional andthree-dimensional ideal grain growth, Phys. Rev. E, 74, (2006), 061605-1-061605-14.
[11] T. Takaki, A. Yamanaka, Y. Tomita, Phase-Field Modeling for Dynamic Recrystallization,
Proceeding of International Conference on Advances and Trends in Engineering Materials and their
Applications, (2007), in press.
[12] F. J. Humphreys, M. Hatherly, Recrystallization and Related Annealing Phenomena, Elsevier,(2004).
[13] D. Raabe, Cellular automata in materials science with particular reference to recrystallization
simulation, Ann. Rev. Mater. Res., 32, (2002), 53-76.
[14] T. Takaki, A. Yamanaka, Y. Higa, Y. Tomita, Development of Phase-Field Model andComputational Procedure During Static Primary Recrystallization, Trans. J. Soc. Mech. Eng., 73,(2007), 482-489. (in Japanese)