Comparative Study and Limits of Different Level-Set Formulations for the Modeling of Anisotropic Grain Growth

In this study, four different finite element level-set (FE-LS) formulations are compared for the modeling of grain growth in the context of polycrystalline structures and, moreover, two of them are presented for the first time using anisotropic grain boundary (GB) energy and mobility. Mean values and distributions are compared using the four formulations. First, we present the strong and weak formulations for the different models and the crystallographic parameters used at the mesoscopic scale. Second, some Grim Reaper analytical cases are presented and compared with the simulation results, and the evolutions of individual multiple junctions are followed. Additionally, large-scale simulations are presented. Anisotropic GB energy and mobility are respectively defined as functions of the mis-orientation/inclination and disorientation. The evolution of the disorientation distribution function (DDF) is computed, and its evolution is in accordance with prior works. We found that the formulation called “Anisotropic” is the more physical one, but it could be replaced at the mesoscopic scale by an isotropic formulation for simple microstructures presenting an initial Mackenzie-type DDF.


Introduction
The study of GB Thermodynamics and Kinetics are two fundamental topics in materials science. The study of thermodynamics provides information about a system at equilibrium; its extrapolation, under the assumption of local equilibrium, provides the basis for kinetic theories. Additionally, kinetics approaches study the evolution of systems out of equilibrium, involving changes in the microstructure. Determining the kinetics of recovery, grain growth (GG), recrystallization, solidification and other metallurgical mechanisms is necessary to predict and optimize material properties [1]. The need for high-performance materials demands a better knowledge and control of the behavior of GBs under thermomechanical loads. This topic became a strong issue of materials science and gave rise to a branch called GB engineering [2].
In the context of GG, the evolution of GB is driven by the reduction of interfacial energy, and its velocity is classically described, at the mesoscopic scale, by the well-known equation v = µP, where µ is the GB mobility and P = −γκ is the curvature flow driving pressure with γ, the GB energy, and κ, the mean curvature (i.e., the trace of the curvature tensor in 3D). This kinetic equation is a simplification of lower scale phenomena in constant discussions [3,4]. In this case, it constitutes the polycrystalline scale, and in metal forms a state-of-the-art kinematically accepted physical framework. In the discussion of whether this kinetic equation is a reasonable approximation [5] and whether the reduced mobility (µγ product) can really be considered as defined by the temperature and macroscopic properties of the interface as mis-orientation and inclination, a clear and univocal answer seems complicated today. First of all, the answer at the few interfaces scale and at the homogenized polycrystal scale can be contradictory as to the statistical effects. Moreover, a bias in the reduced mobility field discussion today lies in the real capacity of full-field methods to take into account a reduced mobility appropriately defined in the 5D space, defined by the mis-orientation and inclination in representative 2D or 3D simulations. As detailed below, such a capacity is typically unclear in the current state-of-the-art. Thus, the discussion between experimental data and anisotropic full-field simulations is to be treated with extreme caution.
If numerical modeling by considering heterogeneous values of GB mobility and GB energy remains a complex discussion, it has in fact been widely studied at the polycrystalline scale with a large variety of numerical approaches: multi phase-field [6][7][8], Monte Carlo [9,10], molecular dynamics [11], orientated tessellation updating method [12], vertex [13], front-tracking Lagrangian or Eulerian formulations in a finite element (FE) context [14][15][16], and level-set (LS) [17][18][19], to cite some examples. During annealing, two properties have been widely studied: the GB energy and mobility. The first models proposed in the literature define the GB mobility and energy as constants, carrying the name of isotropic, [6,9,17,20,21], and this category shows good agreement in terms of mean quantities and distributions; nevertheless, they are restrictive in terms of the grains' morphology and texture predictions. GB energy and mobility were earlier reported as anisotropic by Smith [22] and Kohara [23]. Hence, the models have evolved in order to reproduce more complex microstructures or local heterogeneities, such as twin boundaries. Heterogeneous models were proposed, in which each boundary has its own energy and mobility [10,18,19,[24][25][26][27][28][29][30][31]. For instance, every grain could be related with an orientation, thus the mobility and energy can be computed in terms of the disorientation [7,19], but the mis-orientation axis and inclination dependence are frequently not taken into account. Finally, general frameworks in which the five parameters, mis-orientation, and inclination are discussed have been proposed, and these models could be categorized as fully anisotropic [32][33][34].
However, it must be highlighted that the distinction between 3-parameter and 5parameter full-field frameworks is not straightforward because of unclear terminology. In the literature, heterogeneous values of GB properties have often been categorized as anisotropic. For instance, in [7,12,29,31,35], heterogeneous GB energy and a constant GB mobility to model polycrystal evolution during GG are considered, and the models are categorized as anisotropic even if it is assumed that the GB energy does not depend on the GB normal direction and the GB mobility is not heterogeneous. In [34], the proposed LS formulation in the context of regular grids includes the effect of anisotropic GB energy into the driving force term (P) using both the effect of the mis-orientation and the inclination in a GB energy gradient. However, the GB energy dependence on the normal direction is defined without inquiring if additional torque terms in the solved equations are negligible or not.
Due to the wide variety of formulations, this paper aims to compare four different formulations within an FE-LS approach. The first is an isotropic formulation frequently used in different contexts, such as GG, recrystallization, and GG with second-phase particles [17,[36][37][38][39], referred to as the isotropic model in the following. The second one is a simple extension of the isotropic formulation by considering non-homogeneous values of the reduced mobility. The third formulation was firstly proposed in [26] and extended to polycrystals using different models of GB energy in [19]. The last formulation is based on a more robust thermodynamics and differential geometry framework, but was only applied, as yet, to a bicrystal-like geometry [33]. Another particularity of the discussed approaches is to be usable on unstructured finite element mesh and in the context of large deformations and displacements. The goals of this work are to criticize these existing formulations but also to consider the enrichment of GB mobility in the FE-LS framework. First, some crystallographic definitions, LS treatments, and formulations are introduced in Section 2. In Section 3, simulation results are compared with analytical solutions in the context of simple triple junction geometries. In Section 4, polycrystalline simulations are studied. Mean values and statistical quantities are compared with two different initial textures and using heterogeneous GB energy and mobility. Finally, the last section is dedicated to the inclination dependence discussions.

The Numerical Framework
Before formulating the equations related to GG, the constituents of polycrystalline materials and especially GB structures must be defined.

Crystallographic Definitions
Let us consider a domain Ω of dimension d filled by n grains G i ∈ Ω, being open spaces of Ω and defining the set of grains G = {G i , i = 1, . . . , N G }. The interface between two neighboring grains G i and G j constitutes a GB B ij , and the whole set of boundaries form the GB network Γ. A boundary B ij is characterized by its morphology and its crystallographic properties, which are described using five variables: 2 shape properties, describing the interfaces by the unitary-outward normal direction n ij , and 3 crystallographic properties describing the orientation relationship between the two adjacent grains, O i and O j , known as the mis-orientation M ij . As such, at the mesoscopic scale, each boundary may be characterized by a tuple: The GB space, B, parameterized by the mis-orientation and the normal direction is illustrated in Figure 1. The two quantities of interest, the GB energy γ and the mobility µ, are mapped from B to R + .

FE-LS Formulation
The LS method is a powerful tool firstly proposed by Osher and Sethian [40] to describe curvature flow of interfaces, enhanced later for evolving multiple junctions [41,42], and considered in recrystallization and grain growth problems in [17,36]. The principle for modeling polycrystals is as follows: the grain interfaces are defined through scalar fields called LS functions φ in the space Ω, and more precisely by the zero-isovalue of the φ functions. LS functions to the interfaces are classically initialized as the sign Euclidean distance functions to these interfaces: with d representing the Euclidean distance and φ generally being defined as positive inside the grain and negative outside. The dynamics of the interface is studied by following the evolution of the LS field. The interface may be subjected to an arbitrary velocity field, v, and its movement is described by solving the transport equation: The flexibility of this method lies in the ability to define different physical phenomena encapsulated in the velocity field. This equation is solved to describe the movement of every grain. When the number of grains increases, one may use a graph coloring/recoloring strategy [38] in order to drastically limit the number of involved LS functions: Φ = {φ i , i = 1, . . . , N}, with N N G being N G , the number of grains. Additionally, two more treatments are necessary. Firstly, the LS functions are reinitialized at each time step to keep the metric property of a distance function: Secondly, the evolution may not preserve the impenetrability constraints of the LS functions, leading to overlaps and voids between grain interfaces. These events are corrected after solving the transport equation by resolving Equation (4), as proposed in [41] and classically used in the LS framework [17,43]: Several formulations using the LS framework exist in the literature. The initial GG formulation uses a homogeneous grain boundary energy and mobility, i.e., [17], and the velocity field is thus defined as: where P = −γκ is the capillarity pressure and n is the outward unitary normal to the interface. When dealing with recrystallization, supplemental terms could be added to the velocity, as proposed in [17]. If φ is defined as positive inside the grain and remains a distance function, the mean curvature and the normal may be defined as: then, the velocity in Equation (5) may also be defined as: Four different formulations will be studied. In the first one, an isotropic formulation is considered by introducing Equation (7) into Equation (2), thus the isotropic transport equation may be defined as a pure diffusive problem: This formulation has shown good agreement with experimental data regarding GG predictions concerning the mean grain size and even the grain size distribution (GSD). However, this approach is limited when it comes to reproducing complex grain morphology (non-equiaxed ones), described as special grain boundaries, and with respect to textures. This formulation could be slightly modified in a second one with the introduction of heterogeneous GB properties, leading to a heterogeneous formulation: With this formulation, it is expected to obtain more physical grain shapes. Indeed, some GBs can evolve faster thanks to higher grain boundary mobility values, and triple junctions may have different dihedral angles thanks to different GB energy values. This strategy, classically used in full-field formulations (not only in LS ones), can lead to confusion when it is named as "heterogeneous". Indeed, stricto sensu, the heterogeneity shape of µ and γ can lead to additional terms in the driving pressure of the kinetic equation, Equation (5), but also in the weak formulation derived to solve the GB motion. However, the term "heterogeneous" will be used in the following to distinguish this formulation from the purely isotropic model.
Such discussion is described in [26], where an additional term capturing the local heterogeneity of the multiple junctions is added to the velocity equation, such that: v = µ( ∇γ · ∇φ − γ∆φ) ∇φ. (10) Inserting this term into the transport equation, Equation (2), leads to the, hereafter called, "Heterogeneous with Gradient" formulation [26]: The introduction of the term ∇γ · n only acts at multiple junctions because these are the only places where this term does not vanish. This formulation is equivalent to the isotropic one if no heterogeneity is added.
Finally, in [33,44], a new relation for the velocity was developed using thermodynamics and differential geometry. The five crystalline parameters are taken into account with an intrinsic torque term, which leads to (see Equation (2.43) in [44]): where m αβ is the metric with components α and β of a Riemannian n-manifold, with n the dimension of the space, and∇ the Levi-Civita connection. This equation may be redefined using a flat metric and tensor notations as: and also written as: where I is the unitary matrix, P = Id − n ⊗ n is the tangential projection tensor, and therefore, ∇ n γ = P ∇γ with ∇ n , the surface gradient on the unit sphere of interface normal n, K = ∇ n = ∇ ∇φ, is the curvature tensor. In Equation (12), the term P αβ∇ β γ∇ α φ and its equivalent P ∇γ · ∇φ in Equation (13), i.e., ∇ n γ · ∇φ, should be null in the grain interfaces. However, the front-capturing nature of the LS approach, which requires to solve Equation (13) not only at the GB network but also in its vicinity, needs to consider this term, which could be non-null around the GB interfaces. This stabilization term is then totally correlated to the front-capturing nature of the LS approach and not derived from the GG driving pressure. The resulting tensorial diffusion term, D = ∇ n ∇ n γ + γI [33,44], is also well-known as the GB stiffness tensor Γ( n) in [45,46]. With this formulation, the 5D-GB space B is fully described and is referred to as "Anisotropic-5". If the torque term is neglected, the formulation used could be simplified as: This equation is hereafter called "Anisotropic" and is not equivalent to the "Heterogeneous with Gradient" formulation Equation (11). The strong formulations used in this work are finally the ones defined by the Equations (8), (9), (11) and (15). Moreover, the effect of heterogeneous GB mobility is take into account in the weak formulations in the form of a GB mobility gradient in the Heterogeneous with Gradient and Anisotropic formulations. The weak formulations of Equations (8), (9), (11) and (15), with ϕ ∈ H 1 0 (Ω), can be summarized as: and All the presented formulations are equivalent if the properties are homogeneous, but the main question remains as to the test of their capacity otherwise. In the next sections, a comparative study is presented. In the following, the "Isotropic", "Heterogeneous", "Heterogeneous with Gradient", and "Anisotropic" formulations will be referred to as Iso, Het, HetGrad, and Aniso. It must be highlighted that the formulations proposed in Equations (18) and (19) are slightly more general than those proposed in [26,33] respectively, as here, µ is also considered as heterogeneous.

Description of the Test Case
In this section, simulation results obtained with the Het, HetGrad, and Aniso formulations are compared for a 2D-triple junction configuration proposed in [6] and described in Figure 2. The initial microstructure is a dimensionless T-shape triple junction with L x = 1 and L y = 3. This geometry was chosen because after a transient-state, a quasi-steady-state is reached, where analytical relations, depending on the reduced mobility, are available for the triple junction velocity and equilibrium angles. When the quasi-steady-state is reached, the triple junction moves with a constant velocity towards the bottom of the domain, with a stable triple junction profile which respects the conditions imposed by the Herring's equation [47]: where γ ij is the GB energy and τ ij are the inward pointing tangent vectors of the three boundaries at the triple junction. In the present example, the grain boundary energy is constant per interface (γ(M)) and the above equation may as well be expressed by the Young's law (no torque terms): which may be expressed in terms of the angles ξ i of the grain i, through the Young's equilibrium (see Figure 2): By considering an axially symmetric configuration where γ 13 = γ 23 = γ top and γ 12 = γ bot , and by defining the ratio of grain boundary energies as r = γ top γ bot , an analytical value for the angle ξ 3 can be obtained: Moreover, the stationary transported profile takes the form of the "Grim Reaper" profile, defined as: where v ana T J is the magnitude of the stationary velocity, y 0 is the initial y-value, and (x, y) are the Cartesian coordinates. By using Neumann boundary conditions, the stationary velocity could be related to the x-size of the domain: In order to focus on a considerable level of heterogeneity in the system, r is initially fixed as equal to 10 (γ top = 1 and γ bot = 0.1), and µ is defined as unitary. Several simulations were carried out and compared with the analytical values of ξ ana The dihedral angles are computed using the methodology presented in [26]: one may define, at each time, a circle of radius ε with circumference C ε , and divide it into arcs which pass through grain G i with length L i ε . The angle of the arc, ξ i , could be approximated thanks to the relation ξ i = 2πL i ε /C ε . Hence, these variables are affected by the spatial discretization of the domain and the choice of ε, which must be close enough to the multiple junction while containing a sufficient number of finite elements, as illustrated in Figure 3, where different values of are tested. Here, the value ε = 0.05 is adopted. v T J and ξ i are compared using relative errors which are defined as: where X ana is the analytical value of the variable to be compared. Another discussed quantity is the interfacial energy, calculated using: where T is the set of all elements in the FE mesh, l e is the length of the zero iso-value existing in the element e, and i refers to the number of LS functions, and the 1 2 is necessary due to the duplicity of the LS functions in the interfaces defining a grain boundary. This variable is frequently studied and it may be seen as an energetic measure of how quickly the system reaches equilibrium.

Numerical Strategy
The simulations presented here were carried out with unstructured triangular meshes, a P1 interpolation, and using an implicit backward Euler time scheme for the time discretization. The system is assembled using typical P1 FE elements with a Streamline Upwind Petrov-Galerkin (SUPG) stabilization for the convective term [48]. The boundary conditions (BCs) are classical null-von Neumann BCs applied to all of the LS functions. This choice imposes the orthogonality between the LS functions and the boundary domain (each plane of the boundary domain can be seen as a symmetric plane). By considering a minimal and maximal mesh size (respectively h min and h max ), an optimized anisotropic re-meshing strategy developed by Bernacki et al. [37,49], used in the DIGIMU ® software [50] and illustrated in Figure 4, is adopted here. The mesh is finely and anisotropically refined close to the interfaces (φ < φ min ) and becomes isotropic when φ > φ max , with a linear evolution of the normal mesh size between φ min and φ max . A homogeneous tangential mesh size (h t = h max ) is considered everywhere and the normal mesh size is then defined as: By generalizing this approach at the multiple junctions, a fine isotropic (h n = h t = h min ) re-meshing is automatically performed (see [37] for more details). During grain boundary migration, thanks to a topological mesher/re-mesher, anisotropic re-meshing operations are performed periodically to follow the grain interfaces. Typically, a re-meshing operation is considered each time a LS is about to leave the fine mesh area set by φ min .

Results and Analysis
First, a sensibility analysis for the three formulations was carried out. The values of mesh size and time step used here are: h max = h t = 1 × 10 −2 , h min = {5 × 10 −4 , 1 × 10 −3 , 5 × 10 −3 , 1 × 10 −2 } and ∆t = {1 × 10 −5 , 5 × 10 −5 , 1 × 10 −4 , 5 × 10 −4 }. For all the cases, Φ min and Φ max are fixed respectively to 1 × 10 −2 and 2 × 10 −2 . Figure 3 shows, for the different formulations, the triple junctions at t = 0.25 using h min = 5 × 10 −4 and ∆t = 1 × 10 −5 . One dihedral angle is depicted for different values of ε. In the following, ξ is used to define the converged value of the ξ 3 angle for the k-th value of the h min and ∆t datasets. Indeed, if the results described in Figures 5-7 principally aim to compare the simulations with the quasi-steady-state analytical values, it is also interesting to discuss the obtained converged value of v T J as a function of the converged value of ξ 3 (i.e., if Equation (25) is respected for these values). Figure 5 illustrates the evolution of E Γ , ξ 3 , and v T J using the Het formulation. Two stages appear in E Γ , whereby it initially increases before decreasing. The results illustrate the fact that the approach seems not to converge, in time and space, towards the analytical solutions. However, in terms of the dihedral angle, the results converge towards  (25)). The movement of the Het formulation is mostly influenced by the curvature of the interface, as exposed in Section 2, and one has to keep in mind that there are no additional terms that could influence the movement of the interfaces. These results illustrate that the Het formulation, by considering heterogenous values of reduced mobility and the multiple junction treatment defined by Equation (3), without re-discussing the capillarity driving pressure used in the kinetic equations, is definitively not a good option when a convective/diffusive formulation is solved to model the GG mechanism.       The evolution of the HetGrad formulation is quite different, the interface evolves in the opposite direction (see Figure 8) which explains that E Γ increases during the simulation (see Figure 6). An explanation of this evolution comes from the presence of the grain boundary energy gradient, ∇γ, in the triple junction. The main purpose of this gradient is the correction of the triple junction dihedral angles and velocity. In Figure 6, one can see that The Aniso formulation has an additional term, the projection tensor P, which takes into account the tangential changes of ∇γ. Thanks to this term, the interface evolves in the right direction, with a minimization of the boundary energy. , meaning that the kinetics and topology of the triple junction are well-correlated through Equation (25).
The evolution of the triple junction profile using the Het, HetGrad, and Aniso formulations is illustrated in Figure 8. Both the Het and Aniso formulations produced the Grim Reaper profile, while the profile produced by the HetGrad formulation evolves in the opposite direction. This is reflected in the values of the triple junction velocity shown in Figure 9c. From the comparison of the interfacial energy evolution (Figure 9a) and of the velocities (Figure 9c), one can see that the Aniso formulation has the best energetic behavior and a better approximation of the triple junction velocity. However, the best approximation of dihedral angles is obtained with the HetGrad formulation (Figure 9b). The level of anisotropy defined here is high (r = 10), and this order of value has also been discussed in the literature [26,30,51] and remains necessary to discuss realistic polycrystal aggregates (coherent twin energy, for example). In Figure 10, the effect of the anisotropy level (r value) on the top dihedral angle and the triple junction velocity is illustrated. We have carried out simulations using h = 1 × 10 −3 , ∆t = 1 × 10 −4 and r ∈ {0.55, 0.625, 0.714, 0.833, 1.0, 1.25, 1.66, 2.5, 5, 10}, which are equivalent to γ bot ∈ {1.8, 1.6, 1.4, 1.2, 1.0, 0.8, 0.6, 0.4, 0.2, 0.1}. These results allow us to conclude that the Het methodology is not adapted whatever the r value. Interestingly, the HetGrad formulation seems very good for ξ 3 and v T J for r < 1.5, but the migration direction ends up being reversed for higher r values, while keeping an excellent profile for the equilibrium angles. Finally, if the angle respect is slightly worse for the Aniso formulation, the respect of the triple junction speed is much better as soon as r > 1. In Figure 10b, the three additional dashed lines represent the expected velocity for the ξ 3 values obtained after reaching equilibrium, as illustrated in Figure 10a using Equation (25). One can see that the Het and HetGrad formulations correlate ξ 3 and v T J for r < 1. On the other hand, the Aniso formulation correlates ξ 3 and v T J for every r value. A good correlation could be advantageous if one wants to perform more realistic simulations where correct kinetics and topology of the microstructure are of significant importance. An additional analysis is presented in the Appendix A, where a modified Grim Reaper case with Dirichlet boundary conditions is detailed. The obtained results have similar trends to those presented in Figure 10a.

