An alternative to the Allen-Cahn phase field model for interfaces in solids - numerical efficiency

TL;DR

This paper introduces an alternative phase field model for solid interfaces that is more numerically efficient than the traditional Allen-Cahn model, especially for low interface energy scenarios.

Contribution

It formulates a new phase field equation inspired by sharp interface problems and compares its properties and efficiency to the Allen-Cahn model.

Findings

The alternative model is more effective for low interface energy simulations.

Numerical simulations show improved efficiency with the new model.

The new model maintains key interface evolution properties.

Abstract

The derivation of the Allen-Cahn and Cahn-Hilliard equations is based on the Clausius-Duhem inequality. This is not a derivation in the strict sense of the word, since other phase field equations can be fomulated satisfying this inequality. Motivated by the form of sharp interface problems, we formulate such an alternative equation and compare the properties of the models for the evolution of phase interfaces in solids, which consist of the elasticity equations and the Allen-Cahn equation or the alternative equation. We find that numerical simulations of phase interfaces with small interface energy based on the alternative model are more effective then simulations based on the Allen-Cahn model.

Equations

$$\gamma (t)=\{x\in \Omega \mid S(t,x)=0\},{\gamma }^{\prime }(t)=\{x\in \Omega \mid S(t,x)=1\},$$

$$\varepsilon(\nabla_{x}u)=\frac{1}{2}\big{(}\nabla_{x}u+(\nabla_{x}u)^{T}\big{)}\in{\cal S}^{3}$$

$$-{\mathrm{div}}_{x}T$$

$$T$$

$$u(t,x)$$

$${\sf W}\big{(}\varepsilon(\nabla_{x}u),S\big{)}=\frac{1}{2}\Big{(}D\big{(}\varepsilon(\nabla_{x}u)-\overline{\varepsilon}S\big{)}\Big{)}:\big{(}\varepsilon(\nabla_{x}u)-\overline{\varepsilon}S\big{)},$$

$$\psi^{*}\big{(}\varepsilon(\nabla_{x}u),S,\nabla_{x}S\big{)}={\sf W}\big{(}\varepsilon(\nabla_{x}u),S\big{)}+\hat{\psi}(S)+\frac{1}{2}|\nabla_{x}S|^{2}.$$

$$\hat{\psi }(0)=\hat{\psi }(1)=0,\hat{\psi }(r)>0,\text{mbox}for0<r<1.$$

$$\frac{\partial }{\partial t}{\psi }^{\ast }+{\mathrm{div}}_{x}q\le \mathsf{b}\cdot u_t$$

$$q=-T{u}_t-S_t{\nabla }_{x}S.$$

$$0\geq\frac{\partial}{\partial t}\psi^{*}+{\rm{div}}_{x}q-{\sf b}\cdot u_{t}=\partial_{(\nabla_{x}u)}{\sf W}:\nabla_{x}u_{t}+\big{(}\partial_{S}{\sf W}+\hat{\psi}^{\prime}(S)\big{)}S_{t}+\nabla_{x}S\cdot\nabla_{x}S_{t}\\ -{\rm{div}}_{x}(Tu_{t})-{\rm{div}}_{x}(S_{t}\nabla_{x}S)-{\sf b}\cdot u_{t}=(\partial_{S}{\sf W}+\hat{\psi}^{\prime}(S)-\Delta_{x}S)S_{t}\,.$$

$$\partial_{t}S=-f\Big{(}\partial_{S}{\sf W}\big{(}\varepsilon(\nabla_{x}u),S\big{)}+\hat{\psi}^{\prime}(S)-\Delta_{x}S\Big{)},$$

$$\mathcal{F}={\partial }_S\mathsf{W}+{\hat{\psi }}^{\prime }(S)-{\Delta }_{x}S$$

$$S_{t}(t,x)=-f\big{(}{\cal F}(t,x)\big{)}.$$

$$s(t,x(t))=\frac{\mathrm{d}x(t)}{\mathrm{d}t}\cdot \frac{{\nabla }_{x}S(t,x(t))}{|{\nabla }_{x}S(t,x(t))|}.$$

$$0=\frac{d}{dt}S\big{(}t,x(t)\big{)}=S_{t}\big{(}t,x(t)\big{)}+\frac{dx(t)}{dt}\cdot\nabla_{x}S\big{(}t,x(t)\big{)}\\ =S_{t}\big{(}t,x(t)\big{)}+s\big{(}t,x(t)\big{)}|\nabla_{x}S\big{(}t,x(t)\big{)}|,\qquad x(t_{0})\in\Gamma_{c}(t_{0}),$$

$$S_t+s|{\nabla }_{x}S|=0$$

$$S_{t}(t,x)=-f\big{(}{\cal F}(t,x)\big{)}|\nabla_{x}S(t,x)|,$$

