Application of the level-set method to the analysis of an evolving microstructure

Computational Materials Science 85 (2014) 46–58

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Application of the level-set method to the analysis of an evolving microstructure C.-L. Park a, P.W. Voorhees b, K. Thornton a,⇑ a b

Department of Materials Science and Engineering, University of Michigan, 2300 Hayward Street, Ann Arbor, MI 48109, USA Department of Materials Science and Engineering, Northwestern University, 2220 Campus Drive, Evanston, IL 60201, USA

a r t i c l e

i n f o

Article history: Received 25 September 2013 Accepted 10 December 2013 Available online 11 January 2014 Keywords: Smoothing method Level-set method Interfacial properties Interfacial curvature Microstructural morphology Three-dimensional microstructure

a b s t r a c t Examination of three-dimensional microstructures and their evolution requires accurate quantification of morphologies. We developed a smoothing algorithm, termed ‘‘level-set smoothing,’’ that smoothes a discontinuous or nearly discontinuous function that describes interfaces. The proposed method can be applied to any voxel-based data describing a two-phase structure. The level-set smoothing method is a set of sequential data-processing schemes that consists of first generating the signed distance function for the given microstructure using the level-set method, followed by smoothing via diffusion. This algorithm is designed to substantially reduce the artificial shift in the interfacial positions, as compared to simple smoothing methods, so that interfacial properties, including their rate of change, are accurately calculated. The accuracy of the proposed method is demonstrated by comparisons of the numerically calculated interfacial properties from data sets smoothed by the proposed method with their analytical solutions and of the interfacial locations of a complex microstructure from unsmoothed and smoothed data sets. Higher accuracy in calculating time derivatives of curvatures is achieved by using the advective method and choosing the time difference such that interfacial displacement is about two to three Cartesian grid points. Ó 2014 Elsevier B.V. All rights reserved.

1. Introduction The ‘‘processing–structure–property’’ triangle is the foundational concept in materials science and engineering. During material processing, the microstructure can undergo both morphological and topological changes. This, in turn, directly affects the material properties. One can exploit this relationship by carefully controlling the manufacturing processes that transform the microstructure to induce desired properties. For example, it is a well-known fact that the mechanical properties of metallic alloys are greatly influenced by the shape as well as size and spatial distributions of precipitates [1]. Furthermore, there is an increasing number of claims that indicate the importance of microstructure to the performances of renewable energy devices such as the solid oxide fuel cells [2–4] and Li-ion batteries [5,6]. In order to study microstructures to the full extent, it is often necessary to analyze three-dimensional (3-D) microstructures. Advances in microscopy and the availability of large-scale computing resources have made the studies of 3-D microstructures more accessible [7–9]. Experimental techniques such as mechanicalserial-sectioning optical microscopy [10–13], X-ray tomography

⇑ Corresponding author. Tel.: +1 7346151498. E-mail address: [email protected] (K. Thornton). 0927-0256/$ - see front matter Ó 2014 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2013.12.022

[14–20], and dual-beam focused-ion-beam scanning electron microscopy [2,21–23] can now provide detailed 3-D microstructures. In addition, computational tools such as phase-field method [19,24–27] and level-set method [28–30] can be used to simulate the evolution of 3-D microstructures. As such, the field of 3-D materials science has become increasingly invigorated in recent years. One of the enabling tools for 3-D materials science is new characterization methodologies for morphologically and topologically complex microstructures. As wealth of 3-D data sets of microstructures become available through the improved experimental and computational tools, novel techniques for characterizing local microstructural characteristics such as interfacial curvatures have been developed. For example, an interfacial shape distribution (ISD) can represent the complex morphology of a microstructure by the probability distribution of interfacial area as a function of principal curvatures [13,26,31]. Other characterization methods involve computing spatial correlation for interfacial curvatures [32], calculating genus to quantify the topology of microstructures [33], and determining the channel size distribution of bicontinuous microstructures [34]. In some situations, interfacial curvatures are the preferred measurements for characterizing complex microstructures due to their integral role in microstructural evolution. For example, the kinetics of coarsening is dictated by the curvature of the interface.

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58

Furthermore the time derivatives of curvatures can help elucidate the evolution of morphology and topology. It is thus clear that the ability to numerically calculate the interfacial properties such as curvatures and their time-derivatives (rates of change) from 3-D data sets is essential in developing an understanding of microstructural evolution. Several numerical schemes for computing interfacial curvatures have been developed in the past. These numerical methods can be largely classified into two main groups: surface-mesh-based methods and voxel-based methods (voxels are 3-D counterpart of pixels, which stands for volumetric pixels). The surface-mesh-based methods compute curvatures on a 3-D network of triangulated mesh that represent interfaces using the positional information of the mesh vertices. On the other hand, voxel-based methods compute curvatures using the values of a function at voxels that are close to the interfaces, which are typically defined as regions where the values of the function at the voxels are between the two bulk phase values. Most triangulated-mesh-based methods can be further categorized as either a ‘‘surface fitting’’ method or a ‘‘discrete’’ method [35]. A surface fitting method determines the curvature at each vertex by obtaining a quadratic (or cubic for higher accuracy) equation for a parameterized surface that best fits the interface [36–38]. On the contrary, discrete methods compute curvatures directly using formulas derived from differential geometry, such as Gauss-Bonnet theorem [38–40] and Euler theorem [38,41,42], without analytically fitting the surface. Similarly, voxel-based methods can also be divided into two subgroups: methods that require spatial derivatives across the interface, such those involving the level-set method [43,44], and methods that do not require any spatial derivative [45,46]. In this paper, we focus on the derivative-based approach combined with the level-set method for computing curvatures since the methodology is straightforward and is also computationally efficient, which is particularly advantageous with large 3-D data sets. However, there is one major challenge in computing interfacial properties, such as curvatures, using the level-set method. Numerical calculations of curvatures require the evaluation of high-order derivatives across the interface, which become inaccurate when interfacial resolution is low. Experimental microstructure data often describe the microstructure through a contrast difference between two phases. In such cases, the interfaces between two phases are described by few pixels or voxels and thus the reconstructed interfaces form a ‘‘wedding cake’’ structure (having stepped surface rather than smooth) due to the sharp change of contrast across the interfaces. For the phase-field approach, order parameters with smooth transitions describe the phase boundaries, but they typically contain only two to five Cartesian grid points, which presents problems in the quantitative analysis of the interfacial properties. A simple solution to the sharp transition of the microstructural data across interfaces is to apply the volume smoothing technique. This method simply averages the value of data over the neighboring Cartesian grid points with a set width. However, smoothing the data sufficiently for evaluating high-order derivatives using this approach leads to substantial shift in the interfacial position. As a result, interfacial properties, such as curvatures, calculated from the volume-averaged data may not accurately represent the morphology of the microstructure or its evolution. Therefore, we have developed an advanced yet computationally inexpensive smoothing algorithm, which we term ‘‘level-set smoothing,’’ that can be applied to any voxel-based data describing a two-phase structure. The level-set smoothing method is a set of sequential data-processing schemes that consists of first generating the signed distance function for the given microstructure using the level-set method, followed by smoothing via diffusion. The proposed smoothing method can facilitate accurate calculation of

47

interfacial properties with minimal displacement of interfaces during the smoothing process, enabling quantification of morphology and its evolution for a 3-D microstructure. In this paper, we present the algorithm for the level-set smoothing method and demonstrate its effectiveness in calculating various types of interfacial properties through series of validation tests. We also list numerical algorithms that can be used to calculate various interfacial properties, including time-derivative quantities. For the time derivatives of curvatures, we present two methods, advective and convective, for a comparison. 2. Methods 2.1. Level-set smoothing method The level-set iterations generate the signed distance function near the interfaces for the microstructure. Since the level-set iteration only smoothes out the noise in the first derivative, we apply diffusion smoothing as the second step to remove noise in higher derivatives. The details for each of these processes are described below. From here onwards, a ‘‘grid point’’ refers to a voxel on a Cartesian-grid system. 2.1.1. Level-set method In the level-set approach, an interface is represented by a contour where a level-set function, u, takes a predetermined value, typically zero. The level-set method has been applied to study wide range of microstructural evolution [28–30,47,48]. In many cases, the level-set function is assumed to be a signed distance function, where the magnitude of the function is given by the distance from the nearest interface, while the sign indicates one phase or another. Given a function that is positive in one phase and negative in another phase, the signed distance function can be calculated by iterating the following equation until converged:

