Equation of motion for grain boundaries in polycrystals

Abstract

Grain boundary (GB) dynamics are largely controlled by the formation and motion of disconnections (with step and dislocation characters) along with the GB. The dislocation character gives rise to shear coupling; i.e. the relative tangential motion of two grains meeting at the GB during GB migration. In a polycrystal, the shear coupling is constrained by the presence of other grains and GB junctions, which prevents large-scale sliding of one grain relative to the other. We present continuum equations of motion for GBs that is based upon the underlying disconnection dynamics and accounts for this mechanical constraint in polycrystals. This leads to a reduced-order (zero-shear constrained) model for GB motion that is easily implemented in a computationally efficient framework, appropriate for the large-scale simulation of the evolution of polycrystalline microstructures. We validated the proposed reduced-order model with direct comparisons to full multi-disconnection mode simulations.

Introduction

Grain boundaries (interfaces between domains/grains of different crystallographic orientations) are major components of polycrystalline microstructures. The dynamical properties of grain boundary (GB) networks play essential roles in the mechanical response of materials and influence a wide range of other properties (e.g., electrical resistance, optical transmission, thermoelectric efficiency)¹. Traditional models for isotropic GB dynamics are based on capillarity-driven GB migration; i.e., the driving force for migration is proportional to the local GB mean curvature κ^2,3 such that the normal GB velocity v_n (isotropic limit) is v_n = M_GBγκ, where M_GB is the GB mobility and γ is the GB energy. Theoretical predictions and experimental observations demonstrate that GB migration, in general, is controlled by the motion of discrete line defects constrained to the GBs known as disconnections^{4,5,6,7,8,9,10,11,12,13}. Motion by mean curvature-based continuum models are unable to adequately account for many GB migration phenomena such as stress-driven GB motion associated with disconnection dynamics.

Disconnection is a line defect that is characterised by a dislocation Burgers vector b and a step height H; i.e., (b, H). We recently developed a continuum model for GB migration that is based upon disconnection migration^14,15. These models naturally include GB migration under a wide range of driving forces (internal and applied stress, chemical potential jumps, capillarity), as well as the effect of temperature. We also proposed a continuum model that accounts for the concurrent migration of GB and the junctions along which they meet in a disconnection-based framework^16,17.

One important consequence of GB migration by disconnection motion is shear coupling^18,19,20; the translation of one grain relative to the other in a direction parallel to the GB plane during migration (associated with the disconnection Burgers vector) (Fig. 1a, b). The shear coupling has been observed in bicrystals in both experiments and simulations^{21,22,23,24,25,26,27,28,29,30,31,32,33}. On the other hand, in a polycrystal, constraints imposed by the presence of other grains prevents large-scale sliding of one grain relative to the other (Fig. 1c). Grain boundary migration in such constrained systems requires the simultaneous action of disconnections of multiple modes such that there is zero-net-shear¹⁶. Continuum models^14,15 can, in principle, incorporate this collective behaviour through the long-range stress field generated by all of the disconnections on all of the GBs in the microstructure^16,17. However, such calculations are too complex and computationally inefficient for widespread adoption in microstructure evolution simulations.

Fig. 1: Schematic illustration of shear coupling and constraints during GB migration.

A bicrystal with free top and bottom surface a before and b after GB migration with accompanying shear deformation. The dashed blue line is a fiducial marker (e.g., a scratch) in the bicrystal. The shear-coupling factor β is the inverse of the slope of this line between the initial and final GB positions (red solid and dashed lines in b, respectively). c A grain in a polycrystal with GBs 1–5 and Junctions A–E. The shear of GB ① is constrained by junctions A and B and by the elastic interactions between Grain I and surrounding grains, III–VI. d GB migration in the bicrystal of a where the top and bottom surface do not translate to mimic the effects of Grains IV and V on the shear of GB ①.

In this paper, we propose a reduced-order description of GB migration that accounts for the constraints imposed by other grains in the polycrystalline microstructure (Fig. 1c) that is both disconnection-based and is sufficiently computationally tractable to be broadly applied. The result is a disconnection migration-based equation of motion for a mechanically constrained GB; one for which there is zero-net-shear (see an example in Fig. 1d). To validate our proposed approach, we compare model predictions with a multi-disconnection migration-based approach for GB migration in the presence of different mechanical constraints (as appropriate for microstructural evolution in a polycrystal)¹⁵.

Results

Equation of motion for a constrained grain boundary

On the microscopic level, a high-angle GB may be viewed as being composed of flat sections and disconnections (the stepped, black curve in Fig. 2). The motion of disconnections in one direction translates the GB while motion of a pair of disconnections such that they meet and annihilate with each other changes the GB curvature. Hence, both GB migration and change in GB shape can be characterised by disconnection motion. On a continuum level, a high-angle GB may be modelled as a smoothly varying curve (not resolving the underlying disconnection/step structure—see the red curve in Fig. 2).

Fig. 2: Continuum representation of GB with disconnections.

A GB with disconnections (black curve) and its continuum representation y = h(x, t) (red curve). The GB velocity v (in direction y) results from the glide of disconnections in x, characterised by (b, H) and (−b, −H).