Conclusions
These results highlight that the Aniso formulation seems to be the most physically acceptable approach regarding the velocity of the triple junction and the interfacial energy evolution. Additionally, it also correctly represents the dihedral angles for a wide range of anisotropy levels. Nevertheless, these results must be reinforced with large-scale simulations of polycrystals, which is the subject of the next section.

Effect of the Texture and Heterogeneous GB Properties during GG Simulations for a Polycrystalline Microstructure
In this section, we study a representative GB network in 2D. Figure 11 exhibits the initial characteristics of the microstructure, it consists of a square domain with length L = 1.6 mm and 5000 grains generated using a Laguerre-Voronoi tessellation [52] based on an optimized sphere packing algorithm [53] with a log-normal distribution for the arithmetic mean grain size. The grain size, R, of each grain is defined as √ S/π, with S as its surface (i.e., defined as the radius of the equivalent circular grain of the same surface). Anisotropic re-meshing is used following Equation (27), with a refinement close to the interface, the mesh size in the tangential direction (as well as far from the interface) is fixed at h t = 5 µm and at h n = 1 µm in the normal direction. The time step is fixed at ∆t = 10 s. This section is mainly devoted to studying the heterogeneity of both GB energy and mobility using the four introduced grain growth formulations. Finally, the same study is performed using a different texture in Appendix B.1.

