The Role of Surface Tension and Mobility Model in Simulations of Grain Growth
Tiago Salvador Selim Esedoglu
cond-mat.mtrl-sci physics.comp-ph





Department of Mathematics, Michigan University, 530 Church St. Ann Arbor, MI 48105 (saldanha@umich.edu).



Department of Mathematics, Michigan University, 530 Church St. Ann Arbor, MI 48105 (esedoglu@umich.edu).





We explore the effects of surface tension and mobility models in simulations of grain growth using threshold dynamics algorithms that allow performing large scale simulations, while naturally capturing the Herring angle condition at junctions and automatically handling topological transitions.
The results indicate that in two dimensions, the different surface tension / mobility models considered do not play a significant role in the stationary grain size distribution.
However, in three dimensions, there is a substantial difference between the distributions obtained from the same three models, depending on whether the reduced mobilities are isotropic or anisotropic.
Additional results show that in three dimensions, the misorientation distribution function of a grain network with random orientation texture returns to the close vicinity of the Mackenzie distribution even if started very far from it.








The Role of Surface Tension and Mobility Model
 in Simulations of Grain Growth
    Selim Esedoḡlu
    December 30, 2023
==============================================================================












§ INTRODUCTION




The shapes and sizes of grains constituting polycrystalline materials, such as most metals and ceramics, influence important physical attributes of these materials, for example their conductivity and yield strength.
Therefore, there has been much interest in modeling and simulation of grain boundary motion during common manufacturing processes, such as annealing.
In particular, the evolution of grain size distribution (GSD), the grain boundary character distribution (GBCD), as well as a simpler version of the latter, the misorientation distribution function (MDF), have been a central focus of research in materials science.

In this paper, we report simulation results for grain boundary motion at the mesoscale, using a recent, very concise and flexible numerical algorithm developed for the continuum model of grain boundaries that was introduced by Mullins <cit.> and earlier authors.
According to this well known model, interfaces in the grain network evolve via steepest descent for an energy that penalizes surface areas of the interfaces, often weighted by a factor, known as the surface tension, that depends on the misorientation across each interface.
The resulting dynamics is described by a system of partial differential equations (PDEs), coupled to each other along junctions by what are known as Herring angle conditions <cit.>, and allows specifying in addition a mobility factor for each interface.
The results we present are the first large scale simulations using the new algorithm proposed in <cit.> for this PDE system.

We investigate several questions that have repeatedly come up in previous literature concerning certain statistics of the grain network.
In particular, we study, in both two and three dimensions:


  * The effect of surface tension and mobility model used on statistics such as the GSD. 

  * The asymptotic behavior of the MDF and its dependence on initial configuration.

Specifically, we simulate in both 2D and 3D Mullins' mesoscale model with 


  * Read-Shockley surface tensions <cit.> (and its extension <cit.> to 3D anisotropy) along with isotropic mobilities, 

  * Read-Shockley surface tensions with “reciprocal” mobilities resulting in isotropic reduced mobilities but asymmetric junction angles, and

  * Isotropic surface tensions and mobilities, resulting in isotropic reduced mobilities as well as symmetric (all 120^∘) junction angles.


We compare to experimental data whenever it is available in existing literature.
Some of our results and conclusions differ substantially from earlier studies.
Highlights of our findings are:


  * In 2D simulations, the stationary GSD that is reached appears to not depend on the surface tension and mobility models used: it is essentially the same in all three surface tension / mobility combinations (i), (ii), (iii) used in our simulations.
This is in agreement with a variety of previous studies, e.g. <cit.>.

  * However, in 3D, the stationary GSD reached is substantially different when Read-Shockley surface tensions along with isotropic mobilities (i) are used, compared to (ii) and (iii) that result in isotropic reduced mobilities.
This finding appears to be distinct from earlier studies, e.g. <cit.>.

  * Moreover, the different GSD reached in 3D by models (i) and (ii) that differ only in their choice of mobilities constitutes a contrary point to previous studies that suggest mobilities do not play a significant role in the long time behavior of the GSD.

  * The GSD reached by (ii) and (iii) are almost identical in both 2D and 3D, even though junction angles in these simulations are very different, suggesting that perhaps the isotropy vs. anisotropy of reduced mobilities is what is significant.

  * The GSD reached in our 3D simulations with surface tension / mobility combination (i) is not less consistent with experimental measurements <cit.>, which have substantial variation themselves.

  * Statistics of grain shapes, such as that of isoperimetric ratios (eccentricities), also display considerable dependence on the choice of surface tension / mobility model.

  * Models (i) & (ii) allow us to study the asymptotic behavior of the MDF.
