© D. Weisz-Patrault et al., published by EDP Sciences 2020

1 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 derivatives. 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–4], mobile finite element modeling [5], level set functions [6,7], phase field [8–11] or molecular dynamics [12–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–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> (asymmetrictilt 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.

Fig. 1 FCC and slip directions in <1, 1, 1> plane.

2 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)\in {\left[0,1\right]}^2$ are dimensionless as well as weights. The tessellation is completely determined by the parameter vectors $x=\left(x_1,\cdots ,x_n\right)$, $y=\left(y_1,\cdots ,y_n\right)$, $w=\left(w_1,\cdots ,w_n\right)$. The dimensionless tessellation is scaled by a physical length L₀ 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: $$\begin{array}{lll}{C}_j\hfill & =\hfill & \{\left(\begin{array}{c}x\\ y\end{array}\right)\in {ℝ}^2,\forall k\in \left\{1,\cdots ,n\right\},{\Vert \begin{array}{c}x-x_j\\ y-y_j\end{array}\Vert }^2\hfill \\ \hfill & \hfill & -w_j\le {\Vert \begin{array}{c}x-x_k\\ y-y_k\end{array}\Vert }^2-w_k\}.\hfill \end{array}$$(1)

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: $${\sum }_{j=1}^n{w}_j=1.$$(2)

Thus, weights w lie in an affine hyperplane of dimension n − 1 and denoted by $P_a^{(n-1)}=\left\{w\in {ℝ}^n,{\sum }_{j=1}^n{w}_j=1\right\}$, whose support is the hyperplane denoted by $P^{(n-1)}=\left\{w\in {ℝ}^n,{\sum }_{j=1}^n{w}_j=0\right\}$. 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 $\theta =\left({\theta }_1,\cdots ,{\theta }_n\right)$. Since the crystal lattice is plane hexagonal $\forall j\in \left\{1,\dots ,n\right\},{\theta }_j\in \left[0,\frac{\pi }{3}\right]$. The parameter set ${\mathcal{P}}_{\text{OT}}$ defining the OT (see Fig. 3) reads: $${\mathcal{P}}_{\text{OT}}=\left\{\alpha =(x,y,w,\theta )\in {ℝ}^n\times {ℝ}^n\times P_a^{(n-1)}\times {\left[0,\frac{\pi }{3}\right]}^n\right\}\text{}\text{}.$$(3)

A total energy $\mathcal{E}(T,\alpha )$ 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 ${\mathcal{P}}_{\text{OT}}$ by deriving an evolution law of the form: $$\stackrel{.}{\alpha }=-M(T,\alpha ).\frac{\partial \mathcal{E}(T,\alpha )}{\partial \alpha }.$$(4)

where M(T, α) is a second order tensor to be defined by introducing a dissipative mechanism at microscale and $\partial \mathcal{E}(T,\alpha )/\partial \alpha$ is the driving force.

This paper is limited to the evolution of weights (i.e., seeds and orientations are fixed parameters: $\stackrel{.}{x}=\stackrel{.}{y}=0$ and $\stackrel{.}{\theta }=0$), that is to say that the generalized vector space of speed ${\mathcal{V}}^{*}$ reads: $${\mathcal{V}}^{*}=\left\{{\stackrel{.}{w}}^{*},{\sum }_{j=1}^n{\dot{w}}_j^{*}=0\right\}.$$(5)

Fig. 2 Seed and weight in tessellation.

Fig. 3 Seed and weight in tessellation.

3 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 : $$I_{GB}=\left\{(i,j)\in {\left\{1,\cdots ,n\right\}}^2,j\ge i,C_i\cap C_j\ne \varnothing \right\}.$$(6)

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 (i, j) ∈ I_GB is denoted E_ij. $$E_{ij}=E(T,\Delta {\theta }_{ij},{\phi }_{ij})$$(7)

where Δθ_ij is the crystal lattice disorientation between grains i and j defined by: $$\Delta {\theta }_{ij}=\left|{\theta }_j-{\theta }_i\right|\in \left[0,\frac{\pi }{3}\right].$$(8)

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]: $$E(T,\Delta \theta ,\phi )=E_0(T,\phi )\Delta \theta \left[a(\phi )-\ln \left(\Delta \theta \right)\right].$$(9)

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: $$\{\begin{array}{l}E(T,\Delta \theta ,\phi )=E_M(T,\phi )\sin \left(\frac{\pi }{2}\frac{\Delta \theta -\Delta {\theta }_m}{\Delta {\theta }_M-\Delta {\theta }_m}\right)\hfill \\ \times \left[1-a(\phi )\ln \left(\sin \left(\frac{\pi }{2}\frac{\Delta \theta -\Delta {\theta }_m}{\Delta {\theta }_M-\Delta {\theta }_m}\right)\right)\right]\hfill \\ \left(\Delta {\theta }_m\le \Delta \theta \le \Delta {\theta }_M\right).\hfill \end{array}$$(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: $$E(T,\Delta \theta )=\frac{G(T)}{G(0)}\gamma (\Delta \theta )$$(11)

where G(T) is the shear coefficient. In this paper, data for pure iron [22] are used to calibrate G(T): $$G(T)=a_{G}T+b_G$$(12)

where b_G ≈ 88134 MPa and a_G ≈−24 MPa.K⁻¹. The temperature dependence G(T)∕G(0) is a simple choice similar to [23]. In addition: $$\{\begin{array}{l}\gamma (\Delta \theta )={\gamma }_{\frac{\pi }{6}}\sin \left(3\Delta \theta \right)\left[1-a_1\ln \left(\sin \left(3\Delta \theta \right)\right)\right]\hfill \\ \left(0\le \Delta \theta \le \frac{\pi }{6}\right)\hfill \\ \gamma (\Delta \theta )={\gamma }_{\frac{\pi }{6}}\sin \left(\pi -3\Delta \theta \right)\left[1-a_2\ln \left(\sin \left(\pi -3\Delta \theta \right)\right)\right]\hfill \\ \left(\frac{\pi }{6}\le \Delta \theta \le \frac{\pi }{3}\right)\hfill \end{array}$$(13)

where a₁ = a₂ = 1 and ${\gamma }_{\frac{\pi }{6}}=1$ J.m⁻², 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: $$D_{ij}^{*}=D(T,\Delta {\theta }_{ij},v_{ij}^{*})$$(14)

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}^{*}$: $$D(T,\Delta {\theta }_{ij},v_{ij}^{*})={\tau }_{c}X(\Delta {\theta }_{ij})\left|v_{ij}^{*}\right|$$(15)

where τ_c is the critical shear stress and X : Δθ↦X(Δθ) is the following function: $$X(\Delta \theta )=\frac{6}{\pi }\left(\frac{\pi }{3}+2\sqrt{3}\ln \left(\frac{\sqrt{3}}{2}\right)\right)\min \left\{\Delta \theta ,\frac{\pi }{3}-\Delta \theta \right\}.$$(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 $\left|v_{ij}^{*}\right|$: $${\tau }_c=\frac{\left|v_{ij}^{*}\right|}{m(T)}$$(17)

where the function m(T) is homogenous to a grain mobility (m⁴.J⁻¹.s⁻¹). Hence the dissipated power per unit area: $$D(T,\Delta {\theta }_{ij},v_{ij}^{*})=\frac{X(\Delta {\theta }_{ij})}{m(T)}{\left[v_{ij}^{*}\right]}^2.$$(18)

It should be noted from (10) and (16) that: $$\{\begin{array}{l}\gamma (\Delta \theta )\underset{\Delta \theta \to 0}{~}A\Delta \theta \ln \left(\Delta \theta \right)\hfill \\ X(\Delta \theta )\underset{\Delta \theta \to 0}{~}B\Delta \theta .\hfill \end{array}$$(19)

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 \left(\Delta \theta \right)$).

Fig. 4 Notations.

Fig. 5 Grain boundary energy per unit area at T = 0 K.

4 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: $$\{\begin{array}{l}\mathcal{E}(T,\alpha )=L_0{\sum }_{(i,j)\in I_{GB}}{l}_{ij}{E}_{ij}\hfill \\ \mathcal{D}(T,\alpha ,v^{*})=L_0{\sum }_{(i,j)\in I_{GB}}{l}_{ij}{D}_{ij}^{*}\hfill \end{array}$$(20)

where l_ij the joint length and $v^{*}={\left(v_{ij}^{*}\right)}_{(i,j)\in 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 $n_{GB}=\text{card}\left[I_{GB}\right]$. There exists a bijection between I_GB and $\left\{1,\cdots ,n_{GB}\right\}$ that definesthe numbering of grain boundaries. Thus, for each grain boundary denoted by (i, j) ∈ I_GB there exists a unique $p_{ij}\in \left\{1,\cdots ,n_{GB}\right\}$. Therefore the total dissipated power per unit area reads: $$\mathcal{D}(T,\alpha ,v^{*})=L_0{v}^{*}.\frac{\chi (\alpha )}{m(T)}.v^{*}$$(21)

where χ(α) is a dimensionless diagonal second order tensor of size n_GB × n_GB of diagonal: $$\forall p_{ij}\in \left\{1,\cdots ,n_{GB}\right\},{\chi }_{p_{ij}}(\alpha )=l_{ij}X(\Delta {\theta }_{ij}).$$(22)

As seeds and crystalline orientations are fixed (i.e., $\stackrel{.}{x}=\stackrel{.}{y}=0$ and $\stackrel{.}{\theta }=0$) and considering any virtual variation of weights ${\stackrel{.}{w}}^{*}$, it is straightforward to demonstrate that: $$v_{ij}^{*}=L_0\frac{{\dot{w}}_i^{*}-{\dot{w}}_j^{*}}{2d_{ij}}$$(23)

where d_ij is the dimensionless distance between seeds: $$d_{ij}=\Vert \left(\begin{array}{c}{x}_i\\ y_i\end{array}\right)-\left(\begin{array}{c}{x}_j\\ y_j\end{array}\right)\Vert .$$(24)

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: $$v^{*}=L_{0}K(\alpha ).{\stackrel{.}{w}}^{*}$$(25)

where K(α) is a dimensionless second order tensor of size n_GB × n, which can be evaluated analytically: $$K(\alpha )=p_{ij}\left(\begin{array}{cccccccccccc}& & & & i& & & & j& & & \\ & ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ & 0& \cdots & 0& \frac{1}{2d_{ij}}& 0& \cdots & 0& -\frac{1}{2d_{ij}}& 0& \cdots & 0\\ & ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\end{array}\right)$$(26)

The constraint (2) leads to: $${\sum }_{j=1}^n{\dot{w}}_j^{*}=0$$(27)

which defines a hyperplane of dimension n − 1 denoted by P⁽ⁿ⁻¹⁾, whose normal vector is denoted by $1=(1,\cdots ,1)\in {ℝ}^n$. It is clear that the kernel of K(α) is $\text{ker}\left(K(\alpha )\right)=1ℝ$. The total dissipated power per unit depth reads: $$\mathcal{D}(T,\alpha ,v^{*})=L_0^3{\stackrel{.}{w}}^{*}.\frac{R(\alpha )}{m(T)}.{\stackrel{.}{w}}^{*}$$(28)

where $R(\alpha )=K{(\alpha )}^T.\chi (\alpha ).K(\alpha )$ is a dimensionless second order tensor of size n × n.

The “real” weight variation $\stackrel{.}{w}$ 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): $$\stackrel{.}{\mathcal{E}}+\mathcal{D}=0.$$(29)

Hence: $$\stackrel{.}{w}.\left(\frac{\partial \mathcal{E}}{\partial w}+L_0^3\frac{R(\alpha )}{m(T)}.\stackrel{.}{w}\right)=0.$$(30)

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 constraint (30): $$\stackrel{.}{w}=\{\begin{array}{l}\underset{\stackrel{.}{w}}{\text{argmax}}\left[L_0^3\stackrel{.}{w}.\frac{R(\alpha )}{m(T)}.\stackrel{.}{w}\right]\hfill \\ \text{subjectedto}.\stackrel{.}{w}.\left(\frac{\partial \mathcal{E}}{\partial w}+L_0^3\frac{R(\alpha )}{m(T)}.\stackrel{.}{w}\right)=0.\hfill \end{array}$$(31)

Consider μ the Lagrangian multiplier and L the Lagrangian: $$L(\stackrel{.}{w},\mu )=L_0^3\stackrel{.}{w}.\frac{R(\alpha )}{m(T)}.\stackrel{.}{w}+\mu \stackrel{.}{w}.\left(\frac{\partial \mathcal{E}}{\partial w}+L_0^3\frac{R(\alpha )}{m(T)}.\stackrel{.}{w}\right).$$(32)

The optimality condition reads: $$\{\begin{array}{l}\frac{\partial L}{\partial \stackrel{.}{w}}(\stackrel{.}{w},\mu )=2L_0^3\frac{R(\alpha )}{m(T)}.\stackrel{.}{w}+\mu \left(\frac{\partial \mathcal{E}}{\partial w}+2L_0^3\frac{R(\alpha )}{m(T)}.\stackrel{.}{w}\right)\hfill \\ =0\hfill \\ \frac{\partial L}{\partial \mu }(\stackrel{.}{w},\mu )=\stackrel{.}{w}.\left(\frac{\partial \mathcal{E}}{\partial w}+L_0^3\frac{R(\alpha )}{m(T)}.\stackrel{.}{w}\right)=0.\hfill \end{array}$$(33)

Hence μ = −2 and the evolution law reads: $$R(\alpha ).\stackrel{.}{w}=-\frac{m(T)}{L_0^3}\frac{\partial \mathcal{E}(T,\alpha )}{\partial w}.$$(34)

The evolution law (34) should be inverted. However R(α) is not invertible because of the constrain (27). The kernel of R(α) is $\text{ker}\left(R(\alpha )\right)=1ℝ$. Thus, there are n − 1 strictly positive singular values of R(α) denoted by (λ₁, ⋯, λ_n−1) and there exist two orthogonal tensors U and V such as: $$R(\alpha )=U.\left(\begin{array}{cc}\begin{array}{ccc}{\lambda }_1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\lambda }_{n-1}\end{array}& 0\\ 0& 0\end{array}\right).V$$(35)

The singular value decomposition is computed in Scilab [27] and gives (λ₁, ⋯, λ_n−1), U and V. The Moore-Penrose pseudo-inverse matrix is introduced as follows: $${\Sigma }^{\text{†}}=\left(\begin{array}{cc}\begin{array}{ccc}\frac{1}{{\lambda }_1}& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & \frac{1}{{\lambda }_{n-1}}\end{array}& 0\\ 0& 0\end{array}\right).$$(36)

Hence the evolution law: $$\stackrel{.}{w}=-\frac{m(T)}{L_0^3}M(\alpha ).\frac{\partial \mathcal{E}(T,\alpha )}{\partial w}$$(37)

where: $$M(\alpha )=\left(V^T.{\Sigma }^{\text{†}}.U^T\right).$$(38)

It should be noted that a scale effect is demonstrated in the evolution law (37) as the physical dimension L₀ arises.

5 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: $$\mathcal{E}(T,\alpha )=L_0\frac{G(T)}{G(0)}{\sum }_{(i,j)\in I_{GB}}{l}_{ij}\gamma (\Delta {\theta }_{ij}).$$(39)

Hence: $$\frac{\partial \mathcal{E}(T,\alpha )}{\partial w}=L_0\frac{G(T)}{G(0)}{\sum }_{(i,j)\in I_{GB}}\frac{\partial l_{ij}}{\partial w}\gamma (\Delta {\theta }_{ij}).$$(40)

Consider the set of triple junctions: $$I_3=\left\{(i,j,k)\in {\left\{1,\cdots ,n\right\}}^3,C_i\cap C_j\cap C_k\ne \varnothing \right\}.$$(41)

There is no condition such as k ≥ j ≥ i, thus if (i, j, k) ∈ I₃ then all permutations are also in I₃. Consider the third order tensor δ: $${\tilde{\delta }}_{ijk}=\{\begin{array}{lll}1\hfill & \text{if}\hfill & (i,j,k)\in I_3\hfill \\ 0\hfill & \text{if}\hfill & (i,j,k)\notin I_3\hfill \end{array}$$(42)

Thus, the third order tensor β representing angles at the triple junctions as shown in Figure 6 are defined as follows: $${\beta }_{ijk}=\{\begin{array}{lll}{\beta }_{ijk}\hfill & \text{if}\hfill & (i,j,k)\in I_3\hfill \\ \frac{\pi }{2}\hfill & \text{if}\hfill & (i,j,k)\notin I_3\hfill \end{array}$$(43)

where the following symmetry rule holds β_ijk = β_kji.

Consider a virtual weight variation ${\stackrel{.}{w}}^{*}\in P^{(n-1)}$ (i.e., verifying the constrain (27)). One obtains: $$\begin{array}{ll}{\dot{l}}_{ij}\hfill & =\frac{\tan \left({\beta }_{ijk}-\frac{\pi }{2}\right)+\tan \left({\beta }_{ijm}-\frac{\pi }{2}\right)}{2d_{ij}}{\dot{w}}_i^{*}\hfill \\ \hfill & +\frac{\tan \left({\beta }_{jik}-\frac{\pi }{2}\right)+\tan \left({\beta }_{jim}-\frac{\pi }{2}\right)}{2d_{ij}}{\dot{w}}_j^{*}\hfill \\ \hfill & -\frac{1}{2d_{jk}\cos \left({\beta }_{ijk}-\frac{\pi }{2}\right)}{\dot{w}}_k^{*}-\frac{1}{2d_{jm}\cos \left({\beta }_{ijm}-\frac{\pi }{2}\right)}{\dot{w}}_m^{*}.\hfill \end{array}$$(44)

Moreover: $${\dot{l}}_{ij}=\left(\frac{\partial l_{ij}}{\partial w}\right).{\stackrel{.}{w}}^{*}.$$(45)

Hence by combining (44) and (45) and $\partial l_{ij}/\partial w={\left(\partial l_{ij}/\partial w_q\right)}_{q\in \left\{1,\cdots ,n\right\}}$: $$\begin{array}{lll}\frac{\partial l_{ij}}{\partial w_q}\hfill & =\hfill & {\sum }_{p=1}^n{\tilde{\delta }}_{ijp}(\frac{\tan \left({\beta }_{ijp}-\frac{\pi }{2}\right)}{2d_{ij}}{\delta }_{iq}+\frac{\tan \left({\beta }_{jip}-\frac{\pi }{2}\right)}{2d_{ij}}{\delta }_{jq}\hfill \\ \hfill & \hfill & -\frac{{\delta }_{pq}}{2d_{jp}\cos \left({\beta }_{ijp}-\frac{\pi }{2}\right)})\hfill \end{array}$$(46)

where δ is the Kronecker symbol.

Fig. 6 Triple junctions.

6 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: $$\tilde{F}=\frac{G(T)}{G(0)}\gamma (\Delta \theta )\left(\frac{1}{R_I}+\frac{1}{R_{II}}\right)$$(47)

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: $$\tilde{v}=\tilde{m}(T,\Delta \theta )\tilde{F}$$(48)

where $\tilde{v}$ is the inward normal speed of the considered grain boundary and $\tilde{m}(T,\Delta \theta )$ the classical grain boundary mobility for the curvature driven model (where the symbol $\widetilde{.}$ 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 $\tilde{m}(T,\Delta \theta )$. Thus for the curvature driven model the inward normal speed reads: $$\tilde{v}=\tilde{m}(T,\Delta \theta )\frac{G(T)}{G(0)}\frac{\gamma (\Delta \theta )}{R}.$$(49)

The proposed approach is applied to the hexagonal situation. The partial derivative (46) reads: $$\begin{array}{l}\forall j\in \left\{2,\cdots ,7\right\}\hfill \\ \frac{\partial l_{1j}}{\partial w_1}=\frac{\partial l_{1j}}{\partial w_j}=-\frac{\partial l_{1j}}{\partial w_{j-1}}=-\frac{\partial l_{1j}}{\partial w_{j+1}}=\frac{1}{2\sqrt{3}R}.\hfill \end{array}$$(50)

Hence: $$\begin{array}{ll}\hfill & \frac{\partial l_{12}}{\partial w}=\frac{1}{2\sqrt{3}R}{\left(1,1,-1,0,0,0,-1\right)}^T\hfill \\ \hfill & \frac{\partial l_{13}}{\partial w}=\frac{1}{2\sqrt{3}R}{\left(1,-1,1,-1,0,0,0\right)}^T\hfill \\ \hfill & \frac{\partial l_{14}}{\partial w}=\frac{1}{2\sqrt{3}R}{\left(1,0,-1,1,-1,0,0\right)}^T\hfill \\ \hfill & \frac{\partial l_{15}}{\partial w}=\frac{1}{2\sqrt{3}R}{\left(1,0,0,-1,1,-1,0\right)}^T\hfill \\ \hfill & \frac{\partial l_{16}}{\partial w}=\frac{1}{2\sqrt{3}R}{\left(1,0,0,0,-1,1,-1\right)}^T\hfill \\ \hfill & \frac{\partial l_{17}}{\partial w}=\frac{1}{2\sqrt{3}R}{\left(1,-1,0,0,0,-1,1\right)}^T.\hfill \end{array}$$(51)

The driving force (40) reads: $$\frac{\partial \mathcal{E}(T,\alpha )}{\partial w}=\frac{G(T)}{G(0)}\frac{\gamma (\Delta \theta )}{2\sqrt{3}R}{\left(6,-1,-1,-1,-1,-1,-1\right)}^T.$$(52)

The second order tensor K introduced in (26) reads for this configuration: $$\begin{array}{l}\begin{array}{cccccccc}& 1& 2& 3& 4& 5& 6& 7\end{array}\\ K=\begin{array}{c}{p}_{12}\\ p_{13}\\ p_{14}\\ p_{15}\\ p_{16}\\ p_{17}\end{array}\left(\begin{array}{cccccccc}& -1& 1& 0& 0& 0& 0& 0\\ & -1& 0& 1& 0& 0& 0& 0\\ & -1& 0& 0& 1& 0& 0& 0\\ & -1& 0& 0& 0& 1& 0& 0\\ & -1& 0& 0& 0& 0& 1& 0\\ & -1& 0& 0& 0& 0& 0& 1\end{array}\right)\frac{1}{4R}.\end{array}$$(53)

The second order tensor χ introduced in (22) reads: $$\chi =\frac{2R}{\sqrt{3}}X(\Delta \theta )I$$(54)

where I is the identity of size 6 × 6. The second order tensor R introduced in (28) reads: $$\begin{array}{lll}R\hfill & =\hfill & K^T.\chi .K\hfill \\ \hfill & =\hfill & \frac{X(\Delta \theta )}{8\sqrt{3}R}\left(\begin{array}{ccccccc}6& -1& -1& -1& -1& -1& -1\\ -1& 1& 0& 0& 0& 0& 0\\ -1& 0& 1& 0& 0& 0& 0\\ -1& 0& 0& 1& 0& 0& 0\\ -1& 0& 0& 0& 1& 0& 0\\ -1& 0& 0& 0& 0& 1& 0\\ -1& 0& 0& 0& 0& 0& 1\end{array}\right).\hfill \end{array}$$(55)

Thus, the second order tensor M introduced in (38) reads: $$M=m(T)\frac{8\sqrt{3}R}{X(\Delta \theta )}\frac{1}{49}\left(\begin{array}{ccccccc}6& -1& -1& -1& -1& -1& -1\\ -1& 41& -8& -8& -8& -8& -8\\ -1& -8& 41& -8& -8& -8& -8\\ -1& -8& -8& 41& -8& -8& -8\\ -1& -8& -8& -8& 41& -8& -8\\ -1& -8& -8& -8& -8& 41& -8\\ -1& -8& -8& -8& -8& -8& 41\end{array}\right).$$(56)

Hence the evolution law (37) for the hexagonal situation: $$\stackrel{.}{w}=-m(T)\frac{4}{7}\frac{G(T)}{G(0)}\frac{\gamma (\Delta \theta )}{X(\Delta \theta )}{\left(6,-1,-1,-1,-1,-1,-1\right)}^T.$$(57)

Thus, by using (23) one obtains the inward normal speed according to the proposed model: $$v=\frac{m(T)}{R}\frac{G(T)}{G(0)}\frac{\gamma (\Delta \theta )}{X(\Delta \theta )}.$$(58)

The comparison between the inward normal speed from the curvature driven model (49) and from the proposed approach (58) reads: $$\tilde{m}(T,\Delta \theta )=\frac{m(T)}{X(\Delta \theta )}.$$(59)

For the sake of simplicity it is assumed that $\tilde{m}(T,\Delta \theta )$ depends on disorientation in 1∕X(Δθ) even though more realistic mobilities could be obtained by considering m(T, Δθ) instead of m(T) in (17). $$\tilde{m}(T,\Delta \theta )X(\Delta \theta )\approx \tilde{\tilde{m}}(T).$$(60)

To estimate m(T) one can use the pure iron grain boundary mobility determined by [28]: $$m(T)=\beta \frac{D_{\text{Fe}}(T)V_m\delta }{b^2\text{R}T}$$(61)

where the molar volume of fcc-Fe is V_m = 7.09 × 10⁻⁶ m³, the boundary thickness is δ = 1 nm, the burgers vector is b = 2.48 × 10⁻¹⁰ m, the gas constant is R = 8.3144621 J.mol⁻¹.K⁻¹. In addition, D_Fe(T) is the diffusivity of Fe atoms along the grain boundary (for fcc) and follows an Arrhenius law: $$D_{\text{Fe}}(T)=D_{\text{Fe}}^0\exp \left(-\frac{Q}{\text{R}T}\right)$$(62)

where the pre-factor is $D_{\text{Fe}}^0=1.5\times {10}^{-4}$ m².s⁻¹ the activation energy is Q = 148 kJ.mol⁻¹.

Fig. 7 Validation by comparison to curvature driven spherical grain shrinkage.

7 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: $$N=\Vert M(T,\alpha ).\frac{\partial \mathcal{E}(T,\alpha )}{\partial w}\Vert$$(63) $$\text{d}t=\frac{{\delta }_w}{N}$$(64)

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. 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.)

Grain growth at 800 °C is presented for three OT (denoted by OT1, OT2, OT3) in Figures 10, 11 and 12 (orientations are in degrees). Simulations show grain coarsening and are consistent with experimental observations [30] on an other material (tantalum). Grain growth is heterogeneous depending on the local disorientation and 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.

Fig. 8 Algorithm.

Fig. 9 Disorientation distributions of three OT at 800 °C.

Fig. 10 Grain growth of the OT1 at 800 °C during approximately 4 h.

Fig. 11 Grain growth of the OT2 at 800 °C during approximately 4 h.

Fig. 12 Grain growth of the OT3 at 800 °C during approximately 4 h.

Fig. 13 Grain growth of the OT1 at 700 °C during approximately 4 h.

8 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.

References

All Figures

	Fig. 1 FCC and slip directions in <1, 1, 1> plane.
In the text

	Fig. 2 Seed and weight in tessellation.
In the text

	Fig. 3 Seed and weight in tessellation.
In the text

	Fig. 4 Notations.
In the text

	Fig. 5 Grain boundary energy per unit area at T = 0 K.
In the text

	Fig. 6 Triple junctions.
In the text

	Fig. 7 Validation by comparison to curvature driven spherical grain shrinkage.
In the text

	Fig. 8 Algorithm.
In the text

	Fig. 9 Disorientation distributions of three OT at 800 °C.
In the text

	Fig. 10 Grain growth of the OT1 at 800 °C during approximately 4 h.
In the text

	Fig. 11 Grain growth of the OT2 at 800 °C during approximately 4 h.
In the text

	Fig. 12 Grain growth of the OT3 at 800 °C during approximately 4 h.
In the text

	Fig. 13 Grain growth of the OT1 at 700 °C during approximately 4 h.
In the text