Level-Set Modeling of Grain Growth in 316L Stainless Steel under Different Assumptions Regarding Grain Boundary Properties

Two finite element level-set (FE-LS) formulations are compared for the modeling of grain growth of 316L stainless steel in terms of grain size, mean values, and histograms. Two kinds of microstructures are considered: some are generated statistically from EBSD maps, and the others are generated by the immersion of EBSD data in the FE formulation. Grain boundary (GB) mobility is heterogeneously defined as a function of the GB disorientation. On the other hand, GB energy is considered as heterogeneous or anisotropic, which are, respectively, defined as a function of the disorientation and both the GB misorientation and the GB inclination. In terms of mean grain size value and grain size distribution (GSD), both formulations provide similar responses. However, the anisotropic formulation better respects the experimental disorientation distribution function (DDF) and predicts more realistic grain morphologies. It was also found that the heterogeneous GB mobility described with a sigmoidal function only affects the DDF and the morphology of grains. Thus, a slower evolution of twin boundaries (TBs) is perceived.


Introduction
As most metallic materials exist in the form of polycrystals, determining the kinetics of metallurgical mechanisms such as recovery, grain growth (GG), and recrystallization is crucial, since they determine the final microstructure and properties [1]. Grain boundary (GB) engineering refers to the control of GBs with the aim of obtaining high-performance materials. Thus, numerical models have emerged to help us predict the evolution of microstructures submitted to different thermomechanical loads and the microstructureproperty relationship.
The migration of GBs is classically described at the polycrystalline scale by the wellknown equation v = µP, where v is the GB velocity, µ is the GB mobility, and P is the driving pressure. During GG, the evolution of GBs is driven by the reduction of interfacial energy, and the driving pressure is classically defined as a curvature flow driving pressure P = −γκ, where γ is the GB energy and κ is the mean curvature (i.e., the trace of the curvature tensor in 3D). At the polycrystalline scale, this kinematic equation is widely accepted while largely questioned [2,3]. Moreover, a definition of the reduced mobility (µγ) within the misorientation and inclination 5D space is not straightforward [4][5][6][7].
GG has been widely studied at the polycrystalline scale with various numerical approaches such as Phase-Field [8][9][10], Monte Carlo [11,12], Molecular Dynamics [13], Orientated Tessellation Updating Method [14], Vertex [15], Front-Tracking Lagrangian or Eulerian Section 2, crystallographic definitions, LS treatments, and FE-LS formulations are presented briefly. The methodology to estimate the GB reduced mobility from experimental data is presented in Section 3. Finally, the results using the isotropic and anisotropic formulations are compared using statistically representative Laguerre-Voronoï tessellations (Section 4), immersed microstructures with heterogeneous GB properties (Section 5), and immersed microstructures with anisotropic GB energy using the GB5DOF code [6] (Section 6).

The Numerical Formulation
The LS method was firstly proposed in [62] to describe curvature flow problems, and it was enhanced later for evolving multiple junctions [63,64] and applied to recrystallization and grain growth in [19,57]. The principle for modeling polycrystals is the following: grains are defined by LS functions φ in the space Ω φ(X) = ±d(X, Γ), X ∈ Ω, Γ = ∂G φ(X ∈ Ω) = 0 More precisely, the grain interface Γ is described by the zero-isovalue of the corresponding φ function. In Equation (1), d is the signed Euclidean distance to Γ and φ is classically chosen as positive inside the grain and negative outside. The dynamics of the interface is studied by following the evolution of the LS field. When the interface evolution is characterized by a velocity field v, its movement can be obtained through the resolution of the following transport equation [62]: Classically, one LS function is used to describe one grain and Equation (2) is solved for each grain to describe the grain boundary network evolution. However, when the number of grains, N G , increases, one may use a graph coloring/recoloring strategy [59] in order to limit drastically the number of involved LS functions Φ = {φ i , i = 1, . . . , N} with N N G . Two more treatments are necessary. Firstly, the LS functions are generally reinitialized at each time step in order to keep their initial metric property when they are initially built as distance functions to the grain interface, as proposed in Equation (1): In the proposed numerical framework, the algorithm developed in [65] is used. Secondly, the LS evolutions may not preserve the impenetrability/overlapping constraints leading to potential overlaps/voids between grain interfaces at multiple junctions. The solution proposed in [63] and largely used in an LS context [19] is adopted.
The main interest of this global numerical front-capturing framework lies in its ability to define different physical phenomena when they are encapsulated in the velocity field and to deal easily with topological events such as grain disappearance. In the next section, different formulations of the GB velocity and the subsequent FE resolution are presented.

