An unfitted finite element method for two-phase Stokes problems with slip between phases

M. Olshanskii¹, A. Quaini¹ and Q. Sun¹

¹Department of Mathematics, University of Houston, 3551 Cullen Blvd, Houston TX 77204
[email protected]; [email protected]; [email protected]

Abstract

We present an isoparametric unfitted finite element approach of the CutFEM or Nitsche-XFEM family for the simulation of two-phase Stokes problems with slip between phases. For the unfitted generalized Taylor–Hood finite element pair $\mathbf{P}_{k+1}-P_{k}$ , $k\geq 1$ , we show an inf-sup stability property with a stability constant that is independent of the viscosity ratio, slip coefficient, position of the interface with respect to the background mesh and, of course, mesh size. In addition, we prove stability and optimal error estimates that follow from this inf-sup property. We provide numerical results in two and three dimensions to corroborate the theoretical findings and demonstrate the robustness of our approach with respect to the contrast in viscosity, slip coefficient value, and position of the interface relative to the fixed computational mesh.

Keywords: XFEM, cutFEM, two-phase flow, Stokes problem, finite elements

1 Introduction

The finite element approximation of two-phase problems involving immiscible fluids features several challenging aspects. The first challenge is the presence of a sharp interface between the two phases, that might move and undergo topological changes. A second critical aspect is the presence of surface tension forces that create a jump in the pressure field at the interface. In addition, if one accounts for slip between phases [25], a jump in the velocity field at the interface needs to be captured as well. Finally, lack of robustness may arise when there is a high contrast in fluid densities and viscosities. Tackling all of these challenges has motivated a large body of literature.

One possible way to categorize numerical methods proposed in the literature is to distinguish between diffusive interface and sharp interface approaches. Phase field methods (e.g., [4, 28]) belong to the first category, while level set methods (e.g., [42]), and conservative level set methods (e.g., [38]) belong to the second. Diffusive interface methods introduce a smoothing region around the interface between the two phases to vary smoothly, instead of sharply, from one phase to the other and usually apply the surface tension forces over the entire smoothing region. The major limitation of diffusive interface methods lies in the need to resolve the smoothing region with an adequate number of elements, which results in high computational costs. Sharp interface methods require less elements to resolve the interface between phases. Thus, we will restrict our attention to sharp interface approaches, which can be further divided into geometrically fitted and unfitted methods.

In fitted methods, the discretization mesh is fitted to the computational interface. Perhaps, Arbitrary Lagrangian Eulerian (ALE) methods [16] are the best known fitted methods. In case of a moving interface, ALE methods deform the mesh to track the interface. While ALE methods are known to be very robust for small interface displacement, complex re-meshing procedures are needed for large deformations and topological changes. Certain variations of the method, like the extended ALE [5, 6], successfully deal with large interface displacement while keeping the same mesh connectivity. The price to pay for such improvement is a higher computational cost. Unfitted methods allow the sharp interface to cut through the elements of a fixed background grid. Their main advantage is the relative ease of handling time-dependent domains, implicitly defined interfaces, and problems with strong geometric deformations [8]. The immersed finite element method (e.g., [3]) and front-tracking methods (e.g., [43]) are examples of unfitted approaches. Applied in the finite element framework, these methods require an enrichment of the elements intersected by the interface in order to capture jumps and kinks in the solution. One complex aspect of these methods is the need for tailored stabilization. Popular unfitted methods that embed discontinuities in finite element solvers are XFEM [36] and CutFEM [12]. XFEM enriches the finite element shape functions by the Partition-of-Unity method. To learn more about XFEM applied to two-phase flow problems, we refer the reader to [14, 19, 21, 30, 40]. CutFEM is a variation of XFEM, also called Nitsche-XFEM [23]. CutFEM uses overlapping fictitious domains in combination with ghost penalty stabilization [11] to enrich and stabilize the solution. See [15, 18, 24, 27, 35, 45] for the application of CutFEM or Nitsche-XFEM to approximate two-phase flows. Finally, recently proposed unfitted methods are a hybrid high-order method [10] and an enriched finite element/level-set method [26].

In this paper, we study an isoparametric unfitted finite element approach of the CutFEM or Nitsche-XFEM family for the simulation of two-phase Stokes problems with slip between phases. All the numerical works cited above consider the homogeneous model of two-phase flow, i.e. no slip is assumed between the phases. This assumption is appropriate in three cases: one of the phases has a relatively small volume, one phase forms drops of minute size, or one phase (representing the continuous medium in which droplets are immersed) has high speed [25]. In all other cases, slip between the phases has to be accounted for. In fact, experimentally it is observed that the velocity of the two phases can be significantly different, also depending on the flow pattern (e.g., plug flow, annular flow, bubble flow, stratified flow, slug flow, churn flow) [29]. A variation of our unfitted approach has been analyzed for the homogeneous two-phase Stokes problem in [13], where robust estimates were proved for individual terms of the Cauchy stress tensor. In the present paper, the analysis is done in the energy norm, allowing a possible slip between phases. In particular, we show an inf-sup stability property of the unfitted generalized Taylor–Hood finite element pair $\mathbf{P}_{k+1}-P_{k}$ , $k\geq 1$ , with a stability constant that is independent of the viscosity ratio, slip coefficient, position of the interface with respect to the background mesh, and of course mesh size. This inf-sup property implies stability and optimal error estimates for the unfitted finite element method under consideration, which are also shown. For more details on the isoparametric unfitted finite element, we refer to [31, 32, 34].

Two-phase flow problems with high contrast for the viscosity are known to be especially challenging. While some authors test different viscosity ratios but do not comment on the effects of high contrast on the numerics [15, 26, 46], others show or prove that their method is robust for all viscosity ratios [27, 10, 30, 37, 45]. In other cases, numerical parameters, like the penalty parameters, are adjusted to take into account large differences in the viscosity [18]. Through analysis and a series of numerical tests in two and three dimensions, we demonstrate that our approach is robust not only with respect to the contrast in viscosity, but also with respect to the slip coefficient value and the position of the interface relative to the fixed computational mesh.

For all the simulations in this paper, we have used NGsolve [1, 20], a high performance multiphysics finite element software with a Python interface, and add-on library ngsxfem [2], which enables the use of unfitted finite element technologies.

The remainder of the paper is organizes as follows. In Section 2, we introduce the strong and weak formulations of the two-phase Stokes problem with slip between phases, together with the finite element discretization. We present a stability result in Sec. 3, while in Sec. 4 we prove optimal order convergence for the proposed unfitted finite element approach. Numerical results in 2 and 3 dimensions are shown in Sec. 5. Concluding remarks are provided in Sec. 6.

2 Problem definition

We consider a fixed domain $\Omega\subset\mathbb{R}^{d}$ , with $d=2,3$ , filled with two immiscible, viscous, and incompressible fluids separated by an interface $\Gamma$ . In this study, we assume $\Gamma$ does not evolve with time although our approach is designed to handle interface evolution. We assume that $\Gamma$ is closed and sufficiently smooth. Interface $\Gamma$ separates $\Omega$ into two subdomains (phases) $\Omega^{+}$ and $\Omega^{-}=\Omega\setminus\overline{\Omega^{+}}$ . We assume $\Omega^{-}$ to be completely internal, i.e. $\partial\Omega^{-}\cap\partial\Omega=\emptyset$ . See Fig. 1. Let $\mathbf{n}^{\pm}$ be the outward unit normal for $\Omega^{\pm}$ and $\mathbf{n}$ the outward pointing unit normal on $\Gamma$ . It holds that $\mathbf{n}^{-}=\mathbf{n}$ and $\mathbf{n}^{+}=-\mathbf{n}$ at $\Gamma$ .

Refer to caption — Figure 1: Illustration of a domain $\Omega$ in $\mathbb{R}^{2}$ . On part of the boundary (dashed line) a Neumann boundary condition is imposed, while on the remaining part of the boundary (solid line with three bars) a Dirichlet boundary condition is enforced.

Let $\mathbf{u}^{\pm}:\Omega^{\pm}\to\mathbb{R}^{d}$ and $p^{\pm}:\Omega^{\pm}\to\mathbb{R}$ denote the fluid velocity and pressure, respectively. We assume that the motion of the fluids occupying subdomains $\Omega^{\pm}$ can be modeled by the Stokes equations

	$\displaystyle-\nabla\cdot\boldsymbol{\sigma}^{\pm}$	$\displaystyle=\mathbf{f}^{\pm}$	$\displaystyle\text{ in }\Omega^{\pm},$		(2.1)
	$\displaystyle\nabla\cdot{}\mathbf{u}^{\pm}$	$\displaystyle=0$	$\displaystyle\text{ in }\Omega^{\pm},$		(2.2)

endowed with boundary conditions

	$\displaystyle\mathbf{u}^{+}$	$\displaystyle=\mathbf{g},$	$\displaystyle\text{ on }\partial\Omega_{D},$		(2.3)
	$\displaystyle\boldsymbol{\sigma}^{+}\mathbf{n}^{+}$	$\displaystyle=\mathbf{g}_{N}$	$\displaystyle\text{ on }\partial\Omega_{N}.$		(2.4)

Here, $\overline{\partial\Omega_{D}}\cup\overline{\partial\Omega_{N}}=\overline{\partial\Omega}$ and $\partial\Omega_{D}\cap\partial\Omega_{N}=\emptyset$ . See Fig. 1. In (2.1), $\mathbf{f}^{\pm}$ are external the body forces and $\boldsymbol{\sigma}^{\pm}$ are the Cauchy stress tensors. For Newtonian fluids, the Cauchy stress tensor has the following expression:

\displaystyle\boldsymbol{\sigma}^{\pm}=-p^{\pm}\mathbf{I}+2\mu_{\pm}\mathbf{D}(\mathbf{u}^{\pm}),\quad\mathbf{D}(\mathbf{u}^{\pm})=\frac{1}{2}(\nabla\mathbf{u}^{\pm}+(\nabla\mathbf{u}^{\pm})^{T})\text{ in }\Omega^{\pm},

where constants $\mu_{\pm}$ represent the fluid dynamic viscosities. Finally, $\mathbf{g}$ and $\mathbf{g}_{N}$ in (2.3) and (2.4) are given.

Subproblems (2.1)-(2.2) are coupled at the interface $\Gamma$ . The conservation of mass requires the balance of normal fluxes on $\Gamma$ :

\displaystyle\mathbf{u}^{+}\cdot\mathbf{n}

\displaystyle=\mathbf{u}^{-}\cdot\mathbf{n}\qquad\text{ on }\Gamma.

(2.5)

This is the first coupling condition. We are interested in modelling slip with friction between the two phases. Thus, we consider the following additional coupling conditions:

	$\displaystyle\mathbf{P}{\boldsymbol{\sigma}^{+}\mathbf{n}}$	$\displaystyle=f(\mathbf{P}\mathbf{u}^{+}-\mathbf{P}\mathbf{u}^{-})$	$\displaystyle\text{ on }\Gamma,$		(2.6)
	$\displaystyle\mathbf{P}{\boldsymbol{\sigma}^{-}\mathbf{n}}$	$\displaystyle=-f(\mathbf{P}\mathbf{u}^{-}-\mathbf{P}\mathbf{u}^{+})$	$\displaystyle\text{ on }\Gamma,$		(2.7)

where $f$ is a constant that can be seen as a slip coefficient and $\mathbf{P}=\mathbf{P}(\mathbf{x})=I-\mathbf{n}(\mathbf{x})\mathbf{n}(\mathbf{x})^{T}$ for $\mathbf{x}\in\Gamma$ is the orthogonal projection onto the tangent plane. Finally, the jump of the normal stress across $\Gamma$ is given by:

\displaystyle[\mathbf{n}^{T}\boldsymbol{\sigma}\mathbf{n}]^{-}_{+}

\displaystyle=\sigma\kappa\qquad\text{ on }\Gamma,

(2.8)

where $\sigma$ is the surface tension coefficient and $\kappa$ is the double mean curvature of the interface.

Since the boundary conditions on $\partial\Omega$ do not affect the subsequent discussion, from now on we will consider that a Dirichlet condition (2.3) is imposed on the entire boundary. This will simplify the presentation of the fully discrete problem.

2.1. Variational formulation