$$S_{t}=-f\big{(}\partial_{S}{\sf W}+\hat{\psi}^{\prime}(S)-\Delta_{x}S\big{)}|\nabla_{x}S|,$$

$$-{\mathrm{div}}_{x}T$$

$$T$$

$${\partial }_{t}S$$

$$u(t,x)$$

$${\partial }_{n_{\partial \Omega }}S(t,x)$$

$$S(0,x)$$

$$-{\mathrm{div}}_{x}T$$

$$T$$

$${\partial }_{t}S$$

$$Q=[t_1,t_2]\times \Omega ,$$

$$\Gamma$$

$$\gamma$$

$${\gamma }^{\prime }$$

$$n:\Gamma \to {\mathbb{R}}^3$$

$$w^{(\pm )}(t,x)$$

$$[w](t,x)$$

$$⟨w⟩(t,x)$$

Full text

Dedicated to the memory of Krzysztof Wilmański

1 Introduction

The phase field approach is used to model the evolution of phase interfaces in many different materials and accordingly the resulting models differ widely. However, in spite of all the differences the evolution equations for the order parameter $S$ in the models is almost always formulated by the standard approach to set the time derivative of the order parameter equal to a suitable function of the functional derivative of the Ginzburg-Landau free energy with respect to $S$, which leads to an Allen-Cahn type equation, or equal to the divergence of a suitable function of the gradient of this functional derivative, which leads to a Cahn-Hilliard type equation. Often this function is chosen to be linear. For a thorough discussion of this approach to formulate material models with the Allen-Cahn und Cahn-Hilliard equation we refer to [8].

The leading idea behind these approaches to formulate the evolution equation is that in both cases for the resulting model the Clausius-Duhem inequality is guaranteed to hold. Yet, there are other possibilities to choose the evolution equation such that this inequality holds. Therefore the question arises whether the standard approach is always the best or whether there are situations where other choices of the evolution equation for the order parameter lead to better results.

Of course, this question can only be discussed at a concrete example of an alternative phase field equation in a concrete mathematical material model. We consider here the prototypic model for the evolution of phase interfaces in solids, neglecting temperature effects, which consists of the elasticity equations coupled to an evolution equation for $S$. For this evolution equation one usually inserts the Allen-Cahn equation. We formulate here an alternative phase field equation and compare the mathematical properties of the two different models, which are obtained when we use the Allen-Cahn equation or the alternative equation as the evolution equation. Our main result is that simulations of phase interfaces in solids, which have small or vanishing interface energy density, are numerically more effective when the alternative equation is used instead of the Allen-Cahn equation.

We stress that the alternative phase field equation can replace the Allen-Cahn equation in other models. The properties of the resulting models have as yet to be investigated.

This paper is based on our investigations of phase field equations in the articles [1] – [7]. It summarizes in particular the results obtained in [6] and [7], but adds also some new considerations.

2 The Clausius-Duhem inequality and the Allen-Cahn

equation

To formulate the alternative phase field equation, we must know the form of the Ginzburg-Landau free energy, which appears in the Clausius-Duhem inequality. Therefore we first introduce the physical situation and the elasticity equations, from which the form of the Ginzburg-Landau free energy results.

Let $\Omega\subseteq{\mathbb{R}}^{3}$ be an open bounded set, which represents a solid body. We assume that the atoms of the material can be arranged in crystal lattices of two different types. The crystal type present at a material point $x\in\Omega$ at time $t$ is indicated by the order parameter. The value $S(t,x)=0$ means that type one is present, $S(t,x)=1$ indicates that type two is present. The sets of points

[TABLE]

where crystal type one or crystal type two is present, respectively, are called phase $1$ or phase $2$ of the material at time $t$, respectively. Let $u(t,x)\in{\mathbb{R}}^{3}$ denote the displacement of the material point $x$ at time $t$ and let

[TABLE]

be the linear strain tensor, where ${\cal S}^{3}$ denotes the set of symmetric $3\times 3$–matrices. We assume that only small displacements occur and we consider a quasistatic model. This means that for every given time $t$ the displacement $x\mapsto u(t,x)$ and the Cauchy stress tensor $x\mapsto T(t,x)\in{\cal S}^{3}$ must solve the boundary value problem of linear elasticity posed in the domain $\Omega$, which is given by

[TABLE]

where $\overline{\varepsilon}\in{\cal S}^{3}$ is the given transformation strain, where $D:{\cal S}\rightarrow{\cal S}$ is the elasticity tensor, a linear, symmetric, positive definite mapping, and where ${\sf b}(t,x),{\sf U}(t,x)\in{\mathbb{R}}^{3}$ denote the given volume force and boundary displacement. By (2.2), the material is stree free in phase one if $\varepsilon(\nabla_{x}u)$ is equal to zero, and in phase two if $\varepsilon(\nabla_{x}u)$ is equal to the transformation strain $\overline{\varepsilon}$.

