An efficient and convergent finite element scheme for Cahn--Hilliard equations with dynamic boundary conditions

TL;DR

This paper introduces an efficient, unconditionally energy-stable finite element scheme for the Cahn--Hilliard equation with dynamic boundary conditions, demonstrating convergence and practical effectiveness through simulations.

Contribution

It proposes a novel finite element method for a recent Cahn--Hilliard model with boundary dynamics, establishing convergence and energy stability.

Findings

Scheme is unconditionally energy stable.

Convergence towards weak solutions is proven.

Numerical simulations confirm practicality and convergence order.

Abstract

The Cahn--Hilliard equation is a widely used model that describes amongst others phase separation processes of binary mixtures or two-phase flows. In the recent years, different types of boundary conditions for the Cahn--Hilliard equation were proposed and analyzed. In this publication, we are concerned with the numerical treatment of a recent model which introduces an additional Cahn--Hilliard type equation on the boundary as closure for the Cahn--Hilliard equation in the domain [C. Liu, H. Wu, Arch. Ration. Mech. An., 2019]. By identifying a mapping between the phase-field parameter and the chemical potential inside of the domain, we are able to postulate an efficient, unconditionally energy stable finite element scheme. Furthermore, we establish the convergence of discrete solutions towards suitable weak solutions of the original model. This serves also as an additional pathway to establish existence of weak solutions. Furthermore, we present simulations underlining the practicality of the proposed scheme and investigate its experimental order of convergence.

Tables

Table 1. Table 1. Parameters used in Section 5.1 .

$m$	$\delta$	$\sigma$	$m_{\Gamma }$	${\delta }_{\Gamma }$	$\kappa$	$C_{\text{pen}}$
0.01	0.01	2	0.01	0.02	1	250

Table 2. Table 2. Parameters used in Section 5.2 .

$m$	$\delta$	$\sigma$	$m_{\Gamma }$	${\delta }_{\Gamma }$	$\kappa$	$C_{\text{pen}}$
0.01	0.02	2	0.02	0.02	1	250

Table 3. (a)

$h$	${\text{err}}_h^{\Omega }$	${\text{EOC}}_h^{\Omega }$
$\sqrt{2}\cdot 2^{-6}$	$3.1\cdot {10}^{-2}$	-
$\sqrt{2}\cdot 2^{-7}$	$6.3\cdot {10}^{-3}$	2.3

Table 4. (a)

$h$	${\text{err}}_h^{\Omega }$	${\text{EOC}}_h^{\Omega }$
$\sqrt{2}\cdot 2^{-6}$	$3.1\cdot {10}^{-2}$	-
$\sqrt{2}\cdot 2^{-7}$	$6.3\cdot {10}^{-3}$	2.3

Table 5. (b)

$h$	${\text{err}}_h^{\Gamma }$	${\text{EOC}}_h^{\Gamma }$
$\sqrt{2}\cdot 2^{-6}$	$1.8\cdot {10}^{-1}$	-
$\sqrt{2}\cdot 2^{-7}$	$8.3\cdot {10}^{-2}$	1.1

Table 6. (a)

$\tau$	${\text{err}}_{\tau }^{\Omega }$	${\text{EOC}}_{\tau }^{\Omega }$
$4\cdot {10}^{-5}$	$1.4\cdot {10}^{-2}$	-
$2\cdot {10}^{-5}$	$4.5\cdot {10}^{-3}$	1.6

Table 7. (a)

$\tau$	${\text{err}}_{\tau }^{\Omega }$	${\text{EOC}}_{\tau }^{\Omega }$
$4\cdot {10}^{-5}$	$1.4\cdot {10}^{-2}$	-
$2\cdot {10}^{-5}$	$4.5\cdot {10}^{-3}$	1.6

Table 8. (b)

$\tau$	${\text{err}}_{\tau }^{\Gamma }$	${\text{EOC}}_{\tau }^{\Gamma }$
$4\cdot {10}^{-5}$	$7.3\cdot {10}^{-2}$	-
$2\cdot {10}^{-5}$	$2.4\cdot {10}^{-2}$	1.6

Equations

$${\partial }_t\varphi$$

$$\mu$$

$$\nabla \mu \cdot \mathbf{n}$$

$${\partial }_t\varphi$$

$${\mu }_{\Gamma }$$

$${\partial }_t\varphi$$

$$\mu$$

$$\nabla \mu \cdot \mathbf{n}$$

$$\nabla \varphi \cdot \mathbf{n}$$

$${\mathcal{E}}_{\Omega }(\varphi ):=\sigma \delta {\int }_{\Omega }\frac{1}{2}{\left|\nabla \varphi \right|}^2+\sigma {\delta }^{-1}{\int }_{\Omega }F(\varphi ),$$

$$W_{\log }(\varphi ):=\frac{\theta }{2}(1+\varphi )\log (1+\varphi )+\frac{\theta }{2}(1-\varphi )\log (1-\varphi )-\frac{{\theta }_c}{2}{\varphi }^2$$