Effect of the Heterogeneity
Here, we use a mis-orientation-dependent GB energy and mobility, defined respectively with a Read-Shockley (RS) function [54] and a Sigmoidal (S) function, proposed by Humphreys in [55]: where θ is the disorientation, µ max and γ max are the maximal GB mobility and energy respectively, and θ 0 = 30 • is the disorientation defining the transition from a low-angle grain boundary (LAGB) to a high-angle grain boundary (HAGB). θ 0 is normally considered to be between 15 and 20 • , but here, this parameter is exaggerated to exacerbate the heterogeneity of the system. The maximal values for the GB properties are µ max = 1.379 mm 4 J −1 s −1 , and γ max = 6 × 10 −7 J mm −2 , and are typical for a stainless steel [56]. Figure 11. Initial microstructure (a) with 5000 grains and the grain size distribution, in number (b). Figure 12 shows the orientation field using the vector magnitude O G = ϕ 2 1 + φ 2 + ϕ 2 2 , where (ϕ 1 , φ, ϕ 2 ) are the three Euler angles. The Euler angles defining the crystallographic orientations generated in this case are generated randomly, leading to a Mackenzie-like disorientation distribution function [57]. As the Read-Shockley model is used to define γ, the GB energy is concentrated at high values, as illustrated in Figure 13.