The purpose of this section is to derive the variational formulation of coupled problem (2.1)–(2.8). Let us introduce some standard notation. The space of functions whose square is integrable in a domain $\omega$ is denoted by $L^{2}(\omega)$ . With $L^{2}_{0}(\omega)$ , we denote the space of functions in $L^{2}(\omega)$ with zero mean value over $\omega$ . The space of functions whose distributional derivatives of order up to $m\geq 0$ (integer) belong to $L^{2}(\omega)$ is denoted by $H^{m}(\omega)$ . The space of vector-valued functions with components in $L^{2}(\omega)$ is denoted with $L^{2}(\omega)^{d}$ . $H^{1}(\textrm{div}\ \!,\omega)$ is the space of functions in $L^{2}(\omega)$ with divergence in $L^{2}(\omega)$ . Moreover, we introduce the following functional spaces:

	$\displaystyle V^{-}$	$\displaystyle=H^{1}(\Omega^{-})^{d},~{}V^{+}=\{\mathbf{u}\in H^{1}(\Omega^{+})^{d},\mathbf{u}{\big{\|}_{\partial\Omega_{D}}}=\mathbf{g}\},~{}V^{+}_{0}=\{\mathbf{u}\in H^{1}(\Omega^{+})^{d},\mathbf{u}{\big{\|}_{\partial\Omega_{D}}}=\boldsymbol{0}\},$
	$\displaystyle V^{\pm}$	$\displaystyle=\{\mathbf{u}=(\mathbf{u}^{-},\mathbf{u}^{+})\in V^{-}\times V^{+},\mathbf{u}^{-}\cdot\mathbf{n}=\mathbf{u}^{+}\cdot\mathbf{n}~{}\text{on}~{}\Gamma\},$
	$\displaystyle V^{\pm}_{0}$	$\displaystyle=\{\mathbf{u}=(\mathbf{u}^{-},\mathbf{u}^{+})\in V^{-}\times V^{+}_{0},\mathbf{u}^{-}\cdot\mathbf{n}=\mathbf{u}^{+}\cdot\mathbf{n}~{}\text{on}~{}\Gamma\},$
	$\displaystyle Q^{\pm}$	$\displaystyle=\{p=(p^{-},p^{+})\in L^{2}(\Omega^{-})\times L^{2}(\Omega^{+})\}.$

Notice that space $V^{\pm}$ can be also characterized as $(V^{-}\times V^{+})\cap H^{1}(\textrm{div}\ \!,\Omega)$ . We use $(\cdot,\cdot)_{\omega}$ and $\langle,\rangle_{\omega}$ to denote the $L^{2}$ product and the duality pairing, respectively.

The integral formulation of the problem (2.1)-(2.8) reads: Find $(\mathbf{u},p)\in V^{\pm}\times L^{2}(\Omega)/\mathbb{R}$ such that

	$\displaystyle-(p^{-},\nabla\cdot\mathbf{v}^{-})_{\Omega^{-}}-(p^{+},\nabla\cdot\mathbf{v}^{+})_{\Omega^{+}}+2(\mu_{-}\mathbf{D}(\mathbf{u}^{-}),\mathbf{D}(\mathbf{v}^{-}))_{\Omega^{-}}+2(\mu_{+}\mathbf{D}(\mathbf{u}^{+}),\mathbf{D}(\mathbf{v}^{+}))_{\Omega^{+}}$
	$\displaystyle\quad\quad+\langle f(\mathbf{P}\mathbf{u}^{-}-\mathbf{P}\mathbf{u}^{+}),\mathbf{P}\mathbf{v}^{-}\rangle_{\Gamma}+\langle f(\mathbf{P}\mathbf{u}^{+}-\mathbf{P}\mathbf{u}^{-}),\mathbf{P}\mathbf{v}^{+}\rangle_{\Gamma}$
	$\displaystyle\quad\quad=(\mathbf{f}^{-},\mathbf{v}^{-})_{\Omega^{-}}+(\mathbf{f}^{+},\mathbf{v}^{+})_{\Omega^{+}}+\langle\sigma\kappa,\mathbf{v}^{-}\cdot\mathbf{n}\rangle_{\Gamma}$		(2.9)
	$\displaystyle(\nabla\cdot\mathbf{u}^{-},q^{-})_{\Omega^{-}}+(\nabla\cdot\mathbf{u}^{+},q^{+})_{\Omega^{+}}=0$		(2.10)

for all $(\mathbf{v},q)\in V^{\pm}_{0}\times Q^{\pm}$ . The interface terms in (2.9) have been obtained using coupling conditions (2.6), (2.7), and (2.8) as follows:

	$\displaystyle-\langle\boldsymbol{\sigma}^{-}\mathbf{n},\mathbf{v}^{-}\rangle_{\Gamma}+\langle\boldsymbol{\sigma}^{+}\mathbf{n},\mathbf{v}^{+}\rangle_{\Gamma}$	$\displaystyle=-\langle\mathbf{P}\boldsymbol{\sigma}^{-}\mathbf{n},\mathbf{P}\mathbf{v}^{-}\rangle_{\Gamma}+\langle\mathbf{P}\boldsymbol{\sigma}^{+}\mathbf{n},\mathbf{P}\mathbf{v}^{+}\rangle_{\Gamma}-\langle[\mathbf{n}^{T}\boldsymbol{\sigma}\mathbf{n}]^{-}_{+},\mathbf{v}^{-}\cdot\mathbf{n}\rangle_{\Gamma}$
		$\displaystyle=\langle f(\mathbf{P}\mathbf{u}^{-}-\mathbf{P}\mathbf{u}^{+}),\mathbf{P}\mathbf{v}^{-}\rangle_{\Gamma}+\langle f(\mathbf{P}\mathbf{u}^{+}-\mathbf{P}\mathbf{u}^{-}),\mathbf{P}\mathbf{v}^{+}\rangle_{\Gamma}$
		$\displaystyle\quad-\langle\sigma\kappa,\mathbf{v}^{-}\cdot\mathbf{n}\rangle_{\Gamma}.$

Notice that problem (2.9)-(2.10) can be rewritten as: Find $(\mathbf{u},p)\in V^{\pm}\times L^{2}(\Omega)/\mathbb{R}$ such that

