## Abstract

Diffusion-generated motion is used to perform a very large-scale simulation of normal grain growth in three dimensions with high accuracy. The method is based on the diffusion of signed distance functions and shares similarities with level-set methods. The Herring-angle condition at junctions and topological transitions are naturally captured with this formulation. This approach offers significant advantages over existing numerical methods and allows for accurate computations on scales not previously possible. A fully resolved simulation of normal grain growth, initially containing over 130 000 grains in three dimensions, is presented and analysed. It is shown that the average grain radius grows as the square root of time and the grain-size distribution is self-similar. Good agreement with other theoretical predictions, experimental results and simulation results via other techniques is also demonstrated.

## 1. Introduction

Grain growth is an important process by which the microstructure of a polycrystalline material (including most metals and ceramics) evolves during manufacturing processes. Statistical measures of the resultant microstructure affect important macroscale properties of the material, such as its conductivity and brittleness. Manufacturing processes must typically be tuned to provide for an optimal blend of desired material properties; however, performing such tuning experimentally is costly and time consuming. As a result, simulations of grain growth have been attempted using a variety of numerical techniques. Several common techniques are described in §3.

Grain growth occurs when polycrystalline materials are annealed. The well-known model for grain growth (Harker & Parker 1945; Beck 1952; Mullins 1956) gives the normal velocity (outward from phase *Σ*_{k}) of the interface *Γ*_{kℓ} by
1.1
Here, *μ* denotes the boundary mobility, *γ*_{kℓ} the grain-boundary energy per unit area for the interface *Γ*_{kℓ} and *κ*_{kℓ} the mean curvature of the boundary separating two grains. We use the convention that if phase *Σ*_{k} is a spherical grain of radius *r* surrounded by phase *Σ*_{ℓ}, *κ*_{kℓ}=−2/*r*.

We specialize to the case where all surface tensions are constant and equal: *γ*=*γ*_{kℓ}, often called ‘isotropic grain growth’. We aim to extend to the more general case given by equation (1.1) in future work. The theory for this extension is nearing completion in two dimensions in the absence of topological events. Additional study of these events (including the division of four junctions following grain disappearance) is still required. Here, we non-dimensionalize the normal velocity using the mean initial grain radius 〈*r*_{0}〉, derived from the mean initial grain volume 〈*V* _{0}〉 by 〈*r*_{0}〉=(3〈*V* _{0}〉/(4*π*))^{1/3}. We define the non-dimensionalized curvature as . We further non-dimensionalize the velocity as , with velocity *V* =〈*r*_{0}〉/*T* and time *T*=〈*r*_{0}〉^{2}/(*μγ*), so that
1.2

As shown by Reitich & Soner (1996) and Zhao *et al.* (1996), this non-dimensionalized normal speed arises as a gradient descent for the non-dimensional energy
1.3
We note that the time scale *T* is chosen so that *t*^{★}=1/4 is the time required for an isolated grain of radius 〈*r*_{0}〉 to disappear under pure curvature motion in the non-dimensional system. Hereafter, we drop the ★ notation and refer solely to the non-dimensionalized quantities, e.g. the non-dimensionalized energy *E*^{★} will be referred to as *E*.

The algorithm used in this work is fully described by Elsey *et al.* (2009). For the convenience of the reader, we briefly discuss it in §2. In this work, we present a large-scale three-dimensional simulation of grain growth, far beyond the scale of simulations presented previously. This simulation is performed on a parallelized version of the code executing on a large cluster. The results are analysed in great detail and compared with numerous theoretical predictions, results from experiments and smaller simulations performed with various numerical methods. The initial condition contains over 130 000 fully resolved grains and the evolution is simulated until only about 14 000 remain. We present various statistics collected throughout this simulation. In particular, we demonstrate the anticipated self-similar character of the grain-size distribution as it evolves in time. Furthermore, our results show good agreement with other three-dimensional predictions for grain growth, such as power-law growth of the mean grain volume, the three-dimensional version of the Aboav–Weaire law (Aboav 1970; Weaire 1974) and the Mullins extension of the two-dimensional von Neumann–Mullins relationship (Mullins 1989).

## 2. Algorithm

Our algorithm is an extension of *distance function-based diffusion-generated motion* (Esedoḡlu *et al.* 2010), which is in turn a variant of the *threshold-dynamics* scheme proposed by Merriman *et al.* (1994). The threshold-dynamics scheme reduces the computation of various geometric motions of an interface to the alternation of two highly efficient operations: convolution of a characteristic function representing the interior of the interface with a circularly symmetric kernel, and thresholding the convolution output to return to a characteristic function. On a uniform grid with *n* grid points, the complexity of these operations is just per time step. The method is also unconditionally stable, so that the only restrictions on the size of the time step are owing to accuracy considerations. A major drawback of the threshold-dynamics algorithm is that it is very inaccurate on uniform grids because the interface must be represented as the boundary of a characteristic function, disallowing any possibility of sub-grid accuracy in the absence of adaptive-grid refinement (as is explored in Ruuth 1998*b*).

A signed-distance-function-based algorithm for diffusion-generated geometric motions of the same types attainable by threshold dynamics is proposed by Esedoḡlu *et al.* (2010). This algorithm is a modification of the threshold-dynamics algorithm, replacing the characteristic function of the set used in threshold dynamics with a signed distance function to the boundary of the set, and the thresholding operation by a redistancing operation. The redistancing operation reconstructs the signed distance function to the boundary of the set from the convolution output, much as the thresholding operation reconstructs the characteristic function of the set from convolution output in threshold dynamics.

Unlike the characteristic function of a set, the signed distance function is a Lipschitz continuous function. This allows for sub-grid accuracy in locating the interface on a uniform grid. Furthermore, fast algorithms are known for constructing signed distance functions (e.g. Russo & Smereka 2000; Zhao 2004), so that the computational complexity of the signed-distance-function-based algorithm is still per time step. A more detailed discussion of the differences between Merriman *et al.* (1994) and Esedoḡlu *et al.* (2010) can be found in the latter.

The threshold-dynamics scheme was extended by Ruuth (1998*a*) to multi-phase motion by mean curvature, the case occurring in grain growth, but the same accuracy limitations faced by the original threshold-dynamics scheme still apply. Standard level-set techniques have also been applied to the grain-growth problem by Zhao *et al.* (1996) and Esedoḡlu & Smereka (2008). The level-set method is a sharp-interface method, similar to the signed distance function technique we use here (which may itself be considered a special type of level-set method). Standard level-set methods place only very loose restrictions on the details of the level sets used to implicitly represent the interface. For example, the size of the gradient |∇*ϕ*| of the level-set function *ϕ* is generally restricted from becoming too large or too small. This restriction is insufficient for multi-phase motion, as first noted and fully explained by Esedoḡlu *et al.* (2010). The standard boundary condition for grain growth is the Herring-angle condition (Herring 1951). For equal surface tensions, the condition states that triple junctions must meet at angles of 120^{°}. In level-set methods, an *O*(1) error is, in general, introduced in the motion of triple junctions owing to differences in the profiles of the level sets representing the three phases. This error also prevents the correct angle conditions from being obtained.

In contrast, Esedoḡlu *et al.* (2010) provide analytical justification and extensive numerical convergence studies showing that the correct behaviour is obtained at junctions under the multi-phase version of the signed-distance-function-based algorithm for motion by mean curvature introduced there. The algorithm for multi-phase motion is given in this work, with an enhancement allowing for the simulation of *N* grains using only *M*≪*N* signed distance functions described fully by Elsey *et al.* (2009). This is achieved by representing large collections of individual, well-separated grains with a single signed distance function. Easily implemented precautions are taken to prevent non-physical interactions from occurring between grains contained in a single signed distance function.

The signed-distance-function-based algorithm for multi-phase motion by mean curvature described in Elsey *et al.* (2009) and applied in this work can be thought of as a particular type of level-set method. In order to obtain the correct behaviour at triple junctions, there is a stronger restriction on the level sets than ordinarily enforced; however, other hallmarks of the level-set method are retained: the interfaces are implicitly represented (the algorithm operates only on the values of the signed distance function at grid points), sub-grid accuracy is obtained due to the Lipschitz continuity of the level sets, and the interface is sharp (unlike phase-field methods, in which the interface is represented by a diffuse transition layer). On the other hand, standard level-set methods for evolution by mean curvature require the solution of a degenerate, highly nonlinear PDE. In contrast, the present algorithm is unconditionally stable (allowing much larger time steps to be taken), a feature inherited from its close connection to the diffusion-generated motion of threshold dynamics.

The initial condition is taken to be the Voronoi diagram for *N*_{0} seeds placed uniformly at random in the computational domain *D*. The *N*_{0} initial grains are contained in sets *Σ*_{i}, *i*=1,…,*N*_{0}, with and . We partition the sets *Σ*_{i} into *M*≪*N*_{0} disjoint sets *Ξ*_{k}, such that . We maintain *M* signed distance functions *d*_{k}(*x*), giving the signed distance to the sets *Ξ*_{k}, with *d*_{k}(*x*)>0 for *x*∈*Ξ*_{k}. For example, in the grain-growth simulation described in §4, we have *N*_{0}=133 110 and *M*=64. The signed distance functions are updated using the following algorithm:

For , perform steps 1–4.

Compute for

*k*=1,…,*M*, whereConstruct

*B*_{k}(*x*) for*k*=1,…,*M*to remove overlaps and vacuums from the previous step,Construct the signed distance function for

*k*=1,…,*M*according toIf necessary, swap appropriate grains between signed distance functions to ensure that all the grains associated to a given signed distance function remain well separated. Denote the resulting signed distance functions as .

For a full description of the algorithm, the reader is referred to the work of Elsey *et al.* (2009). We use a simple cubic lattice with periodic boundary conditions to avoid introducing boundary effects into the evolution, except in the case that a single grain grows so large as to be adjacent to itself across a periodic boundary. Such a case cannot arise in normal grain growth until only very few grains remain (e.g. *N*(*t*)≪1000), which does not occur in the results reported here.

The computational complexity of our algorithm is formally , where *M* is the number of level-set functions used, and *n* is the total number of grid points in each set. For a single signed distance function, both the convolution step and the redistancing operation are . The algorithm is second-order accurate in space and first-order accurate in time away from triple points. At triple points, analysis and experiment by Esedoḡlu *et al.* (2010) suggest that the error is .

## 3. Previous work

A number of other numerical methods have been used to simulate grain growth in previous work. One of the best-known techniques is the Monte Carlo Potts model, implemented in Anderson *et al.* (1989). This model approximates curvature motion by a stochastic series of near-interface cell-flipping steps. While the basic Monte Carlo method is quite easy to implement, it is extremely slow and lacks sub-grid accuracy. Furthermore, the stochastic nature of the Monte Carlo evolution ensures that some type of averaging is needed to approximate the true continuum motion. For example, the evolution of a simple circle by mean curvature is very difficult to capture accurately using Monte Carlo methods, even on a well-resolved grid. Beyond these significant accuracy concerns, it is also difficult to connect the Monte Carlo method with some notion of ‘real’ time beyond reorientation attempts.

Front-tracking techniques have also been used to simulate mean curvature motion in both two (Kawasaki *et al.* 1989; Kinderlehrer *et al.* 2004) and three (Wakai *et al.* 2000; Kinderlehrer *et al.* 2004) dimensions. A major advantage of these techniques is computational efficiency, as computational resources are all devoted to the interface region. The fundamental difficulty inherent to this approach is managing the topological changes that abound in grain growth. With explicit representations of the interface, it is difficult to check if curves (in two dimensions) or surfaces (in three dimensions) intersect. In the case of two-dimensional mean-curvature multi-phase motion, it is expected (though not fully proven, see Mantegazza *et al.* 2004) that interfaces interact only through junction–junction collisions. If this conjecture is true, explicitly checking for and handling topological changes may be manageable in this case. However, *no such condition is expected to hold in three dimensions*. Pinch-off can occur with only two phases in three-dimensional mean-curvature motion, as in the standard ‘dumbbell’ example. Topological changes of this type are difficult to detect and manage using front-tracking techniques, especially in three dimensions. In contrast, our algorithm handles topological changes naturally, without any additional computation.

The phase-field technique has also been used extensively in simulations of grain growth (Fan & Chen 1997; Garcke *et al.* 1999; Krill III & Chen 2002; Kim *et al.* 2006; Suwa *et al.* 2008) and is much more similar to our algorithm than Monte Carlo or front-tracking techniques. In the phase-field method, interfaces are implicitly represented. However, there must be a wide transition region representing the interface between grains. For example, in Kim *et al.* (2006), the authors report that at least six grid points are needed in the transition layer to achieve acceptable accuracy. Therefore, a grain needs to be of the order of at least 25 grid points across to be moderately resolved. In contrast, we demonstrate in Elsey *et al.* (2009) that we can simulate evolutions quite accurately with grains approximately 10 grid points in each direction, and that we can track them down to half that length with only a few percent relative error. Though there is no rigorous notion of generalized solutions with uniqueness through topological changes, we also perform a convergence study in that work, showing that our simulations track grains through topological changes quite consistently as well. The large-grain-size requirement imposed by the phase-field model is a serious impediment to performing very large-scale simulations.

Our algorithm has all the major advantages of these various computational approaches. It is closely related to both threshold dynamics and level-set methods. From threshold dynamics, we inherit unconditional stability, allowing the choice of time step to be restricted only by accuracy considerations. From level-set methods, we inherit sub-grid accuracy and graceful handling of topological changes, free from user input. As with front-tracking models, we concentrate most of our computational resources at the interface, as the signed distance function computation (the most computationally intensive part of the algorithm) must be performed accurately only in a tubular neighbourhood of the interface (with width proportional to ). As a discretization of Mullins’s PDE-based continuum description of grain-boundary motion, our algorithm does not suffer from the difficulty of identifying physical time that plagues Monte Carlo-type models.

## 4. Normal grain growth

In this section, a large-scale simulation normal grain growth is described, and a wide variety of statistical measures on the resulting microstructure are reported. The simulation begins with 133 110 grains in the domain *D*=[0,82.306]^{3} with periodic boundary conditions. The evolution runs for time *t*=6.2021, allowing for over 90 per cent of the initial grains to disappear. *D* is discretized as a regular cubic lattice of size 512×512×512, and 300 time steps are taken, ensuring adequate spatial and temporal resolution. The coarsening rate is shown to agree well with theoretical predictions. The grain-size distribution is calculated, and exhibits self-similarity. This distribution is compared with a number of different predictions from the literature. In addition, average numbers of grain edges, faces and corners are computed and compared with other computational approaches and experimental data.

### (a) Qualitative microstructure

We present a three-dimensional simulation with an initial condition containing 133 110 grains. The initial condition was generated as described in §2. Figure 1*a* shows a single grain of average size taken from the simulation at *t*=6.2021. At this time, 〈*r*〉=2.11 and 〈*V* 〉=39.40, indicating that the mean grain radius has more than doubled the initial value and the mean grain volume is over nine times the initial mean grain volume. The grain appears to be very well resolved. Its faces, edges and corners are easy to see. The faces are smooth, and most appear to be concave. Thus, this particular grain, which is of average size at this stage in the evolution, must be growing owing to the curvature of its interfaces. Overall, the grain resembles what is observed in real materials, such as the beta brass grain shown in figure 1*b*.

The coarsening of the grain pattern is demonstrated in figure 2. Here, we display the grains intersecting the *x*=0, *y*=0 and *z*=82.306 surfaces in the initial condition (*a*) and after 300 iterations (*b*). Different colours correspond to different grains. By volume, the average grain at the end of the simulation is nearly 10 times as large as the average grain at initial condition. In figure 3, we show all the grains from 5 of the 64 total set functions *Ξ* at *t*=2.0674 and 6.2021. There is a great variability in the size of grains seen in this figure, from grains contained within a single grid cell (equivalent radius *r*≈0.08) that are about to disappear, to grains with radius *r*≈4.

Cross sections of successive slices at *t*=6.2021 are shown in figure 4*a*. Compare with figure 4*b*, showing results from a fully two-dimensional simulation, in which all three junctions must have 120^{°} angles. The cross-sectional views also feature more grains that are long in one dimension and short in the other as compared with the two-dimensional simulation results, where grains tend to be more regularly shaped.

### (b) Energetics

It is shown in Kinderlehrer *et al.* (2004) that, at least in the absence of topological changes, the surface energy *E* given by equation (1.3) decreases in time under mean curvature motion subject to the Herring-angle condition at junctions. It is natural to expect that *E* would continue to decrease even through topological changes (critical events). We verified that our numerical scheme respects this fundamental behaviour by evaluating the energy at every time step. We note that the energy *E* can be written in terms of the signed distance functions *d*_{k}(*x*) and the Dirac delta function *δ*(*x*) as the following sum of integrals over the computational domain *D*:
4.1
The factor of 1/2 arises as this formula counts each interface twice. We discretize *E* as
4.2
We use a first-order discretization of the delta function , following Smereka (2006). The energy *E* is measured at each time step and is found to decrease monotonically at every time step (see figure 5*a*), even as the evolution naturally handles the topological changes involved in the disappearance of over 100 000 three-dimensional grains through *t*=6.2021. The evolution of the number of grains is shown in figure 5*b*. After a short transition period (approx. *t*=0.4), the number of grains in the system decreases steadily. Note that even during this transition period, the energy of the system is decreasing quickly. The explanation for this transition period seen in the number of grains is that the initial condition is approximately Voronoi, and so there are very few small grains present initially, as demonstrated in the distribution shown in figure 6*a*. The system must evolve significantly before many grains are small enough to disappear.

### (c) Grain-growth rate and grain-size distribution

The average grain size, 〈*R*_{V}〉, and the grain-size distribution function, *f*(*R*_{V}/〈*R*_{V}〉), are probably the most important statistical quantities used to characterize an isotropic polycrystalline material. Texture distributions are also of primary importance for anisotropic polycrystalline materials, but texture is not considered in this work. Here, *R*_{V}=(3*V*/4*π*)^{1/3}, where *V* is the volume of a grain. Analytical approaches (Feltham 1957; Hillert 1965; Louat 1974), experimental results (as reported by Anderson *et al.* 1989) and simulation results (e.g. Anderson *et al.* 1989; Wakai *et al.* 2000) suggest that the average grain radius 〈*R*_{V}〉 exhibits power-law growth as a function of time: 〈*R*_{V}〉≈*Ct*^{n}, for *t* large. Analytically, the prediction *n*=1/2 has been made using a variety of considerations. The experimental results reported by Anderson *et al.* (1989) find 1/4≤*n*≤1/2. In their own simulation, the authors report that *n*=0.48±0.04 for fits to long-time data (obtained by discarding data from the initial transition phase of the simulation). In Wakai *et al.* (2000), the authors show approximately linear long-time dependence of 〈*R*_{V}〉^{2} on *t*. This simulation contains just 1000 grains initially, so the statistical precision of this measure is low. Furthermore, three-dimensional simulations via front-tracking require that explicit assumptions be made on the types of topological changes that can occur.