@u ¼ SðuÞð1  jrujÞ; @s

ð1Þ

where SðuÞ is a smoothed sign function, which will be defined later in this section, and the time, s, is an iteration time that is not associated with any physical evolution. Eq. (1) will be referred to as the level-set equation hereafter. A signed distance function representing the microstructure can be obtained from microstructural data from computer simulations or experiments. To use simulated data, described by the order parameter, / (note the difference in the font between / and u), in the levelset equation, we prepare the data as follows to best preserve interfacial location during the smoothing process. We begin with /, which describes a two-phase structure where the phases are defined by bulk values, such as / ¼ 0 or 1, and the phase boundary defined as the region where / transitions between the bulk values. Since the zero level-set defines the interface, the value of order parameter needs to be adjusted by adding or subtracting a constant if the interfacial value of / is nonzero. Additionally, we found that Eq. (1) converges more quickly and accurately with a minimal shift of the interface, if the range of / is adjusted so that the magnitude of the gradient is approximately one near the interface (which is the steady-state condition of Eq. (1)). For example, if the input order parameter, /, has the bulk phase values defined at 0 and 1 while the phase boundary has a nondimensional length of d, then / needs to be translated and scaled such that it varies from d=2 to d=2. This adjusted order parameter then becomes the initial condition, u0 , and has the property jru0 j ’ 1 near the interface. The level-set equation, Eq. (1), is a modification of re-initialization equation developed by Sussman et al. [49]. The discretized

48

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58

version of Eq. (1) in a 3-D Cartesian grid system, where i; j, and k denote the indices for x; y, and z positions, respectively, is N N N uNþ1 ijk ¼ uijk þ Ds  Sðuijk ÞGðuijk Þ;

ð2Þ

where Sðuijk Þ is the numerical implementation of the smoothed sign function (explained below), and Gðuijk Þ is the numerical approximation for 1  jruj at a grid point (i; j; k). Following Sussman, Gðuijk Þ is calculated using a first-order upwind scheme [49]. This monotone scheme possess improved stability over standard central differencing and also ensures that u near the interfaces converges to the steady-state first. While we adopt the same numerical scheme to calculate 1  jruj as Sussman’s, our numerical implementation of the sign function is modified such that shift in the interfacial position is reduced. The smoothed sign function, SðuÞ, used in our levelset formulation is defined as

Sðuijk Þ ¼ uijk

.qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u2ijk þ c2 ;

ð3Þ

where the numerical parameter, c, is used to regularize the singularity in the numerical implementation of SðuÞ. There are two differences between our smoothed sign function and the sign function adopted by Sussman: the smoothed sign function can vary during the level-set iteration and uses much larger value of c. For a typical sign function, a very small value of c is used so that the sign function returns a value of either approximately 1 or 1. We employ a larger value of c, one that is comparable to the range of u0 , to smooth the sign function; hence we refer to SðuÞ as a ‘‘smoothed sign function.’’ These two modifications ensure that the effective time step, Ds  SðuNijk Þ, actively changes as u approaches the signed distance function near the interface. The effective time step is smaller in the region near the interfaces, where the level-set equation converges faster and higher accuracy is needed, and is larger in the region further away, where the equation converges slower and less accuracy is needed. However, even with the adjustment of the input such that jru0 j ¼ 1 at the interface and the use of the smoothed sign function that reduces the effective time step for grid points near interfaces, the interfacial location can shift slightly in high curvature regions (usually at distances smaller than a voxel, which is still undesired). One way to mitigate this problem is by limiting the number of iterations of the level-set method. Calculation of interfacial properties requires the level-set function in the form of a signed distance function only within a certain number of grid points near the interface. For example, the number of grid points, N, near the interface required for calculating curvature is determined by the width of the stencil used to calculate the high-order derivatives. With the numerical scheme adopted, N ¼ 6Dx (3Dx on both sides of the interface) is sufficient to calculate all of the interfacial properties discussed in this paper, including timederivative quantities. To provide adequate iteration time for the level-set function to converge to the signed distance function in the region within 3Dx of the interfaces, the level-set equation is iterated until the maximum magnitude of the level-set function reaches 6Dx (twice the distance over which the distance function is desired), as shown in Fig. 1a. Once the iteration is complete, the level-set method produces a signed distance function, u, near the interface. 2.1.2. Diffusion smoothing method The signed distance function obtained above contains a significant amount of numerical noise in the second derivative. To remove noise in the second derivative, diffusion smoothing, achieved by evolving the standard diffusion equation, is applied to the level-set function. The dimensionless diffusion equation is given by

@u ¼ r2 u; @s

ð4Þ

where s is a nonphysical iterative time. Over-iteration of diffusion smoothing can significantly alter the location of the interfaces. This is avoided by limiting the diffusion smoothing to the bare minimum required to eliminate noise in the second derivative. The appropriate amount of diffusion smoothing depends on the quantity to be calculated. Calculation of local interfacial properties, such as curvatures, requires a small number of iterations for diffusion smoothing; we use 20 iteration steps with time-step of 0:1ðDxÞ2 . On the other hand, calculation of time-derivative quantities, such as the time derivatives of curvatures, need approximately 30 iterations of diffusion smoothing. These values were determined by the numerical error analysis; the detailed results are omitted due to space constraints. While the interfacial positions can be preserved by limiting the diffusion smoothing, the level-set function near the interfaces can deviate from the signed distance function during this process. We have found that this deviation is nearly uniform over the grid points required to calculate interfacial properties. To remedy this problem, after applying the diffusion smoothing, we compute jruj at all grid points and determine its average, hjruji, where juj < 1Dx, i.e. the region where the level-set function is closest to the distance function after diffusion smoothing, and scale the level-set function with the calculated 1=hjruji. After normalization, the diffusion smoothed level-set function again takes the form of signed distance function. Hereafter, to simplify notations, /; u and uD will be used to denote the order parameter, the level-set function, and the smoothed level-set function that has been processed by the diffusion smoothing and normalized by hjruji. The combination of the level-set function with diffusion smoothing will be referred to as level-set smoothing. 2.2. Structures The three types of structures are employed in this paper to validate the level-set smoothing method: a sphere, a cylinder and a two-phase bicontinuous structure. All structures presented in this paper are constructed in a nondimensional space, and therefore all measurements of length and time are unitless. A grid spacing of Dx ¼ Dy ¼ Dz ¼ 1:0 is employed for simplicity. Order parameters describing spherical and cylindrical geometries can be constructed analytically using a hyperbolic tangent function. For a spherical particle in a Cartesian grid system, the discretized order parameter, /ijk , where i; j, and k denote the indices for the position, can be constructed as

/ijk ¼

   2 rijk  R 1 1 ; þ tanh 2 2 dtanh

ð5Þ

where

rijk ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 ði  xc Þ þ ðj  yc Þ þ ðk  zc Þ ;

ð6Þ

and R is the radius of the spherical particle. The variables xc ; yc ; zc represent the Cartesian coordinate of the center of the sphere, which is placed at the center of the computational domain. The order parameter defined by Eq. (5) will have a bulk value of 0 for the inside of the sphere and 1 for the outside while dtanh represents the nondimensional interfacial thickness, which is defined as the region where the value of / smoothly transitions between the two bulk values. In this paper, we used dtanh ¼ 4:0, which results in about four grid points across the interface and is typical of phasefield simulations. Similarly, a cylindrical tube with an infinite length along the z-axis can be constructed by Eq. (5) with rijk redefined as