\left\{\begin{split}a(\mathbf{u},\mathbf{v})+b(\mathbf{v},p)&=r(\mathbf{v})\\ b(\mathbf{u},q)&=0\end{split}\right.

(2.11)

for all $(\mathbf{v},q)\in V_{0}^{\pm}\times Q^{\pm}$ , where

	$\displaystyle a(\mathbf{u},\mathbf{v})=$	$\displaystyle 2(\mu_{-}\mathbf{D}(\mathbf{u}^{-}),\mathbf{D}(\mathbf{v}^{-}))_{\Omega^{-}}+2(\mu_{+}\mathbf{D}(\mathbf{u}^{+}),\mathbf{D}(\mathbf{v}^{+}))_{\Omega^{+}}+\langle f(\mathbf{P}\mathbf{u}^{-}-\mathbf{P}\mathbf{u}^{+}),\mathbf{P}\mathbf{v}^{-}-\mathbf{P}\mathbf{v}^{+}\rangle_{\Gamma},$
	$\displaystyle b(\mathbf{v},p)=$	$\displaystyle-(p^{-},\nabla\cdot\mathbf{v}^{-})_{\Omega^{-}}-(p^{+},\nabla\cdot\mathbf{v}^{+})_{\Omega^{+}},$
	$\displaystyle r(\mathbf{v})=$	$\displaystyle(\mathbf{f}^{-},\mathbf{v}^{-})_{\Omega^{-}}+(\mathbf{f}^{+},\mathbf{v}^{+})_{\Omega^{+}}+\langle\sigma\kappa,\mathbf{v}^{-}\cdot\mathbf{n}\rangle_{\Gamma}.$

2.2. Finite element discretization

We consider a family of shape regular triangulations $\{\mathcal{T}_{h}\}_{h>0}$ of $\Omega$ . We adopt the convention that the elements $T$ and edges $e$ are open sets and use the over-line symbol to refer to their closure. Let $h_{T}$ denote the diameter of element $T\in\mathcal{T}_{h}$ and $h_{e}$ the diameter of edge $e$ . The set of elements intersecting $\Omega^{\pm}$ and the set of elements having a nonzero intersection with $\Gamma$ are

\displaystyle\mathcal{T}^{\pm}_{h}=\{T\in\mathcal{T}_{h}:T\cap\Omega^{\pm}\neq\emptyset\},\quad\mathcal{T}_{h}^{\Gamma}=\{T\in\mathcal{T}_{h}:\overline{T}\cap\Gamma\neq\emptyset\},

(2.12)

respectively. We assume $\{\mathcal{T}_{h}^{\Gamma}\}$ to be quasi-uniform. However, in practice adaptive mesh refinement is possible. The domain formed by all tetrahedra in $\mathcal{T}_{h}^{\Gamma}$ is denoted by $\Omega^{\Gamma}_{h}:=\text{int}(\cup_{T\in\mathcal{T}_{h}^{\Gamma}}\overline{T})$ . We define the $h$ -dependent domains:

\displaystyle\Omega_{h}^{\pm}=\text{int}\left(\cup_{T\in\mathcal{T}_{h}^{\pm}}\overline{T}\right)

(2.13)

and the set of faces of $\mathcal{T}_{h}^{\Gamma}$ restricted to the interior of $\Omega_{h}^{\pm}$ :

\displaystyle\mathcal{E}_{h}^{\Gamma,\pm}=\{e=\text{int}(\partial T_{1}\cap\partial T_{2}):T_{1},T_{2}\in\mathcal{T}_{h}^{\pm}~{}\text{and}~{}T_{1}\cap\Gamma\neq\emptyset~{}\text{or}~{}T_{2}\cap\Gamma\neq\emptyset\}.

(2.14)

For the space discretization of the bulk fluid problems, we restrict our attention to inf-sup stable finite element pair $\mathbf{P}_{k+1}-P_{k}$ , $k\geq 1$ , i.e. Taylor-Hood elements. Specifically, we consider the spaces of continuous finite element pressures given by:

\displaystyle Q_{h}^{-}=\{p\in C(\Omega_{h}^{-}):q|_{T}\in P_{k}(T)~{}\forall T\in\mathcal{T}^{-}_{h}\}.

(2.15)

Space $Q_{h}^{+}$ is defined analogously. Our pressure space is given by:

\displaystyle Q^{\pm}_{h}=\{p=(p^{-},p^{+})\in Q_{h}^{-}\times Q_{h}^{+}\,:\,\int_{\Omega^{-}}\mu_{-}^{-1}p^{-}+\int_{\Omega^{+}}\mu_{+}^{-1}p^{+}=0\}.

Let

\displaystyle V_{h}^{-}

\displaystyle=\{\mathbf{u}\in C(\Omega_{h}^{-})^{d}:\mathbf{u}|_{T}\in\mathbf{P}_{k+1}(T)~{}\forall T\in\mathcal{T}^{-}_{h}\}.

(2.16)

with the analogous definition for $V_{h}^{+}$ . Our velocity spaces are given by:

\displaystyle V^{\pm}_{h}=\{\mathbf{u}=(\mathbf{u}^{-},\mathbf{u}^{+})\in V_{h}^{-}\times V_{h}^{+}\}

and $V_{0,h}^{\pm}$ , a subspace of $V^{\pm}_{h}$ with vector functions $\mathbf{u}^{+}$ vanishing on $\partial\Omega$ . All above constructions and spaces readily carry over to tessellations of $\Omega$ into squares or cubes and using $\mathbf{Q}_{k+1}-Q_{k}$ elements.

Functions in $Q^{\pm}_{h}$ and $V^{\pm}_{h}$ and their derivatives are multivalued in $\Omega^{\Gamma}_{h}$ , the overlap of $\Omega_{h}^{-}$ and $\Omega_{h}^{+}$ . The jump of a multivalued function over the interface is defined as the difference of components coming from $\Omega_{h}^{-}$ and $\Omega_{h}^{+}$ , i.e. $[\mathbf{u}]=\mathbf{u}^{-}-\mathbf{u}^{+}$ on $\Gamma$ . Note that this is the jump that we have previously denoted with $[\cdot]^{-}_{+}$ . We are now using $[\cdot]$ to simplify the notation. Moreover, we define the following averages:

	$\displaystyle\{\mathbf{u}\}=\alpha\mathbf{u}^{+}+\beta\mathbf{u}^{-},$		(2.17)
	$\displaystyle\langle\mathbf{u}\rangle=\beta\mathbf{u}^{+}+\alpha\mathbf{u}^{-},$		(2.18)

where $\alpha$ and $\beta$ are weights to be chosen such that $\alpha+\beta=1$ , $0\leq\alpha,\beta\leq 1$ . For example, in [15] the setting $\alpha=\mu_{-}/(\mu_{+}+\mu_{-})$ and $\beta=\mu_{+}/(\mu_{+}+\mu_{-})$ is suggested. In [13], the authors choose $\alpha=0$ , $\beta=1$ if $\mu_{-}\leq\mu_{+}$ and $\alpha=1$ , $\beta=0$ otherwise. Below, in (2.22) and (2.25) we will use relationship:

\displaystyle[ab]=[b]\{a\}+\left\langle b\right\rangle[a].

(2.19)

A discrete variational analogue of problem (2.11) reads: Find $\{\mathbf{u}_{h},p_{h}\}\in V^{\pm}_{h}\times Q^{\pm}_{h}$ such that

\left\{\begin{split}a_{h}(\mathbf{u}_{h},\mathbf{v}_{h})+b_{h}(\mathbf{v}_{h},p_{h})&=r_{h}(\mathbf{v}_{h})\\ b_{h}(\mathbf{u}_{h},q_{h})-b_{p}(p_{h},q_{h})&=0\end{split}\right.

(2.20)

for all $(\mathbf{v}_{h},q_{h})\in V_{0,h}^{\pm}\times Q_{h}^{\pm}$ . We define all the bilinear forms in (2.20) for all $\mathbf{u}_{h}\in V_{h}^{\pm}$ , $\mathbf{v}_{h}\in V^{\pm}_{0,h}$ , $p\in Q^{\pm}$ . Let us start with form $a_{h}(\cdot,\cdot)$ :

\displaystyle a_{h}(\mathbf{u}_{h},\mathbf{v}_{h})=

\displaystyle a_{i}(\mathbf{u}_{h},\mathbf{v}_{h})+a_{n}(\mathbf{u}_{h},\mathbf{v}_{h})+a_{p}(\mathbf{u}_{h},\mathbf{v}_{h}),

(2.21)

where we group together the terms that arise from the integration by parts of the divergence of the stress tensors:

	$\displaystyle a_{i}(\mathbf{u}_{h},\mathbf{v}_{h})=$	$\displaystyle 2(\mu_{-}\mathbf{D}(\mathbf{u}_{h}^{-}),\mathbf{D}(\mathbf{v}_{h}^{-}))_{\Omega^{-}}+2(\mu_{+}\mathbf{D}(\mathbf{u}_{h}^{+}),\mathbf{D}(\mathbf{v}_{h}^{+}))_{\Omega^{+}}+\langle f[\mathbf{P}\mathbf{u}_{h}],[\mathbf{P}\mathbf{v}_{h}]\rangle_{\Gamma}$
		$\displaystyle-2\langle\{\mu\mathbf{n}^{T}\mathbf{D}(\mathbf{u}_{h})\mathbf{n}\},[\mathbf{v}_{h}\cdot\mathbf{n}]\rangle_{\Gamma},$		(2.22)

and the terms that enforce condition (2.5) weakly using Nitsche’s method

\displaystyle a_{n}(\mathbf{u}_{h},\mathbf{v}_{h})=\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\frac{\gamma}{h_{T}}\{\mu\}\langle[\mathbf{u}_{h}\cdot\mathbf{n}],[\mathbf{v}_{h}\cdot\mathbf{n}]\rangle_{\Gamma\cap T}-2\langle\{\mu\mathbf{n}^{T}\mathbf{D}(\mathbf{v}_{h})\mathbf{n}\},[\mathbf{u}_{h}\cdot\mathbf{n}]\rangle_{\Gamma}.

(2.23)

We recall that $h_{T}$ is the diameter of element $T\in\mathcal{T}_{h}$ . To define the penalty terms $a_{p}(\mathbf{u}_{h},\mathbf{v}_{h})$ we need $\omega_{e}$ , the facet patch for $e\in\mathcal{E}_{h}^{\Gamma,\pm}$ consisting of all $T\in\mathcal{T}_{h}$ sharing $e$ . Then, we set

	$\displaystyle a_{p}(\mathbf{u}_{h},\mathbf{v}_{h})$	$\displaystyle={\mu_{-}}\mathbf{J}_{h}^{-}(\mathbf{u}_{h},\mathbf{v}_{h})+{\mu_{+}}\mathbf{J}_{h}^{+}(\mathbf{u}_{h},\mathbf{v}_{h}),$
	$\displaystyle\mathbf{J}_{h}^{\pm}(\mathbf{u}_{h},\mathbf{v}_{h})$	$\displaystyle=\gamma^{\pm}_{\mathbf{u}}\sum_{e\in\mathcal{E}_{h}^{\Gamma,\pm}}\frac{1}{h^{2}_{e}}\int_{\omega_{e}}(\mathbf{u}_{1}^{e}-\mathbf{u}_{2}^{e})\cdot(\mathbf{v}_{1}^{e}-\mathbf{v}_{2}^{e})dx,$		(2.24)

where $\mathbf{u}_{1}^{e}$ is the componentwise canonical extension of a polynomial vector function $\mathbf{u}^{\pm}_{h}$ from $T_{1}$ to $\mathbb{R}^{d}$ , while $\mathbf{u}_{2}^{e}$ is the canonical extension of $\mathbf{u}^{\pm}_{h}$ from $T_{2}$ to $\mathbb{R}^{d}$ (and similarly for $\mathbf{v}_{1}$ , $\mathbf{v}_{2}$ ). We recall that $h_{e}$ is the diameter of facet $e\in\mathcal{E}_{h}^{\Gamma,\pm}$ . This version of the ghost penalty stabilization has been proposed in [39]. In [33], it was shown to be essentially equivalent to other popular ghost penalty stabilizations such as local projection stabilization [11] and normal derivative jump stabilization [12]. In the context of the Stokes problem, this stabilization was recently used in [44]. For the analysis in Sec. 3 and 4, we also define $\mathbf{J}_{h}^{\pm}(\mathbf{u},\mathbf{v})$ for arbitrary smooth functions $\mathbf{u},\mathbf{v}$ in $\Omega_{h}^{\pm}$ . In this case, we set $\mathbf{u}_{1}=\left(\Pi_{T_{1}}\mathbf{u}|_{T_{1}}\right)^{e}$ , $\mathbf{u}_{2}=\left(\Pi_{T_{2}}\mathbf{u}|_{T_{2}}\right)^{e}$ , where $\Pi_{T_{i}}$ is the $L^{2}(T_{i})$ -orthogonal projection into the space of degree $k+1$ polynomial vector functions on $T_{i}$ .

The remaining terms coming from the integration by parts of the divergence of the stress tensors are contained in

\displaystyle b_{h}(\mathbf{v}_{h},p_{h})=

\displaystyle-(p_{h}^{-},\nabla\cdot\mathbf{v}_{h}^{-})_{\Omega^{-}}-(p_{h}^{+},\nabla\cdot\mathbf{v}_{h}^{+})_{\Omega^{+}}+\left\langle\{p_{h}\},[\mathbf{v}_{h}\cdot\mathbf{n}]\right\rangle_{\Gamma},

(2.25)

and the penalty terms are grouped together in

	$\displaystyle b_{p}(p_{h},q_{h})$	$\displaystyle=\mu_{-}^{-1}J_{h}^{-}(p_{h},q_{h})+\mu_{+}^{-1}J_{h}^{+}(p_{h},q_{h}),$
	$\displaystyle J_{h}^{\pm}(p_{h},q_{h})$	$\displaystyle={\gamma_{p}^{\pm}}\sum_{e\in\mathcal{E}_{h}^{\Gamma,\pm}}\int_{\omega_{e}}(p_{1}^{e}-p_{2}^{e})(q_{1}^{e}-q_{2}^{e})dx,$		(2.26)

where $p_{1}^{e},p_{2}^{e},q_{1}^{e},q_{2}^{e}$ are canonical polynomial extensions as defined above.

Finally,

\displaystyle r_{h}(\mathbf{v}_{h})=

\displaystyle(\mathbf{f}_{h}^{-},\mathbf{v}_{h}^{-})_{\Omega^{-}}+(\mathbf{f}_{h}^{+},\mathbf{v}_{h}^{+})_{\Omega^{+}}+\langle\sigma\kappa,\left\langle\mathbf{v}_{h}\cdot\mathbf{n}\right\rangle\rangle_{\Gamma}.

We recall that some of the interface terms in $a_{i}(\cdot,\cdot)$ and $b_{h}(\cdot,\cdot)$ have been obtained using relationship (2.19) and interface conditions.

Parameters $\gamma^{\pm}_{\mathbf{u}}$ , $\gamma_{p}^{\pm}$ and $\gamma$ are all assumed to be independent of $\mu_{\pm}$ , $h$ , and the position of $\Gamma$ against the underlying mesh. Parameter $\gamma$ in (2.23) needs to be large enough to provide the bilinear form $a_{h}(\cdot,\cdot)$ with coercivity. Parameters $\gamma^{\pm}_{\mathbf{u}}$ , $\gamma_{p}^{\pm}$ can be tuned to improve the numerical performance of the method.

2.2.1. Numerical integration

It is not feasible to compute integrals entering the definition of the bilinear forms over cut elements and over $\Gamma$ for an arbitrary smooth $\Gamma$ . We face the same problem if $\Gamma$ is given implicitly as a zero level of a piecewise polynomial function for polynomial degree greater than one. Piecewise linear approximation of $\Gamma$ on the given mesh and polygonal approximation of subdomains lead to second order geometric consistency error, which is suboptimal for Taylor–Hood elements. To ensure a geometric error of the same order or higher than the finite element (FE) approximation error, we define numerical quadrature rules on the given mesh using the isoparametric approach proposed in [31].

In the isoparametric approach, one considers a smooth function $\phi$ such that $\pm\phi>0$ in $\Omega^{\pm}$ and $|\nabla\phi|>0$ in a sufficiently wide strip around $\Gamma$ . Next, one defines polygonal auxiliary domains $\Omega^{\pm}_{1}$ given by $\Omega^{\pm}_{1}:=\{\mathbf{x}\in\mathbb{R}^{d}\,:\,\pm I_{h}^{1}(\phi)>0\},$ where $I_{h}^{1}$ is the continuous piecewise linear interpolation of $\phi$ on $\mathcal{T}_{h}$ . Interface $\Gamma_{1}$ between $\Omega^{+}_{1}$ and $\Omega^{-}_{1}$ is then $\Gamma_{1}:=\{\mathbf{x}\in\mathbb{R}^{d}\,:\,I_{h}^{1}(\phi)=0\}.$ On $\Omega^{\pm}_{1}$ and $\Gamma_{1}$ standard quadrature rules can be applied elementwise. Since using $\Omega^{\pm}_{1}$ , $\Gamma_{1}$ alone limits the accuracy to second order, one further constructs a transformation of the mesh in $\mathcal{T}^{\Gamma}_{h}$ with the help of an explicit mapping $\Psi_{h}$ parameterized by a finite element function. The mapping $\Psi_{h}$ is such that $\Gamma_{1}$ is mapped approximately onto $\Gamma$ ; see [31] for how $\Psi_{h}$ is constructed. Then, $\widetilde{\Omega}^{\pm}=\Psi_{h}(\Omega^{\pm}_{1})$ , $\widetilde{\Gamma}=\Psi_{h}(\Gamma_{1})$ are high order accurate approximations to the phases and interface which have an explicit representation so that the integration over $\widetilde{\Omega}^{\pm}$ and $\widetilde{\Gamma}$ can be done exactly. The finite element spaces have to be adapted correspondingly, using the explicit pullback mapping: $\mathbf{v}_{h}\circ\Psi_{h}^{-1}$ .

3 Stability

For the analysis in this and the next section, we assume that the integrals over cut elements in $\Omega^{\pm}$ are computed exactly. In addition, we restrict our attention to the choice $\alpha=0$ and $\beta=1$ for the averages in (2.17)–(2.18), assuming $\mu_{-}\leq\mu_{+}$ .

The key for the stability analysis of the two-phase Stokes problem is an inf-sup stability property of the unfitted generalized Taylor–Hood finite element pair, which extends the classical LBB stability result for the standard $\mathbf{P}_{k+1}-P_{k}$ Stokes element from [7]. There is no similar stability result in the literature for $\mathbf{Q}_{k+1}-Q_{k}$ unfitted elements. However, we expect that the extension, and so the analysis below, can be carried over to these elements as well.

One is interested in the inf-sup inequality with a stability constant that is independent of the viscosity ratio, position of $\Gamma$ with respect to the background mesh and, of course, mesh size $h$ . The result is given in the following lemma.

Lemma 3.1

Denote by $V_{h}$ the space of continuous $P_{k+1}$ finite element vector functions on $\Omega$ , $V_{h}=\{\mathbf{u}\in C(\Omega)^{d}:\mathbf{u}|_{T}\in\mathbf{P}_{k+1}(T)~{}\forall T\in\mathcal{T}_{h}\}$ . There exists $h_{0}>0$ such that for all $h<h_{0}$ and any $q_{h}\in Q_{h}^{\pm}$ there exists $\mathbf{v}_{h}\in V_{h}$ such that it holds

\begin{split}\mu_{-}^{-1}\|q_{h}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{h}^{+}\|_{{\Omega}_{h}^{+}}^{2}&\leq(q_{h}^{-},\nabla\cdot\mathbf{v}_{h})_{{\Omega}^{-}}+(q_{h}^{+},\nabla\cdot\mathbf{v}_{h})_{{\Omega}^{+}}+c\,b_{p}(q_{h},q_{h})\\ \|\mu^{\frac{1}{2}}\nabla\mathbf{v}_{h}\|_{\Omega}^{2}&\leq C\,\left(\mu_{-}^{-1}\|q_{h}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{h}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right).\end{split}

(3.1)

with $h_{0}$ and two positive constants $c$ and $C$ independent of $q_{h}$ , $\mu_{\pm}$ , the position of $\Gamma$ in the background mesh and mesh size $h$ .

Proof:

Consider subdomains $\Omega_{h,i}^{\pm}\subset\Omega^{\pm}$ built of all strictly internal simplexes in each phase:

\overline{\Omega}_{h,i}^{\pm}:=\bigcup\{\overline{T}\,:\,T\in\mathcal{T}_{h},\quad T\subset\Omega^{\pm}\}.

The following two results are central for the proof. First, we have the uniform inf-sup inequalities in ${\Omega}_{h,i}^{-}$ and ${\Omega}_{h,i}^{+}$ [22]: there exist constants $C_{\pm}$ independent of the position of $\Gamma$ and $h$ such that

0<C_{\pm}\leq\inf_{q\in Q_{h}^{\pm}\cap L^{2}_{0}({\Omega}_{h,i}^{\pm})}~{}\sup_{\scriptsize\begin{array}[]{c}\mathbf{v}\in V_{h}\\ {\rm supp}(\mathbf{v})\subset{\Omega}_{h,i}^{\pm}\end{array}}\frac{(q,\nabla\cdot\mathbf{v})_{{\Omega}_{h,i}^{\pm}}}{\|\mathbf{v}\|_{H^{1}({\Omega}_{h,i}^{\pm})}\|q\|_{{\Omega}_{h,i}^{\pm}}}.

(3.2)

The above result can be equivalently formulated as follows: For any $q\in Q_{h}^{\pm}\cap L^{2}_{0}({\Omega}_{h,i}^{\pm})$ there exist $\mathbf{v}_{h}^{\pm}\in V_{h}$ such that ${\rm supp}(\mathbf{v})\subset{\Omega}_{h,i}^{\pm}$ and

\|q^{\pm}\|_{{\Omega}_{h,i}^{\pm}}^{2}=\left(q^{\pm},\nabla\cdot\mathbf{v}_{h}^{\pm}\right)_{{\Omega}^{\pm}_{h}},\quad\|\nabla\mathbf{v}_{h}^{\pm}\|_{\Omega}\leq C_{\pm}^{-1}\|q^{\pm}\|_{{\Omega}_{h,i}^{\pm}}.

(3.3)

The second important results is the simple observation that the $L^{2}$ norm of $q_{h}$ in $\Omega_{h}^{\pm}$ can be controlled by the $L^{2}$ norm in ${\Omega}_{h,i}^{\pm}$ plus the stabilization term in (2.26) (see, [33, 39]):

\|q_{h}\|_{{\Omega}_{h}^{\pm}}^{2}\leq C\,(\|q_{h}\|_{{\Omega}_{h,i}^{\pm}}^{2}+J_{h}^{\pm}(q_{h},q_{h})),

(3.4)

with some constant $C$ independent of the position of $\Gamma$ and $h$ . We note that (3.4) holds also for discontinuous finite elements.

Consider now

q_{\mu}=\left\{\begin{array}[]{r}~{}\mu_{-}|\Omega^{-}|^{-1}\in Q_{h}^{-}\\ -\mu_{+}|\Omega^{+}|^{-1}\in Q_{h}^{+}.\end{array}\right.

Note that $q_{\mu}$ satisfies the orthogonality condition imposed for elements from $Q_{h}^{\pm}$ , and hence $\mbox{span}\{q_{\mu}\}$ is a subspace in $Q_{h}^{\pm}$ . Using a trick from [37], we decompose arbitrary $q_{h}\in Q_{h}^{\pm}$ into a component collinear with $q_{\mu}$ and the orthogonal complement in each phase:

q_{h}=q_{1}+q_{0},\quad\text{with}~{}q_{1}\in\mbox{span}\{q_{\mu}\},\quad\text{and}~{}~{}(q_{0}^{-},1)_{\Omega_{h,i}^{-}}=(q_{0}^{+},1)_{\Omega_{h,i}^{+}}=0.

Thus, $q_{1}$ and $q_{0}$ are orthogonal with respect to $L^{2}$ product in the inner domains $\Omega_{h,i}^{\pm}$ . Next, we let $q^{\pm}=\mu^{-\frac{1}{2}}_{\pm}q_{0}^{\pm}$ in (3.3) and for $\mathbf{v}^{\pm}_{h}\in V_{h}$ given by (3.3) consider $\mathbf{v}_{h}^{0}=\mu_{-}^{\frac{1}{2}}\mathbf{v}^{-}_{h}+\mu_{+}^{\frac{1}{2}}\mathbf{v}^{+}_{h}\in V_{h}$ . Then after applying (3.4) and summing up, the relations in (3.3) become

\begin{split}\mu_{-}^{-1}\|q_{0}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{0}^{+}\|_{{\Omega}_{h}^{+}}^{2}&\leq C\,\left((q_{0}^{-},\nabla\cdot\mathbf{v}_{h}^{0})_{{\Omega}^{-}}+(q_{0}^{+},\nabla\cdot\mathbf{v}_{h}^{0})_{{\Omega}^{+}}+b_{p}(q_{0},q_{0})\right),\\ \|\mu^{\frac{1}{2}}\nabla\mathbf{v}_{h}^{0}\|_{\Omega}&\leq C_{0}\left(\mu_{-}^{-1}\|q_{0}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{0}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)^{\frac{1}{2}},\end{split}

(3.5)

with $C$ from (3.4) and $C_{0}=\max\{C_{-}^{-1},C_{+}^{-1}\}$ , both of which are independent of $\mu_{\pm}$ and how $\Gamma$ overlaps the background mesh. In (3.5), we also used the fact that supports of $\mathbf{v}^{-}$ and $\mathbf{v}^{+}$ do not overlap. Since ${\rm supp}(\mathbf{v}_{h}^{\pm})\subset{\Omega}^{\pm}$ and $q_{1}^{\pm}$ are constant in ${\Omega}^{\pm}$ , integration by parts shows that

(q_{1}^{\pm},\nabla\cdot\mathbf{v}_{h}^{0})_{{\Omega}^{\pm}_{h}}=0.

(3.6)

Next, we need the following result from Lemma 5.1 in [30]: For all $h\leq h_{0}$ there exists $\mathbf{v}_{h}^{1}\in V_{h}$ such that

\begin{split}\mu_{-}^{-1}\|q_{1}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{1}^{+}\|_{{\Omega}_{h}^{+}}^{2}&=(q_{1},\nabla\cdot\mathbf{v}_{h}^{1})_{{\Omega}^{-}}+(q_{1},\nabla\cdot\mathbf{v}_{h}^{1})_{{\Omega}^{+}},\\ \|\mu^{\frac{1}{2}}\nabla\mathbf{v}_{h}^{1}\|_{\Omega}&\leq C_{1}\left(\mu_{-}^{-1}\|q_{1}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{1}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)^{\frac{1}{2}},\end{split}

(3.7)

with $h_{0}>0$ and $C_{1}>0$ independent of $\mu_{\pm}$ and how $\Gamma$ overlaps the background mesh. The above result follows from the classical inf-sup stability condition for $\mathbf{P}_{2}-P_{1}$ Taylor–Hood elements and a simple scaling and interpolation argument. See [30] for details.

As the next step, set $\mathbf{v}_{h}=\tau\mathbf{v}_{h}^{0}+\mathbf{v}_{h}^{1}$ with some $\tau>0$ and proceed with calculations using (3.6), (3.5), (3.7), and the Cauchy-Schwartz inequality:

\begin{split}(q_{h}^{-},&\nabla\cdot\mathbf{v}_{h})_{{\Omega}^{-}}+(q_{h}^{+},\nabla\cdot\mathbf{v}_{h})_{{\Omega}^{+}}\\ &=(q_{1}^{-},\nabla\cdot\mathbf{v}_{h}^{1})_{{\Omega}^{-}}+(q_{1}^{+},\nabla\cdot\mathbf{v}_{h}^{1})_{{\Omega}^{+}}+\tau(q_{0}^{-},\nabla\cdot\mathbf{v}_{h}^{0})_{{\Omega}^{-}}+\tau(q_{0}^{+},\nabla\cdot\mathbf{v}_{h}^{0})_{{\Omega}^{+}}\\ &\qquad+(q_{0}^{-},\nabla\cdot\mathbf{v}_{h}^{1})_{{\Omega}^{-}}+(q_{0}^{+},\nabla\cdot\mathbf{v}_{h}^{1})_{{\Omega}^{+}}\\ &\geq\mu_{-}^{-1}\|q_{1}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{1}^{+}\|_{{\Omega}_{h}^{+}}^{2}+\tau C^{-1}\left(\mu_{-}^{-1}\|q_{0}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{0}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)-\tau b_{p}(q_{0},q_{0})\\ &\qquad-\left(\mu_{-}^{-1}\|q_{0}^{-}\|_{{\Omega}^{-}_{h}}^{2}+\mu_{+}^{-1}\|q_{0}^{+}\|_{{\Omega}^{+}_{h}}^{2}\right)^{\frac{1}{2}}d^{\frac{1}{2}}\|\mu^{\frac{1}{2}}\nabla\mathbf{v}_{h}^{1}\|_{\Omega}\\ &\geq\mu_{-}^{-1}\|q_{1}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{1}^{+}\|_{{\Omega}_{h}^{+}}^{2}+\tau C^{-1}\left(\mu_{-}^{-1}\|q_{0}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{0}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)-\tau b_{p}(q_{0},q_{0})\\ &\qquad-\left(\mu_{-}^{-1}\|q_{0}^{-}\|_{{\Omega}^{-}_{h}}^{2}+\mu_{+}^{-1}\|q_{0}^{+}\|_{{\Omega}^{+}_{h}}^{2}\right)^{\frac{1}{2}}\,C_{1}d^{\frac{1}{2}}\left(\mu_{-}^{-1}\|q_{1}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{1}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)^{\frac{1}{2}}\\ &\geq\frac{1}{2}\left(\mu_{-}^{-1}\|q_{1}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{1}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)+\left(\frac{\tau}{C}-\frac{C_{1}^{2}d}{2}\right)\left(\mu_{-}^{-1}\|q_{0}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{0}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)-\tau b_{p}(q_{0},q_{0}).\end{split}

We set $\tau$ such that $\frac{\tau}{C}-\frac{C_{1}^{2}d}{2}=\frac{1}{2}$ and note that $b_{p}(q_{0},q_{0})=b_{p}(q_{h},q_{h})$ . Using this and the orthogonality condition for $q_{0}$ , we get

\begin{split}(q_{h},&\nabla\cdot\mathbf{v}_{h})_{{\Omega}^{-}}+(q_{h},\nabla\cdot\mathbf{v}_{h})_{{\Omega}^{+}}\\ &\geq\frac{1}{2}\left(\mu_{-}^{-1}\|q_{1}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{1}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)+\frac{1}{2}\left(\mu_{-}^{-1}\|q_{0}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{0}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)-\tau b_{p}(q_{h},q_{h})\\ &\geq\frac{1}{2}\left(\mu_{-}^{-1}\|q_{1}^{-}\|_{{\Omega}_{h,i}^{-}}^{2}+\mu_{+}^{-1}\|q_{1}^{+}\|_{{\Omega}_{h,i}^{+}}^{2}\right)+\frac{1}{2}\left(\mu_{-}^{-1}\|q_{0}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{0}^{+}\|_{{\Omega}_{h,i}^{+}}^{2}\right)-\tau b_{p}(q_{h},q_{h})\\ &=\frac{1}{2}\left(\mu_{-}^{-1}\|q_{h}^{-}\|_{{\Omega}_{h,i}^{-}}^{2}+\mu_{+}^{-1}\|q_{h}^{+}\|_{{\Omega}_{h,i}^{+}}^{2}\right)-\tau b_{p}(q_{h},q_{h})\\ &\geq\frac{1}{2C}\left(\mu_{-}^{-1}\|q_{h}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{h}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)-\left(\tau+\frac{1}{2}\right)b_{p}(q_{h},q_{h}).\end{split}

(3.8)

Using $(q_{0}^{\pm},q_{1}^{\pm})_{{\Omega}_{h,i}^{\pm}}=0$ , $|\Omega^{\pm}_{h}\setminus\Omega^{\pm}_{h,i}|\leq c\,h$ and so $\|q_{1}^{\pm}\|_{\Omega^{\pm}_{h}\setminus\Omega^{\pm}_{h,i}}\leq ch^{\frac{1}{2}}\|q_{1}^{\pm}\|_{\Omega^{\pm}_{h}}$ , we estimate

|\mu_{-}^{-1}(q_{0}^{-},q_{1}^{-})_{{\Omega}_{h}^{-}}+\mu_{+}^{-1}(q_{0}^{+},q_{1}^{+})_{{\Omega}_{h}^{+}}|\\ \leq c\,h^{\frac{1}{2}}\left(\mu_{-}^{-1}\|q_{0}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{0}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)^{\frac{1}{2}}\left(\mu_{-}^{-1}\|q_{1}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{1}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)^{\frac{1}{2}}.

(3.9)

From (3.5), (3.7), and (3.9), we also get the following upper bound for $\mathbf{v}_{h}$ ,

\begin{split}\|\mu^{\frac{1}{2}}\nabla\mathbf{v}_{h}\|_{\Omega}^{2}&\leq 2(\|\mu^{\frac{1}{2}}\tau\nabla\mathbf{v}_{h}^{0}\|_{\Omega}^{2}+\|\mu^{\frac{1}{2}}\nabla\mathbf{v}_{h}^{1}\|_{\Omega}^{2})\\ &\leq 2\tau^{2}C_{0}^{2}\left(\mu_{-}^{-1}\|q_{0}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{0}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)+2C_{1}^{2}\left(\mu_{-}^{-1}\|q_{1}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{1}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)\\ &\leq\frac{2\max\{\tau^{2}C_{0}^{2},C_{1}^{2}\}}{1-c\,h^{\frac{1}{2}}}\left(\mu_{-}^{-1}\|q_{h}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q_{h}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right).\end{split}

(3.10)

The assertion of the lemma follows from (3.8) and (3.10) after simple calculations.
$\square$

The next lemma shows the uniform coercivity of the symmetric form $a_{h}(\mathbf{u}_{h},\mathbf{v}_{h})$ in (2.21) on $V_{h}^{\pm}\times V_{h}^{\pm}$ .

Lemma 3.2

If $\gamma=O(1)$ in (2.23) is sufficiently large, then it holds

a_{h}(\mathbf{u}_{h},\mathbf{u}_{h})\geq C\,\left(\mu_{-}\|\mathbf{D}(\mathbf{u}_{h}^{-})\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}\|\mathbf{D}(\mathbf{u}_{h}^{+})\|_{{\Omega}_{h}^{+}}^{2}+h^{-1}\|\{\mu\}[\mathbf{u}_{h}\cdot\mathbf{n}]\|_{\Gamma}^{2}+f\|[\mathbf{P}\mathbf{u}_{h}]\|_{\Gamma}^{2}\right)

(3.11)

$\forall~{}\mathbf{u}_{h}\in V_{h}^{\pm}$ , with $C>0$ independent of $\mu_{\pm}$ , $h$ , f, and the position of $\Gamma$ with respect to the background mesh.

Proof:

For the proof, we need the local trace inequality in $T\in\mathcal{T}^{\Gamma}_{h}$ (see, e.g. [22, 23]):

\|v\|_{T\cap\Gamma}\leq C(h_{T}^{-\frac{1}{2}}\|v\|_{T}+h_{T}^{\frac{1}{2}}\|\nabla v\|_{T}),\quad~{}~{}\forall~{}v\in H^{1}(T),

(3.12)

with a constant $C$ independent of $v$ , $T$ , how $\Gamma$ intersects $T$ , and $h_{T}<h_{0}$ for some arbitrary but fixed $h_{0}$ . We also need the following estimate

\|\mathbf{D}(\mathbf{v}_{h}^{\pm})\|^{2}_{L^{2}(\Omega_{h}^{\pm})}\leq C(\|\mathbf{D}(\mathbf{v}_{h}^{\pm})\|^{2}_{L^{2}(\Omega^{\pm})}+\mathbf{J}^{\pm}(\mathbf{v}_{h}^{\pm},\mathbf{v}_{h}^{\pm})\,),

(3.13)

which follows from (3.4) by applying it componentwise and further using FE inverse inequality (note $h^{-2}$ scaling in the definition of $\mathbf{J}^{\pm}$ in (2.24)). Applying (3.12), finite element inverse inequalities and (3.13), we can bound the interface term

\begin{split}\langle\{\mu\mathbf{n}^{T}\mathbf{D}(\mathbf{v}_{h})\mathbf{n}\},&[\mathbf{u}_{h}\cdot\mathbf{n}]\rangle_{\Gamma}=\langle\mu_{-}\mathbf{n}^{T}\mathbf{D}(\mathbf{v}_{h}^{-})\mathbf{n},[\mathbf{u}_{h}\cdot\mathbf{n}]\rangle_{\Gamma}\\ &\leq\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\left[\frac{h_{T}\delta}{2}\|\mu_{-}^{\frac{1}{2}}\mathbf{n}^{T}\mathbf{D}(\mathbf{v}_{h}^{-})\mathbf{n}\|^{2}_{T\cap\Gamma}+\frac{1}{2h_{T}\delta}\|\mu_{-}^{\frac{1}{2}}[\mathbf{u}_{h}\cdot\mathbf{n}]\|^{2}_{T\cap\Gamma}\right]\\ &\leq\frac{\delta}{2}\|\mu_{-}^{\frac{1}{2}}\mathbf{n}^{T}\mathbf{D}(\mathbf{v}_{h}^{-})\mathbf{n}\|^{2}_{\Omega_{h}^{-}}+\frac{1}{h_{T}\delta}\{\mu\}|[\mathbf{u}_{h}\cdot\mathbf{n}]|^{2}_{\Gamma},\quad\forall~{}\delta>0,~{}\mathbf{u}_{h},\mathbf{v}_{h}\in V_{h}^{\pm}.\end{split}

This estimate with $\mathbf{v}_{h}=\mathbf{u}_{h}$ and with $\delta>0$ sufficiently small, together with the definition of the bilinear form $a_{h}(\mathbf{u}_{h},\mathbf{u}_{h})$ , allows to show its coercivity. $\square$

We further need the continuity result for the velocity stabilization form contained in the next lemma.

Lemma 3.3

It holds

a_{p}(\mathbf{v}_{h},\mathbf{v}_{h})\leq C\,\left(\mu_{-}\|\mathbf{D}(\mathbf{v}^{-}_{h})\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}\|\mathbf{D}(\mathbf{v}^{+}_{h})\|_{{\Omega}_{h}^{+}}^{2}\right)\qquad\forall~{}\mathbf{v}_{h}\in V^{\pm}_{h},

with $C>0$ independent of $\mu_{\pm}$ , $h$ , and the position of $\Gamma$ in the background mesh.

Proof:

For any $\mathbf{v}=\mathbf{v}_{h}^{-}\in V^{-}_{h}$ , facet $e\in\mathcal{E}_{h}^{\Gamma,-}$ and the corresponding patch $\omega_{e}$ formed by two tetrahedra $T_{1}$ and $T_{2}$ , it holds

\|\mathbf{v}_{1}^{e}-\mathbf{v}_{2}^{e}\|^{2}_{\omega_{e}}=\|\mathbf{v}_{1}-\mathbf{v}_{2}^{e}\|^{2}_{T_{1}}+\|\mathbf{v}_{1}^{e}-\mathbf{v}_{2}\|^{2}_{T_{2}}\leq(1+c)\|\mathbf{v}_{1}-\mathbf{v}_{2}^{e}\|^{2}_{T_{1}},

where the constant $c$ depends only on shape regularity of the tetrahedra, since $\mathbf{v}_{1}^{e}-\mathbf{v}_{2}$ on $T_{2}$ is the canonical polynomial extension of $\mathbf{v}_{1}-\mathbf{v}_{2}^{e}$ from $T_{1}$ .

Now, we need the following local Korn’s inequality:

\|\nabla\mathbf{v}\|_{T}\leq C\|\mathbf{D}(\mathbf{v})\|_{T},\quad\forall~{}\mathbf{v}\in H^{1}(T)^{d},~{}~{}\text{s.t.}~{}\mathbf{v}=0~{}\text{on any face of}~{}T\in\mathcal{T}_{h},

(3.14)

where $C$ depends only on shape regularity of $T$ . The result in (3.14) follows from eq. (3.3) in [9] and the observation that vector fields vanishing on any face $T$ support only zero rigid motions. A simple scaling argument also proves the local Poincare inequality:

\|\mathbf{v}\|_{T}\leq Ch^{2}_{T}\|\nabla\mathbf{v}\|_{T},\quad\forall~{}\mathbf{v}\in H^{1}(T)^{d},~{}~{}\text{s.t.}~{}\mathbf{v}=0~{}\text{on any face of}~{}T\in\mathcal{T}_{h},

(3.15)

where $C$ depends only on shape regularity of $T$ . Applying (3.14), (3.15) and triangle inequalities on $T_{1}$ for $\mathbf{v}_{1}-\mathbf{v}_{2}^{e}$ which vanishes on $e$ (a face of $T_{1}$ ), we obtain:

	$\displaystyle\\|\mathbf{v}_{1}-\mathbf{v}_{2}^{e}\\|^{2}_{T_{1}}$	$\displaystyle\leq C_{p}h^{2}\\|\mathbf{D}(\mathbf{v}_{1}-\mathbf{v}_{2}^{e})\\|^{2}_{T_{1}}\leq 2C_{p}h^{2}(\\|\mathbf{D}\mathbf{v}_{1}\\|^{2}_{T_{1}}+\\|\mathbf{D}\mathbf{v}_{2}^{e}\\|^{2}_{T_{1}})$
		$\displaystyle\leq 2C_{p}h^{2}(\\|\mathbf{D}\mathbf{v}_{1}\\|^{2}_{T_{1}}+c\,\\|\mathbf{D}\mathbf{v}_{2}\\|^{2}_{T_{2}}),$		(3.16)

where for the last inequality we again use shape regularity and the fact that $\mathbf{D}\mathbf{v}_{2}^{e}=(\mathbf{D}\mathbf{v}_{2})^{e}$ . Thus, we see that $\|\mathbf{v}_{1}^{e}-\mathbf{v}_{2}^{e}\|^{2}_{\omega_{e}}\leq c\,h^{2}\|\mathbf{D}\mathbf{v}\|^{2}_{\omega_{e}}$ , with some $c$ depending only on shape regularity. Summing up over all $e\in\mathcal{E}_{h}^{\Gamma,-}$ leads to the required upper bound for $\mathbf{J}_{h}^{-}(\mathbf{v},\mathbf{v})$ : $\mathbf{J}_{h}^{-}(\mathbf{v},\mathbf{v})\leq C\|\mathbf{D}(\mathbf{v})\|_{\Omega_{h}^{-}}$ . Repeating the same argument for the edges in $\mathcal{E}_{h}^{\Gamma,+}$ and summing up the two bounds scaled by viscosity coefficients proves the lemma.
$\square$

The finite element problem (2.20) can be equivalently formulated as follows: Find $\{\mathbf{u}_{h},p_{h}\}\in V^{\pm}_{h}\times Q^{\pm}_{h}$ such that

\mathcal{A}(\mathbf{u}_{h},p_{h};\mathbf{v}_{h},q_{h})=r_{h}(\mathbf{v}_{h}),\quad\forall~{}\{\mathbf{v}_{h},q_{h}\}\in V^{\pm}_{h}\times Q^{\pm}_{h}

(3.17)

with

\mathcal{A}(\mathbf{u}_{h},p_{h};\mathbf{v}_{h},q_{h})=a_{h}(\mathbf{u}_{h},\mathbf{v}_{h})+b_{h}(\mathbf{v}_{h},p_{h})-b_{h}(\mathbf{u}_{h},q_{h})+b_{p}(p_{h},q_{h}).

Lemmas 3.1–3.3 enable us to show the inf-sup stability of the bilinear form $\mathcal{A}$ . The stability result is formulated using the following composite norm:

\|\mathbf{v},q\|^{2}:=\mu_{-}\|\mathbf{D}(\mathbf{v}^{-})\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}\|\mathbf{D}(\mathbf{v}^{+})\|_{{\Omega}_{h}^{+}}^{2}+h^{-1}\|\{\mu\}[\mathbf{v}\cdot\mathbf{n}]\|_{\Gamma}^{2}+f\|[\mathbf{P}\mathbf{v}]\|_{\Gamma}^{2}+\mu_{-}^{-1}\|q^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|q^{+}\|_{{\Omega}_{h}^{+}}^{2}

