Efficient thermal field computation in phase-field models Jing-Rebecca Li

∗†

Donna Calhoun

‡

Lucien Brush

§

February 2, 2008

Abstract We use a fast solver for the heat equation in free-space to compute the thermal field in a simulation of crystal growth at low undercoolings. The solver is based on the efficient direct evaluation of the integral representation of the solution to the constant coefficient, free space heat equation with a smooth source term. The computational cost and memory requirements of the new solver are reasonable compared to standard methods and allow one to solve for the thermal field in a computational domain whose size depends only on the size of the growing crystal and not on the extent of the thermal field. This can result in significant computational savings in situations where the crystals grow slowly relative to the expansion of the thermal field and may make parameter regimes typically encountered in experimental situations more accessible to numerical treatment. We demonstrate the method by solving the phase-field equations for anisotropic, dendritic growth at low undercooling in two dimensions. Key words. phase-field; crystal growth; dendritic solidification; diffusion equation; unbounded domain; integral representation; fast solvers

1

Introduction

Freely growing single crystals exhibit a wide range of crystal-melt interface morphologies, including faceted and dendritic structures, and as a result, have been of great interest to the engineering, scientific, and mathematical communities for decades. Researchers have attempted to understand, predict, and control solidification morphologies which ultimately impact material properties and performance[27]. Crystal growth morphologies have also been a subject of interest to researchers studying nonlinear dynamical systems because freely growing crystals are examples of rich, pattern-forming systems [11, 29]. Numerical simulation of dendritic growth has become an important tool in understanding how crystal morphologies evolve in various experimental and theoretical regimes. It has attracted the attention of computational scientists motivated to develop and improve numerical methods for accurately reproducing the two- and three-dimensional tree-like structures that emerge when crystals form in an undercooled melt. As a result, several effective methods are now available, including those based on solving the physically derived phase-field equations [6, 15, 35, 44], as well as level set and Lagrangian methods which can be used to solve the related Stefan problem [7, 32, 41]. For an overview of phase-field modeling and comtemporary views of crystal growth, see the review papers by Karma [25], Provatas et. al.[39], Boettinger et. al. [3] and Glicksman [17]. Regimes of theoretical and experimental importance most easily accessible by numerical computations are those for which the temperature of the liquid melt is relatively far from the ∗ Corresponding

author. Tel. : +33 1 39 63 53 55; E-mail address : [email protected] Projet POEMS, Domaine de Voluceau - Rocquencourt, 78153 Le Chesnay Cedex, France ‡ DEN/DM2S/SFME/LTMF, Commissariat a l’Energie Atomique, F-91191 Gif-sur-Yvette Cedex, France § Department of Materials Science and Engineering, University of Washington Seattle, WA, 98195 † INRIA-Rocquencourt,

1

equilibrium melting temperature. In such regimes, the non-dimensional undercooling parameter, the Stefan number, is large (typically O(1)) and crystals grow fast relative to the expansion of the thermal field. One can observe interesting results before the extent of the thermal field exceeds the computational domain and adequate resolution of the solid-liquid interface can be achieved at modest computational expense. Many of the early dendritic growth calculations reported in the literature were computed in this regime [4, 34, 42, 43]. In another regime of theoretical importance, one assumes an infinite thermal decay rate and models the heat transport using a free-space Laplace equation. In this case, the computational mesh for the thermal field can be eliminated entirely by using front-tracking methods coupled with boundary integral techniques to directly solve for growth velocities along the solid-liquid interface[5, 8, 9]. By constrast, obtaining accurate numerical solutions in the regime of low undercoolings still proves to be computationally challenging. For low Stefan numbers (see Eqn. (8)) crystals grow very slowly relative to the rate of expansion of the thermal front and the extent of the thermal field in the liquid can quickly exceed by several orders of magnitude the size of the growing crystal. To obtain results which can be used to confirm existing theory or be validated by experiments, long simulation times may be required. In the free growth experiments of Koss et al. carried out in NASA’s USMP-4 (reference [18] in Provatas et. al. [38]), pivalic acid was grown at a variety of undercoolings. The lowest undercooling corresponded to a melt temperature of approximately 0.5K below the melting point, or to a Stefan number of about 0.05. Under these conditions, it can take a long time to establish steady state dendrite tip growth conditions, which has often been the goal of many experimentalists trying to verify theories of dendritic growth. In addition, the length scale representative of the solid-liquid interface thickness is on the order of a few nanometers at most, for the primary, secondary and tertiary dendritic branches it can be up to microns and for the extent of the thermal field, centimeters. A numerical simulation which attempts to reproduce these experimental results would require a computational domain large enough to capture the entire thermal field over long simulation times, while adequately resolving the interface thickness as well as the emerging dendritic features. The need for large computational domains arises from the fact that in order to match experiments, numerical simulations must not exhibit any spurious artifacts introduced by the presence of an artifical computational boundary. During the early stages of actual crystal growth, the experimental container is typically large relative to the size of the crystal. For an initially uniform undercooled melt, disturbances in the thermal field due to the growth of a crystal are damped at distances as far as the container boundary, and the thermal field is fully captured by the container. The experimental container is assumed to have negligible, if any, effect on the evolving solidification front. In fact, this is a critical requirement for experimentalists attempting to validate the steady state theories of dendrite tip velocities and tip curvatures. For slowly growing crystals, long simulations times are required before interesting results can be observed, and as a result, the thermal field will exhibit time dependent behavior far from the cystal interface. Capturing the thermal field numerically places significant demands on computational resources.

1.1

Numerical approaches

In this paper, we focus on the efficient computation of the thermal field during free crystal growth in the low undercooling parameter regime. Unlike typical finite difference, finite volume, or even finite element approaches, our approach does not require artificial boundary conditions or excessively large computational domains. Because we base our approach on efficient, direct evaluation of the analytic solution to the heat equation in the infinite domain, we obtain the correct far-field behavior on reasonable domains and at reasonable computational cost. Before describing our method in detail, however, we describe some other approaches to efficiently simulating crystal growth processes. Many approaches to modeling crystal growth involve finite-difference or finite element schemes. These are easy to implement, and when used to solve the phase-field equations for modeling crystal growth, do not require any additional front tracking for the liquid-melt interface. A general problem encountered by such methods, however, is that artificial boundary conditions must