Consider a slightly curved GB, deviating from a flat reference parallel to x. The GB profile may be described by y = h(x, t), where ∣∂_xh∣ ≪ 1. The disconnection mode is characterised by a Burgers vector and step height (b, H), where the Burgers vectors b = (b, 0) is, here, in x and the step height vector H in y (Fig. 2). There are many possible Burgers vectors b_n and for each a set of possible step heights {H_nk} on each GB, as determined by the misorientation between the grains¹³. Although in principle, the number of possible disconnection modes is infinite, the modes likely to be found are quite limited, depending on thermodynamic variables, bonding, disconnection core structure. The density of disconnections of mode j (b^(j), H^(j)) is denoted ρ^(j) and positive/negative ρ^(j) value represents the density of disconnection (b^(j), H^(j))/(−b^(j), −H^(j)). Evolution equations for the GB profile h(x, t) and disconnection densities ρ^(j)(x, t) on the GB can be described by disconnection fluxes J^(j),

$${h}_{t}=-\mathop{\sum}\limits_{j}{J}^{(j)}{H}^{(j)},$$

(1)

$${\rho }_{t}^{(j)}+{J}_{x}^{(j)}=0.$$

(2)

In a polycrystalline microstructure, two grains meeting at a GB cannot slide with respect to one another without limit. This is because of the mechanical constraint from the other grains surrounding the two grains meeting at the GB in question and/or because of the presence of a nearby GB junction^12,13. (Note that since different GBs have different crystallographically allowed disconnection Burgers vectors and step heights, disconnections can rarely move across GB junctions.) These constraints give rise to a macroscopic zero-net-shear across a GB constraint; i.e.

$$\dot{\epsilon }=\mathop{\sum}\limits_{j}{b}^{(j)}{J}^{(j)}=0.$$

(3)

Since disconnections have both dislocation and step character, this implies a constrained evolution problem: the GB profile and disconnection densities evolve as per (1) and (2) subject to the constraint (3).

The evolution of the GB profile h_t is determined by the disconnection fluxes (1). The flux of disconnections of each mode contributes to the change in GB profile via their step heights H^(j). The disconnection density evolution is consistent with both the Burgers vector and step conservation. The flux J^(j) of the jth disconnection mode is

$${J}^{(j)}={v}^{(j)}(| {\rho }^{(j)}| +2{c}^{(j)}),$$

(4)

where ${c}^{(j)}=\frac{1}{a}{{\rm{e}}}^{-{E}^{* (j)}/({k}_{{\rm{B}}}T)}$ is the thermal equilibrium disconnection density on the GB (which enables the motion of an initially flat GB), a is an atomic spacing, k_B is the Boltzmann constant, T is the temperature, E*^(j) is half the disconnection pair formation energy¹³, and v^(j) is the disconnection glide velocity

$${v}^{(j)}={M}_{{\rm{d}}}^{(j)}[\tau {b}^{(j)}+\Psi {H}^{(j)}-\gamma {h}_{xx}{H}^{(j)}].$$

(5)

Here, ${M}_{{\rm{d}}}^{(j)}$ is the disconnection mobility, τ is the total internal shear stress, and Ψ is the chemical potential-jump across the GB. We assume that the disconnection glide velocity (5) is overdamped, such that it is proportional to the driving force ${f}_{{\rm{d}}}^{(j)}$, i.e. ${v}^{(j)}={M}_{{\rm{d}}}^{(j)}\cdot {f}_{{\rm{d}}}^{(j)}$.

With (4) and (5), the assumption of zero-net-shear across the GB as described in (3) becomes

$$\dot{\epsilon }=\mathop{\sum}\limits_{j}{A}^{(j)}\left[\tau +(\Psi -\gamma {h}_{xx})\frac{1}{{\beta }^{(j)}}\right]=0,$$

(6)

where

$${\beta }^{(j)}={b}^{(j)}/{H}^{(j)},$$

(7)

$${A}^{(j)}={M}_{{\rm{d}}}^{(j)}{{b}^{(j)}}^{2}(| {\rho }^{(j)}| +2{c}^{(j)}).$$

(8)

Here β^(j) is the shear-coupling factor associated with disconnection mode j¹⁹. The zero-net-shear constraint (3) implies the development of a shear stress τ on the GB induced by GB migration:

$$\tau =-(\Psi -\gamma {h}_{xx})\frac{{\sum }_{j}\frac{1}{{\beta }^{(j)}}{A}^{(j)}}{{\sum }_{j}{A}^{(j)}}.$$

(9)

With this effective shear stress, the GB profile and disconnection density evolution, (1) and (2), can be expressed as

$${h}_{t}=-(\Psi -\gamma {h}_{xx})\frac{{\sum }_{k \,{>}\,j}{\left(\frac{1}{{\beta }^{(j)}}-\frac{1}{{\beta }^{(k)}}\right)}^{2}{A}^{(j)}\cdot {A}^{(k)}}{{\sum }_{k}{A}^{(k)}},$$

(10)