$$W_{\mathrm{obst}}(\varphi ):=\{\begin{array}{cc}\theta (1-{\varphi }^2)& \varphi \in [-1,+1]\\ \infty & \text{else}\end{array}$$

$$\displaystyle\left.\kern-1.2pt\mathcal{E}_{\Omega}{(\phi)}\vphantom{\big{|}}\right|_{T}+\int_{0}^{T}\!\!\!\!\int_{\Omega}m\left|{\nabla\mu}\right|^{2}=\left.\kern-1.2pt\mathcal{E}_{\Omega}{(\phi)}\vphantom{\big{|}}\right|_{0}\,,$$

$${\mathcal{E}}_{\Gamma }(\varphi ):=\kappa {\delta }_{\Gamma }{\int }_{\Gamma }\frac{1}{2}{\left|{\nabla }_{\Gamma }\varphi \right|}^2+{\delta }_{\Gamma }^{-1}{\int }_{\Gamma }G(\varphi )$$

$${\partial }_t\varphi$$

$${\mu }_{\Gamma }$$

$$\sigma \delta \nabla \varphi \cdot \mathbf{n}=-\frac{1}{m_{\Gamma }}{\partial }_t\varphi -G^{\prime }(\varphi ),$$

$$\displaystyle\left.\kern-1.2pt\mathcal{E}_{\Omega}{(\phi)}\vphantom{\big{|}}\right|_{T}+\left.\kern-1.2pt\mathcal{E}_{\Gamma}{(\phi)}\vphantom{\big{|}}\right|_{T}+\int_{0}^{T}{\left[\int_{\Omega}m\left|{\nabla\mu}\right|^{2}+\int_{\partial\Omega}\tfrac{1}{m_{\Gamma}}\left|{\partial_{t}\phi}\right|^{2}\right]}=\left.\kern-1.2pt\mathcal{E}_{\Omega}{(\phi)}\vphantom{\big{|}}\right|_{0}+\left.\kern-1.2pt\mathcal{E}_{\Gamma}{(\phi)}\vphantom{\big{|}}\right|_{0}\,$$

$${\partial }_t\varphi$$

$$\mu$$

$$\displaystyle\left.\kern-1.2pt\mathcal{E}_{\Omega}{(\phi)}\vphantom{\big{|}}\right|_{T}+\left.\kern-1.2pt\mathcal{E}_{\Gamma}{(\phi)}\vphantom{\big{|}}\right|_{T}+\int_{0}^{T}{\left[\int_{\Omega}m\left|{\nabla\mu}\right|^{2}+\int_{\partial\Omega}m_{\Gamma}\left|{\mu}\right|^{2}\right]}=\left.\kern-1.2pt\mathcal{E}_{\Omega}{(\phi)}\vphantom{\big{|}}\right|_{0}+\left.\kern-1.2pt\mathcal{E}_{\Gamma}{(\phi)}\vphantom{\big{|}}\right|_{0}\,$$

$${\partial }_t\varphi$$

$$\mu$$

$$\displaystyle\left.\kern-1.2pt\mathcal{E}_{\Omega}{(\phi)}\vphantom{\big{|}}\right|_{T}+\left.\kern-1.2pt\mathcal{E}_{\Gamma}{(\phi)}\vphantom{\big{|}}\right|_{T}+\int_{0}^{T}{\left[\int_{\Omega}m\left|{\nabla\mu}\right|^{2}+\int_{\partial\Omega}m_{\Gamma}\left|{\nabla_{\Gamma}\mu}\right|^{2}\right]}=\left.\kern-1.2pt\mathcal{E}_{\Omega}{(\phi)}\vphantom{\big{|}}\right|_{0}+\left.\kern-1.2pt\mathcal{E}_{\Gamma}{(\phi)}\vphantom{\big{|}}\right|_{0}\,.$$

$${\int }_0^T{\int }_{\Omega }{\partial }_t\varphi \theta$$

$${\int }_0^T{\int }_{\Gamma }{\partial }_t\varphi \tilde{\theta }+$$

$$\begin{array}{rl}{\int }_0^T{\int }_{\Omega }\mu \hat{\theta }+{\int }_0^T{\int }_{\Gamma }{\mu }_{\Gamma }\hat{\theta }=& \delta \sigma {\int }_0^T{\int }_{\Omega }\nabla \varphi \cdot \nabla \hat{\theta }+{\delta }^{-1}\sigma {\int }_0^T{\int }_{\Omega }{F}^{\prime }(\varphi )\hat{\theta }\\ & +\kappa {\delta }_{\Gamma }{\int }_0^T{\int }_{\Gamma }{\nabla }_{\Gamma }\varphi \cdot {\nabla }_{\Gamma }\hat{\theta }+{\delta }_{\Gamma }^{-1}{\int }_0^T{\int }_{\Gamma }{G}^{\prime }(\varphi )\hat{\theta }\end{array}$$

$$X_{\kappa }:=\{\begin{array}{cc}\{v\in H^1(\Omega ):\gamma (v)\in H^1(\Gamma )\}& \text{ifleavevmodenobreak }\kappa >1,\\ H^1(\Omega )& \text{ifleavevmodenobreak }\kappa =0,\end{array}$$

$$\overline{\Omega }\equiv \bigcup _{K\in {\mathcal{T}}_h}\overline{K}$$

$$\forall K^{\Gamma }\in {\mathcal{T}}_h^{\Gamma }\text{leavevmode}\text{leavevmode}\exists !K\in {\mathcal{T}}_h\text{leavevmodenobreak suchleavevmodenobreak thatleavevmodenobreak }\overline{K^{\Gamma }}=\overline{K}\cap \Gamma ,$$

$$\Gamma \equiv \bigcup _{K^{\Gamma }\in {\mathcal{T}}_h^{\Gamma }}\overline{K^{\Gamma }}$$

$$U_h^{\Gamma }=\mathrm{span}\{\gamma ({\theta }_h):{\theta }_h\in U_h^{\Omega }\}.$$

$${\mathcal{I}}_h\left\{a\right\}:=\sum _{k=1}^{\dim U_h^{\Omega }}a({\mathbf{x}}_k){\chi }_{h,k},$$

$${\Vert (I-{\mathcal{I}}_h)\{f_h{g}_h\}\Vert }_{L^p(\Omega )}$$

$${\Vert (I-{\mathcal{I}}_h^{\Gamma })\{f_h{g}_h\}\Vert }_{L^p(\Gamma )}$$

$${\int }_K{\left|(I-{\mathcal{I}}_h)\{f_h{g}_h\}\right|}^p\le C{h}^{2p}{\int }_K{\left|f_h{g}_h\right|}_{W^{2,\infty }(K)}^p.$$

$$\begin{array}{rl}|f_h{g}_h{|}_{W^{2,\infty }(K)}=& \underset{i,j=1,...,d}{\max }{\Vert {\partial }_i{\partial }_j(f_h{g}_h)\Vert }_{L^{\infty }(K)}=\underset{i,j=1,...,d}{\max }{\Vert {\partial }_i{f}_h{\partial }_j{g}_h\Vert }_{L^{\infty }(K)}\\ \le & \underset{i,j=1,...,d}{\max }{\Vert {\partial }_i{f}_h\Vert }_{L^{\infty }(K)}{\Vert {\partial }_j{g}_h\Vert }_{L^{\infty }(K)}\end{array}$$

Full text

1. Introduction

Different approaches to model the hydrodynamics of fluid mixtures have been widely used in literature. In addition to the conventional sharp interface models which consist of separate hydrodynamic systems for each component of the mixture, there are diffuse interface models. In these models, the hyper-surface description of the fluid-fluid interface is replaced by a small transition region, where mixing of the macroscopically immiscible fluids is allowed. This leads to a smooth transition between the pure phases. In this manuscript, we analyze a numerical scheme for a diffuse interface model featuring dynamic boundary conditions, which was derived recently by C. Liu and H. Wu [38]. In an open domain $\Omega$ with boundary $\Gamma=\partial\Omega$ and outer normal vector $\boldsymbol{n}$, this model reads

[TABLE]

with suitable initial conditions for the phase-field parameter $\phi$ in $\overline{\Omega}$. Here, $\Delta_{\Gamma}$ denotes the Laplace–Beltrami operator, the positive parameters $m$ and $m_{\Gamma}$ are the mobility constants in the domain and on the boundary, the constant $\sigma>0$ is related to the surface tension, $\kappa\geq 0$ describes the influence of surface diffusion, and $\delta$ and $\delta_{\Gamma}$ prescribe the width of the transition area in the domain and on the boundary. The potentials $F$ and $G$ which govern the chemical potentials $\mu$ and $\mu_{\Gamma}$ will be discussed below in more detail. In (1.1), it is assumed that $\phi$ is defined on $\overline{\Omega}$ and that its evolution is governed by a chemical potential defined on $\Omega$ and an additional one defined on $\Gamma=\partial\Omega$. In contrast to other approaches (cf. [28]), $\mu$ and $\mu_{\Gamma}$ are distinct quantities which are only coupled via the normal derivative of $\phi$.

In the recent years, various different boundary conditions for Cahn–Hilliard equations have been discussed. The easiest form of a diffuse interface model for two-phase flow reads

[TABLE]

in combination with initial conditions for $\phi$. The chemical potential $\mu$ is given as the first variation of the free energy

[TABLE]

where $F$ is a double-well potential with minima in $\phi=\pm 1$ representing the pure fluid phases. Typical choices for $F$ are the logarithmic double-well potential

[TABLE]

with $0<\theta<\theta_{c}$, the double obstacle potential

[TABLE]

and the polynomial double-well potential $W_{\operatorname{pol}}{(\phi)}:=\tfrac{1}{4}{(\phi^{2}-1)}^{2}$. Cahn–Hilliard equations with polynomial double-well potentials are investigated, e.g. in [20, 58, 46, 5, 29, 49]. For Cahn–Hilliard equations with the singular potentials $W_{\operatorname{log}}$ and $W_{\operatorname{obst}}$, we refer the reader to [7, 18, 1, 13].

The boundary condition (1.2c) states that there is no flux across $\Gamma$, i.e. $\int_{\Omega}\phi$ is conserved. The second boundary condition (1.2d) indicates that the fluid-fluid interface, i.e. the zero level set of $\phi$, intersects the boundary $\Gamma$ at a static contact angle of $\tfrac{\pi}{2}$. This can be interpreted as neglecting the interactions between the fluids and the walls of the surrounding container. Although (1.2) satisfies the energy balance equation

[TABLE]

the boundary condition (1.2d) imposing a static contact angle is considered a major flaw and there are several attempts to improve this boundary condition. For an improved description of the occurring boundary effects, the surface energy functional

[TABLE]

was introduced with $\kappa\geq 0$ and $\delta_{\Gamma}>0$. Here, $\nabla_{\Gamma}$ denotes the surface gradient operator on $\Gamma$ and $G$ denotes a suitable boundary potential which might be chosen similar to $F$. Dynamic boundary conditions were picked in a way that the total energy ${[\mathcal{E}_{\Omega}+\mathcal{E}_{\Gamma}]}$ is decreasing in time. An example for such boundary conditions are Allen–Cahn-type boundary conditions (cf. [21, 22, 48, 55, 14, 27, 42, 13, 37, 16, 15, 41]), where (1.2d) is replaced by

[TABLE]

In the case $\kappa=0$ and $\delta_{\Gamma}=1$, this boundary conditions reduces to

[TABLE]

where the potential $G$ interpolates between the liquid-solid interfacial energies of the two fluid phases and prescribes the stationary contact angle via Young’s formula. This boundary condition was used in [47] to describe dynamic contact angles (see also [52]). Models combining the boundary conditions (1.2c) and (1.8) satisfy an energy equality of the form

[TABLE]

and conserve the mean value of $\phi$ in $\Omega$, i.e. we have \left.\kern-1.2pt\int_{\Omega}\phi\vphantom{\big{|}}\right|_{T}=\left.\kern-1.2pt\int_{\Omega}\phi\vphantom{\big{|}}\right|_{0}.

In [24] a coupled boundary condition replacing (1.2c) and (1.2d) was proposed. Assuming that $\mu_{\Gamma}$ is the trace of $\mu$, boundary conditions of the form

[TABLE]

were introduced. A solution to (1.2a), (1.2b), and (1.11) still minimizes $\mathcal{E}_{\Omega}+\mathcal{E}_{\Gamma}$ in the sense that

[TABLE]

holds true. However, the boundary conditions (1.11) do not allow for conservation of $\phi$.

A third approach was proposed by Goldstein et al. in [28]. In this publication, it was also assumed that $\mu_{\Gamma}$ is the trace of $\mu$. However, (1.11) was replaced by a Cahn–Hilliard-type boundary equations of the form

[TABLE]

Solutions to (1.2a), (1.2b), and (1.13) satisfy the energy equality

[TABLE]

Furthermore, the total mass of $\phi$ is conserved, i.e. $\partial_{t}{[\int_{\Omega}\phi+\int_{\Gamma}\phi]}=0$. A similar model was also discussed in [43].

A fourth approach was proposed by C. Liu and H. Wu. In [38], they derived model (1.1) via a variational approach using different flow maps for $\Omega$ and $\Gamma$. This model also satisfies (1.14). However, the underlying assumptions on the occurring boundary effects are different. While the boundary conditions proposed in [28] allow for mass transfer between $\Omega$ and $\Gamma$ and enforce an instant equilibration of the chemical potentials, i.e. $\mu_{\Gamma}$ has to be the trace of $\mu$, the model derived in [38] allows for differences in the chemical potentials, but prohibits mass transfer, i.e. $\int_{\Omega}\phi$ and $\int_{\Gamma}\phi$ are conserved individually. For boundary conditions interpolating between (1.13) and (1.1c)-(1.1e), which can be interpreted as non-instantaneous adsorption processes, we refer the reader to [34].

A first existence and uniqueness result for weak and strong solutions to (1.1) was provided in [38] by constructing solutions to a regularized system, where (1.1b) and (1.1e) are extended by $\alpha\partial_{t}\phi$, and taking the limit $\alpha\searrow 0$.

A different pathway to proving the existence of solutions to (1.1) was used by Garcke and Knopf in [26]. They interpreted model (1.1) as a gradient flow equation to the total free energy $\mathcal{E}_{\Omega}{(\phi)}+\mathcal{E}_{\Gamma}{(\phi)}$ and used this structure for their proof of existence and uniqueness of weak solutions. A comparable wellposedness result derived from a fully discrete finite element scheme can be found in Theorem 4.4 of this manuscript.

The numerical treatment of the Cahn–Hilliard equation and its variants – often in combination with Navier–Stokes-equations – was intensely discussed through the last years. Consequently, there are various different discretization techniques at hand, which transfer the energy stability (1.6) to a discrete setting. These techniques include approaches based on convex-concave splittings of the energy (cf. [54, 50]) or the polynomial double-well potential (cf. [19] and [32, 30, 25, 31] for an application of Cahn–Hilliard–Navier–Stokes-systems), stabilized linearly implicit approaches (cf. [56, 51]), the method of invariant energy quadratization (cf. [12, 57]) and the recently developed scalar auxiliary variable approach (see [36]).

Although, we will restrict ourselves to non-singular potentials, we do not want to conceal that there are also numerical schemes at hand which are able to deal with the singular potentials $W_{\operatorname{log}}$ and $W_{\operatorname{obst}}$ (see e.g. [17, 8, 6, 3, 4, 23]).

In this publication, we are interested in the numerical treatment of (1.1). A finite difference model for the treatment of the Allen–Cahn-type boundary conditions (1.8) was proposed in [33]. A first finite element scheme for model (1.1) was proposed in the Bachelor’s thesis [53] (see also [26] for numerical results). In this thesis, a straightforward, fully implicit discretization based on continuous, piecewise linear finite element functions was applied to model (1.1), and the arising nonlinear system was solved using Newton’s method. In this publication, we pursue a different approach and investigate the connection between $\phi$ and the chemical potentials.

The peculiarity of (1.1) is the coupling between the chemical potential $\mu$ defined inside of the domain and the $\mu_{\Gamma}$ on the boundary. In the standard Cahn–Hilliard equation (1.2), the chemical potential is merely a definition in terms of $\phi$. This allows us to write (1.2) as a sole, nonlinear, fourth-order equation (see e.g. [32, 31, 25]). In (1.1), however, the chemical potentials $\mu$ and $\mu_{\Gamma}$ are coupled via the normal derivative $\nabla\phi\cdot\boldsymbol{n}$. Therefore, its weak form formally reads

[TABLE]

with sufficiently regular $\theta$, $\tilde{\theta}$, and $\hat{\theta}$. In particular, we have only one equation for $\mu$ and $\mu_{\Gamma}$. Consequently, the chemical potentials have to be determined by solving a system consisting of (1.15c) and the additional assumption that $\phi$ is continuous on $\overline{\Omega}$. The latter one translates to the constraint that (1.15a) and (1.15b) yield compatible results. Deducing a suitable expression for $\mu$ will be key ingredient for the derivation of an efficient numerical scheme, but also for the numerical analysis, as the existence of a unique (discrete) $\mu$ for any given $\phi$ allows us to reuse techniques from the analysis of the standard Cahn–Hilliard equations. As we will discuss in Remark 2.5, this approach also prevents the arising linear system from degenerating for vanishing time increments.

The outline of the paper is as follows. In Section 2, we introduce the discrete function spaces and derive the discrete scheme. In Section 3, we will establish a first a priori estimate which is discrete counterpart of (1.14), and use this estimate to prove the existence of discrete solutions. The main convergence result, Theorem 4.4, which also provides the existence of weak solutions, can be found in Section 4, where we establish improved regularity results and show the convergence of discrete solutions towards weak solutions of (1.1). For uniqueness results for these weak solutions, we refer the reader to Section 5 in [26]. We will conclude this section by briefly discussing the case of Allen–Cahn-type boundary conditions (cf. Remark 4.5). By showing that the presented techniques are also applicable for Allen–Cahn-type boundary conditions, we also cover (1.8) and its special case (1.9) suggested in [47]. In Section 5, we present numerical simulations of phase-separation processes to underline the practicality of the scheme. We shall also validate our scheme in terms of mass conservation, energy dissipation, and compatibility of (1.15a) and (1.15b).

Notation

Given the spatial domain $\Omega\subset\mathds{R}^{d}$ with $d\in{\{2,3\}}$ and a time interval ${(0,T)}$, we denote the space-time cylinder $\Omega\times{(0,T)}$ by $\Omega_{T}$. By $W^{k,p}{(\Omega)}$ we denote the space of $k$-times weakly differentiable functions with weak derivatives in $L^{p}{(\Omega)}$. The symbol $W^{k,p}_{0}{(\Omega)}$ stands for the closure of $C^{\infty}_{0}{(\Omega)}$ in $W^{k,p}{(\Omega)}$. For $p=2$, we will denote $W^{k,2}{(\Omega)}$ by $H^{k}{(\Omega)}$ and $W^{k,2}_{0}{(\Omega)}$ by $H^{k}_{0}{(\Omega)}$. The dual space of $H^{1}{(\Omega)}$ will be denoted by ${(H^{1}{(\Omega)})}^{\prime}$ and the corresponding dual pairing by ${\left\langle.,.\right\rangle}$. For a Banach space $X$ and a time interval $I$, the symbol $L^{p}{(I;X)}$ stands for the parabolic space of $L^{p}$-integrable functions on $I$ with values in $X$. We use a notation similar to the one introduced above for function spaces defined on $\Gamma$. In this publication, we are concerned with domains $\Omega$ having a lipschitzian boundary $\Gamma$. In this case, the spaces $L^{p}{(\Gamma)}$, $W^{1,p}{(\Gamma)}$, and $H^{1}{(\Gamma)}$ are well-defined (cf. [35]). We denote the dual pairing between ${(H^{1}{(\Gamma)})}^{\prime}$ and $H^{1}{(\Gamma)}$ by ${\left\langle.,.\right\rangle}_{\Gamma}$. In addition, we define the space

[TABLE]

where $\gamma$ defines the trace operator. For domains with lipschitzian boundaries, the trace operator is uniquely defined and lies in $\mathcal{L}{(W^{1,p}{(\Omega)},W^{1-1/p,p}{(\Gamma)})}$ (cf. [44]). For brevity, we will sometimes (in particular when the considered function is continuous) neglect the trace operator and write $v$ instead of $\gamma{(v)}$.

2. Derivation of an efficient numerical scheme

We start by introducing the general notation and the discretization techniques used in the considered scheme. Concerning the discretization with respect to time, we assume that

• •

the time interval $I:=[0,T)$ is subdivided in intervals $I_{n}:=[t_{n},t_{n+1})$ with $t_{n+1}=t_{n}+\tau_{n}$ for time increments $\tau_{n}>0$ and $n=0,...,N-1$ with $t_{N}=T$. For simplicity, we take $\tau_{n}\equiv\tau=\tfrac{T}{N}$ for $n=0,...,N-1$.

The spatial domain $\Omega\subset\mathds{R}^{d}$ in spatial dimensions $d\in{\{2,3\}}$ is assumed to be bounded and convex. To avoid additional technicalities, we will asssume that $\Omega$ is polygonal (or polyhedral, respectively). We introduce partitions $\mathcal{T}_{h}$ of $\Omega$ and $\mathcal{T}_{h}^{\Gamma}$ of $\Gamma$ depending on a spatial discretization parameter $h>0$ satisfying the following assumptions:

• •

Let ${\{\mathcal{T}_{h}\}}_{h>0}$ be a quasiuniform family (in the sense of [10]) of partitions of $\Omega$ into disjoint, open simplices $K$, so that

[TABLE]

• •

Let ${\{\mathcal{T}_{h}^{\Gamma}\}}_{h>0}$ be a quasiuniform family of partitions of $\Gamma$ into disjoint, open simplices $K^{\Gamma}$, so that

[TABLE]

and

[TABLE]

• ‣ 2 implies that $\mathcal{T}_{h}^{\Gamma}$ is compatible to $\mathcal{T}_{h}$ in the sense that all elements in $\mathcal{T}_{h}^{\Gamma}$ are edges (or faces) of elements in $\mathcal{T}_{h}$. For the approximation of the phase-field $\phi$ and the chemical potential $\mu$ we use continuous, piecewise linear finite element functions on $\mathcal{T}_{h}$. This space will be denoted by $U_{h}^{\Omega}$ and is spanned by the functions ${\{\chi_{h,k}\}}_{k=1,...,\dim U_{h}^{\Omega}}$ forming a dual basis to the vertices ${\{\boldsymbol{x}_{k}\}}_{k=1,...,\dim U_{h}^{\Omega}}$ of $\mathcal{T}_{h}$, i.e. $\chi_{h,k}{(\boldsymbol{x}_{k})}=\delta_{k,l}$ for $k,l=1,...,\dim U_{h}^{\Omega}$. Analogously, we denote the space of continuous, piecewise linear finite element functions on $\mathcal{T}_{h}^{\Gamma}$ by $U_{h}^{\Gamma}$. This space is spanned by functions ${\{\chi_{h,k}^{\Gamma}\}}_{k=1,...,\dim U_{h}^{\Gamma}}$ forming a dual basis to the vertices ${\{\boldsymbol{x}_{k}^{\Gamma}\}}_{k=1,...,\dim U_{h}^{\Gamma}}$ of $\mathcal{T}_{h}^{\Gamma}$, i.e. $\chi^{\Gamma}_{h,k}{(\boldsymbol{x}^{\Gamma}_{k})}=\delta_{k,l}$ for $k,l=1,...,\dim U_{h}^{\Omega}$. Due to the compatibility condition for $\mathcal{T}_{h}$ and $\mathcal{T}_{h}^{\Gamma}$, we have

[TABLE]

Without loss of generality, we may assume that the first $\dim U_{h}^{\Gamma}$ vertices of $\mathcal{T}_{h}$ are located on $\Gamma$, i.e. ${\{\boldsymbol{x}_{k}^{\Gamma}\}}_{k=1,...,\dim U_{h}^{\Gamma}}={\{\boldsymbol{x}_{k}\}}_{k=1,...,\dim U_{h}^{\Gamma}}$. We define the nodal interpolation operators $\mathcal{I}_{h}\,:\,C^{0}{(\overline{\Omega})}\rightarrow U_{h}^{\Omega}$ and $\mathcal{I}_{h}^{\Gamma}\,:\,C^{0}{(\overline{\Gamma})}\rightarrow U_{h}^{\Gamma}$ by

[TABLE]

For future reference, we state the following estimate for the interpolation operators.

Lemma 2.1**.**

Let $\mathcal{T}_{h}$ and $\mathcal{T}_{h}^{\Gamma}$ satisfy • ‣ 2 and • ‣ 2. Furthermore, let $p\in[1,\infty)$, $1\leq q\leq\infty$, and $q^{*}=\tfrac{q-1}{q}$ for $q<\infty$ or $q^{*}=1$ for $q=\infty$. Then

[TABLE]

holds true for all $f_{h},\,g_{h}\in U_{h}^{\Omega}$.

Proof.

Using the standard error estimates for the nodal interpolation operator (cf. [10]), we compute on each $K\in\mathcal{T}_{h}$:

[TABLE]

As $f_{h},\,g_{h}\in U_{h}^{\Omega}$, they are linear on each $K$, i.e. their second spatial derivatives vanish. Therefore, we obtain

[TABLE]

As the first spatial derivatives of $f_{h}$ and $g_{h}$ are constant, we combine (2.5) and (2.6) and apply Hölder’s inequality to obtain

[TABLE]

Similar computations provide the result for $\mathcal{I}_{h}^{\Gamma}$. ∎

Concerning the potentials $F$ and $G$, we make the following assumptions:

• •

$F,G\in C^{1}{(\mathds{R})}$ are bounded from below, i.e. there exists a constant $C>0$ such that $F{(s)}>-C$ and $G{(s)}>-C$ for all $s\in\mathds{R}$. Furthermore, there exist convex, non-negative functions $F_{+},\,G_{+}\in C^{1}{(\mathds{R})}$ and concave functions $F_{-},\,G_{-}\in C^{1}{(\mathds{R})}$ such that $F\equiv F_{+}+F_{-}$ and $G\equiv G_{+}+G_{-}$.

• •

The convex and concave parts of $F$ and $G$ can be further decomposed into a polynomial part of degree four and an additional part with a globally Lipschitz-continuous first derivative. Moreover, there exists $\beta\geq 0$ such that the concave parts satisfy

[TABLE]

for all $s_{1},\,s_{2}\in\mathds{R}$. In the case $\kappa=0$, we assume that the above assumption holds true for $\beta>0$.

Remark 2.2**.**

The Assumptions • ‣ 2 and • ‣ 2 are in particular satisfied by the polynomial double-well potential $W_{\operatorname{pol}}{(\phi)}:=\tfrac{1}{4}{(1-\phi^{2})}^{2}$ and the penalized double-well potential

[TABLE]

*The latter one is often used in practical computations, as it penalizes $\phi\notin{[-1,+1]}$ but does not introduce singularities (cf. [31]).

The logarithmic potential $W_{\operatorname{log}}$ and the double obstacle potential $W_{\operatorname{obst}}$, which have the advantage of restricting $\phi$ to the interval ${[-1,+1]}$, do not satisfy • ‣ 2 and • ‣ 2 and are therefore not considered in this manuscript.*

Remark 2.3**.**

In this publication, we consider only a convex-concave decomposition of the double-well potential. Other suitable, energy stable discretization techniques can be found in [31]. For a comparison of these techniques, we refer to [40].

Using the notation introduced above and the compatibility condition (2.1), we may write our discrete scheme as follows: For given $\phi_{h}^{n-1}\in U_{h}^{\Omega}$, find ${(\phi_{h}^{n},\,\mu_{h}^{n}\,,\mu_{\Gamma,h}^{n})}\in U_{h}^{\Omega}\times U_{h}^{\Omega}\times U_{h}^{\Gamma}$ satisfying

[TABLE]

for all $\theta_{h}\in U_{h}^{\Omega}$. As discussed on the example of the weak formulation (1.15), (2.9) only provides on equation for both chemical potentials. Consequently, the goal for this section will be to derive an equivalent formulation for (2.9) with unique expressions for $\mu_{h}^{n}$ and $\mu_{\Gamma,h}^{n}$, which allows us to reuse the standard techniques established for (1.2).

We define the lumped mass matrices $\mathbf{M}_{\Omega}$ and $\mathbf{M}_{\Gamma}$ via

[TABLE]

Furthermore, we collect the nodal values of $\phi_{h}^{n}$, $\phi_{h}^{n-1}$, $\mu_{h}^{n}$, and $\mu_{\Gamma,h}^{n}$ in the vectors $\Phi^{n}$, $\Phi^{n-1}$, $P^{n}$, and $P_{\Gamma}^{n}$. In a slight misuse of notation, we will write $F{(\Phi^{n})}$, when we apply a function $F$ to all components of $\Phi^{n}$. With this notation, we are able to rewrite (2.9) as

[TABLE]

Here, we used the extension operator \left.\kern-1.2pt.\vphantom{\big{|}}\right|^{\Omega}\,:\,\mathds{R}^{\dim U_{h}^{\Gamma}}\rightarrow\mathds{R}^{\dim U_{h}^{\Omega}} defined via

[TABLE]

and the restriction operator \left.\kern-1.2pt.\vphantom{\big{|}}\right|_{\Gamma}\,:\,\mathds{R}^{\dim U_{h}^{\Omega}}\rightarrow\mathds{R}^{\dim U_{h}^{\Gamma}}, which restricts a vector its first $\dim U_{h}^{\Gamma}$ entries.

In order to derive a scheme allowing to solve (2.11), we define restriction operators for matrices. In particular, we will split a matrix $\mathbf{A}\in\mathds{R}^{\dim U_{h}^{\Omega}\times\dim U_{h}^{\Omega}}$ into submatrices

[TABLE]

such that

[TABLE]

Hence, the chemical potentials are given as solutions of the ${(\dim U_{h}^{\Omega}+\dim U_{h}^{\Gamma})}\times{(\dim U_{h}^{\Omega}+\dim U_{h}^{\Gamma})}$-system

[TABLE]

[TABLE]

Here, the first two lines are a consequence of (2.11c) and the last line guarantees that (2.11a) and (2.11b) provide the same result for \left.\kern-1.2pt\Phi^{n}\vphantom{\big{|}}\right|_{\Gamma}.

As the (2.11) is nonlinear in $\Phi^{n}$, computing a possible solution requires the application of an iterative scheme (e.g. Newton’s method) and therefore solving (2.15) multiple times per time step. Hence, solving a ${(\dim U_{h}^{\Omega}+\dim U_{h}^{\Gamma})}\times{(\dim U_{h}^{\Omega}+\dim U_{h}^{\Gamma})}$-system each time is not desirable and we have to continue reducing the complexity of the system.

From the second line in (2.15), we immediately get

[TABLE]

while the first line provides

[TABLE]

Using (2.17) and (2.18), we may write the last line in (2.15) as

[TABLE]

and therefore

[TABLE]

As $\mathbf{M}_{\Omega}^{-1}$ is a diagonal matrix, \left.\kern-1.2pt\mathbf{M}_{\Omega}^{-1}\vphantom{\big{|}}\right|_{\Gamma\times\Omega}\left.\kern-1.2pt\mathbf{L}_{\Omega}\vphantom{\big{|}}\right|_{\Omega\times\Gamma}=\left.\kern-1.2pt\mathbf{M}_{\Omega}^{-1}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}\left.\kern-1.2pt\mathbf{L}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma} holds true. This allows us to multiply (2.20) by \left.\kern-1.2pt\mathbf{M}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma} to obtain

