A decoupled, stable, and linear FEM for a phase-field model of variable density two-phase incompressible surface flow

Due to many applications in the physical and biological sciences, binary (two-phase) flows have attracted increasing interest in the past decades. The phase-field method has emerged as a powerful and effective technique for modeling multi-phase problems in general. This method describes the physical system using a set of field variables that are continuous across the interfacial regions separating the phases. The main reasons for the success of the phase-field methodology are the following: it replaces sharp interfaces by thin transition regions called diffuse interfaces, making front-tracking unnecessary, and it is based on rigorous mathematics. The classical diffuse-interface description of binary incompressible fluid with density-matched components is the so-called Model H [22], which is a Navier–Stokes–Cahn–Hilliard (NSCH) system. Several generalization to the case of variable density components have been presented [2, 4, 6, 7, 8, 14, 15, 29, 42]. Some of these models relax the incompressibility constraint to quasi-compressibility and not all models satisfy Galilean invariance, local mass conservation or thermodynamic consistency. In this paper, we focus on the generalization of Model H first presented in [2]. This model is thermodynamically consistent when the density of the mixture depends linearly on the phase-field variable (concentration or volume/surface fraction). We present an extension of the model in [2] that is thermodynamically consistent for a general monotone relation of density and phase-field variable.

Driven by applications in the analysis of lateral organization of plasma membranes (see [45, 46] for more context), we adopt this new model to simulate two-phase flow dynamics on arbitrary-shaped closed smooth surfaces. While there exist many computational studies of multi-component fluid flows in planar and volumetric domains (see, e.g., [19, 28, 31, 44, 47] for some recent works), the number of papers where NSCH systems are treated on manifolds is limited. This can be explained by the fact that solving equations numerically on general surfaces poses additional difficulties that are related to the discretization of tangential differential operators and the approximate recovery of (possibly) complex shapes. In [30], the authors resort to a stream function formulation, decouple the surface Navier–Stokes problem from the surface Cahn-Hilliad problem, and use a parametric finite element approach. Yang and Kim [43] experimented with a finite difference method on staggered marker-and-cell meshes for a NSCH system embedded in the mesh.

In the present paper, we study a geometrically unfitted finite element method for the simulation of two-phase incompressible flow on surfaces. Our approach builds on earlier work on a unfitted FEM for elliptic PDEs posed on surfaces [36] called TraceFEM. Unlike some other geometrically unfitted methods for surface PDEs, TraceFEM employs sharp surface representation. The surface can be defined implicitly and no knowledge of the surface parametrization is required. This method allows to solve for a scalar quantity or a vector field on a surface, for which a parametrization or triangulation is not required. TaceFEM has been extended to the surface Stokes problem in [32, 33] and the surface Cahn–Hilliard problem in [45]. One of the main advantages of TraceFEM is that it works well for PDEs posed on evolving surfaces, including cases with topological changes [27]. Although we will not benefit from this advantage in this paper (as the surfaces are fixed), it will become handy in a future modeling of membrane deformation and fusion for multicomponent lipid bilayers.

Although technical, the numerical analysis of diffuse-interface models for two-phase fluids with matching densities can be largely built on well-established analyses for the incompressible Navier–Stokes equations and Cahn–Hilliad equations alone. See, for example, [12] for the convergence study of a coupled implicit FE method, [25] for the analysis of a decoupled FE scheme, and [21] for the analysis of a linear and second-order in time FE method. The design of energy stable, efficient, and consistent discretizations of diffuse-interface models for two-phase fluids with non-matching densities turns out to be more challenging.

Following the ideas from [18] for the incompressible Navier–Stokes equations with variable density, in [40] Shen and Yang suggested a stable time-stepping method for a NSCH system with the momentum equation modified based on an inconsistent mass conservation equation. In [41], the same authors introduced a time-stepping scheme for the thermodynamically consistent model from [2]. This scheme (which is first order and linear, and it decouples fluid and phase equations) introduces some extra terms in the momentum equation and modifies the advection flux for the order parameter to achieve energy stability. In the same paper, the authors present the stability analysis for a semi-discrete case. In [48], the approach from [41] was extended with the application of a scalar auxiliary variable method [38] and the stability was studied for a finite difference scheme on uniform staggered grids. A fully discrete FE method for the model from [2] was analyzed in [17], subject to the coupled treatment of the phase and fluid equations and a modification of the momentum equation. In [13], this analysis is extended to a posteriori estimates and adaptivity. The question of building provably stable and computationally efficient methods for thermodynamically consistent discretization models of two-phase fluids with variable densities remains largely open.