$${\rho }_{t}^{(j)}=-\frac{\partial }{\partial x}\left((\Psi -\gamma {h}_{xx})\frac{1}{{b}^{(j)}}\cdot \frac{{\sum }_{k}\left(\frac{1}{{\beta }^{(j)}}-\frac{1}{{\beta }^{(k)}}\right){A}^{(j)}\cdot {A}^{(k)}}{{\sum }_{k}{A}^{(k)}}\right),$$

(11)

where β^(j) and A^(j) are given in (7) and (8). This is our reduced-order model for GB motion and disconnection density evolution subject to the zero-net-shear constraint.

Solution of the original dynamical equations (1) and (2) requires the determination of the total, non-local, shear stress τ. However, with the zero-net-shear constraint, the dynamical equations (10) and (11) are local and can be solved much more efficiently than the original ones.

Before concluding our discussion of a GB equation of motion, it is appropriate to note a few limitations. First, the full disconnection description of GB migration explicitly depends on the crystallinity of the grains on either side of the GB. It will not describe GB migration in systems where the GB is a thick amorphous between the crystal grains. Second, the disconnection model described here assumes that the disconnections move in a glissile manner; disconnections with Burgers vector components perpendicular to the GB plane may move by absorption/emission of point defects - disconnection climb. We do not consider climb here since the climb is typically much slower than glide (even on a GB) and, hence, disconnection climb is expected to make a much smaller contribution to GB migration than that associated with glide. Third, we ignore the fact that GBs may exhibit multiple metastable states³⁴, in which case, disconnections may be partial disconnections separated by stacking faults. While such effects may be included in a discrete disconnection model, they are not incorporated here.

Flat grain boundary dynamics

We first apply the GB evolution equations to GB migration driven by a chemical potential-jump Ψ across the GB. Consider the simple case of an infinitely large, initially flat GB h(x, 0) = 0 with no disconnections ρ^(j) = 0 for all j. Here, the GB profile evolves as

$${h}_{t}=-{M}_{{\rm{GB}}}\Psi ,$$

(12)

where

$${M}_{{\rm{GB}}}=\frac{{\sum }_{k \,{>}\,j}{\left(\frac{1}{{\beta }^{(j)}}-\frac{1}{{\beta }^{(k)}}\right)}^{2}{M}_{{\rm{d}}}^{(j)}{{b}^{(j)}}^{2}2{c}^{(j)}\cdot {M}_{{\rm{d}}}^{(k)}{{b}^{(k)}}^{2}2{c}^{(k)}}{\mathop{\sum }\nolimits_{j = 1}^{J}{M}_{{\rm{d}}}^{(j)}{{b}^{(j)}}^{2}2{c}^{(j)}}.$$

(13)

This expresses the GB velocity h_t as the product of the chemical potential-jump Ψ and an effective GB mobility which is determined by the properties/parameters of the disconnections. Here the GB remains flat and moves with a constant velocity and the disconnection densities are zero ρ^(j) = 0 for all j at all time (i.e., the equal density of positive and negative disconnections). The internal stress τ (9) is proportional to the chemical potential-jump Ψ with a negative factor that depends on the disconnection properties; this exerts a back stress that opposes the chemical potential-jump Ψ,

$$\tau =-\Psi \frac{{\sum }_{j}{M}_{{\rm{d}}}^{(j)}{{b}^{(j)}}^{2}{c}^{(j)}\frac{1}{{\beta }^{(j)}}}{{\sum }_{j}{M}_{{\rm{d}}}^{(j)}{{b}^{(j)}}^{2}{c}^{(j)}}.$$

(14)

We now compare the zero-net-shear model (see (12) and (13)) with the predictions of the full continuum disconnection model¹⁵, as represented by (1), (2), (4), and (5) for the migration of a flat GB. More specifically, we assume a bicrystal geometry where the top and bottom surfaces of the bicrystal (flat horizontal, centred GB) are fixed (see Fig. 1d and ref. ¹⁵ for details) to provide no net bicrystal shear. In the full continuum disconnection model, we explicitly evolve two-disconnection modes and evaluate the stresses at every point (x, y) in the bicrystal that arise from all disconnections along the GB and the image stresses generated by the fixed top and bottom surface of the bicrystal boundary conditions.

We choose parameters appropriate for a nominally flat, Σ13[100](015) symmetric tilt GB in Cu driven by a chemical potential-jump of Ψ = −1 meV ⋅ Å⁻³. More specifically, we consider the two lowest formation energy disconnections modes (b⁽¹⁾, H⁽¹⁾) and (b⁽²⁾, H⁽²⁾); i.e. ${b}^{(1)}={b}^{(2)}={a}_{0}/\sqrt{26}$, ${H}^{(1)}=4{a}_{0}/\sqrt{26},{H}^{(2)}=-2.5{a}_{0}/\sqrt{26}$, where the lattice constant a₀ = 3.615 Å^13,15,35. When T = 800 K, the nucleation rate of these disconnections are estimated to be 2c⁽¹⁾ = 0.2188/L₀ and 2c⁽²⁾ = 0.6930/L₀, where L₀ = 100 Å is the length of the GB^13,15,36,37. In these simulations, we set the distance between the top and bottom surfaces as 2L₀ and the GB is initially located in the centre of the simulation cell. These are consistent with a grain size of ~L₀.