Confirming earlier studies, the two drastically different choices of mobility in (i) & (ii) lead to almost identical behavior, both in 2D and 3D.
Distinct from  earlier studies, we exhibit 3D initial data with MDFs perturbed substantially away from the Mackenzie distribution and observe the MDF return to a very close vicinity of the Mackenzie distribution, using either (i) or (ii) as the surface tension / mobility model.


The code for the simulations is publicly available and can be found at <https://github.com/tiagosalvador/statisticsGBM>.






§ PREVIOUS WORK





The exising literature on the statistics of grain shapes and sizes during annealing is absolutely vast, both from the experimental and the modeling and simulation side.
Our discussion will therefore be necessarily limited to a few representative studies closest to the questions addressed in the present paper.

There are a plethora of numerical algorithms developed for the simulation of Mullins' mesoscale model of grain boundary motion.
They include the Monte Carlo Potts, front tracking, phase field, level set, and threshold dynamics methods.
In principle, all these methods should converge to the same continuum model described by the same system of PDEs, provided that various parameters appearing in them are scaled appropriately, at least until topological changes start to take place (even then, one hopes there would be agreement in the vast majority of situations).
Nevertheless, there isn't always good agreement even among simulations.

Very large, well resolved simulations of normal grain growth were carried out by Mason et. al. <cit.> using the front tracking method in both 2D and 3D.
The authors catalogue a large number of grain size and shape statistics,  for the two dimensional model, for the three dimensional model, and for two dimensional slices of the three dimensional model.
However, this work is restricted to the equal surface tension, equal mobility, i.e. the fully isotropic, version of Mullins' model.
That in particular means the exclusion of such statistics as the GBCD or the MDF which are predicated on variation of surface tension between different interfaces in the network.
Furthermore, left unadressed is the question of how sensitively these shape and size statistics depend on the surface tension and mobility model used in conjunction with Mullins' continuum model.

Upmanyu et. al. <cit.> study the dependence of grain boundary statistics on the surface tension and mobility models used via Monte Carlo Potts and phase field simulations in 2D.
The GSDs they report, including when the surface tension, the mobility, and both are anisotropic, show little deviation from that of the isotropic model; they all follow a more or less log-normal distribution.
Their investigation of the MDF in the same vein reveals strong dependence on the surface tension model used, but not so much on mobilities:
The network is driven towards a prevalence of interfaces with low misorientations and therefore low surface tensions.
Other statistics, such as the number of neighbors (or faces), are also observed to have a much more pronounced dependence on surface tensions than mobilities.
Similar MC Potts simulations in <cit.> also show little difference between the GSDs of isotropic vs. Read-Shockley surface tensions.

Gruber et. al. <cit.> carry out MC Potts simulations in 3D, using isotropic as well as Read-Shockley (anisotropic) surface tension models.
Among the grain properties they document are the GSD and the MDF.
They observe no difference between the stationary GSDs reached using isotropic vs. anisotropic surface tension models.
As regards the MDF, this study suggests surface tension anisotropy has a measurable effect on the steady state MDF reached in simulations that start from an initial data with random orientation texture, regardless of the dimension.
When the initial orientation texture is not random, no steady state for the MDF is observed in the case of anisotropic surface tensions: instead, a proliferation of low angle grain boundaries results in a gradual but persistent concentration of the MDF at the origin.

Kim et. al. <cit.> carry out phase field simulations in 2D and 3D using isotropic surface tensions and mobilities. In 3D, the stationary GSD they reach (given in terms of the reduced equivalent radii) follows very closely the Hillert distribution.
On the other hand, the Monte-Carlo Potts simulations of Zoellner et. al. <cit.> in 3D, using isotropic surface tensions and mobilities, and starting from an initial data with random texture, results in a stationary GSD that is quite far from the Hillert distribution.
The more recent very large scale simulations of Miyoshi et. al. <cit.> in 3D using the phase field method also yielded a GSD similar to that of <cit.> that deviates substantially from the Hillert distribution.