The present study contributes to the field with a first order in time, linear, and decoupled FE scheme for our generalization of the model from [2]. One advantage of our scheme is that it neither modifies the momentum equation nor alters the advection velocity. We prove that our scheme is stable under a mild time-step restriction. Moreover, the analysis tracks the dependence of the estimates on $\epsilon$, the critical model parameter that defines the width of the transition region between phases. The fact that the model is posed on arbitrary-shaped closed smooth surfaces and the use of an unfitted finite element method for spatial discretization add some further technical difficulties to our analysis. However, apart from these extra technical details, the general line of arguments can be followed for the simplified to Euclidian setting or extended to other spatial discretization techniques.

Numerical analysis is even less developed for quasi-incompressible diffuse-interface models due to their additional complexity, e.g. a more subtle nonlinearity. Here, we would like to mention the analysis of a semi-discrete method (coupled, but resulting in linear systems at each time step) in [42].

The rest of the paper is organized as follows. In Sec. 2, we state our generalization of the model from [2] and its variational formulations. A decoupled linear FEM based on the application of TraceFEM to our model is presented in Sec. 3, while the analysis is reported in Sec. 4. In Sec. 5, we report several numerical results, including a convergence test and the simulation of the Kevin–Helmholtz and Rayleigh–Taylor instabilities. Sec. 6 collects a few concluding remarks.

On surface $\Gamma$ we consider a heterogeneous mixture of two species with surface fractions $c_{i}=S_{i}/S$, $i=1,2$, where $S_{i}$ are the surface area occupied by the components and $S$ is the surface area of $\Gamma$. Since $S=S_{1}+S_{2}$, we have $c_{1}+c_{2}=1$. Let $c_{1}$ be the representative surface fraction, i.e $c=c_{1}$. Moreover, let $m_{i}$ be the mass of component $i$ and $m$ is the total mass. Notice that density of the mixture can be expressed as $\rho=\frac{m}{S}=\frac{m_{1}}{S_{1}}\frac{S_{1}}{S}+\frac{m_{2}}{S_{2}}\frac{S_{2}}{S}$. Thus,

$\displaystyle\rho=\rho(c)=\rho_{1}c+\rho_{2}(1-c).$

(2.1)

Similarly, for the dynamic viscosity of the mixture we can write

$\displaystyle\eta=\eta(c)=\eta_{1}c+\eta_{2}(1-c),~{}$

(2.2)

where $\eta_{1}$ and $\eta_{2}$ are the dynamic viscosities of the two species.

In order to state the model posed on surface $\Gamma$, we need some preliminary notation. Let $\mathbf{P}=\mathbf{P}(\mathbf{x})\coloneqq\mathbf{I}-\mathbf{n}(\mathbf{x})\mathbf{n}(\mathbf{x})^{T}$ for $\mathbf{x}\in\Gamma$ be the orthogonal projection onto the tangent plane. For a scalar function $p:\,\Gamma\to\mathbb{R}$ or a vector function $\mathbf{u}:\,\Gamma\to\mathbb{R}^{3}$ we define $p^{e}\,:\,\mathcal{O}(\Gamma)\to\mathbb{R}$, $\mathbf{u}^{e}\,:\,\mathcal{O}(\Gamma)\to\mathbb{R}^{3}$, extensions of $p$ and $\mathbf{u}$ from $\Gamma$ to its neighborhood $\mathcal{O}(\Gamma)$. The surface gradient and covariant derivatives on $\Gamma$ are then defined as $\nabla_{\Gamma}p=\mathbf{P}\nabla p^{e}$ and $\nabla_{\Gamma}\mathbf{u}\coloneqq\mathbf{P}\nabla\mathbf{u}^{e}\mathbf{P}$. These definitions are independent of a particular smooth extension of $p$ and $\mathbf{u}$ off $\Gamma$. On $\Gamma$ we consider the surface rate-of-strain tensor [20] given by

$$E_{s}(\mathbf{u})\coloneqq\frac{1}{2}(\nabla_{\Gamma}\mathbf{u}+(\nabla_{\Gamma}\mathbf{u})^{T}).$$

(2.3)

The surface divergence operators for a vector $\mathbf{g}:\Gamma\to\mathbb{R}^{3}$ and a tensor $\mathbf{A}:\Gamma\to\mathbb{R}^{3\times 3}$ are defined as:

$${\mathop{\,\rm div}}_{\Gamma}\mathbf{g}\coloneqq{\rm tr}(\nabla_{\Gamma}\mathbf{g}),\qquad{\mathop{\,\rm div}}_{\Gamma}\mathbf{A}\coloneqq\left({\mathop{\,\rm div}}_{\Gamma}(\mathbf{e}_{1}^{T}\mathbf{A}),\,{\mathop{\,\rm div}}_{\Gamma}(\mathbf{e}_{2}^{T}\mathbf{A}),\,{\mathop{\,\rm div}}_{\Gamma}(\mathbf{e}_{3}^{T}\mathbf{A})\right)^{T},$$