rijk ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ði  xc Þ þ ðj  yc Þ ;

ð7Þ

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58

(a)

49

(b)

Fig. 1. (a) The interfacial profile of the Allen–Cahn structure (details of the structure are discussed in Section 2.2), comparing the input function from phase-field order parameter and the level-set function obtained using the scheme described Section 2.1.1. The solid line shows the phase-field order parameter describing the two phases, scaled such that the magnitude of the gradient at interface is unity. The dotted line shows the level-set function whose value near the interface (3Dx on each side) is the signed distance from the nearest interface. (b) Morphology of the bicontinuous Allen–Cahn structure at the simulation time of t0 ¼ 600.

where xc and yc represent the Cartesian coordinate of the axis of the cylinder. To calculate the time-derivative quantities, pairs of spherical particles and cylindrical tubes with different radii were analytically constructed from the order parameters. While above geometries provide important test cases with well defined analytical solutions to interfacial properties, the true power of the method is it ability to be applied to a wide range of complex microstructures. For demonstration, we choose a bicontinuous structure resulting from phase ordering [50]. Using the phase-field method, a two-phase bicontinuous structure can be simulated via nonconserved Allen–Cahn dynamics. Hereforth, this bicontinuous structure is refered to as the Allen–Cahn structure. The governing equation of Allen–Cahn dynamics is [50]

  @/ @f ¼ L/   2 r2 / ; @t @/

ð8Þ

where the bulk free energy, f ð/Þ, assumes a double-well potential,

f ð/Þ ¼

W 2 / ð1  /Þ2 ; 4

ð9Þ

which is often used for the Cahn–Hilliard equation and has minima at / ¼ 0 and 1 (corresponding to two bulk phases) [51]. Variables L/ ; W;  are the mobility coefficient, the height of the double-well potential, and the gradient energy coefficient, respectively. The value of parameters are assigned as follows: L/ ¼ 1:0; W ¼ 0:8; 2 ¼ 0:4, which leads to an interfacial thickness of d ¼ 4. Starting from a seed of random numbers that varies between / ¼ 0:5  0:1, the bicontinuous Allen–Cahn structure emerges. Fig. 1b shows the morphology of the Allen–Cahn structure at nondimensional simulation time t 0 ¼ 600. 2.3. Calculation of the interfacial properties with convective and advective method In order to study the morphology of a complex microstructure, curvatures of the surface, such as the mean and Gaussian curvatures, need to be evaluated. Furthermore, to study the evolution of microstructures, calculations of time derivatives of relevant quantities are necessary to quantify the evolution of interfaces. These time-derivative quantities include the normal velocity, which is the rate of displacement of interfaces in the direction normal to the interface, and the time derivatives of mean and Gaussian curvatures (calculated following the interface), which quantify the evolution of local curvatures. In this section, we briefly explain how the mean curvature, H, the Gaussian curvature,

K, and the normal velocity, v, can be calculated using the level-set approach. We then introduce two methods, the convective and the advective method, which are used to calculate time derivatives of mean and Gaussian curvatures, DH=Dt and DK=Dt, respectively and defined below. All interfacial properties are calculated from the smoothed level-set function, uD . Since uD takes the form of a signed distance function near the interface, the interfacial properties, H; K, and v can be evaluated from the derivatives of uD using the level-set method [43]. The main advantage of using the smoothed level-set function is that the interfacial properties can be accurately computed on the Cartesian coordinate without parametrizing the surface of the microstructure. Once an interfacial property is calculated at all grid points, the physically meaningful interfacial values are obtained by linear interpolation at the interfacial position where uD ¼ 0. The expressions for H; K, and v can be found in Sethian’s book [43] and are included in A along with the finite differencing methods used to evaluate the spatial derivatives that appear in the expressions. Note that our definition of H varies from [43] by a factor of 1=2. We now present the two numerical methods to evaluate DH=Dt and DK=Dt, each based on a different frame of reference (FOR): Eulerian with fixed FOR and Lagrangian with material FOR. The first approach is the convective method [52]. The convective method uses an Eulerian approach to calculate the Lagrangian derivative D=Dt. For instance, DH=Dt can be expressed in terms of an Eulerian derivative @H=@t and a gradient, rH, as

DH @H ¼ þ v  rH: Dt @t

ð10Þ

Here, v is the normal velocity vector defined as the product of the normal velocity, v (Eq. (A.9)), and the unit normal vector, n (Eq. (A.3)). The convective method calculates each term on the right hand side of Eq. (10) at all grid points, where H and v are computed using the level-set approach. The partial time derivatives of the mean curvature in Eq. (10) is an Eulerian derivative and can be discretized as

@Hðxi ; yj ; zk ; t 1 Þ Hðxi ; yj ; zk ; t2 Þ  Hðxi ; yj ; zk ; t 1 Þ ¼ ; @t Dt

ð11Þ

where t 1 and t 2 are different evolution times, Dt ¼ t 2  t 1 , and xi ; yj , and zk are the coordinates in Cartesian grid, where i; j, and k are the indices associated with grid points along the x; y, and z directions, respectively. Values of Hðxi ; yj ; zk ; t1 Þ and Hðxi ; yj ; zk ; t 2 Þ at each grid point are calculated from the smoothed level-set function at t 1

50

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58

and t 2 , respectively, based on Eq. (A.4). The second term, v  rH, in Eq. (10) is referred to as the convective term and is calculated using a first-order upwind scheme, where the normal velocity vector, v , is calculated as v n, and the gradient of H is evaluated at t 1 ; rHðxi ; yj ; zk ; t 1 Þ. The convective term provides the correction between the Lagrangian derivative ðD=DtÞ, also known as the convective derivative or the material derivative, and the Eulerian derivative ð@=@tÞ, which does not take into account the movement of interfaces. Once the sum of the partial derivative and the convective term is evaluated at all grid points, the interfacial values of DH=Dt are determined by linear interpolation at the interfacial position where uD ðt 1 Þ ¼ 0. The main advantage of the convective method is that the entire calculation is performed at stationary grid points, which circumvents the need to explicitly follow the interfacial movement between time steps. The second approach for calculating the time derivatives of curvatures is the advective method. Unlike the convective method, the advective method requires the explicit determination of the location of interfaces at two consecutive time steps. The advective method can be viewed as directly calculating DH=Dt and DK=Dt on a FOR that is moving with the interface. The time derivatives of the mean curvature, DH=Dt, can be expressed as

DH Hðxðt 2 Þ; t 2 Þ  Hðxðt1 Þ; t1 Þ ðxðt 1 Þ; t 1 Þ ¼ ; Dt Dt

ð12Þ

where

xðt2 Þ ¼ xðt 1 Þ þ v Dt:

ð13Þ