GB Velocity Formulation
The isotropic formulation uses a homogeneous GB energy and mobility [19]; thus, the velocity field is defined as where κ is the mean curvature of the boundary in 2D and the trace of the curvature tensor in 3D, and n is the outward unit normal to the boundary. By verifying Equation (3) and assuming the LS function to be positive inside the corresponding grains and negative outside, the unitary normal and so the mean curvature may be defined as n = − ∇φ and κ = ∇ · n = −∆φ, and then, the velocity in Equation (4) may be written: At the mesoscopic scale, a GB, B ij , between grains G i and G j , is characterized by its morphology and its crystallographic properties, which may be summarized by a tuple B ij = (M ij , n ij ) with two shape parameters describing the interfaces through the unitaryoutward normal direction n ij and three crystallographic parameters describing the orientation relationship between the two adjacent grains known as the misorientation tensor M ij (see Figure 1). The misorientation is frequently defined with the axis-angle parameterization, i.e., M ij (a i , θ), where a i is the misorientation axis and θ is the disorientation [66]. Then, the two quantities of interest, the GB energy γ and GB mobility µ, must be seen as functions from the GB space B to R + . The anisotropic formulation was initially developed using thermodynamics and differential geometry in [37] and was improved in [22] in order to consider heterogeneous GB mobility. Both the GB normal and misorientation are taken into account, and an intrinsic torque term is present: v = µ(M) P ∇γ(M, n) · ∇φ − ∇ n ∇ n γ(M, n) + γ(M, n)I : K ∇φ (6) where I is the unitary matrix, P = I − n ⊗ n is the tangential projection tensor, ∇ n is the surface gradient, and K = ∇ n = ∇ ∇φ is the curvature tensor. The term Γ(M, n) = ∇ n ∇ n γ + γI is a tensorial diffusion term known as GB stiffness tensor [67,68]. The term P ∇γ · ∇φ in Equation (6) should be null in the grain interfaces. However, the frontcapturing nature of the LS approach, which involves solving Equation (2) at the GB network and in its vicinity, requires considering this term, which could be non-null around the interfaces. Then, this stabilization term is totally correlated to the front-capturing nature of the LS approach. The term Γ(M, n) is the subject of recent studies for twin boundaries (TBs) Σ3, Σ5, Σ7, Σ9 and Σ11 [67,69]. In the present work, the torque term ∇ n ∇ n γ is neglected, but the GB energy still depends on the GB misorientation and inclination, and the kinetic equation could be simplified as: Inserting the kinetic Equations (5) and (7) into Equation (2) leads to the weak formulation of the isotropic (Iso) and anisotropic (Aniso) formulation [22] Ω ∂φ ∂t ϕdΩ + Ω µγ ∇ϕ · ∇φdΩ − ∂Ω µγϕ ∇φ · nd(∂Ω) = 0, and respectively. If the properties are homogeneous, then both formulations are equivalent. The main question is the capability of these two numerical models to reproduce experimental evolution assessed through EBSD data. The next section is dedicated to the parameters identification.