Heterogeneous Grain Boundary Energy
In this section, GB energy is defined using Equation (28) and GB mobility is assumed isotropic. Hence, the Het, HetGrad, and Aniso formulations are presented as "Het(µ:Iso)", "HetGrad(µ:Iso)", and "Aniso(µ:Iso)". The results are summarized in Figures 14-16. First, it is noticeable that all the formulations have a similar evolution concerning the total grain boundary energy E Γ , the number of grains N g , and the mean grain size weighted by numberR Nb [%] or by surfaceR S [%] . Additionally, if the grain size distribution weighted by number is normalized (Figure 15), one can recognize that all the formulations have similar distributions and the minima have similar values with respect to the mean radius. Similar results for the "Iso" and "HetGrad" formulations with heterogeneous GB energy defined by the Read-Shockley model were already reported [19].
The slow evolution of the mean values has been reported as a consequence of little local-heterogeneity produced by a Mackenzie-like DDF and/or a low value of θ 0 [31,35,58,59]. If the DDF starts as a Mackenzie distribution, the value of GB mobility and energy is focused at higher values, thus the microstructure cannot easily find a path to minimize its energy faster and the DDF changes slightly from its initial Mackenzie form. In other words, the initial configuration is almost isotropic. Slight differences can be observed after t = 1 h for the different formulations, which may be due to the low final number of grains (N G ≈ 500).
Regarding the morphology of the microstructures at t = 1 h, the grains are equiaxed. If we divide the total group of grains in classes divided by the number of neighbors (defined as the coordination number in the following), n, an interesting analysis regarding the morphology of grains could be performed. In Figure 16, the contribution of every class is depicted, and at t = 0 s, most of the grains verify n = 5. After one hour, one can directly appreciate that the class with n = 6 is the main class using the four formulations. This agrees with theoretical predictions of grain boundary motion with isotropic GB energy, which promotes triple junctions with dihedral angles near 120 • [60]. This aspect again illustrates the limited impact of the considered anisotropy in this configuration. (d) Aniso (µ:Iso) Figure 16. Grain size distribution and contributions from every group of grains of the same coordination number from 3 to 9 at t = 1 h.