The migration velocity of this flat GB (see (12) and (13)) is independent of time, as seen in Fig. 3a for T = 800 K (red line). This constant velocity agrees well with that found at a late time from the simulation based upon the full multi-mode disconnection model¹⁵ (the blue curve). In the full multi-mode disconnection simulation, the velocity starts at a large value but relaxes quickly to the steady-state solution founding the zero-shear model.

Fig. 3: Migration of a flat GB.

Migration of a nominally flat, Σ13[100](015) symmetric tilt GB in Cu driven by a chemical potential-jump across the GB of Ψ = −1 meV ⋅ Å⁻³ using the zero-shear, reduced-order model and the full multi-disconnection mode model¹⁵. The GB migration velocity in a and internal stress in b at T = 800 K predicted by the reduced-order model (red lines) and the full multi-disconnection mode model with fixed boundary conditions (blue curves). The predicted steady state GB migration velocity in c and internal stress in d are shown as a function of T (blue points: full multi-disconnection mode model). Stress is normalised by τ₀ = − ΨH⁽²⁾/b⁽²⁾, GB velocity by the T = 1300 K velocity v₀, and t₀ = L₀/v₀.

This transient in the full multi-mode disconnection simulation can be understood as follows: (i) the GB starts to migrate by the motion of the lowest formation energy disconnection mode, (ii) shear coupling and the constraint set up a back stress that counteracts the chemical jump driving force¹², (iii) continued motion requires the formation of disconnection of opposite shear-coupling factor with higher formation energy and (iv) continued migration is limited by the rate of formation of the more difficult to form secondary (or higher order) disconnections. The concerted motion of these different modes is required for steady-state migration.

This ‘back stress’ picture is confirmed by a consideration of the internal stress during GB migration τ as described by (14) and the full continuum model¹⁵ (Fig. 3b). While the zero-shear model accurately predicts the steady-state internal stress, the full multi-mode disconnection simulation shows the transient development of the internal stress which asymptotically approaches the reduced-order model steady-state stress. Figure 3c, d show the T-dependence of the GB mobility and the internal stress from the zero-shear model and the full multi-mode disconnection simulations. The agreement is excellent. The duration of the transient, not captured by the zero-shear model, depends on both the GB bicrystallography as well as the microstructure (e.g. large grain sizes imply longer transients). It is important to note that since the full multi-mode continuum disconnection model successfully reproduces atomistic simulation results¹⁵ and the reduced-order model agrees with the steady-state of that model, the reduced-order model too is consistent with the steady-state obtained from atomistic simulations¹².

While the flat GB velocity in the zero-shear model is explicitly given by (12), in the full multi-disconnection mode model¹⁵ it depends on the shear stress τ which evolves during disconnection/GB migration during migration (must be determined via elasticity calculations with the appropriate boundary conditions). Calculation of the stress generated by disconnections in the full model has computational complexity ${\mathcal{O}}\left((m+1){N}^{2}\right)$, where m and N are the number of disconnection modes and points in the GB mesh. Calculation of the image stress arising from the boundary conditions in a general full multi-disconnection model is calculated over the entire domain via, for example, the finite element method. (For the simulation geometries/domain considered in the full multi-disconnection mode model, we employ the Airy stress function and Fourier analysis methods.) Elimination of the stress field calculation makes the reduced-order model local and efficient.

Curved grain boundary dynamics

We now consider the dynamics of the capillarity-driven migration of a curved GB. The geometric constraint that connects the disconnection densities with the GB profile is

$$\mathop{\sum }\limits_{j=1}^{m}{\rho }^{(j)}{H}^{(j)}={h}_{x}.$$

(15)

The GB profile and disconnection density evolution equations (1) and (2) imply that $\left(\right.\mathop{\sum }\nolimits_{j = 1}^{m}{\rho }^{(j)}{H}^{(j)}-{h}_{x}{\left)\right.}_{t}=0$. Thus, the geometric constraint (15) applies at all times. The disconnection density evolution and the zero-net-shear constraints, (2) and (3), imply that $\mathop{\sum }\nolimits_{j = 1}^{m}\left(\right.{\rho }^{(j)}{b}^{(j)}{\left)\right.}_{t}=0$; i.e.

$$\mathop{\sum }\limits_{j=1}^{m}{\rho }^{(j)}{b}^{(j)}=0.$$

(16)

This means that the net Burgers vector vanishes at every point on the GB.

We first consider the two-disconnection mode cases: (b⁽¹⁾, H⁽¹⁾) and (b⁽²⁾, H⁽²⁾), with initial disconnection densities chosen to satisfy constraints (15) and (16):

$${\rho }^{(1)}{H}^{(1)}+{\rho }^{(2)}{H}^{(2)}={h}_{x},$$

(17)

$${\rho }^{(1)}{b}^{(1)}+{\rho }^{(2)}{b}^{(2)}=0.$$

