Energy stable discretization of Allen-Cahn type problems modeling the motion of phase boundaries

TL;DR

This paper develops an energy stable numerical scheme for Allen-Cahn type models of phase boundary motion, ensuring energy decay at discrete levels and demonstrating effectiveness through theoretical proofs and numerical tests.

Contribution

It introduces a conforming Galerkin spatial discretization and an implicit time-stepping scheme that guarantee energy decay and well-posedness for Allen-Cahn type problems.

Findings

Proven energy decay for semi-discrete and fully discrete schemes

Numerical tests confirm theoretical energy decay and interface behavior

Scheme effectively captures asymptotic interface velocities

Abstract

We study the systematic numerical approximation of a class of Allen-Cahn type problems modeling the motion of phase interfaces. The common feature of these models is an underlying gradient flow structure which gives rise to a decay of an associated energy functional along solution trajectories. We first study the discretization in space by a conforming Galerkin approximation of a variational principle which characterizes smooth solutions of the problem. Well-posedness of the resulting semi-discretization is established and the energy decay along discrete solution trajectories is proven. A problem adapted implicit time-stepping scheme is then proposed and we establish its well-posed and decay of the free energy for the fully discrete scheme. Some details about the numerical realization by finite elements are discussed, in particular the iterative solution of the nonlinear problems arising in every time-step. The theoretical results are illustrated by numerical tests which also provide further evidence for asymptotic expansions of the interface velocities derived by Alber et al.

Figures

Equations

$${\partial }_t{S}^{ac}=\left({\mu }^{1/2}\lambda {\Delta }_x{S}^{ac}-\frac{1}{{\mu }^{1/2}}{\psi }^{\prime }(S^{ac})-F^{\prime }(S^{ac})\right)\frac{c}{(\mu \lambda {)}^{1/2}}.$$

$${\partial }_{t}S=\Delta S-\frac{1}{{\epsilon }^2}{\psi }^{\prime }(S)$$

$$E^{ac}(S)={\int }_{\Omega }\frac{{\mu }^{1/2}\lambda }{2}|{\nabla }_{x}S{|}^2+\frac{1}{{\mu }^{1/2}}\psi (S)+F(S)dx.$$

$$v^{ac}(\mu ;\lambda )=\frac{c}{c_1}{n}_{\Gamma }\cdot [C]\cdot n_{\Gamma }+c{\lambda }^{1/2}{\kappa }_{\Gamma }+O({\mu }^{1/2}),$$

$$\bar{v}(\lambda )=\frac{c}{c_1}{n}_{\Gamma }\cdot [C]\cdot n_{\Gamma }+c{\lambda }^{1/2}{\kappa }_{\Gamma }.$$

$${\partial }_t{S}^{hyb}=\left(\nu {\Delta }_x{S}^{hyb}-{\psi }^{\prime }(S^{hyb})-F^{\prime }(S^{hyb})\right)|\nabla S^{hyb}|.$$

$$E^{hyb}(S)={\int }_{\Omega }\frac{\nu }{2}|{\nabla }_{x}S{|}^2+\psi (S)+F(S)dx.$$

$$v^{hyb}(\nu )=n\cdot [C]\cdot n+{\nu }^{1/2}{\omega }_1{\kappa }_{\Gamma }+o({\nu }^{1/2}),$$

$$c(|\nabla S|){\partial }_{t}S$$

$$\alpha {\partial }_{n}S$$

$$E(S)={\int }_{\Omega }\frac{\alpha }{2}|\nabla S{|}^2+\beta d(S)dx.$$

$${\partial }_{t}S=-{\text{grad}}_{g_S}E(S),$$

$${\int }_{\Omega }c(|\nabla S|){\partial }_{t}Svdx=:g_S({\partial }_{t}S,v)=-⟨\delta E(S),v⟩:=-{\int }_{\Omega }\alpha \nabla S\nabla v+\beta d^{\prime }(S)vdx$$

$$\left(c(|\nabla S(t)|){\partial }_{t}S(t),v\right)+(\alpha \nabla S(t),\nabla v)+(\beta d^{\prime }(S(t)),v)=0$$

$$\frac{d}{dt}E(S(t))=-{\int }_{\Omega }c(|\nabla S(t)|)|{\partial }_{t}S(t){|}^{2}dx\le 0.$$

$$\frac{d}{dt}E(S(t))$$

$$(\ast )$$

$$=-(c(|\nabla S(t)|){\partial }_{t}S(t),{\partial }_{t}S(t))=-{\int }_{\Omega }c(|\nabla S(t)|)|{\partial }_{t}S(t){|}^{2}dx\le 0.$$

$$\displaystyle\Big{(}c(|\nabla S_{h}(t)|)\partial_{t}S_{h}(t),v_{h}\Big{)}+(\alpha\nabla S_{h}(t),\nabla v_{h})+(\beta d^{\prime}(S_{h}(t)),v_{h})=0$$