with $\mathbf{e}_{i}$ the $i$th standard basis vector in $\mathbb{R}^{3}$.

The classical phase-field model to describe the flow of two immiscible, incompressible, and Newtonian fluids is the so-called Model H [22]. One of the fundamental assumptions for Model H is that the densities of both components are matching. Several extensions have been proposed to account for the case of non-matching densities; see Sec. 1. Here, we restrict our attention to a thermodynamically consistent generalization of Model H first presented in [2]. For surface based quantities the model reads:

	$\displaystyle\rho\partial_{t}\mathbf{u}+\rho(\nabla_{\Gamma}\mathbf{u})\mathbf{u}-\mathbf{P}{\mathop{\,\rm div}}_{\Gamma}(2\eta E_{s}(\mathbf{u}))+\nabla_{\Gamma}p$	$\displaystyle=-\sigma_{\gamma}c\nabla_{\Gamma}\mu+{M\frac{d\rho}{dc}(\nabla_{\Gamma}\mathbf{u})\nabla_{\Gamma}\mu}+\mathbf{f},$		(2.4)
	$\displaystyle{\mathop{\,\rm div}}_{\Gamma}\mathbf{u}$	$\displaystyle=0,$		(2.5)
	$\displaystyle\partial_{t}c+{\mathop{\,\rm div}}_{\Gamma}(c\mathbf{u})-{\mathop{\,\rm div}}_{\Gamma}\left(M\nabla_{\Gamma}\mu\right)$	$\displaystyle=0,$		(2.6)
	$\displaystyle\mu$	$\displaystyle=\frac{1}{\epsilon}f_{0}^{\prime}-\epsilon\Delta_{\Gamma}c.$		(2.7)

Here, $\mathbf{u}$ is the surface averaged tangential velocity $\mathbf{u}=c\mathbf{u}_{1}+(1-c)\mathbf{u}_{2}$, density is given by (2.1) and viscosity by (2.2), $\sigma_{\gamma}$ is line tension and $\mu$ is the chemical potential defined in (2.7). Force vector $\mathbf{f}$, with $\mathbf{f}\cdot\mathbf{n}=0$, is given. In (2.6) $M$ is the mobility coefficient, while in (2.7) $f_{0}(c)$ is the specific free energy of a homogeneous phase and $\epsilon$ is a parameter related to the thickness of the interface between the two phases. We recall that in order to have phase separation, $f_{0}$ must be a non-convex function of $c$. For a comprehensive discussion of the related Navier–Stokes and Cahn–Hilliard equations on surfaces we refer, for example, to [23, 10, 9].

Problem (2.4)–(2.7) has only one additional term with respect to Model H: the middle term at the right-hand side in eq. (2.4). Notice that this term vanishes in the case of matching densities since $\frac{d\rho}{dc}=(\rho_{1}-\rho_{2})$. However, the term is crucial for thermodynamic consistency when the densities do not match and it can be interpreted as additional momentum flux due to diffusion of the components driven by the gradient of the chemical potential.

In practice, we are interested in more general relations than (2.1) between $\rho$ and $c$, since depending on the choice of $f_{0}(c)$ (and because of numerical errors while computing) the order parameter may not be constrained in $[0,1]$ and so $\rho$ and $\eta$ based on (2.1) may take physically meaningless (even negative) values. Since we do not see how to show the thermodynamic consistency of (2.4)–(2.7) for non-linear $\rho(c)$, we propose a further modification. Without the loss of generality let $\rho_{1}\geq\rho_{2}$. For a general dependence of $\rho$ on $c$, it is reasonable to assume that $\rho$ is a smooth monotonic function of $c$, i.e. $\frac{d\rho}{dc}\geq 0$ (for $\rho_{1}\geq\rho_{2}$), and so we can set

$$\frac{d\rho}{dc}=\theta^{2}.$$

(2.8)

Then, we replace (2.4) with

$$\rho\partial_{t}\mathbf{u}+\rho(\nabla_{\Gamma}\mathbf{u})\mathbf{u}-\mathbf{P}{\mathop{\,\rm div}}_{\Gamma}(2\eta E_{s}(\mathbf{u}))+\nabla_{\Gamma}p=-\sigma_{\gamma}c\nabla_{\Gamma}\mu+{M\theta(\nabla_{\Gamma}(\theta\mathbf{u})\,)\nabla_{\Gamma}\mu}+\mathbf{f}.$$

(2.9)