2

be imposed on the thermal field at the computational domain boundary. One hopes that the solution thus obtained mimics the evolution of the transport of heat in the free space. However, for all but the earliest stages of growth, boundary conditions that are most easily, and hence most often, applied, e.g. purely Dirichlet or Neumann (in [4, 34, 40]), are inappropriate and will result in a distorted thermal field and an incorrect solution to the free space crystal growth problem. To reduce the adverse effects caused by imposing inappropriate boundary conditions on finite computational boundaries, one can simply increase the size of the computational domain. However, this only becomes a feasible approach when one includes mesh adaptivity. Provatas et al.[37] have recently developed an adaptive mesh refinement (AMR) technique to successfully model two-dimensional free growth at the low undercoolings described in the pivalic acid experiment above. However, their method and many others like it still require significant investment in developer time and computational resources to gain the full advantage of such an approach [4, 40, 42]. Moreover, parallelizing AMR codes, especially those involving implicit solvers, is considered an active area of research. Another approach, but one which does not appear to be widely used, involves formulating the exact transparent boundary conditions for the heat equation and then coupling these conditions to solvers for the interior of the domain. The idea is that if heat is allowed to escape the computational domain as if the transparent boundary weren’t there, then the domain need only be large enough to contain the crystal interface itself. Such transparant boundary conditions for the heat equation are described in [12, 22–24]. However, we are not aware of any attempts to use these boundary conditions in the problem of modeling crystal growth. Moreover, such boundary conditions typically can only be applied on domain boundaries of particular shape, e.g. circles for two-dimensions and spheres for three dimensions. This may force the use of grids based on polar or spherical coordinates, or unstructured grids, which in turn may require specialized solvers which tend to be less accurate as those available for Cartesian meshes. To our knowledge, artificial boundary conditions for more general domain shapes and which are computationally efficient are not yet available. In our approach, we eliminate the need for a computational boundary for the thermal field by discretizing the integral representation of the solution of the free space heat equation. We will apply a recent numerical method which solves the free space heat equation with a smooth source term, the Fast Recursive Marching (FRM) method [31], which computes these integrals in an efficient way, to solve the thermal field equation in the phase-field model. This will be coupled to a standard finite-difference method for solving the equation for the phase whose level contours describes the interface between the liquid and solid. The main advantage of this approach is that no artificial boundary conditions are required for the thermal field on the computational boundary, and an accurate solution to the evolving thermal field can be obtained on a computational domain whose size is dictated by the size of the crystal, and not by the thermal field. The computational time and memory requirements of the Fast Recursive Marching method are on the same order as for a standard explicit finite difference/element method with the same spatial resolution. The gain comes from the fact that the new approach allows smaller computational domains and hence fewer spatial discretization points in total. The savings increases with increasing simulation time since the diffusion of the √ thermal field occurs on a length scale that is O( t), where t is the time of simulation. We summarize here the most important features of the Fast Recursive Marching method. Details are fully described in [31]. • The numerical solution is based on the efficient evaluation of the integral representation of the solution to the constant coefficient, time dependent heat equation. Hence, no artificial boundary conditions are needed. • The problem of the “burden of history” typically associated with evaluating convolution integrals in time is eliminated. At each new time step, the solution at the current step can be computed using only values from a few previous time steps. Also, there are no numerical stability constraints.

3

• Temporal accuracy depends on the polynomial approximation order for the source term. • The method is spectrally accurate in space. In this paper, we present the method for uniform Cartesian meshes and explain where appropriate how to extend it to locally refined as well as unstructured meshes. The numerical results presented in this paper are for a dimensionless undercooling of 0.1. The refinement of the spatial mesh near the solidification front will need to be included to treat the case of very low undercoolings (e.g. S ∼ 0.001). Incorporating time and spatial adaptivity into the FRM method poses no theoretical difficulties, but at this time, only a simple adaptive enlargement of the computational domain as the crystal approaches the computational boundary has been implemented. The thermal field can be computed on the enlarged computational domain in a totally natural way even though the support of the thermal field has long since expanded outside of the previous computational domain. Such an approach would not be feasible with typical finite-difference based methods. This paper is a follow-up of the work presented in [31]. In that paper, a simplified phase-field model (no anisotropy, no noise) was included in the numerical results section. Here, we include the anisotropic model with noise and more detailed numerical results. More importantly, we discuss how to apply the FRM method in an appropriate way to the solution of the phase-field model. We explain how to choose the various parameters in the FRM method as well as show that one can easily incorporate adaptive enlargement of the computational domain into the phase-field computation. This paper is organized as follows. In Section 2, we briefly describe the phase-field model and our discretization of the coupled thermal field and phase equations. In Section 3 we describe the Fast Recursive Marching method for solving the free space heat equation and explain how it should be applied to solve the thermal field equation. In Section 4 we show numerical results for a case of two-dimensional, anisotropic crystal growth at a dimensionless undercooling of 0.1. Section 5 contains the conclusions.

2

The phase-field model for dendritic crystal growth

We will consider the solidification of a pure material. The model we will use are the phasefield equation, which consist of a coupled system of two diffusion-type equations governing a dimensionless temperature u(x, y, t) (the thermal field) and the phase variable φ(x, y, t).

2.1

The model

We use the following formulation of the model equations in two dimensions, taken from [34]: 30φ2 (1 − φ)2 φt , x = (x, y) ∈ R2 , t > 0, S   ¯2 1 φt − ¯2 ∇2a (θ)(φ) = φ(1 − φ) φ − [1 + AN (x, y, t)] m ¯ 2 ut − ∇2 u = −

2

2

+ 30φ (1 − φ) ¯ α S u,

(1)

(2)

2

x = (x, y) ∈ R , t > 0,

subject to initial conditions: u(x, y, 0) = u0 (x, y), φ(x, y, 0) = φ0 (x, y),

x = (x, y) ∈ R2 , 2

x = (x, y) ∈ R .

(3) (4)