[TABLE]

In order to show that (2.21) provides a well-defined expression for \left.\kern-1.2ptP^{n}\vphantom{\big{|}}\right|_{\Gamma}, we need to prove that the matrix on the left-hand side is indeed invertible.

Lemma 2.4**.**

The matrix {\left(m\left.\kern-1.2pt\mathbf{L}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}+m_{\Gamma}\left.\kern-1.2pt\mathbf{M}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}\mathbf{M}_{\Gamma}^{-1}\mathbf{L}_{\Gamma}\mathbf{M}_{\Gamma}^{-1}\left.\kern-1.2pt\mathbf{M}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}\right)}, that is defined via (2.10) and (2.13), is symmetric, positive definite.

Proof.

It is obvious that m\left.\kern-1.2pt\mathbf{L}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma} and m_{\Gamma}\left.\kern-1.2pt\mathbf{M}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}\mathbf{M}_{\Gamma}^{-1}\mathbf{L}_{\Gamma}\mathbf{M}_{\Gamma}^{-1}\left.\kern-1.2pt\mathbf{M}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma} are symmetric, positive semi-definite matrices. Therefore, it will be sufficient to show that A^{T}\left.\kern-1.2pt\mathbf{L}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}A>0 for all $0\neq A\in\mathds{R}^{\dim U_{h}^{\Gamma}}$ to complete the proof. This is equivalent to showing