$$\frac{d}{dt}E(S_h(t))=-{\int }_{\Omega }c(|\nabla S_h(t)|)|{\partial }_t{S}_h(t){|}^{2}dx\le 0.$$

$$M(y)y^{\prime }(t)+Ay+B(y)=0,y(0)=y_0,$$

$$C^{\prime }:=E(S_{0,h})\ge E(S_h(t))=\frac{\alpha }{2}\Vert \nabla S_h(t){\Vert }_{L^2(\Omega )}^2+\beta \Vert d(S_h(t)){\Vert }_{L^1(\Omega )}.$$

$$C^{\prime }\ge \beta \Vert d(S_h(t)){\Vert }_{L^1(\Omega )}\ge \beta \underline{d}\Vert S_h(t){\Vert }_{L^1(\Omega )}.$$

$${\bar{\partial }}_{\tau }{S}^n=\frac{1}{\tau }(S^n-S^{n-1})$$

$$\displaystyle\Big{(}c(|\nabla S_{h}^{n-1}|)\bar{\partial}_{\tau}S_{h}^{n},v_{h}\Big{)}+(\alpha\nabla S_{h}^{n},\nabla v_{h})+\Big{(}\beta\frac{d(S_{h}^{n})-d(S_{h}^{n-1})}{S_{h}^{n}-S_{h}^{n-1}},v_{h}\Big{)}=0$$

$${\bar{\partial }}_{\tau }E(S_h^n)$$

$$=-{\int }_{\Omega }c(|\nabla S_h^{n-1}|)|{\bar{\partial }}_{\tau }{S}_h^n{|}^{2}dx-\frac{\alpha }{2\tau }\Vert \nabla S_h^n-\nabla S_h^{n-1}{\Vert }_{L^2(\Omega )}^2\le 0.$$

$$\displaystyle E(S_{h}^{n})-E(S_{h}^{n-1})=\frac{\alpha}{2}\Big{(}(\nabla S_{h}^{n},\nabla S_{h}^{n})-(\nabla S_{h}^{n-1},\nabla S_{h}^{n-1})\Big{)}+\beta\Big{(}d(S_{h}^{n})-d(S_{h}^{n-1}),1\Big{)}$$

$$=\left(\alpha \nabla S_h^n,\nabla S_h^n-\nabla S_h^{n-1}\right)-\frac{\alpha }{2}\Vert \nabla S_h^n-\nabla S_h^{n-1}{\Vert }^2+\left(\beta \frac{d(S_h^n)-d(S_h^{n-1})}{S_h^n-S_h^{n-1}},S_h^n-S_h^{n-1}\right).$$

$$\frac{\underline{c}}{\tau }\Vert S_h^n-{\tilde{S}}_h^n{\Vert }_{L^2(\Omega )}^2$$

$$\le \beta \Vert D(S_h^n,S_h^{n-1})-D({\tilde{S}}_h^n,S_h^{n-1}){\Vert }_{L^2(\Omega )}\Vert S_h^n-{\tilde{S}}_h^n{\Vert }_{L^2(\Omega )}.$$

$$D(S,S^{n-1})=\frac{d(S)-d(S^{n-1})}{S-S^{n-1}}$$

Full text

Keywords: Allen-Cahn equation, phase-field models, gradient systems, mean curvature flow, finite elements, implicit time stepping

AMS-classification (2000): 65M60,74N20,35J93,53C44

1. Introduction

We consider the systematic numerical approximation of a class of Allen-Cahn type equations describing, for instance, the motion of anti-phase boundaries in crystalline solids [3], the evolution of interfaces in more general phase field models [6, 23, 36], or the geometric evolution by mean curvature [12, 24]. For motivation of our considerations, let us briefly consider two particular examples which have been studied in detail in [1, 2]. As a first model problem, we consider the generalized Allen-Cahn equation

[TABLE]

Here $S^{ac}$ is an order parameter which takes values close to zero in one phase of the medium and close to one in the other, $\psi(S)=4S^{2}(1-S)^{2}$ is the double well potential, $F^{\prime}(S)$ is a source term due to external loading, and $c$, $\lambda$ are model parameters. The coefficient $\mu$ plays the role of a regularization parameter that allows to widen the transition zone between the phases and thus to alleviate the numerical simulation. For the choice $c=1$, $\lambda=1$, $\mu=\varepsilon^{2}$, and $F^{\prime}(S)=0$, one obtains the standard Allen-Cahn equation

[TABLE]

which has been studied intensively in the literature; see e.g. [18, 26, 33] and the references given there. It is well-known [25] that (1.2), and therefore also (1.1), can be understood as a gradient flow for an associated free energy functional, which for the problem (1.1) has the form

[TABLE]