Here, xðt1 Þ is the interfacial coordinate at reference simulation time t 1 ; and xðt2 Þ is the predicted interfacial coordinate of the point xðt1 Þ at time t2 , which is obtained by displacing xðt 1 Þ by v Dt (see Eq. (13)). While the convective method is more straightforward because the time derivatives of curvatures are computed on stationary grid points, the advective method yields more accurate results as demonstrated below. There are two reasons behind this. As an example, we here consider the calculation of DH=Dt. The first source of error in the convective method is in the term v  rH. The convective term requires calculation of rH, and the higher order derivative introduces additional numerical errors. The second source of error arises from the calculation of Hðxi ; yj ; zk ; t2 Þ near xðt1 Þ. For an accurate calculation of Hðxi ; yj ; zk ; t2 Þ; uD ðt 2 Þ needs to be in the form of the signed distance function. However, if the interfacial displacement between two time steps is large, uD ðt 2 Þ at grid points near xðt1 Þ may not be a signed distance function (since uD ðt 2 Þ is in the form of the signed distance function only within 3Dx from xðt2 Þ). Error in Hðxi ; yj ; zk ; t2 Þ leads to inaccuracy in @H=@t. Hence, the accuracy of the convective method quickly deteriorates as the interfacial displacement between two consecutive time steps increases. A detailed study of the errors associated with the calculation of interfacial properties are presented in the subsequent section. 3. Results and discussion The following section demonstrates the capability of the proposed level-set smoothing method for calculating various types of interfacial properties. As mentioned earlier, we aim to achieve (1) accurate calculations of interfacial properties, (2) minimal displacement of interfaces during the smoothing process, and (3) to demonstrate the method with complex microstructures. As such, we divided this section into four parts:  Validation and uncertainty quantification of the level-set smoothing for simple geometries.  Validation of the calculation of time derivatives of curvatures for simple geometries.

 Comparison of interfacial locations before and after smoothing for a complex microstructure.  Application of the method to a complex microstructure. 3.1. Validation and uncertainty quantification of the level-set smoothing for simple geometries Generally speaking, validation of a numerical method can be accomplished by comparing the numerically calculated quantities with corresponding analytical solutions. To carry out such comparisons for the level-set smoothing method, we must compare the numerically calculated interfacial properties from the smoothed level-set functions with corresponding analytical values. Therefore, we use pairs of spheres and cylinders with different radii (to mimic the morphologies of shrinking objects at two different times). The simple geometries of spheres and cylinders are ideal for obtaining analytical values of interfacial properties and their changes. The construction of these structures using the hyperbolic tangent function, presented in Section 2.2, ensures near perfect spherical and cylindrical geometries, minimizing errors in the geometries themselves. As an initial analysis, we used pairs of spheres and cylinders with R1 ¼ 40 and R2 ¼ 38. Fig. 2a shows the interfacial profiles of the pair of spherical particles. The same smoothing sequence was applied to the four order parameters. Before the smoothing process, we translate the values of the order parameter, /, via

u0 ¼ 4/  2:

ð14Þ

After the translation, the inside and outside the structures are defined as u0 ¼ 2 and u0 ¼ 2, respectively (see Fig. 2b). Throughout p the smoothing sequence, Dx ¼ Dy ¼ Dz ¼ 1:0, Ds ¼ 0:1, ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c ¼ 3 ðDx2 þ Dy2 þ Dz2 Þ in Eq. (3), and periodic boundary condition was used. The level-set equation is iterated until the maximum juj reaches 6Dx. Subsequently, the level-set function is processed via diffusion smoothing and is scaled by 0.9 (hjruji near the interfaces after the diffusion smoothing step) to acquire uD , as shown in Fig. 2b. After obtaining the smoothed level-set functions, we calculated various interfacial properties discussed in Section 2.3. For the calculation of time-derivative quantities, the structures with larger radii are considered the reference point at t1 . For this particular setup, the physical time, Dt, required for the calculation of timederivative quantities can be arbitrarily chosen since it cancels out during the quantification of the fractional error of the numerical calculation with respect to the analytical solution, as discussed below. The numerically calculated interfacial properties are then compared with the analytical solutions to evaluate the fractional errors according to

ðf Þ ¼

fnumerical  fanalytical ; fanalytical

ð15Þ

where f is one of the interfacial properties. The comparison of the numerical calculation with the analytical solution is performed at every vertex that forms the surface mesh of the spheres and the cylinders. Fig. 3 displays the surfaces of the spherical particle colored by the magnitude of fractional errors of various interfacial properties calculated from / and uD . The grid effect, which manifests itself as rapidly varying errors on the surface, appears for both types of data sets. This effect originates from the grid-anisotropy of the underlying Cartesian grid employed. However, the range of the magnitude of percentage errors (as noted by the color bar range) shows that the level-set smoothing method significantly improves the accuracy of the calculated interfacial properties, especially for the time-derivative quantities.

51

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58

(a)

(b)

Fig. 2. (a) A center cut showing the interfacial profiles of the spherical particles of radius 40 and 38 given by /1 (solid line) and /2 (dashed line), respectively. The outside of the sphere is defined as / ¼ 1 and the inside as / ¼ 0. (b) A center cut showing interfacial profiles of the spherical particle with radius of 40 described by u01 (solid line), before smoothing, and uD1 (dashed line), after smoothing. The outside of the sphere is defined as u0 ¼ 2 and the inside as u0 ¼ 2.

The differences in errors from pre-smoothed and postsmoothed data sets do not vary much for the curvatures (as shown in Fig. 3a–d, about a factor of 2). However, for the time-derivative quantities, the magnitudes of errors reduce drastically after smoothing, by as much as a factor of 50, as shown in Fig. 3e–j. The large magnitude of errors in the calculated time derivatives of curvatures from pre-smoothed data sets can be explained by identifying the sources of errors for DH=Dt and DK=Dt. Based on Eq. (15), with f ¼ DH=Dt for an example,

ðDH=DtÞ ¼ ¼

H2 ð1þðH2 ÞÞH1 ð1þðH1 ÞÞ Dt H2 H1 Dt

1  H2DH t

H2 ðH2 Þ  H1 ðH1 Þ ; H2  H 1

ð16Þ

where Dt ¼ t 2  t 1 , and Hi and ðHi Þ are the interfacial mean curvatures and their associated fractional errors at simulation time t i (i ¼ 1 and 2), respectively. Following Eq. (16), the upper-bound error of jðDH=DtÞj; supjðDH=DtÞj , becomes

supjðDH=DtÞj ¼

maxðjðH1 Þj; jðH2 ÞjÞðjH1 j þ jH2 jÞ 2jðH1 Þj  jDHj ; jDHj jH j

ð17Þ

1

where DH ¼ H2  H1 and the right-hand side is based on the assumption that jðH1 Þj ¼ jðH2 Þj and H1  H2 during Dt. Eq. (17) shows that the upper-bound error is observed when ðH1 Þ and ðH2 Þ have opposite signs. Moreover, supjðDH=DtÞj is amplified if the difference in quantities, DH, is much smaller than the quantities from which the differences are calculated, H1  H2 . Without the application of the level-set smoothing method, interfaces of the pre-smooth data have noise as a result of low interfacial resolution. The unsmoothed order parameter across the interfaces produces continuous over/under-estimation of curvatures along the surface. Hence, there is a greater probability that ðH1 Þ and ðH2 Þ have inconsistent signs, which ultimately increases jðDH=DtÞj. However, after applying the level-set smoothing method, the interfaces are represented by the smoothed level-set function with significantly reduced noise. As a result, the magnitudes of ðHÞ and ðKÞ are reduced and, more importantly, the signs of the errors are more consistent between time-steps, which leads to error cancellation when computing DH=Dt and DK=Dt. Tables 1 and 2 summarize the statistics of the errors for all interfacial properties calculated from both /’s and uD ’s describing the spheres and the cylinders. The maximum error, maxðjjÞ, is the magnitude of the largest error on the surface mesh of the sphere. The standard deviation of the error based on area, r , is calculated on the triangulated surface mesh and is defined as

sP ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi N ðf f Þ2 Ai i¼1 numerical analytical PN

rðf Þ ¼

i¼1

Ai

jfanalytical j

;

ð18Þ

where the calculated interfacial property, fnumerical , is interpolated at the center of each triangulated element, N is the number of triangulated elements that make up the surface mesh of the spheres and cylinders, and Ai represent surface area of each triangulated element. The statistics of errors presented in the two tables emphasize the much improved accuracy of the calculation of interfacial properties after applying the level-set smoothing method, especially with the time derivatives of curvatures (c.f. Tables 1 and 2, r for DH=Dt and DK=Dt). As expected, the advective method yielded higher accuracy in the calculation of time derivatives of curvatures. Hence, we suggest the use of advective method for calculating the time derivatives of curvatures when meshing of the surface contour is possible. It should be noted that there exist first-order error associated with linear interpolation of interfacial properties to the location of the interfaces. However, we claim that these errors from interpolation are negligible based on the wide variation of maxðjjÞ and r values for different interfacial properties that were interpolated to the location of interfaces using the same linear interpolation scheme. 3.2. Validation of the calculation of time-derivatives of curvatures for simple geometries Our next step was to identify the effects of the size of the particles and the displacement of interfaces between time-steps on the magnitude of errors associated with the calculation of DH=Dt and DK=Dt. In this section, we considered only structures with spherical geometries. In order to conduct this parametric study, multiple sets of spherical particles with varying R1 values (10, 20, 30, 40 and 50) and R2 ¼ R1  jDRj values, where 1 6 jDRj  4, are created from the hyperbolic tangent function (see Section 2.2). After obtaining the smoothed level-set functions (by using the same smoothing sequence presented in Section 3.1), we calculated DH=Dt and DK=Dt using the advective method and analyzed the statistics of their respective errors. Fig. 4 shows the plots of maximum error, maxðjjÞ, and standard deviation of error based on area, r , associated with calculating DH=Dt and DK=Dt as functions of interfacial displacement, jDRj, for different R1 values. We find that jDRj ¼ 2 or 3 (interfacial displacement of about two or three grid points) is optimal for computing the time-derivative