[TABLE]

From (2.10), we have that $\mathbf{L}_{\Omega}$ is symmetric, positive semi-definite with only constant vectors corresponding to the zero eigenvalue. As the restrictions in (2.22) do not allow for constant vectors, the proof is complete. ∎

Combining (2.21) with (2.17), we obtain an expression for the chemical potential which requires us to solve only a $\dim U_{h}^{\Gamma}$ by $\dim U_{h}^{\Gamma}$ linear system with a sparse, symmetric, positive definite matrix. Having an expression for the chemical potential, we propose the following nonlinear equation for $\Phi^{n}$. For given $\Phi^{n-1}\in\mathds{R}^{\dim U_{h}^{\Omega}}$, we compute $\Phi^{n}\in\mathds{R}^{\dim U_{h}^{\Omega}}$ satisfying

[TABLE]

Here, $\operatorname{R}_{\Gamma}$ and $\operatorname{R}_{\stackrel{{\scriptstyle\circ}}{{\Omega}}}$, which are defined in (2.16), also depend on the known values $\Phi^{n-1}$. As we will show in the next section, solutions $\Phi^{n}$ to (2.23) satisfy the compatibility condition used in (2.15), which allows us to recover (2.9).

Remark 2.5**.**

*At this point, we want to discuss the advantages of formulation (2.23). Although (2.9) can be written in a symmetric form, trying to solve (2.9) directly has one major flaw. As there is no explicit formula for the chemical potentials available, one has to solve for $\phi_{h}^{n}$, $\mu_{h}^{n}$, and $\mu_{\Gamma,h}^{n}$ monolithically. However, this system degenerates for $\tau\searrow 0$, i.e. if we opt for a small time increment to capture rapid changes, we end up with an ill-conditioned system: For $\tau\searrow 0$, (2.9a) and (2.9b) reduce to $\phi_{h}^{n}=\phi_{h}^{n-1}$, i.e. the dependencies on the chemical potentials vanishe leaving only (2.9c) to determine both potentials.

