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

M. Olshanskii1, A. Quaini1 and Q. Sun1

1Department of Mathematics, University of Houston, 3551 Cullen Blvd, Houston
TX 77204
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

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$.

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.

Figure 2: Approximation of exact solution (5.3)-(5.8) for
$\mathbf{c}=\boldsymbol{0}$, $\mu_{-}=1$, $\mu_{+}=10$, and $f=10$, computed
with mesh $h=1/128$: velocity vectors colored with the velocity magnitude
(left) and pressure (right). Figure 3: 2D test with
$\mathbf{c}=\boldsymbol{0}$, $\mu_{-}=1$, $\mu_{+}=10$, and $f=10$: $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$.

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_{-}$.

Figure 4: 2D test with $\mathbf{c}=\boldsymbol{0}$ and $\mu_{-}=1$: $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_{+}$ (left) and against
the value of the slip coefficient $f$ (right).

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.

Figure 5: 2D test with $\mathbf{c}=\boldsymbol{0}$, $\mu_{-}=1$, $\mu_{+}=10$,
and $f=10$: $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).

### 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.

Figure 6: Approximation of exact solution (5.10)-(5.17) computed with the mesh
with $h=0.125$: velocity vectors colored with the velocity magnitude on the
$xz$-section of $\Omega^{+}$ and in $\Omega^{-}$ (left) and pressure in
$\Omega^{-}$ and half $\Omega^{+}$ (right). Figure 7: 3D test: $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$.

##  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\.


is a non-trivial multiple of $\mathbf{M}^{\alpha}$ where
$\alpha=\frac{1-s/2+s/N}{1-s/2+1/N}$.

To prove this theorem, we start with a simple lemma.

######  Lemma 5.2.

Assume that $W$ satisfies 1–5. Then for every $\delta\in(0,1)$, there exists
$c_{\delta}\in(0,+\infty)$ such that for every $\varepsilon>0$ and every
$u\in\mathbb{R}_{+}$,

(5.10) $\delta(u^{s}\wedge c_{\delta}\varepsilon^{-N(1-s)}u)\leq
W_{\varepsilon}(u).$

###### Proof.

Fix $\delta\in(0,1)$. There exists $M>0$ such that $\delta u^{s}\leq W(u)$ for
every $u\geq M$. Besides, the map $w:u\mapsto W(u)/u$ is lower semicontinuous
and positive on $(0,M]$ by 1 and 2, and since $\liminf_{u\to 0}w(u)>0$ by 5,
$w$ is necessarily bounded from below on $(0,M]$ by some contant $c>0$. As a
consequence $W_{\varepsilon}(u)\geq\delta u^{s}$ if $u\geq\varepsilon^{N}M$
and $W_{\varepsilon}(u)\geq c\varepsilon^{N(s-1)}u$ if
$u\leq\varepsilon^{N}M$, hence:

$\forall u\in\mathbb{R},\quad W_{\varepsilon}(u)\geq\delta(u^{s}\wedge
c\varepsilon^{-N(1-s)}u).\qed$

###### Proof of Theorem 5.1.

By 4, there exists a constant $C$ such that $f_{\varepsilon}^{W}\leq Cf_{s}$
for every $\varepsilon$, and since $f_{\varepsilon}^{W}$ does not depend on
the $x$ variable and converges pointwise to $f_{s}$, 1 is satisfied and our
$\Gamma-\limsup$ result stated in Proposition 4.2 yields

$\mathbf{M}^{H_{f_{s}}}\geq\Gamma(\mathscr{C}_{b}^{\prime})-\limsup_{\varepsilon\to
0}\mathcal{E}_{\varepsilon}.$

Fix $\delta\in(0,1)$. By Lemma 5.2, there exists $c_{\delta}$ such that