The updated model (2.9),(2.5)–(2.7) obviously coincides with (2.4)–(2.7) for $\rho(c)$ from (2.1), but exhibits thermodynamic consistency for a general monotone $\rho$–$c$ relation as we show below. The consistency is preserved if $M$ is a non-negative function of $c$ rather than a constant coefficient. Thermodynamic consistent extensions of (2.4)–(2.7) for a generic smooth $\rho(c)$ (no monotonicity assumption) were also considered in [1, 3] for the purpose of well-posedness analysis. Those extensions introduce yet more term(s) in the momentum equation, so for computational needs and numerical analysis purpose we opt for (2.9).

From now on, we will focus on problem (2.9), (2.5)–(2.7). This is a Navier–Stokes–Cahn–Hilliard (NSCH) system that needs to be supplemented with the definitions of mobility $M$ and free energy per unit area $f_{0}$. Following many of the existing analytic studies, as well as numerical studies, we assume mobility to be constant (i.e., independent of $c$). As for $f_{0}$, a classical choice is the Ginzburg–Landau double-well potential

$\displaystyle f_{0}(c)=\frac{\xi}{4}c^{2}(1-c)^{2},$

(2.10)

where $\xi$ defines the barrier height, i.e. the local maximum at $c=1/2$ [11]. We set $\xi=1$ for the rest of the paper.

For the numerical method, we need a weak (integral) formulation. We define the spaces

$$\mathbf{V}_{T}\coloneqq\{\,\mathbf{u}\in H^{1}(\Gamma)^{3}~{}|~{}\mathbf{u}\cdot\mathbf{n}=0\,\},\quad E\coloneqq\{\,\mathbf{u}\in\mathbf{V}_{T}~{}|~{}E_{s}(\mathbf{u})=\mathbf{0}\,\}.$$

(2.11)

We define the Hilbert space $\mathbf{V}_{T}^{0}$ as an orthogonal complement of $E$ in $\mathbf{V}_{T}$ (hence $\mathbf{V}_{T}^{0}\sim\mathbf{V}_{T}/E$), and recall the surface Korn’s inequality [23]:

$$\|\mathbf{u}\|_{H^{1}(\Gamma)}\leq C_{K}\|E_{s}(\mathbf{u})\|,\quad\forall\,\mathbf{u}\in\mathbf{V}_{T}^{0}.$$

(2.12)

In (2.12) and later we use short notation $\|\cdot\|$ for the $L^{2}$-norm on $\Gamma$. Finally, we define $L_{0}^{2}(\Gamma)\coloneqq\{\,p\in L^{2}(\Gamma)~{}|~{}\int_{\Gamma}p\,ds=0\,\}$. To devise the weak formulation, one multiplies eq. (2.9) by $\mathbf{v}\in\mathbf{V}_{T}$, eq. (2.5) by $q\in L^{2}(\Gamma)$, eq. (2.6) by $v\in H^{1}(\Gamma)$, and eq. (2.7) by $g\in H^{1}(\Gamma)$ and integrates all the equations over $\Gamma$. For eq. (LABEL:gracke-1m) and (2.6), one employs the integration by parts identity, which for a closed smooth surface $\Gamma$ reads:

$$\int_{\Gamma}v\textrm{div}\ \!_{\Gamma}\mathbf{g}\,ds=-\int_{\Gamma}\mathbf{g}\cdot\nabla_{\Gamma}v\,ds+\int_{\Gamma}\kappa v\mathbf{g}\cdot\mathbf{n}\,ds,\quad\text{for}~{}\mathbf{g}\in H^{1}(\Gamma)^{3},\,v\in H^{1}(\Gamma),$$

(2.13)

here $\kappa$ is the sum of principle curvatures. Identity (2.13) is applied to the second term in (2.6) (i.e. $\mathbf{g}=c\mathbf{u}$), which leads to no contribution from the curvature term since $\mathbf{u}$ is tangential, and to the third term in (2.6) (i.e. $\mathbf{g}=M\nabla_{\Gamma}\mu$), which makes the curvature term vanish as well. For a similar reason (component-wise), the curvature term vanishes also when identity (2.13) is applied to the diffusion term in (2.9).

