Modeling of grain growth under fully anisotropic grain boundary energy

A parameterized model of the five-parameter grain boundary energy landscape for fcc metals was proposed by Bulatov et al in [41]. The model is based on the idea that the energy variation for general grain boundaries can be modeled approximately by scaffold interpolation from lower dimensional grooves to higher dimensions using combinations of twist and tilt boundaries and their energies. The individual combinations are constructed from idealized rotations about the $\langle 100\rangle$, $\langle 110\rangle$ and $\langle 111\rangle$ families of axes, with two additional parameters to define the orientation of the grain boundary plane. For each idealized rotation, a 'distance' value is evaluated to determine to what extent an idealized rotation agrees with the actual grain orientation. The distance value is used in the model to determine how much the energy of an individual tilt/twist combination will contribute to the total computed grain boundary energy.

The grain boundary energy model in [41] was first considered in conjunction with level sets in [31], in order to obtain the energy for a subset of particular low-energy grain boundary configurations in copper. As part of the work in [31], additional simulations based on density functional theory were also performed to confirm the energy predictions provided by the model in [41]. In the present work, the grain boundary energy model is integrated in the proposed level set simulation framework to provide the energies of arbitrary grain boundaries at any stage of the simulation.

As macroscopic deformation of the solution domain is not in focus of the present study, a finite difference scheme, working on a fixed grid, is employed to obtain a numerical solution to equation (10). It can be noted that also finite elements may be used to obtain the solution, as shown in [22, 30, 31]. Working on a fixed grid does, however, bring benefits in terms of avoiding remeshing and having access to very efficient level set reinitialization schemes, as discussed in section 3.1.3.

In the spatial domain, the level set equations in equation (10) are solved using an upwind scheme while an explicit second-order Runge–Kutta scheme is used for the temporal discretization. When evolving the level set functions by equation (10), there will be a tendency to form voids between adjacent level sets. This is an artifact common to all level set implementations, regardless of physical setting, and is usually mitigated by performing an 'interaction correction' step at regular intervals during time integration of equation (10). The most common approach is to assign grid points in the void regions to the closest non-negative level set, providing a purely geometrical reconstruction of the affected level sets. In most level set implementations this is achieved by updating all level sets according to

$$\begin{eqnarray}&&{\phi }_{i}({\boldsymbol{x}})=\displaystyle \frac{1}{2}\left[{\phi }_{i}({\boldsymbol{x}})-\mathop{\max }\limits_{j\ne i}\left({\phi }_{j}({\boldsymbol{x}})\right)\right]\quad \forall i,j\in \{1\ldots {N}_{\phi }\},\forall {\boldsymbol{x}}\in {\rm{\Omega }}.\end{eqnarray} \tag{ 11 }$$

In the present implementation, however, the method proposed in [43] is used, whereby void regions are closed based on an energy minimization principle which is enforced by a Lagrange multiplier related to a constraint that is added to equation (10). Both approaches to interaction correction reduce void regions to be in the order of one grid spacing (or one element size in a finite element setting). Preliminary tests indicate that both approaches work well in most cases but that equation (11) might give rise to a dragging effect on grain boundary junctions in some settings, for example when the conditions of grain boundary energy anisotropy are such that wetting is imminent and sharp grain boundary wedges are formed. Further implementation details of the numerical scheme are provided in the appendix.

3.1.1. Interface normal definition

The interface normal plays a pivotal part in the present modeling of anisotropic interface energies. To avoid issues with discontinuous interface normals, for example at the sharp corners that may develop at grain boundary junctions, the scheme proposed in [44] is adopted. To illustrate this scheme, the standard backward and forward difference operators along one spatial dimension x can be defined as

$$\begin{eqnarray}&&{D}_{x}^{-}{\phi }_{i}=\displaystyle \frac{{\phi }_{i}-{\phi }_{i-1}}{{\rm{\Delta }}x},\ \ \ {D}_{x}^{+}{\phi }_{i}=\displaystyle \frac{{\phi }_{i+1}-{\phi }_{i}}{{\rm{\Delta }}x},\end{eqnarray} \tag{ 12 }$$