$\forall x,u,\xi,\quad
f_{\varepsilon}^{W}(x,u,\xi)\geq\delta(\lvert\xi\rvert^{2}+(u^{s}\wedge
c_{\delta}\varepsilon^{-N(1-s)}u)\eqqcolon f_{\varepsilon}^{\delta}(x,u,\xi).$

It is easy to check that $f^{\delta}_{\varepsilon}$ satisfies 1, 2 and 5 for
every $\varepsilon>0$. Moreover $f_{\varepsilon}^{\delta}\uparrow\delta f_{s}$
and $(f_{\varepsilon}^{\delta})_{-}^{\prime}(\cdot,0^{+},0)=\delta
c_{\delta}\varepsilon^{-N(1-s)}\uparrow(+\infty)=(\delta
f_{s})_{-}^{\prime}(\cdot,0^{+},0)$ as $\varepsilon\to 0$, thus 8 holds for
the family $(f_{\varepsilon}^{\delta})_{\varepsilon>0}$, and by applying our
$\Gamma-\liminf$ result stated in Proposition 4.1 to the energies
$\mathcal{E}_{\varepsilon}^{\delta}$ induced by $f_{\varepsilon}^{\delta}$ we
get:

$\Gamma(\mathscr{C}_{0}^{\prime})-\liminf\mathcal{E}_{\varepsilon}\geq\Gamma(\mathscr{C}_{0}^{\prime})-\liminf\mathcal{E}_{\varepsilon}^{\delta}\geq\mathbf{M}^{H_{\delta
f_{s}}^{-}}.$

We get the result by taking the limit $\delta\to 1$, noticing that
$(f_{s})_{-}^{\prime}(\cdot,0^{+},0)=+\infty$, so that $H_{\delta
f_{s}}^{-}=H_{\delta f_{s}}=\delta H_{f_{s}}$ and $\mathbf{M}^{H_{\delta
f_{s}}^{-}}=\mathbf{M}^{\delta H_{f_{s}}}=\delta\mathbf{M}^{H_{f_{s}}}$. ∎

######  Remark 5.3.

We recover the $\Gamma$-convergence results of [BDS96] and [Dub98] when
$s\in(-2,1)$ under slightly more general assumptions: besides 2 and 3, the
authors impose the existence of a non-trivial slope $\lim_{u\to
0}\frac{W(u)}{u}\in(0,+\infty)$ and a regularity condition (either $W$ is of
class $\mathscr{C}^{1}$ or continuous and non-decreasing close to $0$), which
are stronger than 1, 4 and 5. Let us stress however that these works also
tackle the cases $s<-2$ in any dimension, where the exponent $\rho$ has to be
fixed to $\rho(-2,N)$, and the case $s=-2$ in dimension one, where a
logarithmic factor must be introduced, replacing $\varepsilon^{-\rho}$ with
$\varepsilon^{-\rho(-2,1)}\lvert\log\varepsilon\rvert^{-1}=\varepsilon^{-1/2}\lvert\log\varepsilon\rvert^{-1}$.
This implies that in our model we get a trivial $\Gamma$-limit when $s\leq-2$,
namely $H_{f_{s}}\equiv+\infty$ on $(0,+\infty)$.

color=green!20,bordercolor=green!20,inline, caption=2do, ]PAP4: Dans le cas 1d
de [BDS96], ils font des hypothèses de régularité plus fortes ($W$ est
continue, strictement croissante au voisinage de $0$ et l’existence d’une
pente non triviale en $0$), alors que dans [Dub98] ils supposent que $W$ est
$\mathscr{C}^{1}$, tout en disant qu’on peut supposer à la place que $W$ est
continue est croissante (pas strictement) au voisinage de $0$. Aussi, ils
disent qu’on peut relâcher l’hypothèse de pente par $\liminf_{t\to
0}\frac{W(t)}{t}>0$ (ça on l’obtient aussi avec notre résultat), et plus loin
qu’on peut s’en passer complètement (ça je suis très dubitatif).