52

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58

Fig. 3. The morphology of the spherical particle with R ¼ 40 described by (left column) the order parameter, /, before the smoothing process and (right column) the smoothed level-set function, uD , after the smoothing process, colored by the magnitude of errors of (a and b) H, (c and d) K, (e and f) v, (g and h) DH=Dt with the advective method and (i and j) DK=Dt also with the advective method.

quantities. The large magnitude of ðDH=DtÞ and ðDK=DtÞ for the smallest particle, R1 ¼ 10, results from the fact that as the radius of the particle approaches to the thickness of the diffuse interface, ðHÞ and ðKÞ values become large. For R1 ¼ 20 (five times the interfacial thickness), the errors were comparable to the cases with larger radius. Therefore, DH=Dt and DK=Dt calculated on interfaces with radii of curvatures greater than five times the interfacial thickness should yield sufficiently accurate results. As a note, further error analysis showed that H; K and v (normal velocity)

can be accurately calculated on interfaces with radii of curvatures greater than four times the interfacial thickness. 3.3. Comparison of interfacial locations before and after smoothing for a complex microstructure For complex morphologies with spatially varying curvatures, such as the Allen–Cahn structure, a direct comparison to known solutions cannot be made. Therefore, an alternative approach to

53

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58 Table 1 Error analysis of the interfacial properties calculated on two spherical particles with different radii. The variables maxðjjÞ and among all errors defined at vertices on the surface mesh and the standard deviation of the error based on area, respectively. Interfacial properties

H at radius = 40 K at radius = 40 H at radius = 38 K at radius = 38 v DH=Dt convective method DK=Dt convective method DH=Dt advective method DK=Dt advective method

r denote the magnitude of the maximum error After smoothing ðuD Þ

Before smoothing ð/Þ maxðjjÞ (%)

r (%)

maxðjjÞ (%)

r (%)

2.2 4.5 2.2 4.5 22 833 795 262 268

1.1 2.5 1.1 2.5 18 498 476 147 149

0.98 2.0 1.0 2.1 1.1 11 13 5.7 5.3

0.45 0.91 0.49 1.0 1.1 6.3 9.1 0.84 0.89

Table 2 Error analysis of the interfacial properties calculated on two cylindrical tubes with different radii. The variables maxðjjÞ and r denote the magnitude of the maximum error among all errors defined at vertices on the surface mesh and the standard deviation of the error based on area, respectively. Error analysis of K and DK=Dt calculations of the cylinder is omitted since their analytical values are zero. Interfacial properties

H at radius = 40 H at radius = 38 v DH=Dt convective method DH=Dt advective method

After smoothing ðuD Þ

Before smoothing ð/Þ maxðjjÞ (%)

r (%)

maxðjjÞ (%)

r (%)

2.7 2.6 20 1658 503

1.6 1.6 18 894 241

1.3 1.3 1.3 11 4.0

0.38 0.39 1.3 7.7 1.8

Fig. 4. The maximum error, maxðÞ, and the standard deviation of error, r , of DH=Dt (left column) and DK=Dt (right column) as a function of interfacial displacement, jDRj, for different particle sizes. The error plots show that jDRj values between two and three are optimal for computing time derivatives of curvatures.

validating a smoothing method is necessary. Here, we will examine the position of the interfaces before and after the smoothing process to determine how the level-set smoothing method alters the surface morphology.

Preservation of the interfacial location during a smoothing process is important, especially for complex morphologies that often have fine features as well as coarse ones. Even if the calculated curvature values appear to be smooth when plotted on the surface of a

54

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58

structure, if the interfacial shifts are significant during data processing, the calculated curvature values will not accurately represent the morphology and the topology of the original structure. In addition, if the shifts of the interfaces caused by a smoothing process is greater than the interfacial displacement due to microstructural evolution, the time-derivative calculations become meaningless. Therefore, an effective smoothing method must demonstrate that it can preserve the interfacial location sufficiently during the smoothing process. The level-set smoothing method was applied to convert the order parameter, /, from the phase-field simulation that describes the Allen–Cahn structure to the smoothed level-set function, uD . The same sequence of smoothing with the same set of parameter values as in Section 3.1 were employed. We then calculated the interfacial properties and plotted the values on the surface of the structure to ensure that the calculated quantities are smooth, as shown in Fig. 5. In order to visualize the effect of smoothing scheme on the interfacial locations, the contour lines (interfacial locations) obtained at / ¼ 0:5 and uD ¼ 0 are compared, as shown in Fig. 6. The overplots reveal a general overlap of the two contours of the interfaces from the pre-smoothed and the post-smoothed data, indicating that the level-set smoothing preserves the location of the interfaces for most interfaces. At surface patches with very high curvatures (where the radii of curvatures are comparable to the interfacial thickness), there are slight shifts in the two contour lines obtained from / and uD (see inset (a) for a magnified image of two regions with notable shifts). However, these regions are confined to small volumes of the entire structure and should not affect the overall results when statistically examined. In addition, these small shifts occur in regions where the interfacial thickness is comparable to the radius of curvature of the interface. Thus the phasefield calculation is under resolved at these points and these areas should not be included in the calculations.

3.4. Application of the method to a complex microstructure We have thus far demonstrated the level-set smoothing method’s capability to smooth the interfaces with minimal motion of the interfaces so that interfacial properties are accurately calculated. We also found that the errors associated with timederivative quantities are minimized when the interfacial displacement between time steps is about two to three grid points. Unlike the shrinking sphere, the Allen–Cahn structure possesses spatially dependent mean curvatures and time-derivatives, where some interfaces are evolving faster than others. Therefore, in order to maintain this range of interfacial displacement, it is necessary to take different Dt that depends on the normal velocity.

For the Allen–Cahn structure, the normal velocity, v, is proportional to the mean curvature [50],

v ¼ MH;

ð19Þ

where M is a proportionality constant as a function of L/ and 2 of the Allen–Cahn equation, Eq. (8). Therefore, to implement a multiple time-differential size to calculate DH=Dt and DK=Dt of the evolving Allen–Cahn structure, one only needs to consider the local mean-curvature values at the interface. First, we group the surface patches of the structure at simulation time t 0 based on their meancurvature values, and then we determine the required time-step, Dt, for each curvature group, such that all surface patches would be displaced two to three grid points except for the low curvature regions (as explained below). If we assume that the individual surface patches of a microstructure can be approximated as spherical patches with varying radii, then we can substitute v ¼ dR=dt and H ¼ 1=R in Eq. (19). After integration, one obtains

Dt ¼

Rðt0 þ DtÞ2  Rðt 0 Þ2 ; 2M

ð20Þ

