A 2D Front-Tracking Lagrangian Model for the Modeling of Anisotropic Grain Growth

Grain growth is a well-known and complex phenomenon occurring during annealing of all polycrystalline materials. Its numerical modeling is a complex task when anisotropy sources such as grain orientation and grain boundary inclination have to be taken into account. This article presents the application of a front-tracking methodology to the context of anisotropic grain boundary motion at the mesoscopic scale. The new formulation of boundary migration can take into account any source of anisotropy both at grain boundaries as well as at multiple junctions (MJs) (intersection point of three or more grain boundaries). Special attention is given to the decomposition of high-order MJs for which an algorithm is proposed based on local grain boundary energy minimisation. Numerical tests are provided using highly heterogeneous configurations, and comparisons with a recently developed Finite-Element Level-Set (FE-LS) approach are given. Finally, the computational performance of the model will be studied comparing the CPU-times obtained with the same model but in an isotropic context.


Introduction
Grain growth phenomenon in polycrystals has been studied for many decades, both from an experimental and numerical point of view [1]. The majority of experimental observations at this scale suggest that the migration of boundaries is, in general, a strongly heterogeneous phenomenon involving complex dynamics and topological transformations of the grain boundary (GB) network. However, in the literature, it is frequently accepted that the microstructure of given materials behave homogeneously enough to ignore their heterogeneities when considering polycrystal modelling. This hypothesis is used in numerical environments to propose Full-Field (FF) models of microstructural evolutions, using homogeneous values in space of the grain boundary energy γ and mobility µ, e.g., isotropic grain growth (GG).
If this hypothesis remains acceptable when low levels of anisotropy are involved, it constitutes, however, a strong approximation when a strong texture with particular γ values are involved, or when special GBs (e.g., twin boundaries) are present [2].
Commonly, in the literature [1], the source of GB anisotropy, i.e., the reduced mobility defined as the µγ product, is considered a function of the crystallographic misorientation and of the inclination of the interface. Typically, the misorientation M lw between two adjacent grains l and w, is computed using the three Euler angles (ϕ e1 , Φ, ϕ e2 ) of these grains and the inclination is considered through the local normal vector n of the corresponding GB. This gives a system with a total of 5 degrees of freedom (DOF) defining GB properties. These kinds of systems at the polycrystal scale need to be modelled through the use of a numerical approach able to take into account this kind of data set. As such, in the same manner as in [3][4][5], here we will differentiate three kinds of numerical models: isotropic, heterogeneous, and anisotropic models. Isotropic models consider constant GB properties in their formulation. On the contrary, Anisotropic models are those using a formulation where, any assumption regarding the invariability of these quantities in space is discarded, being able to use properties dependent on the tuple (M lw , n) (i.e., X (M lw , n) where X is either γ or µ). Of course, anisotropic models are much more complex than those using an isotropic hypothesis, since, in this context, special attention must be given, for example, to the meaning of the surface tension component of interfaces, as one must be aware that torque terms, derived from the variation of the GB energy γ on the parametric space of a surface may appear [6]. As such, deriving a mathematical model ready to use an anisotropic set of GBs properties is a complex task and, historically, authors in the literature have proposed alternatives: heterogeneous models. Heterogeneous models consider within their formulation the existence of a variation of properties, only in function of M lw (X (M lw ) ), neglecting its dependence on n. In this context, each GB is given homogenised intrinsic properties (constant in its parametric space), but different from the ones of other GBs. i.e., GB properties only change at multiple junctions (MJ) (or multiple lines in 3D) when crossing from one GB to another.
Several approaches have been proposed in the literature to model heterogeneous and/or anisotropic GG. Beginning with the Monte Carlo and extending to Phase-Field, Level-Set and Vertex approaches, heterogeneous (X (M lw ) ) [7][8][9][10], and anisotropic (X (M lw , n) ) [3][4][5]11,12] models have been proposed. However, all these methods are constrained by different reasons each, typically: (i.) the use of regular grids [13,14] (which can lead to difficulties to model large deformation), (ii.) high computational cost [4,5,8], and (iii.) the no-discretization of grain interiors [12,15] (which can lead to difficulties when intragranular phenomena are of interests). Additionally to these aspects, in the context of anisotropic boundary properties modelled using Phase-Field models, although being an appropriate numerical environment, showing interesting results in this context, one should be aware of inherent numerical instabilities, especially for high heterogeneous/anisotropic systems [10,16]. Finally, in anisotropic models, the GB energy dependence on the inclination is classically defined without inquiring if additional torque terms in solved equations are needed with the notable exception of the vertex approaches [12,15].
As an alternative to model microstructural evolutions with anisotropic GB properties, we propose the TRM model presented in [17][18][19]; this article will present the needed implementations in order to model fully anisotropic grain properties with the TRM model. Special attention will be given to the development of a robust high order Multiple Junction (MJ) decomposition algorithm and to the reformulation of the velocity equation at triple junctions extending the methodology presented in [15] to an anisotropic context, using the notions used in [12] for its discrete formulation. Finally, the TRM model will be tested in multiple heterogeneous environments identical to the ones presented in [8,9] while the numerical tests in a fully anisotropic environment will be discussed in a forthcoming publication [20].

