Linear Second Order Energy Stable Schemes of Phase Field Model with Nonlocal Constraints for Crystal Growth

TL;DR

This paper introduces linear, energy-stable schemes for a nonlocal Allen-Cahn model that conserves phase mass, compares dynamics with other models, and demonstrates its effectiveness in simulating crystal growth.

Contribution

It develops and verifies a new class of linear, second-order, energy-stable schemes for a nonlocal Allen-Cahn model with mass conservation, offering an alternative to Cahn-Hilliard for crystal growth simulation.

Findings

Schemes are unconditionally energy stable and linearly solvable.

Nonlocal Allen-Cahn model shows intermediate dynamics between Allen-Cahn and Cahn-Hilliard.

Benchmark examples validate the model's predictive capability.

Abstract

We present a set of linear, second order, unconditionally energy stable schemes for the Allen-Cahn model with a nonlocal constraint for crystal growth that conserves the mass of each phase. Solvability conditions are established for the linear systems resulting from the linear schemes. Convergence rates are verified numerically. Dynamics obtained using the nonlocal Allen-Cahn model are compared with the one obtained using the classic Allen-Cahn model as well as the Cahn-Hilliard model, demonstrating slower dynamics than that of the Allen-Cahn model but faster dynamics than that of the Cahn-Hillard model. Thus, the nonlocal Allen-Cahn model can be an alternative to the Cahn-Hilliard model in simulating crystal growth. Two Benchmark examples are presented to illustrate the prediction made with the nonlocal Allen-Cahn model in comparison to those made with the Allen-Cahn model and the Cahn- Hillard model.

Figures

Equations

$$\displaystyle\begin{array}[]{l}F=\int_{\Omega}[\frac{\phi}{2}(-\varepsilon+(\nabla^{2}+1)^{2})\phi+\frac{\phi^{4}}{4}]\mathrm{d{\bf r}},\end{array}$$

$$\displaystyle\begin{array}[]{l}F=\int_{\Omega}[\frac{\phi}{2}(\nabla^{4}+2a\nabla^{2}+\alpha)\phi+\frac{\phi^{4}}{4}]\mathrm{d{\bf r}}\\ =\int_{\Omega}[\frac{1}{2}\|\nabla^{2}\phi\|^{2}+a\phi\nabla^{2}\phi+\frac{\alpha}{2}\phi^{2}+\frac{\phi^{4}}{4}]\mathrm{d{\bf r}}\\ =\int_{\Omega}[\frac{1}{2}(\|\nabla^{2}\phi\|+a\phi)^{2}+\frac{\alpha-a^{2}}{2}\phi^{2}+\frac{\phi^{4}}{4}]\mathrm{d{\bf r}}\\ =\int_{\Omega}[\frac{1}{2}\|\nabla^{2}\phi\|^{2}-a\|\nabla\phi\|^{2}+\frac{\alpha}{2}\phi^{2}+\frac{\phi^{4}}{4}]\mathrm{d{\bf r}},\end{array}$$

$$\mathbf{n}\cdot \nabla \varphi =0,\mathbf{n}\cdot \nabla {\nabla }^2\varphi =0$$

$$\displaystyle\begin{array}[]{l}\frac{\partial\phi}{\partial t}=-M\mu,\quad\text{in}\,\Omega,\\ \mu=\frac{\delta F}{\delta\phi}=(\nabla^{2}+a)^{2}\phi+(\alpha-a^{2})\phi+\phi^{3},\end{array}$$

$$\displaystyle\begin{array}[]{l}\frac{dF}{dt}=-\int_{\Omega}\mu M\mu\mathrm{d{\bf r}}+\int_{\partial\Omega}{\bf n}\cdot[(\frac{\partial f}{\partial\nabla\phi}-\nabla\frac{\partial f}{\partial\nabla^{2}\phi})\phi_{t}+\frac{\partial f}{\partial\nabla^{2}\phi}\nabla\phi_{t}]\mathrm{d}{\bf s}.\end{array}$$

$$\mathbf{n}\cdot \nabla \varphi =0,\mathbf{n}\cdot (\frac{\partial f}{\partial \nabla \varphi }-\nabla \frac{\partial f}{\partial {\nabla }^2\varphi })=0.$$

$$\mathbf{n}\cdot \nabla \varphi =0,\mathbf{n}\cdot \nabla \frac{\partial f}{\partial {\nabla }^2\varphi }=\mathbf{n}\cdot \nabla {\nabla }^2\varphi =0.$$

$$\frac{\partial f}{\partial {\nabla }^2\varphi }=0,\mathbf{n}\cdot (\frac{\partial f}{\partial \nabla \varphi }-\nabla \frac{\partial f}{\partial {\nabla }^2\varphi })=0.$$

$$({\nabla }^2+a)\varphi =0,\mathbf{n}\cdot (\frac{\partial f}{\partial \nabla \varphi }-\nabla \frac{\partial f}{\partial {\nabla }^2\varphi })=0.$$

$$\mathbf{n}\cdot \nabla {\varphi }_t=0,{\varphi }_t=-\beta \mathbf{n}\cdot (\frac{\partial f}{\partial \nabla \varphi }-\nabla \frac{\partial f}{\partial {\nabla }^2\varphi }),$$

$$-{\int }_{\partial \Omega }\beta [\mathbf{n}\cdot (\frac{\partial f}{\partial \nabla \varphi }-\nabla \frac{\partial f}{\partial {\nabla }^2\varphi }){]}^{2}ds.$$

$$\frac{dF}{dt}=-{\int }_{\Omega }\mu M\mu \mathrm{d}\mathbf{r}-{\int }_{\partial \Omega }\beta [\mathbf{n}\cdot (\frac{\partial f}{\partial \nabla \varphi }-\nabla \frac{\partial f}{\partial {\nabla }^2\varphi }){]}^{2}ds.$$

