Fast simulation of grain growth based on Orientated Tessellation Updating Method

. This work is part of a more general idea consisting in developing a macroscopic model of grain growth whose state variables contain for each material point the statistical descriptors of the microstructure (e.g., disorientation, grain size and shape distributions). The strategy is to determine macroscopic free energy and dissipation potentials on the basis of a large number of computations at the scale of the polycrystal. The aim is to determine enriched macroscopic evolution laws. For sake of simplicity, this contribution only deals with grain growth of a single phased metal without diffusion or segregation of alloying elements. In order to test this upscaling strategy it is necessary to establish a simulation tool at the scale of the polycrystal. It should be sufficiently simple and fast to enable a large number of simulations of various microstructures, even if it leads to neglect some phenomena occurring at this scale. Usual grain growth models relying on mobile finite element modeling, level set functions, phase field or molecular dynamics are too computationally costly to be used within the proposed framework. Therefore, this paper focuses on the development of a “toy” model. Tessellation techniques are usually used to approximate polycrystalline microstructures. Therefore, one can approximate the real evolution of the microstructure as a succession of tessellation approximations. It then becomes quite natural to attempt to establish the evolution law of the microstructure directly on the parameters defining the tessellation. The obtained model is very light in terms of computational cost and enables to compute a large number of evolutions within the framework of the proposed statistical upscaling method.


Introduction
This work is part of a more general idea consisting in developing a macroscopic model of grain growth whose state variables contain for each material point the statistical descriptors of the microstructure (e.g., disorientation, grain size and shape distributions).Statistical distributions are too much detailed to be processed at each material point, thus reduced information is considered instead: means and variances at each material point.This ambition arises as very few information about microstructure is usually processed at the macroscopic scale.The general framework is the standard generalized media characterized by a free energy per unit mass and a dissipated power potential.These two potentials arise in the balance energy equation combining the first and second laws of thermodynamics.The balance energy equation should be verified for all possible virtual evolution.The two potentials depend on macroscopic state variables and their time * e-mail: weisz@lms.polytechnique.frderivatives.The determination of the macroscopic free energy and dissipation potentials enables to establish the evolution laws of the state variables characterizing the microstructure evolution at macroscale.
The strategy is to determine macroscopic free energy and dissipation potentials not axiomatically (with parametric functions) and calibration with experiments at macroscale, but on the basis of a large number of computations at the scale of the polycrystal.The aim is to quantify for each microstructure evolution the total grain boundary energy and the total dissipated power as a function of the macroscopic state variables that characterize the statistical distribution of the microstructure in order to determine enriched macroscopic evolution laws.For sake of simplicity, this contribution only deals with grain growth of a single phased metal without diffusion or segregation of alloying elements.In order to test this upscaling strategy it is necessary to establish a simulation tool at the scale of the polycrystal.It should be sufficiently simple and fast to enable a large number of simulations of various microstructures, even if it leads to neglect some phenomena occurring at this scale.Usual grain growth models relying on cellular automaton [1][2][3][4], mobile finite element modeling [5], level set functions [6,7], phase field [8][9][10][11] or molecular dynamics [12][13][14] have been intensively studied and very interesting results have been obtained.However, the computational cost of such approaches is usually incompatible with an intensive use as suggested within the proposed framework.Vertex methods [15][16][17][18] consist in establishing the evolution law directly at the triple junctions and are sufficiently simple in two dimensions to reach short computation time.However, the extension in three dimensions is very difficult.
Evolution of polycrystalline microstructures (such as grain growth) involves coupled mechanisms at different scales: (i) the atomic scale (crystal lattice), (ii) the microscale including mechanisms at grain boundaries (iii) the mesoscale revealing the grain structure ( i.e., polycrystalline Representative Volume Element) and (iv) the macroscopic scale modeled as a continuous medium.Thus, this paper focuses on the development of a "toy" model formulated at mesoscale and enabling to simulate grain growth within short computation time.
Voronoi-Laguerre tessellation techniques are usually used to approximate polycrystalline microstructures at mesoscale.Very efficient algorithms have been developed with the possibility of controlling statistical distributions of grain size and shape.Crystal lattice orientation can also be specified for each grain.The tessellation equipped with such an orientation field is called an Orientated Tessellation (OT).One can approximate the real evolution of the microstructure as a succession of OT approximations.It then becomes quite natural to attempt to establish the evolution law of the microstructure directly on the parameters defining the OT.This modeling strategy is called in this contribution the Orientated Tessellation Updating Method (OTUM).
As grain growth during annealing is essentially viscous, the time derivatives of the OT parameters are obtained as a function of the thermodynamic forces defined as the partial derivatives of the total grain boundary energy with respect to the OT parameters.This mesoscale evolution law is obtained through the energy balance equation by specifying two mechanisms at microscale: (i) the grain boundary energy (disorientation) and (ii) the dissipated power by crystal visco-plasticity through any grain boundary virtual motion.For sake of simplicity, crystal lattice rotation through crystal visco-plasticity is not considered in this contribution.In addition, plane polycrystals are considered with three slip directions in each grain.Thus, the crystal lattice is plane hexagonal.This configuration corresponds to the plane 1, 1, 1 of a face-centered cubic (fcc) crystal, as shown in Figure 1.Disorientation between two neighboring grains (characterized by five parameters in 3D) is characterized only by two parameters in 2D: the disorientation angle ∆θ and the orientation of the grain boundary ϕ.Thus, the grain boundary energy considered in this paper is computed from fcc crystals sharing the same orientation 1, 1, 1 (asymmetric tilt boundaries).
The dissipative mechanism at microscale is calibrated to correspond to the classic notion of grain boundary mobility.The overall mesoscale grain growth model is validated by comparison with the classical curvature driven shrinkage of a spherical grain.The proposed model is very light in terms of computational cost and enables to compute a large number of microstructure evolutions, which is useful for the general statistical upscaling method.

