Polycrystal plasticity with grain boundary evolution: A numerically efficient dislocation-based diffuse-interface model

A body is represented by an open subset $B$ of a Euclidean space $\mathbb{R}^{3}$. A time-dependent motion of $B$, described relative to a reference configuration $B_{0}\subset\mathbb{R}^{3}$, is given by a smooth one-to-one map $\bm{x}(\bm{X},t)$, which maps a material point $\bm{X}$ to its spacial position $\bm{x}$ at time $t$. Let $\bm{F}(\bm{X},t):=\nabla\bm{x}(\bm{X},t)$ denote the gradient of the deformation map at time $t$, where $\nabla$ denotes the gradient with respect to the material coordinate $\bm{X}$.

In finite deformation plasticity, the deformation gradient admits a multiplicative decomposition (Hill, 1966; Lee, 1969; Reina and Conti, 2014)

where $\bm{F}^{\rm L}$ and $\bm{F}^{\rm P}$ denote the lattice¹¹1Also known as the elastic distortion, $\bm{F}^{\rm E}$. The notation $\bm{F}^{\rm L}$ is adopted from Clayton (2010). and plastic distortions, respectively. $\bm{F}^{\rm P}$ maps an infinitesimal line element $d\bm{X}$ to $\bm{F}^{\rm P}d\bm{X}$. If we regard the collection of line elements $\bm{F}^{\rm P}d\bm{X}$ as the lattice configuration, then $\bm{F}^{\rm P}$ can be understood as a linear transformation from the reference configuration to the lattice configuration, whereas $\bm{F}^{\rm L}$ is a linear mapping from the lattice configuration to the deformed configuration (Admal et al., 2018).²²2Strictly speaking, $\bm{F}$, $\bm{F}^{\rm L}$ and $\bm{F}^{\rm P}$ are linear transformations between the tangent spaces of the respective configurations. In this paper, the body $B$ represents a polycrystalline material, and we limit the mechanism of plastic deformation to dislocation slip on the local crystal’s slip systems. Therefore, assuming the crystal has $N_{s}$ slip systems, the evolution of $\bm{F}^{\rm P}$ is given by

where $\bm{L}^{\rm P}$ denotes the plastic velocity gradient given in terms of the Schmid tensors $\bm{S}^{\alpha}$ and the slip rates $\nu^{\alpha}$ ($\alpha=1,\dots,N_{s}$). The Schmid tensor $\bm{S}^{\alpha}:=\bm{s}^{\alpha}\otimes\bm{m}^{\alpha}$ is defined in terms of the $\alpha$-th slip direction and slip plane normal, $\bm{s}^{\alpha}$ and $\bm{m}^{\alpha}$, respectively.

In order to model evolving grain boundaries within a crystal plasticity framework, Admal et al. (2018) proposed a construction of diffuse-interface grain boundaries using GNDs, followed by equipping their model with a GB energy density in terms of the GND density. We will now describe their construction of GNDs that results in a stress-free polycrystal.

Traditional crystal plasticity theories describe the deformation of a polycrystal using a strain-free polycrystal as a reference configuration, along with an initial decomposition ($\bm{F}=\bm{F}^{\rm E}\bm{F}^{\rm P}$) of the deformation gradient given by³³3 Notice that the notation $\bm{F}^{\rm E}$ is used in Eq. (3) instead of $\bm{F}^{\rm L}$. This is because $\bm{F}^{\rm E}$ describes the elastic distortion of the lattice relative to the reference polycrystal, and does not include the orientation of the grain.

Since the reference configuration is a polycrystal, the Schmid tensors, defined in Section 2.1, are piecewise-constant and depend on the orientation of each grain. While the multiplicative decomposition of $\bm{F}=\bm{I}$ in Eq. (3) is a reasonable choice for modeling polycrystal plasticity with fixed grain boundaries, it does not lend itself to a construction of GB energy and modeling evolving GBs.

An alternate approach is to choose a strain-free single crystal as the reference configuration, along with an initial decomposition

where $\bm{R}^{0}(\bm{X})\in SO(3)$ is a piecewise-constant rotation field indicating the crystal orientation of each grain. By construction, $\bm{F}(\bm{X},0)\equiv\bm{I}$, and the initial Green–Lagrangian lattice strain $\bm{E}^{\rm L}(\bm{X},0)=((\bm{F}^{\rm L})^{T}\bm{F}^{\rm L}-\bm{I})/2=\bm{0}$, which implies the resulting polycrystal is strain-free. At the atomistic scale, it is well known that the dislocation cores of the low-angle GBs result in lattice distortion and stress fields, hence the GB is not strain free. However, the stress fields are short-ranged, and the associated elastic strain energy is already accounted for in the form of GB energy. Hence, in the continuum scale, we model GBs with impotent GNDs that do not generate any stresses, which results in a Read–Shockley type GB energy. Moreover, since the reference configuration is a single crystal, the Schmid tensors are now constant, and correspond to the slip systems of the single crystal. For numerical convenience, the piecewise constant field $\bm{R}^{0}(\bm{X})$ is replaced with a smoothed version $\tilde{\bm{R}}^{0}(\bm{X})\in SO(3)$, leading to a diffuse-interface model.

