A crystal symmetry-invariant Kobayashi--Warren--Carter grain boundary model and its implementation using a thresholding algorithm
Jaekwang Kim, Matt Jacobs, Stanley Osher, Nikhil Chandra Admal
AA crystal symmetry-invariant Kobayashi–Warren–Carter grainboundary model and its implementation using a thresholdingalgorithm
Jaekwang Kim a , Matt Jacobs b , Stanley Osher b , Nikhil Chandra Admal a, ∗ a Department of Mechanical Science and Engineering, University of Illinois Urbana–Champaign b Department of Mathematics, University of California Los Angeles
Abstract
One of the most important aims of grain boundary modeling is to predict the evolutionof a large collection of grains in phenomena such as abnormal grain growth, coupled grainboundary motion, and recrystallization that occur under extreme thermomechanical loads.A unified framework to study the coevolution of grain boundaries with bulk plasticity has re-cently been developed by Admal et al. (2018), which is based on modeling grain boundariesas continuum dislocations governed by an energy based on the Kobayashi–Warren–Carter(KWC) model (Kobayashi et al., 1998, 2000). While the resulting unified model demon-strates coupled grain boundary motion and polygonization (seen in recrystallization), it isrestricted to grain boundary energies of the Read–Shockley type, which applies only to smallmisorientation angles. In addition, the implementation of the unified model using finite el-ements inherits the computational challenges of the KWC model that originate from thesingular diffusive nature of its governing equations. The main goal of this study is to gen-eralize the KWC functional to grain boundary energies beyond the Read—Shockley-typethat respect the bicrystallography of grain boundaries. The computational challenges of theKWC model are addressed by developing a thresholding method that relies on a primal dualalgorithm and the fast marching method, resulting in an O ( N log N ) algorithm, where N is the number of grid points. We validate the model by demonstrating the Herring anglerelation, followed by a study of the grain microstructure evolution in a two-dimensional face-centered cubic copper polycrystal with crystal symmetry-invariant grain boundary energydata obtained from the lattice matching method of Runnels et al. (2016a,b). Keywords:
A. phase field model, microstructures, B. motion by curvature, constitutivebehavior, polycrystalline materials ∗ Corresponding author
Email address: [email protected] (Nikhil Chandra Admal)
Preprint submitted to Journal of Computational Materials Science February 5, 2021 a r X i v : . [ c ond - m a t . m t r l - s c i ] F e b . Introduction Most metals and ceramics exist as polycrystals, which are aggregates of single crystalgrains stacked together along grain boundaries. A polycrystal is characterized by the orien-tation distribution of its grains commonly referred to as grain microstructure. The macro-scopic properties of polycrystals, which include yield strength, resistance to creep, fatigue,thermal and magnetic properties, are strongly influenced by the grain microstructure. Grainboundary engineering refers to the strategy of enhancing the properties of polycrystalline ma-terials by transforming the grain microstructure to a desired state using thermomechanicalprocesses. Mapping the microstructure-property relationship, and modeling the evolutionof microstructure under various manufacturing processes are fundamental open problemsrelevant to grain boundary engineering.The evolution of grain boundaries is driven by a long list of thermodynamic forces, ofwhich surface tension plays a central role. For instance, grain boundaries in the isotropicMullins’ model (Mullins, 1956) are driven by their excess surface energy resulting in motionby curvature with velocity given by v = − mγκ, (1.1)where κ and γ denote curvature and misorientation-dependent energy density of the grainboundary respectively, while m represents a constant mobility. More generally, an anisotropicgrain boundary evolution arises from the dependence of m and γ on the grain boundary char-acter defined by the five macroscopic degrees of freedom, which represent the misorientationand the inclination of the grain boundary. Recent advances in the development of accurateinteratomic potentials have resulted in the mapping of the entire grain boundary energy andmobility landscapes as functions of the grain boundary character (Chen et al., 2020; Run-nels et al., 2016a,b; Olmsted et al., 2009; Bulatov et al., 2013). The grain boundary energylandscape reflects the symmetry of a bicrystal, and understanding the role of crystal symme-try contributes enormously towards characterizing grain microstructure evolution. Moreover,precisely identifying the statistical properties of grain microstructures that are responsiblefor phenomena such as abnormal grain growth and recrystallization remains an open prob-lem in materials science. This motivates us to undertake simulations of large ensembles ofappropriately sampled polycrystals to discover lower-order statistical models for grain mi-crostructure evolution.
The goal of this paper is to develop a lightweight model for motion bycurvature in the presence of a misorientation-dependent grain boundary energy density, thatcan be implemented using an ultrafast algorithm.
While motion by curvature is the simplest description of grain boundary evolution, ex-periments (Rollett, 2018; Barmak et al., 2013) and atomistic simulations (Upmanyu et al.,1998; Janssens et al., 2006) demonstrate that surface tension alone is not the dominant force.Molecular dynamics (MD) simulations have revealed that as grain boundaries evolve, theyplastically deform the underlying material resulting in lattice distortions that give rise toadditional forces on the grain boundary. In order to include the effect of grain boundaryplasticity, recent mesoscale models (Admal et al., 2018) have focused on the coevolutiongrain boundaries and deformation. Evidently, such models subsume motion by curvature as Under an anisotropic energy density γ that depends on the inclination n of a grain boundary, (1.1)transforms as v = − mκ ( γ + ∂ γ/∂n ). ,
000 grainson a 256 × ×
256 grid requires ∼ multi phase field (MPF) model (Chen, 2002;Hirouchi et al., 2012; Steinbach, 2009), and the Kobayashi–Warren–Carter (KWC) model(Kobayashi et al., 1998, 2000; Warren et al., 2003) are two examples of diffuse-interfacemodels for grain boundaries. The main advantage of the MPF model lies in the simplicity ofits construction to include anisotropic grain boundary energies and mobilities. On the otherhand, since the number of phase fields in the MPF method is equal to the number of grains,the MPF method, like the level set and the MBO methods, is memory-intensive. In contrast,the KWC model describes a two-dimensional polycrystal using only two order parameters- one for structural order η ranging from 0 (disordered phase) to 1 (crystalline state), andthe other for crystal orientation field θ . The elegance of the KWC model is offset by thesevere restriction it imposes on the grain boundary energy. The energy functional of theKWC model limits the dependence of grain boundary energy on misorientation angle to aRead–Shockley-type (Read and Shockley, 1950) that does not respect the crystal symmetry.In addition, the singular diffusive nature of the KWC model results in stiff equations thatare computationally expensive to solve.Recognizing the elegance of the KWC model, we formulate a generalization of the KWCmodel that can incorporate arbitrary misorientation-dependent grain boundary energies. Inaddition, we design a thresholding method that addresses the challenge of solving the singulardiffusive equation of the KWC model. The resulting model inherits the memory efficiency ofthe original KWC model, while having significantly more computational efficiency comparedto conventional numerical methods such as finite element and finite difference.The paper is structured as follows. In Section 2, we discuss the role of crystal symmetryon the anisotropy of grain boundary energies, and review the original KWC model and itslimitations. Then, we propose a generalization of the KWC model to incorporate grainboundary energies beyond the Read–Shockley type. In Section 3, we design a thresholdingalgorithm to implement the generalized KWC model. Numerical experiments to validateand demonstrate the computational efficiency of the thresholding method are discussed inSection 4. Finally, we summarize and conclude with a description of future directions.
2. Grain boundary energy and the Kobayashi–Warren–Carter model
A grain boundary is characterized by five macroscopic degrees of freedom where three de-grees represent a rotation associated with the misorientation between the two grains, and theremaining two degrees correspond to the inclination of the grain boundary. More precisely, thegrain boundary character space is given by the topological space T = SO (3) × SO (3) /SO (2),where SO ( n ) is the special orthogonal group in n dimensions. Grain boundaries are equippedwith a surface energy density, which is defined as a function on T . An energy density thatis constant is referred to as an isotropic , and anisotropic otherwise. Fig. 1 shows a plot ofgrain boundary energy density as a function of misorientation angle for a [110] symmetric-tilt grain boundary in face-centered cubic (fcc) copper, calculated using molecular dynamicssimulations (Holm et al., 2010; Bulatov et al., 2014). The role of crystal symmetry in grainboundary energy density is evident from Fig. 1, which shows γ satisfies a two-fold symmetrywith respect to the misorientation angle, a characteristic of all fcc lattices. In addition, γ exhibits local minima at certain misorientations, marked as Σ and Σ in Fig. 1, due to an4 igure 1: A plot of grain boundary energy density as a function of misorientation angle of a [110] symmetric-tilt grain boundary in fcc copper, computed using molecular dynamics Holm et al. (2010); Bulatov et al.(2014). Misorientations corresponding to low energy Σ boundaries are marked on the upper axis. enhanced lattice matching (Runnels et al., 2016a,b; Wolf, 1990) between the two adjoininggrains. Recent efforts (Mason and Patala, 2019; Bulatov et al., 2013; Runnels et al., 2016a,b;Olmsted et al., 2009; Kim et al., 2014) by materials scientists in characterizing the grainboundary character space, and parametrizing grain boundary energy using data from atom-istic simulations and experiments, bring us closer to developing an atomistically informedmesoscale model for grain boundaries.The motion of grain boundaries driven by surface tension to decrease the interfacialenergy is a defining characteristic of various grain microstructure models such as the Mullinsmodel (Mullins, 1956), and its diffuse-interface counterparts such as the multiphase field andthe KWC models. The resulting grain boundary motion, commonly referred to as motionby curvature , in the presence of an anisotropic energy density has been shown to have aconsiderable effect on grain statistics leading to changes in the macroscopic properties ofmaterials. In order to explore the structure-property relationship, recent research effortshave focused on developing ultrafast algorithms to simulate grain boundary evolution in largepolycrystals in the presence of atomistically derived anisotropic grain boundary energies.In this paper, we recognize the simplicity of the phase field model of Kobayashi, Warrenand Carter, and show that it can be generalized and implemented using a thresholdingmethod resulting in an ultrafast algorithm for grain boundary evolution. The Kobayashi–Warren–Carter (KWC) model (Kobayashi et al., 1998, 2000; Warrenet al., 2003) is a dual-phase field model to study grain evolution in polycrystalline mate-rials. In this model, an arbitrary polycrystal in two dimensions is described using only twoorder parameters η and θ . This is one of the main advantages of employing the KWC modelas opposed to the multiphase field model, which uses as many order parameter as the numberof grains. The order parameter η ranges from 0, which signifies disorder, to 1 that describescrystalline order. On other hand, the order parameter θ describes the orientation of the5 a) (b) Figure 2: Results on the one-dimensional steady state solution of the original KWC model describing a flatgrain boundary. a) A steady state analytical solution of the KWC model for a given misorientation. b)Variation of the grain boundary energy as a function of misorientation in the KWC model. grains. The KWC model is governed by a free energy functional given by W [ η, θ ] = (cid:90) Ω (cid:20) f ( η ) (cid:15) + (cid:15) |∇ η | + g ( η ) |∇ θ | + (cid:15) |∇ θ | (cid:21) dV, (2.1)where f = (1 − η ) η = 1, and g = − ln(1 − η ) (2.3)is an increasing function. The functional in (2.1) is defined for all functions η and θ inthe Hilbert space H (Ω). The constant (cid:15) > L -norm) of W , are obtained as (cid:15)b η ˙ η = (cid:15) ∆ η − f (cid:48) ( η ) (cid:15) − g (cid:48) ( η ) |∇ θ | , (2.4a) (cid:15)b θ ˙ θ = ∇ · (cid:20) (cid:15) ∇ θ + g ( η ) ∇ θ |∇ θ | (cid:21) , (2.4b) In three dimensions, the order parameter θ is replaced by a rotation tensor. The choice of the logarithmic function for g is supported by the work of Alicandro et al. (1999), whichshows that the KWC functional converges (in the sense of Γ-convergence) to a surface energy function when g (1) = ∞ (See Theorem 4.1 in Alicandro et al. (1999)). H (Ω) denotes the set of all functions on Ω whose first derivatives are square integrable. b η and b θ are the inverse mobilities corresponding to respective order parameters. A one-dimensional steady state solution of (2.4) under Dirichlet boundary conditions isplotted in Fig. 2a. The value of η < θ is constant in the interior of the grains, andhas a non-zero gradient in a finite thickness around the grain boundary. In the limit (cid:15) → f drives η ( x ) towards 1, while the coupled term g ( η ) |∇ θ | tends to decrease η in a neighborhood of thegrain boundary. In addition, the coupled term tends to localize the jump in θ , while |∇ θ | has a tendency to diffuse it, resulting in a regularized step function for θ . It is interestingto note that in the absence of the |∇ θ | term, the steady state solution for θ is a pure stepfunction, resulting in a model with a blend of sharp- and diffuse-interface characteristics, i.e.while θ is sharp, η is diffused. Moreover, grain boundaries cease to evolve in the absenceof |∇ θ | term (Lobkovsky and Warren, 2001). In other words, |∇ θ | in the KWC modelhas a dual role of not only regularizing θ but also rendering non-zero mobility to the grainboundaries.The grain boundary energy γ gb , as a function of misorientation angle, predicted by theKWC model is of the Read–Shockley-type, as shown in Fig. 2b. This is in contrast to theexperimentally observed grain boundary energies shown in Fig. 1. Despite the elegance ofthe KWC model in describing polycrystals with only two order parameters, its restriction toRead–Shockley type grain boundary energies is a major limitation compared to the flexibilityof incorporating arbitrary grain boundary energies into the multiphase field model. Theabove-stated limitation is one of the main motivation for us to seek a new formulation ofthe KWC model to incorporate arbitrary misorientation-dependent grain boundary energiesthat respect the bicrystallography of grain boundaries. In this section, we formulate a new KWC model that can incorporate misorientation-dependent grain boundary energies. We begin with the KWC functional without the |∇ θ | term. From Section 2.1, recall that in the absence of the |∇ θ | term the steady state solu-tion for θ is a step function with the discontinuity occurring at the grain boundary. Thisobservation motivates us to limit the KWC functional to the set of piecewise constant θ , asopposed to θ ∈ H (Ω), and this enables us to simplify the functional as W [ η, θ ] = (cid:90) Ω (cid:20) (1 − η ) (cid:15) + (cid:15) |∇ η | (cid:21) dV − (cid:90) S ln (1 − ¯ η )[[ θ ]] dS, (2.5) The KWC model was originally developed to simultaneously model grain rotation and grain boundarymotion. The model can be specialized to demonstrate only grain boundary motion by enforcing zero mobilityfor θ in the grain interior. This can be achieved by a constant b φ , and a φ -dependent b θ (Dorr et al., 2010). igure 3: A plot (in green) of the core energy J ( (cid:74) θ (cid:75) ) calculated using (2.8) and (2.9) with grain boundaryenergy data (in red) for [110] symmetric-tilt grain boundaries in fcc copper. where ¯ η : S → R is the restriction of η to the jump set S of θ , which represents theunion of all grain boundaries. The steady state solution, given in (D.8), corresponding to aone-dimensional bicrystal governed by (2.5), and the resulting grain boundary energy as afunction of the misorientation, γ ( (cid:74) θ (cid:75) ) = (cid:74) θ (cid:75) − (cid:74) θ (cid:75) ln (cid:32)(cid:114) (cid:74) θ (cid:75) (cid:33) (2.6)are derived in Appendix D. From (2.6), it is clear that the grain boundary energy is of aRead–Shockley-type, which does not respect the crystal symmetry.The above observation leads us to the following generalization of the KWC functional W G [ η, θ ] = (cid:90) Ω (cid:20) (1 − η ) (cid:15) + (cid:15) |∇ η | (cid:21) dV + (cid:90) S g (¯ η ) J ([[ θ ]]) dS, (2.7)which is defined for all η ∈ H (Ω), and piecewise constant functions θ ; and J is an evenfunction of the jump in orientation. Under this new formulation, the grain boundary energyfunction modifies as γ G ([[ θ ]]) = (1 − η ) − log (1 − η ) J ([[ θ ]]) , (2.8)where ¯ η is the value of the stead-state solution on the grain boundary, given implicitly interms of J ( (cid:74) θ (cid:75) ) as − η ) = J ([[ θ ]]) . (2.9)Inspired from the terminology in dislocations, we refer to J ( (cid:74) θ (cid:75) ) as the core energy .From (2.8), it is clear that by appropriately constructing the core energy J , we can arriveat a γ G that faithfully represents the grain boundary energy and symmetry of the bicrystal. The analog of (2.9) in the original KWC model is (D.7), whose derivation is shown in Appendix D. γ cov of a [110] symmetric tilt grain boundary in face-centered cubic (fcc) copper, shown as red data points in Fig. 3. γ cov is computed using the lattice matching method developed by Runnels et al. (2016a,b), wherein it is defined as thecovariance of the two lattices adjoining the grain boundary. For completeness, we describethe lattice matching method in Appendix A, and list the parameters used to arrive at thedata plotted in Fig. 3. We solve for J in (2.8) using the Newton’s method such that theresulting γ G = γ cov at each data point. In other words, the grain boundary energy of theresulting KWC model is identical to γ cov by construction. Fig. 3 shows a plot of the solution J in green, and highlights common features such as the positions of local minimizers andthe two-fold symmetry, shared between J and γ cov .While the modification of the KWC functional from (2.5) to (2.7) allows us to model arbi-trary grain boundary energies, the absence of |∇ θ | term in (2.7) renders the grain boundariesimmobile as mentioned in Section 2.1. In what follows, we address this shortcoming by de-vising a thresholding method to move grain boundaries by evolving the piecewise-constant θ .
3. Grain boundary motion in the new KWC model
In this section, we present our approach in which we alternate between evolving η and θ to evolve a polycrystal governed by W G . The order parameter η is solved in the followingminimization problem for a given θη ∗ = arg min η ∈ H (Ω) ∂η/∂n | ∂ Ω =0 W G [ η, θ ] . (3.1)Next, the orientation field θ is evolved using a thresholding rule described in the next section.In order to solve for η in (3.1), we note that since g ( η ) → ∞ as η →
1, the functional W G is non-smooth in η , which makes the Newton’s method not viable. Therefore, we usea primal-dual method recently developed by Jacobs et al. (2019), which has a O ( e N log N )complexity, where e is the error in the numerical solution to (3.1), and N is the grid size.See Appendix B for a more detailed description of the primal dual method.Next, we develop a thresholding rule to evolve θ for a fixed η ∗ obtained in (3.1). The alter-nate use of the primal dual method and the thresholding rule at every time step constitutesour approach to evolving the grain boundaries. A thresholding method is a sequence of simple rules, executed every time step to reini-tialize the order parameter, such that its evolution describes the motion of a grain boundary.The original idea of using a thresholding method to evolve grain boundaries goes back to thework (Merriman et al., 1992) of Merriman, Bence and Osher (MBO) wherein, similar to themultiphase field model, grains in a polycrystal are described using as many order parameters,with the caveat that the order parameters are piecewise-constant implying a sharp interface. The evolution of θ by the gradient descent of W G results in pure rotation while the position of the grainboundaries remains fixed. G r a i n B o und a r y Γ l ( x ) /(cid:15) = +˜ ll ( x ) /(cid:15) = − ˜ l B r ( x ) θ R θ L Figure 4: Level sets of the distance function l ( x ) in a neighborhood of x . The grain boundary Γ is depictedas a solid curve, and the dashed curves correspond to the level sets l ( x ) /(cid:15) = ± ˜ l . An order parameter in the MBO method is evolved based on a two-step thresholding scheme— a convolution of the order parameter with a Gaussian kernel followed by a trivial thresh-olding — resulting in motion by curvature. The MBO method has recently been generalizedto a variational model by Esedo¯glu and Otto (2015).In this section, we design a thresholding rule for θ that results in motion by curvature. Wefirst recall that θ is a piecewise-constant field with a finite range of orientations. This implies,that a thresholding rule for θ reassigns θ ( x ), for each point x ∈ Ω, to one of the possibleorientations. Our thresholding rule originates from the observation that the asymmetry of η in the neighborhood of a grain boundary characterizes its curvature. Below, we explicitlyidentify this asymmetry before describing our thresholding rule.First, we note that the steady-state solution for η , derived in (D.8) for a flat interface(zero curvature), is symmetric about the grain boundary. Next, we derive an approximateform for η in the presence of a non-zero curvature, and show the dependence of asymmetryon the curvature. Let Γ denote a grain boundary with a non-zero curvature that separatestwo grains with orientations θ L and θ R . We postulate that in a small neighborhood of x ∈ Γ,the solution η ∗ to (2.4a) is approximated as η ∗ ( x ) ≈ u (cid:18) l ( x ) (cid:15) (cid:19) , (3.2)where l ( x ) is the signed distance function from Γ to x , as shown in Fig. 4. In other words,we assume that η only depends on the radial coordinate. Away from the grain boundary, the10olution to the minimization problem in (3.1) satisfies the equation (cid:15) ∆ η ∗ − ( η ∗ − (cid:15) = 0 . (3.3a)The above equation can be simplified by using a local coordinate system x = (˜ l, s ), where˜ l = l ( x ) /(cid:15) is the scaled radial coordinate, and s is the distance measured along Γ between x and the perpendicular projection of x on Γ. In this coordinate system, we note that (cid:52) η ( x ) = u (cid:48)(cid:48) /(cid:15) + κu (cid:48) /(cid:15) , where κ (˜ l ) is the curvature of the coordinate line { x ∈ B r ( x ) : l ( x ) = ˜ l(cid:15) } .Therefore, (3.3) simplifies as u (cid:48)(cid:48) (˜ l ) + (cid:15)κu (cid:48) (˜ l ) − u (˜ l ) + 1 = 0 . (3.4)Assuming κ (˜ l ) = κ (0), we obtain the following closed form solution to (3.4): u (˜ l ) = 1 + C exp (cid:34) − ˜ l (cid:32) (cid:15)κ + √ (cid:15) κ (cid:33)(cid:35) + C exp (cid:34) − ˜ l (cid:32) (cid:15)κ − √ (cid:15) κ (cid:33)(cid:35) , (3.5)where the constants C and C are determined using the boundary conditions u ( ±∞ ) = 1. Subsequently, the solution can be further approximated under the assumption that both (cid:15) and (cid:15)κ are small, resulting in u (˜ l ) = (cid:40) u (0) − e − (1+0 . (cid:15)κ )˜ l if ˜ l > , u (0) − e (1 − . (cid:15)κ )˜ l otherwise. (3.6)The asymmetry of u is apparent from (3.6) by noting that in the presence of a positivecurvature, the rate at which u → l → + ∞ is greater than when ˜ l → −∞ .The asymmetry of u forms the foundation of our thresholding scheme, which is designedto reassign the values of θ in the neighborhood of the grain boundaries resulting in a motionby curvature. To design a thresholding rule, we identify a unique l = l , such that (cid:90) l −∞ (1 − u ( l/(cid:15) )) dl = (cid:90) + ∞ l (1 − u ( l/(cid:15) )) dl. (3.7)The two integrals in (3.7) are depicted as equal areas under the yellow and green regions inFig. 5, which clearly shows that in the presence of a non-zero curvature, the asymmetry of u results in l (cid:54) = 0. A straightforward but tedious calculation (see Appendix C for details)shows that l = − (cid:15) κ + O ( (cid:15) ) . (3.8)By reinitializing the orientations of all x with l ( x ) < l to θ L , and to θ R when l ( x ) > l , wehave a thresholding rule that moves the grain boundary by (cid:15) κ in one time step dt = t (cid:15) / The boundary conditions are interpreted in the limit (cid:15) →
0, which results in the boundary conditions l/(cid:15) = ±∞ for the scaled radial coordinate. Here, we use the approximation √ (cid:15) κ ≈ O ( (cid:15) κ ). igure 5: A plot of η ∗ in a small neighborhood of x (see Fig. 4) is shown in blue, while (1 − η ∗ ) is shown inred. The asymmetry of u around x due to curvature κ is characterized by the position x at which the twoareas shown in yellow and green regions are equal. The position x is given in terms of l (( x − x ) = l /(cid:15) ),which is the solution of (3.7). where t = 1 is a unit conversion factor. Alternating between the η -update using the primal-dual method, and the θ -update using the thresholding rule, results in a grain boundarymotion by curvature with mobility equal to the inverse of the grain boundary energy. The efficiency of the thresholding rule described above rests on the computation of l in(3.7). In the next section, we use the fast marching method to not only compute l in an O ( N log N ) algorithm, but also generalize the above strategy to an arbitrary polycrystal. The fast marching method (FMM), developed by Tsitsiklis (1995), is an algorithm toevolve a surface with a spatially varying normal velocity. A general description of FMM witha stand-alone example is given in Appendix E. Here, we focus on using the fast marchingmethod to implement the thresholding algorithm described in Section 3.1 by solving for l in(3.7).We begin with a description of our implementation of the thresholding scheme for abicrystal consisting of a circular grain, followed by its generalization to a polycrystal. Wefirst recall that the boundary conditions u ( ±∞ ) = 1 used to arrive at (3.6)–(3.7) apply onlyin the limit (cid:15) → l b >
0, and The case where mobility is equal to inverse of surface tension is referred to as reduced mobilties (Salvadorand Esedo¯glu, 2019; Martine La Boissoni`ere et al., 2019). ! 𝐼 " (a) Two interior regions of a bicrystal. (b) The evolution of ∂I and ∂I solved using the fastmarching method. Figure 6: Shrinking of a circular grain simulated using the thresholding method. a) Two interior regions (redand blue) I and I are grown towards the grain boundary with a speed 1 / (1 − η ∗ ) using the fast marchingmethod. b) A closeup of a rectangular region around the grain boundary, depicted in (a), shows the contourlines of the fast marching method, which describe the time it takes for ∂I or ∂I to arrive at a grid point.Therefore, the original grain boundary, shown as a dashed black line in (b), moves to a new position (solidgreen line) where the two grain interiors meet. modify (3.7) as Find l such that (cid:90) l − l b (1 − u ( l/(cid:15) )) dl = (cid:90) l b l (1 − u ( l/(cid:15) )) dl. (3.9)In Appendix C, we show that the error in l due to the introduction of l b exponentiallydecreases as (cid:15) →
0. In order to use FMM to compute l , we interpret the integrand (1 − u ( l/(cid:15) )) in (3.9) as an inverse of the normal velocity of a surface S l := { x : l ( x ) = l } travelingtowards the grain boundary. Under this interpretation, the integrals in (3.9) are a measureof the time it takes for two initial surfaces S − l b and S l b on either side of the grain boundary,to meet at l = l . We use the fast marching method to evolve the surfaces S − l b and S l b ,and implement the thresholding rule described in Section 3.1 by reassigning the orientationof any point x in the region { x ∈ Ω : | l ( x ) | < l b } to θ L if it first encounters the evolvingsurface S − l b , and to θ R otherwise.In practice, however, we do not have access to the signed distance function l ( x ) to identifythe surfaces S − l b and S l b . Instead, we regularize the surface measure ¯ J ( x ) = J ( (cid:74) θ (cid:75) ( x )) dA on the grain boundaries to obtain a bulk field J (cid:63) ( x ) = G ∗ ¯ J , where G denotes the two-dimensional Gaussian G ( x ) = (1 / π(cid:15) ) e − | x | (cid:15) . Using J (cid:63) , we identify the grain interiors I and I defined as I p = { x ∈ Ω : θ ( x ) = θ p , J (cid:63) ( x ) < ξ } , (3.10)13here ξ > In this work, we choose ξ = 2 (cid:15) . Fig. 6a showsthe grain interiors I and I in a bicrystal, and Fig. 6b is a closeup of a rectangular region,marked in yellow, around the grain boundary. The original grain boundary is marked as ablack dashed line in Fig. 6b. By construction, the two surfaces ∂I and ∂I are equidistant,up to O ( (cid:15) ), from the grain boundary, and serve as substitutes for S − l b and S l b . The graininteriors are grown in the outward direction with a velocity (1 − u ( l/(cid:15) )) − using the fastmarching method, and the surface where they meet is the new grain boundary, shown as agreen dashed line in Fig. 6b.We will now generalize the above implementation to an arbitrary polycrystal consistingof N grains, described using a piecewise constant θ with values in { θ , . . . , θ N } . Using (3.10),we identify the N grain interiors, and define I as their union. Next, we grow the graininteriors in their outward unit normal directions until every point (in the almost everywheresense) in the domain is in precisely one grain. We implement this by first collecting allthe boundaries of the interior regions in ∂I = ∂I ∪ · · · ∪ ∂I n , and simultaneously evolvingthem in the outward normal direction with a speed of 1 / (1 − η ∗ ( x )) using the fast marchingmethod. As the grain interiors grow, a point x ∈ Ω − I is reinitialized to an orientation θ q if it encounters ∂I q ⊂ ∂I . At the end of the fast marching method, all points in Ω − I havebeen reinitialized resulting in an updated polycrystal at the end of a time step. Fig. 7 showsthe implementation of the thresholding rule in a tricrystal. Dirichlet boundary conditions on θ are imposed by including all x ∈ ∂ Ω in the grain interiors. On the other hand, periodicboundary conditions are achieved by periodically reinitializing θ for x ∈ ∂ Ω during the fastmarching step.The primal-dual and the fast marching methods are implemented on a regular grid ofresolution, say δx . From (3.8), we know that the resolution of the grid should be largeenough to resolve a grain boundary movement of (cid:15) κ in each time step, i.e. δx (cid:28) (cid:15) κ, (3.11)which is a common requirement of other thresholding methods (Esedo¯glu and Otto, 2015;Merriman et al., 1992). If the above condition is not satisfied, grain boundaries wouldstagnate. Since grain boundary evolution results in an overall decrease in curvature, (3.11)may cease to hold as the simulation progresses. Therefore, we adaptively increase (cid:15) when agrain boundary stagnates, and as a consequence, we obtain a time adaptive algorithm since dt ∝ (cid:15) . Algorithm 1 summarizes our approach that alternates between the primal-dualand the fast marching methods resulting in motion by curvature. A C++ template librarythat implements Algorithm 1 is available at https://github.com/admal-research-group/GBthresholding .Finally, we remark on the computation of ¯ J ( x ) on a grid. Since (cid:74) θ (cid:75) , calculated at agrid point ij in either x - or y -directions using centered-difference, is shared between two grid Although this regularization results in motion by curvature as in (Merriman et al., 1992), the speed scalesas (cid:15) . This is negligible compared to (3.8). Note that the fast marching method is used to evolve all grain interiors in unison as opposed to evolvingthem individually. ! 𝐼 " 𝐼 Figure 7: Movement of a triple junction according to the thresholding algorithm. The triple junction initiallyat ( x, y ) = (0 . , .
5) (black filled circle) moves to a new position (yellow filled circle) where the three graininteriors, evolved using the fast marching method, meet at the same time. points, a factor of 1 / (cid:74) θ (cid:75) ij = 12 (cid:113) ( θ i +1 ,j − θ i − ,j ) + ( θ i,j − − θ i,j +1 ) . (3.12)
4. Numerical results
In this section, we present examples that explore various features of grain boundaryevolution predicted by our model.We begin with a simulation of a one-dimensional bicrystal Ω = [0 ,
1] with a grain boundaryat x = 0 .
5, and J ( (cid:74) θ (cid:75) ) = (cid:74) θ (cid:75) . The purpose of this simulation is to ensure that the resultsof the primal dual algorithm are consistent with the analytical model described in AppendixD. A Neumann boundary condition dη/dx = 0 is enforced at the two ends. In the absenceof a curvature, we expect the grain boundary to remain at x = 0 .
5, and η reach its steadystate. The tolerance e of the primal dual algorithm (B.8) is set to 10 − . Fig. 8a confirmsthe agreement between η obtained from the primal dual algorithm and the analytical formgiven in (D.8). Furthermore, Fig. 8b shows that the grain boundary energies predicted by theprimal-dual algorithm for various misorientation angles are in agreement with the analyticalresult in (D.10). Next, we study the evolution of a triple junction to its steady state. A triple junction is a line where three grains meet, and it is represented as a point intwo dimensions. The equilibrium of a triple junction is guaranteed if it satisfies the Herring15 lgorithm 1:
Thresholding algorithm for the new KWC model
Input :
A polycrystal with N grains with orientations θ , . . . , θ N , grain boundaryenergies γ ( (cid:74) θ (cid:75) ); parameters: (cid:15) , ξ , total time T , and tolerance e . Output:
Time evolution of the polycrystalConstruct the core energy function J ([[ θ ]]) from grain boundary energy dataInitialize t = 0, and the orientation field θ ( x , while t < T do Compute the discrete jump fields (cid:74) θ (cid:75) ( x , t ) and ¯ J := J ( (cid:74) θ (cid:75) ( x , t )) on ΩRegularize the jump field: J (cid:63) = G ∗ ¯ J , where G ( x ) = (1 / π(cid:15) ) e − | x | (cid:15) // Solve for η ( x , t ) using the primal-dual algorithm Initialize η and the dual field ψ : η ( x ) = 0 , ψ ( x ) = 0, and n = 0 do n = n + 1Calculate η n using ψ n − (B.6)Calculate ψ n using η n (B.7) while (cid:107) η n +1 − η n (cid:107) ∞ ≤ e ; η ( x , t ) = η n +1 // Threshold/update the orientation field Identify interiors of grains: I p = { x ∈ Ω : θ ( x , t ) = θ p , ¯ J ( x ) < ξ } , and set I = ∪ np =1 I p Evolve I with speed 1 / (1 − η ( x , t )) using the fast marching method andupdate/threshold the orientations at each point x ∈ Ω − It = t + 0 . (cid:15) (3.8) end relation (Herring, 1951) given by γ sin Θ = γ sin Θ = γ sin Θ , (4.1)where Θ i is the dihedral angle of grain i , and γ ij is the energy density of the grain boundaryshared by grains i and j . While the Herring relation is derived in the sharp-interface frame-work, it is also seen to hold for a triple junction governed by the original KWC model through(2.4). This is not surprising since the KWC model converges to the Mullins model in thesharp-interface limit and the evolution in (2.4) has a variational structure in the form of agradient descent of the functional in (2.1). On the other hand, it is not clear if our approachto evolve the generalized KWC model arises from a variational formulation. Therefore, itis necessary to examine the Herring relation using our thresholding algorithm. We will nowdemonstrate that the Herring relation indeed holds provided the parameter (cid:15) and the interiorregion of the polycrystal are chosen appropriately.16 imulationAnalytic Solution (a) SimulationAnalytic Solution (b)
Figure 8: A comparison of the numerical solution resulting from the thresholding algorithm, implementedwith (cid:15) = 0 . ×
512 grid, with the analytical solution. Plots of a) the steady-state solution η , and b)grain boundary energy as a function of misorientation. We study the evolution of a triple junction in a tricrystal with orientations θ = 0, θ = π/
6, and θ = π/ , × [0 , J = (cid:74) θ (cid:75) ,we note from Fig. 8b that the energy density of the three grain boundaries are γ = 0 . γ = 0 .
62, and γ = 0 .
87. From the Herring relation in (4.1), it follows that the steadystate dihedral angles are Θ = 134 . ◦ , Θ = 90 . ◦ , and Θ = 134 . ◦ . In order to examinethe Herring relation, we consider a tricrystal under periodic boundary conditions, with aninitial orientation distribution given by θ ( x , t = 0) = θ if x ≤ .
25 or x > . ,θ if 0 . < x ≤ .
75 and 0 . ≤ x < . ,θ if 0 . < x ≤ .
75 and x > .
25 or x > . , (4.2)resulting in four triple junctions at ( x , x ) = (0 . , . , (0 . , . , (0 . , . . , . ◦ , 180 ◦ , and 90 ◦ . Fig. 9a shows aplot of the initial orientation distribution of the tricrystal. The thresholding algorithm isimplemented on a 2048 × (cid:15) = 0 .
02, 0 .
03 and 0 . (cid:15) decreases, the angles in the triple junction converge to thosepredicted by the Herring relation. This test confirms that the Herring relation is satisfiedas (cid:15) →
0. Furthermore, we note that if the construction of the interior I of the grains, asdescribed in (3.10), is altered to I p = { x ∈ Ω : θ ( x ) = θ p , η > − ξ } , (4.3)then all the dihedral angles converge to 120 ◦ . This suggests that the choice of I plays an17 a) Initial Condition (b) (cid:15) =0.04(c) (cid:15) =0.03 (d) (cid:15) =0.02 Figure 9: a) The orientation distribution in an initial tricrystal under periodic boundary conditions with atriple junction at ( x, y ) = (0 . . . , Θ , Θ ) = (90 ◦ , ◦ , ◦ ). The polycrystal isevolved using the thresholding algorithm with (cid:15) = 0 .
03, 0 .
02 and 0 .
01. b)-d) Closeups of an evolving triplejunction (red box) clearly show that the dihedral angles converge to (Θ , Θ , Θ ) = (134 . ◦ , . ◦ , . ◦ )predicted by the Herring angle condition (4.1), as (cid:15) converges to zero. important role in the thresholding algorithm. In this section, we compare the evolutions of a polycrystal resulting form our thresholdingmethod and the finite element method. An initial polycrystal consisting of N = 50 grains, asshown in Fig. 10a, is generated using a Voronoi tessellation of uniformly distributed randompoints. The orientations of the grains are randomly chosen from the interval [0 , π/ J = (cid:74) θ (cid:75) .We implement the thresholding algorithm on a 1024 × (cid:15) = 0 . e = 10 − , and ξ = 2 (cid:15) . A snapshot of an evolving grain microstructure at t = 0 .
18, simulatedusing our thresholding scheme, is shown in Fig. 10b. We confirm that the energy W of thesystem decreases with time, as shown in Fig. 11.The finite element method with Lagrange finite elements cannot be used to implement ourmodel since the solution for θ is discontinuous. Therefore, we implement the finite elementmethod on the regularized KWC model (2.4) using second-order quadrilateral Lagrange finite18 a) (b) (c) Figure 10: An initial polycrystal, shown in (a), is evolved using the thresholding and the finite elementmethods resulting in polycrystals shown in (b) and (c) respectively. The two methods are consistent inpredicting the growth (e.g., ○ , ○ ) and shrinkage (e.g., ○ , ○ ) in various grains. The differences in theevolution is attributed to the mobility function introduced in (4.4) to prevent grain rotation.Figure 11: The total free energy W , predicted by the thresholding algorithm, monotonically decreases as afunction of time. elements. Since the regularized model enables grain rotation, we inhibit rotation using thefollowing η -dependent mobility for θb − θ ( η ) = 10 − (cid:15) + (cid:0) − η (10 − η + 6 η ) (cid:1) (1 − − ) (cid:15), (4.4)as suggested by Dorr et al. (2010). On the other hand, ( b η ) − = (cid:15) is chosen to be constant.To address the singularity due to the |∇ θ | term in (2.4b), we use the approximation g ( η ) |∇ θ | ≈ g ( η ) (cid:112) ρ + |∇ θ | , (4.5)where ρ = 2 × − is a constant. The manifestation of ρ on the solution will be discussedbelow. The finite element simulation is implemented on the open-source computing platform Fenics
Alnaes et al. (2015). We take an implicit time step with dt = 0 . a) (b) (c) Figure 12: An initial polycrystal, shown in (a), is evolved using the finite element method. (b) and (c) showthe resulting polycrystals with regularization parameters ρ = 2 × − and ρ = 2 × − respectively. When ρ is not sufficiently small, grains with small misorientation (e.g., ○ , ○ ) blend out and grain boundarieseasily become rounded. same in the thresholding and the finite element simulations. The grain microstructure at t = 17 .
1, simulated using the finite element method, is shown in Fig. 10c.Comparing Figs. 10b–10c, we note that both the methods are consistent in predictinggrowth (see ○ , ○ ) and shrinkage (see ○ , ○ ) in various grains. It is observed that grainboundaries become rounded in the finite element simulation because of the diffusive natureof orientation field. In addition, disparities are more clear for small misorientation grainboundaries, e.g. ○ , which are highly diffused. This is a manifestation of ρ , which resultsin a non-zero gradient in θ in the grain interiors. In Fig. 12, we compare two finite elementsimulations with ρ = 2 × − and 2 × − , which shows that for a smaller ρ , the grainboundaries retain their characteristic width. However, we note that a decreases in ρ increases the stiffness of the equations, which significantly affects the computational time asdiscussed below.Computational time study clearly highlights the advantage of the our thresholding scheme.The two methods are executed on a single 1.6 GHz core with 8 GB RAM. In Fig. 13, we plotthe dependence of the wall-clock time for the two methods as a function of number of degreesof freedom N . The computational complexity of the thresholding method is O ( N log N ), witha dominant contribution from FFT used in the primal dual algorithm to solve (B.7). On theother hand, the cost of a finite element simulation depends on the choice of ρ , and surgesas N increases. It is clear that for a large collection of grains simulated on a sufficiently finegrid, the thresholding scheme can be orders of magnitude faster than finite elements, whichrequires a small value of ρ to resolve small misorientation grain boundaries. In this section, we examine grain growth in a two-dimensional fcc copper polycrystalwith [110]-type grain boundaries simulated using the generalized KWC model, with crystalsymmetry-invariant grain boundary energy. We compare the results with the predictions of Recall that the characteristic width of a grain boundary in the regularized KWC model is a function of (cid:15) . The [110] direction of each grain is is aligned with the z -axis (out of the plane). lope 1Slope 2 FEM =1 10 -4 FEM =5 10 -5 FEM =1 10 -5 Current scheme
Figure 13: A comparison of the complexity of the thresholding and the finite element methods. While bothmethods have a complexity of O ( N log N ), the cost of the finite element method depends on the choice of ρ . the original KWC model.A two-dimensional polycrystal consisting of N = 50 grains, with orientations in therange [0 , . ◦ ] is generated using a Voronoi tessellation of random points. Fig. 15 showsthe initial orientation distribution in the polycrystal. We assume that the grain boundaryenergy density is independent of inclination. We use the core energy J ( (cid:74) θ (cid:75) ) constructed inSection 2.2 (see Fig. 3). In order to compare the generalized model to the original KWCmodel, we scale the function g of the original KWC model in (2.1) to g = − .
93 ln(1 − η )such that the mean of the grain boundary energies as functions of misorientation in the range[0 , . ◦ ] are identical for the two models. Fig. 14 shows a comparison of the grain boundaryenergy densities of the two models.The orientation distributions of the polycrystal at the end of 200 time steps for the gener-alized and the original KWC models are shown in Figs. 15b–15c respectively. Comparing theresulting polycrystals with the initial polycrystal in Fig. 15a, we note that the generalizedKWC model predicts a growth for red grains while the original model results in their shrink-age. This can be attributed to the relatively smaller energy of grain boundary ○ , which hasa misorientation of ≈ . ◦ in the generalized model due to crystal symmetry.On the other hand, we note an opposite trend for light blue grains for which the gen-eralized model predicts shrinkage while the original model results in a growth. This is aresult of relatively larger energy of grain boundary ○ in the generalized model comparedto the original model. The above observations suggest that the generalized model can re-sult in the growth of certain grains with large misorientation, highlighting the importance ofcrystallography in grain growth.
5. Conclusion and future work
In this work, we generalized the two-dimensional KWC model for grain boundaries suchthat it can incorporate arbitrary misorientation-dependent grain boundary energies that21 igure 14: Grain boundary energies used for the polycrystal simulation in Section 4.3. Using a core energy J ([[ θ ]]) designed in Section 2.2, we obtain a crystal symmetry-invariant KWC model with energy that matchesthe covariance model. In order to compare the original and the new KWC models, we scale the function g ofthe original KWC model in (2.1) to g = − .
93 ln(1 − η ) such that the averages of the grain boundary energies(with respect to misorientation) are identical in the two models. In other words, the areas under the aboveplots are equal. respect the bicrystallography of grain boundaries. In addition, we address the computa-tional challenge of solving a singular diffusive equation of the KWC model by developing an O ( N log N ) thresholding algorithm. Below, we summarize the construction of our model, itsimplementation, and research directions for future work.First, we eliminate the |∇ θ | term in the original KWC model, which is responsible forregularizing the orientation order parameter, and rendering a non-zero mobility to the grainboundaries. The lack of a regularizing term results in a discontinuous orientation order pa-rameter with jump across the grain boundaries, and grain boundaries with no mobility. In thepresence of a piecewise-constant orientation order parameter, the modified KWC functionalcan be separated into a bulk contribution that depends on η , and a surface contribution,called the core energy, that depends linearly on the jump (cid:74) θ (cid:75) across the grain boundaries.Next, we show that by generalizing the core energy from a linear function to an arbitraryfunction J of (cid:74) θ (cid:75) , the model can incorporate arbitrary dependence of grain boundary energieson misorientation angles.Since the absence of the regularizing |∇ θ | term renders the grain boundaries immobile,we design an O ( N log N ) thresholding algorithm to evolve grain boundaries by curvature,where N is the number of grid points. The algorithm, which employs a primal-dual and thefast marching methods, is shown to be an order of magnitude faster than the finite elementimplementation of the original KWC model. We validate our implementation by predictingthe Herring angle relation, and simulate a two-dimensional polycrystal consisting of [110]tilt grain boundaries. The computational efficiency and flexibility of our approach opens thedoor to a number of exiting directions for future work. • The present framework will enable us to carry out a statistical study of large scalesimulations of various ensembles of polycrystals to characterize abnormal grain growthin terms of the grain boundary energy landscape and crystal symmetry.22 a) (b) (c)
Figure 15: (a) A polycrystal with N = 50 grains, and an initial orientation distribution. (b) and (c) showevolved polycrystals using the new and the original KWC models respectively. Grains 1 and 2 show oppositegrowth/shrinkage trends in the two models due to the deviation of the grain boundary energy from theRead–Shockley-type in the new formulation. The blue and red colors represent the maximum and minimumorientation angles of 0 ◦ and 70 . ◦ respectively. • While arbitrary grain boundary energies can be incorporated into our model, its imple-mentation is restricted to grain boundary mobility equal to the inverse of the energy.An extension of our algorithm to include mobilities independently will be explored ina future work. • The present algorithm does not allow grain rotation, which is another important phe-nomenon during recrystallization of polycrystalline materials. We plan to augment thecurrent scheme with a step that models grain rotation. • A recent work by Admal et al. (2019) extended the two-dimensional KWC model to athree-dimensional fully anisotropic model, wherein the dependence of grain boundaryenergy on the misorientation angle was restrictive to a Read–Shockley-type. Due tothe high computational cost of the finite element method, the implementation of thethree-dimensional model was restricted to simple bicrystals. It is envisaged that theefficiency of our thresholding algorithm will enable us to explore large three-dimensionalpolycrystals with fully anisotropic grain boundary energy. • Finally, we recall from the introduction that surface tension is not the only dominantdriving force on a grain boundary due to grain boundary plasticity. Adapting ourthresholding algorithm into existing unified frameworks (Admal et al., 2018), whereingrain microstructure and deformation evolve contemporaneously, will enable us to quan-tify the role of grain boundary plasticity, and study phenomena such as dynamic re-crystallization, superplasticity and severe plastic deformation (Thomas et al., 2017; Weiet al., 2020; Runnels and Agrawal, 2020).
Data availability A C++ template library that implements Algorithm 1 is available at https://github.com/admal-research-group/GBthresholding .23
RediT author statementJaekwang Kim:
Formal analysis, Investigation, Software, Validation, Writing-Originaldraft, Visualization.
Matt Jacobs:
Conceptualization, Methodology, Software, Formalanalysis.
Stanley Osher:
Conceptualization.
Nikhil Admal:
Conceptualization, Method-ology, Writing-Review and Editing, Resources, Data Curation, Supervision, Project Manage-ment.
Appendix A. The covariance model of grain boundary energy
The covariance model for grain boundary energy, developed by Runnels et al. (2016a,b),estimates the grain boundary energy as the extent of lattice mismatch between two misori-ented grains. A lattice mismatch is characterized by the covariance of atomic densities of thetwo lattices adjoining a grain boundary.In the covariance model, a lattice density measure ¯ ρ for a given lattice L , defined as aninfinite sum of Dirac measures with support at the lattice points points of L :¯ ρ ( x ) = (cid:88) d ∈L δ ( x − d ) . (A.2)A lattice density field ρ is introduced as the convolution of ¯ ρ with a thermalization function ξ , i.e. ρ ( x ) = ρ ( x ) ∗ ξ ( x ) , (A.3)where ξ ( x ) = 1 σ π / e −(cid:107) x (cid:107) /σ , (A.4)with σ as the dimensionless temperature. The planar covariance of two thermalized lattices L A and L B with their respective density fields ρ A and ρ B , measured on R , is defined as c [ ρ A , ρ B ] = (cid:90) y ∈ R ρ A ( P T y ) ρ B ( P T y ) λ ( y ) dA, (A.5)where λ ( x ) is an appropriately chosen window function (see (A.8), P : R → R is theprojection P = (cid:18) (cid:19) (A.6)on to the plane R . Expressing the functions ρ A and ρ B in Fourier series, the integral in A lattice L is defined using three lattice vectors l , l , and l as L = { n l + n l + n l n i ∈ Z } . (A.1) igure A.16: A plot of the normalized grain boundary energy versus the misorientation angle predicted bythe covariance model for a [110] symmetric-tilt grain boundary in fcc copper, computed using the relaxationalgorithm of Runnels et al. (2016a,b). For comparison, grain boundary energies obtained from experiment(Miura et al., 1994), and MD simulations (Wolf, 1990) are shown in blue and square points respectively. (A.5) simplifies as c [ ρ A , ρ B ] = 1ˆ λ ( ) (cid:88) k A ∈L (cid:48) A (cid:88) k B ∈L (cid:48) B (cid:98) ρ ( k A ) (cid:98) ρ ∗ ( k B )ˆ λ ( P ( k B − k A )) , (A.7)where k A and k B are lattice vectors of the dual lattices L (cid:48) A and L (cid:48) B , and the window functionis defined in terms of its Fourier transform asˆ λ ( k ) = e −(cid:107) k (cid:107) /ω , (A.8)with an adjustable parameter ω . The grain boundary energy in the covariance model isdefined as γ covgb = E (cid:18) − c [ ρ A , ρ B ] c gs (cid:19) , (A.9)where c gs is the ground state covariance defined as the supremum, over all planes, of c [ ρ A , ρ A ].For example, in fcc, c gs corresponds to covariance measured with respect to the [111] plane.Finally, we note that the covariance model has three adjustable parameters { E , σ, ω } thatcan be used to fit γ covgb to data from experiments or molecular dynamics simulations. It isknown that while (A.9) is a good indicator of grain boundary energy, it over-predicts theenergy for low angle grain boundaries as the above model does not account for facet formation.Runnels et al. (2016a,b) have shown that a further relaxation of the grain boundary energy,which signifies the formation of facets, yields necessary corrections to the energy predictedby the model.Fig. A.16 shows a plot of a relaxed γ covgb computed for [110] symmetric-tilt grain boundariesin fcc copper using E = 1 .
45 J / m − , ω = 0 .
5, and σ/a = 0 . a = 3 .
597 is25he lattice constant of copper. From Fig. A.16, it is clear that the grain boundary energypredicted by the covariance model is in good overall agreement with data from moleculardynamics simulations (Wolf, 1990; Miura et al., 1994).
Appendix B. The primal-dual method
Primal-dual methods is a part of a class of first-order algorithms that have a long historyin the context of optimization problems (Powell, 1978; Kuhn, 1955; Komodakis and Pesquet,2015). As the name suggests, primal-dual methods proceed by concurrently solving a primalproblem and a dual problem. The main benefit of primal-dual splitting is that it replacesan original hard problem with a set of two easy sub-problems (primal- and dual-). Becauseof this advantage, the method has been widely used in diverse fields including compressedsensing, image processing, signal processing, and machine learning (Donoho, 2006; Chambolleand Pock, 2011; Shalev-Shwartz and Singer, 2007; Combettes and Pesquet, 2011).The motivation to use a primal-dual algorithm to solve the minimization problem in (3.1)for η arises due to the presence of a highly nonlinear term g ( η ) J ( (cid:74) θ (cid:75) ) along with |∇ η | .Therefore, we adopt a primal-dual method by introducing an auxiliary dual variable whichenables us to cast (3.1) as an equivalent optimization problem. The choice of the dual variableis based on the observation that (cid:15) (cid:107)∇ η (cid:107) L (Ω) = (cid:15) (cid:107)∇ η (cid:107) L (Ω) − (cid:15) (cid:107)∇ η (cid:107) L (Ω) = − (cid:15) (cid:90) Ω η (cid:52) η dV − (cid:90) Ω (cid:15) ∇ η · ∇ η dV, (B.1)where we have used the divergence theorem, and the Neumann boundary condition ∇ η · n = .Introducing an auxiliary variable ψ , and identifying it with − (cid:15) (cid:52) η , we have (cid:15) (cid:107)∇ η (cid:107) = sup ψ ∈ ( ˙ H (Ω)) ∗ (cid:20)(cid:90) Ω η ( x ) ψ ( x ) dV − (cid:15) (cid:107) ∆ − ∇ ψ (cid:107) L (Ω) (cid:21) = sup ψ ∈ ( ˙ H (Ω)) ∗ (cid:20)(cid:90) Ω η ( x ) ψ ( x ) dV − (cid:15) (cid:107) ψ (cid:107)
2( ˙ H (Ω)) ∗ (cid:21) , (B.2)where ˙ H (Ω) denotes the set of all functions in H (Ω) with zero average, and ( ˙ H (Ω)) ∗ is itsdual. Substituting (B.2) into the KWC functional W G , the minimization problem in (3.1)transforms to the following saddle point problem:inf η ∈ L (Ω) sup ψ ∈ ( ˙ H (Ω)) ∗ Φ[ η, ψ ] , (B.3)where Φ[ η, ψ ] = − (cid:15) (cid:107) ψ (cid:107)
2( ˙ H (Ω)) ∗ + (cid:90) Ω ( ηψ + f ( η )) dV + (cid:90) S g ( η ) J ( (cid:74) θ (cid:75) ) dS. (B.4)The problems of minimizing Φ with respect to η , and maximizing it with respect to ψ are An algorithm that only requires the calculation of the gradient of a functional. η and ψ sub-problems respectively. The advantage of using a primal-dualalgorithm is evident from the observation that Ψ does not depend on the gradients of η ,which renders the η sub-problem local, and the nonlinearity in g ( η ) is no longer a concern.The existence and uniqueness of solutions to the sub-problems follows from standard convexanalysis.We solve for the saddle point of Φ using the following primal-dual update scheme (Algo-rithm 2 in (Chambolle and Pock, 2011)): η n +1 = arg min η ∈ L (Ω) (cid:20) Φ( η, ψ n ) + 12 τ n (cid:107) η − η n (cid:107) L (Ω) (cid:21) , (B.5a) ψ n +1 = arg max ψ ∈ ( ˙ H (Ω)) ∗ (cid:20) Φ(˜ η n +1 , ψ ) − σ n (cid:107) ψ − ψ n (cid:107) L (Ω) (cid:21) , (B.5b)where ˜ η n +1 = (1 + µ n ) η n +1 − µ n η n with µ n = 1 / (cid:112) τ n /(cid:15), τ n +1 = µ n τ n , σ n +1 = σ n /µ n . The scalars τ n and σ n are the step sizes of the η - and ψ -update respectively. The stability(Jacobs et al., 2019; Chambolle and Pock, 2011) of the update scheme in (B.5) is guaranteedif τ n σ n ≤
1. We select τ = (cid:15), σ = 1 /(cid:15) . The solution to (B.5) is obtained by solving thefollowing Euler–Lagrange equations corresponding to gradient flows of the two functionals in(B.5) (cid:18) (cid:15) + 1 τ n (cid:19) η ( x ) + (cid:18) ψ n ( x ) − (cid:15) − (1 + η n ) 1 τ n (cid:19) η ( x ) − J (cid:63) ([[ θ ]]) + 1 (cid:15) − ψ n + 1 τ n η n = 0 , (B.6)(1 /(cid:15) − ∆ /σ n +1 ) ψ n +1 = − ∆(¯ η n +1 + ψ n /σ n +1 ) , (B.7)where the surface measure J ( (cid:74) θ (cid:75) ) dS has been replaced by a volume measure J (cid:63) dV = J ( (cid:74) θ (cid:75) ) exp( − x / (cid:15) ) dV that depends on the distance x from the grain boundary. From(B.6), we note that the primal dual algorithm along with the choice g ( η ) = − log(1 − η ) notonly renders the η sub-problem local but also analytically solvable.We solve (B.6) and (B.7) on a uniform grid of size N = N x × N y . Since (B.6) is solvedanalytically at each grid point, its cost remains O ( N ). We solve for ψ n +1 in (B.7) usingthe fast Fourier transform (FFT), resulting in an O ( N log N ) complexity for the primal dualalgorithm. We use the following stopping criterion for the update scheme in (B.5), (cid:107) η n +1 − η n (cid:107) ∞ = max ≤ j ≤ N | ( η n +1 ) j − ( η n ) j | ≤ e , (B.8)where e is the tolerance of the iterative scheme. Finally, we note that the use of FFT tosolve (B.7) necessitates periodic boundary conditions on η . On the other hand, for Neumann In order to obtain (B.7), we note that the constrained gradient in ( ˙ H (Ω)) ∗ of (cid:82) Ω ˜ η n +1 ψ dV with respectto ψ is −(cid:52) ˜ η n +1 . ψ pq = λ p λ q N x − (cid:88) i =0 N y − (cid:88) j =0 ψ (cid:18) iN x , jN y (cid:19) cos (cid:18) π (2 i + 1) p N x (cid:19) cos (cid:18) π (2 j + 1) q N y (cid:19) , ≤ p ≤ N x − ≤ q ≤ N y − , (B.9)with λ p = (cid:40) / (cid:112) N x , p = 0 , (cid:112) /N x , ≤ p ≤ N x − , and λ q = / (cid:112) N y , q = 0 , (cid:113) /N y , ≤ q ≤ N y − . (B.10) Appendix C. Note on the derivations of thresholding scheme
In this section, we describe the steps to obtain (3.8) from (3.7). We begin by separatingthe domain of integration in (3.7) as (cid:90) l −∞ (1 − u ( l/(cid:15) )) dl = (cid:90) l (1 − u ( l/(cid:15) )) dl + (cid:90) + ∞ (1 − u ( l/(cid:15) )) dl. (C.1)Substituting the solution in (3.6) into (C.1), we have( u (0) − ( (cid:15) − κ ) exp (cid:20)(cid:18) (cid:15) − κ (cid:19) l (cid:21) = ( u (0) − ( (cid:15) − κ ) (cid:18) − exp (cid:20)(cid:18) (cid:15) − κ (cid:19) l (cid:21)(cid:19) + ( u (0) − ( (cid:15) + κ ) . (C.2)Dividing both sides by ( u (0) − and collecting the l terms, we obtain2 exp (cid:2)(cid:0) (cid:15) − κ (cid:1) l (cid:3) ( (cid:15) − κ ) = 1( (cid:15) − κ ) + 1( (cid:15) + κ ) = (cid:15) ( (cid:15) − κ ) . (C.3)Taking a logarithm, we have (cid:18) (cid:15) − κ (cid:19) l = log (cid:18) /(cid:15) /(cid:15) + κ (cid:19) = log (cid:18)
11 + ( (cid:15)κ ) / (cid:19) . (C.4)A Taylor expansion of the right-hand-side of (C.4) with respect to (cid:15)κ/ (cid:18) (cid:15) − κ (cid:19) l = − (cid:15)κ (cid:15) κ − (cid:15) κ O ( (cid:15) κ ) . (C.5)Multiplying by (cid:15) on both sides, we have(2 − κ(cid:15) ) l = − (cid:15) κ (cid:15) κ (cid:15) κ O ( (cid:15) κ ) . (C.6)Finally, using the approximation 2 − κ(cid:15) ≈
2, we get (3.8).As mentioned in Section 3.2, in practice, the infinite bounds of the integral in (C.1) are28eplaced by finite bounds of magnitude l b . Under this change, (C.6) modifies as (cid:18) (cid:15) − κ (cid:19) l = log (cid:32) (cid:15) (cid:15) + κ + 12 exp (cid:18) − (cid:18) (cid:15) − κ (cid:19) l b (cid:19) − (cid:0) (cid:15) − κ (cid:1) (cid:0) (cid:15) + κ (cid:1) exp (cid:18) − (cid:18) (cid:15) + κ (cid:19) l b (cid:19) (cid:33) . (C.7)It can be easily shown that the boxed terms resulting from a finite value of d b decay expo-nentially as (cid:15) →
0, which leaves (3.8) unchanged.
Appendix D. Results on the 1D KWC model
In this section, we collect results on the one-dimensional KWC model which describes aninfinite bicrystal with a grain boundary at the origin. In particular, we present the derivationof the steady-state analytical solution under Dirichlet boundary conditions, and the resultinggrain boundary energy as a function of misorientation.Consider the following KWC energy functional without the |∇ θ | regularizing term W [ η, θ ] = (cid:90) ∞−∞ (cid:20) (cid:15) |∇ η | + (1 − η ) (cid:15) + g ( η ) |∇ θ | (cid:21) dV. (D.1)The Euler–Lagrange equation associated with the above functional is (cid:15) (cid:52) η − η − (cid:15) − g, η |∇ θ | = 0 , (D.2)where g, η is used to denote ∂g/∂η . In what follows, we derive a steady-state solution of (D.2)under Dirichlet boundary conditions η ( ±∞ ) = 1 , θ ( ∞ ) = − θ ( −∞ ) = θ/ . (D.3)We begin with the ansatz that θ ( x ) is a step function satisfying (D.3) with a discontinuityat the origin. Multiplying (D.2) by η (cid:48) , and integrating with respect to x in a region awayfrom the origin, we obtain (cid:15) η, x − (1 − η ) (cid:15) = 0 , = ⇒ η, x = ± (1 − η ) (cid:15) , (D.4)On the other hand, multiplying (D.2) with η (cid:48) , and integrating over an arbitrarily smallneighborhood of 0 results in the jump condition (cid:15) (cid:74) η, x (cid:75) = g, η (¯ η ) (cid:74) θ (cid:75) , (D.5)where ¯ η := η (0) is the value of η at the grain boundary. From (D.4) and (D.5), it followsthat (cid:15)η, x = (cid:40) − η if x > , − (1 − η ) otherwise , (D.6)and g, η (¯ η ) (cid:74) θ (cid:75) = 2(1 − η ) , (D.7)29hich relates ¯ η to (cid:74) θ (cid:75) . The analytical solution for η can be obtained by integrating (D.4).With our choice of g = − ln(1 − η ), the result can be explicitly written as a function ofmisorientation (cid:74) θ (cid:75) : (cid:90) η ¯ η (cid:15) − η dη = x, = ⇒ η ( x ) = 1 − (cid:114) (cid:74) θ (cid:75) (cid:18) − | x | (cid:15) (cid:19) . (D.8)The grain boundary energy γ as a function of misorientation is calculated by evaluating W [ η, θ ] using the steady state solution for η derived above. From (D.4), we have γ ( (cid:74) θ (cid:75) ) = W [ η, θ ] = (cid:90) ∞−∞ (cid:20) (cid:15) η ,x + (1 − η ) (cid:15) (cid:21) dx + g (¯ η ) (cid:74) θ (cid:75) = 2 (cid:90) ∞ (1 − η ) (cid:15) dx + g (¯ η ) (cid:74) θ (cid:75) = 2 (cid:90) η (1 − η ) dη + g (¯ η ) (cid:74) θ (cid:75) = (1 − ¯ η ) + g (¯ η ) (cid:74) θ (cid:75) . (D.9)Note that the grain boundary energy γ and ¯ η are independent of (cid:15) , which reinforces that themodel converges to its sharp interface as (cid:15) → γ with the choice of the logarithmic g , we obtain γ ( (cid:74) θ (cid:75) ) = (1 − ¯ η ) − (cid:74) θ (cid:75) ln(1 − ¯ η )= (cid:74) θ (cid:75) − (cid:74) θ (cid:75) ln (cid:32)(cid:114) (cid:74) θ (cid:75) (cid:33) . (D.10) Appendix E. Fast marching method
The fast marching method (FMM), developed by Tsitsiklis (1995) is used to evolve asurface in the outward unit normal direction with a speed V ( x ) >
0. The fast marchingmethod reformulates a time-dependent initial value problem describing the evolution of asurface into an equivalent boundary value formulation. In this section, we summarize theFMM algorithm as described in Sethian (1996). For illustration, let s ( t ) describe a surfaceevolving with speed V from a given initial surface s (0) = Γ. Instead of solving a time-dependent problem for s ( t ), the fast marching method solves for a function ζ ( x ) whichrepresents the time it takes for the surface to reach x . By the definition of ζ , we have ζ ( s ( t )) = t, (E.1)with ζ = 0 on Γ. Differentiating (E.1) with respect to t , and noting that ∇ ζ is normal to thesurface, we arrive at the following boundary value problem |∇ ζ | V = 1 , ζ = 0 on Γ , (E.2)commonly referred to as the Eikonal equation.Next, we describe the algorithm to solve (E.2) on a two-dimensional grid. In order30o compute |∇ ζ | , an operator D − xij , representing the standard backward finite differenceoperation on the grid point ij , is defined as D − xij ζ = ζ ij − ζ ( i − j ∆ x . (E.3)Similarly, D + x , D − y , and D + y denote forward in x , backward and forward in y finite differenceoperators respectively. To guarantee a unique viscosity solution of the evolving surface,one should consider an upwind finite difference scheme to compute the gradient, which isconveniently written as |∇ ζ | ≈ (cid:2) (max( D − xij ζ, + min( D + x ζ ij , + (max( D − yij ζ, + min( D + y ζ ij , (cid:3) / = (cid:2) (max( D − xij ζ, + max( − D + x ζ ij , + (max( D − yij ζ, + max( − D + y ζ ij , (cid:3) / . (E.4)Using (E.4), we rewrite (E.2) in an algebraic form (cid:2) (max( D − xij ζ, D + x ζ ij , + (max( D − yij ζ, − D + y ζ ij , (cid:3) / = 1 V ( x, y ) . (E.5)Note that if the neighboring values of ζ ij are known, then (E.5) is a quadratic equation for ζ ij that can be solved analytically.The fast marching method begins with the following initialization step1. Assign ζ ( x ) = 0 for grid points in the area enclosed by the initial surface, and tag themas accepted .2. Assign ζ ( x ) = + ∞ for the remaining grid points, and tag them as far .3. Among the accepted points, identify the points that are in the neighborhood of pointstagged as far , and tag them as considered .The key step in the fast marching method is to update ζ with a trial value using (E.5) forgrid points tagged as considered , but only accept the update with the smallest value. Inorder to identify the smallest value efficiently, the grid points tagged as considered are storedin a min-heap structure (Sedgewick and Wayne, 2008) borrowed from discrete networkalgorithms. The fast marching method then proceeds as follows.1. Construct a min-heap structure for the considered points.2. Access the root (minimum value) of the heap.3. Find a trial solution ˜ ζ on the neighbors of the root using (E.5). If the trial solution issmaller than the present values, then update ζ ( x ) = ˜ ζ .4. If a point, previously tagged as far , is updated using a trial value, relabel it as consid-ered , and add it to the heap structure.5. Tag the root of the heap as accepted , and delete it from the heap.6. Repeat steps 2 to 5, until every grid point is tagged as accepted . See Sethian (1996) on the reason behind seeking a viscosity solution. A min-heap structure is a complete binary tree with a property that the value at any given node is lessthan or equal to the values at its children. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure E.17: The level sets of the solution to the Eikonal equation (E.2), computed using the fast marchingmethod, describe a surface evolving with outward normal velocity V ( x, y ) = 1. Fig. E.17 demonstrates the fast marching method used to track an initial surface(9 x − − (3 y + 1)(1 − y ) = 0 , (E.6)growing with a uniform outward normal velocity V ( x ) = 1. ReferencesReferences
N. C. Admal, G. Po, J. Marian, A unified framework for polycrystal plasticity with grainboundary evolution, International Journal of Plasticity 106 (2018) 1–30.R. Kobayashi, J. A. Warren, W. C. Carter, Vector-valued phase field model for crystallizationand grain boundary formation, Physica D: Nonlinear Phenomena 119 (1998) 415–423.R. Kobayashi, J. A. Warren, W. C. Carter, A continuum model of grain boundaries, PhysicaD: Nonlinear Phenomena 140 (2000) 141–150.B. Runnels, I. J. Beyerlein, S. Conti, M. Ortiz, An analytical model of interfacial energybased on a lattice-matching interatomic energy, Journal of Mechanics and Physics of Solids89 (2016a) 174–193.B. Runnels, I. J. Beyerlein, S. Conti, M. Ortiz, A relaxation method for the energy andmorphology of grain boundaries and interfaces, Journal of Mechanics and Physics of Solids94 (2016b) 388–408.W. W. Mullins, Two-dimensional motion of idealized grain boundaries, Journal of AppliedPhysics 27 (1956) 900–904. 32. Chen, J. Han, X. Pan, D. J. Srolovitz, The grain boundary mobility tensor, Proceedings ofthe National Academy of Sciences of the United States of America 117 (2020) 4533–4538.D. L. Olmsted, S. M. Foiles, E. A. Holm, Survey of computed grain boundary properties inface-centered cubic metals: I. grain boundary energy, Acta Materialla 57 (2009) 3694–3703.V. V. Bulatov, B. W. Reed, M. Kumar, Anisotropy of interfacial energy in five dimensions,arXiv: Material Science (2013).A. Rollett, One crystal out of many, Science 362 (2018) 996–998.K. Barmak, E. Eggeling, D. Kinderlehrer, R. Sharp, S. Ta’asan, A. D. Rollett, K. R. Coffey,Grain growth and the puzzle of its stagnation in thin films: The curious tale of a tail andan ear, Progress in Materials Science (2013) 987–1055.M. Upmanyu, R. W. Smith, D. J. Srolovitz, Atomicstic simulation of curvature driven grainboundary migration, Interface Science 6 (1998) 41–58.K. G. F. Janssens, D. Olmsted, E. A. Holm, S. M. Foiles, S. J. Plimton, P. M. Derlet,Computing the mobility of grain boundaries, Nature Materials 5 (2006) 124–127.M. P. Anderson, D. J. Srolovitz, G. S. Grest, P. S. Sahni, Computer simulation of graingrowth I, Acta Metallurgica 32 (1958) 783–791.M. P. Anderson, G. S. Grest, D. J. Srolovitz, Computer simulation of grain growth in threedimensions, Philosophical Magazine B 59 (1989) 293–329.M. I. Mendelev, D. J. Srolovitz, Co-segregation effects on boundary migration, InterfaceScience 10 (2002) 191–199.M. Upmanyu, G. N. Hassold, A. Kazaryan, E. A. Holm, Y. Wang, B. Patton, D. J. Srolovitz,Boundary mobility and energy anisotropy effects on microstructural evolution during graingrowth, Interface Science 10 (2002) 201–216.Z. Yang, S. Sista, J. W. Elmer, T. DebRoy, Three dimensional Monte carlo simulation ofgrain growth during GTA welding of titanium, Acta Materialla 48 (2000) 4813–4825.M. Hillert, On the theory of normal and abnormal grain growth, Acta Metallurgica 13 (1965)227–238.S. M. Allen, J. W. Cahn, A microscopic theory for antiphase boundary motion and itsapplication to antiphase domain coarsening, Acta Metallurgica 27 (1979) 1085–1095.H. J. Frost, C. V. Thompson, D. T. Walton, Simulation of thin film grain structures—I.grain growth stagnation, Acta Metallugica et Materialia 38 (1990) 1455–1462.H. J. Frost, C. V. Thompson, C. L. Howe, J. Whang, A two dimensional computer simulationof capillarity-driven grain growth: Preliminary results, Scripta Metallurgica 22 (1988) 65–70. 33. Kinderlehrer, J. H. Lee, I. Livshits, A. D. Rollett, S. Ta’asan, Mesoscale simulation ofgrain growth, Materials Science Forum 467–470 (2004).D. Kinderlehrer, I. Livshits, S. Ta’asan, A variational approach to modeling and simulationof grain growth, SIAM Journal on Scientific Computing 28 (2006) 1694–1715.H.-K. Zhao, T. Chan, B. Merriman, S. Osher, A variational level set approach to multiphasemotion, Journal of Computational Physics 127 (1996) 179–195.J. Fausty, N. Bozzolo, D. P. Munoz, M. Bernacki, A novel level-set finite element formulationfor grain growth with heterogeneous grain boundary energies, Materials & Design 160(2018) 578–590.B. Merriman, J. K. Bence, S. J. Osher, Diffusion generated motion by mean curvature,Proceedings of the Computational Crystal Growers Workshop (1992) 72–83.S. Esedo¯glu, F. Otto, Threshold dynamics for networks with arbitrary surface tensions,Communications on Pure and Applied Mathematics 68 (2015) 808–864.M. Elsey, S. Esedo¯glu, P. Smereka, Diffusion generated motion for grain growth in two andthree dimensions, Journal of Computational Physics 228 (2009) 8015–8033.M. Elsey, S. Esedo¯glu, P. Smereka, Large scale simulations and parameter study for a simplerecrystallization model, Philosophical Magazine 91 (2011) 1607–1642.A. Zaitzeff, S. Esedo¯glu, K. Garikipati, Second order threshold dynamics schemes for twophase motion by mean curvature, Journal of Computational Physics 410 (2020) 109404.A. M. Jokisaari, P. W. Voorhees, J. E. Guyer, J. Warren, O. G. Heinonen, Benchmarkproblems for numerical implementations of phase field models, Computational MaterialsScience 126 (2017) 139–151.L.-Q. Chen, Phase-field models for microstructure evolution, Annual Review of MaterialsResearch 32 (2002) 113–140.T. Hirouchi, T. Tsuru, Y. Shibutani, Grain growth prediction with inclination dependence of[110] tilt grain boundary using multi-phase-field model with penalty for multiple junctions,Computational Materials Science 53 (2012) 474–482.I. Steinbach, Phase-field models in materials science, Modelling and Simulation in MaterialsScience and Engineering 17 (2009) 073001.J. A. Warren, R. Kobayashi, A. E. Lobkovsky, W. C. Carter, Extending phase field modelsof solidification to polycrystalline materials, Acta Materialla 51 (2003) 6035–6058.W. T. Read, W. Shockley, Dislocation models of crystal grain boundaries, Physical Review78 (1950) 275.E. A. Holm, D. L. Olmsted, S. M. Foiles, Comparing grain boundary energies in face-centeredcubic metals: Al, Au, Cu and Ni, Scripta Materialia 63 (2010) 905–908.34. V. Bulatov, B. W. Reed, M. Kumar, Grain boundary energy function for fcc metals, ActaMaterialla 65 (2014) 161–175.D. Wolf, Structure-energy correlation for grain boundaries in fcc metals—III. Symmetricaltilt boundaries, Acta Metallugica et Materialia 38 (1990) 781–790.J. K. Mason, S. Patala, Basis functions on the grain boundary space: Theory, arXiv preprintarXiv:1909.11838 (2019).H.-K. Kim, S. G. Kim, W. Dong, I. Steinbach, B.-J. Lee, Phase-field modeling for 3D graingrowth based on a grain boundary energy database, Modelling and Simulation in MaterialsScience and Engineering 22 (2014) 034004.R. Alicandro, A. Braides, J. M. Shah, Free-discontinuity problems via functionals involvingthe L1-norm of the gradient and their approximations, Interface and Free Boundaries 1(1999) 17–37.A. E. Lobkovsky, J. A. Warren, Sharp interface limit of a phase-field model of crystal grains,Physical Review E 63 (2001) 051605.M. R. Dorr, J.-L. Fattebert, M. E. Wickett, J. F. Belak, P. E. A. Turchi, A numericalalgorithm for the solution of a phase-field model of polycrystalline materials, Journal ofComputational Physics 229 (2010) 626–641.N. C. Admal, J. Segurado, J. Marian, A three-dimensional misorientation axis- andinclination-dependent Kobayashi–Warren–Carter grain boundary model, Journal of theMechanics and Physics of Solids 128 (2019) 32–53.M. Jacobs, F. Leger, W. Li, S. Osher, Solving large-scale optimization problems with aconvergence rate independent of grid size, SIAM Journal on Numerical Analysis 57 (2019)1100–1123.T. Salvador, S. Esedo¯glu, The role of surface tension and mobility model in simulations ofgrain growth, arXiv:1907.11574 (2019).G. Martine La Boissoni`ere, R. Choksi, K. Barmak, S. Esedo¯glu, Statistics of grain growth:Experiment versus the phase-field-crystal and mullins models, Materialia 6 (2019) 100280.J. N. Tsitsiklis, Efficient algorithms for globally optimal trajectories, IEEE Transactions onAutomatic Control 40 (1995) 1528–1538.C. Herring, Surface Tension as a Motivation for Sintering, McGraw Hill, 1951. doi: .M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring,M. E. Rognes, G. N. Wells, The FEniCS Project Version 1.5, Archive of NumericalSoftware 3 (2015).S. Thomas, K. Chen, J. Han, P. K. Purohit, D. J. Srolovitz, Reconciling grain growth andshear-coupled grain boundary migration, Nature communications 8 (2017) 1–12.35. Wei, L. Zhang, J. Han, D. J. Srolovitz, Y. Xiang, Grain boundary triple junction dynamics:a continuum disconnection model, SIAM Journal on Applied Mathematics 80 (2020) 1101–1122.B. Runnels, V. Agrawal, Phase field disconnections: A continuum method for disconnection-mediated grain boundary motion, Scripta Materialia 186 (2020) 6–10.H. Miura, M. Kato, T. Mori, Temperature dependence of the energy of Cu [110] symmetricaltilt grain boundaries, Journal of Material Science Letters 13 (1994) 46–48.M. J. D. Powell, Algorithms for nonlinear constraints that use lagrangian functions, Mathe-matical Programming 14 (1978) 224–248.H. W. Kuhn, The hungarian method for the assignment problem, Naval Research LogisticsQuarterly 2 (1955) 83–97.N. Komodakis, J. Pesquet, Playing with duality: An overview of recent Primal-dual ap-proaches for solving large-scale optimization problems, IEEE Signal Processing Magazine32 (2015) 31–54.D. L. Donoho, Compressed sensing, IEEE Transactions on Information Theory 52 (2006)1289–1306.A. Chambolle, T. Pock, A first-order primal-dual algorithm for convex problems with appli-cations to imaging, Journal of Mathematical Imaging and Vision 40 (2011) 120–145.S. Shalev-Shwartz, Y. Singer, A primal-dual perspective of online learning algorithms, Ma-chine Learning 69 (2007) 115–142.P. L. Combettes, J.-C. Pesquet, Proximal Splitting Methods in Signal Processing, volume 49of
Springer Optimization and Its Applications , Springer, New York, 2011. doi:10.1007/978-1-4419-9569-8_10