Orientated Tessellation Updating Method
A tessellation is defined by n seeds whose Cartesian coordinates are denoted by (x j , y j ) and n weights denoted by w j (where 1 ≤ j ≤ n).Coordinates (x j , y j ) ∈ [0, 1] 2 are dimensionless as well as weights.The tessellation is completely determined by the parameter vectors The dimensionless tessellation is scaled by a physical length L 0 that represents the length of the tessellation edge.
Each cell (or grain) denoted by C j (where 1 ≤ j ≤ n) is defined by Voronoi-Laguerre tessellation as follows: It is clear from the definition (1) that weights are defined up to a constant.Thus, the following constraint is added to obtain a univocal definition: Thus, weights w lie in an affine hyperplane of dimension n − 1 and denoted by P whose support is the hyperplane denoted by P (n−1) = w ∈ R n , n j=1 w j = 0 .In addition, it should be noted that a cell C j may be empty as shown in Figure 2.This property will be intensively used as some grains should disappear during grain growth.

The crystallographic orientation is denoted by
Since the crystal lattice is plane hexagonal ∀j ∈ {1, . . ., n} , θ j ∈ 0, π 3 .The parameter set P OT defining the OT (see Fig. 3) reads: (3) A total energy E(T, α) is associated to the OT, where T is the temperature.In this contribution, the total grain boundary energy per unit depth is considered.The Orientated Tessellation Updating Method consists in computing a time evolution of the OT characterized by P OT by deriving an evolution law of the form: where M (T, α) is a second order tensor to be defined by introducing a dissipative mechanism at microscale and ∂E(T, α)/∂α is the driving force.
This paper is limited to the evolution of weights (i.e., seeds and orientations are fixed parameters: ẋ = ẏ = 0 and θ = 0), that is to say that the generalized vector space of speed V * reads:

Microscale mechanisms
As mentioned in the introduction, the mesoscopic evolution law ( 4) is determined by introducing two local mechanisms at microscale: (i) the grain boundary energy and (ii) the dissipated power by crystal visco-plasticity through any grain boundary virtual motion.Grain boundaries are indexed by pairs (i, j) where i and j denote two neighboring grains.The set of pairs of neighboring grains defining grain boundaries is denoted by I GB : The condition j ≥ i is meant to count each grain boundary only once (otherwise if (i, j) ∈ I GB then (j, i) would also be an element of I GB ).
The grain boundary energy per unit area of each grain where ∆θ ij is the crystal lattice disorientation between grains i and j defined by: In addition ϕ ij is the angle of the grain boundary as shown in Figure 4.
Several analytical formulations have been proposed for E. The simplest and most widely used function is the Read & Shockley formula [19]: However, the range of validity of the Read & Shockley formula is limited to small disorientation angles and do not account for the energy cusps at certain disorientation angles.To overcome this difficulty, Wolf [20] introduced a piecewise interpolation function for the grain boundary energy per unit area of the form: (10) This interpolation function is usually fitted on atomistic simulations.Other interpolations have been proposed within a three dimensional framework [21].Since the main goal of this contribution is to demonstrate the capability of the OTUM and not to present a study on a specific metal, a simplified grain boundary function is directly defined in the form of (10).Since the crystal lattice is plane hexagonal, a cusp is introduced at π/3.In addition it is assumed that the dependance on ϕ is negligible.The grain boundary energy per unit area is defined as follows: where G(T ) is the shear coefficient.In this paper, data for pure iron [22] are used to calibrate G(T ): where b G ≈ 88134 MPa and a G ≈ −24 MPa.K −1 .The temperature dependence G(T )/G(0) is a simple choice similar to [23].In addition: ) where a 1 = a 2 = 1 and γ π 6 = 1 J.m −2 , which roughly corresponds to pure iron.The resulting simplified grain boundary energy per unit area is presented in Figure 5.
In addition to the grain boundary energy per unit area, the dissipated power per unit area through any virtual grain boundary motion should be specified.Consider v * ij a virtual normal velocity of the grain boundary (i, j) ∈ I GB .The dissipated power per unit area reads: where [24,25] used crystal plasticity to determine D analytically (for plane hexagonal crystals) when the boundary between two semi-infinite crystals moves at the virtual velocity v * ij : where τ c is the critical shear stress and X : ∆θ → X(∆θ) is the following function: (16) Since grain growth is viscous, the analytical computation proposed by [24,25] is simply adapted for crystal visco-plasticity by considering that the critical shear stress τ c linearly depends on v * ij : where the function m(T ) is homogenous to a grain mobility (m 4 .J −1 .s−1 ).Hence the dissipated power per unit area: It should be noted from (10) and ( 16) that: A∆θ ln (∆θ) In other words, the driving force (proportional to γ(∆θ)) tends to zero when the disorientation tends to zero.However, the dissipative cost (proportional to X(∆θ)) that resists to the driving force, tends to zero much faster because of the logarithmic term in (19).Thus, for arbitrarily small disorientation, the grain boundary is unstable as the driving force is arbitrarily large in comparison to the resistive mechanism (the ratio is in ln (∆θ)).

Mesoscale evolution law
The mesoscale evolution law is derived by considering the total grain boundary energy per unit depth and the total dissipated power per unit depth through any virtual grain boundary normal velocity field, which are defined as follows: where l ij the joint length and v * = v * ij (i,j)∈I GB .The dependance on α = (x, y, w, θ) comes from the grain boundary lengths l ij that depend on (x, y, w) and E ij , D * ij that depend on θ.The number of grain boundary is There exists a bijection between I GB and {1, • • • , n GB } that defines the numbering of grain boundaries.Thus, for each grain boundary denoted by (i, j) ∈ I GB there exists a unique Therefore the total dissipated power per unit area reads: where χ(α) is a dimensionless diagonal second order tensor of size n GB × n GB of diagonal: As seeds and crystalline orientations are fixed (i.e., ẋ = ẏ = 0 and θ = 0) and considering any virtual variation of weights ẇ * , it is straightforward to demonstrate that: where d ij is the dimensionless distance between seeds: It is clear in (23) that for each grain boundary an arbitrary choice is made for the positive direction of the normal velocity v * ij , which has no consequence as the square of the virtual velocity arises in the dissipated power.Hence: where K(α) is a dimensionless second order tensor of size n GB × n, which can be evaluated analytically:

See equation (26) next page
The constraint (2) leads to: which defines a hyperplane of dimension n − 1 denoted by P (n−1) , whose normal vector is denoted by It is clear that the kernel of K(α) is ker (K(α)) = 1 R.The total dissipated power per unit depth reads: where R(α) = K(α) T .χ(α).K(α) is a dimensionless second order tensor of size n × n.
The "real" weight variation ẇ is determined for isothermal (i.e., Ṫ = 0) and homogenous temperature fields (i.e., ∇T = 0) by using the general energy balance equation (combining the first and second thermodynamic laws), which is established for "real" evolution (and not virtual): Hence: Thus, the energy balance (30) arises as a constraint on the "real" weight variation.The maximum dissipation principle [26], enables to determine the "real" weight variation as the maximum of the dissipation under the (31) Consider µ the Lagrangian multiplier and L the Lagrangian: (32) The optimality condition reads: (33) Hence µ = −2 and the evolution law reads: The evolution law (34) should be inverted.However R(α) is not invertible because of the constrain (27).The kernel of R(α) is ker (R(α)) = 1 R. Thus, there are n − 1 strictly positive singular values of R(α) denoted by (λ 1 , • • • , λ n−1 ) and there exist two orthogonal tensors U and V such as: The singular value decomposition is computed in Scilab [27] and gives (λ 1 , • • • , λ n−1 ), U and V .The Moore-Penrose pseudo-inverse matrix is introduced as follows: Hence the evolution law: where: It should be noted that a scale effect is demonstrated in the evolution law (37) as the physical dimension L 0 arises.

Driving force
The driving force involved in the evolution law (37) should be evaluated.An analytic solution is derived in this section, which contributes to reach short computation time.From ( 20) and ( 11) one obtains: Hence: Consider the set of triple junctions: There is no condition such as k ≥ j ≥ i, thus if (i, j, k) ∈ I 3 then all permutations are also in I 3 .Consider the third order tensor δ: Thus, the third order tensor β representing angles at the triple junctions as shown in Figure 6 are defined as follows: where the following symmetry rule holds β ijk = β kji .Consider a virtual weight variation ẇ * ∈ P (n−1) (i.e., verifying the constrain ( 27)).One obtains: Hence by combining (44) and (45) and ∂l ij /∂w = (∂l ij /∂w q ) q∈{1,••• ,n} : where δ is the Kronecker symbol.

Calibration
The proposed approach is calibrated by comparison with the classic curvature driven grain growth.Within this framework, the driving force acting on any grain boundary reads: where R I and R II are the radii of curvature along the principal directions (in three dimensions).For plane polycrystals R I = R and R II = ∞.A simple linear evolution law accounting for viscosity is usually used: where v is the inward normal speed of the considered grain boundary and m(T, ∆θ) the classical grain boundary mobility for the curvature driven model (where the symbol .refers to the curvature driven model).In the following, a circular grain shrinkage is compared to an hexagonal configuration as shown in Figure 7 in order to determine the mobility m(T ) introduced in the dissipated power (17) as a function of m(T, ∆θ).Thus for the curvature driven model the inward normal speed reads: The proposed approach is applied to the hexagonal situation.The partial derivative (46) reads: Hence: The driving force (40) reads: (52) The second order tensor K introduced in (26) reads for this configuration: The second order tensor χ introduced in ( 22) reads: where I is the identity of size 6 × 6.The second order tensor R introduced in (28) reads: Thus, the second order tensor M introduced in (38) reads: Hence the evolution law (37) for the hexagonal situation: (57) Thus, by using (23) one obtains the inward normal speed according to the proposed model: The comparison between the inward normal speed from the curvature driven model (49) and from the proposed For the sake of simplicity it is assumed that m(T, ∆θ) depends on disorientation in 1/X(∆θ) even though more realistic mobilities could be obtained by considering m(T, ∆θ) instead of m(T ) in (17).
To estimate m(T ) one can use the pure iron grain boundary mobility determined by [28]: where the molar volume of fcc-Fe is V m = 7.09 × 10 −6 m 3 , the boundary thickness is δ = 1 nm, the burgers vector is b = 2.48 × 10 −10 m, the gas constant is R = 8.3144621 J.mol −1 .K −1 .In addition, D Fe (T ) is the diffusivity of Fe atoms along the grain boundary (for fcc) and follows an Arrhenius law: where the pre-factor is D 0 Fe = 1.5 × 10 −4 m 2 .s−1 the activation energy is Q = 148 kJ.mol −1 .