where Rðt 0 Þ is the inverse of the mean curvature of a surface patch at time t 0 and Rðt 0 þ DtÞ is the corresponding value after Dt. Table 3 lists the curvature range, Rðt 0 Þ; Rðt0 þ DtÞ, and the corresponding Dt for each group of surface patches. For each curvature group, Rðt 0 Þ is the reciprocal of the jHj value at the lower end of the curvature range in each group. The time-step, Dt, in each group is determined by Eq. (20) setting Rðt0 þ DtÞ ¼ Rðt0 Þ  2Dx. Note that even though we assigned Dt values for interfaces with very large curvature values (in this case Rðt 0 Þ < 20 as discussed in Section 3.2), in practice the calculated values on these high curvature surfaces (which constitute only about 5% of the total interfacial area of the Allen–Cahn structure) contain significant errors and should still be discarded for quantitative analysis of the evolution of local curvatures. For the surface patches with jHj < 0:025, we use Dt ¼ 100. This is because, as Dt increases, the general morphology of the structure can undergo substantial evolution (the low curvature patches are connected with the rest of the structure) and consequently time derivatives of curvature calculations become meaningless. Therefore, one must strike a balance when choosing the appropriate time steps such that time-derivative computations are accurate and also physically meaningful. It should be noted that, because the interfacial displacements for the low curvature regions are less than the optimal range, the error in the calculated DH=Dt and DK=Dt at these interfaces maybe substantially larger than those discussed earlier. For example, maxðjjÞ and r in DH=Dt for a sphere with R1 ¼ 80 and R2 ¼ 79 are 78% and 7:3%, respectively. However, these large errors are not detrimental to the quantitative analysis

Fig. 5. Surface plots of the Allen–Cahn structure at t0 described by uD ¼ 0 and colored by (a) the mean curvature, (b) the Gaussian curvature, and (c) the normal velocity calculated from the time derivatives of the smoothed level-set functions.

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58

55

Fig. 6. Overplots of the interfacial contour lines of the Allen–Cahn structure described by the order parameter, /, (solid red line) and the smoothed level set function, uD (dotted black line). The contour lines are obtained on the x–y planes of the Allen–Cahn structure at different z values. Two contours are generally coincident except in the very high curvature regions where the radii of curvatures are comparable to the thickness of the diffuse interface, as shown in inset (a).

Table 3 Different evolution times, Dt, are required to ensure displacement of interfaces that yield accurate time derivatives of curvatures of the Allen–Cahn structure (Dx ¼ 1 is used). The left-hand column shows the range of jHj values that each group of surface patches have while Rðt0 Þ represents reciprocal of the H value at the lower end of the curvature range in each group and Rðt 0 þ DtÞ is the corresponding radius value after evolving for Dt. The values of Rðt0 Þ and Rðt0 þ DtÞ are omitted for jHj < 0:0125 because the minimum jHj value in this curvature range is zero. Range of jHj

Rðt0 Þ

Rðt 0 þ DtÞ

Dt

0:2 6 jHj 0:1 6 jHj < 0:2 0:05 6 jHj < 0:1 0:033 6 jHj < 0:05 0:025 6 jHj < 0:033 0:0125 6 jHj < 0:025 jHj < 0:0125

5.0 10.0 20.0 30.0 40.0 80.0 –

3.0 7.8 17.9 27.9 37.9 79.0 –

10.0 25.0 50.0 75.0 100.0 100.0 100.0

of local curvatures since these low curvature regions undergo very small change in curvatures; thus these surfaces do not contribute substantially to the overall microstructural evolution. Fig. 7 shows the calculated DH=Dt and DK=Dt of the Allen–Cahn structure described by the smoothed level-set function, uD , using a multiple time-differential size based on the Dt’s for each curvature

group shown in Table 3. To compare the values of DH=Dt and DK=Dt obtained from the convective and the advective method, the differences in their corresponding values were also plotted. Fig. 7c and f shows that the convective and the advective methods for calculating DH=Dt and DK=Dt agree for most part except in the high curvature regions, where the radii of curvatures are comparable to five times the interfacial thickness. Based on the error analysis performed in Sections 3.1 and 3.2, we have more confidence in the accuracy of the time derivatives of curvatures calculated with the advective method. The implications of the resulting data will be discussed in a future publication. In summary, the above analysis suggests the use of a multiple time-differential size for more accurate calculation of time derivatives of interfacial characteristics of a microstructure with complex morphologies. We note that, while the surface patches of the Allen–Cahn structure were approximated to be spherical for the above analysis, other simple geometries, such as cylinders, can be used for the same analysis. However, the resulting dependence of Dt on H remains the same regardless of the assumed geometry. Furthermore, high curvature regions typically possess geometry similar to a sphere, and therefore the results should be

56

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58

Fig. 7. Surface plots of the Allen–Cahn structure at t 0 described by uD ¼ 0 and colored by (a) DH=Dt using the convective method, (b) DH=Dt using the advective method, (c) the difference between (a) and (b), (d) DK=Dt using the convective method, (e) DK=Dt using the advective method, and (f) the difference between (d) and (e). All DH=Dt and DK=Dt calculations employed a multiple time-differential size, Dt, based on Table 3.

representative of rapidly evolving features in a morphologically complex microstructure.

4. Conclusion We presented an algorithm, ‘‘level-set smoothing,’’ that smoothes voxel-based data describing the interfaces with minimal shift in the interfaces so that interfacial properties can be accurately calculated. The significance of the level-set smoothing method is its ability to be applied to a wide range of complex microstructures obtained from 3-D experimental and computational techniques. The proposed method can be applied to any voxel-based data describing a 3-D two-phase structure so that the interfacial properties, including their rate of change, can be accurately calculated. These interfacial properties can then be used to characterize the morphology of a microstructure and its evolution, which is a crucial aspect of 3-D materials science. The level-set smoothing method is a set of sequential data-processing schemes that consists of first generating the signed distance function for the given microstructure using the level-set method, followed by smoothing via diffusion. We also presented numerical algorithms to calculate various types of interfacial properties, such as the curvatures and their time-derivatives, the latter of which can be calculated by either the convective or the advective methods. The error analysis of interfacial property calculations for simple geometries showed that the method improves the accuracy of the calculated curvatures and their time derivatives and that the magnitude of the errors are within the acceptable range for the analysis of these quantities. For example, in the case of spheres and cylinders with R1 ¼ 40 and R2 ¼ 38 described by the smoothed level-set functions (with Dx ¼ 1 and interfacial thickness of 4), the maximum error level was found to be 1% in H;  2% in K (for a sphere) and 4–7% in time derivatives of the curvatures when using the advective method. The improvement was especially significant for the time derivative quantities. We determined that an

interfacial thickness less than fifth of the radius of curvature ensures accurate values of DH=Dt and DK=Dt. We found that the advective method yielded more accurate results for both DH=Dt and DK=Dt than the convective method. The source of the increased error in the convective method was traced to the calculation of gradient and the partial time derivative, which are not required for the advective method. Furthermore, interfacial displacement of two to three grid points was shown to be optimal for calculating the time derivatives of curvatures. Comparison of the surface contours of the Allen–Cahn structure from the initial order parameter and the post-smoothed level-set function demonstrates that the level-set smoothing method sufficiently preserves the location of the interfaces except for where the curvature is large. Lastly, we have shown that if the evolution of the structure is non-uniform, multiple time differential size is suggested to accurately calculate the time derivatives of curvature. The considerations should be given to balance sufficient motion of interfaces for numerical accuracy while preventing major evolution of the microstructure. Acknowledgment Authors are grateful for the support from the Department of Energy Office of Basic Energy Science (Grant No. DE-FG0299ER45782/A012). Appendix A. Level-set formulation of curvatures and normal velocity The two common mathematical descriptions of the curvatures of the surfaces are the mean and the Gaussian curvatures, denoted by H and K, respectively. The analytical equation for the mean curvature in terms of the normal to the surface is

H¼

1 ðr  nÞ; 2

ðA:1Þ

57

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58

where n is the unit normal vector [43]. The Gaussian curvature can also be expressed in terms of the unit normal vector as

K ¼ n  adjðHeðuD ÞÞn;

ðA:2Þ

where the time derivative can be discretized as

@ uDijk ðt 1 Þ @t

D