As normal grain growth is characterized by the self-similarity of the distribution of *R*_{V}/〈*R*_{V}〉, it follows that . In table 1, we fit 〈*V* 〉 to the function *at*^{b}+*c*, where *c*≈〈*V* _{0}〉 and mollifies the effect of the initial grain-size distribution on the fit. The fits are quite tight, with all reliability factors less than 0.7 per cent. Equating *b*=3*n*, we find that our simulation predicts the 95% CI of 0.501≤*n*≤0.518. The fit of 2.260*t*^{1.515}+3.496 to 〈*V* 〉 is plotted in figure 7.

The grain-size distribution function *f*(*R*_{V}/〈*R*_{V}〉) is defined by
4.3
In figure 6, we show histograms for this distribution at a variety of stages in the simulation. The distribution changes greatly throughout the evolution. The initial condition is approximately the Voronoi diagram for a randomly distributed set of points. The initial distribution of grain sizes is very narrow and sharply peaked. The distribution flattens out rapidly and appears to approach a self-similar state, characteristic of normal grain growth. This self-similar distribution appears to be attained by approximately *t*=4.1347 and is maintained thereafter, through over 10 000 grain disappearance events to the end of the simulation.

Another way to assess the self-similarity of the distribution of the grain-size distribution function across iterations is to look at the evolution of the central moments of the various distributions obtained. For these distributions, the first moment is, by definition, 1, and the first central moment is always 0. The variance and skewness (** E**[(

*X*−

**[**

*E**X*])

^{j}], for

*X*=

*R*

_{V}/〈

*R*

_{V}〉 and

*j*=2 and 3, respectively) are plotted in figure 8. These measures appear to be approximately constant for