The main advantage of the decomposition given in Eq. (4) is that GBs are characterized by the GND density, defined as

where Curl denotes the curl⁴⁴4In indicial notation, $({\rm{Curl}}\bm{A})_{ij}=\epsilon_{ipq}A_{jq,p}$. of a tensor field in material coordinates. The GND density tensor can then be used to construct the grain boundary contribution of the free energy functional, as shown in the next section.

In this section, we begin by describing the constitutive law proposed by Admal et al. (2018) that results in a simultaneous evolution of grain boundary mediated phenomena and bulk plasticity. The construction of the constitutive law for the GB energy is based on the assumption that the GND density that constitutes the GB quantifies its energy. Therefore, the total energy density $\psi$ is given by a sum of a bulk elastic energy density $\psi^{\rm e}(\bm{E}^{\rm L})$, and a GB energy density $\psi^{\rm gb}$, i.e. $\psi=\psi^{\rm e}+\psi^{\rm gb}$. For simplicity, we limit ourselves to isothermal condition.

The elastic energy density is assumed to be of the classical form that depends only on the lattice strain $\bm{E}^{\rm L}$. On the other hand, the construction of the GB energy density is inspired by the KWC model (Kobayashi et al., 1998), a diffuse-interface model for grain microstructure evolution with two order parameters $\phi$ and $\theta$ that describe the degree of crystallinity and crystal orientation, respectively. The GB energy density of the KWC model is given by

where

The parameters $\alpha$, $e_{0}$, $s$, and $\epsilon$ are material parameters that determine the GB energy and width. The construction of $\psi^{\rm KWC}$ by Kobayashi et al. (2000) was centered around the observation that the term $sg(\phi)|\nabla\theta|$ acts to localize the grain boundary, while the term $\frac{\epsilon^{2}}{2}|\nabla\theta|^{2}$ tends to diffuse it. The combined effect of the two gives rise to GBs of finite widths. Moreover, the inclusion of the $|\nabla\theta|^{2}$ term results in motion by curvature. Noting that GNDs measure the gradient of lattice rotation, and the motion of GBs is given by the evolution of GNDs, Admal et al. (2018) have shown that by replacing $\nabla\theta$ with $\bm{G}$, and defining $\psi^{\rm gb}$ as

where $|\bm{G}|$ denotes the standard Euclidean norm of the GND density tensor, the total energy

is a functional of $\phi$, the displacement field $\bm{u}:=\bm{x}(\bm{X},t)-\bm{X}$, and $\bm{F}^{\rm P}$. The resulting model not only inherits all the features of the KWC model, such as GB motion by curvature and grain rotation, but it also describes additional phenomena, such as coupled GB motion, sliding and recovery. This is because the GNDs constituting the GBs respond to stress,⁵⁵5Currently, we limit this response to dislocation glide. and the shear accompanying their glide results in shear-induced grain boundary motion. It is important to recognize that, depending on the availability of slip systems and the loading condition, the slip rates in a GB may also result in a combination of sliding and pure coupling.

We now extend the constitutive law given in Eq. (6) to include energy due to SSDs, which play an important role in the evolution of grain microstructure of cold-formed polycrystals. Our construction of the energy due to SSDs is inspired by the synthetic potential energy for atomistic systems proposed by Janssens et al. (2006), where the motivation was to induce GB motion in bicrystals in the absence of curvature.⁶⁶6Although a flat grain boundary can be driven by shear stress, this approach is not preferred in MD simulations because the limited times scales accessible to MD necessitate the use of extremely high stress and strain rates that do not represent commonly observed deformation rates. Selectively adding potential energy to atoms on one side of the GB, results in an artificial force on the atoms with higher potential which reorients them with atoms on the other side of the GB. We will now present a continuum analog of the synthetic potential, which represents energy due to SSDs. We limit the description to symmetric tilt grain boundaries in two dimensions.

A synthetic potential is incorporated into $\psi^{\rm{gb}}$ by modifying the functional form of $f(\phi)$ in Eq. (7) to