where Heðu Þ is the 3 3 Hessian matrix of the second derivative of the smoothed level-set function, and adjðHeðuD ÞÞ is the adjoint of the Hessian matrix [43]. In the level-set approach, the unit normal vector in Eqs. (A.1) and (A.2) can be written in terms of the smoothed level-set function as

ruD : jruD j

n¼

1 2

2

2

3=2

2ðuDx þ uDy þ uDz Þ n 2 2 2 uDx ðuDyy þ uDzz Þ þ uDy ðuDxx þ uDzz Þ þ uDz ðuDxx þ uDyy Þ o  2ðuDx uDy uDxy þ uDy uDz uDyz þ uDx uDz uDxz Þ ; 1

K¼

uDx þ uDy 2 þ uDz 2

n

2

2

ðA:4Þ

2

2

2

2

2

2

uDx ðuDyy uDzz  uDyz Þ þ uDy ðuDxx uDzz  uDxz ÞuDz ðuDxx uDyy  uDxy Þ

h     þ 2 uDx uDy uDxz uDyz  uDxy uDzz þ uDy uDz uDxy uDxz  uDyz uDxx  io þ uDx uDz uDxy uDyz  uDxz uDyy ; ðA:5Þ where the subscripts of uD denote the spatial derivatives in their respective directions [43]. The first and second derivative of the smoothed level-set function shown in Eqs. (A.4) and (A.5) are calculated with the central differencing scheme. For example, the numerical stencils for the spatial derivatives along the x-direction is defined as

uDx ¼

uDiþ1;j;k  uDi1;j;k

uDxx ¼

; 2 Dx D u þ 16uiþ1;j;k  30uDijk þ 16uDi1;j;k  uDi2;j;k D iþ2;j;k

ð12DxÞ2

ðA:6Þ ;

ðA:7Þ

where i; j, and k denote the indices for x; y, and z positions, respectively. The first derivative is second-order accurate, and the second derivative is fourth-order accurate. The mixed partial derivative stencil is

uDxy ¼

1 ½uD  8uDiþ2;jþ1;k þ 8uDiþ2;j1;k 144DxDy iþ2;jþ2;k  uDiþ2;j2;k  8uDiþ1;jþ2;k þ 64uDiþ1;jþ1;k  64uDiþ1;j1;k þ 8uDiþ1;j2;k þ 8uDi1;jþ2;k  64uDi1;jþ1;k þ 64uDi1;j1;k  8uDi1;j2;k  uDi2;jþ2;k þ 8uDi2;jþ1;k  8uDi2;j1;k þ uDi2;j2;k ;

ðA:8Þ

which is also fourth-order accurate [52]. The use of higher-order schemes can reduce the grid anisotropy effect. On the other hand, the normal velocity, v, can be determined by [43]

v ¼

@ uD @t

uDijk ðt2 Þ  uDijk ðt1 Þ Dt

:

ðA:10Þ

Here, uD ðti Þ represents the smoothed level-set function for the microstructure at time t i and Dt ¼ t2  t 1 . The denominator of Eq. (A.9) is evaluated at t1 ; i.e., jruD ðt1 Þj. It is important to note that the discretization incurs a first order error, which is proportional to Dt.

ðA:3Þ

By substituting Eq. (A.3) into Eqs. (A.1) and (A.2), the mean and the Gaussian curvature can also be calculated by taking the spatial derivatives of uD ,

H¼

¼

jruD j;

ðA:9Þ

References [1] M.A. Meyers, K.K. Chawla, Mechanical Behavior of Materials, second ed., Cambridge University Press, Cambridge, UK, 2009. ISBN 978-0-521-86675-0. [2] J.R. Smith, A. Chen, D. Gostovic, D. Hickey, D. Kundinger, K.L. Duncan, R.T. DeHoff, K.S. Jones, E.D. Wachsman, Solid State Ionics 180 (1) (2009) 90–98, http://dx.doi.org/10.1016/j.ssi.2008.10.017. ISSN 0167-2738. [3] P. Tanasini, M. Cannarozzo, P. Costamagna, A. Faes, J. Van Herle, A. HesslerWyser, C. Comninellis, Fuel Cells 9 (5) (2009) 740–752, http://dx.doi.org/ 10.1002/fuce.200800192. ISSN 1615-6846. [4] H.-Y. Chen, H.-C. Yu, J.S. Cronin, J.R. Wilson, S.A. Barnett, K. Thornton, J. Power Sour. 196 (3, SI) (2011) 1333–1337, http://dx.doi.org/10.1016/ j.jpowsour.2010.08.010. ISSN 0378-7753. [5] R.E. Garcia, Y.-M. Chiang, W.C. Carter, P. Limthongkul, C.M. Bishop, J. Electrochem. Soc. 152 (1) (2005) A255–A263, http://dx.doi.org/10.1149/ 1.1836132. ISSN 0013-4651. [6] S.H. Choi, J.-W. Son, Y.S. Yoon, J. Kim, J. Power Sour. 158 (2, SI) (2006) 1419– 1424, http://dx.doi.org/10.1016/j.jpowsour.2005.10.076. ISSN 0378-7753. [7] E.A. Holm, P.M. Duxbury, Scr. Mater. 54 (6) (2006) 1035–1040, http:// dx.doi.org/10.1016/j.scriptamat.2005.11.048. ISSN 1359-6462. [8] G. Spanos, Scr. Mater. 55 (1) (2006) 3, http://dx.doi.org/10.1016/ j.scriptamat.2006.02.038. ISSN 1359-6462. [9] K. Thornton, H.F. Poulsen, MRS Bull. 33 (6) (2008) 587–595, http://dx.doi.org/ 10.1557/mrs2008.123. ISSN 0883-7694. [10] J. Alkemper, P.W. Voorhees, J. Microsc. 201 (3) (2001) 388–394, http:// dx.doi.org/10.1046/j.1365-2818.2001.00832.x. ISSN 0022-2720. [11] J. Alkemper, R. Mendoza, P.W. Voorhees, Adv. Eng. Mater. 4 (9) (2002) 694– 697. 10.1002/1527-2648(20020916)4:9<694::AID-ADEM694>3.0.CO;2-M. ISSN 1438-1656. [12] D.J. Rowenhorst, A. Gupta, C.R. Feng, G. Spanos, Scr. Mater. 55 (1) (2006) 11– 16, http://dx.doi.org/10.1016/j.scriptamat.2005.12.061. ISSN 1359-6462. [13] J.L. Fife, P.W. Voorhees, Acta Mater. 57 (8) (2009) 2418–2428, http:// dx.doi.org/10.1016/j.actamat.2009.01.036. ISSN 1359-6454. [14] N. Limodin, L. Salvo, M. Suery, M. DiMichiel, Acta Mater. 55 (9) (2007) 3177– 3191, http://dx.doi.org/10.1016/j.actamat.2007.01.027. ISSN 1359-6454. [15] W. Ludwig, S. Schmidt, E.M. Lauridsen, H.F. Poulsen, J. Appl. Crystallogr. 41 (2) (2008) 302–309, http://dx.doi.org/10.1107/S0021889808001684. ISSN 00218898. [16] J.L. Fife, P.W. Voorhees, Scr. Mater. 60 (10) (2009) 839–842, http://dx.doi.org/ 10.1016/j.scriptamat.2008.12.047. ISSN 1359-6462. [17] S. Terzi, L. Salvo, M. Suery, A. Dahle, E. Boller, Trans. Nonferrous Met. Soc. China 20 (3) (2010) S734–S738. ISSN 1003-6326. [18] U. Lienert, S.F. Li, C.M. Hefferan, J. Lind, R.M. Suter, J.V. Bernier, N.R. Barton, M.C. Brandes, M.J. Mills, M.P. Miller, B. Jakobsen, W. Pantleon, JOM 63 (7) (2011) 70–77, http://dx.doi.org/10.1007/s11837-011-0116-0. ISSN 10474838. [19] L.K. Aagesen, J.L. Fife, E.M. Lauridsen, P.W. Voorhees, Scr. Mater. 64 (5) (2011) 394–397, http://dx.doi.org/10.1016/j.scriptamat.2010.10.040. ISSN 13596462. [20] A.P. Mancuso, G.J. Williams, Nat. Photon. 6 (9) (2012) 574–575, http:// dx.doi.org/10.1038/nphoton.2012.209. ISSN 1749-4885. [21] J.R. Wilson, W. Kobsiriphat, R. Mendoza, H.-Y. Chen, J.M. Hiller, D.J. Miller, K. Thornton, P.W. Voorhees, S.B. Adler, S.A. Barnett, Nat. Mater. 5 (7) (2006) 541– 544, http://dx.doi.org/10.1038/nmat1668. ISSN 1476-1122. [22] S.J. Dillon, G.S. Rohrer, J. Am. Ceram. Soc. 92 (7) (2009) 1580–1585, http:// dx.doi.org/10.1111/j.1551-2916.2009.03064.x. ISSN 0002-7820. [23] J.R. Wilson, J.S. Cronin, S.A. Barnett, S.J. Harris, J. Power Sour. 196 (7, SI) (2011) 3443–3447, http://dx.doi.org/10.1016/j.jpowsour.2010.04.066. ISSN 03787753. [24] A. Karma, W.-J. Rappel, Phys. Rev. E 57 (4) (1998) 4323–4349, http:// dx.doi.org/10.1103/PhysRevE.57.4323. ISSN 1063-651X. [25] L.-Q. Chen, Annu. Rev. Mater. Res. 32 (2002) 113–140, http://dx.doi.org/ 10.1146/annurev.matsci.32.112001.132041. ISSN 1531-7331. [26] Y. Kwon, K. Thornton, P.W. Voorhees, Coarsening of bicontinuous structures via nonconserved and conserved dynamics, Phys. Rev. E 75 (2, 1) (2007), http://dx.doi.org/10.1103/PhysRevE.75.021120. ISSN 1539-3755. [27] Y. Kwon, K. Thornton, P.W. Voorhees, Philos. Mag. 90 (1–4) (2010) 317–335, http://dx.doi.org/10.1080/14786430903260701. ISSN 1478-6435. [28] G. Russo, P. Smereka, SIAM J. Sci. Comput. 21 (6) (2000) 2073–2095, http:// dx.doi.org/10.1137/S1064827599351921. ISSN 1064-8275.