The multiplicative term $\frac{c}{(\mu\lambda)^{1/2}}$ on the right hand side of (1.1) yields a scaling of time that leads to a non-trivial behavior for $\mu\to 0$, which is called the sharp interface limit. It is shown in [1, 2] that the velocity of the interface $\Gamma(t)=\{x:S^{ac}(x)=1/2\}$ can be expanded as

[TABLE]

where $c_{1}$ can be expressed expicitly in terms of the problem data. Here $[C]$ denotes the jump of the Eshelby tensor [16, 17] which is related to the external forces $F^{\prime}(S)$ and $n_{\Gamma}$ is the normal vector while $\kappa_{\Gamma}$ is the mean curvature of the interface $\Gamma$. In the limit $\mu\to 0$, the velocity of the interface motion can thus be described by the well-known relation [26, 33]

[TABLE]

Choosing the regularization parameter $\mu>0$ leads to an incorrect speed of the interface motion, and according to (1.4), one should choose $\mu=o(\lambda)$ to keep the consistency error small. For parameter $\mu>0$ the interface is widened to a transition zone of width $O((\mu\lambda)^{1/2})$. In order to resolve the smoothed interface in simulations, the mesh size should therefore be chosen at least as small as $h=O((\mu\lambda)^{1/2})=o(\lambda)$. More restrictive assumptions on the mesh size are required to really prove a good approximation of the interface evolution; see [19, 30] for details. If $\lambda$ is small, then a very small mesh size is thus required to obtain a good approximation.

In order to allow for a larger mesh size and to reduce the computational burden, the following model for the interface motion has been proposed and analyzed in [1, 2]

[TABLE]

This equation, which has been termed hybrid model in [1, 2], describes the motion of the levelsets of the order parameter $S^{hyb}$ with normal speed proportional to the term in parenthesis. Following the arguments of [9, 13, 31], the problem (1.6) can again be interpreted as a gradient flow for an associated energy functional, which here reads

[TABLE]

The gradient here has to be defined with respect to a metric induced by the solution $S$ itselve; we refer to [9, 13] and to Section 2 for details. Using asymptotic analysis, the propagation speed of the interface has been shown in [1, 2] to behave like

[TABLE]

where the parameter $\omega_{1}$ can again be computed from the model data.

A comparison of the two formulas (1.4) and (1.8) reveals that for a proper choice of the model and scaling parameters, the hybrid model leads to the same propagation speed for the interface $\Gamma$ as in the Allen-Cahn model up to higher order terms. To obtain a good agreement of the results, one should choose $\nu\approx\lambda$ and in this case width of the smoothed interface in the hybrid model (1.6) can be shown to behave like $O(\nu^{1/2})=O(\lambda^{1/2})$; see [1, 2] for details. When $\lambda$ is small, the smoothed interface of the hybrid model is thus substantially larger then that of the corresponding Allen-Cahn model and a mesh size of $h=O(\lambda^{1/2})$ should be sufficient to resolve the interface in computations, instead of $h=o(\lambda)$ required for the corresponding Allen-Cahn model. The hybrid model should therefore allow to compute the interface evolution on a much coarser mesh, which was the main motivation for the proposal of this method.

Due to the close connection to mean curvature flow, there has been much interest in the numerical approximation for the Allen-Cahn equation (1.2). Various discretization schemes have been proposed and analyzed, see e.g. [7, 35] and the references given there. Particular aspects of the efficient numerical approximation have been addressed in [20, 21, 34] and a-posteriori error estimates have been derived in [5, 22, 28]. The convergence of numerical approximations of the Allen-Cahn equation in the sharp interface limit $\varepsilon\to 0$ towards the generalized motion by mean curvature has been established in [19, 30]; see also [38] for further numerical tests. The direct numerical approximation of mean curvature flow has also been investigated intensively; see for instance [11, 37] and refer to [12, 25] for a review about available results and further references.

The focus of the current manuscript is somewhat different to these previous works: Here we want to study the systematic numerical approximation of a large class of Allen-Cahn type models which includes (1.1), (1.2), and (1.6) as special cases. The common feature of all these models is the underlying gradient flow structure with respect to an associated energy functional, and we aim to strictly preserve this structure in the approximation process on the semi-discrete and the fully discrete level. With this in mind, we automatically obtain discretization schemes that are uniformly energy stable and thus consistent with the second law of thermodynamics.

In this paper we establish well-posedness and qualitative properties of a discretization strategy appropriate for a large class of Allen-Cahn type problems. The models (1.1) and (1.6) introduced above will be considered as particular examples. As a by-product of our investigations, we also obtain further numerical evidence for the asymptotic expansions (1.4) and (1.8) of the interface velocities. Our main arguments can be extended to elastic Allen-Cahn type models which describe the evolution of phase interfaces in binary alloys [2, 23, 24]. The required coupling to the additional elasticity system will be investigated in a forth-comming publication.