where $\rho$ denotes a user-defined scaling parameter, and $\theta^{\rm{L}}$ denotes the current crystal orientation field, which can be obtained from the lattice rotation $\bm{R}^{\rm L}$ computed using the right polar decomposition of $\bm{F}^{\rm L}$. The term $[1-{\rm{tanh}}(\rho\theta^{\rm{L}})]$ can be interpreted as a grain selector function, which serves to allocate an additional energy of $e_{0}\Delta E^{*}$ to the selected grain (in this case, the grain with a negative $\theta^{\rm{L}}$). In Section 3.2, we demonstrate a numerical implementation of a GB driven by synthetic potential.

One of the main goals of this paper is to formulate a computationally tractable model that can model GB and bulk plasticity. In this section, we derive evolution equations for the unknown fields — displacement $\bm{u}:=\bm{x}-\bm{X}$, plastic distortion $\bm{F}^{\rm P}$, and the order parameter $\phi$ — that are computationally more tractable than those derived by Admal et al. (2018). In Section 3, we demonstrate the computational efficiency of our formulation.

Our derivation is based on the virtual work formulation of Gurtin (2008). The principle of virtual work states that the power expended on an arbitrary part of a body by external forces is equal to the internal power expended within that part by internal forces. We assume that the internal power, denoted by $I(P)$, is expended on a part $P\subset B_{0}$ by a stress tensor $\bm{P}$ power conjugate to $\dot{\bm{F}}$, a stress vector $\bm{p}$ conjugate to $\nabla\dot{\phi}$, microstress vectors $\bm{\xi}^{\alpha}$ conjugate to $\nabla\nu^{\alpha}$, and internal forces⁷⁷7Internal forces contribute neither to the work nor to the second law. $\pi$ and $\Pi^{\alpha}$ conjugate to $\dot{\phi}$ and slip rate $\nu^{\alpha}$, respectively. Therefore, $I(P)$ is given by

Using the divergence theorem, Eq. (11) transforms to

where $\bm{N}$ denotes the outward unit normal to $\partial P$. The boundary integral in Eq. (12) suggests the following form for the power expended on $P$ by external forces $\bm{t}$, $s$ and $\bm{\Xi}$ that are power conjugate to $\dot{\bm{y}}$, $\dot{\phi}$ and $v$, respectively:

Enforcing $I(P)=E(P)$ for all $P\subset B_{0}$, we arrive at the following balance equations

• •

  The standard force balance equation for $\bm{u}$

  $$\begin{array}[]{l}{\rm{Div}}\,\bm{P}=\bm{0}\quad{\rm in}\;B_{0},\\ \bm{t}=\bm{P}\bm{N}\;\;\;\;\;\;{\rm on}\;\partial B_{0}.\end{array}$$

  (14)

• •

  A microscopic force balance for $\phi$

  $$\begin{array}[]{l}{\rm{Div}}\,\bm{p}-\pi=0\;\;\;\;{\rm in}\;B_{0},\\ s=\bm{p}\cdot\bm{N}\;\;\;\;\;\;\;\;\;\;{\rm on}\;\partial B_{0}.\end{array}$$

  (15)

• •

  A microscopic force balance on each slip system for $\nu^{\alpha}$

  $$\begin{array}[]{l}{\rm{Div}}\,\bm{\xi}^{\alpha}-\Pi^{\alpha}=0\;\;\;{\rm in}\;B_{0},\\ \Xi^{\alpha}=\bm{\xi}^{\alpha}\cdot\bm{N}\;\;\;\;\;\;\;\;\;{\rm on}\;\partial B_{0}.\end{array}$$

  (16)

Next, we obtain the following constitutive relations between stresses and kinematic fields by examining the first and second laws of thermodynamics using the Coleman–Noll procedure (see A)

where we split the microstress $\bm{\xi}^{\alpha}$ into an energetic part $\bm{\xi}^{\alpha}_{\rm e}$ and a dissipative part $\bm{\xi}^{\alpha}_{\rm d}$, and $B^{\alpha}>0$ is the inverse mobility for $\nabla v^{\alpha}$. The internal forces $\pi$ and $\Pi^{\alpha}$ are given by

is the resolved shear stress on the $\alpha$th slip system, while $b^{\phi}>0$ and $b^{\alpha}>0$ denote the inverse mobilities for $\phi$ and $\nu^{\alpha}$, respectively.⁸⁸8In Eq. (21), we are using the notation $\bm{A}:\bm{B}$ to denote the inner product between tensors $\bm{A}$ and $\bm{B}$. In indicial notation, $\bm{A}:\bm{B}=A_{ij}B_{ij}$. To conclude, equations (2), (14), (15) and (16) are the governing equations for the unknown fields $\bm{F}^{\rm p}$, $\phi$, $\bm{u}$ and $\nu^{\alpha}$, respectively.