The recent 2D simulations reported in <cit.> using a phase field method treat surface tension / mobility models (i) and (iii), as well as a third model that combines Read-Shockley surface tensions with mobilities that have a sigmoidal dependence on the misorientation angle as proposed by Humphreys <cit.>.
For this third surface tension / mobility model, the stationary GSD they obtain is significantly different from those with models (i) & (iii).
The stationary GSDs of models (i) & (iii) are reported to be very similar.











§ OUR MODEL AND ALGORITHM




Our starting point is the curvature driven grain boundary motion model that goes back at least to Mullins <cit.>.
According to this model, grain boundaries evolve via steepest descent for the energy

E = ∑_i<j σ_ijArea(Γ_ij) 
where Γ_i,j denotes the boundary between two adjacent grains Σ_i and Σ_j. 
In particular we will neglect the normal dependence of the energy density, but allow each grain boundary Γ_i,j to have a distinct surface tension σ_i,j.
The very special case σ_i,j = 1 for all i≠ j is often referred to as the isotropic surface tension model.
We will also be interested in the surface tension model due to Read & Shockley <cit.> and its extension <cit.> to 3D crystallography:

σ_ij = 
θ_ij/θ_*(1-log(θ_ij/θ_*))	    if  θ_ij ∈[0,θ_*],

1	    if  θ_ij ≥θ_*.

where θ_i,j is the misorientation angle between Σ_i and Σ_j and θ_* is the  Brandon angle <cit.> , often taken to be between 15^∘ and 30^∘.
Following the convention in materials science literature, we will refer to (<ref>) as an anisotropic surface tension model, even though there is still no dependence on the normal to the interface.

In two dimensions, each grain Σ_i is modeled as a square lattice and is assigned an orientation angle θ_i of a clockwise rotation about the origin that maps it back to the standard two-dimensional lattice ^2. Thus, the misorientation between two grains Σ_i and Σ_j is given by

θ_i,j = min_k∈ θ_i-θ_j+kπ/2.
In three dimensions, the orientation of a grain with cubic lattice can be described (nonuniquely) by a matrix g ∈ SO(3) that corresponds to the rotation required to obtain the lattice of the grain from the standard integer lattice ^3. The misorientation angle between two grains Σ_i and Σ_j is defined as

θ_i,j = min_r∈ arccos(1/2(( rg_ig_j^T)-1)),
where  denotes the octahedral group (of symmetries of the cube in the three dimensions).

According to Mullins, the dynamics associated with energy (<ref>) is given by L^2 gradient descent for the interfaces, leading to the normal speed

v_i,j = μ_i,j σ_i,j κ_i,j
for interface  Γ_i,j, where κ_i,j denotes the mean curvature of Γ_i,j and μ_i,j>0 is a mobility factor. In addition, a condition known as the Herring angle condition <cit.> holds along triple junctions. At a junction formed by the meeting of the three distinct grains Σ_i, Σ_j and Σ_k, this condition reads

σ_i,j n_i,j + σ_j,k n_j,k + σ_k,i n_k,i = 0.
As a consequence the angles formed by normals n_i,j, n_j,k and n_k,i to the three interfaces Γ_i,j, Γ_j,k and Γ_k,i along the triple junction are determined by their associated surface tensions; this relation is also known as Young's law.

The algorithms used in this study are obtained from a non-local approximation to energy (<ref>)

1/√(δt) ∑_i<j σ_i,j ∫_Σ_i G_√(δt) ∗1_Σ_j where G_t denotes the Gaussian kernel

    G_√(t)(x) = 1/(4π t)^d/2exp(-x^2/4t)

and 1_Σ(x) denotes the characteristic function for a set Σ:

    1_Σ(x) = 
    1    if  x∈Σ,
    
    0    otherwise.

The width δ t of the Gaussian kernel appearing in (<ref>) ends up playing the role of the time step size for our scheme, described below, that approximates gradient decent of (<ref>) in L^2 sense, as described by <cit.>. Energy (<ref>) has been proposed in <cit.> and has been shown to converge in a very precise sense (namely, that of Gamma-convergence) to energy (<ref>) in the limit δ t→ 0^+. Intuitively, we have

    1/√(δ t)∫_Σ_i G_√(δ t)∗1_Σ_j≈1/√(π)Area(Γ_i,j)

