{'fraction_non_alphanumeric': 0.04155161078238001, 'fraction_numerical': 0.015822923515231208, 'mean_word_length': 4.529810954920019, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
New efficient time-stepping schemes for the anisotropic phase-field dendritic crystal growth model

In this paper, we propose and analyze a first-order and a second-order time-stepping schemes for the anisotropic phase-field dendritic crystal growth model. The proposed schemes are based on an auxiliary variable approach for the Allen-Cahn equation and delicate treatment of the terms coupling the Allen-Cahn equation and temperature equation. The idea of the former is to introduce suitable auxiliary variables to facilitate construction of high order stable schemes for a large class of gradient flows. We propose a new technique to treat the coupling terms involved in the crystal growth model and introduce suitable stabilization terms to result in totally decoupled schemes, which satisfy a discrete energy law without affecting the convergence order. A delicate implementation demonstrates that the proposed schemes can be realized in a very efficient way. That is, it only requires solving four linear elliptic equations and a simple algebraic equation at each time step. A detailed comparison with existing schemes is given, and the advantage of the new schemes are emphasized. As far as we know this is the first second-order scheme that is totally decoupled, linear, unconditionally stable for the dendritic crystal growth model with variable mobility parameter.


Introduction
On one side, dendritic growth is a very common phenomenon in nature. We are all familiar with the way how trees grow by spreading branches and roots from the main trunk. This is where the name "dendritic" comes from, although the term "dendrite" itself is used to describe branched projections of neurons. On the other side, dendritic growth phenomena and the shapes of growing crystals are of fundamental interest to physicists and are of practical importance to engineers.
Crystal dendritic growth is one of the most extensively studied topics in the scientific literature.
Crystallization proceeds through the competition between thermodynamics -driven by the local undercooling of the liquid ahead of the solidification front -and the ability of the system to diffuse latent heat of fusion away from the solid-liquid interface. It usually forms natural fractal microstructure, so-called dendrites, which are the ubiquitous crystal form in freezing alloys and supercooled melts. When the molten material is supercooled below the freezing point of the solid, a spherical solid nucleus grows in the undercooled melt initially. Along with some preferred directions of growth, the solid form begins to express some protrusion accompanied by steeper concentration gradients at its end. Dendritic microstructures formed during solidification/freezing play a key role in properties of the final solid material. Understanding these microstructures is therefore considered essential for controlling basic solidification and crystal growth processes.
The first phase-field models were suggested for numerically simulating dendritic growth in 1980s; see, e.g., [3,8,11,17]. This concept has been validated by comparison with theoretical predictions and experimental measurements and is applied to a broad range of investigations in materials science [2]. Nowadays, the phase-field method has emerged as a powerful tool for modelling and simulation of crystal dendritic growth. In contrast to sharp interface approaches with interfaces of zero thickness, the phase-field model introduces a smooth phase-field variable by a diffuse interface profile to distinguish between the solid and liquid phases. In this model the complicated topological changes of a solid-liquid interface can be handled in an easy way without the need of the explicit tracking of the interface. The phase field is considered as an order parameter which is introduced to describe the moving interfacial boundary between unstable and stable phases during phase transformation processes. By asymptotic expansions, it can be shown that the phase-field methods relate to classical sharp interface models such as Hele-Shaw type models and Stefan problems in the limit of zero interfacial thickness, see, e.g. [4].
Another advantage of the phase-field approach is that the governing set of equations in the model can be naturally derived from an energy-based variational principle. The variational framework of phase-field formulations makes them thermodynamically consistent and physically attractive in modeling the general phase-field dendritic crystal growth model [9,10,13,15,[20][21][22]30].
In this paper, we will focus on the numerical approximations for the anisotropic phase-field dendritic crystal growth model proposed in [13,14]. The model is composed of two coupled nonlinear equations. One is the phase-field equation that governs the anisotropy of the crystal. The other is the heat equation that controls heat diffusion of the system. It is shown that this nonlinear coupled system satisfies a thermodynamically consistent energy dissipation law. The main aim of this paper is to design efficient numerical schemes for this nonlinear crystal growth model, which satisfies a discrete version of the energy dissipation law. In fact, constructing schemes that preserve the discrete energy dissipation law for similar models has been subject of many recent papers [1, 5-7, 12, 25, 26, 29, 39]. Although large amounts of works have been devoted to numerical approximation for phase-field dendritic crystal growth models; see, e.g., [18,19,23,40] and the references therein, there is still a need for efforts on developing low-cost, stable, and high order schemes for such models.
The main difficulties in constructing highly efficient schemes for the dendritic crystal growth models come from: 1) the double-well energy potential and the stiffness associated with the interfacial width in the phase equation; 2) the anisotropic coefficient; 3) the nonlinear interaction terms in both the heat equation and the phase field equation. Let's briefly review recent progress in this direction. Firstly, to overcome the difficulty caused by the nonlinearity and the thin interface in the phase field equation, schemes based on the invariant energy quadratization (IEQ) [31] and scalar auxiliary variable approach (SAV) [24,27] have been proposed: a decoupled stable but only first-order scheme in [37] and a second-order stable scheme but fully coupled scheme in [32].
Due to the presence of the nonlinear phase term in the heat equation and the interaction term in the phase equation, it seems not easy to design a fully-decoupled, second-order accurate and energy stable scheme. For example, the splitting method used in [37] is not directly extendable to a second-order discretization for the time derivative of the phase function in the heat equation. In the case of constant mobility parameter, [36] proposes a decoupling, linear, second-order accurate, and unconditionally stable scheme by using multi-auxiliary variables. A similar technique was used in [33][34][35] to deal with some coupling models. However, compared with the traditional SAV approach, multi-auxiliary variable approach means extra computational cost since more equations are to be solved. In particular, the second-order scheme proposed in [36] is a four-step scheme, thus is much more computationally expensive.
The main purpose of the present paper is to propose easy-to-implement, second-order accurate, and unconditionally stable schemes for the anisotropic phase-field dendritic crystal growth model.
First of all, we rewrite the time derivative term of the phase function in the heat equation into an equivalent form, which allows to design a three-step second-order scheme. The idea is to introduce a suitable auxiliary variable to the Allen-Cahn equation and a new technique to treat the coupling terms. Then some carefully chosen stabilization terms are added to result in totally decoupled schemes that satisfy a discrete energy law without losing the convergence order. A careful examination shows that the proposed schemes can be implemented by only solving four linear elliptic equations and a simple algebraic equation. As far as we know this is the first second-order scheme that is totally decoupled, linear, unconditionally stable for the dendritic crystal growth model variable mobility parameter. In the case of constant mobility parameter, compared with [36] (a four-layer scheme that requires solve five linear elliptic equations with constant coefficients and some algebraic equations), our scheme is a three-layer scheme that only needs solve four linear elliptic equations with constant coefficients and a simple algebraic equation.
The rest of the paper is organized as follows. In Section 2, we describe the phase-field dendritic crystal growth model, and present the equivalent reformulation using auxiliary variables. In Section 3 we propose a first-order unconditionally stable time-stepping scheme, and prove the energy decay property of the proposed scheme. Section 4 is devoted to construct and analyze a second-order, linear, decoupled, and unconditionally stable scheme. The implementation detail is also presented to show that the scheme can be efficiently realized through solving a set of decoupled, linear elliptic equations. We give in Section 5 some numerical examples to verify the efficiency of the proposed methods. Finally, the paper ends with some concluding remarks.