Here, u = (T − TM )/∆T ), TM is the equilibrium melting temperature at a planar interface and ∆T is the melt under-cooling (the difference between the melting temperature and the initial temperature in the bulk liquid). Interfacial properties such as the crystal-melt interface energy and the local atomic attachment coefficient depend on interfacial orientation, i.e. they are anisotropic. In (1-4) the anisotropy of the material properties is included through the

4

anisotropic Laplacian, which depends on the angle θ (for example, the angle that the normal to the iso-φ contours makes with the x-axis):    
∂ ∂ ∂ ∂ η(θ)η 0 (θ) + η(θ)η 0 (θ) + ∇ · η 2 (θ)∇ , (5) ∇2a (θ) = − ∂x ∂y ∂y ∂x with θ computed by the formula θ = tan−1



∂y φ ∂x φ

 .

(6)

η(θ) = 1 + γ cos w θ,

(7)

The anisotropy is typically modeled by

where γ is the strength of the anisotropy and w is the crystalline symmetry (e.g., w = 4 for four-fold symmetry). The phase variable varies from φ = 0 to φ = 1, corresponding to the solid and liquid phases, respectively. The value of φ varies smoothly but sharply from 0 to 1 in a narrow diffuse interface region of width O(¯ ). Parameters α and m ¯ are related to specific material properties and the length scale. The function N (x, y, t) is a random noise function, added to promote dendrite growth and the parameter A is the amplitude of the noise. The Stefan number, or dimensionless undercooling, is given by c∆T , (8) L where c is the specific heat per unit volume and L is the latent heat per unit volume. The undercooling is the parameter of primary importance in this paper. Typical simulations reported in the literature assume S to be in the range 0.05 ∼ 1. However, experimental results often report values of S in the range 0.001 ∼ 0.1 [28]. Early work on numerical simulation using phase-field models can be found in [26, 33, 34, 43, 45]. Since we will be using our Fast Recursive Marching method for the thermal field, we do not require any numerical boundary conditions on the thermal field. For the phase-field variable, we use standard finite-difference methods, and so boundary conditions on the computational boundary will be required. However, the phase-field variable φ can be considered constant in the bulk liquid away from the interface, and so, assuming the computational domain is large enough to contain the support of the solid-liquid interface, one can safely impose conditions compatible with φ ≡ 1 on computational boundary (e.g. homogeneous Neumann conditions). Although the equation for φ also appears to be a diffusion-type equation, in fact, it is not a true diffusion equation. This can be seen by writing the term ∇a (θ)∇φ in tensor notation as ∇L(θ)∇φ and noting that the resulting anisotropic tensor L(θ) is skew symmetric, with complex eigenvalues, and so does not behave as a true diffusion tensor [43]. As a consequence, the field for φ does not extend beyond the interface and is constant outside of this region. The field of interest in the phase-field equations is represented by the phase variable. The thermal field is of interest near the dendritic tips but is of only secondary interest away from the crystal interface. For this reason, one would ideally like to use a computational domain for u which is the same size as needed for φ even at low undercoolings. We want the thermal field to diffuse out through the computational boundary as if the boundary were not there. S=

2.2

Discretization of coupled system

We use the following first order time discretization of the model equations:   ¯2 (φm+1 − φm ) − ¯2 ∇2 (φm+1 ) = ¯2 ∇2a (θ)(φm ) − ∇2 (φm ) m ¯ ∆t   1 + φm (1 − φm ) φm − [1 + AN (x, y, t)] 2 + 30 φ2m (1 − φm )2 ¯ α S um ,

5

(9)

Obtain um+1 by solving : ut − ∇2 u = approx to

  30φ2 (1 − φ)2 − φt ≡ f , S

(10)

where the subscript m is used to denote the approximate value of the variable at t = m∆t. The computational domain for both (9) and (10) at t = m∆t will be Ωm , which contains the solid-liquid interface, and may be adaptively enlarged as simulation continues. We treat the coupling of equations (1- 2) by evaluating the right hand side of each equation using previously computed values of the variables. This operator-splitting approach is formally first order accurate, although in practice one can often obtain accuracy rates higher than this. For equation (9) the contribution of u to the right hand side is approximated by its value at t = m∆t. Hence, we also use a first order time discretization of φt . Overall, the discretization of (2) is semi-implicit: the isotropic heat part of (2) is done implicitly, the rest (the difference between the isotropic and the anisotropic Laplacians as well as the nonlinear part) is done explicitly. The boundary conditions on ∂Ωm are homogeneous Neumann. Since φ ≡ 1 in the bulk liquid, these boundary conditions are acceptable. Spatial discretization of (9) can be done in standard ways, for example, via finite differences. We solve (10) using the FRM method, the details of which will be given in Section 3. There is no need for boundary conditions on ∂Ωm because the correct far field behavior of u is automatically satisfied. The right hand side f of (10) will be treated as known, using already computed values of φ, via the following first order in time approximation:
2 30φ2m+1 (x) 1 − φm+1 (x) (φm+1 (x) − φm (x)) fm+1 (x) ≈ − . (11) S ∆t Higher order methods may be formulated with (9) and (10) as the basic components but will not be discussed here because the emphasis of this paper is on the correct treatment of the thermal field equation while avoiding overly large computational domains. The current implementation is done for two dimensions but there is a natural extension to three dimensions.

3

Fast Recursive Marching method

In this section we outline the Fast Recursive Marching (FRM) method and explain how to apply it to solve (10). We will only point out the main features of the FRM method, details can be found in [31]. The FRM method solves the heat equation in free space with a source term: u ˜t (x, t) − ∇2 u ˜(x, t) = f (x, t),

x ∈ Rd , t > 0,

(12)

supplemented by the initial condition u ˜(x, 0) = u ˜0 (x),

x ∈ Rd .

(13)