since the function 1/√(δ t)1_Σ_i G_√(δ t)∗1_Σ_j approximates a delta function concentrating near Γ_i,j as δ t → 0^+. The reason for our interest in this specific - perhaps unusual - approximation of Mullins' energy is that it generates exceptionally simple and efficient algorithms for simulating gradient descent dynamics associated with (<ref>). Indeed, it has been shown in <cit.> to lead, in a systematic way, to the correct multiphase, arbitrary surface tension analogue of a very fast algorithm known as threshold dynamics that was originally proposed in <cit.> for networks with all equal surface tensions (i.e., θ_i,j = 1 for all i and j) and mobilities.

We start by recalling Algorithm <ref>, the simplest version of the generalization of threshold dynamics to arbitrary surface tensions given in <cit.>. Its benefits include its unconditional stability (time step size δ t can be chosen arbitrarily large, constrained only by accuracy considerations), seamless handling of topological changes in any dimension and very low per time step cost (Nlog N) on a uniform grid of N points, using the FFT to compute the convolutions. Yet, despite its extreme simplicity, it handles automatically and correctly all the essential features of the dynamics, including Herring angle conditions (<ref>) along junctions, and countless types of topological changes that may occur during the evolution. Note that Algorithm <ref> can be seen as a type of level set method, where the function

    -ψ_i^n(x) + min_j≠ iψ_j^n(x)

plays the role of the level set function depicting the shape of the i-th grain at the (n+1)-th time step. It turns out that this very simple algorithm results in a misorientation dependent mobility μ_i,j, associated with the boundary Γ_i,j between grains i and j that is given by

    μ_i,j = 1/σ_i,j

leading to the normal speed

    v_i,j = μ_i,jσ_i,jκ_i,j(x) = κ_i,j(x)

at any point x∈Γ_i,j.
In other words, using the Gaussian as the convolution kernel, Algorithm <ref> naturally prefers  equal (isotropic) reduced mobilities: Each interface in the network moves with the same (equal) multiple of its (mean) curvature, and the only difference from the isotropic model are the junction angles.

(in <cit.>)Given the initial grain shapes Σ_1^0, …, Σ_N^0 and a time step size δ t, obtain the grain shapes Σ_1^n+1, …, Σ_N^n+1 at the (n+1)-th time from the grain shapes Σ_1^n, …, Σ_N^n at the end of the n-th times as follows:

	
  * Compute the convolutions:

    ϕ_i^n = G_√(δ t)∗1_Σ_i^k.

	
  * Form the comparison functions:

    ψ^n_i = ∑_j ≠ iσ_i,jϕ_j^n.

	
  * Update the grain shapes:

    Σ_i^n+1 = {x: ψ_i^n(x) < min_j≠ iψ_j^n(x)}.




A more general algorithm was in fact given in <cit.> that allows arbitrary choice of surface tensions and mobilities, but is considerably more complicated and deviates from the simplicity of Algorithm <ref>.
In this paper, we use for the first time on large scale simulations a more recent version, namely Algorithm <ref>, that was proposed in <cit.>. 
It is simpler and closer in spirit to Algorithm <ref> while allowing for both arbitrary surface tensions and arbitrary mobilities.
This is achieved by merely replacing the convolution kernel used in Algorithm <ref> by the positive sum of two distinct Gaussians, thus preserving the efficiency and simplicity of the original version.
As shown in <cit.>, in the special but important and very common case of equal mobilities (i.e., μ_i,j = 1 for all i≠ j), a careful choice of the widths of the two Gaussians, denoted by the parameters α and β in what follows, guarantees the retention of  the major benefits of the original Algorithm <ref>, namely unconditional stability and convergence to the correct Mullins' energy (<ref>), while allowing its extension to e.g. Read-Shockley surface tensions (<ref>).

Given the initial grain shapes Σ_1^0, …, Σ_N^0 and a time step size δ t, obtain the grain shapes Σ_1^n+1, …, Σ_N^n+1 at the (n+1)-th time from the grain shapes Σ_1^n, …, Σ_N^n at the end of the n-th times as follows:

	
  * Compute the convolutions:

    ϕ^n_1,i = G_√(αδ t) * 1_Σ_i^nϕ^n_2,i = G_√(βδ t) * 1_Σ_i^n.

	
  * Form the comparison functions:

    ψ^n_i = ∑_j ≠ i a_i,jϕ_1,j^n + b_i,jϕ_2,j^n

where a_i,j and b_i,j are given by

    a_i,j = √(π)√(α)/α-β(σ_i,j-βμ_i,j^-1)   and   b_i,j = √(π)√(β)/α-β(-σ_i,j+αμ_i,j^-1).

	
  * Update the grain shapes:

    Σ_i^n+1 = {x: ψ_i^n(x) < min_j≠ iψ_j^n(x)}.