for $\mathbf{v}\in V^{\pm}_{h},~{}q\in Q^{\pm}_{h}$ .

Theorem 3.4

There exists $h_{0}>0$ such that for all $h<h_{0}$ it holds

\sup_{\{\mathbf{v}_{h},q_{h}\}\in V^{\pm}_{h}\times Q^{\pm}_{h}}\frac{\mathcal{A}(\mathbf{u}_{h},p_{h};\mathbf{v}_{h},q_{h})}{\|\mathbf{v}_{h},q_{h}\|}\geq C\,\|\mathbf{u}_{h},p_{p}\|,\quad\forall\,\{\mathbf{u}_{h},p_{h}\}\in V^{\pm}_{h}\times Q^{\pm}_{h},

with $h_{0}>0$ and $C>0$ independent of $\mu_{\pm}$ , $h$ , f, and the position of $\Gamma$ in the background mesh.

Proof:

For a given $p_{h}\in Q^{\pm}_{h}$ , Lemma 3.1 implies the existence of such $\mathbf{w}_{h}\in V_{h}$ that

b_{h}(\mathbf{w}_{h},p_{h})+b_{p}(p_{h},q_{h})\geq c\,\left(\mu_{-}^{-1}\|p_{h}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|p_{h}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)

(3.18)

and