To close the system of model equations, we need an evolution equation for $S$. To formulate it, note that according to (2.2), the stored elastic energy is

[TABLE]

which leads to the Ginzburg-Landau free energy

[TABLE]

where $\hat{\psi}:{\mathbb{R}}\rightarrow{\mathbb{R}}$ is a double well potential satisfying

[TABLE]

The second law of thermodynamics requires that there is a flow of the free energy q\big{(}u,u_{t},\varepsilon(\nabla_{x}u),S,S_{t},\nabla_{x}S\big{)}, such that the Clausius-Duhem inequality

[TABLE]

holds for all solutions $(u,T,S)$ of the model equations. We use the flow

[TABLE]

If we insert (2.5) and (2.8) into (2.7) and note (2.1) and the equation $\partial_{(\nabla_{x}u)}{\sf W}=T$, which follows from (2.4), (2.2), then we obtain by a short computation that

[TABLE]

The Clausius-Duhem inequality (2.7) is therefore satisfied, if the evolution equation for $S$ guarantees that the right hand side of (2.9) is non-positive. The simplest possibility to obtain this is to set

[TABLE]

with a function $f:{\mathbb{R}}\rightarrow{\mathbb{R}}$ satisfying $r\cdot f(r)\geq 0$. If for $f$ the linear function $f(r)=cr$ is chosen with a positive constant $c$, then the Allen-Cahn equation results.

(2.1), (2.2), (2.10) form a closed system of partial differential equations. The standard phase field model for the evolution of phase interfaces consists of this system, combined with the boundary condition (2.3) and an initial condition for $S$.

3 Formulation of an alternative phase field equation

By the inequality (2.9), the expression

[TABLE]

and the time derivative $S_{t}$ must have opposite signs, which means that the value of $S_{t}$ at $(t,x)$ cannot be independent of the value ${\cal F}(t,x)$. Instead, there must be a functional relation between both values. Of course, this does not mean that $S_{t}$ must depend on ${\cal F}$ alone as in the ansatz (2.10), it can depend on additional variables as well. The question arises, on which other variables $S_{t}$ should depend.

To discuss this question we start from the usual physical interpretation of the observation, that there must be a functional relation between $S_{t}$ and ${\cal F}$. The interpretation is that ${\cal F}$ is a configurational force, which drives the time evolution of the order parameter $S$. This interpretation is used as an additional justification for the equation (2.10), which we write in the short form

[TABLE]

What one wants to have is that the variation of the order parameter $S$ is confined to a narrow diffuse interface, which moves with a propagation speed, which is a linear or nonlinear function of the configurational force ${\cal F}$. In fact, standard sharp interface models contain an equation, which prescribes the propagation speed of the interface as a function of the configurational driving force. This equation is called kinetic relation. We extend the meaning of this notation also to phase field models.

In a standard sharp interface model the kinetic relation can therefore be explicitly read off from the model equations. It would be of interest to have a phase field model, where the kinetic relation can also be read off directly from the form of the model equations. For the phase field equation (3.2) this is not possible. Instead, the kinetic relation is a hidden property of this equation, which must be determined by a very technical asymptotic analysis of this equation.

Our goal is therefore to formulate a phase field equation, for which the Clausisus-Duhem inequality (2.7) is satisfied, and which allows to read off the kinetic relation directly from the form of the equation. To formulate such an equation, assume that $S$ is an order parameter, whose transition from [math] to $1$ defines a diffuse phase interface moving in time. We say that the speed of the diffuse interface at $(t,x_{0})$ is equal to the normal speed $s(t,x_{0})$ of the level set $\Gamma_{c}(t)=\{x\in\Omega\mid S(t,x)=c\}$, which contains $x_{0}$. The normal speed of $\Gamma_{c}(t)$ at $x\in\Gamma_{c}(t)$ can be defined as follows: If $\tilde{t}\mapsto x(\tilde{t})\in{\mathbb{R}}^{3}$ is a function defined for all $\tilde{t}$ from a neighborhood of $t$ and if $x(\tilde{t})\in\Gamma_{c}(\tilde{t})$ holds for all $\tilde{t}$, then the normal speed $s(t,x)$ of $\Gamma_{c}(t)$ at $x=x(t)\in\Gamma_{c}(t)$ is the component of the velocoity $x^{\prime}(t)$ in the direction of the unit normal vector $n(t,x)$ to $\Gamma_{c}(t)$ at $x$. Since $n(t,x)=\frac{\nabla_{x}S(t,x)}{|\nabla_{x}S(t,x)|}$, we obtain

[TABLE]