The functions f (x, t) and u ˜0 (x) are assumed to be compactly supported for any fixed t. Because of the choice of the normalization in the phase-field model, we have u = 1 as the base line value for the thermal field in the bulk liquid in the initial configuration. To obtain the baseline value of 0 and compact support for the initial configuration, we have defined the new variable u ˜ = u − 1. We will concentrate on the case d = 2 but will use the notation d if the statement applies to other dimensions as well. 2 2 φt , and that we will use the apWe recall the right hand side of (10) is f = − 30φ (1−φ) S proximation of it given in (11). Clearly, the support of f (·, t) is on a narrow band around the crystal-melt interface (φ = 0 in solid, φ = 1 in liquid) at time t. Given ∆t, we will compute the approximate solution at t = ∆t, 2∆t, · · · , m∆t, · · · , using an increasing set of computational domains, Ω1 ⊂ · · · ⊂ Ωm · · · , where each Ωi contains the support of f (·, t), for all t ≤ i∆t. The domain will be enlarged adaptively after the support of the phase variable φ, which is computed via (9), is determined. In addition, we will assume that all the computational domains contain the support of u0 .

6

The FRM method is applicable to problems where the source term f is k-times differentiable. Let tmax be an estimate of the longest simulation time of interest, the actual simulation time Q may be much shorter, and let B = di=1 [−Ri /2, Ri /2] be a rectangle which contains the support of f (·, t), for all t ≤ tmax . The values tmax , max{R1 , · · · , Rd } will be used to choose certain parameters in the FRM method. We assume that f (x, t) ∈ C k (B × [0, tmax ]). This is satisfied because φ is smooth in the diffuse interface region. The fact that φ varies sharply in this narrow region will affect the parameters we choose in the algorithm and we will address this when we come to the relevant parameters. ¿From standard potential theory [21], the solution to (12-13) can be written as Z u ˜(x, t) = k (x − y, t) u ˜0 (y) dy + Rd

Zt Z k (x − y, t − τ ) f (y, τ ) dy dτ,

(14)

0 Rd

where the fundamental solution for the heat equation in Rd is 2

k (x, t) :=

e−kxk

/4t d

, t > 0.

(15)

(4πt) 2

We will refer to the first integral in (14) as an initial potential and the second integral as a volume potential. Straightforward discretization of the above integrals leads to an prohibitively expensive numerical scheme - the solution is dependent on the full space-time history of the diffusion process. If M is the number of time steps, and the computational domain is not adaptively enlarged during these M time steps and contains N spatial discretization points, it is easy to see that O(N 2 M + N 2 M 2 ) work is required for those M time steps with a naive discretization. The Fast Recursive Marching method is the accelerated computation of these integrals. There are three components to the algorithm; • We accelerate the time convolution by doing the time stepping in the Fourier domain. • To compute the inverse Fourier transform, we use an efficient quadrature of the inverse Fourier integral that involves many fewer Fourier nodes than a naive trapezoidal rule. • These Fourier nodes are not equi-spaced, so to transform to and from Fourier space in an efficient manner, we use the Non-uniform FFT (NUFFT) [18, 30]. In this way, we are able to reduce the total computational cost of evaluating (14) to O(M · (Nf +Np ) log(Nf +Np )), where Np is the number of discretization points in the physical domain, and Nf is the number of nodes in the Fourier domain (which depends on the desired tolerance and an estimate of the smoothness of f ). In our numerical experiments, Nf is typically on the same order as Np . For the same Np , this algorithmic complexity is comparable to finite difference or finite elements methods. The advantage of this approach is that many fewer points in the physical domain will be needed for long time simulation of slow solidification because of the smaller computational domain.

3.1

Fourier representation of the solution

We propose time stepping using the Fourier transform of u ˜: Z u ˆ(s, t) = eis·x u ˜(x, t) dx, Rd

7

(16)

where the inverse Fourier transform is given by Z 1 u ˜(x, t) = e−is·x u ˆ(s, t) ds. (2π)d

(17)

Rd

It is straightforward to show that u ˆ(s, t) satisfies the ordinary differential equation dˆ u (s, t) = −ksk2 u ˆ(s, t) + fˆ(s, t), dt where fˆ(s, t) =

Z

eis·x f (x, t) dx.

(18)

(19)

Rd

Thus, time stepping for u ˆ(s, t) can take the form 2

u ˆ(s, t) = e−ksk where

∆t

u ˆ(s, t − ∆t) + Φ(s, t, ∆t),

Zt

2

e−ksk

Φ(s, t, ∆t) =

(t−τ )

fˆ(s, τ ) dτ .

(20)

(21)

t−∆t

Thus, in the Fourier domain, time convolution can be accelerated from quadratic complexity to linear complexity: u ˆ(s, t) is simply damped and updated by a local (in time) computation at each time step, as indicated in (20). We will refer to Φ(s, t, ∆t) in (21) as the update integral. The observation in (20) is not new, but for it to be useful, one must be able to choose a small set of Fourier modes which can adequately represent the function u ˜ and one must be able to move easily between Fourier and physical domains.

3.2

Choice of Fourier nodes

The choice of Fourier nodes s depends on our choice of quadrature for evaluating Z 1 e−is·x u ˆ(s, t) ds. u ˜(x, t) = (2π)d

(22)

Rd

An optimal choice for the Fourier nodes is one that results from a quadrature for (22) which meets the desired accuracy requirements with as few nodes as possible. Let us recall that the source is smooth: f (·, t) ∈ C k (B). Thus, fˆ(s, t) = O(ksk−k ) for large s. A straightforward calculation shows that if  is the error tolerance, then the evaluation of 1 fˆ(s, t) in (19) needs to be done only for ksk ≤ P , where P = O(1/) k . We will refer to P as the high-frequency cutoff. For the phase-field application, f is a function of φ. Since we will be using methods which are of low order accuracy in time, we choose an error tolerance of  = 10−4 because we do not expect to get more than 4 digits of accuracy. We make an estimate for the smoothness of φ. 1 1 For example (1/) 2 = 100, (1/) 3 ≈ 22. For our simulations we typically set P = 50 but this can be modified if it is found to be inadequate. To evaluate the finite Fourier integral Z 1 u ˜(x, t) ≈ e−is·x u ˆ(s, t) ds, (23) (2π)d ksk 0 be the desired precision, let Lmin = − log(1/) and let Lmax = dlog P e, where P is the high-frequency cutoff. Further, let {sj,1 , . . . , sj,n(j) } and {wj,1 , . . . , wj,n(j) } be the nodes and weights for the n(j)-point Gauss-Legendre quadrature on the interval [2j , 2j+1 ], where n(j) = max(Lx 2j+3/2 , 8 log(1/)). Then, in one space dimension, Z LX max n(j) X isj,k x −isj,k x −is·x (24) (e +e )u ˆ(sj,k , t) wj,k = O() u ˆ(s, t)ds − e j=Lmin k=1 |s| 0, which is a significantly larger number of nodes, especially when  is small. The parameter Lx is the largest magnitude of the evaluation points x. Since we have supQ posed that the largest computational domain needed will be contained in B = di=1 [−Ri /2, Ri /2], we have the estimate Lx ≤ max{R1 /2, · · · , Rd /2}. In the code we have written we in fact use a more efficient quadrature than indicated by Theorem 1 because we take into account tf , the longest needed simulation time. We also verify the error bound in (24) numerically and adaptively choose the smallest n(j), the number of nodes in [2j , 2j+1 ], which has a quadrature error smaller than 3(2j+1 −2j ) . These modifications reduce the number of s nodes required for a given error tolerance .

