Phase-field modeling of stored-energy-driven grain growth with intra-granular variation in dislocation density

To account for the variation in intra-granular dislocation densities observed in the HEDM experiments, we introduce an individual dislocation density field, ${\rho _i}$, for each grain $i$, which stores not only the initial dislocation density value for that grain but also what the dislocation density value would be if the grain boundary of the grain were to reach that position. While the exact mechanism by which the reduction in the dislocation density occurs has yet to be investigated, there are multiple mechanisms that would lead to an elimination of dislocations as a grain boundary migrates, including annihilation of dislocation pairs and absorption of dislocations onto grain boundaries [42, 43]. The enhancement of the rate by which these processes occur (vs. recovery in the bulk, which appears to be slow on the time scale of the experiment) is expected from the fact that grain boundary migration requires atomic rearrangements and consequently changes the strain fields. Below, we derive a model for dislocation density evolution that accounts for such a process, which is used to define ${\rho _i}$.

Mathematically, it is reasonable to assume that the dislocation density is reduced at a rate, $\frac{{{\text{d}}\rho }}{{{\text{d}}x}}$, that is proportional to how many dislocations per unit volume exist at the grain boundary position as the grain boundary migrates:

$$\begin{align}\frac{{{\text{d}}\rho }}{{{\text{d}}x}} = - \frac{\rho }{{{l_e}}},\end{align} \tag{ 1 }$$

where $x$ is the distance traversed by the grain boundary and the inverse of ${l_e}$ provides the rate constant for the decay. This differential equation can be solved analytically, yielding

$$\begin{align}\rho = {\rho _{{\text{gb}}}}\exp \left( { - \frac{x}{{{l_e}}}} \right),\end{align} \tag{ 2 }$$

where ${\rho _{{\text{gb}}}}$ is the dislocation density just behind the initial grain boundary and the constant ${l_e}$ can now be identified as the length scale at which the dislocation density decays by a factor of $e$ in the wake of the migrating grain boundary. We assume that grain boundaries migrate in their normal direction, which is a simplifying assumption, and therefore the value of $x$ can be determined from a signed distance function using the level-set method (detailed in section S1 of the supplementary material).

We now describe how this model is implemented. As in figure 2, we denote the grain with index $i$ as ${G_i}$. For each grain ${G_i}$, we employ a time- and position-dependent order parameter field ${\eta _i}\left( {\boldsymbol{r},t} \right)$ to track its evolution. The value of ${\eta _i}\left( {\boldsymbol{r},t} \right)$ is 1 within the grain, 0 outside the grain, and varies smoothly from 1 to 0 across its interface. In this work, we omit the effect of recovery, and therefore the dislocation density at a certain position, $r$, retains its initial value until the grain boundary from another grain reaches this point, at which time a reduced dislocation density is assigned, according to equation (2). Since we do not allow dislocation density to evolve otherwise, this model introduces discontinuities in the gradient of dislocation density, which persist; if we consider recovery, the magnitude of such discontinuities would gradually reduce and eventually vanish. We precompute the reduced dislocation density for each grain in regions not initially occupied by the grain and store these values into its individual dislocation density field ${\rho _i}\left( \boldsymbol{r} \right)$, along with its initial dislocation density values within the grain. Specifically, ${\rho _i}$ is defined as