(18)

This yields an explicit expression for the disconnection densities,

$${\rho }^{(1)}=\frac{{b}^{(2)}}{{b}^{(2)}{H}^{(1)}-{b}^{(1)}{H}^{(2)}}{h}_{x},\quad {\rho }^{(2)}=\frac{{b}^{(1)}}{{b}^{(1)}{H}^{(2)}-{b}^{(2)}{H}^{(1)}}{h}_{x}.$$

(19)

Substituting these into the GB profile evolution equation (10), we obtain an explicit expression for the GB evolution:

$${h}_{t}=-B(\Psi -\gamma {h}_{xx})\left(| {h}_{x}| +D\right),$$

(20)

where

$$\begin{array}{lll}B\,=\,{\left(\frac{1}{{\beta }^{(1)}}-\frac{1}{{\beta }^{(2)}}\right)}^{2}\frac{{M}_{{\rm{d}}}^{(1)}{M}_{{\rm{d}}}^{(2)}{{b}^{(1)}}^{2}{{b}^{(2)}}^{2}}{{M}_{{\rm{d}}}^{(1)}| {b}^{(1)}| +{M}_{{\rm{d}}}^{(2)}| {b}^{(2)}| }\frac{1}{| {b}^{(2)}{H}^{(1)}-{b}^{(1)}{H}^{(2)}| },\\ D\,=\,\frac{2\left({M}_{{\rm{d}}}^{(2)}{c}^{(1)}+{M}_{{\rm{d}}}^{(1)}{c}^{(2)}\right)}{{M}_{{\rm{d}}}^{(1)}| {b}^{(1)}| +{M}_{{\rm{d}}}^{(2)}| {b}^{(2)}| }| {b}^{(2)}{H}^{(1)}-{b}^{(1)}{H}^{(2)}| .\end{array}$$

Note, we assumed that the equilibrium disconnection densities c⁽¹⁾ and c⁽²⁾ are small and only keep their leading order terms.

We performed a series of numerical simulations to study the evolution of a curved GB with two-disconnection modes using the reduced-order (20) and full continuum models¹⁵. We again focus on the Σ13[100](015) (near) symmetric tilt GB in Cu (shear modulus μ = 45 GPa and Poisson ratio ν = 0.36) in the two-mode limit (parameters above); the symmetric GB has energy γ = 0.878 J ⋅ m^−2 19. The GB has an initially sinusoidal profile $h=0.03{L}_{0}\sin (2{\rm{\pi }}x/{L}_{0})$ with period L₀ = 100 Å. While the initial disconnection densities in the full multi-disconnection mode model are arbitrary chosen (we examine two cases below), in the reduced-order model these satisfy (19).

First, we examine the case where the initial disconnection densities in the reduced-order and full model are identical (19) (Fig. 4a–c). As shown in Fig. 4a, b, the reduced-order model predictions for the evolution of both the GB profile and densities of the two-disconnection modes are in excellent agreement with the predictions from the full continuum model with the same number of modes. To provide a more quantitatively comparison of the two sets of predictions, we also calculated the differences of the GB migration velocities and the differences in the rates of changes of the disconnection densities between the reduced-order model and the full model, as shown in Fig. 4c. The differences between the two models decay rapidly to zero; indicating again that the errors intrinsic to the reduced-order model are associated with small, short-lived transients.

Fig. 4: Comparison of the reduced-order model and the full multi-disconnection mode model.

The temporal evolution of the GB profile in a, d and density in b, e of each disconnection mode using the reduced-order model (solid lines) and the full multi-disconnection mode model (dashed lines with circles) for two-disconnection modes. The arrows indicate the direction of the evolving curves. c, f The difference between the predictions of the reduced-order model and the full multi-disconnection mode model for the GB velocity and rates of change of the disconnection densities during the evolution of the initially sinusoidal profile GB. These data are calculated at the positions of the arrows in a, b, d and e. Panels a–c correspond to the case where the initial disconnection densities in both models are the same (19), and d–f correspond to a case where the initial densities are different—as described in the text.

In the second comparison, we choose the initial disconnection densities in the full multi-disconnection mode model to be different from those in (19) in the reduced-order model; i.e. the initial disconnection densities in the full multi-disconnection mode model was ρ⁽¹⁾ = h_xc⁽¹⁾/(c⁽¹⁾H⁽¹⁾ − c⁽²⁾H⁽²⁾), ρ⁽²⁾ = h_xc⁽²⁾/(c⁽²⁾H⁽²⁾ − c⁽¹⁾H⁽¹⁾). The comparisons, shown in Fig. 4d–f, indicate that the GB profiles evolution are well-predicted by the reduced-order model and the full multi-disconnection mode model disconnection densities rapidly decay to those predicted from the reduced-order model. While the differences decay rapidly, the magnitude of the initial differences in the evolution rates (Fig. 4f) is an order of magnitude larger than when the initial disconnection densities are the same, as expected. These results demonstrate that there are two sources of error associated with the reduced-order model. The first is associated with the zero-shear assumption (the heart of the reduced-order formalism) and the second is associated with the initial disconnection densities which are inevitably unknown (in most applications). Fortunately, both of these only affect transient behaviour, such that the reduced-order and full multi-disconnection mode models match after a short-time transient.