Material Characterisation
The chemical composition of the 316L stainless steel is reported in Table 1. The samples were machined in the form of rectangular parallelepipeds of 8.5 mm × 8.5 mm × 12 mm. Then, the samples were annealed at 1050 • C during 30 min, 1 h, and 2 h. Afterwards, the samples were prepared for EBSD characterization. The preparation consisted of mechanical polishing followed by fine polishing and finally electrolitic polishing; the details of the polishing are listed in Table 2. Microstructures were analyzed at the center of the sample using a TESCAN FERA 3 Field Emission Gun Scanning Electron Microscope (FEGSEM). It is equipped with several detectors including Symmetry and C-Nano EBSD detectors from the Oxford company. In this work, the Symmetry EBSD detector was used. Post-processing was conducted using the MTEX toolbox in a MATLAB environment [70]. The EBSD map at t = 0 h has a size of 1.138 mm × 0.856 mm and was acquired with a constant step size of 1.5 µm. The other three EBSD maps at t = 30 min, 1 h, and 2 h have a size of 1.518 mm × 1.142 mm and were acquired with a constant step size of 2 µm. Grain boundaries have a disorientation above 5 degrees (θ > 5 • ) and Σ3 twin boundaries have a misorientation axis <111> ±5 • and θ = 60 ± 5 • .
The main properties of the initial microstructure are reported in Figure 2. Figure 2b illustrates the grain size and disorientation distribution ignoring Σ3 TBs. The grain size is defined as an equivalent radius, R = √ S/π, where S is the grain area. The microstructure consists of equiaxed grains with an arithmetic mean radius of 15 µm, and few bigger grains with a radius around 60 µm. Additionally, the microstructure presents a Mackenzie-like DDF typical of random grain orientations [71]. On the other hand, if Σ3 TBs are considered, the DDF presents an additional sharp peak at a disorientation angle θ = 60 • , which comes from the TBs and then constitutes a strong source of anisotropy with regard to GB properties (see Figure 2c).   Figures 3 and 4, the evolution of the microstructure seems to mostly proceed by normal grain growth (NGG), but the surface grain size distribution shows that the microstrucuture has a bimodal population of grains (see Figure 5). However, some of the grains can reach an equivalent diameter above 0.1 mm, which is much larger than the average grain size.