The proposed scheme (2.23), on the other hand, is based on another $\tau$-independent relation between the chemical potentials, which prevents the system from becoming ill-conditioned for vanishing time increments. For an illustration of the dependence of the condition numbers on the size of the time increment based on practical computations, we refer the reader to Section 5.

The downside of (2.23) is that it requires us to solve an additional smaller linear system with the matrix {\left(m\left.\kern-1.2pt\mathbf{L}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}+m_{\Gamma}\left.\kern-1.2pt\mathbf{M}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}\mathbf{M}_{\Gamma}^{-1}\mathbf{L}_{\Gamma}\mathbf{M}_{\Gamma}^{-1}\left.\kern-1.2pt\mathbf{M}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}\right)} repeatedly. However, as this matrix is symmetric and positive definite, it can be tackled efficiently using a conjugate gradient method.*

3. Stability and existence of discrete solutions

In this section, we analyze the discrete scheme (2.23) proposed in the previous section. As we derived explicit expressions for $\mu_{h}^{n}$ and $\mu_{\Gamma,h}^{n}$ in the previous section, we could return to the variational form (2.9) and derive stability and existence results from there. Namely, testing (2.9a) by $\mu_{h}^{n}$, (2.9b) by $\mu_{\Gamma,h}^{n}$, and (2.9c) by ${(\phi_{h}^{n}-\phi_{h}^{n-1})}$ will provide a discrete version of (1.14). However, as (2.23) is the formulation we suggest to implement, we establish first stability and existence results based on this formulation. By doing so, we shall verify that all information from (2.9) are preserved in (2.23) and shed light on the structure of (2.23).

Although (2.23) is entirely written in terms of the unknown quantity $\Phi^{n}$, we will continue using $P^{n}$ and $P_{\Gamma}^{n}$, which are defined in (2.21), (2.17), and (2.18), to simplify the notation. For the ease of representation, we will set $\sigma=\delta=\delta_{\Gamma}=1$ for the remainder of this publication. As a first step, we shall verify that (2.23) indeed satisfies the compatibility constraint m\left.\kern-1.2pt\mathbf{M}_{\Omega}^{-1}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}\left.\kern-1.2pt\mathbf{L}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Omega}P^{n}=m_{\Gamma}\mathbf{M}_{\Gamma}^{-1}\mathbf{L}_{\Gamma}P_{\Gamma}^{n}. This auxiliary result allows us derive an a priori stability result for (2.23) which serves as the corner stone for proving the existence of discrete solutions without additional restrictions on $h$ or $\tau$.

Lemma 3.1**.**

Let $P^{n}$ and $P_{\Gamma}^{n}$ be defined via (2.21), (2.17), and (2.18). Then the identity

[TABLE]

holds true.

Proof.

Using (2.17), we compute

[TABLE]

Recalling (2.21) and (2.18), we obtain

[TABLE]

which completes the proof. ∎

This result allows us to show that the phasefield parameter is conserved in $\Omega$ and on $\Gamma$. Multiplying (2.23) by $\boldsymbol{1}^{T}\mathbf{M}_{\Omega}$ and by \left.\kern-1.2pt\boldsymbol{1}_{\Gamma}^{T}\mathbf{M}_{\Gamma}\vphantom{\big{|}}\right|^{\Omega} proves the following corollary.

Corollary 3.2**.**

Let $\Phi^{n}$ be a discrete solution of (2.23). Then

[TABLE]

with $\boldsymbol{1}:={(1,...,1)}^{T}\in\mathds{R}^{\dim U_{h}^{\Omega}}$ and \boldsymbol{1}_{\Gamma}:=\left.\kern-1.2pt\boldsymbol{1}\vphantom{\big{|}}\right|_{\Gamma}.

Using the above auxiliary results, we are now able to state a first stability result which is a discrete version of the energy equality (1.14).

Lemma 3.3**.**

Let the assumptions • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 2, and • ‣ 2 hold true and let $\Phi^{n-1}\in\mathds{R}^{\dim U_{h}^{\Omega}}$ be given. Then a solution $\Phi^{n}\in\mathds{R}^{\dim U_{h}^{\Omega}}$ to (2.23), if it exists, satifies

[TABLE]

with $\boldsymbol{1}:={(1,...,1)}^{T}\in\mathds{R}^{\dim U_{h}^{\Omega}}$, \boldsymbol{1}_{\Gamma}:=\left.\kern-1.2pt\boldsymbol{1}\vphantom{\big{|}}\right|_{\Gamma}, and $P^{n}$ and $P_{\Gamma}^{n}$ defined in (2.17), (2.21), and (2.18).

Proof.

We multiply (2.23) by ${\left(\mathbf{M}_{\Omega}P^{n}+\begin{pmatrix}\mathbf{M}_{\Gamma}P_{\Gamma}^{n}\\ \boldsymbol{0}\end{pmatrix}\right)}$ and use Lemma 3.1 to obtain

[TABLE]

As $III$ and $IV$ provide the dissipative parts of the desired estimate, we have show to that $I$ and $II$ yield the time difference of the energy. Recalling (2.18), we compute

[TABLE]

Consequently, we obtain from (2.21)

[TABLE]

As $\mathbf{M}_{\Omega}$ and $\mathbf{M}_{\Gamma}$ are diagonal matrices, we may combine ${(F^{\prime}_{+}{(\Phi^{n})}+F^{\prime}_{-}{(\Phi^{n-1})})}$ and ${(\Phi^{n}-\Phi^{n-1})}$, and {(G^{\prime}_{+}{(\left.\kern-1.2pt\Phi^{n}\vphantom{\big{|}}\right|_{\Gamma})}+G^{\prime}_{-}{(\left.\kern-1.2pt\Phi^{n-1}\vphantom{\big{|}}\right|_{\Gamma})})} and \left.\kern-1.2pt{(\Phi^{n}-\Phi^{n-1})}\vphantom{\big{|}}\right|_{\Gamma} componentwise. In combination with $s_{1}{(s_{1}-s_{2})}=\tfrac{1}{2}s_{1}^{2}+\tfrac{1}{2}{(s_{1}-s_{2})}^{2}-\tfrac{1}{2}s_{2}^{2}$, this provides the result. ∎

Using the a priori estimate from Lemma 3.3, we are able to prove the existence of discrete solutions.

Lemma 3.4**.**

Let the assumptions • ‣ 2, • ‣ 2, • ‣ 2, and • ‣ 2 hold true and let $\Phi^{n-1}\in\mathds{R}^{\dim U_{h}^{\Omega}}$ be given. Then, there exists at least one vector $\Phi^{n}\in\mathds{R}^{\dim U_{h}^{\Omega}}$ solving (2.23).

Proof.

We will prove the existence of discrete solutions by contradiction. Let $\left|\!\left|\!\left|.\right|\!\right|\!\right|$ denote the discrete $L^{2}$-norm which is derived from the inner product ${(A,B)}:=A^{T}\mathbf{M}_{\Omega}B$. According to Corollary 3.2, the mean-value of the phase-field is conserved in $\Omega$. This allows us to assume w.l.o.g. that $\boldsymbol{1}^{T}\mathbf{M}_{\Omega}\Phi^{n}=\boldsymbol{1}^{T}\mathbf{M}_{\Omega}\Phi^{n-1}=0$. Therefore, $\sqrt{{\Phi^{n}}^{T}\mathbf{L}_{\Omega}\Phi^{n}}$ is also a norm of $\Phi^{n}$. Under the assumption that (2.23) has no solution in