*t*≥4.1347, agreeing with the visual impression of self-similarity obtained from figure 6

*c*.

Many closed-form distributions have been suggested as appropriate fits for the distribution *f*(*R*/〈*R*〉), including the Louat (1974) distribution, the Hillert (1965) distribution, the Rios (a modification of the Hillert distribution, Rios *et al.* 2006) distribution, the Weibull distribution (for two-dimensional grain growth, Fayad *et al.* 1999) and the lognormal distribution (for the distribution of grain radii in cross sections of three-dimensional experiments, Feltham 1957). These distributions are compared with the distribution of *R*_{V}/〈*R*_{V}〉 in figure 9*a*. The Rios distribution, with *ν*=3.34, appears to fit our simulation data the best. The lognormal and Louat distributions fit quite poorly, showing the wrong behaviour near *R*_{V}/〈*R*_{V}〉=0, for large *R*_{V}/〈*R*_{V}〉, and also peaking at *R*_{V}/〈*R*_{V}〉<1, all in disagreement with the simulation results. The Weibull and Hillert distributions show a better fit, but can be seen both visually and by reliability factor (table 2) to be inferior to the fit of the Rios distribution.

We also fit these distributions to data from cross sections of the three-dimensional simulation. This is of interest as, experimentally, it is difficult to slice materials thinly enough for the experiments to be two dimensional in nature, though carefully conducted thin film experiments are possible for polycrystalline grains of sufficient size. Measuring grain volumes is also difficult, though recent progress in X-ray and focused ion-beam techniques have made these measurements more feasible. However, it is still easiest to take cross sections of three-dimensional grains and measure areas and effective radii in cross section. Defining , where *A* is the area of a grain in cross section, we generate the distribution of *R*_{A}/〈*R*_{A}〉 from the simulation data at *t*=6.2021. We take 512 cross sections of constant *z*-value and aggregate the grain-slice area data across all these cross sections to create the simulation distribution. These cross sections contain a total of 368 138 two-dimensional grain slices. In figure 9*b*, we fit this distribution to the closed-form distributions discussed previously. None of these distributions fit the cross-sectional data as well as the Rios distribution fit the fully three-dimensional data taken from grain volumes. The Louat distribution fits the data the best with *α*=0.685, but with a reliability factor of *χ*=0.152. For comparison, the Rios distribution fits the three-dimensional data with *χ*=0.031. The distribution of *R*_{A}/〈*R*_{A}〉 is seen in figure 9 to be much flatter and wider than the distribution of *R*_{V}/〈*R*_{V}〉, reemphasizing the importance of interpreting these distributions separately.