§ LARGE SCALE SIMULATIONS




In this section, we use Algorithms <ref> and <ref> to study the evolution of large networks of grains, in both 2D and 3D.
The 2D simulations are initialized with approximately 100,000 well-resolved grains and are repeated 3 times, while the 3D simulations are initialized with 10,000 grains and repeated 10 times.
Recall from the Introduction the following surface tension and mobility models that will be explored in both cases:


  * Read-Shockley surface tensions with constant (isotropic) mobilities: σ_i,j given by (<ref>) and μ_i,j = 1.

  * Read-Shockley surface tensions and reciprocal mobilities: σ_i,j given by (<ref>) and μ_i,j = σ_i,j^-1.
This results in equal (isotropic) reduced mobilities.

  * Equal (isotropic) surface tensions and mobilities: σ_i,j = μ_i,j = 1.
This also results in equal (isotropic) reduced mobilities.

Whenever Read-Shockley surface tensions are used, we take the Brandon angle θ_* in (<ref>) to be θ_* = 30^∘ to agree with simulations performed in Refs. <cit.> and to lie within the experimentally observed range <cit.>.

In all simulations, the initial grain configuration is taken to be the Voronoi diagram of a set of points chosen uniformly at random from a periodic box.




 §.§ Two spatial dimensions



In these 2D simulations, the initial grain configuration of around 100,000 well-resolved grains is assigned a random fiber texture: All grain orientations are obtained from a reference configuration by a rotation about the axis normal to the 2D domain, and the angle of the rotation is taken from [0,2π) uniformly at random.
We present results at two distinct times for each model: model (i) at t_i = 1.907 × 10^-5 and t_f = 8.345 × 10^-5, model (ii) at t_i =  1.550 × 10^-5 and t_f = 6.259 × 10^-5 and model (iii) at t_i =  1.788 × 10^-5 and t_f = 7.033 × 10^-5. The different stopping times ensure that at the intermediate time t_i, around 30% of all grains remain in each one of our runs, while at the final time t_f around 10% of all grains remain.
The simulations are performed on a 4096 × 4096 grid discretizing [0,1]^2.


Algorithm <ref> requires choosing the parameters α and β, the widths of the two Gaussian kernels. 
They need to satisfy

    α > max_i≠ jσ_ij  and  β < min_i≠ jσ_ij

to guarantee the properties described in <ref>.
However, from a practical point of view,  we also need the quotient α/β to remain moderate for accuracy considerations so that the Gaussians have comparable widths and thus may be sampled well without requiring unduly high spatial resolution.
Since both α and β depend on the surface tensions which in turn depend on the misorientation angle, this can be accomplished by imposing a minimum misorientation angle between any two grains to be 1^∘. 
This reduces the effective size of the surface tension matrix, while keeping the same number of distinct grains, and ultimately results in a small quotient α/β. 

We achieve this in the following way. Once the grain orientations are assigned, we partition the domain [0,2π) into equally spaced bins of 1^∘ length. Then all grains with orientations in the same bin are reassigned a new grain orientation, the midpoint of the bin. This is equivalent to assigning grain orientations from the uniform distribution on {0.5^∘, 1.5^∘,…,359.5^∘}. 





Figure <ref> shows the time evolution of the GSD for the three different models, while Figure <ref> shows the corresponding distributions of reduced effective radii.
For models (i) & (ii) with Read-Shockley surface tensions, the distributions appear stationary with no significant difference between the two.
The distributions for the isotropic model (iii) are also very similar, except at final time t_f, we see some evidence of the two distinct peaks observed previously in simulations of <cit.> with the same isotropic model, but using the very different algorithm of front tracking.
This two peak distribution has not been observed in other recent 2D simulations of the isotropic (constant) surface tension and mobility model (i) in e.g. <cit.>, whose results look closer to the ones we obtained with Read-Shockley surface tensions.
On the other hand, the recent study <cit.>, which includes simulations with models (i) & (iii), also has some evidence of the double peak with model (iii); our results are in close agreement with theirs for both of these models.





The similarity in the GSDs of the three surface tension / mobility models considered here notwithstanding, it is natural to expect some differences when we move on to statistics of grain shapes.
After all, it is well known that the microstructure of anisotropic models can be visually quite different than that of the isotropic model (see Figures <ref> and <ref>).