The function $t\mapsto x(t)$ satisfies $x(t)\in\Gamma_{c}(t)$ if and only if t\mapsto S\big{(}t,x(t)\big{)}=c holds, and this last equation holds if and only if for a fixed time $t_{0}$ the function $x(t)$ satisfies the initial value problem

[TABLE]

with $s$ defined by (3.3). From this we conclude that if $t_{1}<t_{2}$ are given times and if $s:[t_{1},t_{2}]\times\Omega\rightarrow{\mathbb{R}}$ is a given function, then $S$ satisfies the partial differential equation

[TABLE]

in the domain $[t_{1},t_{2}]\times\Omega$, if and only if every level set $\Gamma_{c}(t)$ moves with normal speed $s(t,x)$ at $x\in\Gamma_{c}(t)$.

This suggests to combine the equations (2.1) – (2.3) with the evolution equation

[TABLE]

with the driving force ${\cal F}$ defined by (3.1) and with a given linear or nonlinear function $f:{\mathbb{R}}\rightarrow{\mathbb{R}}$. If we compare (3.4) and (3.5), then we see that the propagation speed of the diffuse interface defined by (3.5) is equal to s=f\big{(}{\cal F}(t,x)\big{)}, whence the kinetic relation is given by $f$ and can be read off directly from the evolution equation (3.5). From (2.9) we immediately see that every solution $(u,T,S)$ of the equations (2.1), (2.2), (3.5) satisfies the Clausius-Duhem inequality (2.7) if $f$ satisfies $r\cdot f(r)\geq 0$ for all $r\in{\mathbb{R}}$. The evolution equation (3.5) has therefore the desired properties.

(3.5) has the form of a Hamilton-Jacobi equation. However, if one inserts the definition (3.1) of ${\cal F}$ into (3.5), one obtains the phase field equation

[TABLE]

which is degenerate parabolic. (3.6) has therefore mixed hyperbolic–parabolic properties. This is why we call (3.6) hybrid phase field equation.

4 The Allen-Cahn and the hybrid models

We have now two different phase field models for the evolution of phase interfaces in solids: If we combine the equations (2.1), (2.2) with the phase field equation (2.10) of Allen-Cahn type we obtain the system

[TABLE]

which must be solved in the domain $[0,\infty)\times\Omega$. As boundary and initial conditions we choose, for example,

[TABLE]

To obtain (4.3) from (2.10) we specialized the function $f$ in (2.10) to be $f(r)=cr$ with a positive constant $c$ and we introduced two scaling parameters $\mu>0$ and $\lambda>0$, whose meaning will become clear later. To have a short name, we call the system (4.1) – (4.3) the Allen-Cahn phase field model.

The second model is obtained by combination of (2.1), (2.2) with the hybrid phase field equation (3.6). If we specialize the function $f$ in (3.6) to be $f(r)=cr$ with a constant $c>0$ and introduce a scaling parameter $\nu>0$, the resulting system is

[TABLE]

These equations must be solved in the domain $[0,\infty)\times\Omega$. For the boundary and initial conditions we can again take (4.4) – (4.6). We call the system (4.7) – (4.9) the hybrid phase field model.

Several questions arise immediately. (4.9) is a quasilinear, degenerate parabolic equation. Little is known about equations of the form (4.9). The first question therefore concerns existence and uniqueness of solutions to the system (4.7) – (4.9). Moreover, if solutions $(u,T,S)$ exist, does the function $S$ have the properties required from an order parameter? If both questions can be answered positively, what is then the difference between the Allen-Cahn model and the hybrid model? We have studied these questions in recent years. To the first two questions only partial answers can be given, whereas the answer to the third question is quite well known.

In [3] it is proved that weak solutions of the hybrid model (4.7) – (4.9), (4.4) – (4.6) exist in the case of one space dimension. The proof is based on the observation that the one-dimensional version of the evolution equation (4.9) has some monotonicity properties. In higher space dimensions no rigorous existence proof is available. We must therefore rely on extensive numerical tests and on formal asymptotic analysis. The numerical test computations seem to indicate quite clearly, that solutions $(u,T,S)$ exist and that the function $S$ in these solutions has the properties required from an order parameter. In fact, the test computations converge in higher space dimensions better then in one space dimension. A part of the test computations is documented in [6].

The last question on the difference of the models is answered in the remainder of this paper. Of course, to answer the question we need to have more information on the properties of the models. This information is collected in Sections 5 and 6. The information is obtained by asymptotic analysis of the models, more precisely by construction of approximate solutions to the Allen-Cahn and the hybrid models. The answer to the comparison question is finally given in Section 7.

5 Model error and asymptotics

To compare the models we need to define what we understand under the model error. In this section we first give this definition and subsequently state in Theorems 5.2 and 5.3 some results on approximate solutions, which have been obtained in [6] and [7].