The weak (integral) formulation of the surface NSCH problem (2.9), (2.5)–(2.7) reads: Find $(\mathbf{u},p,c,\mu)\in\mathbf{V}_{T}\times L_{0}^{2}(\Gamma)\times H^{1}(\Gamma)\times H^{1}(\Gamma)$ such that

	$\displaystyle\int_{\Gamma}\left(\rho\partial_{t}\mathbf{u}\cdot\mathbf{v}+\rho(\nabla_{\Gamma}\mathbf{u})\mathbf{u}\cdot\mathbf{v}+2\eta E_{s}(\mathbf{u}):E_{s}(\mathbf{v})\right)\,ds-\int_{\Gamma}p\,{\mathop{\,\rm div}}_{\Gamma}\mathbf{v}\,ds=-\int_{\Gamma}\sigma_{\gamma}c\nabla_{\Gamma}\mu\cdot\mathbf{v}\,ds$
	$\displaystyle\quad\quad+\int_{\Gamma}M(\nabla_{\Gamma}(\theta\mathbf{u}))(\nabla_{\Gamma}\mu)\cdot(\theta\mathbf{v})\,ds+\int_{\Gamma}\mathbf{f}\cdot\mathbf{v}\,ds,$		(2.14)
	$\displaystyle\int_{\Gamma}q\,{\mathop{\,\rm div}}_{\Gamma}\mathbf{u}\,ds=0,$		(2.15)
	$\displaystyle\int_{\Gamma}\partial_{t}c\,v\,ds-\int_{\Gamma}c\mathbf{u}\cdot\nabla_{\Gamma}v\,ds+\int_{\Gamma}M\nabla_{\Gamma}\mu\cdot\nabla_{\Gamma}v\,ds=0,$		(2.16)
	$\displaystyle\int_{\Gamma}\mu\,g\,ds=\int_{\Gamma}\frac{1}{\epsilon}f_{0}^{\prime}(c)\,g\,ds+\int_{\Gamma}\epsilon\nabla_{\Gamma}c\cdot\nabla_{\Gamma}g\,ds,$		(2.17)

for all $(\mathbf{v},q,v,g)\in\mathbf{V}_{T}\times L^{2}(\Gamma)\times H^{1}(\Gamma)\times H^{1}(\Gamma)$.

Indeed, the first two terms in eq. (2.14) tested with $\mathbf{v}=\mathbf{u}$ can be handled as follows:

$$\begin{split}\int_{\Gamma}\left(\rho\partial_{t}\mathbf{u}\cdot\mathbf{u}+\rho(\nabla_{\Gamma}\mathbf{u})\mathbf{u}\cdot\mathbf{u}\right)ds&=\int_{\Gamma}\left(\frac{\rho}{2}\partial_{t}|\mathbf{u}|^{2}-\frac{1}{2}{\mathop{\,\rm div}}_{\Gamma}(\rho\mathbf{u})|\mathbf{u}|^{2}\right)ds\\ &=\int_{\Gamma}\left(\frac{1}{2}\partial_{t}(\rho|\mathbf{u}|^{2})-\frac{1}{2}\left(\partial_{t}\rho+{\mathop{\,\rm div}}_{\Gamma}(\rho\mathbf{u})\right)|\mathbf{u}|^{2}\right)ds.\end{split}$$

Using generic $\rho=\rho(c)$, (2.8), and (2.5) we get

$$\partial_{t}\rho+{\mathop{\,\rm div}}_{\Gamma}(\rho\mathbf{u})=\frac{d\rho}{dc}\left(\partial_{t}c+{\mathop{\,\rm div}}_{\Gamma}(c\mathbf{u})\right)=\theta^{2}\left(\partial_{t}c+{\mathop{\,\rm div}}_{\Gamma}(c\mathbf{u})\right).$$

With the help of (2.6) , this yields

	$\displaystyle\int_{\Gamma}\left(\partial_{t}\rho+{\mathop{\,\rm div}}_{\Gamma}(\rho\mathbf{u})\right)|\mathbf{u}|^{2}ds$	$\displaystyle=-\int_{\Gamma}M(\nabla_{\Gamma}\mu)\cdot(\nabla_{\Gamma}|\theta\mathbf{u}|^{2})ds$
		$\displaystyle=-2\int_{\Gamma}M(\nabla_{\Gamma}(\theta\mathbf{u}))(\nabla_{\Gamma}\mu)\cdot(\theta\mathbf{u})ds,$

which will cancel with the second term at the right-hand side in eq. (2.14) tested with $\mathbf{v}=\mathbf{u}$. Thus, from eq. (2.14) tested with $\mathbf{v}=\mathbf{u}$ we obtain the following equality:

$$\frac{1}{2}\frac{d}{dt}\int_{\Gamma}\rho|\mathbf{u}|^{2}ds+\int_{\Gamma}2\eta|E_{s}(\mathbf{u})|^{2}ds=-\int_{\Gamma}\sigma_{\gamma}c\nabla_{\Gamma}\mu\cdot\mathbf{u}\,ds+\int_{\Gamma}\mathbf{f}\cdot\mathbf{u}\,ds,$$

(2.19)

