Research paper
An accurate and parallel method with post-processing boundedness control for solving the anisotropic phase-field dendritic crystal growth model

https://doi.org/10.1016/j.cnsns.2022.106717Get rights and content

Highlights

  • Direction splitting method for simulating anisotropic dendritic crystal growth.

  • Stabilizing and post-processing for the boundedness control of numerical solution.

  • The method has the advantages of accuracy, stability, and speed.

  • Theoretical analyses and numerical simulations are carried out.

Abstract

A fast, accurate, and stable numerical algorithm is proposed to solve the anisotropic phase-field dendritic crystal growth model. The algorithm combines the first-order direction splitting method and the linear stabilization technique, which preserve the energy stability while ensuring fast computing. The related stability analysis is shown. To obtain high accurate solution, the solution of first-order direction splitting method is extrapolated to be fourth-order, and the fourth-order difference method is applied to spatial discretization. To further improve the stability and prevent the blow-up of numerical solution, two post-processing methods, the bound preserving least-distance modification and the cut-off approach, are developed to control the boundedness of numerical solution, and their theoretical proofs are given. The proposed novel algorithm can be performed in parallel in each time step, significantly improving computing efficiency. Several numerical examples demonstrate the effectiveness of the proposed algorithm, including convergence and stability tests, the effectiveness of two post-processing techniques, and a series of 2D and 3D dendritic crystal growth simulations.

Introduction

Snowflake patterns, frost patterns, metal solidification, crystals generated in supersaturated fluids are examples of dendritic crystal growth [1], [2], [3]. The formation of dendritic crystals and the evolution of the structure are critical aspects of the solidification process. By controlling the dendritic crystal growth, the organization can be modified to optimize the final properties of the cured product. For decades, to understand and simulate crystal growth, several methods have been developed, including integral boundary [4], [5], cellular automata [6], front tracking [7], [8], level set [9], [10], Monte-Carlo space [11] and phase-field [12], [13], [14], [15] methods. Because of its simplicity and versatility, the phase-field approach has been frequently employed in this area for mathematical modeling and numerical simulation of the dendritic crystal growth.

The present paper aims to design an efficient numerical algorithm for the anisotropic phase-field dendritic crystal growth model proposed in [16]. The model is a second-order system of equations coupled by Allen–Cahn type equations with gradient-dependent anisotropy coefficients and heat conduction equations. In recent years, the efficient numerical algorithm of the isotropic Allen–Cahn type equations have been extensively developed. Also, significant progress has been made in dealing with the Allen–Cahn equation’s difficulties, which includes stiff, high computational complexity and energy stability, etc. [17], [18], [19], [20]. Solving anisotropically coupled systems faces more challenges than the isotropic case: (1) how to deal with nonlinear gradient anisotropy coefficients; (2) how to preserve computational stability; (3) how to realize fast computing and high precision simultaneously. Therefore, there have been a limited number of successful attempts to build efficient algorithms.

Now we briefly introduce the anisotropic phase-field dendritic crystal growth model proposed in [16] by energy variational approach and present some notations which will be used in the later analysis. The total free energy takes the form: E(ϕ,u)=Ω12|κ(ϕ)ϕ|2+F(ϕ)4ϵ2+λ2ϵKu2dx,where xΩRd,d=2,3, ϕ(x,t) is the density field where ϕ=1 for the solid and ϕ=1 for the fluid. These two values are connected by a smooth transitional layer with the thickness ϵ. Among them, ϵ, λ and K are all positive parameters, u(x,t) is the temperature, F(ϕ)=ϕ212 is the double well potential. κ(ϕ) is a function describing the anisotropic property that depends on the direction of the outer normal vector n which is the interface normal defined as n=ϕ|ϕ|. For the 2D system, the anisotropy coefficient κ(ϕ) is usually given by κ(ϕ)=1+ϵ4cos(mΘ),where m is the number of folds of anisotropy, ϵ4 is the parameter for the anisotropy strength, and Θ=arctan ϕyϕx. When m=4 (i.e., fourfold anisotropy), for instance, κ(ϕ) can be easily reformulated in terms of the phase-field variable ϕ, namely, for 2D κ(ϕ)=13ϵ41+4ϵ413ϵ4ϕx4+ϕy4|ϕ|4,and for 3D, κ(ϕ)=13ϵ41+4ϵ413ϵ4ϕx4+ϕy4+ϕz4|ϕ|4.By adopting the Allen–Cahn type (L2gradient flow) relaxation dynamics for the dendritic crystal growth, one obtains the governing dynamical equations via the variational approach, which reads as follows: τ(ϕ)ϕt=δEδϕλϵp(ϕ)u=κ2(ϕ)ϕ+κ(ϕ)|ϕ|2H(ϕ)f(ϕ)ϵ2λϵp(ϕ)u,ut=DΔu+Kp(ϕ)ϕt,where H(ϕ) is the variational derivative of κ(ϕ). In 2D, it reads as H(ϕ)=δκ(ϕ)δϕ=4ɛ44|ϕ|6ϕxϕx2ϕy2ϕy4,ϕyϕx2ϕy2ϕx4,whereas in 3D, H(ϕ)=δκ(ϕ)δϕ=4ɛ44|ϕ|6ϕxϕx2ϕy2+ϕx2ϕz2ϕy4ϕz4,ϕyϕy2ϕz2+ϕx2ϕy2ϕx4ϕz4,ϕzϕx2ϕz2+ϕy2ϕz2ϕx4ϕy4.In the system Eq. (1.5), τ(ϕ)>0 is the mobility constant that is chosen either as a constant or a function of ϕ [16], [21], constant D is the diffusion rate of the temperature, δEδϕ is the variational derivative of the total energy with respect to ϕ, the function p(ϕ) accounts for the generation of latent heat and it is a phenomenological functional taking the form preserving the minima of ϕ at ±1 independently of the local value of u. For p(ϕ), we choose p(ϕ)=15ϕ523ϕ3+ϕ and p(ϕ)=1ϕ22. The periodic boundary condition is considered in this study. The model equations Eq. (1.5) follow the dissipative energy law. By taking the L2 inner product of the first equation of Eq. (1.5) with ϕt, and of the second equation of Eq. (1.5) with λϵKu, using the integration by parts and combining the obtained two equalities, we obtain ddtE(ϕ,u)=τ(ϕ)ϕt2λDϵKu20.

In recent years, many numerical methods have been studied. In the explicit scheme [22], [23], all non-linear terms are explicitly processed. Although it is easy to construct, it has a strict step size limit, and the solution becomes unstable under a large time step. For this reason, the authors in [21] suggested δt<min{h2/4D,ϵ2} for stability of explicit methods. The stabilized semi-implicit scheme [24], [25], [26], [27] adds appropriate stability items to ensure unconditional energy stability. An artificial dissipation term is added to the semi-implicit scheme, resulting in an error term. The convex splitting approach [18], [28], [29], [30], [31], [32] has been widely used and analyzed in isotropic and anisotropic models. The energy potential’s convex part is treated implicitly, whereas the concave part is processed explicitly. Although it can be demonstrated to be energy stable, it frequently results in a non-linear system, making implementation difficult and computationally expensive. However, whether the product of the gradient potential and the anisotropy coefficient can be separated into convex and concave components for some dendritic models is still debatable. Recently, Yang studied the invariant energy quantization (IEQ) and scalar auxiliary variable (SAV) method to solve the model Eq. (1.5) [33], [34], [35], [36], [37]. By incorporating appropriate auxiliary variables, the IEQ and SAV techniques write the energy function into a modified quadratic form. As a result, a new system of equations can be generated that is equivalent to the original equation. In the new system, the equation terms are handled implicitly or explicitly. The resulting scheme has a high time order, an energy that is unconditionally stable, and an extraordinarily high efficiency. The operator splitting approach [19], [38], [39], [40], [41], [42], [43], [44], [45] decouples the model equations’ different physical processes (diffusion and reaction) to produce multiple subproblems. Then solve the sub-problems using appropriate numerical methods, which can be linear and quick. Nonetheless, proving its discrete energy stability is challenging. In [38], the author uses the operator splitting method to split the dendritic crystal growth model according to different physical processes, which greatly reduces the difficulty in computing. However, the computational complexity is still high since no direction splitting is performed.

In this paper, our research focuses on efficient and accurate computation. Direct spatial discretization in high-dimensional space will result in high computing costs and storage space. Therefore, we employ the first-order sequential operator splitting approach to decompose high-dimensional issues into a sequence of one-dimensional subproblems that may be performed in parallel, avoiding complex and time-consuming calculations. In order to avoid the more significant influence of anisotropic coefficients and nonlinear terms, when solving the subproblems of each direction, two stability terms are added to enhance stability. However, stability comes at the expense of accuracy, and the stiff of the equation will also cause the first-order method to be inaccurate. Therefore, the first-order operator splitting method is extrapolated to the fourth-order. For the accurate spatial discretization, we use the long-stencil fourth-order central difference scheme which has been successfully applied to various nonlinear partial differential equations, such as fluid equations, phase field equations [46], [47], [48], [49]. Obviously, our algorithm is very accurate due to the fourth-order precision in both space and time. It is worth mentioning that the considered extrapolation approach can still be performed in parallel, which significantly speeds up the calculation.

Due to the double-well potential F(ϕ), the conventional Allen–Cahn type equations satisfy the maximum principle: if the initial value is bounded by 1ϕ0(x)1,then the entire solution is also bounded, i.e., 1ϕ(x,t)1,t>0.The maximum principle is a stronger stability under the infinity norm has a great effect on avoiding numerical oscillations. For the anisotropic Allen–Cahn equation in model Eq. (1.5), although it is difficult to derive such a theoretical maximum principle, a large number of numerical experiments show that the its numerical solution is bounded by 1 and 1. During numerical modeling, this property of boundedness is vital to prevent the numerical solution from blowing up. For the conventional Allan–Cahn equations, there have been many successful schemes that can preserve the discrete maximum principle unconditionally, such as operator splitting, convex splitting and stabilized semi-implicit methods, but they are also limited to first-order schemes [18], [26], [27]. The exponential time difference (ETD) method [20], [50], [51], [52], [53] is a temporal integration method that has become popular in recent years for solving Allen–Cahn type equations. lt precisely computes the integral of the linear part of the equation and interpolates the nonlinear part. With appropriate spatial discretizations, the resultant schemes can be linear, high-order, unconditional stable, and maximum principle preserving. For the anisotropic Allen–Cahn equation, due to nonlinear diffusion term, it is hard to guarantee the boundedness of the numerical solution obtained by the direction splitting method. In addition, even if there is a boundedness control of direction splitting method, the proposed fourth-order extrapolation method is a non-convex combination of several first-order solutions, the boundedness of the numerical solution cannot be maintained. To this end, we develop two post-processing techniques to force the boundedness of the numerical solution: bound preserving least-distance modification (BPLDM) [54] and cut-off approach [55]. The BPLDM is a least-squares optimization problem, which can unconditionally force the boundedness of the numerical solution. The global correction can be simplified to a directional correction. Different modifications in the same direction do not affect each other, so they can be performed in parallel. This is an improvement of its original version proposed in [54], which dramatically speeds up the computing. The cut-off approach cuts off the solution at each moment. In other words, values greater than 1 are forced to 1, and values less than −1 are forced to −1. Both methods maintain the boundedness of the numerical solution while improving accuracy and stability.

In short, our novel algorithm has the advantages of fast, accurate, stable and boundedness control. The rest of this article is organized as follows: In Section 2, we detail the spatiotemporal discretization scheme of the model. In Section 3, an extrapolation scheme is used to increase the order of convergence. In Section 4, two post-processing methods are proposed to improve the accuracy and stability of the numerical method. In Section 5, the convergence experiment and numerical simulation of dendritic crystal growth are carried out. Finally, some concluding comments are given in Section 6.

Section snippets

Numerical schemes

In this section, we introduce the time discretization and spatial discretization formats of the model Eq. (1.5) in detail.

Improve accuracy

In order to improve accuracy, the first-order operator splitting method is extrapolated to fourth-order. We take the interval (tn,tn+1] as an example. By the proposed stable first-order operator splitting method, the solution at tn+1 by δt, 12δt, 13δt and 14δt are denoted by ϕn+1,1, ϕn+1,2, ϕn+1,3 and ϕn+1,4, respectively. Similarly, we can define un+1,1, un+1,2, un+1,3 and un+1,4. Proceeding linear combinations ϕn+1=16ϕn+1,1+4ϕn+1,2272ϕn+1,3+323ϕn+1,4, un+1=16un+1,1+4un+1,2272un+1,3+323un+

Two post-processing approaches

The boundedness of the numerical solution is vital to numerical stability. With the above proposed method, even if there is a boundedness control of direction splitting method with linear stabilization, the proposed fourth-order extrapolation method is a non-convex combination of several first-order solutions, the boundedness of the numerical solution may not be maintained. Therefore, the long-time numerical simulation may fail with solution’s blow up. To force the boundedness of the high-order

Accuracy and stability test

We firstly tested the time convergence of the fourth-order operator splitting method. The exact solution is ϕ=sin(x)cos(y)cos(t),u=cos(x)sin(y)cos(t).Other initial conditions are set to Ω=[π,π]2,K=0.5,λ=D=1,τ=100,ϵ=0.06,ϵ4=0.05,S1=6,S2=4.We discretize the space with 1292 grid points so that the spatial discretization error is negligible relative to the time discretization error. To determine the order of convergence in time, we run a mesh refinement test in real time. By taking a linear

Conclusion

Based on a stabilized direction splitting method, combined with extrapolation and post-processing techniques, we propose an efficient parallel numerical method to solve the anisotropic phase-field dendritic crystal growth model. There related stability and error analyses are provided. Numerical experiments show that the algorithm can simulate 2D and 3D dendritic crystal growth processes quickly, accurately, and stably. Future work is to design an efficient direction-splitting scheme for a

CRediT authorship contribution statement

Yan Wang: Investigation, Software, Validation, Writing – original draft. Xufeng Xiao: Conceptualization, Methodology, Formal analysis, Software, Writing – review & editing. Xinlong Feng: Formal analysis, Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References (57)

  • ChenC. et al.

    Efficient numerical scheme for a dendritic solidification phase field model with melt convection

    J Comput Phys

    (2019)
  • YangX.

    Efficient linear, stabilized, second-order time marching schemes for an anisotropic phase field dendritic crystal growth model

    Comput Methods Appl Mech Engrg

    (2019)
  • ZhangJ. et al.

    A novel decoupled and stable scheme for an anisotropic phase-field dendritic crystal growth model

    Appl Math Lett

    (2019)
  • LiY. et al.

    A fast robust and accurate operator splitting method for phase field simulations of crystal growth

    J Cryst Growth

    (2011)
  • LiuC. et al.

    Operator splitting scheme for reaction–diffusion equations with detailed balance

    J Comput Phys

    (2021)
  • ZhaiS. et al.

    Fast explicit operator splitting method and time-step adaptivity for fractional non-local Allen–Cahn model

    Appl Math Model

    (2016)
  • ZhangC. et al.

    A second order operator splitting numerical scheme for the “good” Boussinesq equation

    Appl Numer Math

    (2017)
  • ChengK. et al.

    An energy stable fourth order finite difference scheme for the Cahn-Hilliard equation

    J Comput Appl Math

    (2019)
  • BurdakovO. et al.

    Monotonicity recovering and accuracy preserving optimization methods for postprocessing finite element solutions

    J Comput Phys

    (2012)
  • ChalmersB.

    Principles of solidification

  • HuangS. et al.

    Fundamentals of dendritic solidification, I and II

    Acta Metall

    (1981)
  • WeeksJ.D. et al.

    Dynamics of crystal growth

    Adv Chem Phys

    (1979)
  • ZhuM. et al.

    A modified cellular automaton model for the simulation of dendritic growth in solidification of alloys

    Isij Int

    (2001)
  • IhleT.

    Competition between kinetic and surface tension anisotropy in dendritic growth

    Eur Phys J B

    (2000)
  • GibouF. et al.

    A level set approach for the numerical simulation of dendritic growth

    J Sci Comput

    (2002)
  • KimY.T. et al.

    Computation of dendritic microstructures using a level set method

    Phys Rev E

    (2000)
  • JeongJ. et al.

    Phase-field model for three-dimensional dendritic growth with fluid flow

    Phys Rev E

    (2001)
  • KarmaA. et al.

    Three-dimensional dendrite-tip morphology at low undercooling

    Phys Rev E

    (2000)
  • Cited by (13)

    • An efficient numerical method for the anisotropic phase field dendritic crystal growth model

      2024, Communications in Nonlinear Science and Numerical Simulation
    • An efficient and fast adaptive numerical method for a novel phase-field model of crystal growth

      2024, Communications in Nonlinear Science and Numerical Simulation
    View all citing articles on Scopus

    This work is supported by the NSF of China (No. 12001466, No. U19A2079 and No. 12071406) the scientific research plan of universities in Xinjiang, PR China (No. XJEDU2020Y001 and No. XJEDU2020I001), the Research Fund from Key Laboratory of Xinjiang Province, PR China (No. 2020D04002), and the Graduate Student Research Innovation Program of Xinjiang, PR China (No. XJ2022G068).

    View full text