[TABLE]

for any $R>0$, the function $\mathcal{H}$ defined via

[TABLE]

has no root and is continuous on $B_{R}$. This allows us to define a function $\mathcal{A}\,:B_{R}\rightarrow\partial B_{R}\subset B_{R}$ as

[TABLE]

As $\mathcal{A}$ is continuous and maps a closed set onto itself, Brouwer’s fixed point theorem provides the existence of at least one fixed point $\Phi^{*}$. In the following, we will show

[TABLE]

for a suitable $\Psi\in B_{R}$ and $R$ large enough. This contradiction shows that the initial assumption of (2.23) not having solutions in $B_{R}$ is wrong. To prove the contradiction (3.9), we choose $\Psi=\tilde{\Psi}_{1}+\tilde{\Psi}_{2}-\boldsymbol{1}^{T}\mathbf{M}_{\Omega}{(\tilde{\Psi}_{1}+\tilde{\Psi}_{2})}{(\boldsymbol{1}^{T}\mathbf{M}_{\Omega}\boldsymbol{1})}^{-1}\boldsymbol{1}$ with

[TABLE]

and

[TABLE]

i.e. the test vector is the sum of the chemical potentials deprived of their mean values. The computations from the proof of Lemma 3.3 provide

[TABLE]

where the constant $C$ depends on $\Phi^{n-1}$ and the lower bound from • ‣ 2, but not on the fixed point $\Phi^{*}$ or $R$. Since all norms on finite dimensional spaces are equivalent, there exists $c>0$ such that $\tfrac{1}{2}{\Phi^{*}}^{T}\mathbf{L}_{\Omega}\Phi^{*}\geq c{\Phi^{*}}^{T}\mathbf{M}_{\Omega}\Phi^{*}$ and we obtain

[TABLE]

for $R$ large enough. This provides the second inequality in (3.9). In order to establish the first inequality we again use the computations from the proof of Lemma 3.3 to show

[TABLE]

with $0<\varepsilon,\tilde{\varepsilon}<\!\!<1$. For every fixed $h$, there is a constant $C_{h}>0$ such that \left.\kern-1.2pt\Phi^{*}\vphantom{\big{|}}\right|_{\Gamma}^{T}\mathbf{M}_{\Gamma}\left.\kern-1.2pt\Phi^{*}\vphantom{\big{|}}\right|_{\Gamma}\leq C_{h}{\Phi^{*}}^{T}\mathbf{M}_{\Omega}\Phi^{*}. Hence, we have

[TABLE]

with $C_{\varepsilon,\tilde{\varepsilon}}>0$ independent of $\Phi^{*}$ and $R$. Choosing $\varepsilon$ and $\tilde{\varepsilon}$ small enough provides ${(c-\varepsilon-C_{h}\tilde{\varepsilon})}>0$. Hence, we obtain the first inequality in (3.9) for $R$ large enough, which completes the proof. ∎

Remark 3.5**.**

The existence result in Lemma 3.3 implies no constraints on the time increment $\tau$. Therefore, we have the existence of discrete solutions for arbitrary time increments.

4. Convergence of the discrete scheme

In this section, we show that the discrete solutions established in the last section converge towards suitable weak solutions of (1.1). This requires some assumptions on the initial data. In particular, we will assume that

• •

the initial data $\phi_{0}\in X_{\kappa}$ and its projection $\phi_{h}^{0}$ onto $U_{h}^{\Omega}$ satisfies

[TABLE]

with some $C>0$ independent of $h$ and $\tau$.

Furthermore, the regularity results provided in this section require additional assumptions on $h$ and $\tau$. In particular, we will need

• •

that $\tfrac{h^{4}}{\tau}\searrow 0$ for ${(h,\tau)}\searrow 0$ when $\kappa>0$ and that $\tfrac{h^{2}}{\tau}\searrow 0$ for ${(h,\tau)}\searrow 0$ when $\kappa=0$.

Assumption • ‣ 4 allows us to state our first regularity result.

Corollary 4.1**.**

*Let the assumptions • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 2, and • ‣ 4 hold true and let $h>0$ be small enough. Then a solution ${(\phi_{h}^{n},\mu_{h}^{n},\mu_{\Gamma,h}^{n})}_{n=1,...,N}$ to (2.9) satisfies *

[TABLE]

with a constant $C>0$ independent of $h$ and $\tau$.

Proof.

After summing the result of Lemma 3.3 over all time steps and recalling Corollary 3.2 and • ‣ 4, it remains to show that we have indeed control over the complete $H^{1}$ norm of $\mu_{h}^{n}$ and $\mu_{\Gamma,h}^{n}$. To establish this result, we will follow the lines of [26]. Testing (2.9c) by $\mathcal{I}_{h}{\left\{\eta\right\}}$ with $\eta\in C_{0}^{\infty}{(\Omega;{[0,1]})}$, which is not identically zero, we obtain

[TABLE]

From • ‣ 2 and • ‣ 4, we obtain

[TABLE]

Hence, there exists a constant $\tilde{C}{(\eta)}$ independent of $h$ and $\tau$ such that $\left|{\int_{\Omega}\mathcal{I}_{h}{\left\{\mu_{h}^{n}\eta\right\}}}\right|\leq\tilde{C}{(\eta)}$. We now define

[TABLE]

From standard error estimates for the interpolation operator $\mathcal{I}_{h}$ (cf. [10]), we derive the existence of $c{(\eta)}>0$ such that $\int_{\Omega}\mathcal{I}_{h}{\left\{\eta\right\}}\geq c{(\eta)}$ for $h$ small enough. Therefore, we may use the generalized Poincaré inequality (cf. [2]), which we cite in the appendix as Lemma A.1, with $u_{0}\equiv 0$ and $C_{0}:=\tilde{C}{(\eta)}/c{(\eta)}$ to obtain

[TABLE]

To obtain the $L^{2}$-bound for $\mu_{\Gamma,h}^{n}$, we test (2.9c) by $\theta_{h}\equiv 1$ and obtain

[TABLE]

Considerations similar to (4.2) show that the last term on the right-hand side is also bounded by a constant independent of $h$ and $\tau$. Therefore, we may use Poincaré’s inequality to complete the proof. ∎

In a second step, we derive uniform bounds for the time difference quotient $\partial_{\tau}^{-}\phi_{h}^{n}:=\tau^{-1}{(\phi_{h}^{n}-\phi_{h}^{n-1})}$ of the phase-field parameter on $\Omega$ and $\Gamma$.

Lemma 4.2**.**

Let the assumptions • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 4, and • ‣ 4 hold true. Furthermore, let $h>0$ be small enough such that Corollary 4.1 holds true. Then a solution ${(\phi_{h}^{n})}_{n=1,...,N}$ to (2.9) satisfies

[TABLE]

with $C>0$ independent of $h$ and $\tau$.

Proof.

We take $\theta\in H^{1}{(\Omega)}$ and test (2.9a) by $\theta_{h}:=\mathcal{P}_{U_{h}^{\Omega}}\theta$, where $\mathcal{P}_{U_{h}^{\Omega}}$ is the orthogonal $L^{2}$-projection onto $U_{h}^{\Omega}$. We decompose the first term in (2.9a) into

[TABLE]

The first term will be used to obtain a norm on the dual space of $H^{1}{(\Omega)}$. The second term can be controlled via Lemma 2.1 and the $H^{1}$-stability of $\mathcal{P}_{U_{h}^{\Omega}}$ (cf. [9]). Using these considerations and Hölder’s inequality, we obtain

[TABLE]

Dividing by $\left\|{\theta}\right\|_{H^{1}{(\Omega)}}$, taking the second power on both sides, multiplying by $\tau$, and summing over all time steps provides

[TABLE]

Applying the already established regularity results and • ‣ 4 completes the proof of the left inequality in (4.5). For the case $\kappa>0$, the right inequality in (4.5) can be established using similar computations. In the case $\kappa=0$, we combine Lemma 2.1 with an inverse estimate and obtain

[TABLE]

Again, the already established regularity results and • ‣ 4 complete the proof. ∎

In order to pass to the limit ${(h,\tau)}\searrow 0$, we define time-interpolants of time-discrete functions $a^{n}$, $n=0,...,N$, and introduce some time-index-free notation as follows.

[TABLE]

We want to point out that the time derivative of $a^{\tau}$ coincides with the difference quotient, i.e.

[TABLE]

If a statement is valid for $a^{\tau}$, $a^{\tau,+}$, and $a^{\tau,-}$, we use the abbreviation $a^{\tau,(\pm)}$. With this notation, system (2.9) reads as follows.

[TABLE]

for all $\theta_{h}\in L^{2}{(0,T;U_{h}^{\Omega})}$. Similarly, we can rewrite the regularity results obtained in Corollary 4.1 and Lemma 4.2 as

[TABLE]

These regularity results can be used to identify converging subsequences.

Lemma 4.3**.**

Let the assumptions • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 4, and • ‣ 4 hold true. Furthermore, let ${(\phi_{h}^{\tau,(\pm)},\,\mu_{h}^{\tau,+},\,\mu_{\Gamma,h}^{\tau,+})}$ be a solution to (4.12). Then there exists a subsequence (again denoted by ${(\phi_{h}^{\tau,(\pm)},\,\mu_{h}^{\tau,+},\,\mu_{\Gamma,h}^{\tau,+})}$) and functions