We note that Eqs. (14)–(22) are identical to those derived by Admal et al. (2018), wherein the primary unknowns in Eq. (16) are the slip rates $\nu^{\alpha}$, which have to be solved for each corresponding slip system. This is a major disadvantage, as correctly pointed out by Ask et al. (2019), in materials with a large number of slip systems — for example, body-centered-cubic materials like $\alpha$-Fe have up to 48 slip systems. To this end, we note that the slip rate on a specific slip system is seldom of interest since the motion of GBs, through the evolution of $\bm{F}^{\rm P}$ via Eq. (2), is a combined effect of slip on all available slip systems. This motivates us to seek an alternative formulation, which captures the combined effects of all slip rates, while avoiding the need to solve for each slip rate explicitly. To this end, we substitute Eq. (21) into Eq. (16) and rearrange as

From Eq. (2), $\bm{L}^{\rm P}$ can be expressed as

Substituting Eq. (23) into Eq. (24), and rearranging, we have

Before proceeding to simplifying Eq. (25), we consider explicit functional forms for the inverse mobilities $B^{\alpha}$ and $b^{\alpha}$ appearing in Eq. (19) and Eq. (21) respectively.

In this work, since we focus primarily on GB plasticity, we assume inverse mobilities $b^{\alpha}$s are of the form

where $\overline{b}$ denotes a nominal inverse mobility, and is defined as

The quantities $m_{\rm min}$ and $m_{\rm max}$ are the minimum and maximum mobilities, respectively. In addition, we assume $B^{\alpha}=c_{B}b^{\alpha}$, where the scaling coefficient $c_{B}$ is a constant with dimension of length. The functional form of $\overline{b}$ in Eq. (27) ensures $\overline{b}=m_{\rm min}^{-1}$ in the grain interiors where $\phi\approx 1$, and by choosing $m_{\rm min}^{-1}\gg m_{\rm max}^{-1}$, plastic slip in the grain interiors is suppressed. The factor $c^{\alpha}$ is a dimensionless slip system-dependent availability factor — $c^{\alpha}=0$ disables the $\alpha^{th}$ slip system, while $c^{\alpha}=1$ ensures the slip system is available with inverse mobility $\overline{b}$.

Recalling the definitions of the energetic and dissipative components of $\bm{\xi}^{\alpha}$ introduced in Eq. (19), we note that except for the dissipative part of the left-hand-side of Eq. (25), i.e.

all terms in Eq. (19) are in terms of $\bm{F}^{\rm P}$. Therefore, we now draw our attention to (28). Using Eq. (24), the definition of $B^{\alpha}$ described above, and noting that $\bm{S}^{\alpha}$ is a constant tensor, the expression in (28) simplifies as

where $\nabla\bm{A}\bullet\bm{a}$ is a second-order tensor whose $ij$-component is defined as $A_{ij,k}a_{k}$.⁹⁹9We have used the identity $\mathrm{Div}(a\bm{A})=a\mathrm{Div}\bm{A}+\bm{A}\nabla a$ for some scalar field $a$ and tensor field $\bm{A}$ to arrive at Eq. (29) Finally, substituting Eq. (29) into Eq. (25), we obtain

which is a system of equations for $\bm{F}^{\rm P}$ that does not involve any slip rates. The transformation of Eq. (16) into Eq. (30) gives rise to a new formulation wherein the unknown fields are $\phi$, $\bm{u}$ and $\bm{F}^{\rm P}$, and their governing equations are (14), (15) and (30), respectively. In other words, slip rates are no longer solved for in the new formulation. Compared to Eq. (16), which is a system of $N_{s}$ equations, Eq. (30) is a system of $d^{2}$ equations — where $d$ denotes the space dimension — our formulation results in a computational advantage whenever $d^{2}<N_{s}$. As noted earlier, this is indeed the case for most crystalline materials in three dimensions. In addition, the boundary conditions (see Eq. (16)) for $\nu^{\alpha}$, in the formulation by Admal et al. (2018), transform into Dirichlet/Neumann boundary conditions for $\bm{F}^{\rm P}$. Finally, we note that although Eq. (16) and Eq. (30) are not equivalent — satisfying Eq. (30) is a necessary but not sufficient condition for satisfying Eq. (16) on each slip system — the two equations should, in principle, result in identical evolution for $\bm{F}^{\rm P}$.

In this example, we simulate grain boundary motion induced by an applied synthetic driving force. In addition, we implemented the same synthetic driving force in the formulation by Admal et al. (2018), and compare the results of the two formulations.

For the model problem, we consider a bicrystal with a size 20 by 3.33 nm² centered at $\bm{X}=(10\,\mathrm{nm},10\,\mathrm{nm})$ with a symmetric tilt grain boundary at $X_{1}=10\,\mathrm{nm}$. The initial condition for the crystal orientation is given by $\bm{F}^{\rm L}(\bm{X},0)=\bm{R}(\theta^{0}(\bm{X}))=(\bm{F}^{\rm P}(\bm{X},0))^{T}$, where