\|\mu^{\frac{1}{2}}\nabla\mathbf{w}_{h}\|_{\Omega}^{2}\leq C\,\left(\mu_{-}^{-1}\|p_{h}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|p_{h}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right),

(3.19)

with some positive $c$ , $C$ independent of $\mu$ and how $\Gamma$ overlaps the background mesh. Next, we extend the finite element function $\mathbf{w}_{h}\in V_{h}$ to the element of the product space $\widehat{\mathbf{w}}_{h}\in V_{h}^{\pm}$ by setting $\widehat{\mathbf{w}}_{h}^{\pm}=\mathbf{w}_{h}|_{\Omega_{h}^{\pm}}\in V_{h}^{\pm}$ . We let $\mathbf{v}_{h}=\mathbf{u}_{h}+\tau\widehat{\mathbf{w}}_{h}$ for some $\tau>0$ and $q_{h}=p_{h}$ . Using the definition of the form $\mathcal{A}$ and (3.18), we calculate

\begin{split}\mathcal{A}(\mathbf{u}_{h}&,p_{h};\mathbf{v}_{h},q_{h})=a_{h}(\mathbf{u}_{h},\mathbf{u}_{h})+\tau a_{h}(\mathbf{u}_{h},\widehat{\mathbf{w}}_{h})+\tau b_{h}(\widehat{\mathbf{w}}_{h},p_{h})+b_{p}(p_{h},p_{h})\\ &\geq\frac{1}{2}a_{h}(\mathbf{u}_{h},\mathbf{u}_{h})-\frac{\tau^{2}}{2}a_{h}(\widehat{\mathbf{w}}_{h},\widehat{\mathbf{w}}_{h})+\min\{\tau,1\}\,c\,\left(\mu_{-}^{-1}\|p_{h}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|p_{h}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right),\end{split}