Numerical Method
This section will introduce the TRM model's necessary concepts and new implementations to model GBM using anisotropic GB properties. The topological changes that may occur in this context have the same level of complexity as the ones produced under the influence of stored energy, presented in a previous work [19]. Additionally, in [17], the decomposition of high-order multiple junctions (MJs) was simplified for isotropic GB properties. A more developed algorithm is needed to obtain valid predictions in an anisotropic context. This section will cover these notions.
Hereafter, we will consider γ as a function of (M lw , n), while the mobility term µ will be considered constant in space. Before considering a misorientation in the computation of grain boundary properties, each grain requires an orientation. In this work, these orientations are generated at random using a uniform distribution for each of Euler's angles (ϕ 1 , Φ, ϕ 2 ). Figure A1 (bottom-right) gives un example of the disorientation angle distribution obtained with this approach and compares it to the Mackenzie distribution [21] for disorientation angles in a cubic sample.

Grain Boundary Motion by Capillarity: Anisotropic Context for the TRM Model
In [12], a formulation for the computation of the velocity of GBs and triple junctions using anisotropic GB properties was proposed in the context of the Vertex model. This formulation uses the tensile character of the capillarity forces exerted at every node based on a discrete analysis, similar to the one used in [19] for the computation of a velocity from a stored energy field at triple junctions. The model in [12] writes for the velocity at MJs: where the index i denotes the node representing the MJ N i and j their connection to node N j , µ i is the mobility of node N i , γ ij , t ij and n ij are, respectively, the boundary energy, the unit tangent vector and the normal of the segment N i N j . Note that γ ij = γ ji but t ij = − t ji and the direction of the normal n ij is arbitrary. Finally, note the apparition of the term τ ij , which corresponds to the torque experienced by the segment N i N j due to the change of the GB energy given by its dependence on the inclination angle ω [6]. This torque term is defined as follows: In [15], three formulations were given for the computation of the velocity at MJs in the context of isotropic GB properties, from which we have used the so-called model II to find our velocity at MJs in previous works [17,19]. This formulation can be rewritten in the context of heterogeneous grain boundary properties (hence, in the absence of torque terms) and for MJs of arbitrary order, in a very similar way as in Equation (1): where c is the number of connections of the MJ, and where the only difference with Equation (1) is that the terms in the numerator contribute all in the same amount to the summation, instead of being escalated each by the term |N i N j | −1 . Indeed, in our experience, the homogenization of the contributions of the numerator term by the separated summation ∑ j |N i N j | has proven to be more stable than the one given in Equation (1), especially when the value of any |N i N j | approaches to zero (or when |N i N j | << |N i N k | for all k = j).
For this reason, the use of Equation (3) is preferred here but maintaining the torque terms of Equation (1): Finally, note that torque effects also need to be considered at GBs. For this purpose, the analytical model introduced in [22] for single surfaces can be used: where n ix is the projection of the variable normal vector n, onto the tangent vector to the interface at node i. In practice, we have found that applying Equation (5) might produce oscillatory effects on the computation of velocity v c given by the second derivatives of γ.
To avoid such instabilities, we use a combination of the discrete approach given in [12] and the standard approach of the TRM model: where the terms κ i and n i are computed using the numerical approximation by splines at node i. Note that n i = n ij , as n ij denotes the normal of the segment N i N j and n i the normal to the numerical approximation at node i.