Governing equations and auxiliary variable reformulation
2.1. Governing equations. We are interested in numerically solving the following anisotropic phase-field dendritic crystal growth model equations in the domain Ω ⊂ R 2 : where φ(x, t) is the phase function to label the liquid and solid phases, (φ) > 0 is the mobility parameter that can be chosen either as a constant [36], or as a function of φ [32]. ε > 0 is a parameter used to control the interface width, λ is the linear kinetic coefficient. In Equation (2.2), T (x, t) is the scaled temperature, D is the constant diffusion rate of the temperature, and K is the latent heat parameter that controls the speed of heat transfer along with the interface. It is worth noting that the efficiency of the schemes we propose below covers the case D is a function of φ [28].
The function h(φ) is defined by which represents a generation of latent heat. Following the phenomenological free energy used in [13], we consider here where F (φ) = 1 4 φ 2 − 1 2 is the double-well type Ginzburg-Landau potential. κ(·) in (2.3) is a function describing the anisotropic property, which takes the form [13,14]: κ(∇φ) = 1 + σ cos(mθ), (2.4) where m is a model number of anisotropy, σ is the parameter for the anisotropy strength, and θ = arctan φy φx . The variational derivative of E with respect to φ is: where H(φ) is the variational derivative of κ(∇φ), and f (φ) = F (φ). In the case m = 4, a direct calculation shows For convenience, we only consider that the equations (2.1) and (2.2) are subject to the Neumann boundary conditions ∂φ ∂n although other boundary conditions such as the periodic conditions are possible.
Now we briefly recall main property of (2.1)-(2.6). A key property of the model is that it satisfies an energy law, which can be derived by taking the inner product of (2.1) with −φ t and using integration by parts: Then taking the inner product of (2.2) by − λ εK T gives: Combining the above two equalities gives the following energy law where · denotes the standard L 2 (Ω) norm. This means that the energy E(φ, T ) decays in time during the crystal-growing process.