where we have also used eq. (2.15) tested with $q=p$. From eq. (2.16) tested with $v=\mu$ and multiplied by $\sigma_{\gamma}$, we get:

$\displaystyle\int_{\Gamma}\sigma_{\gamma}(\mathbf{u}\cdot\nabla_{\Gamma}c)\,\mu\,ds=-\int_{\Gamma}\sigma_{\gamma}(\mathbf{u}\cdot\nabla_{\Gamma}\mu)\,c\,ds=-\int_{\Gamma}\sigma_{\gamma}\partial_{t}c\,\mu\,ds-\int_{\Gamma}\sigma_{\gamma}M|\nabla_{\Gamma}\mu|^{2}ds.$

(2.20)

The first term on the right-hand side can be handled using (2.7) as follows:

$$\int_{\Gamma}\partial_{t}c\,\mu\,ds=\int_{\Gamma}\partial_{t}c\left(\frac{1}{\epsilon}f_{0}^{\prime}-\epsilon\Delta_{\Gamma}c\right)ds=\frac{d}{dt}\int_{\Gamma}\left(\frac{1}{\epsilon}f_{0}+\frac{\epsilon}{2}|\nabla_{\Gamma}c|^{2}\right)ds.$$

Plugging this into (2.20) and then using what we get in (2.19), we obtain the energy balance (2.18).

For the discretization of the variational problem (2.14)–(2.17) we apply an unfitted finite element method called trace finite element approach (TraceFEM). To discretize surface equations, TraceFEM relies on a tessellation of a 3D bulk computational domain $\Omega$ ($\Gamma\subset\Omega$ holds) into shape regular tetrahedra untangled to the position of $\Gamma$.

We start with required definitions to set up finite element spaces and variational form. A few auxiliary results will be proved. Then we proceed to the fully-discrete method by introducing a splitting time discretization, which decouples (2.14)–(2.17) into one linear fluid problem and one phase-field problem per every time step.

Surface $\Gamma$ is defined implicitly as the zero level set of a sufficiently smooth (at least Lipschitz continuous) function $\phi$, i.e. $\Gamma=\{\mathbf{x}\in\Omega\,:\,\phi(\mathbf{x})=0\}$, such that $|\nabla\phi|\geq c_{0}>0$ in a 3D neighborhood $U(\Gamma)$ of the surface. The vector field $\mathbf{n}=\nabla\phi/|\nabla\phi|$ is normal on $\Gamma$ and defines quasi-normal directions in $U(\Gamma)$. Let $\mathcal{T}_{h}$ be the collection of all tetrahedra such that $\overline{\Omega}=\cup_{T\in\mathcal{T}_{h}}\overline{T}$. The subset of tetrahedra that have a nonzero intersection with $\Gamma$ is denoted by $\mathcal{T}_{h}^{\Gamma}$. The grid can be refined towards $\Gamma$, however the tetrahedra from $\mathcal{T}_{h}^{\Gamma}$ form a quasi-uniform tessellation with the characteristic tetrahedra size $h$. The domain formed by all tetrahedra in $\mathcal{T}_{h}^{\Gamma}$ is denoted by $\Omega^{\Gamma}_{h}$.

On $\mathcal{T}_{h}^{\Gamma}$ we use a standard finite element space of continuous functions that are piecewise-polynomials of degree $k$. This bulk (volumetric) finite element space is denoted by $V_{h}^{k}$:

$$V_{h}^{k}=\{v\in C(\Omega^{\Gamma}_{h})\,:\,v\in P_{k}(T)~{}\text{for any}~{}T\in\mathcal{T}_{h}^{\Gamma}\},$$

In the trace finite element method formulated below, the traces of functions from $V_{h}^{k}$ on $\Gamma$ are used to approximate the surface fraction and the chemical potential. Our bulk velocity and pressure finite element spaces are either Taylor–Hood elements on $\Omega^{\Gamma}_{h}$,

$$\mathbf{V}_{h}=(V_{h}^{m+1})^{3},\quad Q_{h}=V_{h}^{m}\cap L^{2}_{0}(\Gamma),$$

(3.1)

or equal order velocity–pressure elements

$$\mathbf{V}_{h}=(V_{h}^{m})^{3},\quad Q_{h}=V_{h}^{m}\cap L^{2}_{0}(\Gamma),\quad m\geq 1.$$

(3.2)

These spaces are employed to discretize the surface Navier–Stokes system. Surface velocity and pressure will be represented by the traces of functions from $\mathbf{V}_{h}$ and $Q_{h}$ on $\Gamma$. In general, approximation orders $k$ (for the phase-field problem) and $m$ (for the fluid problem) can be chosen to be different.