(3.20)

where we used the Cauchy-Schwartz inequality:

\tau a_{h}(\mathbf{u}_{h},\widehat{\mathbf{w}}_{h})\leq\tau|a_{h}(\mathbf{u}_{h},\mathbf{u}_{h})|^{\frac{1}{2}}|a_{h}(\widehat{\mathbf{w}}_{h},\widehat{\mathbf{w}}_{h})|^{\frac{1}{2}}\leq\frac{1}{2}a_{h}(\mathbf{u}_{h},\mathbf{u}_{h})+\frac{\tau^{2}}{2}a_{h}(\widehat{\mathbf{w}}_{h},\widehat{\mathbf{w}}_{h}).

Note that it holds $[\widehat{\mathbf{w}}_{h}\cdot\mathbf{n}]=0$ and $[\mathbf{P}\widehat{\mathbf{w}}_{h}]=0$ on $\Gamma$ . Since all Nitsche and ‘friction’ terms in $a_{h}(\widehat{\mathbf{w}}_{h},\widehat{\mathbf{w}}_{h})$ vanish, the results of the Lemma 3.3 and estimate (3.19) imply the upper bound

a_{h}(\widehat{\mathbf{w}}_{h},\widehat{\mathbf{w}}_{h})\leq C\,\|\mu^{\frac{1}{2}}\nabla\widehat{\mathbf{w}}_{h}\|_{\Omega}^{2}\leq C\,\left(\mu_{-}^{-1}\|p_{h}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|p_{h}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right).

Using it in (3.20) and choosing $\tau>0$ small enough, but independent of all problem parameters, leads us to the lower bound

\mathcal{A}(\mathbf{u}_{h},p_{h};\mathbf{v}_{h},q_{h})\geq\frac{1}{2}a_{h}(\mathbf{u}_{h},\mathbf{u}_{h})+c\,\left(\mu_{-}^{-1}\|p_{h}^{-}\|_{{\Omega}_{h}^{-}}^{2}+\mu_{+}^{-1}\|p_{h}^{+}\|_{{\Omega}_{h}^{+}}^{2}\right)\geq c\,\|\mathbf{u}_{h},p_{h}\|^{2},

(3.21)

with some $c>0$ independent of $\mu_{\pm}$ , $h$ , and the position of $\Gamma$ in the background mesh. For the last inequality, we used (3.11).

Finally, by the construction of $\mathbf{v}_{h}$ and thanks to (3.19) it is straightforward to see the upper bound:

\|\mathbf{v}_{h},q_{h}\|\leq c\,\|\mathbf{u}_{h},p_{h}\|.

This combined with (3.21) proves the theorem.
$\square$

The stability of the finite element solution in the composite norm immediately follows from (3.17) and Theorem 3.4:

\|\mathbf{u}_{h},p_{h}\|\leq C\,\sup_{\{\mathbf{v}_{h},q_{h}\}\in V^{\pm}_{h}\times Q^{\pm}_{h}}\frac{|r_{h}(\mathbf{v}_{h})|}{\|\mathbf{v}_{h},q_{h}\|},

where on the right-hand side we see the dual norm of the functional $r_{h}$ and constant $C$ , which is independent of the mesh size $h$ , the ratio of the viscosity coefficients $\mu_{\pm}$ , and the position of $\Gamma$ in the background mesh.

4 Error analysis

The stability result shown in Sec. 3 and interpolation properties of finite elements enable us to prove optimal order convergence with uniformly bounded constants.

We assume in this section that the solution to problem (2.1)–(2.8) is piecewise smooth in the following sense: $\mathbf{u}^{\pm}\in H^{k+2}(\Omega^{\pm})^{d}$ and $p^{\pm}\in H^{k+1}(\Omega^{\pm})$ . For the sake of notation, we define the following semi-norm

\|\mathbf{u},p\|_{\ast}=\left(\mu_{-}|\mathbf{u}^{-}|_{H^{k+2}(\Omega^{-})}^{2}+\mu_{+}|\mathbf{u}^{+}|_{H^{k+2}(\Omega^{+})}^{2}+\mu_{-}^{-1}|p^{-}|_{H^{k+1}(\Omega^{-})}^{2}+\mu_{+}^{-1}|p^{+}|_{H^{k+1}(\Omega^{+})}^{2}\right)^{\frac{1}{2}}.

(4.1)

Since we assume $\Gamma$ to be at least Lipschitz, there exist extensions $\mathcal{E}\mathbf{u}^{\pm}$ and $\mathcal{E}p^{\pm}$ of the solution from each phase to $\mathbb{R}^{d}$ such that $\mathcal{E}\mathbf{u}^{\pm}\in H^{k+2}(\mathbb{R}^{d})^{3}$ , $\mathcal{E}p^{\pm}\in H^{k+1}(\mathbb{R}^{d})$ . The corresponding norms are bounded as follows

\|\mathcal{E}\mathbf{u}^{\pm}\|_{H^{k+2}(\mathbb{R}^{d})}\leq C\,\|\mathbf{u}^{\pm}\|_{H^{k+2}(\Omega^{\pm})},\quad\|\mathcal{E}p^{\pm}\|_{H^{k+1}(\mathbb{R}^{d})}\leq C\,\|p^{\pm}\|_{H^{k+1}(\Omega^{\pm})}

(4.2)

See [41]. Denote by $I_{h}\mathbf{u}^{\pm}$ the Scott-Zhang interpolants of $\mathcal{E}\mathbf{u}^{\pm}$ onto $V^{\pm}_{h}$ and $I_{h}\mathbf{u}:=\{I_{h}\mathbf{u}^{-},I_{h}\mathbf{u}^{+}\}$ . Same notation $I_{h}p^{\pm}$ will be used for the Scott-Zhang interpolants of $\mathcal{E}p^{\pm}$ onto $Q^{\pm}_{h}$ . For the pressure interpolants, we can always satisfy the orthogonality condition of $Q^{\pm}_{h}$ by choosing a suitable additive constant in the definition of $p$ .

Applying trace inequality (3.12), standard approximation properties of $I_{h}$ , and bounds (4.2), one obtains the approximation property in the product norm:

\|\mathbf{u}-I_{h}\mathbf{u},p-I_{h}p\|\leq C\,h^{k+1}\|\mathbf{u},p\|_{\ast}.

(4.3)

The following continuity result is the immediate consequence of the Cauchy–Schwatz inequality:

\mathcal{A}(\mathbf{u}-I_{h}\mathbf{u},p-I_{h}p;\,\mathbf{v}_{h},q_{h})\leq C\,\|\mathbf{u}-I_{h}\mathbf{u},p-I_{h}p\|\|\mathbf{v}_{h},q_{h}\|\\ +|\langle\{\mu\mathbf{n}^{T}\mathbf{D}(\mathbf{v}_{h})\mathbf{n}\},[(\mathbf{u}-I_{h}\mathbf{u})\cdot\mathbf{n}]\rangle_{\Gamma}+\langle\{\mu\mathbf{n}^{T}\mathbf{D}(\mathbf{u}-I_{h}\mathbf{u})\mathbf{n}\},[\mathbf{v}_{h}\cdot\mathbf{n}]\rangle_{\Gamma}|,

(4.4)

for all $\{\mathbf{v}_{h},q_{h}\}\in V_{h}^{\pm}\times Q_{h}^{\pm}$ . The last term on the right-hand side in (4.4) needs a special treatment. Applying the Cauchy–Schwatz, inequalities (3.12) and (3.13), FE inverse inequalities and approximation properties of the interpolants, we get

\begin{split}|\langle\{\mu\mathbf{n}^{T}\mathbf{D}(\mathbf{v}_{h})\mathbf{n}\},[(\mathbf{u}-I_{h}\mathbf{u})\cdot\mathbf{n}]\rangle_{\Gamma}|\leq C\,h^{k+1}\|\mathbf{u},0\|_{\ast}\|\mathbf{v}_{h},0\|,\\ |\langle\{\mu\mathbf{n}^{T}\mathbf{D}(\mathbf{u}-I_{h}\mathbf{u})\mathbf{n}\},[\mathbf{v}_{h}\cdot\mathbf{n}]\rangle_{\Gamma}|\leq C\,h^{k+1}\|\mathbf{u},0\|_{\ast}\|\mathbf{v}_{h},0\|.\end{split}

(4.5)

The consistency of the stabilization term is formalized in the estimates that follow from [33, lemma 5.5]: For $p^{-}\in H^{k+1}(\Omega^{-})$ , $\mathbf{u}^{-}\in H^{k+2}(\Omega^{-})^{d}$ , it holds

J_{h}^{-}(p^{-},p^{-})\leq C\,h^{2k+2}\|p^{-}\|_{H^{k+1}(\Omega^{-})}^{2},\quad\mathbf{J}_{h}^{-}(\mathbf{u}^{-},\mathbf{u}^{-})\leq C\,h^{2k+2}\|\mathbf{u}^{-}\|_{H^{k+2}(\Omega^{-})}^{2}.

(4.6)

The above estimates and the stability of the interpolants also imply

\begin{split}J_{h}^{-}(p^{-}-I_{h}p^{-},p^{-}-I_{h}p^{-})&\leq C\,h^{2k+2}{|p^{-}|}_{H^{k+1}(\Omega^{-})}^{2},\\ \mathbf{J}_{h}^{-}(\mathbf{u}^{-}-I_{h}\mathbf{u}^{-},\mathbf{u}^{-}-I_{h}\mathbf{u}^{-})&\leq C\,h^{2k+2}{|\mathbf{u}^{-}|}_{H^{k+2}(\Omega^{-})}^{2}.\end{split}