Minimal-State Energy of High-Order MJs
The grain growth phenomenon is driven by the minimization of GB energy. In an isotropic context, this means that a GB network will continuously reduce its total surface, producing numerous topological changes (interface destruction/creation and grain disappearance) in their structure over time. In [17], it was explained how, when using anisotropic GB properties, the topological changes can diverge from the isotropic behaviour, mainly during the decomposition of multiple junctions occurring after an interface destruction. This section provides an insight of the decomposition of high-order MJs when considering anisotropic GB properties.
The main challenge here is to explore all possible configurations that may proceed after a decomposition process. The size of the possibilities set P(z) is only dependent on the MJ's order z to be decomposed. A fourth-order MJ (i.e., four grain boundaries meeting in a point) can be decomposed only in two ways. However, the set of possibilities P(z) grows much higher when the MJ's order increases. Consider the configuration given in Figure 1, here we provide a fifth-order MJ, as well as the firsts five possible decompositions given by the separation of consecutive interface pairs. Note, however, that each possibility regroups one fourth-order and one third-order MJ, from which the fourth-order one might decompose in two third-order MJs. In total, for a fifth-order MJ, five final possible decompositions are allowed when decomposing all MJ with z > 3 (P 3 (5) = 5 where the upper script means that all final MJs are z = 3, see Figure 2). P 3 (z) increase rapidly with the MJ's order z: for z = (2, 3, 4, 5, 6, 7, 8, . . . ), P 3 (z) = (1, 1, 2, 5, 14, 42, 132, . . . ). In general, the number of possible combinations in this context is given by the Catalan numbers [23] formula C n : Of course, the probability of encountering a MJ of order z decreases as z increases, as for a MJ of order z to form, it would require that all P (z−1) possible decompositions were stable. This notion of stability is related to the total minimum energy state able to be reproduced for a given initial configuration. Note that this notion also suggests that one could obtain a total minimal energy state for a MJ with z > 3, in which case this MJ should not be decomposed [12]. As such, not only the configurations given by P 3 (z) need to be considered, but also those in between (e.g., the ones given in Figure 1), hence producing P (z) >> P 3 (z) . We have simplified this aspect by accepting configurations presenting local minimal energy states and by not testing all possible configurations P (z) . Details of the decomposition algorithm are given below.