$$\displaystyle M=\left\{\begin{array}[]{ll}M_{0},&\hbox{Allen-Cahn},\\ -\nabla\cdot M_{0}\nabla,&\hbox{Cahn-Hilliard},\end{array}\right.$$

$$\displaystyle\begin{array}[]{l}\frac{\partial\phi}{\partial t}=-M\mu,\quad\text{in}\,\Omega,\\ \frac{\partial\phi}{\partial n}=0,\frac{\partial\nabla^{2}\phi}{\partial n}=0,\quad\text{in}\,\partial\Omega,\\ \phi|_{t=0}=\phi(0,\mathbf{x}),\end{array}$$

$$\displaystyle\begin{array}[]{l}\mu=\frac{\delta F}{\delta\phi}=(\nabla^{2}+a)^{2}\phi+(\alpha-a^{2})\phi+\phi^{3}.\\ \end{array}$$

$$\displaystyle\begin{array}[]{l}\frac{dF}{dt}=\int_{\Omega}\frac{\delta F}{\delta\phi}\phi_{t}\mathrm{d{\bf r}}=-\int_{\Omega}\mu(M\mu)\mathrm{d{\bf r}}\leq 0,\end{array}$$

$$\displaystyle\begin{array}[]{l}\frac{\partial\phi}{\partial t}=\nabla\cdot(M\nabla\mu),\quad\text{in}\,\Omega,\\ \frac{\partial\phi}{\partial n}=0,\frac{\partial\nabla^{2}\phi}{\partial n}=0,\frac{\partial\mu}{\partial n}=0\quad\text{in}\,\partial\Omega,\\ \phi|_{t=0}=\phi(0,\mathbf{x}),\end{array}$$

$$\displaystyle\begin{array}[]{l}\mu=\frac{\delta F}{\delta\phi}=(\nabla^{2}+a)^{2}\phi+(\alpha-a^{2})\phi+\phi^{3}.\\ \end{array}$$

$$\displaystyle\begin{array}[]{l}\frac{dF}{dt}=\int_{\Omega}\frac{\delta F}{\delta\phi}\phi_{t}\mathrm{d{\bf r}}=-\int_{\Omega}\nabla\mu M\nabla\mu\mathrm{d{\bf r}}\leq 0,\end{array}$$

$$\displaystyle\begin{array}[]{l}F=\int_{\Omega}[\frac{1}{2}\|\nabla^{2}\phi\|^{2}-a\|\nabla\phi\|^{2}+\frac{\alpha}{2}\phi^{2}+\frac{\phi^{4}}{4}]\mathrm{d{\bf r}}+\frac{\eta}{2}(\int_{\Omega}\phi(t)\mathrm{d{\bf r}}-{\bf M_{0}})^{2},\end{array}$$

$$\displaystyle\begin{array}[]{l}\frac{\partial\phi}{\partial t}=-M\tilde{\mu},\\ \frac{\partial\phi}{\partial n}=0,\frac{\partial\nabla^{2}\phi}{\partial n}=0,in\;\partial\Omega,\\ \phi|_{t=0}=\phi(0,\mathbf{x})\end{array}$$

$$\displaystyle\begin{array}[]{l}\tilde{\mu}=\mu+\sqrt{\eta}\zeta,\quad\zeta=\sqrt{\eta}(\int_{\Omega}\phi(t)\mathrm{d{\bf r}}-{\bf M_{0}}).\\ \end{array}$$

$$\displaystyle\begin{array}[]{l}\frac{dF}{dt}=\int_{\Omega}\frac{\delta F}{\delta\phi}\phi_{t}\mathrm{d{\bf r}}=\int_{\Omega}\tilde{\mu}(-M\tilde{\mu})\mathrm{d{\bf r}}\leq 0,\end{array}$$

$$\tilde{F}=F-L({\int }_{\Omega }\varphi (t)\mathrm{d}\mathbf{r}-{\mathbf{M}}_0).$$

$$\displaystyle\begin{array}[]{l}\frac{\partial\phi}{\partial t}=-M\tilde{\mu},\\ \frac{\partial\phi}{\partial n}=0,\frac{\partial\nabla^{2}\phi}{\partial n}=0,in\;\partial\Omega,\\ \phi|_{t=0}=\phi(0,\mathbf{x}),\end{array}$$

$$\displaystyle\begin{array}[]{l}\tilde{\mu}=\mu-L,\quad L=\frac{1}{\int_{\Omega}M\mathrm{d{\bf r}}}\int_{\Omega}[M\tilde{\mu}]\mathrm{d{\bf r}}.\end{array}$$

$$\displaystyle\begin{array}[]{l}\frac{dF}{dt}=\int_{\Omega}\frac{\delta F}{\delta\phi}\phi_{t}\mathrm{d{\bf r}}=\int_{\Omega}\tilde{\mu}(-M\tilde{\mu})\mathrm{d{\bf r}}\leq 0,\end{array}$$

$$\displaystyle\begin{array}[]{l}\phi=\phi^{ss}+\delta v.\end{array}$$

$$\displaystyle\begin{array}[]{l}\frac{\partial\delta v}{\partial t}=-M((\nabla^{2}+a)^{2}\delta v-a^{2}\delta v+(\alpha+3{(\phi^{ss})}^{2})\delta v).\end{array}$$

$$\delta v=\sum _{k=0,l=0}^{\infty }{a}_{kl}cos(kx)cos(ly),$$

Full text

Keywords: Allen-Cahn equation with nonlocal constraints, Phase field model,Crystal growth, Energy stable schemes, Energy quadratization.

1 Introduction

Phase field crystal (PFC) growth model, developed as an extension to the phase field formalism [12, 11, 35, 25, 27], has been successfully applied to various applications in materials science across different time scales [13, 34, 1], capturing the interaction between material defects [12] and modeling the microstructure evolution [12, 11, 13, 2, 3, 42, 22, 4, 29, 33]. It is a challenge to develop efficient and stable numerical algorithms to faithfully simulate dynamics described by the PFC model. The PFC model is thermodynamically consistent in that the free energy of the thermodynamic model is dissipative. Numerical algorithms that respect the free energy dissipation property at the discrete level are known as energy stable schemes.

The Cahn-Hilliard equation is a popular phase field model for crystal growth because of its mass (or volume) preserving property. However, the Cahn-Hilliard equation for the crystal growth problem is of up to the 6th order spatial derivative. Searching for a lower order phase field model that can also preserve mass and free energy dissipation properties is therefore a viable alternative. Allen-Cahn equation is a popular phase field model which normally has lower spatial derivatives than the Cahn-Hilliard model. It describes relaxation dynamics of the thermodynamical system to equilibrium. However, in the case of a phase field description, when the phase variable represents the mass fraction or the volume fraction of a material component, this model does not warrant the conservation of mass of that component. In order to conserve mass, the free energy functional has to be augmented by a mass preserving penalty term or with a Lagrange multiplier [30, 49, 36, 28, 10]. This thus modifies the Allen-Cahn equation into a nonlocal equation. We call this the nonlocal Allen-Cahn model or the Allen-Cahn model with a nonlocal constraint. Rubinstein and Sternberg studied the Allen-Cahn model with a mass constraint analytically and compared it with the Cahn-Hilliard model [36]. Their result seems to favor using the Allen-Cahn model with a mass constraint in place of the Cahn-Hilliard model when studying interfacial dynamics of incompressible, immiscible multi-component material systems.

For the classical Allen-Cahn equation as well as the Cahn-Hilliard equation, there have been several popular numerical approaches to construct energy stable schemes for the equations, including the convex splitting approach [14, 17, 16, 47, 37, 45, 23] , the stabilizing approach [40, 9, 7], the energy quadratization (EQ) approach [57, 52, 21, 20] and the scalar auxiliary variable approach [39, 8, 46]. Recently, the energy quadratization (EQ) and its reincarnation in the scalar auxiliary variable (SAV) method have been applied to a host of thermodynamic and hydrodynamic models owing to their simplicity, ease of implementation, computational efficiency, linearity, and most importantly their energy stability property [53, 52, 21, 20, 51, 54, 55, 6, 48, 58, 56, 59, 8, 46]. Ones have shown that these strategies are general enough to be useful for developing energy stable numerical approximations to any thermodynamically consistent models, i.e., the models satisfy the second law of thermodynamics or are derived from the Onsager principle [31, 32, 57]. The convex splitting approach, scalar auxiliary variable, energy quadratization approach and other methods have been applied to the Cahn-Hillard model for crystal growth [44, 47, 38, 57, 50, 24, 15, 19, 18, 41, 43].

In this paper, we develop a set of linear, second order, unconditionally energy stable schemes using the energy quadratization (EQ) and scalar auxiliary variable (SAV) approach to solve the nonlocal Allen-Cahn equation numerically. The numerical schemes for the Allen-Cahn and the Cahn-Hilliard model are recalled in the paper simply for comparison purposes. In some of these schemes, both EQ and SAV methods are combined to yield linear, energy stable schemes. We note that when a nonlocal Allen-Cahn model is discretized, it is inevitable to yield an integral which has to be treated with a scalar auxiliary variable. When multiple integrals are identified as SAVs in the free energy functional of the Allen-Cahn model with nonlocal constraints, new solution procedures are developed to solve the subproblems in which elliptic equations can be solved efficiently. All these schemes are linear and second order accurate in time. On the other hand,when the EQ strategy is coupled with the discretized integrals, the Sherman-Morrison formula can lead us to an efficient numerical scheme as well. In fact, this can be equivalently dealt with using the SAV method, which will be discussed in the Appendix. The numerical schemes developed in this study for the Allen-Cahn equation with nonlocal constraints preserve not only mass but also the energy dissipation rate at the discrete level.

In the end, we conduct two numerical experiments to assess the performance of the schemes. The results based on EQ and those based on SAV methods perform equally well in preserving mass and the energy dissipation rate. In addition, the computational efficiency of the schemes is comparatively studied in one of the benchmark examples as well. To simplify the presentation, we present the temporal discretization of the models using EQ and SAV approaches in detail. Then, we only briefly discuss the strategy to obtain fully discrete schemes by discretizing the semi-discrete schemes in space later. We refer readers to our early publications in [21, 20] for details. We show that the linear systems resulting from the schemes are all solvable uniquely if the time step size is suitable so that the solution existence and uniqueness in the full-discrete system is warranted.

The rest of the paper is organized as follows. In $\S 2$, we present the mathematical models for the classical Allen-Cahn, the Cahn-Hilliard, and the Allen-Cahn model with nonlocal constraints. In $\S 3$, we compare their near equilibrium dynamics. In $\S 4$, we present a set of second order, linear, energy stable numerical schemes for the models. In $\S 5$, we conduct mesh refinement tests on all the schemes and carry out two simulations on crystal growth as well as its grain-boundary effects using the models. Finally, we give the concluding remark in section $\S 6$.

2 Phase Field Models for Crystal Growth

We consider a phase field model for modeling crystal growth in solids with a focus on resolving the detail of transient dynamics. The free energy of the phase field model for crystal growth is given by [11, 25, 19]

[TABLE]

where $\phi$ represents an atomistic density field, which is the deviation of the density from the average density and is a conserved field variable, $\epsilon$ is a parameter related to the temperature, that is, higher $\epsilon$ corresponds to a lower temperature, and $\nabla$ is the gradient operator ($\nabla^{2}$ denotes the Laplacian). In this study, we use a more general form of the free energy given by

[TABLE]

where $a=1,\alpha=1-\varepsilon$ recovers equation (2.2) and the boundary conditions

[TABLE]

are assumed.

Based on the Onsager linear response theory[31, 32], transient dynamics of such a system is customarily governed by a time-dependent partial differential equation given by

[TABLE]

subject to appropriate boundary and initial conditions, where $M$ and $\mu=\frac{\delta F}{\delta\phi}$ are the mobility matrix and the chemical potential, respectively.

The time rate of change of the free energy is given by

[TABLE]

The following boundary conditions will annihilate the boundary terms in the energy dissipation functional

[TABLE]

In the crystal growth model, $\frac{\partial f}{\partial\nabla\phi}=0$. So, the boundary conditions reduce to

[TABLE]

This set of boundary conditions is consistent with the derivation of the free energy functional.

An alternative set of boundary conditions is given by

[TABLE]

This is equivalent to

[TABLE]

This is different from the previous one.

We can also assign dissipative boundary conditions to the model as follows

[TABLE]

where $\beta>0$ is inversely proportional to a relaxation time. The boundary contribution to the energy dissipation is then given by

[TABLE]

The total energy dissipation rate or energy dissipation functional is given by

[TABLE]

If $\beta\to\infty$, we recover (2.15).

Two well-known phase field models are the Allen-Cahn and the Cahn-Hilliard equation, whose mobility is given respectively by

[TABLE]

where $M_{0}$ is a prescribed mobility coefficient matrix, which can be a function of $\phi$. The Allen-Cahn equation does not conserve the total mass $\int_{\Omega}\phi d\mathbf{r}$ if $\phi$ is the mass-fraction while the Cahn-Hilliard equation does. However, these two models predict similar near equilibrium dynamics. On the other hand the Allen-Cahn equation is an equation of lower spatial derivatives, and presumably costs less when solved numerically. Thus, one can impose the mass conservation as a constraint to the Allen-Cahn equation for it to be used to describe dynamics in which the mass is conserved. Next, we will briefly recall several ways to enforce mass conservation to dynamics described by the Allen-Cahn equation.

2.1 Allen Cahn model

The classical Allen-Cahn equation with the non-flux Neumann boundary conditions is given by

[TABLE]

where $M$ is the mobility coefficient and $\mu$ is the chemical potential

[TABLE]

The energy dissipation rate of the Allen-Cahn equation is given by

[TABLE]

provided nonnegative $M$. The Allen-Cahn model does not conserve the mass if the mass is denoted as $\int_{\Omega}\phi\mathrm{d}{\bf r}$.

2.2 Cahn-Hillard model

The Cahn-Hilliard equation with the non-flux Neumann boundary condition is given by

[TABLE]

where $M$ is the mobility coefficient and $\mu$ is the chemical potential

[TABLE]

The energy dissipation rate of the equation is given by

[TABLE]

provided nonnegative $M$. The Cahn-Hillard model conserves the mass. We next discuss the nonlocal Allen-Cahn equation that conserves the mass.

2.3 Nonlocal Allen-Cahn models

we present two methods to impose mass conservation. One is called the Allen-Cahn model with a penalizing potential and the other is called the Allen-Cahn model with a Lagrange multiplier.

2.3.1 Allen-Cahn model with a penalizing potential

In the Allen-Cahn model with a penalizing potential model, a penalizing term is augmented to the free energy to enforce the mass conservation by the model as follows

[TABLE]

where $\eta$ is penalizing parameter, ${\bf M_{0}}=\int_{\Omega}\phi(0)\mathrm{d{\bf r}}$ is the initial mass.

The transport equation for $\phi$ is given by the Allen-Cahn equation

[TABLE]

where $M$ is the mobility coefficient and $\tilde{\mu}$ is the chemical potential given by

[TABLE]

We calculate the energy dissipation rate as follows

[TABLE]

provided $M\geq 0$. The modified Allen-Cahn equation is approximately mass preserving depending on the size of $\eta$ and nonlocal. We next discuss another approaches to obtain a nonlocal Allen-Cahn model.

2.3.2 Allen-Cahn model with a Lagrange multiplier

In this model, the free energy is augmented by a penalty term with a Lagrange multiplier $L$ as follows.

[TABLE]

The transport equation for $\phi$ is given by the Allen-Cahn equation

[TABLE]

where $M$ is the mobility coefficient and $\tilde{\mu}$ is the chemical potential given by

[TABLE]

We calculate the energy dissipation rate as follows

[TABLE]

provided $M\geq 0$.

3 Near equilibrium dyamics of the models

We study dynamics of the models near their equilibrium solution $\phi^{ss}$. We consider a small perturbation of the steady state given by $\delta v(t,\mathbf{r})$:

[TABLE]

For the Allen-Cahn model, substituting (3.2) into (2.27), we get

[TABLE]

We seek the solution of the linearized partial differential equation system given by

[TABLE]

in the domain $\Omega=[-\pi,\pi]^{2}$. Then, we have

[TABLE]

Instability can emerge if $(k^{2}+l^{2})^{2}-2a(k^{2}+l^{2})+(\alpha+3(\phi^{ss})^{2})<0$, for some wave numbers $k$, $l$.

For the Allen-Cahn model with a penalizing potential, substituting (3.2) into the transport equation, we obtain the linearized system as follows

[TABLE]

Using the ansatz (3.5), we have

[TABLE]

If $(k^{2}+l^{2})^{2}-2a(k^{2}+l^{2})+\alpha+3(\phi^{ss})^{2}+4\pi^{2}\eta\delta_{k0}\delta_{l0}<0$, instability will occur. In comparison, the Allen-Cahn model with a penalizing potential is more stable than the classical Allen-Cahn model.

For the Allen-Cahn model with a Lagrange multiplier, substituting equation 3.2 into the transport equation, we get the linearized system

[TABLE]

where $g(\phi)=M((\nabla^{2}+a)^{2}\phi+(\alpha-a^{2})\phi+\phi^{3})$. Solving the linear system using ansatz (3.5), we have

[TABLE]

If $(k^{2}+l^{2})^{2}-2a(k^{2}+l^{2})+\alpha+3{(\phi^{ss})}^{2}-\frac{\int_{\Omega}(\alpha+3{(\phi^{ss})}^{2})\mathrm{d{\bf r}}}{\int_{\Omega}M\mathrm{d{\bf r}}}\delta_{k0}\delta_{l0}<0$, instability may ensure. The contribution of the Lagrange multiplier is to introduce a destabilizing mechanism depending on steady state solution $\phi^{ss}$.

For the Cahn-Hillard model, repeating the above analysis, we have the dynamical equation for the Fourier coefficients:

[TABLE]

The window of instability in the Cahn-Hillard model is identical to that in the Allen-Cahn model. However, the growth rates differ.

The linear stability results dictate initial transient dynamics of the solution towards or away from the given steady state. We will resort to numerical computations for long time transient behavior of the solution.

4 Numerical Approximations to the Phase Field Models

We design numerical schemes to solve the above nonlocal phase field equations to ensure that the energy dissipation property as well as mass conservation are respected. We do it by employing the energy quadratization (EQ) and the scalar auxiliary variable method (SAV) developed recently [54, 58, 39, 48]. Both methods depend on a reformulation of the models into equivalent ones with a quadratic energy. From the latter, we have effective ways to design linear numerical schemes. For a full review on EQ methods for thermodynamical models, readers are referred to a recent review article [57]. All schemes presented below are firstly given as semi-discretized ones in time and then the full discretization in space will be discussed. In fact, we have shown recently that BDF and Runge-Kutta methods can be used to design energy stable schemes for thermodynamical systems up to arbitrarily high order in time [53]. For comparison purposes, we also present analogous schemes for the classical Allen-Cahn and the Cahn Hilliard model as well.

4.1 Temporal discretization

4.1.1 Numerical methods for the Allen-Cahn model by EQ methods

We reformulate the free energy density by introducing an intermediate variable:

[TABLE]

Then, the free energy recast into

[TABLE]

We rewrite (2.27) as

[TABLE]

where

[TABLE]

We now discretize it using the linear Crank-Nicolson method in time to arrive at a second order semi-discrete scheme.

Scheme 4.1**.**

Given initial conditions $\phi^{0},q^{0}$, we first compute $\phi^{1},q^{1}$ by a first order scheme. Having computed $\phi^{n-1},q^{n-1}$, and $\phi^{n},q^{n}$, we compute $\phi^{n+1},q^{n+1}$ as follows.

[TABLE]

where

[TABLE]

The numerical implementation can be done as follows

[TABLE]

where $A=I+\Delta t\overline{M}^{n+1/2}[\frac{\nabla^{4}}{2}+a\nabla^{2}+\frac{\alpha}{2}+\frac{{(\overline{q^{\prime}}^{n+1/2})}^{2}}{4}],b^{n}=\phi^{n}-\Delta t\overline{M}^{n+1/2}[\frac{\nabla^{4}}{2}\phi^{n}+a\nabla^{2}\phi^{n}+\frac{\alpha}{2}\phi^{n}+\frac{1}{2}q^{n}\overline{q^{\prime}}^{n+1/2}-\frac{{(\overline{q^{\prime}}^{n+1/2})}^{2}}{4}\phi^{n}]$. So, $\phi^{n+1}$ is solved independent of $q^{n+1}$.

We define the discrete energy as follows

[TABLE]

4.1.2 Numerical methods for the Allen-Cahn model by SAV methods

Introducing an intermediate variable$r=\sqrt{\int_{\Omega}\frac{\phi^{4}}{4}\mathrm{d}{\bf r}+C_{0}}$ as the scalar auxiliary variable, the free energy recast into

[TABLE]

We rewrite (2.27) as

[TABLE]

We then discretize it using the linear Crank-Nicolson method in time to arrive at a second order semi-discrete scheme.

Scheme 4.2**.**

Given initial conditions $\phi^{0},r^{0}$, we first compute $\phi^{1},r^{1}$ by a first order scheme. Having computed $\phi^{n-1},r^{n-1}$, and $\phi^{n},r^{n}$, we compute $\phi^{n+1},r^{n+1}$ as follows.

[TABLE]

where

[TABLE]

We define the discrete energy as follows

[TABLE]

The numerical scheme can be recast into

[TABLE]

where $A=I+\Delta t\overline{M}^{n+1/2}[\frac{\nabla^{4}}{2}+a\nabla^{2}+\frac{\alpha}{2}],c=\overline{g}^{n+1/2},d=\Delta t\overline{M}^{n+1/2}\overline{g}^{n+1/2}$ and $b^{n}=\phi^{n}-\Delta t\overline{M}^{n+1/2}[\frac{\nabla^{4}}{2}\phi^{n}+a\nabla^{2}\phi^{n}+\frac{\alpha}{2}\phi^{n}+2\overline{g}^{n+1/2}r^{n}-\overline{g}^{n+1/2}\int_{\Omega}\overline{g}^{n+1/2}\phi^{n}\mathrm{d}{\bf r}]$. Multiplying the inverse of $A$ firstly and taking the inner product of the equation with c secondly, we have

[TABLE]

Then, the solution in the scheme is solved in the following steps,

[TABLE]

4.1.3 Numerical methods for the Cahn-Hillard model by EQ methods

The free energy density is reformulated by introducing an intermediate variable:

[TABLE]

Then, the free energy recast into

[TABLE]

We rewrite (2.35) as

[TABLE]

where

[TABLE]

We then discretize it using the linear Crank-Nicolson method in time to arrive at a second order semi-discrete scheme.

Scheme 4.3**.**

Given initial conditions $\phi^{0},q^{0}$, we first compute $\phi^{1},q^{1}$ by a first order scheme. Having computed $\phi^{n-1},q^{n-1}$, and $\phi^{n},q^{n}$, we compute $\phi^{n+1},q^{n+1}$ as follows.

[TABLE]

The discrete energy is defined as follows

[TABLE]

The numerical implementation can be done as follows

[TABLE]

where $A=I-\Delta t\nabla\cdot(\overline{M}^{n+1/2}\nabla[\frac{\nabla^{4}}{2}+a\nabla^{2}+\frac{\alpha}{2}+\frac{{(\overline{q^{\prime}}^{n+1/2})}^{2}}{4}]),b^{n}=\phi^{n}+\Delta t\nabla\cdot(\overline{M}^{n+1/2}\nabla[\frac{\nabla^{4}}{2}\phi^{n}+a\nabla^{2}\phi^{n}+\frac{\alpha}{2}\phi^{n}+\frac{1}{2}q^{n}\overline{q^{\prime}}^{n+1/2}-\frac{{(\overline{q^{\prime}}^{n+1/2})}^{2}}{4}\phi^{n}])$.

4.1.4 Numerical methods for the Cahn-Hillard model by SAV methods

An intermediate variable $r=\sqrt{\int_{\Omega}\frac{\phi^{4}}{4}\mathrm{d}{\bf r}+C_{0}}$ are introduced to reformulate the free energy, which is

[TABLE]

We rewrite (2.35) as

[TABLE]

Linear Crank-Nicolson method is used in time to arrive at a second order semi-discrete scheme.

Scheme 4.4**.**

Given initial conditions $\phi^{0},r^{0}$, we first compute $\phi^{1},r^{1}$ by a first order scheme. Having computed $\phi^{n-1},r^{n-1}$, and $\phi^{n},r^{n}$, we compute $\phi^{n+1},r^{n+1}$ as follows.

[TABLE]

where

[TABLE]

The discrete energy is defined as follows

[TABLE]

The numerical scheme can be rewritten into

[TABLE]

where $A=I-\Delta t\nabla\cdot(\overline{M}^{n+1/2}\nabla[\frac{\nabla^{4}}{2}+a\nabla^{2}+\frac{\alpha}{2}]),c=\overline{g}^{n+1/2},d=\Delta t\overline{M}^{n+1/2}\overline{g}^{n+1/2}$ and $b^{n}=\phi^{n}+\Delta t\nabla\cdot(\overline{M}^{n+1/2}\nabla[\frac{\nabla^{4}}{2}\phi^{n}+a\nabla^{2}\phi^{n}+\frac{\alpha}{2}\phi^{n}+2\overline{g}^{n+1/2}r^{n}-\overline{g}^{n+1/2}\int_{\Omega}\overline{g}^{n+1/2}\phi^{n}\mathrm{d}{\bf r}])$. Multiplying the first equation by the inverse of $A$ firstly and taking the inner product of the equation with c secondly, we have

[TABLE]

So we solve the solution in the scheme in the following steps,

[TABLE]

4.1.5 Numerical method for the Allen-Cahn model with a penalizing potential by EQ methods

In the Allen-Cahn model with a penalizing potential, we reformulate the free energy density by introducing two intermediate variables

[TABLE]

Then, the free energy recast into

[TABLE]

We rewrite the nonlocal Allen-Cahn equation as follows:

[TABLE]

where

[TABLE]

We then discretize it using the linear Crank-Nicolson method in time to arrive at a new scheme as follows.

Scheme 4.5**.**

Given initial conditions $\phi^{0},q^{0}$, we first compute $\phi^{1},q^{1}$ by a first order scheme. Having computed $\phi^{n-1},q^{n-1}$, and $\phi^{n},q^{n}$, we compute $\phi^{n+1},q^{n+1}$ as follows.

[TABLE]

where

[TABLE]

The discrete energy is defined as follows

[TABLE]

From the scheme, it follows that

[TABLE]

It can be written into a compact form,

[TABLE]

where $A=I+\Delta t\overline{M}^{n+1/2}[\frac{\nabla^{4}}{2}+a\nabla^{2}+\frac{\alpha}{2}+\frac{1}{2}{(\overline{q^{\prime}}^{n+1/2})}^{2}],c=1,d=\frac{\Delta t\overline{M}^{n+1/2}}{2}$ and $b^{n}=(I-\Delta t\overline{M}^{n+1/2}[\frac{\nabla^{4}}{2}+a\nabla^{2}+\frac{\alpha}{2}+q^{n}\overline{q^{\prime}}^{n+1/2}-\frac{1}{2}{(\overline{q^{\prime}}^{n+1/2})}^{2}])\phi^{n}-\sqrt{\eta}\zeta^{n}+\Delta t\overline{M}^{n+1/2}\frac{\eta}{2}\int_{\Omega}\phi^{n}\mathrm{d{\bf r}}.$

Then, the solve is solved in the following steps,

[TABLE]

4.1.6 Numerical method for the Allen-Cahn model with a penalizing potential by SAV methods

In the penalizing nonlocal Allen-Cahn model, we reformulate the free energy density by introducing two intermediate variables

[TABLE]

Then, the free energy recast into

[TABLE]

We rewrite the nonlocal Allen-Cahn equation as follows:

[TABLE]

where

[TABLE]

We then discretize it using the linear Crank-Nicolson method in time to arrive at a new scheme as follows.

Scheme 4.6**.**

Given initial conditions $\phi^{0},r^{0}$, we first compute $\phi^{1},r^{1}$ by a first order scheme. Having computed $\phi^{n-1},r^{n-1}$, and $\phi^{n},r^{n}$, we compute $\phi^{n+1},r^{n+1}$ as follows.

[TABLE]

where

[TABLE]

The discrete energy is defined as follows

[TABLE]

The scheme can be recast into

[TABLE]

where

[TABLE]

It implies that

[TABLE]

We solve for $(\phi^{n+1},c_{1})$ and $(\phi^{n+1},c_{2})$ from the above equation after we obtain

[TABLE]

So, the solution is solved in the following steps,

[TABLE]

4.1.7 Numerical method for the Allen-Cahn model with a Lagrange multiplier by EQ methods

We reformulate the free energy density by introducing an intermediate variable

[TABLE]

Then, the free energy recast into

[TABLE]

We rewrite (2.54) as

[TABLE]

where

[TABLE]

We then discretize it using the linear modified Crank-Nicolson method in time as follows.

Scheme 4.7**.**

Given initial conditions $\phi^{0},q^{0}$, we first compute $\phi^{1},q^{1}$ by a first order scheme. Having computed $\phi^{n-1},q^{n-1}$, and $\phi^{n},q^{n}$, we compute $\phi^{n+1},q^{n+1}$ as follows.

[TABLE]

where

[TABLE]

Remark: ${L}^{n+1/2}\neq\frac{L^{n}+L^{n+1}}{2}$.

Then, we have the following theorem

Theorem 4.1**.**

The mass of each phase is conserved, i.e.,

[TABLE]

Proof: Substituting the ${L}^{n+1/2}$ into the equation below, we have

[TABLE]

This implies the masse-conservation property.

We define the discrete energy as follows

[TABLE]

The solution is solved in the following steps

[TABLE]

where

[TABLE]

4.1.8 Numerical method for the Lagrangian models by SAV methods

We reformulate the free energy density by introducing an intermediate variable

[TABLE]

Then, the free energy recast into

[TABLE]

We rewrite (2.54) as

[TABLE]

where

[TABLE]

We then discretize it using the modified Crank-Nicolson method in time as follows.

Scheme 4.8**.**

Given initial conditions $\phi^{0},r^{0}$, we first compute $\phi^{1},r^{1}$ by a first order scheme. Having computed $\phi^{n-1},r^{n-1}$, and $\phi^{n},r^{n}$, we compute $\phi^{n+1},r^{n+1}$ as follows.

[TABLE]

where

[TABLE]

Then, we have the following theorem

Theorem 4.2**.**

The mass of each phase is conserved, i.e.,

[TABLE]

Proof: The proof is similar to that of theorem 4.1 and is thus omitted.

We define the discrete energy as follows

[TABLE]

This scheme can be recast into

[TABLE]

So we have

[TABLE]

We solve for $(\phi^{n+1},c_{1})$ and $(\phi^{n+1},c_{2})$ from the above equations after we obtain

[TABLE]

where

[TABLE]

Finally, the solution is obtained as follows

[TABLE]

4.2 Spatial discretization

We use the finite difference method to discretize the equations with the Neumann boundary condition in space for all models. The linear spatially dependent PDE systems resulting from all semi-discrete schemes are spatially discretized by compact second order finite difference methods at the cell center. We divide the 2D domain $\Omega=[0,L_{x}]\times[0,L_{y}]$ into rectangular meshes with mesh sizes $h_{x}=L_{x}/N_{x}$ and $h_{y}=L_{y}/N_{y}$, where $L_{x}$, $L_{y}$ are two positive real numbers and $N_{x}$, $N_{y}$ are the number of meshes in each direction. After this, the sets of the cell center points $C_{x}$ and $C_{y}$ for the uniform partition are defined as follows

[TABLE]

where $x_{i}=(i-\frac{1}{2})h_{x}$ and $y_{j}=(j-\frac{1}{2})h_{y}$.

We define the east-west-edge-to-center and center-to-east-west-edge difference operators $d_{x}$ and $D_{x}$ as follows

[TABLE]

Similarly, we define the north-south-edge-to-center and center-to-north-south-edge difference operators $d_{y}$ and $D_{y}$ as follows

[TABLE]

The fully discrete Laplacian and fourth order gradient operator are given by

[TABLE]

The discrete inner product is defined as follows

[TABLE]

where $f_{ij}$ and $g_{ij}$ are given at the cell center. In particular,

[TABLE]

With the notations, we replace the continuous differential operators in the semi-discrete schemes by the discrete difference operators to arrive at the corresponding fully discrete schemes. To save space, we will not enumerate them here.

4.3 Energy dissipation property and solvability of the linear systems resulting from the schemes

We summarize the energy dissipation law and unique solvability for all linear systems resulting from the semidiscrete schemes presented in this section into two theorems. The proofs of the energy dissipation property and solvability for the fully discrete schemes are similar, we only prove the theorems for the Allen-Cahn model with a Lagrangian multiplier discreized using EQ methods below and omit others for simplicity. Based on scheme 4.7 and the spatial discretization, the fully discrete scheme corresponding to scheme 4.7 is summarized below.

Scheme 4.9**.**

Given initial conditions $\phi^{0},q^{0}$, we first compute $\phi^{1},q^{1}$ by a first order scheme. Having computed $\phi^{n-1},q^{n-1}$, and $\phi^{n},q^{n}$, we compute $\phi^{n+1},q^{n+1}$ as follows.

[TABLE]

where

[TABLE]

We define the discrete energy as follows

[TABLE]

Theorem 4.3**.**

The fully discrete scheme obeys the following energy dissipation law

[TABLE]

Proof: Taking inner product of (4.222) with $-u^{n+1/2}$, we obtain

[TABLE]

Taking inner product of $\mu^{n+1/2}$ with $\frac{\phi^{n+1}-\phi^{n}}{\Delta t}$, we have

[TABLE]

Taking inner product of $q^{n+1}-q^{n}$ with $\frac{q^{n+1}+q^{n}}{\Delta t}$, we obtain

[TABLE]

Combining the above equations, we obtain

[TABLE]

Substituting the expression of ${L}^{n+1/2}$ into the equation, we have

[TABLE]

Remark: (i). This proof applies to the semi-discrete schemes as well. (ii). When the linear schemes of the nonlocal Allen-Cahn model involve integrals discretized by a composite Trapezoidal rule, efficient numerical methods can be devised to solve the resulting linear systems. Such methods are derived from the Sherman-Morisson formula (See Appendix).

Theorem 4.4**.**

The linear system resulting from the above fully discrete scheme admits a unique solution.

Proof: Note that the solution in scheme 4.222 is solved via the following steps

[TABLE]

where

[TABLE]

From the Sherman-Morrison formula, we notice that the solution uniqueness of $A\phi^{n+1}+(\phi^{n+1},c)d=b^{n}$ depends on the uniqueness of the corresponding linear system $A\phi^{n+1}=b^{n}$. Now we only need to prove the uniqueness of the solution for

[TABLE]

If $\Delta t$ is small enough, we have

[TABLE]

So, $A\phi=0$ has only zero solution.

Remark: One can prove the uniqueness of the solution for any time step size if ${\mu}^{n+1/2}=\nabla_{h}^{4}{\phi}^{n+1/2}+2a\nabla_{h}^{2}\overline{\phi}^{n+1/2}+\alpha{\phi}^{n+1/2}+\frac{1}{2}{q}^{n+1/2}\overline{q^{\prime}}^{n+1/2}$ in equation (4.222).

5 Numerical Results and Discussions

In this section, we conduct mesh refinement tests to validate the accuracy of the proposed schemes and then present some numerical examples to assess the schemes for the nonlocal Allen-Cahn models against those for the Cahn-Hilliard model. For convenience, we refer the numerical schemes designed by EQ methods for the Allen-Cahn model, the Cahn-Hilliard model, the Allen-Cahn model with a penalizing potential and the Allen-Cahn model with a Lagrangian multiplier as AC-EQ, CH-EQ, AC-P-EQ and AC-L-EQ, respectively. Similarly, we name the numerical schemes obtained using SAV approaches for the models as AC-SAV, CH-SAV, AC-P-SAV, AC-L-SAV, respectively. In the following, we set the constant in the free energy at $C_{0}=1\times 10^{4}$ in all computations.

5.1 Accuracy test

We confirm the convergence rates of the proposed schemes for the PFC models through mesh refinement tests. The computational domain is set as $\Omega=[0,1]^{2}$ . The model parameter values are chosen as $a=1,\varepsilon=0.1,M=1\times 10^{-3}$. We solve the equations with the initial condition given by

[TABLE]

We choose the space step size $h_{x}=h_{y}=\frac{1}{256}$. By taking a linear refinement path $\Delta t=\frac{0.05}{2^{k}}$, $k=0,1,\cdots,5$, we calculate the $L^{2}$ errors of the phase variable with adjacent k at $t=1$. The tables show the schemes are second order accurate in time numerically.

We also compare the computational efficiency of all schemes designed by EQ and SAV methods with $M=1\times 10^{-3},1$ and $1\times 10^{2}$ in table 3. The AC-P-EQ/SAV schemes perform the best among the schemes for nonlocal Allen-Cahn models. Besides this, AC-P-EQ/SAV schemes also perform better than CH-EQ/SAV schemes in most test cases. In fact, the accuracy of the Cahn-Hillard model relies on the mobility coefficient $M$ more sensitively than the Allen-Cahn models do. Hence, the accuracy of the schemes for the nonlocal Allen-Cahn models is better than that for the Cahn-Hillard model if $M$ is large. We will discuss it in more details next (see Figure 5.5).

Table 1 Mesh refinement tests for the proposed schemes using EQ methods.

Scheme AC-EQ CH-EQ AC-P-EQ AC-L-EQ

Coarse $\Delta t$ Fine $\Delta t$ $L^{2}$ error order $L^{2}$ error order $L^{2}$ error order $L^{2}$ error order

5.00E-02 2.5E-2 1.21E-06

1.23E-05

1.21E-06

1.21E-06

2.5E-2 1.25E-2 3.02E-07 2.00 3.12E-06 1.98 3.02E-07 2.00 3.02E-07 2.00

1.25E-2 6.25E-3 7.56E-08 2.00 7.84E-07 1.99 7.56E-08 2.00 7.56E-08 2.00

6.25E-3 3.125E-3 1.89E-08 2.00 1.96E-07 2.00 1.89E-08 2.00 1.89E-08 2.00

3.125E-3 1.5625E-3 4.72E-09 2.00 4.90E-08 2.00 4.72E-09 2.00 4.72E-09 2.00

Table 2 Mesh refinement tests for the proposed schemes using SAV methods.

Scheme AC-SAV CH-SAV AC-P-SAV AC-L-SAV

Coarse $\Delta t$ Fine $\Delta t$ $L^{2}$ error order $L^{2}$ error order $L^{2}$ error order $L^{2}$ error order

5.00E-02 2.5E-2 1.24E-06

1.22E-05

1.24E-06

1.24E-06

2.5E-2 1.25E-2 3.09E-07 2.00 3.10E-06 1.98 3.09E-07 2.00 3.09E-07 2.00

1.25E-2 6.25E-3 7.72E-08 2.00 7.78E-07 1.99 7.72E-08 2.00 7.72E-08 2.00

6.25E-3 3.125E-3 1.93E-08 2.00 1.95E-07 2.00 1.93E-08 2.00 1.93E-08 2.00

3.125E-3 1.5625E-3 4.83E-09 2.00 4.87E-08 2.00 4.82E-09 2.00 4.82E-09 2.00

Table 3 Efficiency of the schemes at $1000$ time steps with respect to $M=0.001,1$ and $100$ (From top to below).

Scheme AC-EQ AC-SAV CH-EQ CH-SAV AC-P-EQ AC-P-SAV AC-L-EQ AC-L-SAV

Time (s) 45.3 77.9 66.0 107.8 68.3 90.0 110.0 103.2

Time (s) 44.4 76 54.6 104 38.8 84.2 53.9 100.7

Time (s) 53.4 75.0 66.8 111.9 61.4 83.2 103.3 99.5

5.2 Assessment of the numerical schemes

To further assess the numerical schemes, we numerically solve the model equations using the schemes with respect to two benchmark problems. Firstly we simulate the phase transition of crystal growth in 2D. We use time step $\Delta t=1\times 10^{-3}$ and $256\times 256$ space mesh in the 2D simulation. A solid crystallite with Hexagonal ordering in 2D is initially placed in the centre of the domain, which is assigned an average density $\overline{\phi}$. The initial condition is given by

[TABLE]

where

[TABLE]

[TABLE]

$\bf r_{0}$ is the center coordinate of the domain, and $\bf d_{0}$ is $\frac{1}{6}$ of the domain length in the x-direction. The domain is given by $\Omega=[0,\frac{2\pi}{q}a]\times[0,\frac{\sqrt{3}\pi}{q}b]$, $a=10$ and $b=12$. The other values are $\varepsilon=0.325,\overline{\phi}=\frac{\sqrt{\varepsilon}}{2},A=\frac{4}{5}(\overline{\phi}+\frac{\sqrt{15\varepsilon-36\overline{\phi}^{2}}}{3})$ and $q=\frac{\sqrt{3}}{2}.$

Figure 5.1 shows time evolution of the crystal growth process computed by AC-EQ, CH-EQ, AC-P-EQ and AC-L-EQ schemes, respectively. In Figure 5.1-(a), the crystal growth simulated by the Allen-Cahn model can’t preserve the Hexagonal ordering, different from the results simulated by the Cahn-Hillard model in Figure 5.1-(b) and the nonlocal Allen-Cahn models in Figure 5.1-(c) and (d). The time evolution of mass and free energy are shown in Figure 5.2-(a) and (b) respectively. The results computed by the EQ and SAV schemes for the same model (Allen-Cahn model, Cahn-Hillard model and nonlocal Allen-Cahn model) are identical. We don’t see any differences between the results of the Allen-Cahn model with a penalizing potential and the Allen-Cahn model with a Lagrangian multiplier either. The mass decays in the Allen-Cahn model to nearly zero in finite time. In contrast, the mass in the Cahn-Hillard model and the nonlocal Allen-Cahn models is conserved in the simulations. Meantime, the free energies of the Cahn-Hillard model and the nonlocal Allen-Cahn models are larger than that of the Allen-Cahn model. Figure 5.2-(b) shows that the free energy computed by the nonlocal Allen-Cahn model reaches the steady state faster than that of the Cahn-Hillard model.

Secondly we simulate another case of polycrystalline growth involving the grain boundary effect, where the two initial crystallites with a hexagonal configuration oriented in different direction (or misorientation) are put in the domain. Grain boundaries appear when the two crystals meet during the growth, which yields some orientation mismatch. The initial condition is given by

[TABLE]

where

[TABLE]

${\bf r_{1}}=\sqrt{(x-\frac{1}{2}x_{0})^{2}+(y-\frac{1}{2}y_{0})^{2}}$, ${\bf r_{2}}=\sqrt{(x-\frac{3}{2}x_{0})^{2}+(y-\frac{3}{2}y_{0})^{2}}$ and $(x_{0},y_{0})$ is the center of the domain. The other parameters and the domain are the same as in the first example. By doing an affine transformation of the Cartesian coordinates $(x,y)$ to produce a rotation $\theta$ in the domain, the modified coordinates $(x_{\theta},y_{\theta})$ can be used to generate the crystallites in different directions,

[TABLE]

We put two crystallites in the domain, the first one is defined as equation (5.15) with $\theta=0$, the other is with $\theta=\frac{\pi}{8}.$

Figure 5.3 depicts the grain boundary effect during polycrystalline growth computed by AC-SAV, CH-SAV, AC-P-SAV and AC-L-SAV schemes, respectively. The snapshots of the phase transitions computed by the Cahn-Hillard and the nonlocal Allen-Cahn models show the Hexagonal ordering are broken at the center at $T=400$. The time evolution of mass and free energy are shown in Figure 5.4-(a) and (b), respectively. The mass in the Allen-Cahn model decays in time. In contrast, the mass in the Cahn-Hillard model and the nonlocal Allen-Cahn models are conserved during the simulation. Meanwhile, the free energies in the Cahn-Hillard model and the nonlocal Allen-Cahn models are larger than that of the Allen-Cahn model. In Figure 5.4-(b), the nonlocal Allen-Cahn models and the Cahn-Hillard model predict comparable time evolution of the free energy but the nonlocal Allen-Cahn models reach the steady state first. The results computed by the EQ and SAV schemes for the same model (the Allen-Cahn model, Cahn-Hillard model and nonlocal Allen-Cahn model) are nearly identical. We don’t see any differences between the results of the Allen-Cahn model with a penalizing potential and the Allen-Cahn model with a Lagrangian multiplier either.

We note that the results computed by the Allen-Cahn model with a penalizing potential and the Allen-Cahn model with a Lagrangian multiplier in the above two examples at $\eta=1\times 10^{3}$ are the same. However, the choice of $\eta$ can certainly affect the outcome. In principle, $\eta$ should be chosen as large as possible. However, when $\eta$ is too large, the governing equation becomes very stiff, which forces one to use extremely small time-step size in order to resolve the detail correctly. On the other hand, as we have shown in the two examples, $\eta=1\times 10^{3}$ is good enough to produce the results that conserve mass very well.

Both Figure 5.2 and Figure 5.4 show that the free energy computed by the nonlocal Allen-Cahn models reach the steady state faster than that of the Cahn-Hillard model, which is different from our previous results on a different free energy functional [26]. Especially, the simulations computed by the models with mobility coefficient $M=1$ and time steps $1\times 10^{-2},1\times 10^{-4}$ show the same time evolution behavior. When the mobility is large, say $M=10$, and $\Delta t=0.01$, the CH-SAV produce an erroneous result while the others produce comparable ones. In practice, if the steady state is more important than the transition dynamics, one can enlarge the mobility coefficient of the models to accelerate the convergence to steady state. But, CH-SAV seems to have some accuracy issues with this approach. The results of the nonlocal Allen-Cahn models in 5.5-(a) and (b) show better performance than the Cahn-Hillard model dose in Figure 5.5-(c) and (d).

6 Conclusions

We have developed a set of linear, second order, energy stable schemes for the Allen-Cahn equation with nonlocal constraints that conserve mass and compared them with the energy stable, linear schemes for the Allen-Cahn and the Cahn-Hilliard model. These schemes are devised based on the energy quadratization strategy in the form of EQ and SAV formulation, respectively. We show that they are unconditionally energy stable and uniquely solvable. All schemes can be solved using efficient numerical methods, making the models alternatives to the Cahn-Hilliard model to describe interface dynamics of immiscible materials while conserving mass. The nonlocal Allen-Cahn models show a faster coarsening rate than that of the Cahn-Hilliard model at the same mobility, but one can enlarge the mobility coefficient of the nonlocal Allen-Cahn model to accelerate their dynamics in case only steady states are of interest. In addition, we have compared the two Allen-Cahn models with nonlocal constraints numerically. The computational efficiency of the Allen-Cahn model with a penalizing potential is slightly better than the one with a Lagrange multiplier, but the accuracy of the former depends on a suitable choice of model parameter $\eta$. If the steady state is desired rather than the transient dynamical behavior, large time step size and mobility coefficient can be used to accelerate the simulation. In the end, we show that the nonlocal Allen-Cahn models perform better that the Cahn-Hilliard model in the case of a large time step and mobility coefficient.

Acknowledgements

Xiaobo Jing and Qi Wang’s research is partially supported by NSFC awards #11571032, #91630207 and NSAF-U1530401.

Appendix

Appendix A Shermann-Morrison formula and its application to solving integro-differential equations

Here we give a brief review over the Sherman-Morrison formula [5] and explain its applications in the practical implementation of our various relevant schemes.

Suppose $A$ is an invertible square matrix, and $u$,$v$ are column vectors. Then $A+uv^{T}$ is invertible iff $1+v^{T}A^{-1}u\neq 0$. If $A+uv^{T}$ is invertible, then its inverse is given by

[TABLE]

So if $Ay=b$ and $Az=u$, $(A+uv^{T})x=b$ has the solution given by

[TABLE]

For the integral term(s) in the semi-discrete schemes in this study such as (4.150), we need to discretize it properly. $\forall f$, we discretize $\int_{\Omega}f\mathrm{d{\bf r}}$ using the composite trapezoidal rule and adding all the elements of the new matrix $w_{1}w_{2}^{T}f$, where $w_{1}=\frac{h_{x}}{2}S$, $w_{2}=\frac{h_{y}}{2}S$, $h_{x}$, $h_{y}$ are the spatial step sizes and $S={[1,2,2,...,2,2,1]}^{T}$. For convenience, we use $w_{1}w_{2}^{T}f$ to represent the integral discretized by the composite trapezoidal rule.

To solve equation (4.150), we discretize the integral or the inner product of functions $(c,\phi^{n+1})d$ as $u{v}^{T}\phi^{n+1}$. The scheme is recast to $A\phi^{n+1}+u{v}^{T}\phi^{n+1}=b^{n}$. After using the Sherman-Morrison formula, we get

[TABLE]

In the inner product of vectors, (4.150) can be rewritten into

[TABLE]

So, indeed the approach we take in the study using the discrete inner product is essentially equivalent to applying the Sherman-Morrison formula.