(4.7)

Similar estimates to (4.6), (4.7) hold for $J_{h}^{+}$ and $\mathbf{J}_{h}^{+}$ with $p^{+}\in H^{k+1}(\Omega^{+})$ , $\mathbf{u}^{+}\in H^{k+2}(\Omega^{+})^{d}$ , which can be combined with suitable weights to yield

b_{p}(p-I_{h}p,p-I_{h}p)+a_{p}(\mathbf{u}-I_{h}\mathbf{u},\mathbf{u}-I_{h}\mathbf{u})\leq C\,h^{2k+2}\|\mathbf{u},p\|_{\ast}^{2}.

(4.8)

Denote the error functions by $\mathbf{e}_{u}=\mathcal{E}\mathbf{u}-\mathbf{u}_{h}$ and $e_{p}=\mathcal{E}p-p_{h}$ . Galerkin orthogonality holds up to the consistency terms

\mathcal{A}(\mathbf{e}_{u},e_{p};\,\mathbf{v}_{h},q_{h})=b_{p}(p-I_{h}p,q_{h})+a_{p}(\mathbf{u}-I_{h}\mathbf{u},\mathbf{v}_{h}),

(4.9)

for all $\mathbf{v}_{h}\in V_{h}^{\pm}$ and $q_{h}\in Q_{h}^{\pm}$ .

The result of Lemma 3.2, (4.8) and the trivial bound $b_{p}(q_{h},q_{h})\leq C\|\mathbf{0},q_{h}\|^{2}$ imply the following estimate for the consistency term on the right-hand side of (4.9):

\begin{split}|b_{p}(p&-I_{h}p,q_{h})+a_{p}(\mathbf{u}-I_{h}\mathbf{u},\mathbf{v}_{h})|\\ &\leq|b_{p}(p-I_{h}p,p-I_{h}p)|^{\frac{1}{2}}|b_{p}(q_{h},q_{h})|^{\frac{1}{2}}+|a_{p}(\mathbf{u}-I_{h}\mathbf{u},\mathbf{u}-I_{h}\mathbf{u})|^{\frac{1}{2}}|a_{p}(\mathbf{v}_{h},\mathbf{v}_{h})|^{\frac{1}{2}}\\ &\leq C\,h^{k+1}\|\mathbf{u},p\|_{\ast}\|\mathbf{v}_{h},q_{h}\|,\end{split}

(4.10)

The optimal order error estimate in the energy norm is given in the next theorem.

Theorem 4.1

For sufficiently regular $\mathbf{u},p$ solving problem (2.1)–(2.8) and $\mathbf{u}_{h},p_{h}$ solving problem (2.20), the following error estimate holds:

\|\mathbf{u}-\mathbf{u}_{h},p-p_{h}\|\leq\,Ch^{k+1}\|\mathbf{u},p\|_{\ast},

(4.11)

with a constant $C$ independent of $h$ , the values of viscosities $\mu_{\pm}$ , slip coefficient $f\geq 0$ , and the position of $\Gamma$ with respect to the triangulation $\mathcal{T}_{h}$ .

Proof:

This result follows by standard arguments (see, for example, section 2.3 in [17]) from the inf-sup stability results of Theorem 3.4, continuity estimates (4.4) and (4.5), Galerkin orthogonality and consistency (4.9)–(4.10), and approximation properties (4.3). $\square$

Remark 4.2

If we consider using isoparametric elements to handle numerical integration over cut cells (see section 2.2.1), then the Sobolev seminorms in the definition of $\|\mathbf{u},p\|_{\ast}$ on the right-hand side in (4.11) should be replaced by the full Sobolev norms of the same order; see the error analysis of the isoparametric unfitted FEM in [34].

5 Numerical results

The aim of the numerical results collected in this section is twofold: (i) support the theoretical results presented in Sec. 4 and (ii) provide evidence of the robustness of the proposed finite element approach with respect to the contrast in viscosity, slip coefficient value, and position of the interface relative to the fixed computational mesh.

For the averages in (2.17)-(2.18), we set $\alpha=0$ and $\beta=1$ for all the numerical experiments since we have $\mu_{-}\leq\mu_{+}$ . Recall that this is the choice for the analysis carried out in Sec. 3 and 4. In addition, we set $\gamma^{\pm}_{\mathbf{u}}=0.05$ , $\gamma_{p}^{\pm}=0.05$ , and $\gamma=40$ . The value of all other parameters will depend on the specific test.

For all the results presented below, we will report the $L^{2}$ error and a weighted $H^{1}$ error for the velocity defined as

\left(2\mu_{-}\|D(\mathbf{u}-\mathbf{u}_{h}^{-})\|_{L^{2}(\Omega^{-})}^{2}+2\mu_{+}\|D(\mathbf{u}-\mathbf{u}_{h}^{+})\|_{L^{2}(\Omega^{+})}^{2}\right)^{\frac{1}{2}},

(5.1)

and a weighted $L^{2}$ error for the pressure defined as

\left(\mu^{-1}_{-}\|p-p_{h}^{-}\|_{L^{2}(\Omega^{-})}^{2}+\mu^{-1}_{+}\|p-p_{h}^{+}\|_{L^{2}(\Omega^{+})}^{2}\right)^{\frac{1}{2}}.

(5.2)

5.1. 2D tests

First, we perform a series of tests in 2D. For all the tests, the domain $\Omega$ is square $[-1,1]\times[-1,1]$ and interface $\Gamma$ is a circle of radius $2/3$ centered at $\mathbf{c}=(c_{1},c_{2})$ . Let $(x,y)=(\tilde{x}-c_{1},\tilde{y}-c_{2})$ , $(\tilde{x},\tilde{y})\in\Omega$ . The exact solution we consider is given by:

	$\displaystyle p^{-}$	$\displaystyle=(x-c_{1})^{3},\hskip 68.28644ptp^{+}=(x-c_{1})^{3}-\frac{1}{2},$		(5.3)
	$\displaystyle\mathbf{u}^{-}$	$\displaystyle=g^{-}(x,y)\left[\begin{array}[]{c}-y\\ x\end{array}\right],\qquad\mathbf{u}^{+}=g^{+}(x,y)\left[\begin{array}[]{c}-y\\ x\end{array}\right],$		(5.8)

where

\displaystyle g^{+}(x,y)=\frac{3}{4\mu_{+}}(x^{2}+y^{2}),\quad g^{-}(x,y)=\frac{3}{4\mu_{-}}(x^{2}+y^{2})+\frac{\mu_{-}-\mu_{+}}{3\mu_{+}\mu_{-}}+\frac{1}{f}.

The forcing terms $\mathbf{f}^{-}$ and $\mathbf{f}^{+}$ are found by plugging the above solution in (2.1). The surface tension coefficient $\sigma$ is set to -0.5. The value of the other physical parameters will be specified for each test.

We impose a Dirichlet condition (2.3) on the entire boundary, where function $\mathbf{g}$ is found from $\mathbf{u}^{+}$ in (5.8).

Spatial convergence. First, we check the spatial accuracy of the finite element method described in Sec. 2.2. The aim is to validate our implementation of the method and support the theoretical findings in Sec. 4. For this purpose, we consider exact solution (5.3)-(5.8) with $\mathbf{c}=\boldsymbol{0}$ (i.e., interface $\Gamma$ is a circle centered at the origin of the axes), viscosities $\mu_{-}=1$ and $\mu_{+}=10$ , and $f=10$ .

We consider structured meshes of quads with six levels of refinement. The initial triangulation has a mesh size $h=1/2$ and all the other meshes are obtained by halving $h$ till $h=1/128$ . We choose to use finite element pairs $\mathbf{Q}_{2}-Q_{1}$ . Fig. 2 shows the velocity vectors colored with the velocity magnitude and the pressure computed with mesh $h=1/128$ . Fig. 3 shows the $L^{2}$ error and weighted $H^{1}$ error (5.1) for the velocity and weighted $L^{2}$ error (5.2) for the pressure against the mesh size $h$ . For the range of mesh sizes under consideration, we observe close to cubic convergence in the $L^{2}$ norm for the velocity and quadratic convergence in the weighted $L^{2}$ norm for the pressure and in the weighted $H^{1}$ norm for the velocity.

Robustness with respect to the viscosity contrast. As mentioned in Sec. 1, the case of high contrast for the viscosities in a two-phase problem is especially challenging from the numerical point of view. To test the robustness of our approach, we consider exact solution (5.3)-(5.8) and fix $\mu_{-}=1$ , while we let $\mu_{+}$ vary from 1 to $10^{8}$ . We set $\mathbf{c}=\boldsymbol{0}$ and $f=10$ .

We consider one of the meshes adopted for the previous sets of simulations (with $h=1/64$ ) and use again $\mathbf{Q}_{2}-Q_{1}$ finite elements. Fig. 4 (left) shows the $L^{2}$ error and weighted $H^{1}$ error (5.1) for the velocity and weighted $L^{2}$ error (5.2) for the pressure against the value of $\mu_{+}$ . We observe that all the errors quickly reach a plateau as the $\mu_{+}/\mu_{-}$ ratio increases, after initially decreasing. These results show that our approach is substantially robust with respect to the viscosity contrast $\mu_{+}/\mu_{-}$ .

Robustness with respect to the slip coefficient. For the next set of simulations, we consider exact solution (5.3)-(5.8) and let the slip coefficient $f$ in (2.6)-(2.7) vary from 1/256 to 256. Notice that the larger $f$ becomes, the closer the two-phase problem gets to the homogeneous model. The other parameters are set as follows: $\mathbf{c}=\boldsymbol{0}$ , $\mu_{-}=1$ , and $\mu_{+}=10$ .

We consider again the structured mesh with mesh size $h=1/64$ and $\mathbf{Q}_{2}-Q_{1}$ finite elements. Fig. 4 (right) shows the $L^{2}$ error and weighted $H^{1}$ error (5.1) for the velocity scaled by the $H^{3}$ norm of $\mathbf{u}$ and weighted $L^{2}$ error (5.2) for the pressure against the value of $f$ . We observe that the scaled weighted $H^{1}$ error for the velocity does not vary substantially as $f$ varies, while the other two errors increase as $f$ decreases. When $f$ goes to zero, the external phase loses its control over tangential motions in the internal fluid on $\Gamma$ , thus allowing for purely rigid rotations in the perfectly circular $\Omega^{-}$ ; see the definition of $\mathbf{u}^{-}$ in (5.8). While the seminorm $\|\mathbf{u},p\|_{\ast}$ appearing on the right-hand side in (4.11) remains the same, the full Sobolev norm $\|\mathbf{u}^{-}\|_{k+2}$ grows as $O(f^{-1})$ . Since we use isoparametric unfitted FE, we indeed see the uniform error bound with respect to $f\to 0$ if we normalize the error by the full Sobolev norm of the solution. See Remark 4.2. Summarizing, the approach proves to be robust in the energy norm as the physical parameter $f$ varies.

Robustness with respect to the position of the interface. We conclude the series of the 2D tests with a set of simulations aimed at checking that our approach is not sensitive to the position of the interface with respect to the background mesh. For this purpose, we vary the center of the circle that represents $\Gamma$ :

\mathbf{c}=(c_{1},c_{2}),~{}c_{1}=\frac{h}{20}k\cos\left(\frac{k}{10}\pi\right),~{}c_{2}=\frac{h}{20}k\sin\left(\frac{k}{10}\pi\right),\quad k=1,2,...,20,

(5.9)

where $h$ is the mesh size. We set $\mu_{-}=1$ , $\mu_{+}=10$ and $f=10$ .

Just like the two previous sets of simulations, we consider the mesh with mesh size $h=1/64$ and the $\mathbf{Q}_{2}-Q_{1}$ pair. Fig. 5 shows the $L^{2}$ error and weighted $H^{1}$ error (5.1) for the velocity and weighted $L^{2}$ error (5.2) for the pressure against the value of $k$ in (5.9). We see that all the errors are fairly insensitive to the position of $\Gamma$ with respect to the background mesh, indicating robustness.

5.2. 3D tests

For the 3D tests, the domain $\Omega$ is cube $[-1.5,1.5]\times[-1.5,1.5]\times[-1.5,1.5]$ and interface $\Gamma$ is the unit sphere, centered at origin of the axes. We characterize $\Gamma$ as the zero level set of function $\phi(\mathbf{x})=||\mathbf{x}||^{2}_{2}-1$ , with $\mathbf{x}=(x,y,z)$ . We consider the exact solution given by:

	$\displaystyle p^{+}$	$\displaystyle=\frac{1}{2}x,\hskip 79.6678ptp^{-}=x,$		(5.10)
	$\displaystyle\mathbf{u}^{-}$	$\displaystyle=g^{-}(x,y)\left[\begin{array}[]{c}-y\\ x\\ 0\end{array}\right],\qquad\mathbf{u}^{+}=g^{+}(x,y)\left[\begin{array}[]{c}-y\\ x\\ 0\end{array}\right],$		(5.17)