In the present two-disconnection-mode examples, the constraints (17) and (18) in the reduced-order model uniquely determine the densities of the two modes of disconnections (19) and the GB evolution equation (20) is an explicit expression for the GB profile evolution, without requiring the tracking of the evolution of the disconnection densities. However, in the general m-disconnection modes (m ≥ 3) case, the two constraints (15) and (16) are not sufficient to uniquely determine the disconnection densities, i.e. there is an infinite set of possible disconnection densities for any given GB profile h. In this case, we must solve both the equations for the GB motion (10) and the evolution of the disconnection densities (11) in the reduced-order model; here the disconnection density and GB profile evolution are coupled; see Supplementary Fig. 1 for simulation results for GB evolution with three disconnection modes.

While the results presented above demonstrated that the reduced-order model provides accurate predictions for the evolutions of arbitrarily curved GBs and the disconnection densities on them, a major advantage of the reduced-order model (as compared with the full continuum multi-disconnection mode model) is computational efficiency. As discussed in the previous section, the reduced-order model avoids the time-consuming ${\mathcal{O}}({N}^{2})$ calculation (N is the number of numerical grid/mesh points on GBs) of singular integrals for the stress generated by disconnections and the auxiliary elasticity problem over the entire computational domain associated with elastic boundary conditions. Hence, the reduced-order model provides a simple, efficient and accurate method to track GB motion based upon a mechanistically accurate disconnection mechanism.

Discussion

We presented a reduced-order, local continuum model for GB migration appropriate for microstructure evolution in polycrystals that accounts for mechanical constraints through the zero-net-shear condition. Our reduced-order model is based on the underlying microscopic mechanism of GB dynamics, in which GB migration and shear are associated with disconnection migration along with a GB. Our model leads to explicit evolution equations for the GB profile and disconnection densities. The strong zero-net-shear constraint implies that the reduced-order model does not capture short-time, transient effects but accurately describes long-time (steady-state) GB dynamics. The zero-net-shear assumption implies that the reduced-order model does not capture GB sliding or grain rotation.

When the GB is flat, Eq. (12) implies that the GB migration velocity is proportional to the driving force (chemical potential-jump Ψ). The migration of a general curved GB is described by h_t = − B(Ψ − γh_xx)(∣h_x∣ + D), where B and D are functions of disconnection parameters (20). D is associated with the equilibrium disconnection densities c^(j) and is small at low T. At high T, with large equilibrium disconnection densities c^(j), D is large and ∣h_x∣ ≪ D. In this limit, the GB equation of motion reduces to classical curvature flow h_t = − M_GB(Ψ − γh_xx) (M_GB = BD).

The reduced-order model implies that two or more disconnection modes must be active on a GB to satisfy the zero-net-shear constraint (unless the lowest formation energy disconnection is a pure step). If there is only one disconnection mode active on a GB, the constraint (3) reduces to $\dot{\epsilon }=bJ=0$. This constraint directly implies that GB motion is described by ${h}_{t}=-HJ=-\frac{1}{\beta }bJ=0$; in other words, GBs in polycrystals cannot migrate with a single disconnection mode. This conclusion is consistent with the observation of GB stagnation in bicrystal atomistic simulations¹² and the observations that most GBs do not move at low temperature, where only one disconnection mode is activated. The internal stress on the GB can be obtained from the constraint,

$$\dot{\epsilon }={M}_{{\rm{d}}}{b}^{2}\left[\tau +(\Psi -\gamma {h}_{xx})\frac{1}{\beta }\right](| \rho | +2c)=0,$$

(21)

which is satisfied only when βτ + (Ψ − γh_xx) = 0. When the GB is flat, we have h_xx = 0 and τ = −Ψ/β, which implies that if a chemical potential-jump Ψ is applied to the GB, its migration will generate a back stress τ. Their relationship τ = −Ψ/β is in agreement with molecular dynamics simulations¹².

The character of different disconnection modes active on a GB influences GB dynamics. Figure 5 compares the GB evolution when there are two and three disconnection modes active. At T = 800 K, the GB profiles evolve almost identically (Fig. 5a). On the other hand, at the relatively high temperature of T = 1200 K, the GB profile evolves more quickly when three modes are included (relatively to the two-mode case) (Fig. 5b). The disconnection formation energies increase and formation rates decrease with increasing mode number^13,15,36,37. The thermally nucleated disconnections have a non-negligible effect at high T, while disconnection migration dominates GB evolution at low T.

Fig. 5: Comparison of the GB evolution with different active disconnection modes.

Comparison of the GB profile evolution in the case of two and three disconnection modes. The arrows indicate the direction of the evolving curves. The three disconnections modes share the same Burgers vector ${b}^{(i)}=b={a}_{0}/\sqrt{26}$, and H⁽¹⁾ = 4b, H⁽²⁾ = −2.5b, H⁽³⁾ = − 9b. The initial disconnection densities for the first two modes are as per (19), while in the three mode calculations in a, b, we also set the initial values of ρ⁽³⁾ = 0. In c, the initial disconnection densities are as per (15) and (16) with ρ⁽¹⁾:ρ⁽²⁾:ρ⁽³⁾ = 1:−0.65:−0.35.

