Research paperAn accurate and parallel method with post-processing boundedness control for solving the anisotropic phase-field dendritic crystal growth model☆
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: where , is the density field where for the solid and for the fluid. These two values are connected by a smooth transitional layer with the thickness . Among them, , and are all positive parameters, is the temperature, is the double well potential. is a function describing the anisotropic property that depends on the direction of the outer normal vector which is the interface normal defined as . For the system, the anisotropy coefficient is usually given by where is the number of folds of anisotropy, is the parameter for the anisotropy strength, and . When (i.e., fourfold anisotropy), for instance, can be easily reformulated in terms of the phase-field variable , namely, for and for 3D, By adopting the Allen–Cahn type (gradient flow) relaxation dynamics for the dendritic crystal growth, one obtains the governing dynamical equations via the variational approach, which reads as follows: where is the variational derivative of . In , it reads as whereas in 3D, In the system Eq. (1.5), is the mobility constant that is chosen either as a constant or a function of [16], [21], constant is the diffusion rate of the temperature, is the variational derivative of the total energy with respect to , the function accounts for the generation of latent heat and it is a phenomenological functional taking the form preserving the minima of at independently of the local value of . For , we choose and . The periodic boundary condition is considered in this study. The model equations Eq. (1.5) follow the dissipative energy law. By taking the inner product of the first equation of Eq. (1.5) with , and of the second equation of Eq. (1.5) with , using the integration by parts and combining the obtained two equalities, we obtain
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 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 , the conventional Allen–Cahn type equations satisfy the maximum principle: if the initial value is bounded by then the entire solution is also bounded, i.e., 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 and . 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 as an example. By the proposed stable first-order operator splitting method, the solution at by , , and are denoted by , , and , respectively. Similarly, we can define , , and . Proceeding linear combinations
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 Other initial conditions are set to We discretize the space with 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)
- et al.
Nonlinear morphological control of growing crystals
Physica D
(2005) - et al.
Crystal growth and dendlritic solidification
J Comput Phys
(1992) - et al.
Numerical simulation of dendritic solidification with convection: Two-dimensional geometry
J Comput Phys
(2002) - et al.
Multiscale finite-difference-diffusion-Monte-Carlo method for simulating dendritic solidification
J Comput Phys
(2000) - et al.
Thermodynamically-consistent phase-field models for solidification
Physica D
(1993) - et al.
Linear energy stable and maximum principle preserving semi-implicit scheme for Allen–Cahn equation with double well potential
Commun Nonlinear Sci Numer Simul
(2021) - et al.
Prediction of dentric growth and microsegregation patterns in a binary alloy using the phase field method
Acta Metall Mater
(1995) - et al.
A pseudo-front tracking technique for the modelling of solidification microstructures in multi-component alloys
Acta Mater
(2002) - et al.
Efficient energy stable schemes for isotropic and strongly anisotropic Cahn-Hilliard systems with the Willmore regularization
J Comput Phys
(2018) - et al.
A weakly nonlinear, energy stable scheme for the strongly anisotropic Cahn-Hilliard equation and its convergence analysis
J Comput Phys
(2020)
Efficient numerical scheme for a dendritic solidification phase field model with melt convection
J Comput Phys
Efficient linear, stabilized, second-order time marching schemes for an anisotropic phase field dendritic crystal growth model
Comput Methods Appl Mech Engrg
A novel decoupled and stable scheme for an anisotropic phase-field dendritic crystal growth model
Appl Math Lett
A fast robust and accurate operator splitting method for phase field simulations of crystal growth
J Cryst Growth
Operator splitting scheme for reaction–diffusion equations with detailed balance
J Comput Phys
Fast explicit operator splitting method and time-step adaptivity for fractional non-local Allen–Cahn model
Appl Math Model
A second order operator splitting numerical scheme for the “good” Boussinesq equation
Appl Numer Math
An energy stable fourth order finite difference scheme for the Cahn-Hilliard equation
J Comput Appl Math
Monotonicity recovering and accuracy preserving optimization methods for postprocessing finite element solutions
J Comput Phys
Principles of solidification
Fundamentals of dendritic solidification, I and II
Acta Metall
Dynamics of crystal growth
Adv Chem Phys
A modified cellular automaton model for the simulation of dendritic growth in solidification of alloys
Isij Int
Competition between kinetic and surface tension anisotropy in dendritic growth
Eur Phys J B
A level set approach for the numerical simulation of dendritic growth
J Sci Comput
Computation of dendritic microstructures using a level set method
Phys Rev E
Phase-field model for three-dimensional dendritic growth with fluid flow
Phys Rev E
Three-dimensional dendrite-tip morphology at low undercooling
Phys Rev E
Cited by (13)
Second-order accurate and unconditionally stable algorithm with unique solvability for a phase-field model of 3D volume reconstruction
2024, Journal of Computational PhysicsFast and stable dimension splitting simulations for the hydrodynamically coupled three-component conserved Allen–Cahn phase field model
2024, International Journal of Multiphase FlowNumerical simulation for the conserved Allen–Cahn phase field model of two-phase incompressible flows by an efficient dimension splitting method
2024, Communications in Nonlinear Science and Numerical SimulationAn efficient numerical method for the anisotropic phase field dendritic crystal growth model
2024, Communications in Nonlinear Science and Numerical SimulationAn efficient and fast adaptive numerical method for a novel phase-field model of crystal growth
2024, Communications in Nonlinear Science and Numerical SimulationA highly efficient variant of scalar auxiliary variable (SAV) approach for the phase-field fluid-surfactant model
2023, Computer Physics Communications
- ☆
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).