describes a GB with misorientation of $30^{\circ}$. The bicrystal is fixed on its left end but is free to deform on the right end, yielding a free-end condition. Periodic boundary conditions are applied to the top and bottom faces of the computational domain to mimic an infinite bicrystal. We assume the material is equipped with $12$ slip systems with equal inverse mobilities ($b^{\alpha}$s). To construct the slip systems, 12 uniformly spaced angles $\beta^{\alpha}$ in the range $[0,90^{\circ}]$ are generated, and the corresponding slip directions and normals are computed as

The GB is driven to the right by adding a synthetic potential of strength $\Delta E=0.147\,\rm{GPa}$ to the right-hand grain using $f(\phi)$ given in Eq. (10). The model parameters used in the simulations discussed in this section are listed in Table 1.

The unknown fields $\phi$, $\bm{u}$, and $\bm{F}^{\rm P}$ are interpolated using Lagrange quadratic finite elements. At this point, it is emphasized that the initial condition for $\bm{F}^{\rm P}$ is a smooth rotation field satisfying the orthogonality condition $(\bm{F}^{\rm P})^{T}\bm{F}^{\rm P}=\bm{I}$ everywhere in the domain. Such is not the case if $\bm{F}^{\rm P}$ is interpolated component-wise using finite element shape functions. To ensure that the interpolation of initial $\bm{F}^{\rm P}$ is an orthogonal tensor field, we express $\bm{F}^{\rm P}$ using its right polar decomposition $\bm{F}^{\rm P}=\bm{R}^{\rm P}\bm{U}^{\rm P}$, and in 2D, the four components of $\bm{F}^{\rm P}$ are replaced by $\theta^{\rm P}$, $U^{\rm P}_{11}$, $U^{\rm P}_{12}$ and $U^{\rm P}_{22}$ which are interpolated using Lagrange quadratic finite elements. This treatment is adopted in all simulations. When solving the equations from both the current formulation and the framework by Admal et al. (2018), a time-dependent solver with a constant step size of $\Delta t=1\times 10^{-7}$s was used for time stepping, and both systems were evolved for $5\times 10^{-6}$s.

First, the effect of the synthetic potential on the order parameter $\phi$ is discussed. The contour of $\phi$ at the end of the simulation, extracted from the current framework, is plotted in Fig. 1a, along with that obtained from a simulation with no added potential. A time history plot of the value of $\phi$ on the right-hand grain is shown in Fig. 1b.

From Fig. 1a, we see that minimum of $\phi$ has clearly shifted to the right, indicating rightward GB migration, as expected. We also notice that the added potential manifests as an asymmetry in $\phi$ with $\phi<1$ in the right grain. As the GB migrates, the value of $\phi$ in the right grain decreases gradually and converges to its steady-state value, as evidenced in Fig. 1b. The lower value of $\phi$ in the right grain is a consequence of the synthetic potential, interpreted as energy due to SSDs, added to the right grain. While at the continuum scale SSDs do not result in lattice distortions, at the atomic scale, they contribute to disorder. Since $\phi<1$ describes disorder, the prediction that $\phi<1$ in the right grain is in line with the physical interpretation of $\phi$. The synthetic potential constructed in this paper allows us to induce migration of GBs surrounding a particular grain, identified by its crystal orientation, and can be used as a surrogate for many forms of driving forces due to energy differences, either coming from differences in SSD densities, stored strain energy, or plastic anisotropy.

Next, we leverage this example to compare the results from the current formulation to that by Admal et al. (2018). The time history plots of the GB displacements and measured coupling factor are shown in Fig. 2.

As evidenced by Fig. 2a, some differences in GB displacements are visible, but they remain below 5%. The fitted GB migration velocities reveal a 6.07% difference. However, the simulated coupling factor show good agreement, as shown in Fig. 2b. These observed differences are within acceptable ranges, and the difference in migration velocity can be compensated for in model parameter calibration. In Section 2.4, it was hypothesized that the solutions to Eq. (30) and Eq. (16) should yield identical evolutions for $\bm{F}^{\rm P}$, we now check the validity of this hypothesis with simulation data. Recalling that GB plasticity (which is what $\bm{F}^{\rm P}$ measures) is due to GB migration, we plot the mean percent difference in components of $\bm{F}^{\rm P}$ as a function of the GB displacement, see Fig. 2c. The plot of the mean differences in components of $\bm{F}^{\rm P}$ shows good agreement between the evolution histories computed from the two formulations, thus confirming our hypothesis.