Indeed, the similarities finally cease when we look at the isoperimetric ratio, also called circularity, which is given by 4π A/P^2 where A and P denote the area and perimeter of the grain, respectively.
Figure <ref> shows that all three surface tension / mobility models lead to peak distributions but with a higher amplitude for the isotropic model. From time t_i to time t_f, while 20,000 disappear, the isoperimetric ratio retains its peaked distribution but with the peak moving towards 1 which indicates that the grains become more rounded.


























Finally, Figure <ref> shows the MDF for the two unequal surface tension models: Read-Shockley with equal and reciprocal mobilities.
Due to its random fiber texture, the MDF of the initial data is uniform.
As time evolves, it concentrates at small misorientations corresponding to small surface tensions.
This is in agreement with previous studies, e.g. 
<cit.>.
There is a slight but noticeable difference between the MDFs of the two models, which may in any case keep concentrating at the origin with further evolution: Model (i) has a more concentrated MDF compared to (ii), perhaps because low angle boundaries move slower in model (i) and hence persist longer.








 §.§ Three spatial dimensions




For the 3D simulations, we initialize the microstructure with approximately 10,000 well-resolved grains, chosen as the approximate Voronoi diagram to points placed uniformly at random in a periodic box that is discretized into a 256 × 256 × 256 uniform grid.
This results in approximately 12 grid points across each grain in every direction at time t=0, which is roughly the same resolution as in our 2D simulations.
Each grain is assigned a fully random orientations from SO(3).
We present results at two distinct times for each model: model (i) at t_i = 4.578 × 10^-3 and t_f = 7.629 × 10^-3, model (ii) and (iii) at t_i =  1.831 × 10^-3 and t_f = 4.578 × 10^-3. 
Just like in the 2D simulations, the different stopping times ensure that at the intermediate time t_i, around 30% of all grains remain in each one of our runs, while at the final time t_f around 10% of all grains remain.

As in our 2D simulations, we impose a minimum misorientation angle, this time of 5^∘, between the grains once again in order to keep the quotient α/β that appears in our algorithm at a moderate value.
Given that the minimum misorientation angle imposed here is larger than usual, a natural question arises: if it were taken to be smaller, would the different statistics change significantly? In order to investigate this, we run our simulations for model (ii) without doing any merging and observed that the statistics remained essentially the same.






Unlike in our 2D simulations, we see substantial differences in the eventual GSD reached by surface tension / mobility model (i) vs. models (ii) & (iii).
Figures <ref> and Figure <ref> show the time evolution of the GSD in terms of the reduced volume and corresponding effective radii distributions, respectively.
We note that while the effective radii distributions of (ii) & (iii) remain very close to each other throughout the evolution (and agree with previous 3D simulations of the same model via different numerical methods, e.g. <cit.>), that of (i) is more spread out, and its peak is shifted in the direction of smaller grains – a difference not seen in the 2D simulations with the same algorithms, and at very similar resolutions.
Figure <ref> suggests the different eventual GSD observed with surface tension / mobility model (i) may not be inconsistent with experimental data:
Two experimental distributions available in existing literature, namely <cit.> and <cit.>, display substantial deviation from the GSD seen in simulations of the isotropic model (iii) in our and previous simulations such as <cit.>, in the same direction as our simulations with model (i), namely in the direction of smaller grains and more spread out distribution.







Figure <ref> shows the distribution of the isoperimetric ratio, namely 36π V^2/S^3, where V denotes the volume and S the surface area of the grain.
Remarkably, we see less difference between surface tension / mobility models (ii) and (iii) in 3D than we did in 2D simulations: their stationary isoperimetric ratio distributions are almost identical.
The distribution for model (i), however, is broader, indicating a greater presence of eccentric grains compared to the isotropic model (iii), as in our 2D results show in Figure <ref>.














We now turn to statistics of grain shapes and sizes in 2D slices of 3D data, often called 3DX statistics in the literature.
To this end, we took a series of cross sections parallel to the faces of the cubic computational domain [0,1]^3.
They were spaced roughly five average grain diameters apart to reduce correlations between neighboring cross sections. Figure <ref> shows the distribution of GSD in terms of the reduced effective radii distributions.
The results are similar to the fully 3D grains: the surface tension / mobility models (ii) and (iii) are close to each other, while the distribution for model (i) is more spread out and its peak is achieved at smaller grains. 
In Figure <ref> we compare the distribution of the isoperimetric ratio, namely 4π A/P^2, where A denotes the area and P the perimeter of the grain in the 2D slice. All surface tension / mobility models present almost identical peak distributions, with models (ii) and (iii) attaining a higher peak.