Heterogeneous Grain Boundary Energy and Mobility
In this section, both GB energy and mobility are heterogeneous, respectively defined with Equations (28) and (29), and for that reason, the names introduced above are replaced by "Het (µ:S)", "HetGrad (µ:S)", and "Aniso (µ:S)". In order to compare the results presented above, the same initial microstructure and crystallographic orientations are used. The mean values' evolution and distributions remain similar among the four formulations and retain similar values as presented before. The heterogeneous GB mobility may affect the morphology of the microstructure due to a retarding effect from boundaries with disorientation lower than θ 0 . There is similarity between the four microstructures shown in Figure 17, showing mostly equiaxed grains. Two important aspects of these microstructures are that the microstructure obtained by the "Het" formulation is the most dissimilar, with a lower number of boundaries with disorientation inferior to θ 0 . Second, the presence of lowangle boundaries (θ < 30 • ) looks higher using the Anisotropic formulation. Nevertheless, this is not reflected in the interfacial energy evolution nor the DDF (see Figure 18). Finally, Figure 19 shows the disorientation distribution function using both an isotropic and heterogeneous mobility at t = 1 h. As stated before, the initial Mackenzie-type distribution evolves slowly, and a slow preference of low-angle boundaries is found. Using heterogeneous mobility slightly affects the DDF, and one can see that the Anisotropic formulation (Aniso (µ:S)) exacerbates low values of disorientation reflected in higher values in the distribution at 0 < θ < 10 • . Due to the Mackenzie-like DDF, the GB energy distribution is concentrated around γ max , leading to microstructures with triple junction angles around 120 • (see Figure 17). These results are in accordance with prior works [58,60,61]. (b) µ sigmoidal Figure 19. Disorientation distribution function at t = 1h using an isotropic (a) and a heterogeneous (b) mobility.
At this point, one can see that for a non-textured polycrystal with an initial Mackenzielike DDF, the evolution of the GB network and of the GB energy and mobility fields are similar to an isotropic case. That is the fundamental reason explaining the weak differences among the results of the different formulations. The results exhibit similar evolution of mean values, distributions, and grains' morphology. In order to study the behavior of the different formulations for a wider spectrum of GB properties, Appendix B.1 is devoted to studying the effect of a strong texture using the four formulations with isotropic and heterogeneous mobility. Under the effect of a textured microstructure, the Anisotropic formulation seems to be the more physical by promoting a higher percentage of boundaries with lower values of disorientation.