In this example, we compare the computational cost of the current formulation to that by Admal et al. (2018) using a similar model problem as discussed in the previous section, where the scaling of total number of DOFs and simulation time are measured. We limit our comparison to the framework by Admal et al. (2018) since it is closely related to our current work and includes identical physics. It is true that many other models for GB migration exist, but a comparison of numerical efficiency is not fair unless the models being compared include similar physics. To this end, it is appropriate to compare our work with some dislocation-based models such as Zhang and Xiang (2018) and Zhang and Xiang (2020), and disconnection-based models such as Runnels and Agrawal (2020) and Gokuli and Runnels (2021). However, as those models are still under active development and are not formulated for numerical efficiency, further comparison is not pursued here.

In order to explore the scaling of our formulation, we change the size of the bicrystal to a square with an edge length of 20 nm and consider a sequence of uniform meshes with $N^{2}$ elements, generated using the sequence $N=\{25,50,100,110,150,200\}$. All other simulation settings are identical to that described in Section 3.1. A plot comparing the number of DOFs in each mesh, and the corresponding simulation time is shown in Fig. 3.

The reduction in the number of DOFs, and its dependence on $N$ is shown in Fig. 3a. In this example, the framework in Admal et al. (2018) requires a total of 19 DOFs (1 for $\phi$, 2 for $\bm{u}$, 12 for $\nu_{\alpha}$ and 4 for $\bm{F}^{\rm P}$) per node. While using the current formulation, the number of DOFs per node is reduced to 7 (no DOFs related to $\nu_{\alpha}$), which is a 63.2% decrease. The reduced number of DOFs in turns leads to a huge decrease in simulation time, as evidenced in Fig. 3b. Apart from the direct advantage that the new formulation yields fewer number of DOFs for a given mesh, it is interesting to note that even when the number of DOFs are identical (e.g., two blue points in Fig. 3a), the new formulation has a shorter run time (e.g., comparing the two blue points in Fig. 3b, shows that the simulation time of the formulation by Admal et al. (2018) is about 6.6 times longer). We hypothesize that this additional gain in performance is due to the change in sparsity pattern (hence the bandwidth) of the system matrix associated with fewer DOFs per node, as the solution time for sparse matrix solvers depend both on the system size and its sparsity pattern. It should be noted that although the scaling test was performed on a simplified geometry, similar performance enhancements are expected in general. The reason for this claim is that the performance gain roots from the reduced number of DOFs per node, which is independent of the underlying geometry or the polycrystal system. This significant increase in computational efficiency, resulting from fewer DOFs per node and improved sparsity pattern, is beneficial in many ways. For more complex crystal configurations, a refined mesh near the GBs is desired to resolve high local gradients. The current formulation relaxes the requirements on computational cost, giving room for the use of a more locally refined mesh. In addition, reduced number of DOFs entails less memory usage during the simulations.

In this example, we study the evolution of a copper tricrystal, originally constructed by Trautt and Mishin (2014), that is formed by embedding a circular grain into a bicrystal, as shown in Fig. 4. In particular, we make a quantitative comparison with predictions from MD, and explore the effect of shear stress on the GB motion and grain rotation. All results are are compared to the MD simulations by Trautt and Mishin (2014).

The three grains are oriented such that their $[001]$ crystal directions are all aligned along the $z$-axis pointing vertically out of the page. A circular grain of radius 8.33 nm with an orientation of 0^∘ is placed in the center of a square bicrystal of side length 30 nm, whose upper and lower grains have orientations of 73.74^∘ and 16.26^∘, respectively. The $[100]$ crystal direction for the center grain is parallel to the $x$-axis. The computational domain is discretized with quadratic Lagrange elements. This configuration corresponds to Case 1 in the MD study by Trautt and Mishin (2014). We note that due to the four-fold symmetry of the cubic crystal, the orientations of the grains are not uniquely defined, and this manifests as non-uniqueness of the GND density. For example, the orientations described above result in a larger GND density for the upper semi-circular grain boundary relative to its lower counterpart, while changing the orientation of the top grain to an equivalent $-16.26^{\circ}$ results in equal misorientations and GND densities. Since the two semi-circular GBs have identical atomistic structures, and therefore equal energies, we proceed with the latter choice. The tricrystal is equipped with three slip systems with equal inverse mobilities, and oriented along directions given by Eq. (32) with $\beta^{\alpha}=0^{\circ}$, $45^{\circ}$, and $-45^{\circ}$. To match the initial grain shrinkage rate as observed in MD by Trautt and Mishin (2014) (cf. Figure 8b therein), the inverse mobility for $\phi$, $b^{\phi}$, and the maximum slip mobility, $m_{\rm max}$, were calibrated as $b^{\phi}=4.165\times 10^{8}\;\;\rm{kg}/(m\cdot s)$ and $m_{\rm max}=1.2\times 10^{-7}\rm{(m\cdot s)}/kg$, respectively. All other simulation parameters are identical to those listed in Table 1, except that no synthetic potential is applied.