Finally, we report asymptotic behavior of the MDF in our 3D simulations with surface tension / mobility models (i) & (ii).
When the initial data is generated by assigning orientations at random sampled from the uniform distribution on SO(3) to grains from a Voronoi construction, it is well known that the initial MDF is given by the Mackenzie distribution <cit.>.
With this initial data, both models (i) & (ii) behave very similarly: Their MDFs stay in close vicinity of the Mackenzie distribution; see the first row of plots in Figure <ref>.
This observed behavior is in close agreement with previous simulation results, e.g. <cit.>.
The very slight deviation from the Mackenzie distribution (observed in these previous studies as well) appears to be present mainly in model (ii), and not so much in (i).

We also explore the asymptotic behavior of the MDF when the initial configuration deviates substantially from the Mackenzie distribution.
To that end, we devised and implemented an algorithm to modify the randomly assigned initial orientations so that the initial MDF is perturbed substantially away from the Mackenzie distribution.
We leave the shapes of the initial grains unchanged, as the random Voronoi diagram.
To be specific, consider the function




    MDF_δ (x) = ∑_i<jArea(Γ_ij) b_δ(x-θ_ij)

 for δ>0 with
 
    b_δ(x) = 1/δ√(2π)e^-x^2/2δ^2.

 Without loss of generality assume that
 
    ∑_i<jArea(Γ_ij) = 1.

Then MDF_δ is a regularized approximation of the MDF of the grain network.
Let T(x) denote the desired initial MDF.
Based on (<ref>), MDF_δ can be treated as a function of g_1,…,g_N ∈ SO(3), the orientations of each grain, by expressing the misorientation angles θ_ij in terms of these.
We apply a steepest descent procedure on the L^2 norm of the difference between MDF_δ and T with respect to the orientations g_i of the grains; this updates the initial orientations g_j until the initial MDF closely matches the desired distribution T(x).
As this procedure has no preference for a particular orientation, we expect that the initial orientation distribution remains uniform.
In other words, out approach allows perturbing the MDF without introducing texture.

The second and third rows of Figure <ref> show the evolution of the MDF when it is initially perturbed away from the Mackenzie distribution by this procedure, by two different sized (and quite large) perturbations.
In each case, the MDF evolves back to a very close vicinity of the Mackenzie distribution, which is quite different from earlier results in <cit.>, where no steady state MDF is observed in simulations with initial data whose MDFs differ from the Mackenzie distribution.
We expect the difference is due to the different ways of generating non-Mackenzie initial data: The method of <cit.> introduces texture (non-uniform orientation distribution), while ours does not. This is illustrated in Figure <ref> where we plot the orientation angle distribution. Here, by orientation angle of a grain we mean its misorientation angle when measured against the sample reference grain, the one that is completely aligned with the axes. In the presence of a random texture the orientation angle distribution will match the Mackenzie distribution.

 
 
 


 





§ CONCLUSIONS





We investigated, via large scale simulations, the role that surface tension and mobility models play in continuum descriptions of grain growth.
This was accomplished by applying in the context of very large numbers of grains the new, extremely simple and streamlined threshold dynamics algorithm recently introduced in <cit.> together with the algorithm previously introduced in <cit.> to perform 2D and 3D simulations.
Grain statistics, such as the GSD and the MDF, for three distinct choices of surface tension / mobility models were compared to each other and to available experimental data in existing literature.
The 2D simulations did not reveal a strong dependence of the GSD on the surface tension / mobility model used.
However, in 3D, a clear distinction emerged between the stationary GSDs of the model with anisotropic reduced mobilities and the other two models considered that had isotropic reduced mobilities.
The asymptotic behavior of the MDF was also studied, both in 2D and 3D, using 2D and 3D crystallography, respectively.
In 3D, it was observed that as long as the orientation texture of the initial data is random, even an initial MDF that is very far from the Mackenzie distribution eventually enters a very close vicinity of it.





§ ACKNOWLEDGEMENTS



The authors gratefully acknowledge support from the NSF grant DMS-1719727.


amsalpha