CPU Time
All the performed polycrystalline simulations were considered on 20 cores with the same mesh size, h n = 1 µm, in the normal direction of the interface, and h t = 5 µm in the tangential direction of the interface and far from the interface. As expressed earlier, both heterogeneous formulations and the Anisotropic formulation have additional terms which can be synonymous of more complex resolutions. This aspect, if not significant when moderated by anisotropy, is considered, as illustrated in the first line of Table 1.
However, the CPU time changes significantly for the textured case presented in Appendix B.1. The HetGrad and Aniso formulations present, respectively, an increase of 35% and of 74% of the calculation time, in comparison to the Isotropic formulation.

Accounting for Mis-Orientation and Inclination Dependence
The formulations presented above have dealt with heterogeneous GB properties. However, we know that the nature of the GB is described in a 5D space generated by the inclination and the mis-orientation. The effect of the normal direction has been described by Herring in [62] as a torque term. Hence, a triple junction should respect a condition frequently known as Herring's equation, Equation (20).
Due to the high dimensional space of GBs, many researchers have attempted to propose metrics that properly represent symmetries [63][64][65][66][67][68][69]. With these metrics, one can compare and compute the shortest paths (geodesics) between GBs. As the mis-orientation and the normal orientation can change during the microstructure evolution due to grain rotation or grain disappearance/appearance, the evolution of the metric could reveal important information about the structure-property relationship. Recent works by Chesser et al. and Francis et al. have proposed new metrics using octonions [70,71], revealing good predictions of GB energy of the data published by Olmsted in [72]. To the authors' knowledge, the effect of the GB normal orientation is not clear, and more experimental, numerical, and theoretical works are needed. Here, we define the effect of the normal orientation using a model of GB energy proposed for fcc metals by Bulatov et al. [73] and available in the GB5DOF code. When γ is defined using the GB5DOF code, both the effect of the mis-orientation and inclination are taken into account using the crystallographic orientations of the two adjacent grains and the local coordinate system of the corresponding GB [73].

Triple Junction
This case again consists of a triple junction, as described by Figure 20. We performed simulations with a constant GB mobility set to µ = 1 × 10 6 mm 4 J −1 s −1 , taken from [74], a domain of 1 × 1 mm 2 , and a time step of ∆t = 5 × 10 −5 s. The Aniso formulation is used by considering γ as only initially defined by the mis-orientation and then also dependent on the inclination (obtained through the GB5DOF code and denoted as Aniso-GB5DOF). The Iso, Het, and HetGrad are not presented here because they evolve in the wrong direction (the expected movement should reduce the length of the interface between grain G 1 and G 2 , depicted in yellow). The evolution of the interfaces shown in Figure 21 presents similar tendencies to the cases presented by Garcke in [75] and Hallberg in [34]. If both evolutions (without or with the inclination dependence) seem to promote similar triple junction evolution, the Aniso-GB5DOF case exhibits a much faster evolution, which illustrates the importance of accounting for the inclination in the reduced mobility description.