A schematic showing the initial crystal orientation and the applied boundary conditions is shown in Fig. 4. Two cases are considered, one with no external loading, and the other with an external shear of 500 MPa acting on the top surface, and the system is evolved for 19 ns and 14 ns, respectively. Selected snapshots of crystal orientations during their evolution for the two load cases are shown in Fig. 5 and Fig. 7, where the latter corresponds to the case of applied shear. Fig. 6 compares the variations of grain area and rotations as functions of time for the two load cases.

From Fig. 5, we note that in the absence of external loads the two semi-circular GBs, driven by curvature, evolve at identical speeds preserving the symmetry of the grain microstructure. Moreover, the orientation of the grains are preserved throughout the simulation as shown in Fig. 6b. As noted in Section 2.3, the GB motion observed in this simulation may be interpreted as a superposition of coupled GB motion and GB sliding that produces no grain rotation. These features, which are also observed in MD simulations of Trautt and Mishin (2014), can be attributed to coupling factors that are equal in magnitude but opposite in sign for the two GBs, resulting in no grain rotation. In the absence of grain rotation, the equivalence (energetic and kinetic) of the two GBs is preserved, and they continue to evolve with equal speeds. In addition, the evolution history predicted by the continuum simulation is in close agreement with that observed by MD simulation, which proves that the framework can produce accurate results after material parameter calibration.

In the case of shear-driven tricrystal, the GBs experience a driving force due to shear stress in addition to the force due to curvature. Fig. 7 shows that the planar GB migrates upward due to the external driving force. For the lower semi-circular GB, both the curvature and shear stress drive the GB upwards, whereas for the upper semi-circular GB, the two driving forces have a tendency to drive the GB in different directions (curvature flow drives the GB downwards, while shear stress drives it upwards). This leads to a difference in migration velocities, with the lower GB migrating faster than the upper one. This difference in migration velocities, along with GB coupling, leads to a nonzero net rotation rate for the center grain, in the direction such that its crystal orientation is decreased. When comparing the continuum simulation results with those by MD, we see that the two agree qualitatively, i.e., the planar GB migrates upward, the lower GB migrates at a faster rate than the upper one, and the inner grain decreases in orientation.

The qualitative agreements between continuum and MD simulations discussed above provide confidence for the model to be used in more complex scenarios. While using GB energies and mobilities obtained from atomistics would yield a quantitative comparison with MD simulations, we do not pursue such a study due to certain limitations of our framework, which we will now discuss. In our framework, since the construction of GB energy is in terms of the GND density, and the norm of the GND density increases monotonically with GB misorientation, our model is only applicable for small angle GBs. In the tricrystal example, although the misorientations computed using the originally prescribed orientations ($0^{\circ}$, $16.26^{\circ}$ and $73.74^{\circ}$) are large angles, the cubic symmetry enabled us to identify them with small angles by noting the equivalence of grain orientations $\theta$ and $\theta\pm 90^{\circ}$, resulting in the correct GND density that accurately reflects the GB energy. While such an unambiguous choice of GND density is possible for small angle GBs, it is not possible for large angle GBs. For example, consider a tricrystal with grain orientations $0^{\circ}$, $30^{\circ}$ and $60^{\circ}$. Using cubic symmetry, we note that all GBs are equivalent with misorientation angle of $30^{\circ}$. On the other hand, no choice of orientations among their equivalent classes will result in equal misorientations for the three GBs. Therefore, our framework limits us to small angle GBs with a Read–Shockley type GB energy. In Section 4, we discuss how our dislocation-based GB model can complement a disconnection-based framework to model GBs with arbitrary misorientations.

Finally, we remark that the mobility of a GB depends on $\phi$ through the relation given in Eq. (26), which results in the qualitatively correct trend, where GBs with larger misorientation have higher mobilities Trautt and Mishin (2014).¹⁰¹⁰10MD simulations from Trautt and Mishin (2014) highlight the importance of misorientation-dependent GB mobilities in the shape evolution of the inner grain, especially when the two semi-circular GBs have vastly different misorientations. However, we note that the mobility of the GB also depends on orientations of the slip systems relative to the GB. Therefore, in order to make detailed quantitative comparison between our model and MD, further material model parameters (e.g., $b^{\phi}$, $m_{\rm min}$, $m_{\rm max}$ and the functional form of $m$ as a function of $\phi$) calibration is necessary.

In this last simulation, we study the evolution of a square polycrystal under two different loading cases. In the first case, no external load is applied. In the second case, a displacement boundary condition is imposed to intentionally induce grain rotation. The evolution of texture in the two cases are compared, and the importance of external loads in GB evolution is highlighted.