### (d) Topology

Interesting topological characteristics of the grain network include the number of faces, corners and edges of individual grains in three dimensions, and the number of edges of grains viewed in cross section. Such characteristics have been the subject of numerous experimental studies (e.g. Desch 1919; Williams & Smith 1952; Rhines & Patterson 1982; Hull 1988; Zhang *et al.* 2004). Here, we compare the topological measures extracted from our large three-dimensional simulation to those obtained from experimental data as well as to those from other simulations. In all the following results, we take data from *t*=6.2021. At this time, 14 150 grains remain. In the 512 cross sections of constant *z*-value, there are a total of 368 138 grain slices.

Unlike front-tracking methods, methods using implicit representation of surfaces do not explicitly track topological features. The locations of these features are still well defined: a location *x* is on a face, edge or corner if, for every *ε*>0, there exist *m*=2,3 or 4, respectively, distinct subsets *c*_{1},…,*c*_{m} and locations *x*_{1},…,*x*_{m} satisfying |*x*_{i}−*x*|<*ε* and *d*_{ci}(*x*_{i})>0. The numerical implementation that allows association of topological descriptors to individual grains is described below. In order to count faces, corners and edges of individual grains at any fixed time *T* in the evolution, each grid point in the discretization is assigned a value from the set {1,…,*N*(*T*)} corresponding to the grain at that location. The number of faces of grain *i* is then the number of unique identifiers different from *i* contained in a 1-neighbourhood of the set of grid points that have identifier *i*. Counting corners is more challenging. In three dimensions, corners are characterized as being locations where four or more grains come together. We denote the set of all such locations as *C*. Because adjacent grid locations may, as part of a highly resolved corner, be marked as each being such a location, we take the number of connected components of *C* (as opposed to simply the number of points in *C*) within a single grain to be the number of corners possessed by that grain. However, this procedure will cause two corners connected by a short edge to be counted as one. To alleviate this problem, we subdivide the grid twice before applying the above procedure (so that a grid of size *n*×*n*×*n* is subdivided to size 4*n*×4*n*×4*n* before counting vertices). Having thus counted the number *f* of faces and the number *c* of corners as described above, we appeal to the well-known formula *c*−*e*+*f*=2 of Euler to infer the number *e* of edges of each grain. This formula holds for all polyhedra that are topologically equivalent to the sphere, which appears by inspection to be true for all the grains in our simulations of grain growth.