Coherent and Incoherent Twin Boundary
The main advantage of the GB5DOF code is that it is possible to characterize coherent and incoherent twin boundaries. These special GBs play an important role on polycrystalline microstructures, and their modeling is not frequently discussed at the mesoscopic scale. The next example was firstly proposed by Brown and Ghoniem in [76] and also reproduced at the mesoscopic scale in [34]. It consists of two grains composed of two coherent twin boundaries (CTB) and one incoherent boundary (ICB). Figure 22 shows the crystallographic orientation, the initial GB energy, and the variation of the GB energy as a function of the GB inclination. The Iso and Aniso formulations were used to model the GB movement. For the Aniso formulation, the GB5DOF code was used to compute the GB energy all along the simulation. On the other hand, the GB energy of the Iso case is constant and set to γ = 0.65969 J m −2 . The evolution of the GB is shown in Figure 23. The time step was set to ∆t = 0.1 ns and GB mobility was set to µ = 1.3 × 10 7 µm 4 J −1 ns −1 in order to reproduce the velocity of the ICB found by Brown and Ghoniem in [76], v ICB = 1.2 m s −1 . The movement of the ICB should be uniform and it should respect the flatness of the CTB. The Aniso-GB5DOF simulation enables to respect the expected behavior.

Conclusions
Different FE-LS formulations to simulate grain growth were presented and compared in this text, with the isotropic formulation, in which the grain boundary mobility and energy are assumed constants, being the most used framework in the literature. The isotropic formulation is able to reproduce mean grain size and grain size distribution evolutions when a moderated anisotropy is involved.
From the results presented using the triple junction cases, the Anisotropic formulation was the most accurate. The triple junction velocity predictions were the closest to the theoretical values while predicting accurate dihedral angles. In addition, the interfacial energy was always minimized and faster compared to the other approaches.
Additionally to these academic configurations, simulations using two different polycrystalline microstructures were performed. First, the initial orientations were generated using a uniform distribution, producing an initial Mackenzie-like disorientation distribution. Finally, another example with a textured orientation was considered. It was then illustrated that for a simple microstructure with initial random orientation, an isotropic formulation can be used, and that for a textured configuration, the Anisotropic formulation presents the best behavior in terms of grain morphology, DDF, and interfacial energy evolution predictions, while keeping a reasonable efficiency, compared to the isotropic formulation.
It was also illustrated that the Anisotropic approach is the most versatile approach, enabling to take into account the inclination dependence. Future works will be focused on the use of 2D and 3D experimental results concerning 304L and 316L, which are currently capitalized. These experimental results will be used to validate the Anisotropic formulation with more complex datasets.

Data Availability Statement:
The raw data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study. The processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A. Grim Reaper Case: Effect of the Boundary Conditions
In [26], the authors proposed the HetGrad formulation and performed several simulations for different values of r. The authors compared the dihedral angles against the analytical Grim Reaper values, see Equation (23), and found a very good estimation of the dihedral angles. A triangular domain was used with an initial triple junction equilibrium at 120°and Dirichlet boundary conditions (fixing the GB in the border domain). In other words, a final configuration respecting the Young's equilibrium is attained without the possibility to describe the transient state with an analytical solution. In order to study the Aniso formulation behavior, the same case is presented here. An isotropic mesh is used with a local adaptation around the triple junction where the mesh is refined in a circle of radius ε = 0.05, allowing the simulation to be more computationally efficient in terms of CPU time and memory storage. Figure A1 illustrates the mesh around the triple junction, where one can see the change of the mesh size close to the triple junction. Figure A1. Initial configuration of the triangular case with a refined isotropic mesh around the triple junction and a coarse mesh outside the triple junction with r = 10.
Multiple simulations were carried in order to study the effect of r. The constant parameters are the GB mobility µ = 1, the GB energy of the top interfaces γ top = 0.1, the mesh size at the triple junction h T J = 0.001, the mesh size outside the triple junction h = 0.01, and the time step ∆t = 1 × 10 −4 . As for the case presented before, the GB energy of the bottom interface is changed to obtain r ∈ {0.55, 0.625, 0.714, 0.833, 1.0, 1.25, 1.66, 2.5, 5, 10} (γ bot ∈ {0.18, 0.16, 0.14, 0.12, 0.1, 0.08, 0.06, 0.04, 0.02, 0.01}). In Figure A2, one can see the same tendencies as in Figure 10, with the HetGrad formulation being the best option in terms of dihedral angles' prediction.  Figure A3 shows the interface evolution with r = 10. The evolution is similar for the Grim Reaper example in Figure 8. The Het and Aniso formulations exhibit a Grim Reaperlike profile, while the HetGrad formulation evolves in the upward direction. This may seem wrong, however, for this particular geometry an upward movement is expected for r > 1 in order to match the analytical angles and as the initial angles are fixed to 120°. Thus, one can say that the interface obtained with the Het formulation evolves in the wrong direction. On the other hand, the movement obtained by the HetGrad formulation exaggerates the expected displacements and the interface is highly curved. Another illustration of the interface movement is shown in Figure A4 for r = 1.66, where the HetGrad and Aniso formulations have a correct evolution of the interface and the dihedral angle is closer to the analytical value, as shown in Figure A2.
In Figure A5, one can see the evolution of E Γ . The trends are similar to the previous test case. The HetGrad and Aniso formulations have a better energetic behavior, and the Aniso formulation remains the best option for high anisotropy levels.