where

	$\displaystyle g^{+}(x,y)=\frac{1}{2\mu_{+}}(x^{2}+y^{2}+z^{2}),$
	$\displaystyle g^{-}(x,y)=\frac{1}{2\mu_{-}}(x^{2}+y^{2}+z^{2})+\frac{\mu_{-}-2\mu_{+}\mu_{-}-\mu_{+}}{2\mu_{+}\mu_{-}}.$

The forcing terms $\mathbf{f}^{-}$ and $\mathbf{f}^{+}$ are found by plugging the above solution in in (2.1). We set $f=1$ , $\mu_{-}=1$ , and $\mu_{+}=100$ . The surface tension coefficient is set to $\sigma=-0.5x$ .

Just like for the 2D tests, we impose a Dirichlet condition (2.3) on the entire boundary, where function $\mathbf{g}$ is found from $\mathbf{u}^{+}$ in (5.17).

To verify our implementation of the finite element method in Sec. 2.2 in three dimensions and to further corroborate the results in Sec. 4, we consider structured meshes of tetrahedra with four levels of refinement. The initial triangulation has mesh size $h=1$ and all the other meshes are obtained by halving $h$ till $h=0.125$ . All the meshes feature a local one-level refinement near the corners of $\Omega$ . We choose to use finite element pair $\mathbf{P}_{2}-P_{1}$ . Fig. 6 shows a visualization of the solution computed with mesh $h=0.125$ . Fig. 7 shows the $L^{2}$ error and weighted $H^{1}$ error (5.1) for the velocity and weighted $L^{2}$ error (5.2) for the pressure against the mesh size $h$ . For the small range of mesh sizes that we consider, we observe almost cubic convergence in the $L^{2}$ norm for the velocity, quadratic convergence in the weighted $L^{2}$ norm for the pressure and in the weighted $H^{1}$ norm for the velocity.

6 Conclusions

In this paper, we focused on the two-phase Stokes problem with slip between phases, which has received much less attention than its homogeneous counterpart (i.e. no slip between the phases). For the numerical approximation of this problem, we chose an isoparametric unfitted finite element approach of the CutFEM or Nitsche-XFEM family. For the unfitted generalized Taylor–Hood finite element pair $\mathbf{P}_{k+1}-P_{k}$ , we prove stability and optimal error estimates, which follow from an inf-sup stability property. We show that the inf-sup stability constant is independent of the viscosity ratio, slip coefficient, position of the interface with respect to the background mesh and, of course, mesh size.

The 2D and 3D numerical experiments we used to test our approach feature an exact solution. They have been designed to support the theoretical findings and demonstrate the robustness of our approach for a wide range of physical parameter values. Finally, we show that our unfitted approach is insensitive to the position of the interface between the two phases with respect to the fixed computational mesh.

Acknowledgments

This work was partially supported by US National Science Foundation (NSF) through grant DMS-1953535. M.O. also acknowledges the support from NSF through DMS-2011444. A.Q. also acknowledges the support from NSF through DMS-1620384.

References

[1] Netgen/NGSolve. https://ngsolve.org/.
[2] ngsxfem. https://github.com/ngsxfem/ngsxfem/tree/49205a1ae637771a0ed56d4993ce99008f3a00e0.
[3] Slimane Adjerid, Nabil Chaabane, and Tao Lin. An immersed discontinuous finite element method for stokes interface problems. Computer Methods in Applied Mechanics and Engineering, 293:170 – 190, 2015.
[4] D. M. Anderson, G. B. McFadden, and A. A. Wheeler. Diffuse-interface methods in fluid mechanics. Annual Review of Fluid Mechanics, 30(1):139–165, 1998.
[5] S. Basting and M. Weismann. A hybrid level set/front tracking approach for finite element simulations of two-phase flows. Journal of Computational and Applied Mathematics, 270:471–483, 2014.
[6] Steffen Basting, Annalisa Quaini, Suncica Canic, and Roland Glowinski. Extended ALE method for fluid-structure interaction problems with large structural displacements. Journal of Computational Physics, 331:312 – 336, 2017.
[7] Michel Bercovier and Olivier Pironneau. Error estimates for finite element method solution of the stokes problem in the primitive variables. Numerische Mathematik, 33(2):211–224, 1979.
[8] S.P. Bordas, E. Burman, M.G. Larson, and eds. M. A. Olshanskii. Geometrically Unfitted Finite Element Methods and Applications, volume Lecture Notes in Computational Science and Engineering 121. Springer, Berlin, 2018.
[9] Susanne C Brenner. Korn’s inequalities for piecewise H1 vector fields. Mathematics of Computation, pages 1067–1087, 2004.
[10] E. Burman, G. Delay, and A. Ern. An unfitted hybrid high-order method for the Stokes interface problem. hal-02519896v3, 2020.
[11] Erik Burman. Ghost penalty. C. R. Math. Acad. Sci. Paris, 348(21-22):1217–1220, 2010.
[12] Erik Burman, Susanne Claus, Peter Hansbo, Mats G Larson, and André Massing. Cutfem: Discretizing geometry and partial differential equations. International Journal for Numerical Methods in Engineering, 104(7):472–501, 2015.
[13] Ernesto Cáceres, Johnny Guzmán, and Maxim Olshanskii. New stability estimates for an unfitted finite element method for two-phase stokes problem. SIAM Journal on Numerical Analysis, 58(4):2165–2192, 2020.
[14] J. Chessa and T. Belytschko. An extended finite element method for two-phase fluids. ASME Journal of Applied Mechanics, 70:10–17, 2003.
[15] Susanne Claus and Pierre Kerfriden. A cutfem method for two-phase flow problems. Computer Methods in Applied Mechanics and Engineering, 348:185 – 206, 2019.
[16] Jean Donea, Antonio Huerta, J.-Ph. Ponthot, and A. Rodríguez-Ferran. Arbitrary Lagrangian–Eulerian Methods, chapter 14. John Wiley & Sons Ltd.,, 2004.
[17] A. Ern and J.-L. Guermond. Theory and practice of finite elements, volume 159. Springer, New York, 2013.
[18] Thomas Frachon and Sara Zahedi. A cut finite element method for incompressible two-phase Navier–Stokes flows. Journal of Computational Physics, 384:77 – 98, 2019.
[19] T. P. Fries. The intrinsic XFEM for two-fluid flows. International Journal for Numerical Methods in Fluids, 60(4):437–471, 2009.
[20] P. Gangl, K. Sturm, M. Neunteufel, and J. Schöberl. Fully and semi-automated shape differentiation in NGSolve. arXiv:2004.06783, 2020.
[21] S. Groß, V. Reichelt, and A. Reusken. A finite element based level set method for two-phase incompressible flows. Comp. Visual. Sci., 9:239–257, 2006.
[22] Johnny Guzmán and Maxim Olshanskii. Inf-sup stability of geometrically unfitted stokes finite elements. Mathematics of Computation, 87(313):2091–2112, 2018.
[23] A. Hansbo and P. Hansbo. An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems. Comput. Methods Appl. Mech. Engrg., 191:5537–5552, 2002.
[24] Peter Hansbo, Mats G. Larson, and Sara Zahedi. A cut finite element method for a Stokes interface problem. Applied Numerical Mathematics, 85:90 – 114, 2014.
[25] Jerzy Hapanowicz. Slip between the phases in two-phase water–oil flow in a horizontal pipe. International Journal of Multiphase Flow, 34(6):559 – 566, 2008.
[26] Mohammad R. Hashemi, Pavel B. Ryzhakov, and Riccardo Rossi. An enriched finite element/level-set method for simulating two-phase incompressible fluid flows with surface tension. Computer Methods in Applied Mechanics and Engineering, 370:113277, 2020.
[27] X. He, F. Song, and W. Deng. Stabilized nonconforming Nitsche’s extended finite element method for Stokes interface problems. https://arxiv.org/abs/1905.04844, 2019.
[28] David Jacqmin. Calculation of Two-Phase Navier-Stokes Flows Using Phase-Field Modeling. Journal of Computational Physics, 155(1):96–127, October 1999.
[29] Mohammad J. Kermani and John M. Stockie. The effect of slip velocity on saturation for multiphase condensing mixtures in a pem fuel cell. International Journal of Hydrogen Energy, 36(20):13235 – 13240, 2011. 3rd Iranian Fuel Cell Seminar.
[30] Matthias Kirchhart, Sven Gross, and Arnold Reusken. Analysis of an XFEM discretization for Stokes interface problems. SIAM Journal on Scientific Computing, 38(2):A1019–A1043, 2016.
[31] Christoph Lehrenfeld. High order unfitted finite element methods on level set domains using isoparametric mappings. Computer Methods in Applied Mechanics and Engineering, 300:716–733, 2016.
[32] Christoph Lehrenfeld. A higher order isoparametric fictitious domain method for level set domains. In Stéphane P. A. Bordas, Erik Burman, Mats G. Larson, and Maxim A. Olshanskii, editors, Geometrically Unfitted Finite Element Methods and Applications, pages 65–92, Cham, 2017. Springer International Publishing.
[33] Christoph Lehrenfeld and Maxim Olshanskii. An Eulerian finite element method for PDEs in time-dependent domains. ESAIM: Mathematical Modelling and Numerical Analysis, 53(2):585–614, 2019.
[34] Christoph Lehrenfeld and Arnold Reusken. Analysis of a high-order unfitted finite element method for elliptic interface problems. IMA Journal of Numerical Analysis, 38(3):1351–1387, 2018.
[35] A Massing, M.G. Larson, A. Logg, and M.E. Rognes. A stabilized Nitsche overlapping mesh method for the Stokes problem. Numer. Math., 128:73 – 101, 2014.
[36] Nicolas Moës, John Dolbow, and Ted Belytschko. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 46(1):131–150, 1999.
[37] Maxim A Olshanskii and Arnold Reusken. Analysis of a Stokes interface problem. Numerische Mathematik, 103(1):129–149, 2006.
[38] Elin Olsson and Gunilla Kreiss. A conservative level set method for two phase flow. Journal of Computational Physics, 210(1):225 – 246, 2005.
[39] J. Preuß. Higher order unfitted isoparametric space-time FEM on moving domains. Master’s thesis, NAM, University of Göttingen, 2018.
[40] Henning Sauerland and Thomas-Peter Fries. The stable XFEM for two-phase flows. Computers & Fluids, 87:41 – 49, 2013. USNCCM Moving Boundaries.
[41] Elias M Stein. Singular integrals and differentiability properties of functions, volume 30. Princeton university press, 1970.
[42] Mark Sussman, Peter Smereka, and Stanley Osher. A level set approach for computing solutions to incompressible two-phase flow. Journal of Computational Physics, 114(1):146 – 159, 1994.
[43] S O Unverdi and G Tryggvason. A front-tracking method for viscous, incompressible, multi-fluid flows. Journal of Computational Physics; (United States), 100, 3 1992.
[44] Henry von Wahl, Thomas Richter, and Christoph Lehrenfeld. An unfitted Eulerian finite element method for the time-dependent Stokes problem on moving domains. arXiv preprint arXiv:2002.02352, 2020.
[45] N. Wang and J. Chen. A nonconforming Nitsche’s extended finite element method for Stokes interface problems. J Sci Comput, 81:342–374, 2019.
[46] Qiuliang Wang and Jinru Chen. A new unfitted stabilized Nitsche’s finite element method for Stokes interface problems. Computers & Mathematics with Applications, 70(5):820 – 834, 2015.

	$\displaystyle\\|\mathbf{v}_{1}-\mathbf{v}_{2}^{e}\\|^{2}_{T_{1}}$	$\displaystyle\leq C_{p}h^{2}\\|\mathbf{D}(\mathbf{v}_{1}-\mathbf{v}_{2}^{e})\\|^{2}_{T_{1}}\leq 2C_{p}h^{2}(\\|\mathbf{D}\mathbf{v}_{1}\\|^{2}_{T_{1}}+\\|\mathbf{D}\mathbf{v}_{2}^{e}\\|^{2}_{T_{1}})$
		$\displaystyle\leq 2C_{p}h^{2}(\\|\mathbf{D}\mathbf{v}_{1}\\|^{2}_{T_{1}}+c\,\\|\mathbf{D}\mathbf{v}_{2}\\|^{2}_{T_{2}}),$		(3.16)