[TABLE]

such that $\gamma{(\phi)}=\psi$ almost everywhere on $\Gamma_{T}$ and for ${(h,\tau)}\searrow 0$

[TABLE]

Proof.

The weak and weak∗ convergence expressed in (4.15a), (4.15b), (4.15g), and (4.15h) follows directly from the bounds in (4.13a) and (4.13b). The strong convergence in (4.15c) then follows from the bounds for $\phi_{h}^{\tau,(\pm)}$ in $L^{\infty}{(0,T;H^{1}{(\Omega)})}$, the bounds on $\partial_{t}\phi_{h}^{\tau}$ in $L^{2}{(0,T;{(H^{1}{(\Omega)})}^{\prime})}$, the Aubin–Lions theorem, and the fact that $\phi_{h}^{\tau,+}$, $\phi_{h}^{\tau,-}$, and $\phi_{h}^{\tau}$ converge towards the same limit function due to the bound on $\tau^{-1}\left\|{\nabla\phi_{h}^{\tau,+}-\nabla\phi_{h}^{\tau,-}}\right\|_{L^{2}{(0,T;L^{2}{(\Omega)})}}^{2}$.

Similar arguments provide (4.15d)-(4.15f) in the case $\kappa>0$. In the case $\kappa=0$, we use the uniform bound on $\left\|{\phi_{h}^{n}}\right\|_{H^{1}{(\Omega)}}$ to deduce a uniform bound for $\left\|{\gamma{(\phi_{h}^{n})}}\right\|_{H^{{(1/2)}}{(\Gamma)}}$. As $H^{1/2}{(\Gamma)}$ is compactly embedded in $L^{\tilde{s}}{(\Gamma)}$ for $\tilde{s}\in[1,\tfrac{2{(d-1)}}{d-2})$ (cf. [45]), we verify (4.15d)-(4.15f) for $\kappa=0$. It remains to show that $\psi$ can be identified with $\gamma{(\phi)}$. We choose $\boldsymbol{\theta}\in L^{2}{(0,T;{(C^{\infty}{{(\Omega)}})}^{d})}$ and compute

[TABLE]

∎

Theorem 4.4**.**

Let $d\in{\{2,3\}}$ and let the assumptions • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 2, • ‣ 4, and • ‣ 4 hold true. Then a tuple ${(\phi,\mu,\mu_{\Gamma})}$ satisfying

[TABLE]

can be obtained from discrete solutions to (2.9) by passing to the limit ${(h,\tau)}\searrow 0$. This tuple solves (1.1) in the following weak sense:

[TABLE]

Proof.

We start by passing to the limit in (4.12a). Choosing $\theta_{h}:=\mathcal{I}_{h}{\left\{\theta\right\}}$ for $\theta\in L^{2}{(0,T;C^{\infty}{(\overline{\Omega})})}$, we have $\theta_{h}\rightarrow\theta$ in $L^{2}{(0,T;H^{1}{(\Omega)})}$ (cf. [10]). We decompose the first term as

[TABLE]

This allows us to combine the results from (4.6) and (4.7) with (4.15b) and (4.15g) to derive (4.18a) for $\theta\in L^{2}{(0,T;C^{\infty}{(\overline{\Omega})})}$. Noting that $L^{2}{(0,T;C^{\infty}{(\overline{\Omega})})}$ is dense in $L^{2}{(0,T;H^{1}{(\Omega)})}$ yields the result. Similar arguments allow us to pass to the limit in (4.12b) to obtain (4.18b).

In order to pass to the limit in (4.12c), we choose $\theta_{h}\!:=\!\mathcal{I}_{h}{\left\{\theta\right\}}$ with $\theta\in L^{2}{(0,T;C^{\infty}{(\overline{\Omega})})}$ and assume that $\gamma{(\phi_{h}^{\tau,(\pm)})}\in L^{\infty}{(0,T;H^{1/2}{(\Gamma)})}$, which is the case for $\kappa>0$ and $\kappa=0$. While the convergence of the left-hand side of (4.12c) and the gradient terms is straightforward, the convergence of the terms including the derivative of the potential functions $F$ and $G$ require more finesse. We will showcase the convergence of $\int_{\Gamma_{T}}\mathcal{I}_{h}^{\Gamma}{\left\{G_{+}^{\prime}{(\phi_{h}^{\tau,+})}\theta_{h}\right\}}$. Then, the convergence of the remaining parts can be obtained in an analogous manner. According to • ‣ 2, $G_{+}^{\prime}$ can be written as the sum of a polynomial of degree three and a globally Lipschitz-continuous component ${G_{+}^{L}}^{\prime}$. We start with the decomposition

[TABLE]

The convergence of the first term on the right-hand side follows directly from (4.15f) and the strong convergence of $\theta_{h}\rightarrow\theta$. Therefore, it remains to show that the remaining terms vanish when passing to the limit. Recalling that $H^{1/2}{(\Gamma)}$ is continuously embedded in $L^{4}{(\Gamma)}$ (cf. [45]), the estimates in Lemma 2.1 and the standard inverse estimates (cf. [10]) provide

[TABLE]

Therefore, the last term in (4.20) vanishes. Furthermore, we derive the estimates

[TABLE]

and

[TABLE]

As $\phi_{h}^{\tau,+}\in L^{p}{(0,T;L^{4}{(\Gamma)})}$, we obtain the convergence of the polynomial part of $G^{\prime}_{+}$. To deal with the Lipschitz-continuous part ${G_{+}^{L}}^{\prime}$, we start with the decomposition

[TABLE]

Combining Lemma 2.1 with a standard inverse estimate, we compute

[TABLE]

Using the Lipschitz-continuity of ${G_{+}^{L}}^{\prime}$, we deduce

[TABLE]

with a constant $C$ depending on the Lipschitz-constant of ${G_{+}^{L}}^{\prime}$. Furthermore, the Lipschitz-continuity provides on every $K^{\Gamma}\in\mathcal{T}_{h}^{\Gamma}$

[TABLE]

Consequently, an inverse estimate yields

[TABLE]

which proves that $II$ will also vanish when passing to the limit. From the strong convergence (4.15f), we deduce ${G_{+}^{L}}^{\prime}{(\phi_{h}^{\tau,+})}\rightarrow{G_{+}^{L}}^{\prime}{(\gamma{(\phi)})}$ almost everywhere. Recalling ${G_{+}^{L}}^{\prime}{(\phi_{h}^{\tau,+})}\in L^{\infty}{(0,T;L^{4}{(\Gamma)})}$, we may use Vitali’s convergence theorem (see e.g. [2]) to show ${G_{+}^{L}}^{\prime}{(\phi_{h}^{\tau,+})}\rightarrow{G_{+}^{L}}^{\prime}{(\gamma{(\phi)})}$ in $L^{\infty}{(0,T;L^{\tilde{s}}{(\Gamma)})}$ for $\tilde{s}<4$. The convergence of derivatives of the concave parts of $G$ follows from the same arguments. The uniform bounds of $\phi_{h}^{\tau,(\pm)}$ in $L^{\infty}{(0,T;H^{1}{(\Omega)})}$ provide enough regularity, to adapt the previously presented arguments to three spatial dimensions, which proves the convergence of the remaining terms. As $C^{\infty}{(\overline{\Omega})}$ is dense in $X_{\kappa}$, this concludes the proof. ∎

Remark 4.5**.**

The results presented in the preceding sections carry over to the case of Allen–Cahn-type dynamic boundary conditions (cf. (1.8)), where we use

[TABLE]

instead of (2.9b). The resulting scheme reads

[TABLE]

*and is well defined, as {\left(m\left.\kern-1.2pt\mathbf{L}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}+m_{\Gamma}\left.\kern-1.2pt\mathbf{M}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}\mathbf{M}_{\Gamma}^{-1}\left.\kern-1.2pt\mathbf{M}_{\Omega}\vphantom{\big{|}}\right|_{\Gamma\times\Gamma}\right)} is obviously a symmetric, positive definite matrix.

Although, $\int_{\Gamma}\phi_{h}^{n}$ is not conserved when using Allen–Cahn-type boundary conditions, testing (4.29) by $1$ shows that $\left|{\int_{\Gamma}\phi_{h}^{n}}\right|$ is bounded. Consequently, the energy estimate still provides control over $\left\|{\phi_{h}^{n}}\right\|_{H^{1}{(\Gamma)}}$.

Testing (4.29) by $\partial_{\tau}^{-}\phi_{h}^{n}$ shows $\tau\sum_{n=1}^{N}\left\|{\partial_{\tau}^{-}\phi_{h}^{n}}\right\|_{L^{2}{(\Gamma)}}^{2}\leq C$, i.e. we obtain a slightly better regularity result for the discrete time derivative than we obtained for Cahn–Hilliard-type boundary conditions. Using the time-index-free notation introduced in (4.10), the bounds read*

[TABLE]

*with $C>0$ independent of $h$ and $\tau$. Based on these uniform bounds, we are able to identify converging subsequences and pass to the limit.

Variants of these boundary conditions are used to describe dynamic contact angles. To recover the boundary condition suggested in [47], we choose $\kappa=0$, $\delta_{\Gamma}=1$, and $G{(\phi)}=\tfrac{\gamma}{2}\sin{(\tfrac{\pi}{2}\min{\{\max{\{\phi,-1\}},1\}})}+\tfrac{|{\gamma}|}{2}$, where the parameter $\gamma$ prescribes the static contact angle via Young’s formula. As $G$ satisfies • ‣ 2 and • ‣ 2 the previous results are also valid for the boundary condition (1.9).