Data for the mean number of edges per face, mean number of faces and mean number of corners are presented in table 3 and compared with other simulations, with data reported for some regular polyhedrally based grain models and experimental results. The summary statistics vary to some extent with the simulation technique. Ours are well within the range of values found with other simulation techniques (though the other simulations were smaller and must be less statistically valid whether due to a smaller number of grains or the potential effects of ensemble averaging). Regular polyhedra such as the pentagonal dodecahedron and the tetrakaidecahedron have been proposed as space-filling approximations for grain shapes (Kelvin 1887; Smith 1952; Meijering 1953; Williams 1968), though experimentally it is well known that grains come in a variety of shapes and sizes. The tetrakaidecahedron matches the mean values we found well, but cannot explain more complex features of grain growth, such as the grain-size distribution function (4.3). The Voronoi model is generated by distributing seeds uniformly at random and growing crystals simultaneously and isotropically from these seeds. The Johnson–Mehl model grows crystals isotropically but allows for varying nucleation times (Meijering 1953). Both these models ignore grain-boundary motion due to interface curvature, holding grain boundaries stationary once crystals meet. These are in fact models for primary recrystallization, a different annealing phenomenon occurring when cold-worked metals are annealed. The experimental data contains a wide range of values, clearly demonstrating the difficulty of computing these measures in three dimensions and also suggesting that other higher-order effects (such as variable surface tension and mobility due to grain-boundary misorientation and inclination) play an important role in the evolution of polycrystalline grain systems. In future work, we will investigate extending our algorithm so that such effects can be simulated.

