A fast, robust, and accurate operator splitting method for phase-field simulations of crystal growth Yibao Li, Hyun Geun Lee, Junseok Kim à Department of Mathematics, Korea University, Seoul 136-701, Republic of Korea article info Article history: Received 30 August 2010 Received in revised form 27 October 2010 Accepted 25 February 2011 Communicated by J.J. Derby Available online 4 March 2011 Keywords: A2. Crystal growth A1. Phase-field simulation A1. Operator splitting A1. Multigrid method abstract In this paper we propose a fast, robust, and accurate operator splitting method for phase-field simulations of dendritic growth in both two- and three-dimensional space. The proposed method is based on operator splitting techniques. We split the governing phase-field equation into three parts: the first equation is calculated by using an explicit Euler’s method. The second is a heat equation with a source term and is solved by a fast solver such as a multigrid method. The third is a nonlinear equation and is evaluated using a closed form solution. We also present a set of representative numerical experiments for crystal growth simulation to demonstrate the accuracy and efficiency of the proposed method. Our simulation results are also consistent with previous numerical experiments. & 2011 Elsevier B.V. All rights reserved. 1. Introduction Crystal growth is a classical example of phase transformations from the liquid phase to the solid phase via heat transfer. For several decades, to understand and simulate crystal growth, several meth- ods have been developed including boundary integral [1–4], cellular automaton [5–8], front-tracking [9–13], level-set [14–17], Monte- Carlo [18,19], and phase-field [20–37] methods. Among these various methods, the phase-field method is popular and widely used. Its advantage is that the explicit tracking of the interface is unnecessary by introducing an order parameter, i.e., a phase-field variable. In this paper, we focus on the phase-field method for crystal growth problems which avoids difficulties associated with tracking the interface and computes complex crystal shapes. We consider here the solidification of a pure substance from its supercooled melt in both two- and three-dimensional space. A great challenge in the simulation with various supercoolings is the large difference in time and length scales. In order to overcome this, many numerical methods have been proposed such as explicit [11,27, 31,36,38], mixed implicit-explicit [30,35,37], and adaptive [21,22,29, 32,33] methods. In explicit methods, which are widely used, the solutions become unstable for large time steps. For this reason the authors in [11,36] suggested Dt oh 2 =ð4DÞ for stability of explicit methods. Here, Dt is the time step, h is the mesh size, and D is the thermal diffusivity. In [11], the time step is also restricted to Dt rh=ð10jV max jÞ, where jV max j is the magnitude of the maximum value of the interface velocity. Also the authors in [36] showed that Dt ¼ h 2 =ð5D L Þ works well for not too large choices of the anisotropy by numerical experiments, where D L ¼ M f e 2 , M f is the kinetic mobility, and e is the interface energy anisotropy. Implicit methods allow relatively larger time steps, however, they are computation- ally more expensive per step than explicit ones. Another classical method is a multiple time-step algorithm that uses a larger time step for the flow-field calculations while reserving a finer time step for the phase-field evolution [34]. The use of mesh adaptivity is a natural choice to overcome this problem. However, explicit adaptive technology also suffers the time step restriction. Therefore, we need a scheme that allows the use of a sufficiently large time step without the technical limitations. In this paper we present a new, computationally efficient, and robust operator splitting algorithm for solving phase-field simulations of dendritic growth and demon- strate the accuracy and efficiency of the method by a set of representative numerical experiments. This paper is organized as follows: in Section 2 the governing equations for crystal growth based on the phase-field method are given. In Section 3 we describe the computationally efficient operator splitting algorithm. In Section 4 we present numerical results for solving the crystal growth simulation both in 2D and 3D. Finally, conclusions are given in Section 5. 2. The phase-field model The basic equations of the phase-field model are derived from a single Lyapounov functional [39]. We model the solidification in Contents lists available at ScienceDirect journal homepage: www.elsevier.com/locate/jcrysgro Journal of Crystal Growth 0022-0248/$ - see front matter & 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.jcrysgro.2011.02.042 à Corresponding author. Tel.: þ82 2 3290 3077; fax: þ82 2 929 8562. E-mail address: cfdkim@korea.ac.kr (J. Kim). URL: http://math.korea.ac.kr/~cfdkim/ (J. Kim). Journal of Crystal Growth 321 (2011) 176–182