The remainder of the manuscript is organized as follows: In Section 2, we first introduce a unified formulation of the problems discussed in the introduction and then present a variational characterization of smooth solutions. In addition, we verify that the respective energies decrease along solution trajectories which is due to the underlying gradient flow structure. In Section 3, we then discuss the Galerkin approximation of the underlying variational principles. We establish the well-posedness of the resulting semi-discrete schemes and verify the energy decay on the discrete level. In Section 4, we then discuss the subsequent discretization in time by a problem adapted implicit time stepping scheme. We establish the well-posedness of the fully discrete problem and again prove the decay in energy. A particular space discretization by finite elements is discussed in Seciton 5, and some further aspects of the implementation are presented. For illustration of our theoretical results, we report in Section 6 about some numerical tests,in which we also verify the asymptotic expansions for the interface velocities given in (1.4) and (1.8). We close with a short discussion of our results and open problems for future research.

2. A unified variational framework

For the rest of the manuscript, we consider the following general model problem:

[TABLE]

Here $\Omega\subset\mathbb{R}^{d}$, $1\leq d\leq 3$ is a bounded Lipschitz domain, $\alpha,\beta>0$ are model parameters, and the functions $c(\cdot)$, $d(\cdot)$ are given a-priori. For our analysis later on, we will assume that

• (A1)

$c\in W^{1,\infty}(\mathbb{R})$ and $0<\underline{c}\leq c(x)\leq\overline{c}<\infty$;

• (A2)

$d\in W^{2,\infty}(\mathbb{R})$ with $d(x)\geq 0$, $d(x)\geq\underline{d}|x|-d_{1}$, and $|d^{\prime\prime}(x)|\leq\overline{d}$ for all $x\in\mathbb{R}$.

Some of these technical assumptions are made only for convenience and could further be relaxed without problems. This will become clear from our analysis below and we will make some remarks in this direction. To complete the model description, we additionally assume that

[TABLE]

Other boundary conditions can also be treated without much difficulty. The above problem is of quasilinear parabolic type and existence of solutions for appropriate initial values $S(0)=S_{0}$ can be proven with standard arguments; see e.g. [4, 29, 32].

Notation 2.1**.**

A function $S:C^{1,2}([0,T]\times\overline{\Omega})$ is called smooth solution of (2.1)–(2.2), if it satisfies the respective equations in a pointwise sense.

Here $C^{1,2}([0,T]\times\overline{\Omega})$ denotes the space of continuous functions that are continuously differentiable with respect to $t$ and twice continuously differentiable with respect to $x$. Note that acording to our definition, smooth solutions and their spatial derivatives are uniformly bounded.

Remark 2.2**.**

The Allen-Cahn model (1.1) can be phrased in the form (2.1) by setting $\alpha=\mu^{1/2}\lambda$, $\beta=\mu^{-1/2}$, $c(x)=(\mu\lambda)^{1/2}/c$, and defining the potential as $d(S)=\psi(S)+\mu^{1/2}F(S)$. Note that the second derivative $d^{\prime\prime}(\cdot)$ of the potential is unbounded here. Using maximum principles, the solution $S$ can however be shown to remain bounded and therefore $d(S)$ could be modified for large $S$ in order to satisfy (A2) without changing the solution; see [19] for other appropriate conditions on the potential. The hybrid model (1.6) can also be cast into the form (2.1) with $\alpha=\nu$, $\beta=1$, $d(s)=\psi(s)+F(S)$, and $c(x)=1/x$. To avoid problems for vanishing gradient, one can set $c(x)=1/\max(\delta,\min(1/\delta,x))$ for some $\delta>0$ instead, and thus obtain a regularized version of (1.6). The limit $\delta\to 0$ may be studied in the framework of $\Gamma$-convergence [10, 27].

As a next step, let us define an associated energy functional for the problem (2.1)–(2.2) by

[TABLE]

Equation (2.1) can then be interpreted as a gradient flow for this energy with respect to the metric $g_{S}(v,w)=\int_{\Omega}vw\;c(|\nabla S|)dx$. Problem (2.1) can thus be formally expressed as

[TABLE]

which, following [9, 13, 25], has to be understood in the sense that

[TABLE]

for all $v\in C_{0}^{\infty}(\Omega)$ and $t>0$. The last term is the negative directional derivative of the energy functional in direction $v$. These formal arguments motivate the following variational principle.

Lemma 2.3** (Variational characterization).**

Let $S$ be a smooth solution of (2.1)–(2.2). Then

[TABLE]

for all $v\in H^{1}(\Omega)$ and all $t>0$. Here $(v,w)=\int_{\Omega}vw\;dx$ is the scalar product of $L^{2}(\Omega)$.

As a direct consequence of the gradient flow structure, we obtain

Lemma 2.4** (Energy decay).**