In practice, $\Gamma$ has to be approximated by a (sufficiently accurate) “discrete” surface $\Gamma_{h}\approx\Gamma$ in such a way that integrals over $\Gamma_{h}$ can be computed accurately and efficiently. For first order finite elements ($m=1$, $k=1$), a straightforward polygonal approximation of $\Gamma$ ensures that the geometric approximation error is consistent with the finite element interpolation error; see, e.g., [36]. For higher order elements, numerical approximation of $\Gamma$ based on, e.g., isoparametric trace FE can be used to recover the optimal accuracy [16] in a practical situation when parametrization of $\Gamma$ is not available explicitly.

There are two well-known issues related to the fact that we are dealing with surface and unfitted finite elements:

1. 1.

   The numerical treatment of condition $\mathbf{u}\cdot\mathbf{n}$. Enforcing the condition $\mathbf{u}_{h}\cdot\mathbf{n}=0$ on $\Gamma$ for polynomial functions $\mathbf{u}_{h}\in\mathbf{V}_{h}$ is inconvenient and may lead to locking (i.e., only $\mathbf{u}_{h}=0$ satisfies it). Instead, we add a penalty term to the weak formulation to enforce the tangential constraint weakly.

2. 2.

   Possible small cuts of tetrahedra from $\mathcal{T}_{h}^{\Gamma}$ by the surface. For the standard choice of finite element basis functions, this may lead to poorly conditioned algebraic systems. We recover algebraic stability by adding certain volumetric terms to the finite element formulation.

To make the presentation of the finite element formulation more compact, we introduce the following finite element bilinear forms related to the Cahn–Hilliard part of the problem:

	$\displaystyle a_{\mu}(\mu,v)\coloneqq\int_{\Gamma}M\nabla_{\Gamma}\mu\cdot\nabla_{\Gamma}v\,ds+\tau_{\mu}\int_{\Omega^{\Gamma}_{h}}(\mathbf{n}\cdot\nabla\mu)(\mathbf{n}\cdot\nabla v)\,d\mathbf{x}$		(3.3)
	$\displaystyle a_{c}(c,g)\coloneqq\epsilon\int_{\Gamma}\nabla_{\Gamma}c\cdot\nabla_{\Gamma}g\,ds+\tau_{c}\int_{\Omega^{\Gamma}_{h}}(\mathbf{n}\cdot\nabla c)(\mathbf{n}\cdot\nabla g)\,d\mathbf{x}.$		(3.4)

Forms (3.3)–(3.4) are well defined for $\mu,v,c,g\in H^{1}(\Omega^{\Gamma}_{h})$. The second term in (3.3) and (3.4) take care of the issue of the small element cuts (the above issue 2) [16]. These terms are consistent up to geometric errors related to the approximation of $\Gamma$ by $\Gamma_{h}$ and $\mathbf{n}$ by $\mathbf{n}_{h}$ in the following sense: any smooth solution $\mu$ and $c$ can be extended off the surface along (quasi)-normal directions so that $\mathbf{n}\cdot\nabla\mu=0$ and $\mathbf{n}\cdot\nabla c=0$ in $\Omega^{\Gamma}_{h}$. We set the stabilization parameters as follows:

$$\tau_{\mu}=h,\quad\tau_{c}=\epsilon\,h^{-1}.$$

For the stability of a numerical method, it is crucial that the computed density and viscosity stay positive. The polynomial double-well potential does not enforce $c$ to stay within $[0,1]$ interval and hence $\rho$ and $\eta$ may eventually take negative values, if one adopts the linear relation between $c$ and $\rho$ in (2.1) or between $c$ and $\eta$ in (2.2). Numerical errors may be another reason for the order parameter to depart significantly from $[0,1]$. Therefore, assuming without loss of generality that $\rho_{1}\geq\rho_{2}$ and $\eta_{1}\geq\eta_{2}$, we first replace (2.1) and (2.2) with the following cut-off functions that ensure density and viscosity stay positive:

$$\rho(c)=\left\{\begin{array}[]{cc}\rho_{2}&c\leq 0\\ c\rho_{1}+(1-c)\rho_{2}&c>0\end{array}\right.\qquad\eta(c)=\left\{\begin{array}[]{cc}\eta_{2}&c\leq 0\\ c\eta_{1}+(1-c)\eta_{2}&c>0\end{array}\right.$$

(3.5)

Note that unlike some previous studies, we clip the linear dependence (2.1) and (2.2) only from below. The resulting convexity of $\rho(c)$ plays a role in the stability analysis later. Nevertheless, (3.5) is not completely satisfactory since $\theta^{2}=\frac{\delta\rho}{\delta c}$ is discontinuous, while we need $\theta$ to be from $C^{1}$. To this end, we approximate $\rho(c)$ from (2.1) by a smooth monotone convex and uniformly positive function. In our implementation we let $\theta^{2}=\frac{\rho_{1}-\rho_{2}}{2}\left(\tanh(c/\alpha)+1\right)$, with $\alpha=0.1$, and $\rho(c)=\int_{0}^{c}\theta^{2}(t)dt+\rho_{2}$.

Later we make use of the decomposition of a vector field on $\Gamma$ into its tangential and normal components: $\mathbf{u}=\overline{\mathbf{u}}+(\mathbf{u}\cdot\mathbf{n})\mathbf{n}.$ For the Navier–Stokes part, we introduce the following forms:

	$\displaystyle a(\eta;\mathbf{u},\mathbf{v})\coloneqq\int_{\Gamma}2\eta E_{s}(\overline{\mathbf{u}}):E_{s}(\overline{\mathbf{v}})\,ds+\tau\int_{\Gamma}(\mathbf{n}\cdot\mathbf{u})(\mathbf{n}\cdot\mathbf{v})\,ds+\beta_{u}\int_{\Omega^{\Gamma}_{h}}[(\mathbf{n}\cdot\nabla)\mathbf{u}]\cdot[(\mathbf{n}\cdot\nabla)\mathbf{v}]\,d\mathbf{x},$		(3.6)
	$\displaystyle c(\rho;\mathbf{w},\mathbf{u},\mathbf{v})\coloneqq\int_{\Gamma}\rho\mathbf{v}^{T}(\nabla_{\Gamma}\overline{\mathbf{u}})\mathbf{w}\,ds+\frac{1}{2}\int_{\Gamma}\widehat{\rho}({\mathop{\,\rm div}}_{\Gamma}\overline{\mathbf{w}})\overline{\mathbf{u}}\cdot\overline{\mathbf{v}}\,ds,$		(3.7)
	$\displaystyle b(\mathbf{u},q)=\int_{\Gamma}\mathbf{u}\cdot\nabla_{\Gamma}q\,ds,$		(3.8)
	$\displaystyle s(p,q)\coloneqq\beta_{p}\int_{\Omega^{\Gamma}_{h}}(\mathbf{n}\cdot\nabla p)(\mathbf{n}\cdot\nabla q)\,d\mathbf{x}.$		(3.9)

with $\widehat{\rho}=\rho-\frac{d\rho}{d\,c}c$. Forms (3.6)–(3.9) are well defined for $p,q\in H^{1}(\Omega^{\Gamma}_{h})$, $\mathbf{u},\mathbf{v},\mathbf{w}\in H^{1}(\Omega^{\Gamma}_{h})^{3}$. In (3.6), $\tau>0$ is a penalty parameter to enforce the tangential constraint (i.e., to address the above issue 1), while $\beta_{u}\geq 0$ in (3.6) and $\beta_{p}\geq 0$ in (3.9) are stabilization parameters to deal with possible small cuts. They are set according to [24]:

$$\tau=h^{-2},\quad\beta_{p}=h,\quad\beta_{u}=h^{-1}.$$

(3.10)

If one uses equal order pressure-velocity trace elements instead of Taylor–Hood elements, then form (3.9) should be replaced by

$\displaystyle s(p,q)\coloneqq\beta_{p}\int_{\Omega^{\Gamma}_{h}}\nabla p\cdot\nabla q\,d\mathbf{x}.$

The second term in (3.7) is consistent since the divergence of the true tangential velocity is zero. To avoid differentiation of the projector operator, one may use the identity $\nabla_{\Gamma}\overline{\mathbf{u}}=\nabla_{\Gamma}\mathbf{u}-(\mathbf{u}\cdot\mathbf{n})\mathbf{H}$ to implement the $a$-form and $c$-form. For the analysis, we need the identity for the form (3.7) given in the following elementary lemma.

For the analysis in this section, we assume no forcing term, i.e. $\mathbf{f}^{n+1}=\boldsymbol{0}$ for all $n$. To avoid technical complications related to handling Killing vector fields (see, e.g. discussion in [5]), we shall also assume that the surface does not support any tangential rigid motions, and so $\mathbf{V}_{T}^{0}=\mathbf{V}_{T}$.

We further use the $a\lesssim b$ notation if inequality $a\leq c\,b$ holds between quantities $a$ and $b$ with a constant $c$ independent of $h$, $\Delta t$, $\epsilon$, and the position of $\Gamma$ in the background mesh. We give similar meaning to $a\gtrsim b$, and $a\simeq b$ means that both $a\lesssim b$ and $a\gtrsim b$ hold.

The following lemma will be useful in the proof of the main result, which is reported in Theorem 4.2.