In figure 10, we plot the frequency with which grains with *f* faces occur. The distribution is skewed towards grains with many faces. The peak occurs at *f*=12 faces and the mean number of faces is 〈*f*〉=13.79. It is natural to expect that larger grains will have more faces, on average. However, the exact nature of this relationship is unknown. Figure 11 shows the relationship between the mean value of *R*_{V}/〈*R*_{V}〉 for grains with *f* faces and *f*, as determined from our simulation data. We also compare with measurements made by Rhines & Patterson (1982) on aluminium, by Zhang *et al.* (2004) on *α*-iron, and with simulation data generated by Anderson *et al.* (1989), using a Potts model and kinetic Monte Carlo techniques. The fit, particularly to the data for aluminium, is quite good and appears to describe the experimental data better than the linear fit posited by Anderson *et al.* (1989). The simulation results of Anderson *et al.* (1989) do appear to fit the measurements of Zhang *et al.* (2004) well for small *f*, but poorly for large *f*.

Stable corners occur where three triple lines come together on the surface of a grain. Under the assumption that every corner is stable, 3*c*=2*e*. Together with Euler’s formula, we can then calculate the number of corners and edges as a function of the number of faces *f* as *c*(*f*)=2(*f*−2) and *e*(*f*)=3(*f*−2). This prediction of a linear relationship between *c* and *f* is plotted in figure 12*a* against the values obtained from our simulation data, suggesting acceptable accuracy in our algorithm for counting corners and that our method does produce stable corners. Note that figure 10 illustrates that very few grains have less than 4 or more than 30 faces, so small inaccuracies in the count or the presence of only a few unstable corners will cause the small deviations from the prediction shown.