To define the model error we must first specify the type of material interfaces, which we want to model. Of great current interest are phase interfaces in functional materials. Very often such interfaces are thin and consist only of a few atomic layers. A large number of phase field models to simulate the time evolution of such interfaces have been devised and more are developed. It is therefore of interest to study how well the Allen-Cahn and the hybrid models are adapted to the simulation of thin interfaces in solids. More precisely, it is of interest to study how large the difference between the propagation speed of a thin phase interface in the real material and of the interface in the respective phase field model is. This difference is the model error.

To give a precise definition of the model error, we must approximately know the propagation speed of the real phase interface. For very thin interfaces mathematical models with sharp interface are appropriate. We therefore base the following considerations on the hypothesis that the propagation speed of the interface in the sharp interface model is a good approximation to the propagation speed of the interface in the real material. The model error of a phase field model is then the difference of the propagation speed of the sharp interface and the propagation speed of the diffuse interface in the phase field model.

To formulate the sharp interface model to be used we must introduce some notations. The asymptotic solution is constructed in the bounded domain

[TABLE]

where $0\leq t_{1}<t_{2}<\infty$ are given times. $\Gamma(t)$ denotes the sharp interface at time $t$. We assume that the phase sets $\gamma(t)$, $\gamma^{\prime}(t)$ introduced in Section 2 are open, disjoint subsets of $\Omega$, whose common boundary is $\Gamma(t)$, such that $\Omega=\gamma(t)\cup\gamma^{\prime}(t)\cup\Gamma(t)$. We set

[TABLE]

Let

[TABLE]

be the continuous vector field, for which $n(t,x)$ is the unit normal vector to $\Gamma(t)$ at $x\in\Gamma(t)$, which points into the domain $\gamma^{\prime}(t)$. For a function $w$ defined in a neighborhood of $\Gamma$ and $(t,x)\in\Gamma$ we set

[TABLE]

Now we can formulate the sharp interface model. Let $\hat{S}:Q\rightarrow\{0,1\}$ be a piecewise constant function, which only takes the values [math] and $1$ with a jump across $\Gamma$, such that

[TABLE]

The sharp interface model consists of a transmission problem for the elasticity equations and of a kinetic relation. The transmission problem is given by

[TABLE]

To determine the kinetic relation we proceed as in Section 2. We use the Clausius-Duhem inequality

[TABLE]

with the free enery and the flux

[TABLE]

where $c_{1}\geq 0$ is an arbitrarily chosen constant. The last term on the right hand side of (5.7) is the interface energy, hence $\lambda^{1/2}c_{1}$ is the interface energy density. It is well known that if $(\hat{u},\hat{T})$ is a solution of the transmission problem (5.1) – (5.5) and if the interface $\Gamma(t)$ in this problem moves with the normal speed $s_{\rm sharp}(t,x)$ at $x\in\Gamma(t)$, then the Clausius-Duhem inequality (5.6) holds if and only if the inequality

[TABLE]

is satisfied at every point $x\in\Gamma(t)$. Here $\kappa_{\Gamma}(t,x)$ denotes twice the mean curvature of the surface $\Gamma(t)$ at $x\in\Gamma(t)$.

A proof of this well known result is given in [1], however only for the case where $c_{1}=0$ in (5.7). The proof can be readily generalized to the case $c_{1}>0$.

A simple linear kinetic relation, for which (5.8) obviously holds, is

[TABLE]

with a positive constant $\hat{c}$. The sharp interface problem thus consists of the transmission problem (5.1) – (5.5) combined with the kinetic relation (5.9).

We can now define the model error. To this end note that solutions of the Allen-Cahn model depend on the parameters $\mu$ and $\lambda$, whereas solutions of the hybrid model depend on the parameter $\nu$. Therefore we record these parameters in the notation. For a solution $(u^{(\mu\lambda)}_{\rm AC},T^{(\mu\lambda)}_{\rm AC},S^{(\mu\lambda)}_{\rm AC})$ of the Allen-Cahn model and for a solution $(u^{(\nu)}_{\rm hyb},T^{(\nu)}_{\rm hyb},S^{(\nu)}_{\rm hyb})$ of the hybrid model consider the level sets

[TABLE]

Let $s_{\rm AC}^{(\mu\lambda)}(t,x)$ denote the normal speed of $\Gamma_{\rm AC}^{(\mu\lambda)}(t)$ at $x\in\Gamma_{\rm AC}^{(\mu\lambda)}(t)$, and let $s_{\rm hyb}^{(\nu)}(t,x)$ denote the normal speed of $\Gamma_{\rm hyb}^{(\nu)}(t)$ at $x\in\Gamma_{\rm hyb}^{(\nu)}(t)$. These normal speeds are approximately equal to the propagation speeds of the diffuse phase interfaces defined by the solutions of the Allen-Cahn and hybrid models.