We also compare the evolution of the GB profile and disconnection densities in the two and three disconnection mode cases for arbitrary initial disconnection densities (see Fig. 5c and the Supplementary Information for details). Again, we see that the observation that GB migration is faster when more disconnections modes are included at high temperature is a general phenomenon. We reiterate that with a single disconnection mode GB migration does not occur, at low T GB migration with two or three (or more) modes are nearly identical and adding more modes increases the migration rate at high T.

The present paper presents a reduced-order model for GB motion in a two-dimensional continuum dynamics framework based upon the underlying (discrete) disconnection dynamics mechanism; our formulation is computationally efficient and appropriate for large-scale simulations. Here, we described the model and showed examples specifically for the case of GB migration in a two-dimensional system. This two-dimensional model can be extended to three dimensions, appropriate for application to polycrystalline microstructures in typical materials. While the main ideas and model present here easily translate to three dimensions, there are some important differences. First and foremost, disconnections in three dimensions are lines rather than points and hence the model must include disconnection topology and its evolution. Second, disconnection nucleation in three dimensions includes a wider range of option than considered here in two dimensions; including disconnection bowing, disconnection loop formation, etc. The three-dimensional continuum framework for the dynamics of low angle grain boundaries incorporating dislocation structures^38,39 can be adopted. Alternatively, the three-dimensional model can also be implemented in a discrete dislocation dynamics-type of framework, albeit constrained to a plane (similar to dislocation dynamics where dislocations are constrained to a slip plane; see e.g. refs. ^{40,41,42,43,44}). The constraint of disconnections lines to a GB plane rather than dislocation lines in three dimensions means that such calculations need only be codimension 1 rather than 2 in the three-dimensional dislocation dynamics case; hence such three-dimensional simulations are much easier to implement than full, three-dimensional dislocation dynamics simulations. Such discrete dislocation-dynamic like calculations may be implemented in either sharp or diffuse dislocation/disconnection dynamics formulations.

Methods

Derivation of disconnection velocity

The disconnection glide velocity expression (5) is applicable to the case where the motion is overdamped; i.e. the velocity is proportional to the driving force ${f}_{{\rm{d}}}^{(j)}$, i.e., ${v}^{(j)}={M}_{{\rm{d}}}^{(j)}\cdot {f}_{{\rm{d}}}^{(j)}$, where ${M}_{{\rm{d}}}^{(j)}$ is the disconnection mobility. The driving force on a jth mode disconnection is

$${f}_{{\rm{d}}}^{(j)}=\tau {b}^{(j)}+\Psi {H}^{(j)}-\gamma {h}_{xx}{H}^{(j)},$$

(22)

which has contributions from the stress τ (both applied stress and that originating from all other disconnections), the jump in the energy/chemical potential across the GB Ψ, and capillarity where γ is the GB energy density. Positive/negative disconnection density ρ^(j) implies disconnections of type (b^(j), H^(j))/(−b^(j), −H^(j)) with velocities v^(j)/ − v^(j). The disconnection flux is J^(j) = v^(j)(∣ρ^(j)∣ + c^(j)), as per (4).

Numerical methods of the grain boundary evolution

In the numerical simulations, we employ a forward Euler and finite difference methods in time and space. The GB length L₀ is discretised into N equally spaced intervals [x_i−1, x_i], where x_i = x₀ + iΔx, Δx = L₀/N, i = 0, 1, … , N. h(x) and ρ(x) are evaluated at spatial discrete points and h_x is calculated using a backward difference scheme and h_xx using a centre difference scheme. ∣ρ^(j)∣ in ${A}^{(j)}={M}_{{\rm{d}}}^{(j)}{{b}^{(j)}}^{2}(| {\rho }^{(j)}| +2{c}^{(j)})$ is calculated using an upwind scheme. The discretised temporal points are t_n = nΔt, n = 0, 1, 2, …. The numerical value of a function g(x, t) at the discrete spacial point x_i and temporal point t_n is denoted as ${g}_{i}^{n}$.

Specifically, the GB evolution in Eq. (20) is numerically implemented as