Let $S$ be a smooth solution of (2.1)–(2.2). Then

[TABLE]

Proof.

Since $S$ is smooth, we can apply the chain rule to get

[TABLE]

By Lemma 2.3, any smooth solution of (2.1)–(2.2) also satisfies the variational principle (2.6), and testing this variational equation with $v=\partial_{t}S(t)$ immediately yields

[TABLE]

For the last step we only had to use that $c(\cdot)$ is non-negative. ∎

3. Galerkin approximation in space

Motivated by the results of the previous section, we now consider the Galerkin approximation of the variational principle (2.6) to construct semi-discretizations in space. Let us choose some finite dimensional subspace $V_{h}\subset H^{1}(\Omega)$ and consider the following discrete variational problem.

Problem 3.1** (Galerkin semi-discretization).**

Let $S_{h,0}\in V_{h}$ and $T>0$ be given. Find $S_{h}\in C^{1}([0,T];V_{h})$ such that $S_{h}(0)=S_{h,0}$ and

[TABLE]

for all test functions $v_{h}\in V_{h}$ and for all $0\leq t\leq T$.

Since we employ a conforming Galerkin approximation, the gradient flow structure and thus also the energy decay estimates are inherited automatically by the semi-discrete problem.

Lemma 3.2**.**

Let $S_{h}$ denote a solution of Problem 3.1. Then

[TABLE]

Proof.

The assertion follows with literally the same arguments as in Lemma 2.4. ∎

As a next step, let us establish the well-posedness of the semi-discrete problem.

Lemma 3.3** (Well-posedness).**

For any $S_{h,0}\in V_{h}$ and $T>0$, Problem 3.1 has a unique solution.

Proof.

After choosing a basis for $V_{h}$, we can rewrite the discrete problem as initial value problem

[TABLE]

where $y(t)$ here denotes the coordinate vector when expanding $S_{h}(t)$ in the chosen basis. Due to assumptions (A1)–(A2) on the coefficients and since $V_{h}$ is finite dimensional, the matrix functions $M(\cdot)$ and $B(\cdot)$ can be seen to be Lipschitz continuous. From the lower bounds on $c(\cdot)$ we further obtain that $M(y)$ is symmetric positive definite. Existence of a unique local solution then follows from the Picard-Lindelöf theorem. From (2.3) and (3.2), we can see that

[TABLE]

Here we used that $d(\cdot)\geq 0$ due to assumption (A2). This yields $\alpha\|\nabla S_{h}(t)\|_{L^{2}\Omega)}^{2}\leq 2C^{\prime}$ and by the growth condition and the positivity assumption for $d$, we further obtain

[TABLE]

By the Poincaré inequality and positivity of $\alpha,\beta$, we thus get $\|S_{h}(t)\|_{H^{1}(\Omega)}\leq C^{\prime\prime}$ with constant $C^{\prime\prime}$ only depending on the domain, on the bounds in the assumptions, and on the energy of the initial value. This implies that $y(t)$ remains uniformly bounded and hence the solution can be extended uniquely to arbitrarily large $t>0$. ∎

Remark 3.4**.**

Under assumptions (A1)–(A2), we thus obtain global existence of a unique semi-discrete solution $S_{h}$. The constant $C^{\prime\prime}$ in the a-priori bound $\|S_{h}(t)\|_{H^{1}(\Omega)}\leq C^{\prime\prime}$ only depends on the initial energy $E(S_{h,0})$, on the values of $\alpha$ and $\beta$, and on the lower bound $\underline{d}$ for the growth rate of the potential $d(\cdot)$. Let us note that the energy decay estimate also provides a uniform bound for the time derivative of the solution which might be useful for a full convergence analysis for the Galerkin approximation; we refer to [7, 19] for results in this direction.

4. Time discretization

As a second step in the approximation process, let us now discuss the discretization in time. Given some time step $\tau>0$, we define $t^{n}=n\tau$ for $n\geq 0$, and we denote by

[TABLE]

the backward difference quotient at time $t^{n}$. As will become clear from our analysis, non-uniform time steps could be considered as well. For the time discretization of the Galerkin approximation defined in Problem 3.1, we then consider the following implicit time-stepping scheme.

Problem 4.1** (Fully discrete scheme).**

Set $S_{h}^{0}=S_{h,0}$ with initial value as in Problem 3.1. For any $n\geq 1$, find $S_{h}^{n}\in V_{h}$ such that

[TABLE]

for all $v_{h}\in V_{h}$. We set $\frac{d(S_{h}^{n})-d(S_{h}^{n-1})}{S_{h}^{n}-S_{h}^{n-1}}=d^{\prime}(S_{h}^{n})$ wherever $S_{h}^{n}-S_{h}^{n-1}=0$.

Let us note that a similar treatment of the nonlinear term has been proposed in the context of mean-curvature flow in [37]. The particular of the discrete variational problem (4.1) allows us to establish a decay estimate for the energy of the fully discrete approximation with very similar arguments as used for the analysis on the continuous level.

Lemma 4.2**.**

Let $\{S_{h}^{n}\}_{n\geq 0}\subset V_{h}$ denote a solution of Problem 4.1. Then

[TABLE]

Proof.

Using definition (2.3), we can decompose the discrete energy as

[TABLE]

We now utilize the fully discrete variational principle with $v_{h}=S_{h}^{n}-S_{h}^{n-1}$ to replace the first and last term together by $-\tau(c(|\nabla S_{h}^{n-1}|)\bar{\partial}_{\tau}S_{h}^{n},\bar{\partial}_{\tau}S_{h}^{n})$ which directly yields the result. ∎

The second term on the right hand side of the above energy estimate stems from numerical dissipation of the implicit time-stepping scheme. The particular time discretization thus enhances the stability of the numerical solution. As further theoretical backup, we next establish the well-posedness of the fully discrete scheme.

Lemma 4.3**.**

Let $S_{h}^{n-1}\in V_{h}$ be given. Then (4.1) has at least one solution $S_{h}^{n}\in V_{h}$, and the solution is unique for all time steps $0<\tau\leq\tau_{0}$ with some $\tau_{0}$ sufficiently small.

Proof.

Existence of a solution $S_{h}^{n}$ follows from the energy estimate given in the previous section and the Brouwer fixed point theorem, and we know that $\|S^{n}\|_{H^{1}(\Omega)}\leq C^{\prime\prime}$. Now let $S_{h}^{n}$ and $\tilde{S}_{h}^{n}$ be two solutions of (4.1). Then testing with $v_{h}=S_{h}^{n}-\tilde{S}_{h}^{n}$ and using assumption (A1) leads to

[TABLE]

For ease of notation, we utilized the symbol

[TABLE]

here to abbreviate the difference quotient. By the fundamental theorem of calculus and using the bounds of assumption (A2) for the coefficient $d(\cdot)$, we further get

[TABLE]

A combination of the two estimates then directly leads to

[TABLE]

For $\tau>0$ small enough, the right hand side can be absorbed by the terms on the left side and we thus also obtain the uniqueness of the solution. ∎

Before we proceed, let us briefly indicate at this point how the assumption (A1)–(A2) could be possibly relaxed, which can be deduced by inspection of the previous proof.

Remark 4.4**.**

The bounds for $|d^{\prime\prime}(S)|$ are not required for all $S\in\mathbb{R}$ here, but only for values $S$ that are actually attained during the evolution. If a bound $|S_{h}^{n}|\leq C_{h}^{\prime\prime}$ is available, which may be proved by a discrete comparison principles [34], then one could replace the constant $\overline{d}$ in the above proof by $\overline{d}^{\prime}=\sup\{|d^{\prime\prime}(S)|:|S|\leq C_{h}^{\prime\prime}\}$. If the space $V_{h}$ admits an inverse inequality

[TABLE]

then the a-priori estimate $\|S_{h}^{n}\|_{H^{1}(\Omega)}\leq C^{{}^{\prime\prime}}$ implies $|S_{h}^{n}|\leq C_{h}^{{}^{\prime\prime}}$ with $C_{h}^{\prime\prime}=C_{h}C^{{}^{\prime\prime}}$. The constant $C_{h}$ can be estimated explicitly for many discretization schemes, e.g., for finite element methods which discussed in the following section. In one space dimension, the above inverse inequality holds true for any choice of $V_{h}$ due to the continuous embedding of $H^{1}$ into $L^{\infty}$. More general conditions on the potential $d(\cdot)$ have also been considered in [19].

The considerations of this section show that we can expect unique solutions and global well-posedness for the numerical approximations of the model problems (1.1) and (1.6), if the time step $\tau$ is chosen sufficiently small. This is also observed in our computations.

5. Details on the implementation

Before we turn to numerical tests, let us briefly discuss some details of the implementation that will be used in our computations later on. This particularly involves the choice of the approximation spaces $V_{h}$ and the solution of the nonlinear systems (4.1) in every time step.

5.1. A finite element method

Let $T_{h}=\{T\}$ denote a non-overlapping partition of the domain $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ into triangles or tetrahedra. We assume that the partition is conforming and uniformly shape-regular in the usual sense [8, 15]. We then choose

[TABLE]

as the space of continuous piecewise polynomials of degree $k\geq 1$. In our numerical tests we will actually only utilize the lowest order case $k=1$, but the methods can in principle be extended to higher order without difficulty. Similar approximations in space have been considered in [7, 14].

5.2. Nonlinear potential

Without loss of generality, one may assume that the potential $d(\cdot)$ is a polynomial function. Otherwise, $d(\cdot)$ can be replaced by a polynomial and due to the Weierstrass theorem, the error of this approximation can be made arbitrarily small on bounded intervals. In our numerical tests, we will utilize a function of the form

[TABLE]

The first part is the usual double-well potential and the linear term amounts to external forces. Note that the potential $d(\cdot)$ does not strictly satisfy the assumptions (A2), in particular, the upper bounds on the second derivative are violated. As noted in Remark 4.4, this assumption was made for convenience and is not required as long as $S$ stays uniformly bounded. For the above choice of the potential $d(S)$, the difference quotient arising in (4.1) can be expressed as

[TABLE]

The nonlinear problem (4.1) that has to be solved in every time step of the fully discrete scheme thus has a polynomial nonlinearity and is non-degenerate. Standard iterative methods can therefore be utilized for an efficient solution.

5.3. Iterative solver

The following fixed-point iteration will be used in our numerical tests for the solution of the nonlinear problem that arises in every time step of the fully discrete scheme.

Problem 5.1** (Fixed-point iteration).**

Let $S_{h}^{n-1}\in V_{h}$ and $\gamma\geq 0$, $\tau>0$ be given and define $c^{n-1}=c(|\nabla S_{h}^{n-1}|)$. Further set $\tilde{S}_{h}^{n,0}=S_{h}^{n-1}$ and for $k\geq 1$, find $S_{h}^{n,k}\in V_{h}$ such that

[TABLE]

for all $v_{h}\in V_{h}$. If convergence is reached after iteration $k=k^{*}$, set $S_{h}^{n+1}=S_{h}^{n,k^{*}}$.

Choosing $\gamma\approx\beta$ allows to improve the convergence behavior also in case $\beta\gg 0$; see [20] for a similar argument. As an immediate consequence of the particular construction, we obtain

Lemma 5.2**.**

For any $\alpha,\beta,\gamma\geq 0$, the above fixed-point iteration is well-defined. If $\tau>0$ is sufficiently small, then the iteration converges to the unique solution of Problem 4.1.

Proof.

First note that the problem that has to be solved in every step of the fixed-point iteration of Problem 5.1 is linear. Moreover, the left hand side in the above iteration defines a positive definite quadratic form. This implies the existence of a unique solution $S_{h}^{n,k}$ in every step of the iteration. For $\tau>0$ sufficiently small, one can show that the map $\phi:S_{h}^{n,k-1}\to S_{h}^{n,k}$ is contractive, and convergence follows by the Banach fixed-point theorem. ∎

Remark 5.3**.**

For time step $\tau>0$ sufficiently small, the contraction constant in the above iteration becomes small and convergence will typically be observed within a few iterations. To speed up convergence, one may alternatively also utilize Newton-type iterations.

6. Numerical illustration

For our computations, we consider the Allen-Cahn equation (1.1) and the hybrid model (1.6). The model parameters are set to $\lambda=\mu=\nu=0.1$ in all tests, and the double-well potential is defined as $\psi(S)=4S^{2}(1-S)^{2}$. As potential for the external forces, we choose $F(S)=CS$ with constant $C$ to be specified below. The constant $c$ in the Allen-Cahn equation is chosen as $c=\int_{0}^{1}\sqrt{2\psi(S)}dS\approx 0.4714$. For this setting, the interface velocities of Allen-Cahn and the hybrid model given in (1.4) and (1.8) agree up to higher order terms; see [1, 2] for details.

For the space discretization, we utilize the finite element method described in Section 5 with polynomial degree $k=1$. The time integration is performed by the implicit time-stepping scheme (4.1) with constant time step $\tau>0$, and the nonlinear systems arising in every time step are solved with the fixed-point iteration outlined above.

6.1. A quasi one-dimensional problem

As a first test case, we consider a quasi one-dimensional geometric setting. We choose $\Omega=(-3,3)^{2}$ and set

[TABLE]

where $\chi_{D}$ is the characteristic function of the set $D$. The interface $\Gamma$ between the two phases at time $t=0$ is a straight line here and thue the curvature $\kappa_{\Gamma}$ is zero. This can be shown to remain true for all $t>0$ by a symmetry argument. The evolution of the interface is therefore only driven by external forces in this example. For our computations, we choose $F(S)=CS$ with $C=\frac{1}{2}$ which yields a constant driving force $F^{\prime}(S)=\frac{1}{2}$ acting on the whole domain.

As mentioned above, the solution $S(x,y;t)$ can be shown to be independent of $y$ for all $t\geq 0$. We thus obtain $S(x,y;t)=\widehat{S}(x;t)$ with $\widehat{S}$ defined by

[TABLE]

Here $\omega=(-3,3)$ is the one-dimensional cross-section at arbitrary $y$ and $\partial\omega=\{-3,3\}$. The initial condition for the one-dimensional problem is $\widehat{S}(x;0)=\chi_{\{|x|<3/2\}}(x)$. In Figure 6.1, we display a few snapshots of the solution $S^{ac}(x,y;t)$ obtained with the Allen-Cahn equation (1.1) in two dimensions and the corresponding solutions $\widehat{S}^{ac}$ and $\widehat{S}^{hyb}$ for the Allen-Cahn problem and the hybrid model in one space dimension.

The results clearly illustrate the relation $S^{ac}(x,y;t)=\widehat{S}^{ac}(x;t)$ between the solutions of the one- and the two-dimensional problem. Moreover, the propagation of interface $\widehat{\Gamma}=\{x:\widehat{S}(x)=1/2\}$ for the one-dimensional Allen-Cahn and the corresponding hybrid model take place at the same speed, which is in perfect agreement with formulas (1.4) and (1.8). In Figure 6.2, we display the area $A(S(t))=|\{x:\widehat{S}(x;t)>1/2\}|$ surrounded by the interface $\widehat{\Gamma}$ and we depict the decay of the energy functional for the two models illustrating the underlying gradient flow structure.

6.2. Shrinking of a circle by mean curvature

As a second test case, we consider the evolution of a circular interface $\Gamma$. We again choose $\Omega=(-3,3)^{2}$ and define the initial values as

[TABLE]

Here we set $F(S)=0$ which yields $[C]=0$ and thus the motion is only driven by the curvature of the interface. The solution of the sharp interface limit can be shown to be radially symmetric and according to (1.5), the propagation of the interface is governed by

[TABLE]

here $r(t)$ denotes the radial position of the interface $\Gamma(t)$ at time $t$. This allows to compute the area $A(t)$ surrounded by the interface $\Gamma(t)$ analytically, i.e.,

[TABLE]

In Figure 6.3, we display a few snapshots of the interface $\Gamma(t)$ for the numerical solutions of the Allen-Cahn model (1.1) and the hybrid model (1.6).

A brief visual inspection shows that the motion of the interface is rather similar for both models. One can further deduce from the plots that the width of the transition zone in the hybrid model is somewhat larger than for the Allen-Cahn model; compare with our discussion in the introduction. A coarser mesh might therefore be used for the simulation of the interface evolution by the hybrid model. For a better comparison of the two models, we again display in Figure 6.4 the evolution of the area $A(t)$ and of the energy $E(S(t))$ representing one phase of the system.

In accordance with formulas (1.4) and (1.8) for our choice of parameters, the evolution of the interface here occurs at approximately the same speed. Moreover, the energy $E(S(t))$ is monotonically decreasing for both models which is in perfect agreement with the assertions of Lemma 4.2.

6.3. A non-smooth and non-convex geometry

As a last test case, we now consider the evolution of an interface which initially is non-smooth and encloses a non-convex region at time $t=0$. Some snapshots of the evolution are depicted in Figure 6.5.

Like in the previous test, we set the force potental to $F(S)=0$ which implies $[C]=0$, and hence the motion is again driven only by the curvature. Note that in the vicinity of convex corners, the interface moves towards the interior of the enclosed area while at the non-convex corners, the interface propagates in the other direction, which is caused by the change of sign in the curvature. In Figure 6.6, we again depict the area surrounded by the interface and the uniform decay of the energy for the numerical solutions obtained with the two models.

Like in the previous examples, the velocity of the interface motion in the Allen-Cahn and the hybrid model are very similar; compare with the formulas (1.4) and (1.8). Due to the gradient flow structure of the fully discrete evolution, a strict decrease in the energy $E(S(t))$ can again be observed. The results are very similar, and the slight discrepancies can be explained by approximation errors.

7. Discussion

In this paper, we investigated the systematic numerical approximation of a general class of Allen-Cahn type equations. The common feature of these models was a gradient flow structure with respect to an associated energy. We proposed and analyzed semi-discrete and fully discrete numerical schemes which strictly preserve the underlying gradient flow structure and which therefore automatically yield uniformly energy stable discrete approximations. Well-posedness of the numerical schemes and energy decay was established theoretically and illustrated by numerical tests. Our computational results also provide further numerical evidence for asymptotic expansions of interface motion obtained [1, 2].

In our analysis, we used some rather strong structural assumptions on the coefficients arising in the equations, which should allow to conduct a full a-priori convergence analysis of the semi- and fully discrete scheme by following the arguments of [7, 19]. As indicated in remarks, a detailed inspection of the proofs allows to further relax these assumptions. Our arguments and methods of proof are rather general and seem to be applicable also to phase-field models for the dynamics of solid-solid phase transition which are governed by coupled Allen-Cahn and elasticity equations. These topics and further extensions are left for future research.

Acknowledgments

This work was supported by the German Research Foundation (DFG) via grants IRTG 1529, TRR 154, and Eg-331/1-1 and by the German Excellence Initiative via grant GSC 233.