Time discretization and numerical results
The temperature kinetics is imposed in the following, thus T (t) is a known function.Voronoi-Laguerre tessellations (x, y, w) are generated with the free software Neper [29] and orientations θ are set in Scilab [27].Weights are initially set to a common value and are updated at each time step according to the evolution law (37) and the driving force (46).The numerical discretization consists in a simple explicit scheme with adaptive time step in order to capture accurately growth rate variations.The time step is determined as follows: where δ w > 0 is a constant value.As the simple model proposed in this paper is purely viscous, even though the temperature is maintained constant, grain growth never stops strictly speaking.The evolution becomes only slower and slower.Thus, a stopping criterion needs to be defined.The algorithm is simply presented in Figure 8. Different Orientated Tessellations with different distributions of sphericity and grain size are generated.The orientation field correspond to a disorientation distribution very similar for each OT as shown in Figure 9    grain size arrangement leading to heterogeneous grain size distribution during the transient state.Of course, grain size heterogeneity eventually reduce as bigger grain tends to grow at the expense of smaller grains.In addition, the disorientation distribution tends to favor small disorientations as shown in Figure 9. Sphericity tends to increase as shown in Figure 11 but this may be due to the fact that weights w are the only state variables.Future developments will include seeds x and y as state variables, which will enable to consider additional degrees of freedom to approximate more accurately grain shape evolution.In addition, orientations θ can also be added as state variables with an associated dissipative mechanism related to crystal plasticity in the grain bulk.Thus, reorientation is likely to occur for small grain sizes as suggested by molecular dynamics [14] because reorientation involves dissipated energy per unit volume although grain boundary migration involves dissipated energy per unit area.
Grain growth is a thermally activated process as captured by the Arrhenius law (62).This is illustrated by simulating the same OT as those presented in Figure 10 at 700 • C. The resulting evolution is much slower than at 800 • C as shown in Figure 13.

Conclusion
In this paper the Orientated Tessellation Updating Method has been introduced to simulate grain growth at mesoscale (i.e., at the scale of the polycrystalline structure).The method relies on Voronoi-Laguerre tessellation techniques and energetic contributions at microscale (i.e., at the scale of the grain boundary), namely the grain boundary energy due to the crystal disorientation and the dissipated power by crystal plasticity for any virtual motion of the grain boundary.The method is very light in terms of computation cost.Even though 2D structures have been studied in this contribution, the extension in 3D is straightforward and future works will focus on this aspect.In addition, the only state variables are the weights of the tessellation.Thus, further developments including seeds and orientations as state variables are needed to approximate more accurately structural evolution of polycrystals.Furthermore, tessellation techniques include options to generate subgrain structures that could also be introduced in the framework presented in this paper.
. The temperature is set to T = 800 • C and grain growth is simulated by Orientated Tessellation Updating Method during approximately 4 h.The code is not optimized as Scilab has been used to compute the evolution law and Neper to generate tessellations.(Computation time is roughly 15 min for 1000 time steps and is mainly due to writing and reading text files between Scilab and Neper.Thus, a unifed code not relying on text files would significantly reduce computation time.)