The computational domain is a square polycrystal of side length 100 nm with 13 grains, generated using DREAM.3D (Groeber and Jackson, 2014).¹¹¹¹11DREAM.3D is an architecture for computational microstructure tools that allows users to create ’recipes’ or pipelines for processing digital instances of microstructure. The domain is discretized by quadratic Lagrange elements with a uniform mesh size of 1.6 nm. The crystal orientations are drawn at random from a normal distribution with a mean of 22.5^∘ and a standard deviation of 10.5^∘. Additional checks are performed to limit the GB misorientations to the range 5^∘ to 20^∘, which ensures that all GBs are within the small angle misorientation range to which our model is applicable. The global X- and Y-axes are assumed to coincide with the $[\overline{1}10]$ and $[001]$ directions of a single face-centered cubic (fcc) crystal, such that the in-plane effective slip directions are given by Eq. (32) using $\beta^{\alpha}=0^{\circ}$, $54.7^{\circ}$, and $125.3^{\circ}$ (Kysar et al., 2010). Other simulation parameters are identical to those listed in Table 1, except that no synthetic potential is applied.

A schematic of the polycrystal with a color density plot of initial grain orientations along with applied boundary conditions for the two loading cases is shown in Fig. 8a. The GND density corresponding to the initial grain orientation is shown in Fig. 8b. Since all GBs in the system are low-angle GBs, the norm of the GND tensor serves as a valid measure of GB misorientation. The two loading cases are simulated using boundary conditions imposed on the top edge of the polycrystal. In the first case, the top edge is traction-free, while in the second case, a vertical displacement of 4 nm is imposed on a part of the boundary, shown in red in Fig. 8a, at a rate of $4\times 10^{4}$ nm/s. In both cases, the bottom edge of the domain is held fixed. The duration of the simulation is set to be $2\times 10^{-4}$ s in both cases.

Fig. 9 and Fig. 11 show snapshots (taken at the same time instances) of the time evolution of grains under no external load and uniaxial load cases, respectively. Clearly, in both cases, GB evolution is accompanied by grain rotation. However, we note that uniaxial loading has a noticeable impact on the rate of grain rotation. For instance, a comparison of Fig. 9d and Fig. 11d shows that while grain 1 rotates to decrease its misorientation with the neighbouring grains, the rate of lattice rotation is lower in the presence of external loads. This can be seen in Fig. 9d wherein grain 1 ceases to exist by merging with its neighbor by $t=3\times 10^{-5}$s, while it is still visible in Fig. 11d. The above observation extends to grain 4 as well.

The effect of external load on grain rotation is more evident in Fig. 10, which compares grain lattice rotations for the two load cases at time $t=3\times 10^{-5}$s. The arrows in Fig. 10 depict displacement vectors. From Fig. 10a, we note that in the absence of external load, lattices of grains 1 and 3 rotate clockwise, while the lattice of grain 2 rotates counterclockwise. On the other hand, Fig. 10b shows that the presence of an external load on grain 2, superimposes a counterclockwise rotation in grain 1 and a clockwise rotation in grain 3. This results in a reduced lattice rotation for grain 1 while the lattice rotation is enhanced in grain 3. Moreover, lattice rotation of grain 2 is partially suppressed in the presence of external loads.

In addition, additional rotations are introduced to grains 4 and 5, causing them to slightly rotate in the negative and positive direction, respectively.

In this simulation, equal mobilities are assigned to all available slip systems for simplicity. The work by Thomas et al. (2017) highlights the importance of GB-specific migration behaviors in polycrystal evolution by the MD simulation of a simple toy polycrystal. As shown in Admal et al. (2018), the migration behavior of a GB (e.g., geometric coupling, sliding, or a mix of both) within this framework can be controlled numerically by the available slip systems. For example, it was shown that having two slip systems instead of three, would result in motion by curvature with no grain rotation. Within our framework, various GB mobilities can be explored by having a spatially varying field $c^{\alpha}$ (see Eq. (26)), as guided by experimental measurements or individual MD simulations Thomas et al. (2017), thereby allowing different migration behaviors for GBs in the system.

It is worth noting that the polycrystal simulation reveals a drawback of the diffuse-interface KWC-based models. During the evolution, as the misorientation of the GB decreases, the crystal orientation becomes more diffused. This can be seen in later stages of the two simulations (e.g., Fig. 9d and Fig. 11d). At these time instances, the crystal orientation field become diffused, which renders the identification of the location of GBs difficult. Note that similar diffused GBs can be seen in the polycrystal simulations by Admal et al. (2019). This can be mitigated by increasing the value of $\gamma$ in Table 1, which is used to approximate the singular-diffusive term $|\bm{G}|$ in Eq. (15). However, doing so stiffens the resulting PDEs and worsens the conditioning of the global system, hence it is not pursued here.