5. Numerical simulations

In this section, we present simulations to underline the practicality of the proposed scheme (2.23), which we implemented in the C++ framework EconDrop (cf. [32, 11, 31, 40, 39]). This framework allows for adaptivity in space and time using the ideas presented in [32], i.e. we are able to use meshes with a high resolution in the evolving interfacial area and a lower resolution in the bulk phases where $\phi\approx\pm 1$. Similarly, the time increments can be varied such that they are small, when the solution changes rapidly and larger when the solution is almost stationary. In the presented simulations, we used Newton’s method to linearize (2.23), a biconjugate gradient stabilized method to solve the arising linear system, and a preconditioned conjugate gradient method to tackle the smaller auxiliary problem.

5.1. Scenario 1

As a first test case, we consider a phase-separation scenario in $\Omega={(0,1)}^{2}$ starting from the initial condition

[TABLE]

To keep $\phi$ close to the physical meaningful interval ${[-1,+1]}$, we use the penalized double-well potential $W_{\operatorname{pen}}$ (cf. (2.8)) with penalty parameter $C_{\text{pen}}=250$ for $F$ and $G$. In this simulation, we use adaptive mesh refinement based on the criteria proposed in [32] using triangles with diameters between $\sqrt{2}\cdot 2^{-6}$ and $2^{-8}$. The dimensions of the corresponding finite element spaces are depicted in Fig. 5a. Examples of the used triangulations can be found in Fig. 2. The size of the time increment is also chosen adaptively based on the ideas presented in [32], which leads to time increments between $6.3\cdot 10^{-7}$ and $6.4\cdot 10^{-4}$. The remaining parameters are listed in Tab. 1.

The evolution of the phase-separation process is depicted in Figure 1. Thereby, the light color represents the phase $\phi=+1$ and the dark color represents $\phi=-1$ in all pictures except the first one. As the initial data only contains values in ${[-0.1,+0.1]}$, we rescaled Fig. 1a such that the light color corresponds to $\phi=+0.1$, while the dark color corresponds to $\phi=-0.1$.

The phase-separation starts with the typical wave pattern (see Figs. 1b-1d), which then transforms into two entangled spirals (cf. Figs. 1e-1f). In the course of the simulation, the spirals retreat (cf. Figs. 1f-1l) reducing the fluid-fluid contact area. The corresponding decrease in energy of depicted in Fig. 3b.

Fig. 3a shows evolution of $\int_{\Omega}\phi_{h}^{\tau}$ and $\int_{\Gamma}\phi_{h}^{\tau}$. In theory, these quantities should be conserved. In our simulation, the deviation of these quantities from their initial values is of order $10^{-4}$ and originates from the adaptivity of the used triangulation.

To compare (2.23) with the straightforward approach of solving (2.9) directly in terms of practicality, we investigate the dependence of the condition numbers on the size of the time increment. For this purpose, we computed the matrices needed in the first Newton step of both schemes for artificial time increments $\tau\in\{1\cdot 10^{-4},\,4\cdot 10^{-5},\,2\cdot 10^{-5},\,1\cdot 10^{-5},\,1\cdot 10^{-6},\,1\cdot 10^{-7}\}$ and estimated their condition number using the Matlab function condest. Averages of the condition numbers based on 51 equidistant points in time are plotted in Fig. 4. As expected, the condition numbers in the straightforward approach grow for vanishing $\tau$. Applying a Jacobi preconditioner to the system reduces the condition number drastically. However, for vanishing $\tau$ the condition number of the preconditioned system still increases rapidly (cf. blue triangles in Fig. 4). On the contrary, the condition number for (2.23), drops for decreasing $\tau$, as the matrix in (2.23) converges towards the identity.

As discussed in Rem. 2.5, solving (2.23) requires us to solve a smaller auxiliary problem repeatedly. In the discussed simulation, the average number of needed cg-iterations remained below 20 (cf. Fig. 5b).

When developing our scheme, we boiled (2.9) down to (2.23) using the fact that (2.9a) and (2.9b) have to provide identical values for the trace of $\phi_{h}^{n}$. In order to validate our scheme, we use (2.23) to compute the new phase-field values, recover $\mu_{\Gamma,h}^{n}$ via (2.18), and check whether (2.9b) still holds true. In this simulation, the $L^{2}{(\Gamma)}$-norm of the deviation averages out to $3.3\cdot 10^{-9}$.

5.2. Scenario 2

In the second scenario, we are interested in the experimental order of convergence (EOC) for our scheme. Similar to the last section, we choose $W_{\operatorname{pen}}$ (cf. (2.8)) with penalty parameter $C_{\text{pen}}=250$ for $F$ and $G$ and use $\Omega:={(0,1)}^{2}$. The remaining parameters are collected in Tab. 2. We consider the initial configuration

[TABLE]

As this configuration is unstable, the two phases will separate and form areas where $\phi$ is close to $\pm 1$. For $h=\sqrt{2}\cdot 2^{-7}$ and $\tau=2\cdot 10^{-5}$, the evolution of $\phi$ in the time interval ${[0,1]}$ is depicted in Fig. 6. Thereby, the light color represents the phase $\phi=+1$ and the dark color represents the phase $\phi=-1$. In order to visualize the initial condition, we rescaled Fig. 6a so that the dark color corresponds to $\phi=0$ and the light color to $\phi=0.1$. Again, the separation process reduces the energy of the system (cf. Fig. 7a).

To compute an experimental order of convergence, we repeat the above simulation with $\tau=k\cdot 10^{-5}$, $k\in{\{1,2,4\}}$, and $h=\sqrt{2}\cdot 2^{-l}$, $l\in{\{6,7,8\}}$. These spatial discretizations correspond to $\dim U_{h}^{\Omega}\in{\{4225,\,16641,\,66049\}}$ and $\dim U_{h}^{\Gamma}\in{\{256,\,512,\,1024\}}$. Fixing $\tau=2\cdot 10^{-5}$, we use the solution obtained for $h=\sqrt{2}\cdot 2^{-8}$ as reference solution $\overline{\phi}$ and define

[TABLE]

where the integration in time is approximated by a trapezoidal rule with step size $2\cdot 10^{-4}$. The experimental order of convergence w.r.t. $h$ is then defined as

[TABLE]

Analogously, we define $\text{err}_{h}^{\Gamma}$ and $EOC_{h}^{\Gamma}$ using the $L^{2}{(0,1;L^{2}{(\Gamma)})}$-norm. As shown in Tab. 3, we obtain order 2.3 for the convergence of $\phi_{h}^{\tau}$ w.r.t. $h$ on $\Omega$ and $1.1$ for the convergence of $\gamma{(\phi_{h}^{\tau})}$ on $\Gamma$.

In a similar manner, we fix $h=\sqrt{2}\cdot 2^{-7}$, use the solution corresponding to $\tau=1\cdot 10^{-5}$ as reference solution, and define $\text{err}_{\tau}^{\Omega}$, $\text{err}_{\tau}^{\Gamma}$, $\text{EOC}_{\tau}^{\Omega}$, and $\text{EOC}_{\tau}^{\Gamma}$ accordingly. As it can be seen in Tab. 4, the computed convergence order w.r.t. the time increment is 1.6 for $\phi_{h}^{\tau}$ and its trace $\gamma{(\phi_{h}^{\tau})}$.

Again, we conclude this section by evaluating the reliability and the efficiency of our scheme based on the conservation of mass, validity of (2.9b), the average condition number, and the average number of cg-iterations needed to solve the auxiliary problem. As illustrated in Fig. 7b, our scheme conserves $\int_{\Omega}\phi_{h}^{\tau}$ and $\int_{\Gamma}\phi_{h}^{\tau}$ perfectly for all considered values of $h$ and $\tau$, if the triangulation is fixed. As explained in the last section, we may use (2.18) to recover $\mu_{\Gamma,h}^{n}$ and compute the trace of $\phi_{h}^{n}$ using (2.9b). For $h=\sqrt{2}\cdot 2^{-7}$ and $\tau=2\cdot 10^{-5}$ the $L^{2}{(\Gamma)}$-norm of the deviation averages to $1.6\cdot 10^{-9}$.

As it was done for the last scenario, we computed the matrices used in the first Newton iteration in the different schemes for various artificial time increments and estimated their condition number using the Matlab function condest. The average condition numbers based on 51 equidistant points in time are plotted in Fig. 8. As before, the matrices arising when solving (2.9) directly are ill-conditioned for small time increments. A Jacobi preconditioner is able to reduce the condition number drastically, but can not overcome the structural problem of this approach. Fig. 7c shows the average number of cg-iterations needed to solve the auxiliary problem for $\tau=2\cdot 10^{-5}$ and different values of $h$. Again, the auxiliary problem can be solved with very few iterations.

Appendix A Appendix

For the reader’s convenience, we provide the generalized Poincaré inequality which can be found in [2].

Lemma A.1**.**

Let $\Omega\subset\mathds{R}^{d}$ be open, bounded and connected with Lipschitz boundary $\partial\Omega$. Moreover, let $1<p<\infty$ and let $\mathcal{M}\subset W^{1,p}{(\Omega)}$ be nonempty, closed and convex. Then the following items are equivalent for every $u_{0}\in\mathcal{M}$:

1. (1)

There exists a constant $C_{0}<\infty$ such that for all $\xi\in\mathds{R}$

[TABLE] 2. (2)

There exists a constant $C<\infty$ with

[TABLE]