3.3

The forward transform

To compute fˆ(s, t) from the data f (x, t), since we have assumed that the data is supported in the rectangle B, we need to compute the finite integral: Z R1 /2 Z Rd /2 fˆ(s, t) = ··· eis·x f (x, t) dx, (25) −R1 /2

−Rd /2

where x = (x1 , ..., xd ), s = (s1 , ..., sd ). 1 Let us recall that P = O(1/) k is the high-frequency cutoff. This bounds the oscillatory is·x behavior of the term e in the integrand of (25). Together with the fact that f (x, t) is smooth, it follows that the trapezoidal rule applied to each dimension in (25) with N points will yield O(N −k ) accuracy [10]. In dimension i, the error will begin decaying rapidly once Ni is of the order O(P Ri ), meaning that the integrand is well-resolved. If f (x, t) is given on a uniform mesh with Ni points in each dimension, the trapezoidal rule yields fˆ(s, t) ≈

 N1 d  Y Ri X i=1

Ni

n1 =1

···

Nd X

eis·xn f (xn , t),

(26)

nd =1


where xn = − R1 /2 + n1 (R1 /N1 ), · · · , −Rd /2 + nd (Rd /Nd ) . We have used the uniform trapezoidal rule in (26) to obtain our numerical results. But we emphasize here that it is desirable to use a locally refined spatial discretization of (25)

9

instead, putting a finer discretization of x points near the crystal-melt interface and a coarser discretization away from it. Given a code which computes φ on an adaptive mesh, for example on locally refined rectangles or triangles, it is easy to generate the quadrature weights to get a modified rule for the Fourier transform in (25). The only restriction is that the density of the discretization in physical space still needs to satisfy ∆xi ≤ O( P 1Ri ) globally to resolve the oscillations in eis·x . Spatial refinement is then used to resolveQ the other term, f (x, t). We denote the number of points x in physical space by Np (with Np = di=1 Ni if refinement is not used).

3.4

The nonuniform FFT

Let us denote by s1 , · · · , sQ and w1 , · · · , wQ the full set of one dimensional quadrature nodes and weights described in the preceding section. Then, we need to compute the discrete transform u ˜(x, t) ≈

Q Q X 1 X · · · e−isj ·x u ˆ(sj , t) wj1 · · · wjd , (2π)d j =1 j =1 1

(27)

d

where sj = (sj1 , . . . , sjd ). The summation in (27) must be carried out for each evaluation point x. Direct evaluation at each of the Np points x defined in Section 3.3 would require O(Nf · Np ) work. Likewise, the sum (26) must be evaluated at each of the Nf spectral nodes, also requiring O(Nf · Np ) work. Since the {sj } are not uniformly spaced, the classical Fast Fourier Transform (FFT) cannot be applied to accelerate this computation. In the last decade, however, suitable fast algorithms have been developed that we refer to as Non-Uniform Fast Fourier Transforms (NUFFTs). Detailed descriptions and additional references can be found in [2, 13, 14, 16, 18, 30, 36]. Like the FFT, these algorithms evaluate sums of the form (26) and (27) in O((Nf + Np ) log (Nf + Np )) work, for any fixed precision. When the x points are uniformly spaced, the sums in (26) and (27) are computed by what are called type-1 and type-2 NUFFTs [30]. When spatial refinement in physical space is used, the target and source points in (26) and (27) are both non-uniformly spaced, and the sums are computed by what is called a type-3 NUFFT.

3.5

Computing the update integral

The update integral (21) can be computed to high order accuracy using a standard product integration approach [10]. We first approximate fˆ(s, t) by a polynomial in time, and then integrate the resulting polynomial analytically. Linear approximation of fˆ(s, t) yields second order accuracy and the rule [20] Φ(s, t, ∆t) = W0 (s, ∆t)fˆ(s, t) + W1 (s, ∆t)fˆ(s, t − ∆t), −z

W0 (s, ∆t) =

e

−1+z ∆t, z2

−z

W1 (s, ∆t) =

1−e

− ze z2

(28)

−z

∆t ,

(29)

where z = ksk2 ∆t. For a cubic approximation of fˆ(s, τ ), yielding fourth order accuracy in time, we have Φ(s, t, ∆t) = W0 (s, ∆t)fˆ(s, t) + W1 (s, ∆t)fˆ(s, t − ∆t) + W2 (s, ∆t)fˆ(s, t − 2∆t) + W3 (s, ∆t)fˆ(s, t − 3∆t),

10

(30)

where (6 + 2z 2 − 6z)e−z + (6z 3 − 11z 2 − 6 + 12z) ∆t, 6z 4 (−z 2 + 2z 3 + 6 − 4z)e−z + (−6z 2 − 6 + 10z) ∆t, W1 (s, ∆t) = −2z 4 (2z + 2z 2 − 6)e−z + (−8z + 3z 2 + 6) W2 (s, ∆t) = ∆t, −2z 4 (z 2 − 6)e−z + (6 + 2z 2 − 6z) W3 (s, ∆t) = ∆t. 6z 4 Some care needs to be taken in computing the weights above. If z is small, the formulae above are subject to significant cancellation error. In that case, the weights can be computed by Taylor expansion of the exponentials, carried out to a sufficient number of terms. A reasonable compromise is to switch from the analytic expression to the Taylor series if z < 10−3 . Four terms in the Taylor series are then sufficient for twelve digit accuracy. W0 (s, ∆t) =