Appendix B. Large-Scale Simulation: Effect of a Strong Texture
Here, the crystallographic orientations are defined differently: one Euler angle is generated randomly with a uniform distribution function and the two others are constants. As a result, the final disorientation distribution is more uniform, as seen in Figure A6. Properties are defined using Equations (28) and (29), and the transition disorientation angle is set to 30 • , as previously used. The main effect of the wider resulting GB energy distribution (GBED) is the increase of local anisotropy at triple junctions, as illustrated in Figure A6, compared to the previous test case (Mackenzie-like DDF).

Appendix B.1. Heterogeneous Grain Boundary Energy
The results described in Figure A7 illustrate that the Iso formulation predicts the fastest evolution. Additionally, one can see that the interfacial energy is better minimized using the Aniso formulation. From these results, one can infer that the isotropic formulation seems not adapted in this context. For a wider range of anisotropy levels, such as the one used in this test case, a particular coordination number with n = 4, 5 may be more present [60,61]. However, the Iso formulation promotes equiaxed grains (n = 6). Once again, this tendency discredits the Isotropic approach for highly heterogeneous interfaces.
Regarding both heterogeneous formulations (Het and HetGrad), the evolution of mean values and distributions are similar, as illustrated in Figures A7 and A8. First, both predicted distributions have similar groups, with n = 4, 5, 6, and second, the predicted microstructures show mostly equiaxed grain with a similar distribution of GB disorientation. In Figure A9, one can see similar clusters of GBs, with high values of disorientation depicted in red. From the morphology of grain boundaries ( Figure A9), the formulation that respects the most, on average, the triple junction angles is the Anisotropic one. This is illustrated in Figure A10, where the dihedral angles of a triple junction formed by GBs with low and high disorientation angles are shown. For this particular example, in Figure A10, blue and red boundaries have values of γ of about 0.25 × 10 −7 J mm −2 and 6 × 10 −7 J mm −2 , respectively. One can estimate an approximated value of the dihedral angle opposite the blue interface using Equation (23), which is about 177 • , with r = 6/0.25 = 24. The results described in Figures A7 and A11 show that while promoting a slower evolution of the microstructure, the Aniso (µ:Iso) formulation exhibits a better behaviour concerning the decreasing GB energy.  (d) Aniso (µ:Iso) Figure A8. Grain size distribution and contributions from every group of grains of the same coordination number from 3 to 9 at t = 1 h. Figure A9. Disorientation of the boundaries using the four formulations with homogeneous grain boundary mobility at t = 1 h. Boundaries with a disorientation higher than 30 • are colored in red. Figure A10. TJ dihedral angles among boundaries with high (red) and low (blue) GB energy. The disorientation of the boundaries is also depicted for the four formulations and homogeneous grain boundary mobility at t = 2 h. Boundaries with a disorientation higher than 30 • are colored in red.

Appendix B.2. Heterogeneous Grain Boundary Energy and Mobility
If heterogeneous GB mobility is added, the evolution of the microstructures can vary significantly. The results presented in Figure A12 show two regimes for the Het(µ:S) formulation. First, one can infer that the Het formulation presents issues to reduce the interfacial energy and presents a peak which is the result of an evolution dominated by curvature flow, without any effect of the heterogeneity. If we compare the results shown in Figures A7 and A12, one can see the retarding effect of using a heterogeneous GB mobility. This effect is stronger on the HetGrad and Anisotropic formulations due to the gradients introduced by heterogeneous fields, and is completely natural because technically, the effect of the crystallography is taken into account twice in the µγ product. Thus, the isotropic case evolves faster than the other formulations.
The results presented in Figures A13 and A14 show that, at t = 1 h, the Het and HetGrad formulations present more grains within the classes n = 4 and n = 7. On the other hand, the Anisotropic case did not evolve enough to compare it to the other cases.
Interestingly, the DDF of the Het formulation disagrees with the results presented in [35]. Here, the evolution of the DDF evolves in the opposite way to the expected results (see Figure A15). Indeed, the DDF tends to increase the percentage of interfaces with θ > θ 0 and decrease those with θ < θ 0 , which clearly seems nonphysical. Moreover, the Iso and HetGrad formulations do not exacerbate a particular disorientation. Finally, the Anisotropic formulation seems to exhibit a more physical behavior by promoting a higher percentage of boundaries with lower values of disorientation. (d) Aniso (µ:S) Figure A13. Equivalent radius distribution and contribution from every group of grains with coordination number from 3 to 9 at t = 1 h. Figure A14. Disorientation of the boundaries using the four formulations with heterogeneous grain boundary mobility at t = 1 h. Boundaries with a disorientation higher than 30 • are colored in red. (b) Distribution of γ Figure A15. Grain boundary characteristics' distributions at t = 1 h.