## References

  * [AFP00] Luigi Ambrosio, Nicola Fusco and Diego Pallara  “Functions of Bounded Variation and Free Discontinuity Problems”, Oxford Mathematical Monographs  Oxford, New York: Oxford University Press, 2000
  * [AG99] Patricio Aviles and Yoshikazu Giga  “On Lower Semicontinuity of a Defect Energy Obtained by a Singular Limit of the Ginzburg–Landau Type Energy for Gradient Fields”  In _Proceedings of the Royal Society of Edinburgh Section A: Mathematics_ 129.1 Royal Society of Edinburgh Scotland Foundation, 1999, pp. 1–17  DOI: 10.1017/S0308210500027438
  * [BB90] Guy Bouchitté and Giuseppe Buttazzo  “New Lower Semicontinuity Results for Nonconvex Functionals Defined on Measures”  In _Nonlinear Analysis: Theory, Methods & Applications_ 15.7, 1990, pp. 679–692  DOI: 10.1016/0362-546X(90)90007-4
  * [BB93] Guy Bouchitté and Giuseppe Buttazzo  “Relaxation for a Class of Nonconvex Functionals Defined on Measures”  In _Annales de l’Institut Henri Poincaré C, Analyse non linéaire_ 10.3, 1993, pp. 345–361  DOI: 10.1016/S0294-1449(16)30216-5
  * [BBH17] Fabrice Bethuel, Haim Brezis and Frederic Helein  “Ginzburg-Landau Vortices”, Modern Birkhäuser Classics  Cham: Springer International Publishing, 2017  DOI: 10.1007/978-3-319-66673-0
  * [BCM09] Marc Bernot, Vicent Caselles and Jean-Michel Morel  “Optimal Transportation Networks: Models and Theory”, Lecture Notes in Mathematics  Berlin Heidelberg: Springer-Verlag, 2009  DOI: 10.1007/978-3-540-69315-4
  * [BDS96] Guy Bouchitté, Christophe Dubs and Pierre Seppecher  “Transitions de phases avec un potentiel dégénéré à l’infini, application à l’équilibre de petites gouttes”  In _Comptes Rendus de l’Académie des Sciences. Série I_ 323.9 Elsevier, Paris, 1996, pp. 1103–1108
  * [BPP12] Patricia Bauman, Jinhae Park and Daniel Phillips  “Analysis of Nematic Liquid Crystals with Disclination Lines”  In _Archive for Rational Mechanics and Analysis_ 205.3, 2012, pp. 795–826  DOI: 10.1007/s00205-012-0530-7
  * [Bra02] Andrea Braides  “Gamma-Convergence for Beginners”, Oxford Lecture Series in Mathematics and Its Applications  Oxford: Oxford University Press, 2002  DOI: 10.1093/acprof:oso/9780198507840.001.0001
  * [But89] Giuseppe Buttazzo  “Semicontinuity, Relaxation and Integral Representation in the Calculus of Variations”  Harlow: Longman Scientific &| Technical; New York: John Wiley &| Sons, 1989
  * [Col+17] Maria Colombo, Antonio De Rosa, Andrea Marchese and Salvatore Stuvard  “On the Lower Semicontinuous Envelope of Functionals Defined on Polyhedral Chains”  In _Nonlinear Analysis_ 163, 2017, pp. 201–215  DOI: 10.1016/j.na.2017.08.002
  * [DH03] Thierry De Pauw and Robert Hardt  “Size Minimization and Approximating Problems”  In _Calculus of Variations and Partial Differential Equations_ 17.4, 2003, pp. 405–442  DOI: 10.1007/s00526-002-0177-6
  * [Dub98] Christophe Dubs  “Problemes de Perturbations Singulieres Avec Un Potentiel Degenere a l’infini”, 1998  URL: http://www.theses.fr/1998TOUL0005
  * [Fed59] Herbert Federer  “Curvature Measures”  In _Transactions of the American Mathematical Society_ 93.3, 1959, pp. 418–491  DOI: 10.1090/S0002-9947-1959-0110078-1
  * [Fle66] Wendell H. Fleming  “Flat Chains over a Finite Coefficient Group”  In _Transactions of the American Mathematical Society_ 121.1, 1966, pp. 160–186  DOI: 10.1090/S0002-9947-1966-0185084-5
  * [Mar14] Mihai Mariş  “Profile Decomposition for Sequences of Borel Measures”, 2014  arXiv: http://arxiv.org/abs/1410.6125
  * [MM77] Luciano Modica and Stefano Mortola  “Un esempio di $\Gamma^{-}$-convergenza”  In _Bollettino della Unione Matematica Italiana. Series V. B_ 14 Zanichelli, Bologna, 1977, pp. 285–299
  * [Mon15] Antonin Monteil  “Elliptic Approximations of Singular Energies under Divergence Constraint”, 2015  URL: https://tel.archives-ouvertes.fr/tel-01326231
  * [Mon17] Antonin Monteil  “Uniform Estimates for a Modica–Mortola Type Approximation of Branched Transportation”  In _ESAIM: Control, Optimisation and Calculus of Variations_ 23.1 EDP Sciences, 2017, pp. 309–335  DOI: 10.1051/cocv/2015049
  * [Nir59] L. Nirenberg  “On elliptic partial differential equations”  In _Annali della Scuola Normale Superiore di Pisa - Classe di Scienze_ 13.2, 1959, pp. 115–162  URL: http://www.numdam.org/item/ASNSP_1959_3_13_2_115_0/
  * [OS11] Edouard Oudet and Filippo Santambrogio  “A Modica-Mortola Approximation for Branched Transport and Applications”  In _Archive for Rational Mechanics and Analysis_ 201.1, 2011, pp. 115–142  DOI: 10.1007/s00205-011-0402-6
  * [RV73] A.W. Roberts and D.E. Varberg  “Convex Functions” 57, Pure and Applied Mathematics; a Series of Monographs and Textbooks  Academic Press, 1973
  * [San15] Filippo Santambrogio  “Optimal Transport for Applied Mathematicians” 87, Progress in Nonlinear Differential Equations and Their Applications  Cham: Springer International Publishing, 2015  DOI: 10.1007/978-3-319-20828-2
  * [Wir19] Benedikt Wirth  “Phase Field Models for Two-Dimensional Branched Transportation Problems”  In _Calculus of Variations and Partial Differential Equations_ 58.5, 2019  DOI: 10.1007/s00526-019-1615-z