The three-dimensional version of the Aboav–Weaire law (Aboav 1970; Weaire 1974), proposed by Edwards & Pithia (1994), provides a relationship between the number of faces *f* exhibited by a grain and the mean number of faces of its neighbouring grains, *M*_{f},
4.4
where *μ*_{f} is the variance of *f*. Following Wakai *et al.* (2000), we plot the mean value of 〈*fM*_{f}〉 against *f* and find the linear relationship predicted by Edwards & Pithia (1994), but find the best linear fit to be *fM*_{f}=13.6*f*+25.4. This is in good agreement with the results of Wakai *et al.* (2000) who found *fM*_{f}=13.3*f*+23.4. Based on their experimental data, Zhang *et al.* (2004) found *fM*_{f}=13.97*f*+12.61. Equation (4.4) predicts *fM*_{f}=12.8*f*+37.7, using the values of 〈*f*〉 and *μ*_{f} determined by our simulation data. Thus, simulation, experiment and theory for the three-dimensional Aboav–Weaire law agree well up to an additive constant. See figure 12*b* for simulation data and best-fit line.

In two dimensions, the well-known von Neumann–Mullins relationship (Mullins 1956) states that grains with more than six sides grow, and grains with fewer than six sides shrink,
4.5
where *n* is the number of sides of the grain. Mullins (1989) proposed the following relationship for three dimensions, relating the mean growth rate of three-dimensional grains to their number of faces, *f*:
4.6
where
4.7
and
4.8
In figure 12*c*, we plot the simulation results for 〈(*dV*/*dt*)/*R*_{V}〉, taken from 5.9953≤*t*≤6.2021. For *t*=6.0367, 6.0780, 6.1194 and 6.1607 and *δt*=0.0413, we approximate *dV*/*dt*=(*V* (*t*+*δt*)−*V* (*t*−*δt*))/(2*δt*). The simulation results follow the same curve as the predictions, but appear to differ by a constant additive value of approximately 2.2. Our simulation results agree well with those of Wakai *et al.* (2000) (using Surface Evolver, a front-tracking software package), and of Weygand & Bréchet (1999) (via vertex dynamics). Other topological generalizations have been proposed by Weaire & Glazier (1993) and Hilgenfeldt *et al.* (2001). The Weaire relationship gives a linear relation between *f* and 〈(*dV*/*dt*)/*R*_{V}〉, which does not appear to fit the data presented here or in other simulations well. The Hilgenfeldt relationship agrees closely up to a scaling constant with the von Neumann–Mullins extensions for *f*≥10 and is thus not shown. Recently, MacPherson & Srolovitz (2007) gave a generalization of the von Neumann–Mullins relationship to three dimensions; however the quantities involved in their formula (mean width and total edge length) are not topological in nature, unlike the two-dimensional von Neumann–Mullins relation. Furthermore, mean width is quite difficult to calculate for grains. Simplifications are known for convex polyhedra (Cahn 1967) and for regular polyhedra (Glicksman *et al.* 2009), but grains are irregular and may possess both convex and concave faces. We elect to compare only with the Mullins generalization, which is a topological relationship depending only on the number of faces *f*.

## 5. Summary and conclusions

We apply the algorithm developed by Esedoḡlu *et al.* (2010) and extended by Elsey *et al.* (2009) based on diffusion-generated motion of signed distance functions to a three-dimensional simulation of grain growth. This approach naturally captures the Herring condition at triple junctions. In addition, numerical evidence unequivocally shows that the energy of the simulated system decays, even through topological changes. The efficiency of this algorithm allows us to compute the *accurate* evolution of over 130 000 grains until less than 15 000 grains remain. To the best of our knowledge, this evolution contained at least twice as many grains as any other currently published to date. In the next largest simulation (Suwa *et al.* 2008), the authors implement a phase-field model initially containing 50 000 grains on a 512×512×512 grid. Grains have an average initial size of approximately 14×14×14 grid points, with a diffuse interface width *ϵ* of three grid points. This suggests that the initial resolution of their system is quite low. We are able to verify, with greater confidence, that the coarsening rate for normal grain growth is 〈*r*〉∼*t*^{1/2} and that the grain-size distribution function is self-similar. We are also able to provide accurate average values of the number of edges, corners and faces of individual grains. We observe that in many cases, these are in agreement with experimental results. This provides further validation that approximately normal grain growth is present in experimental settings.

## Acknowledgements

This work was supported, in part, by grants from the National Science Foundation: DMS-0748333 and DMS-0810113. S.E. was also supported by an Alfred P. Sloan Foundation fellowship.

- Received April 9, 2010.
- Accepted June 21, 2010.

- This journal is © 2010 The Royal Society