3.6

Adaptive enlargement of the computational domain

The computational domain for u ˜ must contain the crystal-melt interface. As this interface grows in size, we can enlarge the computational domain. Suppose at t = m∆t, the computational domain for both u ˜ and φ is Ωm and at the next time step, we would like to enlarge it to Ωm+1 ⊃ Ωm . Since we compute the Fourier coefficients of u ˜ at t = (m + 1)∆t, which are the values u ˆ(s, (m + 1)∆t), to obtain u ˜(x, (m + 1)∆t) on the enlarged domain Ωm+1 , we simply perform the inverse Fourier transform in (27) at the set of points x ∈ Ωm+1 instead of only in Ωm . All the information about u ˜ (hence u) is stored in the Fourier coefficients and no information is lost even though the thermal field has expanded far outside of Ωm . This is in contrast to a finite difference/element method where even if exact artificial boundary conditions were imposed on ∂Ωm , it is not possible to accurately reconstruct the solution u outside of the domain Ωm .

3.7

The numerical scheme

With the full complement of tools in place, we can now provide an informal description of the Fast Recursive marching (FRM) method. For lth order accuracy in time, we refer to the method as FRM(l). The method for a second order accurate scheme is described in Figure 1. The modifications necessary for the fourth order accurate FRM(4) method are straightforward. The total computational cost is of the order O(M · (Nf + Np ) log(Nf + Np )), where Nf is the number of Fourier nodes and Np is the number of x points in the physical domain. Two NUFFT calculations are required per time step. This is true for both FRM(2) and FRM(4), so that their execution times are nearly identical. FRM(4) does, however, require twice as much storage since it is a four-level marching scheme.

4

Numerical results

We show numerical results for the phase-field model equations shown in (1-2). The problem parameters we have chosen are m ¯ = 0.035, ¯ = 0.005, α = 400, A = 0.025, with anisotropy parameters w = 4, γ = 0.015, and the undercooling parameter S = 0.1, which is at the low end of the undercoolings which have been simulated in literature [4, 34]. To compute at even lower undercoolings, local spatial refinement should be used in complement to our approach and this is currently under study.

11

The FRM(2) Method Step 1: Initialization a) Select the time step ∆t. b) Select the error tolerance  for the quadrature of the inverse Fourier integral. c) Obtain Q quadrature nodes and weights in Fourier space. d) Compute weights W0 , W1 from (29). Step 2: Transform initial data a) Initialize u ˆ(s, 0) at all Nf = Qd spectral nodes by computing the forward NUFFT of u ˜0 (x). ˆ b) Compute f (s, 0) at all Nf spectral nodes by computing the forward NUFFT of f (x, 0). Step 3: For m = 1, . . . a) Compute fˆ(s, m∆t) at all Nf spectral nodes by computing the forward NUFFT of f (x, m∆t) using the data on Ωm−1 . b) Update u ˆ(s, m∆t) according to (20): u ˆ(s, m∆t) = e−ksk

2

∆t

u ˆ(s, (m − 1)∆t) + W0 (s, ∆t)fˆ(s, m∆t) + W1 (s, ∆t)fˆ(s, (m − 1)∆t)

c) Compute u ˜(x, m∆t), x ∈ Ωm , from u ˆ(s, m∆t) using the (inverse) NUFFT. The domain Ωm may be larger than Ωm−1 .

Figure 1: The Fast Recursive Marching algorithm FRM(2) for solving the heat equation in the infinite domain. We use the time discretization of the coupled thermal field and phase equations in (9-10). The thermal field equation in (10) is solved using the FRM(2) method, the source term being the first order approximation in (11). The cost of solving (10) is that of performing two NUFFTs. The discretization of the phase equation in (9) will be denoted SIE(1) in the figures which follow, which stands for first order semi-implicit Euler. We discretize the spatial derivatives in (9) on a uniform grid and solve the resulting Helmholtz equation using the software package FISHPACK [1], which utilizes cyclic reduction for separable PDEs. Because we have a uniform grid, the cost of solving (9) is that of performing one uniform FFT. Numerical results are shown in Figure 2 with the time step ∆t = 0.002, and ∆x = ∆y = 0.01. In Figure 2, we show the contour plots of u and φ at various times. The background color (or grayscale) corresponds to the field values of u on the computational domain. The black lines are contour lines for φ, corresponding to the values 0.01, 0.1, 0.5, 0.99, respectively. These lines indicate the limits of the crystal-melt interface. In Figure 2(a) we show the initial seed configuration. Up to t = 0.5 we used the computational domain [−0.5, 0.5]2 , which is marked as a black frame in Figures 2(b) and 2(c). At t = 0.5, φ is approaching the boundary of the computational domain, so we enlarge Ωm to [−1, 1]2 . Since all the information about u(·, m∆t) is stored in the Fourier coefficients, u ˆ(·, m∆t), enlarging Ωm just involves computing the inverse Fourier transform at a larger set of x points at the current time step. This is in contrast to finite difference/elements methods where one must enlarge Ωm as soon as the support of um approaches ∂Ωm ; in our approach, we do not have to enlarge Ωm

12

until the support of φm approaches ∂Ωm , a great advantage at low undercoolings. In Figures 2(d) and 2(e) we show the simulation up to t = 2, on the computational domain [−1, 1]2 (marked as a black frame). At t = 2 we enlarge Ωm again, to [−2, 2]2 , which we use up to t = 12. The last picture shown in Figure 2(f) is for t = 25 computed on a grid of [−3, 3]2 . It is easy to see that when the thermal field is computed with the Fast Recursive Marching method, heat diffuses out of the computational domain correctly and simulation can continue until the solidification front reaches the computational boundary. On the largest computational domain [−3, 3]2 , the side length Lx = 3, and with our choice of the high frequency cutoff P = 50, we have Q = 120, meaning there are Nf = 2 ∗ Q2 = 28, 800 nodes in the Fourier domain. By contrast, there are Np = 600 ∗ 600 = 360, 000 nodes in the physical domain. In this example, both FRM(2) (for u) and SIE(1) (for φ) use the same number of (uniformly spaced) discretization points in physical space. Due to the fact that one FFT is computed in SIE(1) and two NUFFTs are computed in FRM(2), and typically, a type 1 or type 2 NUFFT takes three times as long as one uniform FFT, we expect that for the same Np , one iteration of FRM(2) will take six times as long as one iteration of SIE(1). Indeed, this is what we observe numerically. For this example, the computational time of FRM(2) is between 4 to 8 times that of SIE(1). We emphasize that for the case of low undercoolings, this increase is small compared to the additional work that would have been required if we were forced to compute u on a much larger computational domain due to the√wrong artificial boundary conditions. In that case, the domain would have had side length O( t), where t is the simulation time, which translates to √d Np = O( t ) discretization points.