Let $t\in[t_{1},t_{2}]$ be a given, fixed number. As initial conditions for the sharp interface problem we can choose

[TABLE]

Definition 5.1

We call the functions ${\cal E}^{(\mu\lambda)}(t):\Gamma(t)\rightarrow{\mathbb{R}}$ and ${\cal E}^{(\nu)}(t):\Gamma(t)\rightarrow{\mathbb{R}}$, respectively, which are defined by

[TABLE]

the error of the Allen-Cahn model or the error of the hybrid model at time $t$, respectively.

We next state some results for the Allen-Cahn and hybrid models obtained by asymptotic analysis.

By $B_{\rm AC}^{(\mu\lambda)}>0$ and $B_{\rm hyb}^{(\nu)}>0$ we denote the widths of the diffuse interfaces defined by the order parameter in solutions of the Allen-Cahn model and by the order parameter in solutions of the hybrid model. Here we do not define the interface width precisely. If $S$ is an order parameter, one could define the interface width to be the maximal distance between the level surfaces $\{x\in\Omega\mid S(t,x)=0.1\}$ and $\{x\in\Omega\mid S(t,x)=0.9\}$, for example. We are interested in the limits $\mu\rightarrow 0$, $\lambda\rightarrow 0$, $\nu\rightarrow 0$ and assume therefore that $\mu\in(0,\mu_{0}]$, $\lambda\in(0,\lambda_{0}]$, $\nu\in(0,\nu_{0}]$, with suitably chosen fixed constants $\mu_{0},\lambda_{0},\nu_{0}>0$.

Theorem 5.2

Let $(u^{(\mu\lambda)}_{\rm AC},T^{(\mu\lambda)}_{\rm AC},S^{(\mu\lambda)}_{\rm AC})$ be a solution of the Allen-Cahn model (4.1) – (4.5), let $t\in[t_{1},t_{2}]$ be a given time, and let $(\hat{u}(t),\hat{T}(t))$ be the solution of the transmission problem (5.1) – (5.5) with the interface given by $\Gamma(t)=\Gamma_{\rm AC}^{(\mu\lambda)}(t)$. Then

[TABLE]

where s_{0}=s_{0}\big{(}\Gamma_{\rm AC}^{(\mu\lambda)}(t)\big{)}, s_{10}=s_{10}\big{(}\Gamma_{\rm AC}^{(\mu\lambda)}(t)\big{)} and s_{11}\big{(}\Gamma_{\rm AC}^{(\mu\lambda)}(t)\big{)} are nonlocal functions of $\Gamma_{\rm AC}^{(\mu\lambda)}(t)$. In particular, we have

[TABLE]

with the constant

[TABLE]

For the remainder term $R_{\rm AC}(\mu,\lambda,t,x)$ there is a function $\mu\rightarrow{\cal C}_{\cal E}(\mu)$ with the property that $\lim_{\mu\rightarrow 0}{\cal C}_{\cal E}(\mu)=0$, such that for all $0<\mu\leq\mu_{0}$, $0<\lambda\leq\lambda_{0}$ and all $(t,x)\in\Gamma_{\rm AC}^{(\mu\lambda)}$ the inequality

[TABLE]

holds. Moreover, there is a constant $C_{1}>0$ such that for all $0<\mu\leq\mu_{0}$, $0<\lambda\leq\lambda_{0}$

[TABLE]

These results are contained in [7]. We stress here the fact, that the results are obtained by formal asymptotic analysis. No rigorous mathematical proof of these statements is given in [7]. The asymptotic analysis with respect to $\mu\rightarrow 0$ uses mathematical methods, which are standard in the analysis of phase field models. This is different for the estimate (5.15), which says that the remainder term $R_{\rm AC}$ tends to zero for $\mu\rightarrow 0$, uniformly with respect to $\lambda$. This uniformity estimate is obtained by a second asymptotic analysis with respect to $\lambda\rightarrow 0$. The formal derivation of this estimate is a novelty introduced in [7].

Theorem 5.3

Let $(u^{(\nu)}_{\rm hyb},T^{(\nu)}_{\rm hyb},S^{(\nu)}_{\rm hyb})$ be a solution of the hybrid model (4.7) – (4.9), (4.4), (4.5), let $t\in[t_{1},t_{2}]$ be a given time, and let $(\hat{u}(t),\hat{T}(t))$ be the solution of the transmission problem (5.1) – (5.5) with the interface given by $\Gamma(t)=\Gamma_{\rm hyb}^{(\nu)}(t)$. Then

[TABLE]