Algorithm for the Decomposition of High-Order Multiple Junctions
Algorithm 1 summarizes the TRM implementation of MJ decomposition. Here we have used the function E b (B) which gives the total surface energy of the internal boundary segments B, of a given element patch e p , obtained thanks to the function Boundaries(e p ). Furthermore, we use the function GL(i, j, N) to return a list of size i with the j th set of consecutive (consecutiveness is measured in this context by following the angular coordinate in the polar coordinates system) boundary segments attached to Node N (i.e., for the case given in Figure 1, if z (P i ) > 3 then 3: for all number of connections of N i : j do 11:

18:
S min ← tuple (e pj , B j ) 19: if E min < E 0 then 20: Replace S 0 by S min 21: else if i < z (P i ) /2 then 22: i ← i + 1 23: goto 10: Algorithm 1 first searches between each pair of consecutive segments, the one that would reduce the boundary energy the most if it is separated from the MJ (just as depicted in Figure 1), and selects this configuration. Then, if this configuration reduces the initial GB energy given by the initial state, the initial configuration is replaced, and the algorithm continues to the next MJ. If not, instead of searching between pairs, the algorithm will re-iterate between consecutive triplets of lines if the order z of the MJ is sufficiently high (at least z = 6) and so on. Finally, if no configuration tested has lower energy than the initial configuration, the algorithm considers the MJ as stable and continues to the next. Note that the decompositions are made one at a time for a given call of Algorithm 1 over a given MJ. This means that a MJ of order z = 5 might be entirely decomposed in two increments and one with z = 6 in three.
This procedure accepts configurations with higher energy than the total minimal energy state (the configuration with the total minimum boundary energy), especially for MJs of high order (z > 6). However, in practice, such configurations have a very low probability of appearance in real microstructures.

Numerical Results
In this section, the TRM model will be tested in a heterogeneous context, meaning that the influence of the inclination angle ω over the value of γ will be ignored. Then, γ depends only on the disorientation angle θ and is equal for all segments defining the boundary between two given grains but different from all other boundaries. In such a context, the torque term τ is equal to 0 for all boundary segments, and the velocity of all nodes can be computed using Equation (3).
All tests performed in this section have been inspired by the ones presented in [8,9] in the same heterogeneous context. In [8], the classical FE-LS formulation of [24][25][26] has been reformulated with the primary objective of taking into account the gradients terms produced by a variation of γ, that were otherwise neglected in a homogeneous context (where γ is constant in Ω). Given that this formulation considers all variational terms relevant in this context, it will be named hereafter the heterogeneous FE-LS formulation. The numerical testing of this approach was divided into two parts: firstly, in [8], the numerical analysis is focused on the evolution of multiple junctions as a means to test the heterogeneous FE-LS formulation presented in the same publication. Secondly, in [9], the same heterogeneous FE-LS formulation was tested in the context of heterogeneous GG, using different formulations for the computation of the grain boundary energy γ, as a function of the disorientation angle. The approach used for the computation of the misorientation and disorientation angles can be found in Appendix A.
We reproduce these studies in the following with the TRM model:

Triple Junction Test Case
The first test corresponds to an academic triple junction test. Here, the motion of MJs is dictated by the GB energies of the interfaces meeting at the central node. Figure 3 (left) illustrates this aspect, where the term γ ij denotes the GB energy between grains G i and G j and φ i is the angle measured at the junction between the interfaces of grain G i and the other two grains. For this test, γ 13 = γ 23 , this will provoke a vertical movement of the junction for any value of γ 12 = γ 23 , until it arrives at its equilibrium position. As such, we will study the motion and the equilibrium of the junction in function of the ratio r = γ 23 /γ 12 . This test used dimensionless simulations, the value of the grain boundary energies γ 13 = γ 23 = 0.1 and the mobility term µ = 1 were held constant; moreover, for practical reasons (for r < 1, the MJ moves downwards, which when using Neuman type boundary conditions, induces a global movement of the interfaces in the same direction, and eventually leads to the contact of the junction with the base of the triangle. This behaviour is not wanted in this context), Dirichlet boundary conditions with v i = 0 ∀ N i ∈ ∂Ω were imposed, hence impeding the movement of the nodes at the intersection of the GBs and the edges of the triangular domain. Finally, the mesh size parameter h trm = 0.006 and the time step dt = 5 × 10 −5 will be used for all tests. These values were selected correspondingly to the limit of the stability region of the TRM model when using piece-wise polynomials (splines) as a means to obtain values of curvature, and normal [17]. The initial mesh is illustrated in Figure 3. While there is not an analytical formulation for the movement of the triple junction during its transient state in this context, triple junctions present stationary dihedral angles relying on the energies of the grain boundaries meeting at the junction [22]. In the absence of torque terms, i.e., when the energy of each interface is maintained constant, the dihedral angles φ 1 , φ 2 and φ 3 (see Figure 3) verify the Young's equilibrium, leading to the relation: Similarly to [8], we tested ratios in the range of r = [0.53, 10], and the obtained equilibrium angles were compared to the analytical equilibrium state obtained thanks to Equation (8). Figure 4 (left) illustrates the evolution of the ϕ 3 angle for different values of r obtained with the TRM model. These values are compared to the ones obtained in [8] (see Figure 4 (right)), where we have found that the TRM model evolves faster to its equilibrium state than the heterogeneous FE-LS method for values of r > 1.67. Furthermore, the TRM model can reproduce more accurately the analytical values of ϕ 3 for r < 1.0. Figure 5 also illustrates this aspect, where the final value of ϕ 3 is plotted against the grain boundary energy ratio r and compared to the analytical equilibrium value via an L2-Error plot.   Figure 6 illustrates the final interface states for both approaches at the end of the simulation. In [8], it was found that, while the equilibrium angles of φ 3 were accurately described for values of r > 2.5 near the junction, the behaviour of the interfaces was strongly affected by the boundary conditions applied to the FE resolution, inducing nonminimal energy configurations. This behaviour was not found nor expected with the TRM model as boundary conditions only affect the velocity of nodes belonging to the boundary, and as a result, the TRM model reduces in all cases (until the equilibrium) the total energy of the system. This result can be found in Figure 7 where the evolution of the normalised GB energy E Γ (each curve was scaled to start from a value equals to 1), has been plotted as a function of time.

2D GG with Heterogeneous GB Properties
Similarly to the triple junction test, in this section, we reproduce the same testing approach of [9] for heterogeneous FE-LS simulations of 2D-GG.
The first set of simulations measures the accuracy of the TRM model to reproduce results using different sets of mesh size and time step (h trm , dt). Results of these simulations are given in Appendix B. These simulations used a Read-Shockley (RS) type function [27] for the determination of the GB energy γ as a function of the disorientation angle θ: where γ max is the maximal grain boundary energy equals to 1.012 J m −2 , θ max corresponds to a threshold angle of 30 • and the mobility term µ has been held constant and equal to 0.1 mm 4 J −1 s −1 . These values are identical to the ones used in [9] for pure Nickel at 1400 K, where the authors have explained that contrary to the common usage of the RS function (using values for θ max in the range of [10,15] • ) a value of θ max = 30 • enables to numerically increase the system's heterogeneity. However, even with this choice, the heterogeneity level using a RS type function remains minimal. Indeed, only a narrow percent of the grain boundaries present a disorientation angle in the "variational" zone of the RS function (see Figure A1 (bottom-right) for the values with a disorientation angle θ < 30 • .) while the majority of the interfaces present a disorientation angle θ ≥ 30 • , and thus they acquire a value of γ = γ max . In [9], as a means to extend the representativity of the heterogeneous LS-FE formulation, multiple functions were used to compute the value of γ as a function of θ. This section will test the TRM model using two of the proposed functions. As such, the results presented here can be directly compared to those detailed in [9]. These functions correspond to the Read-Shockley+ and the Gaussian functions, which produce the most heterogeneous configurations. These functions are defined as follows: where γ max = 1.1 Jm −2 and θ thresh = 55 •

Gaussian
where γ q = 1.54 Jm −2 , θ µ = 40 • and θ σ = 10 • . These formulations were used along with the RS function and a homogeneous formulation (γ = 1.012 Jm −2 ) in the TRM model for the full-field modelling of annealing. These simulations were performed over four different initial Laguerre-Voronoi tessellations [28] based on an optimized sphere packing algorithm [29] and representative of the same statistical grain size distribution given in [9] with approximately 40,000 grains each. One example of initial tessellation is given in Figure 8. Hereafter, all results will contain data taken from the results of the four initial states and mean quantities will be averaged.  Figure 9 shows the evolution of this tessellation in time and for different γ functions, here, the colours are representative of the GB energy of each interface. These figures illustrate how the RS formulation is too "homogeneous", presenting just a few variations in the GB properties (even at the end of the simulation), while the RS+ and Gaussian are more heterogeneous. Additionally, in the RS+ and Gaussian cases, interfaces with a high GB energy seem to be eliminated during the early stages of the simulations, giving a higher predominance to low-energy GBs, which is not the case for the RS configuration. Another essential aspect observed in the final states of the RS+ and Gaussian cases is the appearance of stable high-order multiple junctions as predicted in Section 2.2. Figure 9. Examples of the evolution of the microstructure, from top to bottom as a function of time, and using the functions from right to left: RS, RS+, and Gaussian. A much higher heterogeneity is found in the cases using the Gaussian and RS+ functions. These figures illustrate the presence of stable high-order multiple junctions. Figure 10 illustrates the normalised GB disorientation distributions of the heterogeneous configurations for every hour of annealing. Results show how the RS maintains its shape near the Mackenzie plot, hence not giving almost any preference to low energy GBs. Contrarily, the Gaussian and RS+ cases tend to avoid the disorientation angles with high energy. The plot shows maximum values at disorientation angles with low energy (e.g., for the RS+ configuration, finds one maximum at θ > θ thresh = 55 • ). These results can also be observed in terms of the normalised grain boundary energy distributions given in Figure 11. In only one hour of annealing, the Gaussian and RS+ configurations tend to dissipate high energy GBs, giving a much higher predominance to low energy GBs and promoting their permanency (or their appearance) as time advances. In contrast, for the RS configuration, the changes in the distribution of energy remain negligible. The Gaussian configuration is a perfect example of how the TRM algorithm responds to grain boundary energy minimisation when opposed to a highly heterogeneous configuration.  Low energy GBs predominance may produce a deceleration of the evolution of the grain size in the domain. Figure 12 illustrates the grain size distribution of the different test cases showing how the RS configuration produces a grain size distribution with larger sizes while the RS+ and Gaussian configurations promote smaller grains. Figure 13 gives the evolution of the mean grain size, the total number of grains, the total GB energy, and the total grain boundary length. The minimisation of the GB energy is much higher for the most heterogeneous cases (RS+ and Gaussian), even though their number of grains and mean size appears to have a "slower" evolution than the RS and homogeneous cases. Furthermore, the responses of the homogeneous and the RS cases are very similar. Figure 14b gives the total GB length plotted against the number of grains of the simulation showing how the Gaussian and RS+ cases have a higher value than the RS and the homogeneous cases. This result can not be anticipated as one could have guessed the contrary, by seeing the evolution curves of the total GB energy given in Figure 13 (bottomleft) as a function of time and as illustrated in Figure 14a as a function of the number of grains. This result is a product of the preference of the higher heterogeneous cases for grain boundaries of low energy, but also by the more frequent apparition of high-order multiple junctions that decelerate the reduction of the total GB length.
. Evolution of different parameters as a function of the number of grains, for the 2D heterogeneous GG test case simulated with the TRM model: (a) total grain boundary energy E Γ and (b) total lenght of GBs. Table 1 gives the CPU-time of each simulation, showing how the computational cost of the TRM model in this context may be more related to the length of boundaries than to the number of grains (see Figure 14b). Additionally, the differences between the computational cost of the homogeneous and the heterogeneous cases are very high. This can be explained by the fact that for the homogeneous case, it is not necessary to compute the misorientation at GBs nor the lowest energy configuration in the event of a separation of MJs. These operations are very demanding as both rely on an iterative computation of the lowest rotation angle between two orientations in a set of 24 possible rotations, where all of them have to be tested. Table 1 also shows that the higher the heterogeneity of the case, the higher its computational cost. This behaviour can be anticipated by seeing the evolution of the number of grains and the total length of boundaries (Figure 14), as the more homogeneous cases reduce these quantities much faster. The results presented here are very similar to the ones obtained in [9] in the context of the heterogeneous FE-LS formulation. This suggests that both methodologies are valid to predict microstructural states in a full-field context, as even though the mechanisms behind their evolution are the same, they have been modelled using a completely different numerical scheme and still produce a very similar outcome. It mush be highlighted that simulations in [9] where performed with initial states with around 5000 grains while here we performed simulations eight times larger (for simulations using the same initial grains as in [9] see [30]). Moreover, the computational power needed to produce these results using the heterogeneous FE-LS may be much higher than the one needed by the TRM. The TRM model performed all sequential simulations in less than eight hours for the heterogeneous configurations and in less than one hour for the homogeneous ones, using an AMD Ryzen 73,600× processor.

Discussion, Conclusions and Perspectives
This article has provided the necessary implementation for modelling grain growth using heterogeneous grain boundary properties with the TRM model. These implementations consist of: (i.) a numerical framework on top of the TRM base code to measure neighbors' misorientation. The algorithm only takes these measurements at grain interfaces, namely, L-Nodes and PP-Connections. (ii.) A decomposition algorithm for high-order multiple junctions, which searches for the lowest energy configuration among all possible decompositions. These decompositions are obtained by the separation of pairs of interfaces from the MJ, storing for each, the total GB energy change ∆E Γ and applying the one with the lowest ∆E Γ only if it is negative (as for events with a minimum value of ∆E Γ > 0 the original configuration should remain stable). Finally, (iii.) a formulation for the computation of the velocity using anisotropic data was established using a discrete formulation inspired by the literature [12,15]. The new methodology implemented for the TRM model was tested in the context of heterogeneous grain boundary properties, using identical test cases like the ones presented in [8,9]. Results show that the TRM model can produce more accurate results regarding the equilibrium angles of triple junctions compared to the analytical values given by Young's equilibrium. Additionally, the TRM model ensures at all times low-energy stable configurations contrary to the heterogeneous LS-FE model presented in [8] which may produce stable configurations with non-minimal energy states. Furthermore, the TRM model was tested in a GG context using heterogeneous grain boundary properties in function of the disorientation angle. The initial configurations of all tests were statistically identical to the one presented in [9] with around 40,000 initial grains. Sensitivity analyses were performed, resulting in a tendency of the model to converge to a fixed solution when decreasing the set of parameters (h trm , dt), controlling the mesh size and the time step, respectively. Then, multiple formulations for the determination of the grain boundary energy γ as a function of the disorientation angle θ were used, namely the Read-Shockley (RS) [27], the modified Read-Shockley (RS+), and the Gaussian formulations used in [9]. Results showed a similar statistical behaviour to the results presented in [9] in a LS-FE context, hence validating both approaches at this scale.
Results also show that the CPU-time depends on to the total length of GBs. Additionally, the computational needs of the heterogeneous cases are higher than when using a homogeneous configuration. This result is strongly related to the computation of the disorientation angle which, even if it is only performed at the interfaces, it remains a brute force algorithm, which in [30] showed a poor performance, taking up to 60% of the total CPU-time in the heterogeneous configurations.
The proposed test cases use properties which are representative of pure Nickel at 1400 K (as in [9] ), however, the scope of this article remains purely academic, validating the approach through comparisons with well-known numerical models (i.e., the LS-FE model) and not with existing experimental data. Furthermore, although the LS-FE approach has been validated using experimental data in other contexts [31], it is not the case concerning the impact of anisotropic grain boundary properties. This is given by the fact that it is experimentally extremely difficult to measure local dynamics of grain boundaries and to correlate these data to their anisotropic properties (this is in fact an active field of research). Of course, It remains a perspective of the present work to test the TRM model in a fully anisotropic environment, where the influence of the inclination of the interface on the computation of γ is taken into account, producing variations of properties over curved GB and torque terms. Such a study will be presented in a forthcoming publication [20] on statistical behaviour of GB when using anisotropic GB properties that is more in accordance with the way that in situ experimental data at the polycrystalline scale are classically presented [32], paving the way to experimental-numeric correlations.
Finally, the TRM model (and its implementation for isotropic GB properties) can be extended to 3D by following the same set of rules defined for 2D modeling in [17]. These rules remain principally, the selective remeshing procedure that needs to carefully reconnect mesh entities (nodes, facets and elements) that belong to GBs. Then, a 3D surface approximation able to compute curvature and normal at GBs needs to be stablished. This constitutes a particularly difficult task since local gradients of curvature are needed to establish proper GBs' dynamics; however, the literature offers some workarounds [33].
Finally, the discrete methods presented to compute DRX dynamics [19] and anisotropic behavior need to be developed, the latter, being addressed in the past in [22,34].

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available as the data forms part of an ongoing study. the one with the lowest θ. In this work, we will consider only cubic-type crystals, hence 24 symmetric representations must be iterated.
TRM model performs these operations during two stages of the algorithm: First, a misorientation computation is held before the computation of the nodes' velocities v i , as all boundary properties must be defined at this stage. The computation is performed once per Line, which attributes all misorientations and disorientation angles for all L-Nodes and segments of the Lines entities. Note that a Line can only compute one misorientation, hence, it is unnecessary to compute it for every node lying on the Line. Then, some edges still need to define their orientation: the ones describing a PP-Connection, namely, the edges defining a connection between two Point entities (see [17] for more information about the data structure used by the TRM model). This computation is necessary, as, even though the notion of grain boundary energy does not hold at MJs in the same way as for normal boundaries, the GB properties of all interfaces attached to the MJ are needed.
Secondly, the TRM model performs a misorientation computation during the decomposition of MJs for all possible new interfaces (see line 13-15: of Algorithm 1). This could be a very demanding procedure as, for instance, each possible decomposition seeks the minimal disorientation angle among all 24 equivalent symmetries defined for the crystal. We study the relative cost of the misorientation computation at the end of Section 3.

Appendix B. Sensitivity Analysis on GG Simulations: Mesh Size and Time Step
This section uses a squared RVE domain of 1.5 mm of side length to model annealing. Figure A1 (top-left) illustrates the initial state of the polycrystal used in this sensitivity analysis. This polycrystal contains exactly 5089 initial grains and its grain size distribution (pondered by surface) is given in Figure A1 (top-right). Additionally, in all cases, the mobility term has been held constant and equal to 0.1 mm 4 J −1 s −1 . Figure A1 (bottom-left) shows the initial microstructure colored following the orientation magnitude e = ϕ 2 1 + Φ 2 + ϕ 2 2 of each grain. Disorientation angles (θ) have been computed for each Line, L-Node and PP-Connection (i.e., for all nodes and segments belonging to the GBs) of the interface using the methodology presented in Appendix A. Figure A1 (bottom-right) gives the disorientation angle distribution of the initial microstructure, which shows a good agreement with the Mackenzie plot. Figure A2 (top) gives the evolution of various parameters for the first 3 h of simulated time using a constant mesh size of h trm = 0.003 mm and for various time steps dt. The mean grain size, the number of grains, and the total GB energy have been plotted, showing a tendency to converge to a fixed solution when the time step decreases. Figure A2 (bottom) gives the L2-difference of each iteration taking as a reference the curve using dt = 5 s, confirming these results. Similarly, the simulations were repeated using a constant time step dt = 50 s and for various mesh sizes. Similarly, as before, decreasing the mesh size produces a tendency to converge to a fixed evolution (see Figure A3), reducing the L2difference to the reference curve (here the one using h trm = 0.002 mm) with every iteration. This study shows that one can expect good accuracy when using a set of parameters (h trm , dt) in the surroundings of (0.003 mm, 50 s). These values will be used in all other polycrystal simulations.