Auxiliary variable reformulation.
The main purpose of this paper is to develop novel efficient schemes for the anisotropic crystal growth model (2.1)-(2.6). We start with an auxiliary variable approach, which will be used later to construct time-stepping schemes for the phase field equation (2.1). We define the variable where S 1 and S 2 are two positive constants, 0 < S 1 < (1 − σ) 2 , B is a positive constant used to make E 1 positive. Notice κ 2 (∇φ) ≥ (1 − σ) 2 and F (φ) is a quartic polynomial, one can verify that )dx is bounded from below. Therefore such a constant B exists. The introduction of the constants S 1 and S 2 is inspired by the work [32,36,38]. We will see that these constants help in ensuring the H 1 -stability of the phase function.
Using the auxiliary variable R(t), the total free energy (2.3) can be rewritten as and the original equations (2.1)-(2.2) can be reformulated into the following equivalent form: where M (φ) = 1 (φ) , and Obviously, the equation (2.10c) can be obtained by taking the time derivative of the auxiliary variable R(t). The initial conditions for (2.10) take Remark 2.1. To separate the computation of different unknown functions, one may think about the splitting method used in [37] or explicit treatment of φ t as in [36]. However the former is unlikely to lead to a second-order scheme, while the latter may result in expensive four-layer computation [36]. To construct more efficient second-order scheme, our idea here is to split ( Now we are in a position to construct and analyze our schemes for the anisotropic phase-field dendritic crystal growth model (2.1)-(2.6). To better follow the main argument, let's start with a first-order scheme.

A first order scheme and stability analysis
Let τ > 0 be the time step size, t n = nτ, 0 ≤ n ≤ N , T = N τ . We propose the following scheme: , S 3 and S 4 are two extra positive stabilization parameters.
Before carrying out the stability analysis, the scheme (3.1) is worthy of some explanation. First, we notice that the coupling terms in both the phase field equation and the temperature equation are treated explicitly. This is for ease of calculation. The implicit treatment of the coupling terms in the auxiliary variable equation, i.e., eq.(3.1c), may make the implementation difficult. However, as we will see in the next section, the extra terms added to the scheme, i.e., terms involving the parameters S i , play a dual role. On one side, some of extra terms are useful in decoupling the calculation of different unknowns. On the other side, some other extra terms help in enhancing the stability. For example, the term S 3 ε 2 φ n+1 − φ n is used to balance the explicit treatment of 1 ε 2 f (φ) in the phase field equation, and the term S 4 ∆ φ n+1 − φ n has purpose to balance the explicit treatment of the gradient term. The last point we want to emphasize is that the extra terms have the same order as the approximation to the time derivatives, thus do not affect the overall accuracy.
For example, we can check that the term . Therefore, formally the convergence of the scheme (3.1) is first order.
3.1. Stability analysis. In the following theorem, we establish the stability result for the scheme (3.1). That is, we prove that a discrete "energy" decays in time, and consequently the discrete solution remains bounded during the time stepping.
Theorem 3.1. Let {φ n , T n , R n } be the solution of the discrete problem (3.1). Then the following discrete energy law holds:
To summarize, the algorithm involves the solution of two variable coefficient linear elliptic equations, two constant coefficient linear elliptic equations, and one algebraic equation. Furthermore, it is a three-layer scheme, compared to the four-layer scheme proposed in [36]. In the case of constant mobility, our scheme is reduced to solve four elliptic equations with constant coefficients and one algebraic equation.

A second order scheme
4.1. Construction of the scheme. For ease of notation we will useφ n+1 to denote 2ϕ n − ϕ n−1 .
The second order scheme we propose reads: Obviously, to start the calculation the scheme (4.1) must be accompanied by a suitable one-step scheme to compute φ 1 , T 1 , R 1 . This can be done, for example, by employing the first step of the scheme (3.1).
Intuitively this is a second-order scheme since all involved terms are approximated with secondorder precision. Although rigorous proof is not available for the time being, the convergence order of the scheme will be confirmed through a series of numerical tests.

Stability analysis.
Below we prove the stability of the scheme (4.1). The stability analysis will make use of the following well-known identities: Theorem 4.1. Let {φ n , T n , R n } be the solution of the discrete problem (4.1). Then for n ≥ 1 it satisfies the discrete energy law: where E n is defined by and Q n+1 is defined by Proof First, we take the inner product of (4.1a) with 2 M (φ n+1 ) (3φ n+1 − 4φ n + φ n−1 ), and (4.1b) with 2(3φ n+1 − 4φ n + φ n−1 ). Then we sum up the resulting equations and use the identities (4.3) and (4.4) to obtain By virtue of (4.1c) and (4.1a), we can rewrite (4.1c) as Multiplying the both sides by 8τ R n+1 , we have Furthermore, it follows from taking the inner product of (4.1d) with 4λτ εK T n+1 , and using (4.3): Finally we sum up (4.7), (4.8), and (4.9) to conclude. This completes the proof.
The second order scheme (4.1) can also be implemented in an efficient way. We are not going to describe the implementation details of this scheme since it is very similar to the first order scheme (3.1) explained in the previous section. Nevertheless we give here the algebraic equation to compute the auxiliary variable ξ n+1 :

Numerical experiments
In order to illustrate the performance of the proposed numerical method and confirm our analysis results, several numerical tests are carried out and presented in this section. For the sake of convenience, we consider the numerical examples with the square domain Ω = (−1, 1) 2 and fourfold anisotropy, i.e., m = 4 in (2.4). We consider constant mobility (φ) in the calculation. The spatial discretization is a Legendre-Galerkin spectral method. The approximation space for the phase function φ and the temperature T is IP N (Ω), where IP N (Ω) denotes the space of polynomials of degree ≤ N at each space direction. Since the focus of our numerical tests is the verification of the stability and convergence order of the time-stepping schemes, we will take N = 128, which is large enough so that the spatial discretization errors are negligible compared with the temporal one.
Example 5.1. (Accuracy test) As the first example, we test the convergence order of the schemes (3.1) and (4.1). In the case one, we fabricate two forcing functions in (2.1) and (2.2) so that the exact solution to (2.1)-(2.6) is Case-I φ(x, y, t) = sin(t) cos(πx) cos(πy), T (x, y, t) = sin(t) cos 2 (πx) cos 2 (πy). In the case two, we choose the initial conditions: We know from (2.9) that E n defined in (4.6) can be regarded as a discrete version of the original energy functional E(φ n , T n ) defined in (2.3). According to (2.7) and Theorem 4.1, both E(φ n , T n ) and E n should be monotonically decreasing with the time step n ≥ 1. Notice that E n is usually called as modified energy functional, which is not necessarily an approximation to the original energy functional E(φ n , T n ).
The modified energy functional E n , original energy functional E(φ n , T n ), and ξ n computed by the scheme (4.1) with different time step sizes are presented in Figure         However, the dissipation of E(φ n , T n ) loses monotonicity during some time period(τ ≥ 10) for the time step sizes bigger than 10. This is most likely caused by imprecise calculation with large time step sizes. This guess is supported by the observation from Figure 5.2(c), in which the time evolution of ξ n is given. Remember that ξ n , defined in (4.2), should be close to 1 if the approximation is good enough. The results presented in Figure 5.2(c) demonstrate good accuracy on ξ n when τ is not very large, say less than 10. However when τ = 10, ξ n becomes oscillatory and error becomes visible. When τ increases to 100, the computed auxiliary variable ξ n is not any more close to 1.
This implies that the numerical solution is not accurate enough for τ ≥ 10, leading to a violation of the monotonic dissipation of the original energy. Thus in practice, it is not recommended to use time step sizes too large, although the calculation can always be stable.
Next test concerns impact of the parameters S 1 , S 2 , S 3 , and S 4 on the quality of numerical solutions. Since the auxiliary variable is indicative of the quality of the computed solutions, we only report the computed values of ξ n . We consider a benchmark problem, which has been extensively studied; see, e.g., [14,16,34,37].
We set the initial conditions as follows: We let the latent heat parameter K vary. The spatial spectral discretization uses the polynomial space of degree 512 at each direction.
In starting with the same small circle. Moreover, it is seen that the parameter K affects the width of the branch: larger is K, thinner is the width of the branches, and sharper are the tips.
The isocontours of the temperature T at the last moment of each simulation is plotted in the figures 5.4 (e). It is observed that the contours of the temperature T take similar dendrite crystal shape as the phase field. This is due to the fact that the heat is propagating only at the interface.
The dissipation behavior of the modified energy E n in time is shown in Figure 5.5(a). The monotonic decay feature of E n for all tested K reflects good stability property of the scheme used in the calculation. Finally we give in Figure 5.5(b) evolution of the area of the crystal, defined by the quantity Ω 1+φ 2 dx, for several values of K. We see that the area of the crystal keeps increasing during the simulation. This is in a good agreement with the existing results; see, e.g., [14,32,36].

Concluding remarks
We have proposed a class of new time-stepping schemes for the anisotropic phase-field dendritic crystal growth model. The proposed schemes were constructed based on an auxiliary variable approach for the Allen-Cahn equation and sophisticated treatment of the terms coupling the Allen-Cahn equation and temperature equation. In particular, the new reformulation of the model introduced in the paper plays a key role in developing efficient schemes. Thanks to the carefully chosen extra terms added to the time discretization, we were able to construct a second-order scheme, which is linear, decoupled, uniquely solvable, and unconditionally stable. A detailed comparison with existing schemes is given, and the advantage of the new schemes are emphasized. The stability property of the proposed schemes was rigorously established, while the convergence rate was carefully examined through a series of numerical tests. Our analysis and numerical experiments demonstrated the efficiency of the proposed method. It seems to us that the approach proposed in this paper is extendable to more complex models such as those studied in [33][34][35].