where $c>0$ is the constant from (4.9). For the remainder term $R_{\rm hyb}(\nu,t,x)$ there is a constant $C_{2}$ such that for all $0<\nu\leq\nu_{0}$ and all $(t,x)\in\Gamma_{\rm hyb}^{(\nu)}$ the inequality

[TABLE]

holds. Moreover, there is a constant $C_{3}>0$ such that for all $0<\nu\leq\nu_{0}$

[TABLE]

These results are obtained in [6], again by formal asymptotic analysis.

6 Characteristic equations

From the results on the asymptotic behavior of the models stated in Theorems 5.2 and 5.3 we derive in this section for both models some relations between parameters of the models. We call these relations the characteristic relations of the models. The comparison of the models in Section 7 is based on these relations.

We first consider the Allen-Cahn model. For $c_{1}$ in the free energy (5.7) we choose the value given by (5.14), With this value we adapt the interface energy density $\lambda^{1/2}c_{1}$ to the value in the real material by varying $\lambda$. In (4.3) we choose $c=\hat{c}c_{1}$. By (5.9) and (5.13) we then have

[TABLE]

hence (5.10) and (5.12) together imply

[TABLE]

This equation and (5.15) together yield

[TABLE]

with a constant $C$, which can be chosen independently of $\lambda$. By this inequality, $\mu^{1/2}$ controls the model error. Therefore we write $F=\mu^{1/2}$ and call $F$ the error parameter. Moreover, since $\lambda^{1/2}c_{1}$ is the interface energy density, we call $E=\lambda^{1/2}$ the interface energy parameter. Also, since by (5.16) the interface width is bounded by a constant, which is proportional to $(\mu\lambda)^{1/2}$, we call ${\cal W}=(\mu\lambda)^{1/2}$ the interface width parameter. These three parameters and the propagation speed $s_{\rm AC}=s_{\rm AC}^{(\mu\lambda)}$ are connected by the fundamental relations

[TABLE]

where we use the notation ${\cal E}[E,F]={\cal E}^{(\mu\lambda)}$. The first equation is an immediate consequence of the definition of the parameters, the second is obtained by insertion of (5.9) into (5.10), and the last inequality is just a restatement of (6.2).

Now assume that we want to use a phase field model to numerically simulate the propagation of a phase interface. In such a simulation the numerical effort is proportional to $h^{-p}$, where $h$ denotes the grid spacing and where the power $p>1$ depends on whether we want to simulate a problem in $2$–d or in $3$–d and it depends on the numerical scheme we use. In order for the simulation to be precise, we must guarantee that the model error and the numerical error are small. To make the numerical error small, we must choose the grid spacing $h$ small enough to resolve the transition of the order parameter across the interface, which means that we must choose $h<{\cal W}$, hence we have $h^{-p}>{\cal W}^{-p}$. Therefore we see that the numerical effort of a simulation based on a phase field model is measured by the number ${\cal W}^{-p}$. We call the number

[TABLE]

the parameter of numerical effort. For a simulation based on the Allen-Cahn model we see from (6.3) that the numerical effort is

[TABLE]

We call the relations (6.3) – (6.6) characteristic relations for the Allen-Cahn model.

Next we derive the characteristic relations for the hybrid model. In the free energy (5.7) we choose $c_{1}=0$, and in (4.9) we set $c=\hat{c}$. By (5.9) and (5.17) we then have

[TABLE]

We insert this equation into (5.11) and obtain for the model error

[TABLE]

From this equation and from (5.18) we infer that

[TABLE]

By this equation, $\nu^{1/2}$ controls the model error. In the case of the hybrid model we therefore choose $F=\nu^{1/2}$ as the error parameter. By (5.19), the interface width is bounded by a constant, which is proportional to $\nu^{1/2}$, whence the interface width parameter is ${\cal W}=\nu^{1/2}$. For the hybrid model we therefore have the characteristic relations

[TABLE]

where we used the notations $s_{\rm hyb}=s_{\rm hyb}^{(\nu)}$ and ${\cal E}[F]={\cal E}^{(\nu)}$. The first of these relations follows from the definitions of $F$ and ${\cal W}$, the second one is obtained by combination of (6.7) and (6.8), noting (5.9), the third one is just a restatement of (6.9), and the last one follows from the definition $e_{\rm num}={\cal W}^{-p}$ of the parameter of numerical effort and from (6.10).

7 Comparison of the models, numerical efficiency

From (6.4) we see that the Allen-Cahn model can describe the evolution of a phase interface with propagation speed $\hat{c}\,(-\overline{\varepsilon}:\langle\hat{T}\rangle+c_{1}\kappa_{\Gamma}E)$, which by (5.9) is the propagation speed of an interface with interface energy density $c_{1}\lambda^{1/2}=c_{1}E$. The interface energy density is always positive, since we cannot set $\lambda=0$ in the Allen-Cahn equation (4.3). Varying of the parameter $E$ to adjust the interface energy density does not change the model error; this error can be adjusted to a desired value by choosing the parameter $F=\mu^{1/2}$ suitably. Varying of $F$ does not change the interface energy density. From (6.6) we see that if the interface energy density parameter $E$ is fixed, then the effort of a numerical simulation grows with $F^{-p}$, where the power $p>1$ depends on the numerical method employed and on the space dimension of the problem, which we want to simulate.