58

C.-L. Park et al. / Computational Materials Science 85 (2014) 46–58

[29] Y. Xiang, L.-T. Cheng, D.J. Srolovitz, E. Weinan, Acta Mater. 51 (18) (2003) 5499–5518, http://dx.doi.org/10.1016/S1359-6454(03)00415-4. ISSN 13596454. [30] L. Tan, N. Zabaras, J. Comput. Phys. 221 (1) (2007) 9–40, http://dx.doi.org/ 10.1016/j.jcp.2006.06.003. ISSN 0021-9991. [31] R. Mendoza, I. Savin, K. Thornton, P.W. Voorhees, Nat. Mater. 3 (6) (2004) 385– 388, http://dx.doi.org/10.1038/nmat1138. ISSN 1476-1122. [32] A.L. Genau, P.W. Voorhees, Acta Mater. 57 (20) (2009) 6226–6233, http:// dx.doi.org/10.1016/j.actamat.2009.08.049. ISSN 1359-6454. [33] R. Mendoza, K. Thornton, I. Savin, P.W. Voorhees, Acta Mater. 54 (3) (2006) 743–750, http://dx.doi.org/10.1016/j.actamat.2005.10.010. ISSN 1359-6454. [34] V.W.L. Chan, K. Thornton, Acta Mater. 60 (6–7) (2012) 2509–2517, http:// dx.doi.org/10.1016/j.actamat.2011.12.042. ISSN 1359-6454. [35] M. Zhihong, C. Guo, M. Yanzhao, K. Lee, Comput.-Aided Des. 43 (12) (2011) 1561–1566, http://dx.doi.org/10.1016/j.cad.2011.06.006. ISSN 0010-4485. [36] B. Hamann, Curvature approximation for triangulated surfaces, in: G. Farin, H. Noltemeier, H. Hagen, W. Knödel (Eds.), Geometric Modelling, Computing Supplementum, vol. 8, Springer Vienna, Vienna, Austria, 1993, pp. 139–153. [37] J. Goldfeather, V. Interrante, ACM Trans. Graph. 23 (1) (2004) 45–63, http:// dx.doi.org/10.1145/966131.966134. ISSN 0730-0301. [38] E. Magid, O. Soldea, E. Rivlin, Comput. Vis. Image Understand. 107 (3) (2007) 139–159, http://dx.doi.org/10.1016/j.cviu.2006.09.007. ISSN 1077-3142. [39] D.S. Meek, D.J. Walton, Comput. Aided Geomet. Des. 17 (6) (2000) 521–543, http://dx.doi.org/10.1016/S0167-8396(00)00006-6. ISSN 0167-8396. [40] N. Meyer, M. Desbrun, P. Schroder, A.H. Barr, Discrete differential-geometry operators for triangulated 2-manifolds, in: H.C. Hege, K. Polthier (Eds.), Visualization and Mathematics III, Springer-Verlag Berlin, Berlin, Germany, 2003, pp. 35–57. ISBN 3-540-01295-8.

[41] X. Chen, F. Schmitt, Intrinsic surface properties from surface triangulation, in: G. Sandini (Ed.), Computer Vision — ECCV’92, Lecture Notes in Computer Science, vol. 588, Springer Berlin Heidelberg, Berlin, Germany, 1992, pp. 739– 743. ISBN 978-3-540-55426-2ISSN 0302-9743. [42] K. Watanabe, A.G. Belyaev, Comput. Graph. Forum 20 (3) (2001) C385–C392, http://dx.doi.org/10.1111/1467-8659.00531. ISSN 0167-7055. [43] J. Sethian, Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science, second ed., Cambridge University Press, Cambridge, UK, 1999. [44] S. Osher, J.A. Sethian, J. Comput. Phys. 79 (1) (1988) 12–49, http://dx.doi.org/ 10.1016/0021-9991(88)90002-2. ISSN 0021-9991. [45] J.W. Bullard, E.J. Garboczi, W.C. Carter, E.R. Fuller Jr., Comput. Mater. Sci. 4 (2) (1995) 103–116, http://dx.doi.org/10.1016/0927-0256(95)00014-H. ISSN 0927-0256. [46] O.I. Frette, G. Virnovsky, D. Silin, Comput. Mater. Sci. 44 (3) (2009) 867–875, http://dx.doi.org/10.1016/j.commatsci.2008.06.006. ISSN 0927-0256. [47] Y.-T. Kim, N. Goldenfeld, J. Dantzig, Phys. Rev. E 62 (2, B) (2000) 2471–2474, http://dx.doi.org/10.1103/PhysRevE.62.2471. ISSN 1063-651X. [48] K.A. Smith, F.J. Solis, L. Tao, K. Thornton, M.O. de la Cruz, Phys. Rev. Lett. 84 (1) (2000) 91–94, http://dx.doi.org/10.1103/PhysRevLett.84.91. ISSN 0031-9007. [49] M. Sussman, P. Smereka, S. Osher, J. Comput. Phys. 114 (1) (1994) 146–159, http://dx.doi.org/10.1006/jcph.1994.1155. ISSN 0021-9991. [50] S.M. Allen, J.W. Cahn, Acta Mater. 27 (6) (1979) 1085–1095, http://dx.doi.org/ 10.1016/0001-6160(79)90196-2. ISSN 0001-6160. [51] J.W. Cahn, J.E. Hilliard, J. Chem. Phys. 28 (2) (1958) 258–267, http://dx.doi.org/ 10.1063/1.1744102. ISSN 0021-9606. [52] L.K. Aagesen, Phase-Field Simulation of Solidification and Coarsening in Dendritic Microstructures, Ph.D. Thesis, Northwestern University, 2010.