5

Conclusions

We used the Fast Recursive Marching method, a Fourier-based method for the solution of the heat equation in free space with a smooth source term, to compute the thermal field in the phasefield model of dendritic solidification. In the case of low undercoolings, this algorithm facilitates efficient and accurate long-time simulations because it allows the use of smaller computational domains for the thermal field. Future work needs to be done to incorporate spatial adaptivity into the code, so that the mesh is locally refined to resolve the micro-structures that emerge at the solid-liquid interface. The FRM method can also be used on more general types meshes: logically Cartesian, triangular, or unstructured meshes. There is also a natural extension to three dimensions. While we have demonstrated our solver on the phase-field equations, it can also be incorporated into solvers for the sharp interface model. This would involve solving an initial boundary value heat equation instead of the free space heat equation with a source term.

References [1] J. Adams, P. Swarztrauber, and R. Sweet. FISHPACK: A Package of FORTRAN Subprograms for the Solution of Separable Elliptic Partial Differential Equations. The National Center for Atmospheric Research, Boulder, Colorado (80307) U.S.A., version 3.2 edition, November 1988. [2] G. Beylkin. On the fast fourier transform of functions with singularities. Appl. Comput. Harmonic Anal., 2:363–383, 1995. [3] W. Boettinger, J. A. Warren, and C. Beckermann. Phase-field simulation of solidification. Ann. Rev. Mat. Res, 32:163–194, 2002. [4] R. J. Braun and B. T. Murray. Adaptive phase-field computations of dendritic crystal growth. J. Crystal Growth, 174(1-4):41–53, April 1997.

13

FRM(2)/SIE(1), S = 0.1, time = 0

FRM(2)/SIE(1), S = 0.1, time = 0.25

2

0 −0.1

1.5

FRM(2)/SIE(1), S = 0.1, time = 0.5

2

0 −0.1

1.5

−0.2 1

−0.4

0

−0.5 −0.6

−0.5

−0.2 1

−0.3 0.5

−0.4

0

−0.6

−0.5 −1

−0.9

u: colorbar φ: contours at {0.01, 0.1, 0.5, 0.99} −1

0

1

2

−1

−0.4

0

−0.5 −0.6

−0.5

−0.7

−0.7 −1

−0.8 −1.5

−0.3 0.5

−0.5

−0.7 −1

−0.1

1.5

1

0.5

0

−0.2

−0.3

−2 −2

2

−0.8 −1.5 −2 −2

−0.9

u: colorbar φ: contours at {0.01, 0.1, 0.5, 0.99} −1

(a) t=0

0

1

2

−1

−0.8 −1.5 −2 −2

−1

(b) t=0.25

0

0

1

2

2

FRM(2)/SIE(1), S = 0.1, time = 25

0 3

−0.1

1.5

1.5

−0.2

2

−0.2

1 −0.3

−0.3 0.5

−0.4

0

−0.5 −0.6

−0.5

1

0.5

−0.4

0

−0.5

−0.4 0

−0.6

−0.5

−0.7

−0.7

−0.6 −1

−1

−1

−0.8

−0.8 −1.5 −2 −2

0

−0.1

−0.2 1

−1

(c) t=0.5

FRM(2)/SIE(1), S = 0.1, time = 2

FRM(2)/SIE(1), S = 0.1, time = 1.5 2

−0.9

u: colorbar φ: contours at {0.01, 0.1, 0.5, 0.99}

−0.9

u: colorbar φ: contours at {0.01, 0.1, 0.5, 0.99} −1

0

(d) t=1.5

1

2

−1

−1.5 −2 −2

−1

0

1

(e) t=2

−0.8

−2

−0.9

u: colorbar φ: contours at {0.01, 0.1, 0.5, 0.99} 2

−1

u: colorbar φ: contours at {0.8, 0.85, 0.9, 0.99} −3 −3

−2

−1

0

1

2

3

−1

(f) t=25

Figure 2: The thermal field is computed via the second order Fast Recursive Marching method (FRM(2)). The phase is computed via first order semi-implicit Euler (SIE(1)). When the thermal field is computed with the Fast Recursive Marching method heat diffuses out of the computational domain correctly and simulation can continue until the solidification front reaches the computational boundary. [5] L. N. Brush and R. F. Sekerka. A numerical study of two-dimensional crystal-growth forms in the presence of anisotropic growth-kinetics. J. Crystal Growth, 96(2):419–441, 1989. [6] G. Caginalp and P. Fife. Phase-field methods for interfacial boundaries. Phys. Rev. B, pages 7792–7794, 1986. [7] S. Chen, B. Merriman, S. Osher, and P. Smereka. A simple level set method for solving Stefan problems. J. Comput. Phys., 135(1):8–29, 1997. [8] V. Cristini and J. Lowengrubb. Three-dimensional crystal growth I: linear analysis and self-similar evolution. J. Crystal Growth, 240:267–276, 2002. [9] V. Cristini and J. Lowengrubb. Three-dimensional crystal growth II: nonlinear simulation and control of the Mullins Sekerka instability. J. Crystal Growth, 266:552–567, 2004. [10] P. J. Davis and P. Rabinowitz. Methods of numerical integration. Academic Press, San Diego, 1984.

14