From (6.11) we see that the hybrid model, on the other hand, can describe the evolution of a phase interface with propagation speed $-\hat{c}\>\overline{\varepsilon}:\langle\hat{T}\rangle$, which by (5.9) is the propagation speed of an interface with interface energy density $c_{1}\lambda^{1/2}=0$. The model error can be adjusted to a desired value by choosing the parameter $F=\nu^{1/2}$ suitably. By (6.13), also for this model the effort of a numerical simulation grows with $F^{-p}$, where the power $p>1$ depends on the numerical method employed and on the space dimension of the problem, which we want to simulate.

These observations suggest the following rule:

Simulations of phase interfaces with positive interface energy density should be based on the Allen-Cahn model, simulations of interfaces with zero or small interface energy density should be based on the hybrid model.

One can object to this rule by arguing that the Allen-Cahn model can also be used to simulate interfaces with zero interface energy density by choosing the interface energy density parameter positive, but very small. However, because of the presence of the factor $E^{-p}$ in the formula (6.6) the numerical effort will become very large.

To be more specific, we consider an interface with vanishing interface energy density, hence $c_{1}\lambda^{1/2}=0$, which by (5.9) means that the propagation speed of the sharp interface is

[TABLE]

For the Allen-Cahn model it follows from this equation and from (6.4) that in this case the total model error, which we denote by ${\cal E}_{\rm total}$, is

[TABLE]

This means that the term $\hat{c}c_{1}\kappa_{\Gamma}E$ is now part of the total model error.

If we prescribe the maximal value ${\cal E}_{\max}$ of the total model error $|{\cal E}_{\rm total}|$, we must therefore choose the parameters $E$ and $F$ such that

[TABLE]

where the second condition is imposed by the requirement to make the numerical effort $e_{\rm num}=(EF)^{-p}$ as small as possible. To discuss this optimization problem, we assume first that the term $s_{10}$ in the asymptotic expansion (5.12) of the propagation speed $s^{(\mu\lambda)}_{\rm AC}$ is not identically equal to zero. In this case we conclude from (6.1) and (5.15) that for sufficiently small $\lambda^{1/2}=E$ and for sufficiently small $\mu^{1/2}=F$ the error ${\cal E}[E,F]={\cal E}^{(\mu\lambda)}$ satisfies

[TABLE]

This inequality and (7.1) imply that the solution $(E,F)$ of the optimization problem (7.1), (7.2) satisfies

[TABLE]

From this result we obtain

Corollary 7.1

Let ${\cal E}_{\rm max}$ denote the total model error of the Allen-Cahn model in the simulation of an interface without interface energy. If the term $s_{10}$ in the asymptotic expansion (5.12) of the propagation speed $s^{(\mu\lambda)}_{\rm AC}$ is not identically equal to zero, then the interface width $B_{\rm AC}$ satisfies

[TABLE]

In a numerical simulation of an interface without interface energy based on the Allen-Cahn model the parameter of numerical effort satisfies

[TABLE]

with a power $p>1$ depending on the space dimension and the numerical method used.

For the hybrid model we have by (6.11) and (6.12) that ${\cal E}_{\max}=\max_{\Gamma}|{\cal E}[F]|\leq CF$. From (6.13) and from (7.4) we thus see that in a simulation of an interface without interface energy or with small interface energy the numerical efforts behave like

[TABLE]

Since the time step in a simulation must be decreased when the grid spacing $h$ in $x$–direction is decreased, the number $p$ can be larger than $4$ in a three dimensional simulation. From (7.5) we thus see that the numerical effort for the Allen-Cahn model grows much faster for the Allen-Cahn model than for the hybrid model when the required accuracy is increased. This confirms the rule stated above for the usage of both models in simulations.

This picture does not change essentially when the term $s_{10}$ vanishes identically. In this case the same considerations show that instead of (7.3) and (7.4) we would have $B_{\rm AC}=O({\cal E}_{\rm max}^{3/2})$ and $e_{\rm num}^{\rm AC}\geq C{\cal E}_{\rm max}^{-\frac{3}{2}p}$, hence the numerical effort for the Allen-Cahn model would still grow faster than for the hybrid model. However, a close investigation of the terms, which constitute $s_{10}$ and which are computed in [7], shows that only in very exceptional situations one can expect that $s_{10}$ vanishes identically.