where i is the grid position along the x-axis. In two and three dimensions, the difference operators along the other axes are defined in the same way as in equation (12). In two dimensions, denoted by x and y and with grid positions ($i,j$), the local interface normal ${\boldsymbol{n}}$ is evaluated as

$$\begin{eqnarray}\begin{array}{rcl}{\boldsymbol{n}} & = & \left[\displaystyle \frac{{\phi }_{,x}}{{\left({\phi }_{,x}^{2}+{\phi }_{,y}^{2}\right)}^{1/2}},\displaystyle \frac{{\phi }_{,y}}{{\left({\phi }_{,x}^{2}+{\phi }_{,y}^{2}\right)}^{1/2}}\right]\\ & \approx & \displaystyle \frac{1}{4}\left[\displaystyle \frac{\left({D}_{x}^{+}{\phi }_{{ij}},{D}_{y}^{+}{\phi }_{{ij}}\right)}{{\left[{\left({D}_{x}^{+}{\phi }_{{ij}}\right)}^{2}+{\left({D}_{y}^{+}{\phi }_{{ij}}\right)}^{2}\right]}^{1/2}}+\displaystyle \frac{\left({D}_{x}^{+}{\phi }_{{ij}},{D}_{y}^{-}{\phi }_{{ij}}\right)}{{\left[{\left({D}_{x}^{+}{\phi }_{{ij}}\right)}^{2}+{\left({D}_{y}^{-}{\phi }_{{ij}}\right)}^{2}\right]}^{1/2}}\right.\\ & & +\left.\displaystyle \frac{\left({D}_{x}^{-}{\phi }_{{ij}},{D}_{y}^{+}{\phi }_{{ij}}\right)}{{\left[{\left({D}_{x}^{-}{\phi }_{{ij}}\right)}^{2}+{\left({D}_{y}^{+}{\phi }_{{ij}}\right)}^{2}\right]}^{1/2}}+\displaystyle \frac{\left({D}_{x}^{-}{\phi }_{{ij}},{D}_{y}^{-}{\phi }_{{ij}}\right)}{{\left[{\left({D}_{x}^{-}{\phi }_{{ij}}\right)}^{2}+{\left({D}_{y}^{-}{\phi }_{{ij}}\right)}^{2}\right]}^{1/2}}\right].\end{array}\end{eqnarray} \tag{ 13 }$$

Equation (13) provides an average of the four limiting normals at a grid point and if one of the denominators vanish, that component is neglected and the remaining normals are given new weights. The extension of equation (13) to three dimensions is straightforward and for additional details, the reader is referred to [44].

3.1.2. Neighbor identification scheme

In addition to the handling of discontinuous interface normals by equation (13), a second pivotal point in the present level set implementation lies in the determination of neighbors at junctions between multiple level sets, for example at grain boundary triple junctions or along triple lines in 3D. Identification of the correct neighbor is critical when evaluating the grain boundary energy between adjacent crystals. In most level set representations of interfaces, the spatial discretization—by finite element node points or finite difference grid points—is not adapted to the interfaces and the selection of the neighbor crystal at grain boundary junctions is not straightforward. One exception is the finite element-based level set 'interface reconstruction scheme' proposed in [22], whereby two neighboring crystals always share a single unique interface segment and can be unambiguously defined.

In the present framework, a method is proposed whereby the neighbor level set is identified by looking for the minimum angular deviation between the interface normal of a level set with respect to the normals of all other level sets (crystals) that are present at a certain location. This conveniently simple approach turns out to provide a neighbor definition that result in junction configurations that are in agreement with theoretical models, as discussed in section 4. The neighbor identification scheme is illustrated in figure 2 and consists of the following steps:

• 1.  
  For a level set ${\phi }_{p}$ and grid locations $(i,j)$ in 2D, or $(i,j,k)$ in 3D, the set ${ \mathcal S }$ is identified which contains grid points in the vicinity of the zero level, by evaluating:

  $$\begin{eqnarray}&&\,{ \mathcal S }=\left\{(i,j):\left(| \mathrm{sign}({\phi }_{i+1,j})-\mathrm{sign}({\phi }_{i-1,j}| +| \mathrm{sign}({\phi }_{i,j+1})-\mathrm{sign}({\phi }_{i,j-1}| \right)\gt 0\right\}.\end{eqnarray} \tag{ 14 }$$
• 2.  
  For each point in the set ${ \mathcal S }$, all other level sets are checked one at the time using the condition in equation (14). If only a single neighbor is identified in this step, the neighbor identification algorithm is terminated since all neighbors have been identified.
• 3.  
  At junction points or, for example, along triple lines in a 3D setting, more than one neighbor will be identified and a selection has to be performed. To this end, the interface normals ${{\boldsymbol{n}}}_{i}$ of all identified neighbors are evaluated at the present grid point using equation (13), also see figure 2.
• 4.  
  The angle between the normal of ${\phi }_{p}$ and each of the neighbor's normals is evaluated and the neighbor is picked as the level set that provides the minimum angular deviation between the normal directions.

It is emphasized that the proposed neighbor identifications scheme is identical in both 2D and 3D. In the latter case, neighbors at quadruple or other higher-order junctions or along lines where multiple crystals meet are efficiently identified without modifying the algorithm.

Zoom In Zoom Out Reset image size

3.1.3. Level set reinitialization

When tracing the temporal evolution of a level set function by repeatedly solving equation (10), there is usually a tendency for the gradient of ϕ to drift away from satisfying $| {\rm{\nabla }}\phi | =1$ at the interface. This is commonly mitigated by performing a redistancing or reinitialization operation on the level set field with some regularity. In the present implementation, the very efficient reinitialization scheme proposed in [45] is adopted. By this scheme, each level set is first reinitialized using the fast sweeping method (FSM) described in [46]. The FSM is followed by a small number of iterations—only a few iterations are required due to the previous FSM step—using the method presented in [47]. In the latter step

$$\begin{eqnarray}&&\left\{\begin{array}{l}\displaystyle \frac{\partial {\phi }_{i}}{\partial t}-\mathrm{sign}\left({\phi }_{i}^{0}\right)\left(1-| {\rm{\nabla }}{\phi }_{i}| \right)=0\\ {\phi }_{i}({\boldsymbol{x}},t=0)={\phi }_{i}^{0}({\boldsymbol{x}})\end{array}\right.\quad \forall \,i\in \{1...{N}_{\phi }\}\end{eqnarray} \tag{ 15 }$$

is solved iteratively up to a (fictitious) time τ. The obtained solution provides $\phi ({\boldsymbol{x}},\tau )$ as the signed distance function for all points within a distance τ from the interface.

3.1.4. Evaluation of grain boundary energy

In [48], the local equilibrium configuration of a grain boundary triple junction is found to satisfy the Herring equation

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{3}\left({\gamma }_{i}{{\boldsymbol{t}}}_{i}+\displaystyle \frac{\partial {\gamma }_{i}}{\partial \beta }{{\boldsymbol{n}}}_{i}\right)=0,\end{eqnarray} \tag{ 16 }$$

where ${{\boldsymbol{t}}}_{i}$ and ${{\boldsymbol{n}}}_{i}$ are the tangent and normal vector to boundary i, respectively. These quantities are illustrated in figure 3(a).

Zoom In Zoom Out Reset image size

The first component in equation (16) will tend to move the junction to shorten the boundary having the highest energy and the second, so-called torque term, will act to rotate the junction in order to minimize the energy. Any dependence of the grain boundary energy on the boundary plane inclination enters the equilibrium through the torque term in equation (16). However, as noted in [49], the torque term is usually small for general grain boundaries, far away from any 'special' boundary type. In such cases, when neglecting the torque term, equation (16) simplifies into

$$\begin{eqnarray}&&{\gamma }_{i}{{\boldsymbol{t}}}_{i}+{\gamma }_{j}{{\boldsymbol{t}}}_{j}+{\gamma }_{k}{{\boldsymbol{t}}}_{k}=0\end{eqnarray} \tag{ 17 }$$

which, using the law of sines, can be stated as

$$\begin{eqnarray}&&\displaystyle \frac{{\gamma }_{i}}{\sin {\chi }_{i}}=\displaystyle \frac{{\gamma }_{j}}{\sin {\chi }_{j}}=\displaystyle \frac{{\gamma }_{k}}{\sin {\chi }_{k}},\end{eqnarray} \tag{ 18 }$$

where the dihedral angles ${\chi }_{i,j}$ that separate the boundaries i and j, as illustrated in figure 3(b), were introduced. If the grain boundary energy is isotropic, i.e. if ${\gamma }_{1}={\gamma }_{2}={\gamma }_{3}$, the dihedral angles will be ${\chi }_{\mathrm{1,2}}={\chi }_{\mathrm{1,3}}={\chi }_{\mathrm{2,3}}=120^\circ$.

As noted in [15], the dihedral angles between grain boundaries meeting at a junction are always measured in a plane perpendicular to the line junction, see figure 3(a). Hence, a verification of the dihedral angles obtained in 2D simulations is sufficient, as illustrated in figure 3(b). To evaluate the dihedral angles obtained from the present level set model, the idealized triple junction in figure 4(a) is considered. The grain boundary energies are chosen such that ${\gamma }_{2}={\gamma }_{3}$ and an 'anisotropy ratio' R is defined as

$$\begin{eqnarray}&&R=\displaystyle \frac{{\gamma }_{2}}{{\gamma }_{1}}=\displaystyle \frac{{\gamma }_{3}}{{\gamma }_{1}}.\end{eqnarray} \tag{ 19 }$$

The dihedral angle ${\chi }_{2,3}$ can now be obtained from equation (18) as

$$\begin{eqnarray}&&{\chi }_{2,3}=2\arccos \left(\displaystyle \frac{1}{2R}\right).\end{eqnarray} \tag{ 20 }$$

Note that $R\,=\,1$ corresponds to isotropic grain boundary energies and $R\,=\,1/2$ is the special case where wetting will occur and the two boundaries 2 and 3 will tend to merge into one. Figure 4 shows the initial triple junction and the different configurations obtained in the simulations for different choices of the anisotropy ratio R. The results are also shown in figure 5 where the numerically evaluated values of ${\chi }_{\mathrm{2,3}}$ for different values of R are compared to the theoretical values provided by equation (20).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The simulations are performed on a 200 × 200 point finite difference grid, discretizing a domain of $1\times 1$ mm². Neumann boundary conditions are used throughout this paper. The grain boundary mobility is also held constant at $m=1\times {10}^{-5}$ m⁴ J⁻¹ s⁻¹ in the simulations, except in section 4.4 where the mobility is adjusted to comply with the atomistic simulations that are used for comparison.

In the level set simulations, the dihedral angles are evaluated by placing a circle centered around the triple point, as illustrated by the red circle in figure 4. Tracing the perimeter of the circle, the points where each of the three level sets change sign can be interpolated to identify where the circle crosses the individual grain boundaries. Knowing these crossing points, the dihedral angles can be obtained directly. The radius of the circle is set to $5h$, where h is the finite difference grid spacing, similar to the choice made in [39] for evaluation of dihedral angles in a finite element mesh.

As illustrated in figure 4, the initial triple junction configuration is a perfect 'T', but it only takes a few simulation steps before a stable junction configuration is obtained for each value of R. The individual configurations are then maintained throughout the remainder of the simulation and the junctions simply move downwards with approximately constant dihedral angles.

Although it is sufficient to evaluate dihedral angles at triple junctions in 2D, figure 6 provides results from 3D simulations of a junction corresponding to the 2D junction in figure 4, just to illustrate that the proposed level set scheme works equally well in 3D as it does in 2D. Figure 6(a) shows the initial junction configuration between four grains and figure 6(b) shows the configuration when running the simulation with isotropic grain boundary energies, i.e. with $R=1$. In figure 6(c), an anisotropy ratio of $R\,=\,4$ is used, such that the vertical boundary between the two smaller grains (colored in blue and orange), in the bottom half of the domain, has a lower energy than the other grain boundaries. The mobility m is the same as in the 2D case and the $1\times 1\times 1$ mm³ 3D domain is discretized by $150\times 150\times 150$ grid points.

Zoom In Zoom Out Reset image size

In a polycrystal, the grain boundaries will meet at junctions and in most cases form triple lines or quadruple junctions in 3D, or triple points in 2D. Higher-order junctions are energetically unfavorable and will tend to split into multiple junctions of lower energy. Such splitting of unfavorable junctions is an example of a topological transform that frequently takes place during grain growth. To illustrate this transform, the energetically unfavorable 2D quadruple junction in figure 7(a) is considered. If all grain boundary energies are exactly equal, the configuration will be fixated in space. In real polycrystals, however, there will be local heterogeneities present in the microstructure that will initiate and drive the topological transform. Figure 7(b) shows results from a simulation in which the energy ${\gamma }_{\mathrm{1,4}}$ of a grain boundary separating grains 1 and 4, is twice as high as for the other boundaries, which are assumed to have the same energy. Over the course of the simulation, the system will try to avoid creation of such a high-energy boundary between grains 1 and 4 which causes the split shown in figure 7(b). In figure 7(c), the energies are changed such that a boundary between grains 2 and 3 would have twice as high energy as any other boundary in the system and this causes the split shown. It is emphasized that no explicit triggering of a particular topological transform is required in the proposed model. The transforms shown in figure 7 come as a direct consequence of considering grain boundary energy anisotropy in the model.

Zoom In Zoom Out Reset image size

As discussed in section 2, not only the crystallographic orientation across a grain boundary but also the boundary plane orientation will have an impact on the local grain boundary energy variation. These combined dependencies are not considered in the vast majority of mesoscale simulations of polycrystal aggregates, but can be conveniently addressed in the present framework. To provide an illustration, the triple junction configuration in figure 8(a) is considered. Using the same mobility and finite difference grid setup as before, the triple junction evolution is simulated using three different implementations of grain boundary energy. First, an isotropic and constant grain boundary energy of $\gamma =0.625$ J m⁻² is used, which is represented by the solid black lines in figures 8(b)–(d). Under such conditions, the triple junction forms equal dihedral angles and moves downwards without any sideways deflections. Next, GB5DOF is used to identify the energy of each of the initial grain boundary configurations in figure 8(a). This provides γ_1,2 = 0.863, γ_1,3 = 0.685 and γ_2,3 = 0.313 J m⁻². These values are then used in the second simulation case to provide anisotropic, but constant, energies assigned to each grain boundary. Still, no influence of the grain boundary orientation is considered and the outcome is shown by the red lines in figures 8(b)–(d). The introduction of anisotropic grain boundary energies triggers the system to evolve such that the high-energy boundary Γ_1,2 is shortened, compared to the first isotropic case. Finally, also the influence of grain boundary plane inclination is included by evaluating the grain boundary energy by GB5DOF throughout the course of the simulation. The results are represented by the blue lines in figures 8(b)–(d). As illustrated by figures 8(b)–(d), the boundary inclination dependence of the grain boundary energy provides additional driving forces to consume the high-energy boundary Γ_1,2, compared to when constant anisotropic energies are used.

Zoom In Zoom Out Reset image size

The example shown in figure 8 can also be viewed as an illustration of the errors that can be expected when adopting simplified descriptions of the grain boundary energy variation in numerical simulations. As discussed in sections 1 and 2, the most common assumption is to adopt a constant isotropic grain boundary energy which corresponds to the 'isotropic' case in figures 8(b)–(d). In a much more limited number of investigations, anisotropic grain boundary energies are considered by assigning different constant energies to the individual boundaries, corresponding to the 'anisotropic' case in figures 8(b)–(d). As one of very few studies—and probably the only one using level sets—the present work also adds the boundary inclination dependence which provides the 'GB5DOF' case, seen in figures 8(b)–(d).

As discussed in section 4.1, the variation of the grain boundary energy with boundary plane inclination emerges as the torque terms in equation (16). In the vicinity of low-energy boundary orientations, these torque terms will strive to rotate a triple junction into an energetically more favorable position. This causes the additional displacement of the triple junction, shown by blue lines in figures 8(b)–(d), when using GB5DOF to evaluate the grain boundary energy as the system evolves. The impact of the torque terms can be estimated by considering the individual grain boundary energies as functions of the angle of grain boundary plane inclination. This is shown in figure 9, where (a) indicates the normal directions and inclination angle of the individual grain boundaries that meet at the triple junction. Figure 9(b) shows the grain boundary energies as functions of inclination angle, evaluated using GB5DOF. The torque terms in equation (16) correspond to the derivative of γ(θ), which means that the magnitude and direction of the torques are indicated by the slopes of the curves in figure 9(b). Consequently, it can be expected that the torques on the boundaries Γ_1,2 and Γ_1,3 provide a clockwise rotation as the inclination angle increases, while the torque on boundary Γ_2,3 strives in the other direction. An indication of this interaction is visible in figures 8(b)–(d) where the dihedral angle between the boundaries Γ_1,3 and Γ_2,3 is smaller when the boundary plane inclination, i.e. the torque, is considered (blue lines), than when it is not (red or black lines).

Zoom In Zoom Out Reset image size

The interaction and migration of twin boundaries has attracted significant research interest, not least in the field of grain boundary engineering and in relation to nano-crystalline materials. In such studies, the stability of twin grain boundaries and the interaction between coherent (CTB) and incoherent (ICB) twin boundaries is of central interest. The junctions between CTBs and ICBs have been shown to, for example, have a significant effect on the twin boundary concentration in Cu [50]. A range of studies have focused on the dislocation dynamics of such twin grain boundary configurations, revealing that twin boundary migration occurs in Cu under emission of partial dislocations [51]. In [52], it is shown that twin migration enhances plastic deformation in Al. The thermal grain structure stability in nano-crystalline Cu is related to the presence of twin boundaries in [53]. In passing, it can be noted that not only the energy but also the mobility of CTBs and ICBs depend on the local crystal structure, see [54]. The mobility anisotropy is, however, beyond the scope of the present study that focuses on the role of anisotropic grain boundary energy.

To provide yet another example of the applicability of the proposed model, the case of an ICB pinned between two coherent twin boundaries (CTB) is considered. This corresponds to the setup investigated by MD simulations in [55] and is illustrated in figure 10(a). From the present GB5DOF implementation, the grain boundary energies for the ICB and CTB shown in figure 10(a) are obtained as 0.664 J m⁻² and 0.024 J m⁻², respectively. These values are in reasonable agreement with the corresponding values of 0.711 J m⁻² and 0.022 J m⁻², obtained from molecular dynamics simulations in [55]. These values indicate that a considerable energy anisotropy is present in the system, by a factor of nearly 30 between the maximum and minimum energies. Even such a significant anisotropy ratio is handled by the present simulation framework.

Zoom In Zoom Out Reset image size

To comply with the simulation settings used in [55], a 20 × 30 nm² simulation domain is used, as shown in figure 10(a). The domain is discretized using 200 × 300 grid points and the length of the ICB is 15.7 nm, as in [55]. The ICB is initially at a distance of 24 nm from the bottom of the simulation domain. Figure 10(b) shows how the initial grain boundary configuration (black lines) gradually evolves (red lines) under the influence of anisotropic grain boundary energy, evaluated from GB5DOF. The ICB is first arched due to the curvature component of the level set formulation trying to alleviate the influence of the sharp corners between the ICB and the CTBs, and then the ICB moves downwards with its shape maintained. For comparison, figure 10(c) shows the same system evolving under an isotropic grain boundary energy, taken as γ = 0.625 J m⁻². In this case, no one of the three grain boundaries is favored or penalized and the grain boundary triplet collapses as a smooth arc. In the last stage, shown at the bottom of figure 10(c), the equal-energy grain boundaries form a circular arc, before disappearing entirely.

In the atomistic simulations performed in [55], the ICB remains nearly flat and maintains a constant length. Looking at figures 10(b) and (c), it is apparent that arching of the ICB will take place in the present simulations. However, as seen in figure 10(c), the assumption of an isotropic grain boundary energy provides a grain boundary topology that quite drastically differ from a straight ICB of constant length. Although a direct comparison between the present mesoscale modeling and the atomistic simulations done in [55] might be questionable, it is of interest to compare grain boundary migration rates. In [55], the ICB is found to move downwards at approximately 1.2 m s⁻¹ (it is worth noting that the atomistic simulations are performed at a temperature of 1000K). Figure 11 shows the vertical position of the center point of the ICB as function of simulation time. The plots in figure 11 correspond to the two cases of the grain boundary being isotropic and being evaluated continuously from GB5DOF, see figures 10(b) and (c). As seen in figure 11, the center of the ICB stays at the initial position for a short period of time as the initial arching of the ICB occurs. Thereafter, the ICB moves downwards with constant velocity until the very end of the simulations where the Neumann boundary conditions at the bottom of the simulation domain start to influence the shape of the ICB and alter the migration rate. The results in figure 11 are obtained using a grain boundary mobility of m = 1.8 × 10⁻⁸ m⁴ J⁻¹ s⁻¹, which provides a migration velocity of the ICB of about 1.2 m s⁻¹ when using GB5DOF, in agreement with [55]. It can be noted from figure 11 that the assumption of isotropic grain boundary energy significantly increases the migration velocity of the ICB, approximately by a factor of two, compared to when full grain boundary energy anisotropy is considered, using GB5DOF. In comparing our level set-based simulations with the MD simulations in [55], we note that not only grain boundary energy but also grain boundary mobility can be anisotropic and inclination dependent. Although potentially important, anisotropy of grain boundary mobility is beyond the scope of the present study.

Zoom In Zoom Out Reset image size

To provide another 3D example, a spherical grain embedded in a matrix grain is considered. The setup is illustrated in figure 12(a) and the 1 × 1 × 1 mm³ domain is again discretized by 150 × 150 × 150 grid points. The simulations are ran in two scenarios: first, a constant grain boundary energy of γ = 0.625 J m⁻² is used and, second, the grain boundary energy is continuously evaluated using GB5DOF. As an illustration of the grain boundary energy variation, figure 12(b) shows the grain boundary energy for the initial spherical grain. It can be noted that low-energy CTBs are present with boundary normals aligned with the x-axis.

Zoom In Zoom Out Reset image size

The shape and size evolution of the initially spherical grain is shown in figure 13 at three consecutive simulation steps. Figure 13(a) shows the results obtained when using an isotropic grain boundary energy. As seen, the spherical grain maintains its shape throughout the simulation. In figure 13(b), the grain boundary energy is evaluated using GB5DOF and the spherical grain evolves into an ellipsoid shape that is flattened in the y- and z-directions. The presence of low-energy twin boundaries normal to the x-axis means that there is little energy to be lost by changing the curvature along the x-axis. This induces the elliptical arching along the x-axis.

Zoom In Zoom Out Reset image size

Related to grain growth is also the nucleation of new grains, for example as a key stage during recrystallization. The origin of such nucleation events is usually described in terms of either discontinuous subgrain growth or strain-induced grain boundary migration. Regardless of the initiation mechanism at work, a tug-of-war between bulk energy and interface energy driving forces sets the conditions for a critical nucleus in terms of the minimum size required for subsequent growth. The nucleation event is most likely to take place at local microstructure heterogeneities, such as in the presence of orientation gradients, that provide favorable conditions. Nucleation sites are usually found along grain boundaries, within shear bands or near particle inclusions. Solid state nucleation has been the subject for numerous investigations and the related literature is extensive. A recent review on the matter is provided in [56].

One aspect of the growth of a critical nucleus is the energy of the boundaries that the nucleus form with respect to its neighbors. Energy variations can induce directional growth of the nucleus and also cause faceting of the nucleus boundary. As an additional example of application of the proposed modeling framework, the growth and shape evolution of a critical nucleus at a grain boundary is considered. The example is a close analog to the cases studied in the very recent investigation in [57], where a phase field model is used to characterize the faceting of the nucleus surface that appear due to an artificial low-energy boundary plane. The shape of the nucleus obtained from the phase field model is compared to the geometrical shape constructions done in [58, 59], with good agreement. Similar faceting during solid state nucleation has also been shown by atomistic simulations, for example in [60].

In the present example, 3D simulations of the growth of an initially lenticular nucleus, situated at a grain boundary, are performed using the crystal orientations defined in figure 14. As arbitrary grain boundary energy anisotropy and influence of boundary plane inclination is considered, the restrictions in assuming a single artificial low energy plane as is done in [57] are lifted. As can be seen from the crystal orientations in figure 14, a low-energy coherent Σ 3 twin boundary, with its normal in the negative y-direction, exists between the nucleus and the lower grain (grain 2) in figure 14(a). In figure 14(b), the nucleus and the lower grain (grain 2) are tilted such that the twin boundary is inclined 30° from the yz-plane. As in the previous example, see figure 12(b), significant grain boundary energy anisotropy is also present in the polycrystal considered here and is handled by the proposed modeling framework without further modifications of the implementation.

Zoom In Zoom Out Reset image size

To replicate the typical conditions involved in nucleation during recrystallization, a bulk energy difference is introduced in both scenarios in figure 14 such that the nucleus has a zero bulk energy while the upper and lower grains (grains 1 and 2) have equal and positive bulk energies. This bulk energy could for example be due to a higher dislocation density in the cold-worked matrix material, compared to the interior of a recrystallization nucleus and enters the model through equation (8). The precise magnitude of the bulk energy difference is quite arbitrary in the present case and a value of $[[{E}^{b}]]=2750$ J m⁻³ is tentatively set to induce expansion of the nucleus, rather than the shrinkage that would take place under purely curvature-driven grain boundary migration. The simulation domain is 1 × 1 × 1 mm³ in size, although also the domain dimensions are quite arbitrary in the present case, and it is discretized using 100 × 100 × 100 grid points. The radius of curvature of the nucleus surfaces is 0.1 mm.

The evolution of the nucleus morphology in the two polycrystals is illustrated in figure 15 at different stages in time. In figure 15(a), it can be seen how the low-energy boundary—with a normal in the negative y-direction—causes the development of a pointed tip that protrudes straight downwards. In contrast, in figure 15(b) where the low-energy twin boundary is inclined with respect to the vertical plane, the lower tip of the nucleus extends along the normal direction of the twin boundary plane.

Zoom In Zoom Out Reset image size

To illustrate the evolving nucleus morphology more clearly, figure 16 shows cross-sections through the center of the simulation domain, parallel to the coordinate planes.

Zoom In Zoom Out Reset image size

In [57], the co-deformation of the grain boundary plane—on which the nucleation takes place—with the growth of the nucleus is presented as a particular novel finding. A corresponding co-deformation of the grain boundary plane can be observed also in the present results, as seen in figures 16(a), (c), (d) and (f). In addition, the asymmetric shape of the nucleus seen in figure 16(d) is very similar to the nucleus morphology obtained in [57], as well as by the geometric constructions done in [58, 59]. It is, however, emphasized that while an idealized low-energy boundary plane was introduced in [57], no such simplified construction is adopted in the present simulations in which full grain boundary energy anisotropy is considered.

As noted in [57], the co-deformation of the grain boundary plane can be attributed to the influence of the torque terms in the Herring equation, see equation (16), as also manifested in the previous example in section 4.3. The influence of the torque terms permit the grain boundary structure to attain configurations of lower energy than what would be possible without the torques.