$${h}_{i}^{n+1}=\left\{\begin{array}{ll}{h}_{i}^{n}-{{\Delta }}t\cdot B\left(\Psi -\gamma \frac{{h}_{i+1}^{n}-{h}_{i}^{n}+{h}_{i-1}^{n}}{{({{\Delta }}x)}^{2}}\right)\left(\left|\frac{{h}_{i}^{n}-{h}_{i-1}^{n}}{{{\Delta }}x}\right|+D\right),&\,\text{if}\,B\left(\Psi -\gamma \frac{{h}_{i+1}^{n}-{h}_{i}^{n}+{h}_{i-1}^{n}}{{({{\Delta }}x)}^{2}}\right)\left(\frac{{h}_{i+1}^{n}-{h}_{i-1}^{n}}{2{{\Delta }}x}\right)>\,0\,\\ {h}_{i}^{n}-{{\Delta }}t\cdot B\left(\Psi -\gamma \frac{{h}_{i+1}^{n}-{h}_{i}^{n}+{h}_{i-1}^{n}}{{({{\Delta }}x)}^{2}}\right)\left(\left|\frac{{h}_{i+1}^{n}-{h}_{i}^{n}}{{{\Delta }}x}\right|+D\right),&\,\text{if}\,B\left(\Psi -\gamma \frac{{h}_{i+1}^{n}-{h}_{i}^{n}+{h}_{i-1}^{n}}{{({{\Delta }}x)}^{2}}\right)\left(\frac{{h}_{i+1}^{n}-{h}_{i-1}^{n}}{2{{\Delta }}x}\right)<\,0\,\end{array}.\right.$$

(23)

The GB profile and disconnection density evolution equations (10) and (11) are numerically implemented as

$${h}_{i}^{n+1}={h}_{i}^{n}-{{\Delta }}t\left(\Psi -\gamma \frac{{h}_{i+1}^{n}-{h}_{i}^{n}+{h}_{i-1}^{n}}{{({{\Delta }}x)}^{2}}\right)\frac{{\sum }_{k \,{>}\,j}{\left(\frac{1}{{\beta }^{(j)}}-\frac{1}{{\beta }^{(k)}}\right)}^{2}{A}^{(j)}{| }_{i}^{n}\cdot {A}^{(k)}{| }_{i}^{n}}{{\sum }_{k}{A}^{(k)}{| }_{i}^{n}},$$

(24)

$${\rho }_{{\rm{back}}}^{(j)}{| }_{i}^{n+1}={\rho }_{{\rm{back}}}^{(j)}{| }_{i}^{n}-{{\Delta }}t\frac{{J}^{(j)}{| }_{i}^{n}-{J}^{(j)}{| }_{i-1}^{n}}{{({{\Delta }}x)}^{2}},$$

(25)

$${\rho }_{{\rm{forw}}}^{(j)}{| }_{i}^{n+1}={\rho }_{{\rm{forw}}}^{(j)}{| }_{i}^{n}-{{\Delta }}t\frac{{J}^{(j)}{| }_{i+1}^{n}-{J}^{(j)}{| }_{i}^{n}}{{({{\Delta }}x)}^{2}},$$

(26)

$${\tau }_{i}^{n+1}=-\left(\Psi -\gamma \frac{{h}_{i+1}^{n}-{h}_{i}^{n}+{h}_{i-1}^{n}}{{({{\Delta }}x)}^{2}}\right)\frac{{\sum }_{j}\frac{1}{{\beta }^{(j)}}{A}^{(j)}{| }_{i}^{n}}{{\sum }_{j}{A}^{(j)}{| }_{i}^{n}};$$

(27)

where

$${J}^{(j)}{| }_{i}^{n}=\left(\Psi -\gamma \frac{{h}_{i+1}^{n}-{h}_{i}^{n}+{h}_{i-1}^{n}}{{({{\Delta }}x)}^{2}}\right)\frac{1}{{b}^{(j)}}\cdot \frac{{\sum }_{k}\left(\frac{1}{{\beta }^{(j)}}-\frac{1}{{\beta }^{(k)}}\right){A}^{(j)}{| }_{i}^{n}\cdot {A}^{(k)}{| }_{i}^{n}}{{\sum }_{k}{A}^{(k)}{| }_{i}^{n}},$$

(28)

$${A}^{(j)}{| }_{i}^{n}=\left\{\begin{array}{ll}{M}_{{\rm{d}}}^{(j)}{{b}^{(j)}}^{2}(| {\rho }_{{\rm{back}}}^{(j)}{| }_{i}^{n}+2{c}^{(j)}),&\,\text{if}\,{v}^{(j)}{| }_{i}^{n}({\rho }_{{\rm{back}}}^{(j)}{| }_{i}^{n}+{\rho }_{{\rm{forw}}}^{(j)}{| }_{i}^{n})>\,0\,\\ {M}_{{\rm{d}}}^{(j)}{{b}^{(j)}}^{2}(| {\rho }_{{\rm{forw}}}^{(j)}{| }_{i}^{n}+2{c}^{(j)}),&\,\text{if}\,{v}^{(j)}{| }_{i}^{n}({\rho }_{{\rm{back}}}^{(j)}{| }_{i}^{n}+{\rho }_{{\rm{forw}}}^{(j)}{| }_{i}^{n})<\,0\,\end{array},\right.$$

(29)

$${v}^{(j)}{| }_{i}^{n}={M}_{{\rm{d}}}^{(j)}\left[{\tau }_{i}^{n}{b}^{(j)}+\Psi {H}^{(j)}-\gamma \frac{{h}_{i+1}^{n}-{h}_{i}^{n}+{h}_{i-1}^{n}}{{({{\Delta }}x)}^{2}}{H}^{(j)}\right].$$

(30)

Periodic boundary conditions are applied at the endpoints of the GB.