[11] S. H. Davis. Theory of Solidification. Cambridge University Press, Cambridge, England, 2001. Cambridge Monographs on Mechanics. [12] E. Dubach. Artificial boundary conditions for diffusion equations: numerical study. J. Comput. Appl. Math., 70(1):127–144, 1996. [13] A. Dutt and V. Rokhlin. Fast fourier transforms for nonequispaced data. SIAM J. Sci. Comput., 14:1368–1383, 1993. [14] J. A. Fessler and B. P. Sutton. Nonuniform fast fourier transforms using min-max interpolation. IEEE Trans. Signal Proc., 51:560–574, 2003. [15] G. Fix. Phase field models for free boundary problems. In A. Fasano and M. Primicerio, editors, Free Boundary Problems, Theory and Application, pages 580–589, 1983. [16] K. Fourmont. Non-equispaced fast fourier transforms with applications to tomography. IEEE Trans. Signal Proc., 9:431–450, 2003. [17] M. E. Glicksman and A. O. Lupulescu. Dendritic crystal growth in pure materials. J. Crystal Growth, 264:541–549, 2004. [18] L. Greengard and J. Y. Lee. Accelerating the nonuniform fast fourier transform. SIAM Rev., 46:443–454, 2004. [19] L. Greengard and P. Lin. Spectral approximation of the free-space heat kernel. Appl. Comput. Harmon. Anal., 9(1):83–97, 2000. [20] L. Greengard and J. Strain. A fast algorithm for the evaluation of heat potentials. Comm. Pure Appl. Math., 43(8):949–963, 1990. [21] R. B. Guenther and J. W. Lee. Partial differential equations of mathematical physics and integral equations. Prentice Hall, Inglewood Cliffs, New Jersey, 1988. [22] L. Halpern and J. Rauch. Absorbing boundary conditions for diffusion equations. Numer. Math., 71(2):185–224, 1995. [23] H. Han and Z. Huang. Exact and approximating boundary conditions for the parabolic problems on unbounded domains. Comput. Math. Appl., 44(5-6):655–666, 2002. [24] H. Houde and Z. Huang. A class of artificial boundary conditions for heat equation in unbounded domains. Comput. Math. Appl., 43(6-7):889–900, 2002. [25] A. Karma. Phase field modeling. In S. Yip, editor, Handbook of Materials Modeling, Volume I: Methods and Models, pages 2087–2103. Springer, Netherlands, 2005. [26] A. Karma and W. J. Rappel. Phase-field method for computationally efficient modelling of solidification with arbitrary interface kinetics. Phys. Rev. E, pages R3017–R3020, 1996. [27] W. Kurz. Fundamentals of Solidification. Trans Tech, Aedermannsdorf, Switzerland, 1992. [28] J. C. LaCombe, M. B. Koss, D. C. Corrigan, A. O. Lupulescu, L. A. Tennenhouse, and M.E. Glicksman. Implications of the interface shape on steady-state dendritic crystal growth. J. Crystal Growth, 206:331–344, 1999. [29] J. S. Langer. Instabilities and pattern formation in crystal growth. Rev. Mod. Phys., 52:1–28, 1980. [30] J. Y. Lee and L. Greengard. The type 3 nonuniform FFT and its applications. J. Comput. Phys., 206(1):1–5, 2005. [31] J. R. Li and L. Greengard. On the numerical solution of the heat equation i: Fast solvers in free space. J. Comput. Phys., 226:1891–1901, 2007. [32] Z. Li and B. Soni. Fast and Accurate Numerical Approaches for Stefan Problems and Crystal Growth. Numerical Heat Transfer, Part B: Fundamentals, 35:461–484, 1999. [33] B. T. Murray, W. J. Boettinger, G. B. McFadden, and A. A. Wheeler. Computation of dendritic solidification using a phase-field model. ASME Heat Transfer Melting, Solidification Cryst. Growth, pages 67–76, 1993.

15

[34] B. T. Murray, A. A. Wheeler, and M. E. Glicksman. Simulations of experimentally observed dendritic growth behavior using a phase-field model. J. Crystal Growth, 154:386–400, 1995. [35] O. Penrose and P. Fife. Thermodynamically consistent models of phase-field type for the kinetics of phase transitions. Physica D, pages 44–62, 1990. [36] D. Potts, G. Steidl, and M. Tasche. Fast fourier transforms for nonequispaced data: A tutorial. In J.J. Benedetto and P. Ferreira, editors, Modern Sampling Theory: Mathematics and Applications, pages 249–274. Birkhauser, Boston, 2001. [37] N. Provatas, N. Goldenfeld, and J. Dantzig. Efficient computation of dendritic microstructures using adaptive mesh refinement. Physical Review Letters, 1998. [38] N. Provatas, N. Goldenfeld, J. Dantzig, J. C. LaCombe, A. Lupulescu, M. Koss, M. Glicksman, and R. Almgren. Cross-over scaling in dendrtic evolution at low undercooling. Physical Review Letters, 1999. [39] N. Provatas, M. Greenwood, B. Athreya, N. Goldenfeld, and J. Dantzig. Multiscale modeling of solidification: Phase-field methods to adaptive mesh refinement. Intern. J. Mod. Phys. B, 19(31):4525–4565, 2005. [40] J. Rosam, P. K. Jimack, and A. Mullis. A fully implicit, fully adaptive time and space discretisation method for phase-field simulation of binary alloy solidification. J. Comput. Phys., 225(2):1271–1287, 2007. [41] J. A. Sethian and J. Strain. Crystal growth and dendritic solidification. J. Comput. Phys., 98(2):231–253, 1992. [42] R. Tonhardt and G. Amberg. Phase-field simulation of dendritic growth in a shear flow. J. Crystal Growth, 194:406–425, 1998. [43] S. L. Wang and R. F. Sekerka. Algorithms for Phase Field Computation of the Dendritic Operating State at Large Supercoolings. J. Comput. Phys., 127:110–117, 1996. [44] S. L. Wang, R. F. Sekerka, A. A. Wheeler, B. T. Murray, S. R. Coriell, R. J. Braun, and G. B. McFadden. Thermodynamically-consistent phase-field models for solidification. Physica D, pages 189–200, 1993. [45] A. A. Wheeler, B. T. Murray, and R. J. Schaefer. Computation of dendrites using a phase field model. Physica D, pages 243–262, 1993.

16