$$\begin{align}{\rho _i}\left( \boldsymbol{r} \right) =\left \{ \begin{array}{*{20}{l}} {{\rho _{0,i}}\left( \boldsymbol{r} \right),\kern105.4pt{\text{within}}\;{\text{Grain}}{\text{ }}{G_i}} \\ {\max \left( {{\rho _{{\text{base}}}},{\rho _{\text{gb,i}}}\exp \left( { - \frac{{\left| {{d_i}\left( \boldsymbol{r} \right)} \right|}}{{{l_e}}}} \right)} \right),\;\;{\text{otherwise}}} \end{array}\right..\end{align} \tag{ 3 }$$

Here, ${\rho _{0,i}}\left( \boldsymbol{r} \right)$ is the initial dislocation density within grain ${G_i}$, ${\rho _{{\text{gb}},i}}$ is the dislocation density just behind ${G_i}$'s initial boundary, and ${d_i}\left( \boldsymbol{r} \right)$ is the signed distance between the point $r$ and the nearest point to the ${G_i}$'s initial boundary; note that the sign of ${d_i}\left( \boldsymbol{r} \right)$ is taken to be positive within the grain and negative outside, and therefore ${d_i}\left( \boldsymbol{r} \right) < 0$ in equation (3). The term ${\rho _{{\text{base}}}}$ in equation (3) denotes a baseline dislocation density, i.e., the lower bound for the dislocation density measured in the experiment, that remains even after a long period of annealing. In general, ${\rho _{{\text{gb}},i}}$ is a spatially varying quantity, but we here assume it to be a constant for each grain, as discussed in section 3.4.2. We found that simulations with ${l_e}$ = 20 $\mu {\text{m}}$ yield results that closely resemble the experimental observations, and thus we employ this value throughout this work. We estimate the baseline dislocation density ${\rho _{{\text{base}}}}$ to be $3 \times {10^{12}}$ m⁻² by taking the average of the dislocation densities at state S₂ for the regions outside grain ${G_1}$'s initial position (indicated in figure 2(c)). While this may slightly overestimate the base value because it may include some enhanced dislocation density near the initial grain boundaries, we deemed this assumption to be reasonable given the statistical variation and uncertainties in the measurements [44]. In addition, the determined value of ${\rho _{{\text{base}}}}$ is of the same order of magnitude as dislocation densities in well-annealed metals ($2 \times {10^{11}}$ m⁻² to $3 \times {10^{12}}$ m⁻² [45]), which supports the validity of the assumption. Further details are provided in section 3.4.2.

It is necessary to note that the values of ${\rho _i}\left( \boldsymbol{r} \right)$ outside the boundaries of the initial grains should not be considered as physical values, other than the fact that it sets the value of the dislocation density if the grain boundary of ${G_i}$ sweeps past the point. Rather, at any point in time, the actual dislocation density is given by constructing it from $\left\{ {{\rho _i}\left( \boldsymbol{r} \right)} \right\}$, along with the set of order parameters $\left\{ {{\eta _i}\left( {\boldsymbol{r},{\text{ }}t} \right)} \right\}$, using the following interpolation function:

$$\begin{align}\rho \left( {\boldsymbol{r},t} \right) = \frac{{\mathop \sum \nolimits_{i = 1}^N \eta _i^2\left( {\boldsymbol{r},t} \right){\rho _i}\left( \boldsymbol{r} \right)}}{{\mathop \sum \nolimits_{i = 1}^N \eta _i^2\left( {\boldsymbol{r},t} \right)}}.\end{align} \tag{ 4 }$$

This formula is identical to those used in [16, 25, 37, 46], except that in these previous works, ${\rho _i}$ was a constant associated with the $i^{\text{th}}$ grain, and therefore these models could not treat intra-granular inhomogeneity nor a decay in the dislocation density as the grain boundaries migrate.

To facilitate our understanding of dislocation density evolution based on the above formulation, we provide schematics in figure 4 showing the change in the dislocation density field as grain ${G_i}$ consumes its surrounding grains. The grains that are adjacent to ${G_i}$ are omitted in figure 4(a) for clarity. Grain ${G_i}$, indicated in blue, initially possesses dislocation density ${\rho _{{\text{gb}},i}}$ at the position of its grain boundaries. As it progressively consumes its surrounding grains with higher dislocation density (denoted by ${\rho _{{\text{surr}}}}$ in figure 4) and reaches the point $A$ at time ${t_1}$, the dislocation density reduces to a value given by ${\rho _{{\text{gb}},i}}\exp \left( { - \frac{{\left| {{d_i}\left( A \right)} \right|}}{{{l_e}}}} \right)$. As the grain continues to grow to point $B$ at time ${t_2}$, where by ${\rho _{{\text{gb}},i}}\exp \left( { - \frac{{\left| {{d_i}\left( B \right)} \right|}}{{{l_e}}}} \right) = {\rho _{{\text{base}}}}$, the dislocation density value behind the migrating grain boundary reduces to the baseline dislocation density ${\rho _{{\text{base}}}}$, which is the lower bound of the dislocation density. Beyond this point, the dislocation density behind the grain boundary remains at ${\rho _{{\text{base}}}}$; for example, when the grain boundary reaches point $C$ at time ${t_3}$, the dislocation density at point $C$ reduces to ${\rho _{{\text{base}}}}$.

Zoom In Zoom Out Reset image size

A total free energy functional that encodes the thermodynamics of the polycrystalline material is given by the volume integral of the sum of three energy density terms:

$$\begin{align}\mathcal{F}\left( {\left\{ {{\eta _i}} \right\}} \right) = \int\limits_V \left( {{f_{{\text{bulk}}}} + {f_{{\text{gradient}}}} + {f_{{\text{stored}}}}} \right){\text{d}}V.\end{align} \tag{ 5 }$$

We utilize the functional form proposed by Moelans et al [47] to calculate the bulk energy density, ${f_{{\text{bulk}}}}$:

$$\begin{align}{f_{{\text{bulk}}}} = W\left[ {\mathop \sum \limits_{i = 1}^N \left( { - \frac{1}{2}\eta _i^2 + \frac{1}{4}\eta _i^4} \right) + \alpha \mathop \sum \limits_{i = 1}^N \mathop \sum \limits_{j > 1}^N \left( {\eta _i^2\eta _j^2} \right) + \frac{1}{4}} \right],\end{align} \tag{ 6 }$$

where $W$ is the bulk energy density coefficient, and $\alpha = 1.5$ is set to ensure a symmetrical order parameter profile across an interface, as discussed in [47]. The gradient energy density [20, 48], ${f_{{\text{gradient}}}}$, is employed to penalize a sharp interface,

$$\begin{align}{f_{{\text{gradient}}}} = \frac{\kappa }{2}\mathop \sum \limits_{i = 1}^N {\left| {\nabla {\eta _i}} \right|^2},\end{align} \tag{ 7 }$$

where $\kappa$ is the gradient energy coefficient. The values of $W$ and ${{\kappa }}$ are determined based on the grain boundary energy, ${\gamma _{{\text{gb}}}}$, a material parameter, and the grain boundary width, ${l_{{\text{gb}}}}$, a model parameter, using the formula derived by Moelans et al [47]:

$$\begin{align}W& = \frac{{6{\gamma _{{\text{gb}}}}}}{{{l_{{\text{gb}}}}}},\end{align} \tag{ 8 }$$

$$\begin{align}\kappa& = \frac{3}{4}{\gamma _{{\text{gb}}}}{l_{{\text{gb}}}}.\end{align} \tag{ 9 }$$

The grain boundary energy for Cu–Zn [49], ${\gamma _{{\text{gb}}}} \approx 0.595$ J m⁻², is employed here for Cu–Al–Mn, due to a lack of data on the latter system. An approximate form [17], which neglects the energy of the dislocation core and assumes isotropic elasticity, is utilized to describe the stored energy density, ${f_{{\text{stored}}}}$, due to the presence of the dislocation density, $\rho \left( {\boldsymbol{r},t} \right)$, as

$$\begin{align}{f_{{\text{stored}}}} = \frac{1}{2}G{b^2}\rho \left( {\boldsymbol{r},t} \right),\end{align} \tag{ 10 }$$

where $G = 28.4$ GPa is the shear modulus for the $\beta$ phase of Cu–Al–Mn alloy at 650 $^\circ$C, calculated by [14] based on [50]. The magnitude of the Burger's vector for the FCC copper [51], $b = 0.255$ nm, is used here as an estimate value for the copper-rich Cu–Al–Mn alloy. The evolution of the order parameter ${\eta _i}$ is obtained via steepest descent, using Allen–Cahn dynamics [52]:

$$\begin{align}\frac{{\partial {\eta _i}}}{{\partial t}} = - L\left( {\frac{{\delta F}}{{\delta {\eta _i}}}} \right),\end{align} \tag{ 11 }$$

where $L$ is the mobility parameter. Substituting equations (5)–(7) and (10) into equation (11) yields

$$\begin{align}\frac{{\partial {\eta _i}}}{{\partial t}} = - L\left[ {W\left( { - {\eta _i} + \eta _i^3 + 2\alpha {\eta _i}\mathop \sum \limits_{\begin{array}{*{20}{c}} {j = 1} \\ {j \ne i} \end{array}}^N \eta _j^2{\text{ }}} \right) - \kappa {\nabla ^2}{\eta _i} + G{b^2}\frac{{{\eta _i}}}{{\mathop \sum \nolimits _{j = 1}^N\eta _j^2}}\left( {{\rho _i}\left( \boldsymbol{r} \right) - \rho } \right)} \right].\end{align} \tag{ 12 }$$

We chose 40 $\mu {\text{m}}$ as the length scale, which is on the order of the radius of a typical grain in our simulations. We divide all the lengths by ${l_0}$ to obtain their dimensionless forms. We select the bulk energy density coefficient $W$ in equation (6) as the energy scale and divide all the energies by $W$ to nondimensionalize them. We denote ${\bar \rho _{{\text{init}}}}$ as the spatial average of $\rho (r,t = 0){\text{ }}$ over the entire computational domain (excluding the non-solid region defined by the dummy order parameter in section 4.1) and utilize this value to scale $\left\{ {{\rho _i}( \boldsymbol{r})} \right\}$ and $\rho \left( {\boldsymbol{r},t} \right)$. After these scaling quantities are chosen, we set the time scale to

$$\begin{align}{t_0} = \frac{1}{{WL}}.\end{align} \tag{ 13 }$$

We divide the time by this scale, and rearrange the resulting equation, which results in the dimensionless form of equation (12). The dimensional and dimensionless quantities are summarized in table 1. The dimensionless quantities are denoted by their counterpart accompanied by an asterisk ($*$) in their superscripts. We employ the dimensionless quantities for the study, which can be converted into dimensional quantities based on the scaling factors.

Table 1. Dimensional and dimensionless phase-field variables.

Variable	Nondimensionalized variable
$W$	${W^*} = \frac{W}{W} = 1$
$\kappa$	${\kappa ^*} = \frac{\kappa }{{Wl_0^2}} = \frac{1}{8}{\left( {\frac{{{l_{{\text{gb}}}}}}{{{l_0}}}} \right)^2}$
$L$	${L^*} = \frac{L}{L} = 1$
${\boldsymbol{r}}$	${\boldsymbol{r}^*} = \frac{\boldsymbol{r}}{{{l_0}}}$
$t$	${t^*} = \frac{t}{{{t_0}}} = WLt$
${\rho _i}$	$\rho _i^* = \frac{{{\rho _i}}}{{{{\bar \rho }_{{\text{init}}}}}}$
$\rho$	${\rho ^*} = \frac{\rho }{{{{\bar \rho }_{{\text{init}}}}}}$

The initial condition for each simulation is composed of two parts: an initial microstructure described by a set of order parameters $\left\{ {{\eta _i}\left( {\boldsymbol{r},t = 0} \right)} \right\}$, and a set of corresponding individual dislocation densities, $\left\{ {{\rho _i}\left( \boldsymbol{r} \right)} \right\}$. In this work, we perform two sets of 2D simulations. The first set of simulations directly imports a grain map from the experimental cross-sectional data containing a relatively small number of grains, while the second set of simulations leverages the GSD from the experimental data to initialize a synthetic microstructure that contains a large number of grains.

3.4.1. Initialization of order parameters.

First, the microstructure at state ${S_1}$ captured in the experiment is processed to produce a set of binary fields $\left\{ {{n_i}\left( {\boldsymbol{r}} \right)} \right\}$. For the first set of simulations, the 2D grain map, ${{\Omega }}\left( \boldsymbol{r} \right)$, of the confined region at ${S_1}$, which marks the region of each grain with a unique but consecutive integer, is used to initialize a set of ${n_i}\left( \boldsymbol{r} \right)$ fields according to

$$\begin{align}{n_i}\left( \boldsymbol{r} \right) = \left\{ \begin{array}{*{20}{l}} {1,\quad \,{\text{if}}{\text{ }}\Omega \left( \boldsymbol{r} \right) = i} \\ {0, \quad {\text{otherwise}}} \end{array}\right..\end{align} \tag{ 14 }$$

For the second set of simulations, the GSD calculated from the entire 3D volume of the sample at the ${S_1}$ state (figure 3(a)) is converted to a 2D GSD by implementing an established method in stereology [53, 54] in a reversed manner, as detailed in section S2 of the supplementary material. The 2D GSD is then used to inform the selection of the locations of $K$ seeds of grains, as detailed in section S3 of the supplementary material. A set of binary fields $\left\{ {{n_i}\left( \boldsymbol{r} \right)} \right\}$ is then constructed from the locations of $K$ seeds of grains, using the method discussed in section 3.5.

For both cases, the set of binary fields $\left\{ {{n_i}\left( \boldsymbol{r} \right)} \right\}$ is evolved via Allen–Cahn dynamics without the stored-energy for a dimensionless time of $t_{{\text{ini}}}^*$, just long enough to obtain the order parameters with diffuse interfaces, $\left\{ {{\eta _i}\left( \boldsymbol{r} \right)} \right\}$. The parameters employed for this smoothing process for the small- and large-scale simulations are summarized in tables 2 and 3, respectively. As mentioned earlier, the grain boundary energy for Cu–Zn [49], ${\gamma _{{\text{gb}}}} \approx 0.595$ J m⁻², is employed here for Cu–Al–Mn, which was consistent with Kusama's study [14]. For interfacial width, domain size, and grid spacing, both dimensionless and dimensional values are presented. Since bulk energy coefficient and gradient energy coefficient are model parameters that are calculated from the dimensionless values of grain boundary energy and grain boundary width, only dimensionless values of them are presented. The conversion of dimensionless time to dimensional time is complicated because it requires the mobility parameter $L$, which depends on the temperature and was not reported for Cu–Al–Mn alloy under a non-isothermal annealing condition. Thus, this study focuses on qualitative agreement between the dynamics observed in simulations and experiments. Note that the initial microstructure for the small-scale simulations is derived from the experimental measurement and contained tiny grains that cannot be resolved by the PF model. Therefore, it was necessary to apply the smoothing operation longer (about twice as long) until these grains are eliminated during the smoothing process. In addition, slight adjustments were needed in the grid spacing and related parameters due to the domain size change.

Table 2. Parameters employed for the small-scale simulations to generate the order parameters from the binary fields. Dimensional values in parentheses denote numerical parameters (not physical parameters).

Parameters	Variable	Value (dimensionless)	Value (dimensional)
Grain boundary energy	${\gamma _{{\text{gb}}}}$	1.67 $\times {10^{ - 2}}$	$0.595$ J m−2
Grain boundary width	${l_{{\text{gb}}}}$	0.1	(4 $\mu$m)
Domain size	$l$	$12$	4.80 $\times$ 102 $\mu$m
Grid spacing	$\Delta x = \Delta y$	$2.50 \times {10^{ - 2}}$	(1 $\mu$m)
Bulk energy coefficient	$W$	$1$	Not applicable
Gradient energy coefficient	$\kappa$	$1.25 \times {10^{ - 3}}$	Not applicable
Time step size	${{\Delta }}t$	$0.1$	Not available
Duration	${t_{{\text{ini}}}}$	100	Not available

Table 3. Parameters employed for the large-scale simulations to generate the order parameters from the binary fields. Dimensional values in parentheses denote numerical parameters (not physical parameters).

Parameters	Variable	Value (dimensionless)	Value (dimensional)
Grain boundary energy	${\gamma _{{\text{gb}}}}$	1.46 $\times {10^{ - 2}}$	$0.595$ J m−2
Grain boundary width	${l_{{\text{gb}}}}$	$8.75 \times {10^{ - 2}}$	(3.5 $\mu$m)
Domain size	$l$	$69.43$	2.78 $\times$ 103 $\mu$m
Grid spacing	$\Delta x = \Delta y$	$2.19 \times {10^{ - 2}}$	(0.876 $\mu$m)
Bulk energy coefficient	$W$	$1$	Not applicable
Gradient energy coefficient	$\kappa$	$9.57 \times {10^{ - 4}}$	Not applicable
Time step size	${{\Delta }}t$	$0.1$	Not available
Duration	${t_{{\text{ini}}}}$	50	Not available

3.4.2. Construction of individual dislocation density fields.

As seen in equation (3), the initial dislocation density ${\rho _{0,i}}\left( \boldsymbol{r} \right)$ and the distance function ${d_i}\left( \boldsymbol{r} \right)$ for each grain ${G_i}$ are two required components to construct an individual dislocation density field ${\rho _i}\left( \boldsymbol{r} \right)$. Although in reality a cycle of heat treatment involves continuous dislocation generation and microstructure evolution, we model this process by injecting all the generated dislocations, ${{\Delta }}{\rho _i}$, at the start of the thermal treatment cycle; thus, this term is added to the baseline value before the microstructure evolution is simulated. A more sophisticated approach that considers the dynamic generation of dislocations during precipitation will be left for future work. Given the experimental observations that a majority of the second-phase particles is found on or near the grain boundaries and the fact that a small grain has a higher line length per unit area in two dimensions (equivalent to surface area per unit volume in three dimensions) than a large grain, we propose that the density of second-phase particles is higher in a small grain. Additionally, it is reported in [16] that the relationship between the grain-averaged dislocation density and the density of second-phase particles can be roughly described by a positive linear correlation. We therefore postulate that the value of ${{\Delta }}{\rho _i}$ is inversely proportional to the equivalent radius (assuming a circular shape), ${R_i}$, of grain ${G_i}$, according to

$$\begin{align}{{\Delta }}{\rho _i} = \frac{{{k_\rho }}}{{{R_i}}}.\end{align} \tag{ 15 }$$

The value of ${k_\rho }$ is estimated by substituting the known values of ${{\Delta }}{\rho _1}$ and ${R_1}$ for grain ${G_1}$ within the confined region (see figure 2) into equation (15). To do so, ${{\Delta }}{\rho _1}$ is estimated by subtracting the baseline dislocation density from the average dislocation density of ${G_1}$ at ${S_2}$ within its initial position at ${S_1}$. The value of the equivalent radius ${R_1}$ is determined for ${G_1}$ based on its area at ${S_1}$. Accordingly, the value of ${k_\rho }$ is found to be $8 \times {10^8}$ m⁻¹.

When considering multiple cycles, each cycle must be initialized based on either the initial state of the sample or the state from the previous cycle. Hereafter in this section, we add the cycle number on the superscript of the dislocation density to indicate the cycle number. The individual dislocation density at positions initially within grain ${G_i}$ for the first cycle, $\rho _{0,i}^1\left( \boldsymbol{r} \right)$, is set by adding ${{\Delta }}{\rho _i}$ to the baseline dislocation density ${\rho _{{\text{base}}}}$:

$$\begin{align}\rho _{0,i}^1\left( \boldsymbol{r} \right) = {{\Delta }}{\rho _i} + {\rho _{{\text{base}}}} = \frac{{8 \times {{10}^8}{\text{ }}{{\text{m}}^{ - 1}}}}{{{R_i}}} + 3 \times {10^{12}}{\text{ }}{{\text{m}}^{ - 2}}.\end{align} \tag{ 16 }$$

To set the individual dislocation density at positions initially outside of grain ${G_i},$ we first apply the level-set method discussed in section S1 of the supplementary material to obtain the distance function, ${d_i}\left( \boldsymbol{r} \right)$, based on the order parameter ${\eta _i}\left( \boldsymbol{r} \right)$. Since $\rho _{0,i}^1\left( \boldsymbol{r} \right)$ is a constant for the first cycle, as indicated by equation (16) and illustrated in figure 4(b) (see the blue line), the individual dislocation density at the grain boundary, $\rho _{{\text{gb}},i}^1$, is equal to this constant value given by

$$\begin{align}\rho _{gb,i}^1 = {{\Delta }}{\rho _i} + {\rho _{{\text{base}}}} = \frac{{8 \times {{10}^8}{\text{ }}{{\text{m}}^{ - 1}}}}{{{R_i}}} + 3 \times {10^{12}}{\text{ }}{{\text{m}}^{ - 2}}.\end{align} \tag{ 17 }$$

Equation (3) is then employed to generate the individual dislocation density field $\rho _i^1\left( \boldsymbol{r} \right)$ for the first cycle.

To reinject dislocations at the beginning of each subsequent cycle $C{\text{ }}(C > 1)$, the value of the individual dislocation density at positions initially within ${G_i}$ for the subsequent cycle, $\rho _{0,i}^C\left( \boldsymbol{r} \right)$, is set by adding ${{\Delta }}{\rho _i}$ to its individual dislocation density at the end of the previous cycle, $\rho _{f,i}^{C - 1}\left( \boldsymbol{r} \right)$:

$$\begin{align}\rho _{0,i}^C\left( \boldsymbol{r} \right) = {{\Delta }}{\rho _i} + \rho _{f,i}^{C - 1}\left( \boldsymbol{r} \right) = \frac{{8 \times {{10}^8}{\text{ }}{{\text{m}}^{ - 1}}}}{{{R_i}}} + \rho _{f,i}^{C - 1}\left( \boldsymbol{r} \right).\end{align} \tag{ 18 }$$

We assume that after each cycle of grain growth, the value of dislocation density at the grain boundary of each remaining grain is reduced to a value that is sufficiently close to the baseline dislocation density. Thus, for simplicity, the individual dislocation density at the grain boundary for the cycle $C$ (where $C > 1$), $\rho _{{\text{gb}},i}^C$, is set again using equation (17). Then, equation (3) is employed to generate the initial dislocation density field $\rho _i^C\left( \boldsymbol{r} \right)$, where $C > 1$.

Although the implementation is simple when an order parameter is used to track one grain, the computation becomes expensive when the simulation involves a large number of grains, as in the second set of simulations. To enable large-scale simulation of microstructure evolution in an efficient manner, we require the capability to track multiple grains with one order parameter, which reduces the total number of order parameters required for the simulations. However, it is necessary to ensure that the grains tracked by the same order parameter do not overlap during the microstructure evolution to avoid grain coalescence (and, for our model, that the decaying individual dislocation density fields do not overlap, which would introduce numerical artifacts). The grain remapping schemes [32, 55, 56] have been proposed in the literature for PF simulations of grain growth, which reassign a grain to a different order parameter when the grain is too close to another grain tracked by the same order parameter. However, the standard remapping scheme cannot be applied without modifications because of the necessity to track the decay in the dislocation density outside of the grains. Therefore, we take a simpler approach in which grains are placed sufficiently far apart initially and are not dynamically remapped. In practice, we assign up to nine grains to one order parameter to reduce the number of required order parameters (and thus the number of equations to be evolved) and employ a corresponding individual dislocation density field to track the decay in dislocation densities for these grains. The detailed procedure is described in section S4 of the supplementary material.

In this section, we apply the PF model to study the evolution of microstructure and dislocation density within a confined region observed in the experiment (see figures 2(c)–(f)), as discussed in section 2. We select $4{\text{ }}\mu {\text{m}}$ as the interfacial width, which is ∼1/10 of the radius of a typical grain in our simulations and is sufficient to resolve the grains in the confined region. The order parameters are initialized using the method discussed in section 3.4.1. Since the boundary of the confined region is not rectangular, we utilize a dummy order parameter, which is not evolved throughout the microstructure evolution, to describe the area surrounding the confined region. The individual dislocation densities are initialized using equation (3), with ${\rho _{0,i}}\left( \boldsymbol{r} \right)$ given by equation (16) and ${\rho _{{\text{gb}},i}}$ given by equation (17). The dimensionless parameters employed for the simulation are summarized in table 4.

Table 4. Parameters employed in the simulation for microstructure evolution of the confined region. Dimensional values in parentheses denote numerical parameters (not physical parameters).

Parameters	Variable	Value (dimensionless)	Value (dimensional)
Grain boundary energy	${\gamma _{{\text{gb}}}}$	1.67 $\times {10^{ - 2}}$	$0.595$ J m−2
Grain boundary width	${l_{{\text{gb}}}}$	0.1	(4 $\mu$m)
Domain size	$l$	$12$	4.80 $\times$ 102 $\mu$m
Grid spacing	$\Delta x = \Delta y$	$2.5 \times {10^{ - 2}}$	(1 $\mu$m)
Bulk energy coefficient	$W$	$1$	Not applicable
Gradient energy coefficient	$\kappa$	$1.25 \times {10^{ - 3}}$	Not applicable
Time step	${{\Delta }}t$	$0.05$	Not available

The initial microstructure is shown in figure 5(a). Each grain is assigned to a unique number from 1 to 8. The color indicates the dislocation density within the grains. As shown in figure 5(b), grain 2 is rapidly consumed by the neighboring grains 3, 4, and 7. This behavior is expected since grain 2 possesses a higher dislocation density and a smaller size than the surrounding grains, leading to a large combined driving force due to capillarity and stored energy. As microstructure evolution proceeds, grains 3, 4, 5, and 6 are continuously consumed, as shown in figures 5(c) and (d). It is worth noting that the traversed regions that are far away from the initial boundary of the growing grains, such as the left regions of grains 1 and 8, have a nearly uniform deep blue color, implying the reduction of dislocation density to a baseline value. In addition, some of the grain boundaries, such as the one between grains 1 and 8 in figure 5(d), become curved during microstructure evolution. Such curved boundaries, often associated with island grains or peninsula grains, have been reported in samples that have undergone mechanical deformation and subsequent (isothermal) annealing leading to AGG [57, 58]. The simulation results thus indicate that such curved grains could also be induced by heat treatment alone without mechanical processing. Since grain 8 (with radius ∼92 $\mu$m) is slightly larger than grain 1 (with radius ∼90 $\mu$m), the assigned initial dislocation density for grain 8 (∼1.17 × 10¹³ m⁻²) is lower than that for grain 1 (∼1.19 × 10¹³ m⁻²). Consequently, grain 8 eventually consumes grain 1, as shown in figures 5(e) and (f), which does not agree with the experimental observation that grain 1 consumes grain 8, as shown in figures 2(c) and (d).

Zoom In Zoom Out Reset image size

The disagreement could arise from a few potential reasons. First, the simulation is confined to 2D, and the evolution of a cross section within a 3D volume may differ from that of a 2D simulation. Moreover, although the density of the second-phase particles is in general expected to be larger within a smaller grain, there may be statistical variations due to the stochastic nature of precipitate nucleation, growth, and subsequent dissolution that produce dislocations. Finally, it has been known that grain boundary energies are anisotropic, depending on not only the crystalline orientations of the pair of grains but also the relative orientation of the grain boundary plane [59]. The resulting variation in grain boundary energies would result in different probabilities of precipitating second-phase particles on or near grain boundaries, leading to a spread in the density of dislocations produced by the precipitates. Therefore, in reality, the dislocation density within a grain will likely not follow equation (15) exactly, and it is plausible that grain 1 has a lower dislocation density than grain 8. By conducting a parametric sweep, we determined that when the dislocation density of grain 1 is reduced by ∼16% to 10¹³ m⁻², grain 1 consumes grain 8 during the simulated evolution. The results of this simulation are shown in figure 6. The early evolution of the two simulations are very similar, as can be seen by comparing figures 5(a)–(c) and 6(a)–(c). On the other hand, during the later stage, grain 1 progressively consumes grain 8, as shown in figures 6(d) and (e), and it is the only grain that remains at the end of the simulation, as shown in figure 6(f). The high dislocation density at grain 1's initial position and the low dislocation density in the surrounding region closely resemble the experimental observations from figure 2(f). Although a quantitative agreement is not expected because of the 2D nature of the simulations and is therefore not a goal of this work, the observed similarities provide support for the proposed model.

Zoom In Zoom Out Reset image size

In this section, we perform a set of large-scale simulations to examine the effect of CHT on the microstructure evolution. To maintain computational feasibility, the smallest resolvable grain in the simulations is chosen to have a radius of $12{\text{ }}\mu {\text{m}}$. This grain radius requires the interfacial width (defined by order parameter values from ∼0.1 to ∼0.9) to be no larger than 3.5 $\mu {\text{m}}$ such that the bulk region within the smallest grain occupies at least 50% of the area, assuming a circular shape of the grain. We initialize the order parameters using the method discussed in sections 3.4.1 and 3.5. The individual dislocation densities are initialized using equation (3). For this equation, the values of ${\rho _{0,i}}\left( \boldsymbol{r} \right)$ are given by equation (16) when the dislocations are injected for the first time and by equation (18) when the dislocations are injected for the second and the third times. The values of ${\rho _{{\text{gb}},i}}$ are given by equation (17). The dimensionless parameters employed for the simulations are summarized in table 5. Periodic boundary conditions are utilized throughout the simulations.

Table 5. Parameters employed in the large-scale simulations. Dimensional values in parentheses denote numerical parameters (not physical parameters).

Parameters	Variable	Value (dimensionless)	Value (dimensional)
Grain boundary energy	${\gamma _{{\text{gb}}}}$	1.46 $\times {10^{ - 2}}$	$0.595$ J m−2
Grain boundary width	${l_{{\text{gb}}}}$	$8.75 \times {10^{ - 2}}$	(3.5 $\mu$m)
Domain size	$l$	$69.43$	2.78 $\times$ 103 $\mu$m
Grid spacing	$\Delta x = \Delta y$	$2.19 \times {10^{ - 2}}$	(0.876 $\mu$m)
Bulk energy coefficient	$W$	$1$	Not applicable
Gradient energy coefficient	$\kappa$	$9.57 \times {10^{ - 4}}$	Not applicable
Time step (with stored energy)	${{\Delta }}{t_{\text{s}}}$	$0.05$	Not available
Time step (without stored energy)	${{\Delta }}{t_{{\text{ns}}}}$	$0.1$	Not available

The initial microstructure, which contains 1778 grains, is shown in figure 7(a). It can be observed from the figure that a number of small grains are clustered in a chain-like structure. This arrangement is reasonable in 2D given the facts that a large portion of small grains was observed in the experiment at state ${S_1}$ (see figure 3(a)) and that these small grains would have to cluster in order to avoid an unphysical (non-equiaxed) grain morphology. The high dislocation density of these small grains yields a high stored energy as compared to their surrounding grains, and consequently, these grains are rapidly consumed, as shown in figure 7(b). The microstructure at the end of the first cycle is shown in figure 7(c). Most grains with high dislocation densities (i.e. yellow grains) are consumed by their surrounding grains at this stage and the large grains next to these small grains become larger. The corresponding GSDs are shown in figure 8. Note that the probability is scaled to yield an area of unity under the curve. It can be seen that GSD broadens during this cycle (figures 8(a) and (b)), primarily due to the increase in the average grain size. A second cycle is then simulated, which consists of reinjections of dislocations to the microstructure and subsequent microstructure evolution. The increase of dislocation densities is inversely proportional to the grain radius (equation (15)), and therefore the smaller grains experience a higher increase in dislocation densities, as shown in figure 7(d). Microstructures at an intermediate stage and at the end of the second cycle are shown in figures 7(e) and (f), respectively. The GSD continues to broaden, as shown in figure 8(c). A second peak in the GSD appears after the third cycle, as shown in figure 8(d). The comparison of scaled GSDs for each cycle is provided in figure S4(a) of the supplementary material, in which all the radii are divided by the average radius. The scaled GSD for the simulation with stored energy is clearly not self-similar. Instead, it continues to evolve throughout the three simulation cycles, exhibiting a slight broadening especially at the later cycles, and a bimodal distribution may start to develop at the end of the third cycle. It has been suggested that a bimodal GSD is one of the signatures of AGG [2, 3].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that the GSD calculated from the simulated data (figure 8(b)) do not match the experimental result in the ${S_3}$ state (see figure 3) on the extreme ends of the distribution. Besides the difference in the dimensionality between the simulation and experiment, there are additional reasons for the discrepancies. On the low end of the radii, the simulation misses some of the smallest grains due to the diffuse interface effect, while the experiment may overestimate the number of small grains due to the noise in the data. On the high end of the radii, the experiment had temperature gradient over the sample, which led to faster grain growth in the hotter end, where the large grains were identified. Furthermore, the timing of the end of the cycle is approximated by the increase in the radius, which introduces additional uncertainty.

To confirm that dislocation injections are necessary to yield persistent AGG [60], a simulation that omits the driving force due to the stored energy difference is performed. To allow for a direct comparison, an identical initial microstructure is employed, as depicted in figure 9(a). The microstructure at times ${t^*} = 800$, $1700$, and $3300$ during the evolution is shown in figures 9(b)–(d), respectively, which correspond to the final net times of the three cycles for the stored-energy-driven grain growth (see figures 7(c), (f), and (i)). In this case, the large grains still progressively consume the smaller grains surrounding them, due to capillarity. To quantify and compare the rate of the microstructure evolution between the simulations with and without the stored energy, we show the average grain radius as a function of time for the two cases in figure 10. The average grain radius during the microstructure evolution with stored energy for cycles 1, 2, and 3 are indicated by red, blue, and black dots, respectively. The average grain radius during the microstructure evolution without stored energy is indicated by green dots. A much faster increase in the average radius is observed when stored energy is considered, which indicates that the microstructure evolution proceeds more rapidly when dislocations are injected via CHT. Additionally, at the start of each cycle, there is a rapid increase in the average radius, which is due to the large variation in initial dislocation densities. Then, the growth rate decreases as the stored-energy driving force reduces and the system shifts toward NGG. An additional abrupt increase in the average radius in the middle of cycle 3, observed from time ${t^*} = 3000$ to $3100$, is caused by a significant fraction of grains (14 out of 246 grains) that disappeared between these two times. Unlike the apparent broadening of the GSD observed in the evolution of microstructure with dislocations, the GSD for capillary-driven grain growth only slightly broadens (see figures 8(e)–(g)). In addition, none of the grains are observed to grow abnormally faster than the others. Consequently, the GSD remains singly peaked from ${t^*} = 800$ to $3300$. The convergence to self-similarity is observed in the scaled GSD, as shown in figure S4(b) of the supplementary material.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To further examine the necessity of recurrent dislocation injections to continuously provide stored-energy driving force, we performed an additional simulation in which the dislocations are only injected once at the beginning. The average grain radius as a function of time from this simulation is shown in figure S3 of the supplementary material, in which red crosses indicate the continued portion, along those of the three-cycle simulation. The microstructure evolution with one-time dislocation injection gradually turns slower, and eventually the rate reduces to a value that is close to that of capillary-driven grain growth, suggesting a progressive transition from stored-energy-driven grain growth to capillary-driven grain growth. In order to prevent such transition, dislocation reinjections are necessary, which provide continuous stored-energy driving force to maintain the rate of microstructure evolution.