Estimation of the Average Grain Boundary Mobility Based on the Burke and Turnbull GG Method
In order to compute the average mobility necessary to run full-field simulations, the evolution of the arithmetic mean grain radiusR Nb must be known. Figure 6 shows the evolution ofR Nb as a function of the annealing time. Using the methodology discussed in [72,73], one can obtain an average reduced mobility µγ using the Burke and Turnbull model [74]. This model, where topological and neighboring effects are neglected, is based on five main assumptions: the driving pressure is proportional to the mean curvature, grains are equixaed, the GB mobility and energy are isotropic, the annealing temperature is constant, and no second-phase particles are present in the material. In this context, one can obtain a simplified equation describing the mean radius evolution: This methodology has been used in [72,73,[75][76][77] assuming general grain boundaries with homogeneous GB energy and mobility. From the evolution ofR Nb in Figure 6 (excluding the Σ3 TBs), one can then obtain a first approximation of the product µγ for the general boundaries at 1050 • C. This approximation will be used for the µγ definition in isotropic simulations. Nevertheless, as illustrated by the second orange curve in Figure 6, when Σ3 TBs are considered in the analysis, grains are of course smaller, but they also grow much slower, with a direct impact on the apparent reduced mobility. This slow evolution can be produced by the strong anisotropy brought by special GBs in the global grain boundary network migration. Different ways to improve the description of the reduced mobility and their impacts in the results are discussed in the following.
'woS3' 'All' Figure 6. Mean grain radius evolution at 1050 • C from experimental data measured in 2D sections by taking into account all boundaries (in orange) and without TBs (in green). The outer grains that share a boundary with the image borders are not taken into account in the analysis.
In the following sections, general boundaries make reference to the case without Σ3 TBs, and the case with Σ3 is referred to as all boundaries.

Statistical Cases
In this section, the 2D GB network is created from the initial experimental grain size distribution shown in Figure 7. The square domain has a length L = 2.0 mm, and grains are generated using a Laguerre-Voronoi tessellation [61] based on an optimized sphere packing algorithm [78]. Anisotropic remeshing is used with a refinement close to the interfaces, the mesh size in the tangential direction (and far from the interface) is set to h max = 5 µm and in the normal direction h min = 1 µm, with transition distances set to φ min = 1.2 µm and φ max = 5 µm (see [22,58,79] for more details concerning the remeshing procedure and parameters). The time step is set to ∆t = 10 s. The orientation field was generated randomly from the grain orientations measured by EBSD in the initial microstructure ( Figure 2a). The first part studies grain boundaries without Σ3 TB, (Section 4.1). In the second part, Σ3 TBs are included in the analysis (Section 4.2). The interfacial energy and average GB properties are computed as where T is the set of elements in the FE mesh, l e is the length of the LS zero iso-values existing in the element e, i refers to the number of LS functions, L Γ is the total length of the GB network Γ, and x e is the GB property of the element e.

Statistical Case with General Boundaries
The first case with general boundaries is composed of N G = 4397 initial grains. Figures 8a,b show the initial GB disorientation and the initial DDF distribution. Most of the interfaces have a disorientation higher than 15 • due to the random generation of orientations that leads to a Mackenzie-like DDF [71]; see Figure 8b. Then, the GB energy and mobility are defined as being disorientation dependent, using a Read-Shockley (RS) [55] and a Sigmoidal (S) function [56], respectively: and where θ is the disorientation, µ max and γ max are the GB mobility and energy of general high-angle grain boundaries (HAGBs). θ 0 = 15 • is the disorientation defining the transition from a low-angle grain boundary (LAGB) to a HAGB. The maximal value of GB energy is set to γ max = 6 × 10 −7 J·mm −2 , which is typical for stainless steel [72,80]. The value of general HAGB mobility was computed using the methodology presented in Section 3.2 and is fixed at µ max = 0.476 mm 4 ·J −1 ·s −1 for both isotropic and anisotropic formulations. The simulations carried out using the anisotropic formulation consider heterogeneous GB energy defined by Equation (12) and two descriptions of the mobility. If GB mobility is isotropic, the formulation is referred as "Aniso(µ:Iso)", and in the cases where GB mobility is heterogeneous (i.e., defined by Equation (13), the formulation is referred to as "Aniso(µ:S)". Figure 9 shows the evolution of average quantities: normalized total GB energy E Γ /E Γ (t = 0), normalized number of grains N G /N G (t = 0), arithmetic mean grain radiusR Nb , and normalized average GB disorientationθ/θ(t = 0). One can see that the mean grain radius evolution agrees with the experimental data and that the evolution of the other mean values are close to each other when using the different formulations and present reasonable variations from the experimental data. As stated in [22], the effect of a heterogeneous GB mobility does not affect the evolution of the mean values and distributions when orientations are generated randomly, and the DDF is similar to a Mackenzie distribution. One can see that the only mean value affected by the GB mobility is the average GB disorientationθ, being the "Aniso(µ:S)" formulation the one that is closer to the experimental evolution.     Figure 9. Mean values time evolution for the isotropic (Iso) formulation, anisotropic formulations with isotropic GB mobility (Aniso(µ:Iso)) and heterogeneous GB mobility (Aniso(µ:S)) and the experimental data (EBSD). Numerical results obtained from the initial microstructure shown in Figure 8a. Figures 10 and 11 show a good match of GSD and DDF between simulation results and experimental data after one and two hours of annealing at 1050 • C. One can see that the three cases are alike. The initial Mackenzie-like DDF evolves slowly for cases with random orientations and low anisotropy (as in [22,35,[81][82][83]). Finally, in Figure 12, one can see the similarity between the microstructures obtained in the different simulations with most of the grains being equiaxed and few LAGBs.    Length Figure 11. Disorientation distribution obtained excluding TBs at (a) t = 1 h and (b) t = 2 h for the isotropic (Iso) formulation, anisotropic formulations with isotropic GB mobility (Aniso(µ:Iso)), and heterogeneous GB mobility (Aniso(µ:S)) and the experimental data (EBSD). The Y-axis represents the GB length percentage. The results presented here show that the evolution of an untextured polycrystal with an initial Mackenzie-like DDF could be simulated using an isotropic formulation or an anisotropic formulation with heterogeneous GB energy or both heterogeneous GB mobility and energy. This methodology has been used in different contexts under different annealing processes [72,[75][76][77] and with academic microstructures in [22]. In the next section, the same analysis is performed by considering the same domain but by introducing special grain boundaries through an update of the µ and γ fields.

Statistical Case with an Improved Description of the γ and µ Fields
The microstructure used in this section, described in Figure 13, was generated using the same domain with L = 2.0 mm and the GSD shown in Figure 7b. The initial number of grains is N G = 14,956 and their orientation is also generated randomly from the EBSD orientations. The initial number of grains is more important comparatively to the previous test case as the GSD described in Figure 2c, (left side) where Σ3 TBs are taken into account is used to generate the Laguerre-Voronoï polycrystal. The same mesh and time step as the previous simulations are used in order to be able to fairly compare the obtained results. In order to define the behavior of special GBs with properties close to Σ3 TBs, the µ and γ fields are updated as follows: with θ 0 = 15 • and θ Σ3 = 60 • and a value of GB energy and mobility set to γ max * 0.1 and µ max * 0.1 for GBs with θ ≥ θ Σ3 . The value of µ max is estimated using the evolution ofR Nb [%] considering all GBs (see Figure 6). The estimated GB mobility is µ max = 0.069 mm 4 ·J −1 ·s −1 , and it is one order of magnitude lower from the GB mobility estimated without Σ3 TBs. The decrease of the GB mobility is proportional to the decrease of grain size (see Figure 6) due to the high number of TBs. From the results shown in Figure 14, one can see that all formulations minimize the energy with the same efficiency ( Figure 14a) and the microstructures evolve at the same rate, leading to a good fit of mean grain size and number of grain evolutions comparatively to the experimental data. On the other hand, the anisotropic formulation shows a better agreement in terms of mean disorientation evolution.  'Aniso(µ:S)' 'EBSD' (d)θ/θ(t = 0) = f (t) Figure 14. Mean values time evolution for the isotropic (Iso) formulation, anisotropic formulations with isotropic GB mobility (Aniso(µ:Iso)) and heterogeneous GB mobility (Aniso(µ:S)), and the experimental data (EBSD). Numerical results obtained from the initial microstructure shown in Figure 13a. Figure 15 shows the evolution of the grain size distribution at t = 1 h and 2 h. GSDs present a good match ( Figure 15) with small differences between the Iso and anisotropic formulations. However the DDF predictions ( Figure 16) are quite bad for all formulations even if the anisotropic calculations tend to be better. This result can easily be explained by the use of statistics (GSD and orientations) from EBSD data which are not sufficient to accurately describe the real microstructure. First, the initial topology is simplified, but above all, even if the orientation data come from the EBSD measurements, the resulting initial DDF is not accurate, as a Mackenzie-like distribution is obtained, as illustrated in Figure 13b. Hence, the effect of the heterogeneous GB energy and the mobility anisotropy of the real microstructure are underestimated. A way of improvement of the proposed statistical generation methodology will be to modify the algorithm dedicated to the assignment of the orientation of each virtual grain by imposing also a better respect of the experimental DDF [84].
The GB energy of the microstructure at t = 2 h is shown in Figure 17, a higher number of blue GBs, which correspond to TBs (low value of γ), are obtained using the anisotropic formulations, and the effect of heterogeneous GB mobility seems negligible. The next simulations are carried out by immersing the EBSD data in order to overcome the limits discussed above.

Immersion of EBSD Data
In this section, a digital twin microstructure obtained by immersion of the EBSD map acquired on the initial microstructure ( Figure 2) is discussed. Figure 18 shows the band contrast map of the microstructure and its numerical twin. The dimensions of the domain are L x = 0.856 mm and L y = 1.138 mm, and it contains 3472 crystallites. The time step is fixed at ∆t= 10 s, and the domain is discretized here using an unstructured static triangular mesh with a mesh size of h= 1 µm. This microstructure is more appropriate to compare simulations and experimental data. The evolution of the numerical microstructure is compared to EBSD maps obtained at three different times: t = 30 min, 1 h, and 2 h (see Figure 3). The GB energy and mobility are defined using Equations (14) and (15), respectively. The maximal value of GB energy is set to γ max = 6 × 10 −7 J·mm −2 , and the maximal value of GB mobility is set to fit the mean grain size evolution. The maximal value of GB mobility for the Aniso(µ:Iso) and Aniso(µ:S) formulations are µ max = 0.146 mm 4 ·J −1 ·s −1 and µ max = 0.272 mm 4 ·J −1 ·s −1 , respectively. Regarding the isotropic formulation, the value of GB reduced mobility remains constant µγ = 0.414 × 10 −7 mm 2 ·s −1 . The changes in µ max are due to the more complex geometry and the higher quantity of special boundaries that produce additional gradients of GB energy and GB mobility (see Equation (9)). As stated before, the additional jump at θ Σ3 is set to define special boundaries similar to Σ3 TBs. Figure 19 confirms the good match between the TBs colored in red in the EBSD band contrast map (left side) and TBs colored in blue corresponding to a low GB energy in the numerical twin microstructure (right side). With the immersed data, one can obtain a close digital twin of the real microstructure with the initial GB distributions presented in Figure 20. The initial GB energy distribution (GBED) is shown in Figure 20b. With this particular distribution, several junctions will have a high anisotropy level, and as stated in [22,83], one can expect different behaviors using the different formulations. As illustrated in Figure 21, the three simulations predict similar trends concerning the mean grain size and the grain number evolution. Comparatively to experimental EBSD data, these predictions are very good concerning the mean grain size prediction, but all of them tend to predict, at the beginning, a faster disappearance of the small grains. Concerning the total energy, mean GB disorientation, and mean GB energy evolutions, the anisotropic formulation is closer to the EBSD data. This means that the Aniso formulation is more physical and promotes GBs with low GB energy. Figure 22 illustrates the topology of grains at t = 2 h. One can notice the higher quantity of GBs with low GB energy using the anisotropic formulation and its similarity to the EBSD band contrast map even if one can notice from the EBSD data that the real microstructure contains more TBs that create small grains as reflected in the GSD in Figure 23. Another advantage of the anisotropic formulation is the better reproduction of the DDF compared to the Iso formulation, which contains a lower percentage of GBs with θ ≈ θ Σ3 and tends to promote a Mackenzie-like DDF (see Figure 24). One can also see in  that the heterogeneous GB mobility improves the morphology of grain, the GSD, and the DDF.
In this section, it has been shown that the immersed data give a better insight of the real microstructure evolution. In terms of mean values, there are small differences between the results from the Laguerre-Voronoï tesselation and the immersed microstructure. However, the GSD and DDF distributions are better reproduced for the immersed case when the anisotropic formalism is adopted. The heterogeneous GB mobility affects the grain topology, the GSD, and the DDF due to its additional retarding effect. Regarding the GSD, the Aniso(µ:S) formulation can reproduce the peak at low values of grain size and the peak around θ Σ3 for the DDF, which are due to the TBs. Nevertheless, the behavior of TBs is still not perfectly reproduced by the proposed simulations, being the anisotropic formulation the one that seems more physical. In the next section, this issue is addressed using the GB5DOF code, which allows to define the GB energy in terms of misorientation and GB inclination (in 2D) in order to better characterize the evolution of TBs. (e)γ = f (t) Figure 21. Mean values time evolution for the isotropic (Iso) formulation, anisotropic formulations with isotropic GB mobility (Aniso(µ:Iso)) and heterogeneous GB mobility (Aniso(µ:S)), and the experimental data (EBSD). Numerical results obtained from the initial microstructure shown in Figure 18.  formulation, anisotropic formulations with isotropic GB mobility (Aniso(µ:Iso)) and heterogeneous GB mobility (Aniso(µ:S)), and the experimental data (EBSD); N G refers to the number. Numerical results obtained from the initial immersed microstructure shown in Figure 18 and the RS and Sigmoidal model to define GB energy and mobility. formulation, anisotropic formulations with isotropic GB mobility (Aniso(µ:Iso)) and heterogeneous GB mobility (Aniso(µ:S)), and the experimental data (EBSD). The Y-axis represents the GB length percentage. Numerical results obtained from the initial immersed microstructure shown in Figure 18 and the RS and sigmoidal model to define GB energy and mobility.

Simulation Results
In this section, the immersed polycrystalline microstructure and the FE mesh presented in Section 5 are used. Anisotropic GB energy values are defined using the GB5DOF code [6], and the heterogeneous GB mobility is described using Equation (15). It means that the GB energy can vary with the GB misorientation and inclination even if the torque terms are neglected. Note that the GB inclination is measured in 2D and not 3D; in other words, the GB is supposed to be perpendicular to the EBSD map. The GB energy of the microstructure and its GBED are shown in Figure 25. The initial microstructure is shown in Figure 25a, and one can see that the maximum value of GB energy is set around γ max ≈ 7 × 10 −7 J·mm −2 . As discussed in [6], incoherent Σ3 TBs have a GB energy defined as γ Σ3 ≈ 0.6 * γ max , meaning that the modified Read-Shockley model described by Equation (14) seems exaggerated. In Figure 25b, the GBED is concentrated within the values 4 × 10 −7 J·mm −2 ≤ γ ≤ 7 × 10 −7 J·mm −2 , which means that the level of heterogeneity is low and the different formulations are expected to promote similar trends, as stated in [22]. The GB mobility was set to fit the evolution of the mean grain size, the maximal GB mobility for the Aniso(µ:Iso) and Aniso(µ:S) formulations are, respectively, set to µ max = 0.0767 mm 4 ·J −1 ·s −1 and µ max = 0.1423 mm 4 ·J −1 ·s −1 . The difference between the µ max values is generated by the higher values of GB energy produced by the GB5DOF code; note that with the RS model, all the TBs are defined as coherent, while the GB5DOF code can distinguish between coherent twin boundaries and incoherent twin boundaries, as pointed out in [22,38]. First, the mean grain size and grain size distribution evolutions are well reproduced by the different simulations (see Figures 26 and 27). The mean GB disorientation is not well represented by any of the formulations, the EBSD data show a stable value around 50 • while all numerical results exhibit a decreasing trend (see Figure 26b). This effect is due to the TBs and illustrates the inability for the numerical formulations to preserve or generate them. Additionally, the DDF from both formulations are similar and do not correspond to the experimental DDF (see Figure 28). The similarity between the isotropic and anisotropic simulations is due to the low anisotropy level, which may be produced by the lack of information of the GB inclination (see Figure 29). As stated in [22], when the GBED is concentrated around a specific value, both formulations can present a similar trend. This is confirmed with Figure 30, where a zoom on the GB network is shown at four different times, and one cannot see any obvious difference among the obtained microstructures with the three different simulations. The main difference of these results lies in the ability of the anisotropic formulation to keep more Σ3 TBs when the sigmoid description of µ is used. formulation, anisotropic formulations with isotropic GB mobility (Aniso(µ:Iso)) and heterogeneous GB mobility (Aniso(µ:S)), and the experimental data (EBSD). The Y-axis represents the GB length percentage. Numerical results obtained from the initial immersed microstructure shown in Figure 25a and the GB5DOF code to define the GB energy.

Current State of the Modeling of 3D Anisotropic Grain Growth
A final question regarding the anisotropy of GB properties is still open: do the 3D descriptions of GB properties affect the microstructure evolution? Until now, most of the studies of GG in 3D have presented simulations of polycrystalline microstructures using different textures and a mathematical description of GB properties [33,[50][51][52] or using databases of GB energy [53,54]. The following conclusions are pointed out: • The effect of the heterogeneity is stronger when the material is textured or the disorientation transition between LAGBs and HAGBs, θ 0 , is high [51,52,54]; • The individual effect of GB energy and mobility is small on the GG [52].
Note that similar conclusions were presented in the first part of this work [22]. In [33,[50][51][52], GB properties are defined as heterogeneous and not as anisotropic. The inclination dependency can have an important impact; hence, a complete description of the GB properties is necessary; i.e., µ(M(θ, a), n) and γ (M(θ, a), n), as well as 3D non-destructive in situ characterization [45,47,48] in order to obtain more realistic values of GB mobility, must be considered as a key perspectives concerning the full-field modeling of GG.
In the simulations presented in this section, the GB inclination is simplified as it is projected in the observation plane. In other words, the description of the GB properties is simplified, and only a slice of the GB energy is considered. For instance, Figure 31 shows the GB energy and mobility of a Σ3 TB. One can see that the 3D surface of the TB properties has a complex geometry. On the other hand, the anisotropy of the GB mobility is simplified to a sigmoidal model with a cusp at θ Σ3 . Unluckily, the GB mobility data are not available for the complete GB space and at different temperatures.

Summary and Conclusions
Different FE-LS formulations to study the GG of 316L stainless steel were compared. The isotropic formulation is able to reproduce, for statistically generated or immersed polycrystals, the average grain size and grain size distribution for a wide range of anisotropy levels.
The results obtained using representative Laguerre-Voronoï polycrystals show that the heterogeneous GB mobility values do not affect the response of the different formulations and that the anisotropic formulation is more physical being the only formulation that enables promoting GBs with low energy. However, the anisotropy level was largely underestimated because of the initial Mackenzie-like DDF, which could not be controlled during the polycrystal generation.
Two additional cases were presented with a twin numerical microstructure obtained directly from EBSD data. The main advantage is that the initial DDF and topology are accurately defined. First, GB energy and mobility were defined using the modified Read-Shockley and sigmoidal models already tested for the virtual polycrystals. Then, the model was coupled with an anisotropic model of GB energy that takes into account the GB misorientation and inclination [6]. However, it is highlighted again that the GB inclination is not well defined, as the GB is supposed to be perpendicular to the observation plane. The proposed RS model seemed to exaggerate the anisotropy level comparatively to the GB5DOF code. However, the predictions are clearly better with the proposed RS and sigmoid models associated to the anisotropic formulation while not allowing to be predictive concerning the DDF whatever the method chosen.
These results illustrate that the prediction of grain growth at the polycrystal scale can be ambiguous depending the aimed attributes and the available data. First of all, 3D simulations should be considered. Of course, this aspect is, firstly, essential to improve the representation of the considered polycrystals but is also essential to describe correctly the γ dependence on the inclination. Indeed, the proposed 2D model/data context limits the actual use of the inclination, as this parameter is described here with one degree of freedom and not in a 3D framework with 3D experimental data. This aspect can explain the low anisotropy level obtained using the GB5DOF code. Finally, this objective must also be correlated to the fact to integrate the torque effects and thus the GB stiffness tensor in the simulations and analysis. It should be highlighted that this conclusion is common to all existing works of the state of the art involving anisotropic 2D GG simulations and 3D simulations where the inclination dependence, or torque terms, or both are not taken into account. Finally, it must also be emphasized that the influence of impurities segregation and/or oxides can be important for the considered alloy but was not discussed in the numerical framework; therefore, this constitutes another perspective of this work.

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

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

Abbreviations
The following abbreviations are used in this manuscript: