# A finite element method for two-phase flow with material viscous interface

Maxim Olshanskii∗ Annalisa Quaini∗ and Qi Sun Department of Mathematics,
University of Houston, 3551 Cullen Blvd, Houston TX 77204
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
This work has been partially supported by US National Science Foundation (NSF)
through grants DMS-1620384, DMS-1953535, DMS-2011444.

###### Abstract

This paper studies a model of two-phase flow with an immersed material viscous
interface and a finite element method for numerical solution of the resulting
system of PDEs. The interaction between the bulk and surface media is
characterized by no-penetration and slip with friction interface conditions.
The system is shown to be dissipative and a model stationary problem is proved
to be well-posed. The finite element method applied in this paper belongs to a
family of unfitted discretizations. The performance of the method when model
and discretization parameters vary is assessed. Moreover, an iterative
procedure based on the splitting of the system into bulk and surface problems
is introduced and studied numerically.

###### keywords:

two-phase flow, surface Stokes equations, finite elements, cutFEM, TraceFEM

##  1 Introduction

We consider two immiscible, viscous, and incompressible fluids separated by a
viscous inextensible material interface modeled as a Boussinesq–Scriven
surface fluid. The following coupling conditions are prescribed between the
bulk two-phase flow and the surface fluid: (i) the immiscibility condition,
i.e. the bulk fluid does not penetrate through the interface; (ii) slip with
friction between the bulk fluid and the viscous interface; and (iii) the load
exerted from the bulk fluid onto the surface fluid defined by the jump of the
normal stress across the interface. To the best of our knowledge, it is the
first time that the slip with friction conditions are imposed in this context.
This extends the Navier boundary condition to the interface case and
interpolates between no-slip of the bulk phase and embedded material layer, as
considered in [4, 10, 38, 45], and uncoupled lateral dynamics of fluidic
surfaces as studied, for example, in [46, 48, 51]. The model contributes to
characterizing the interaction between a rather viscous and dense interface
material and less viscous and dense bulk fluid, as is the case of, e.g., a
lipid vesicle or cell membrane with its content and surrounding extracellular
fluid.

We first state the problem in terms of incompressible two-phase Navier–Stokes
flow in the bulk coupled to surface fluid equations and show the energy
balance for this system. This balance demonstrates that the system with the
proposed coupling conditions is thermodynamically consistent, i.e. in the
absence of external forces and with no energy inflow through the boundary the
system is dissipative. Then, we introduce a simplifying assumption: the
coupled system has reached a steady state and inertia terms can be neglected.
Although the resulting surface–bulk Stokes problem is a strong simplification
of the original problem, it still retains interesting features for the purpose
of numerical analysis.

For the simplified model, we show well-posedness following the abstract
framework for generalized saddle point problems proposed in, e.g., [5, 37].
Next, we present a sharp interface, geometrically unfitted finite element (FE)
method. Unfitted methods allow the sharp interface to cut through the elements
of a fixed background grid. Their main advantage is the relative ease in
handling time-dependent domains, implicitly defined interfaces, and problems
with strong geometric deformations [9]. For the bulk two-phase Stokes problem,
we apply an isoparametric unfitted finite element approach [32] of the CutFEM
(or Nitsche-XFEM) family[14, 17, 25, 26, 35, 52]. For a (bulk only) interface
Stokes problem with slip between phases, the approach was studied in [40].
CutFEM uses overlapping fictitious domains in combination with ghost penalty
stabilization [11] to enrich and stabilize the solution. In this paper, we
consider the unfitted generalized Taylor–Hood finite element pair
$\mathbf{P}_{k+1}-P_{k}$, $k\geq 1$. For more details on the isoparametric
unfitted finite element, we refer to [33, 34]. For the discretization of the
surface Stokes problem, we apply the trace finite element method [41, 29],
which uses traces of the bulk finite element functions. In a set of numerical
experiments, we address the dependence of the discretization error of the
method on variations in the friction coefficient, ratio of bulk fluid
viscosity, surface fluid viscosity, mesh size, and position of the interface
relative to the fixed computational mesh.

One option to solve the coupled bulk-surface problem is to adopt a monolithic
approach. However, such an approach would quickly lead to high computational
costs as the mesh gets refined. Instead, we introduce a partitioned scheme
based on fixed-point iterations [44]. With this algorithm, the bulk and
surface flow problems are solved separately and sequentially and the coupling
conditions are enforced in an iterative fashion. We choose this algorithm for
its simplicity of implementation and study its performance when model and
discretization parameters vary. Our numerical experiments provide evidence of
the robustness of the proposed approach with respect to the contrast in
viscosity in the bulk fluid, surface fluid viscosity, and position of the
interface relative to the background mesh. At the same time, convergence
slowdown was observed for certain values of the slip coefficients.

The outline of the paper is as follows. In Sec. 2, we introduce the strong
formulation of the coupled problem and the associated energy balance. Sec. 3
presents the strong and weak formulations of the simplified problem, together
with the finite element discretization. In Sec. 4, we propose a partitioned
algorithm for the numerical solution of the coupled problem. Numerical results
in 3 dimensions are reported in Sec. 5. For all the simulations in this paper,
we have used NGsolve [2, 19], a high performance multiphysics finite element
software with a Python interface, and add-on library ngsxfem [3].

##  2 A two-phase fluid with material viscous interface

Consider a fixed volume $\Omega\subset\mathbb{R}^{3}$ filled with two
immiscible, viscous, and incompressible fluids separated by an interface
$\Gamma(t)$ for all $t\in[0,T]$. We assume that $\Gamma(t)$ stays closed and
sufficiently smooth (at least $C^{2}$) for all $t\in[0,T]$. Surface
$\Gamma(t)$ separates $\Omega$ into two phases (subdomains) $\Omega_{+}(t)$
and $\Omega_{-}(t):=\Omega\setminus\overline{\Omega_{+}(t)}$. We assume
$\Omega_{-}(t)$ to be completely internal, i.e.
$\partial\Omega_{-}(t)\cap\partial\Omega=\emptyset$ for all times. See Fig. 1.

Fig. 1: Fluid domains and interface in $\mathbb{R}^{2}$.

Denote by $\mathbf{n}^{\pm}$ the outward normals for $\Omega_{\pm}(t)$ and
$\mathbf{n}$ the normal on $\Gamma$ pointing from $\Omega_{-}(t)$ to
$\Omega_{+}(t)$: it holds that $\mathbf{n}^{-}=\mathbf{n}$ and
$\mathbf{n}^{+}=-\mathbf{n}$ at $\Gamma$. For ease of notation, from now on we
will drop the dependance on $t$ for $\Gamma$, $\Omega_{+}$, and $\Omega_{-}$.

To describe the fluid motion, consider the bulk fluid velocities
$\mathbf{u}^{\pm}(t):\Omega_{\pm}\to\mathbb{R}^{3}$ and pressures
$p^{\pm}(t):\Omega_{\pm}\to\mathbb{R}$. The motion of the fluids occupying
subdomains $\Omega_{\pm}$ is governed by the incompressible Navier–Stokes
equations

$\displaystyle\rho^{\pm}\accentset{\mbox{\large.}}{\partial}_{t}\mathbf{u}^{\pm}$
$\displaystyle=\mathop{\rm div}\boldsymbol{\sigma}^{\pm}+\mathbf{f}^{\pm}$
$\displaystyle\text{ in }\Omega_{\pm},$ (2.1) $\displaystyle\mathop{\rm
div}{}\mathbf{u}^{\pm}$ $\displaystyle=0$ $\displaystyle\text{ in
}\Omega_{\pm},$ (2.2)

for all $t\in(0,T)$. In (2.1), constants $\rho^{\pm}$ represent the fluid
density, $\accentset{\mbox{\large.}}{\partial}_{t}$ denotes the material
derivative, $\mathbf{f}^{\pm}$ are the external body forces, and
$\boldsymbol{\sigma}^{\pm}$ are the Cauchy stress tensors. For Newtonian
fluids in both phases, the Cauchy stress tensors have 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},$ (2.3)

where constants $\mu^{\pm}$ represent the fluid dynamic viscosity.

We assume interface $\Gamma$ to be a thin material layer with possibly
different material properties from the bulk fluid. Motivated by applications
in cell biology, we consider a viscous inextensible interface modeled as an
“incompressible” surface fluid. The evolution of the material interface can be
described in terms of the velocity of this surface fluid denoted by
$\mathbf{U}$. Later, we will need the decomposition of $\mathbf{U}$ into
tangential and normal components: $\mathbf{U}=\mathbf{U}_{T}+U_{N}\mathbf{n}$,
with $\mathbf{U}_{T}\cdot\mathbf{n}=0$, $U_{N}=\mathbf{n}\cdot\mathbf{U}$. The
surface Navier-Stokes equations governing the motion of a fluidic deformable
layer appear in several work [53, 49, 28, 36]. Here, we adopt the formulation
in terms of tangential differential operators from [28], where the equations
have been derived from conservation principles. To introduce these equations,
we need some further definitions. Let
$\mathbf{P}(\mathbf{x}):=\mathbf{I}-\mathbf{n}(\mathbf{x})\mathbf{n}(\mathbf{x})^{T}$
for $\mathbf{x}\in\Gamma$ be the orthogonal projection onto the tangent plane.
For a scalar function $\pi:\,\Gamma\to\mathbb{R}$ or a vector function
$\mathbf{U}:\,\Gamma\to\mathbb{R}^{3}$ we define
$\pi^{e}\,:\,\mathcal{O}(\Gamma)\to\mathbb{R}$,
$\mathbf{U}^{e}\,:\,\mathcal{O}(\Gamma)\to\mathbb{R}^{3}$, smooth extensions
of $\pi$ and $\mathbf{U}$ from $\Gamma$ to its neighborhood
$\mathcal{O}(\Gamma)$. The surface gradient and covariant derivatives on
$\Gamma$ are then defined as $\nabla_{\Gamma}\pi=\mathbf{P}\nabla\pi^{e}$ and
$\nabla_{\Gamma}\mathbf{U}:=\mathbf{P}\nabla\mathbf{U}^{e}\mathbf{P}$. The
definitions of $\nabla_{\Gamma}\pi$ and $\nabla_{\Gamma}\mathbf{U}$ are
independent of the particular smooth extension of $\pi$ and $\mathbf{U}$ off
$\Gamma$. On $\Gamma$, we consider the surface rate-of-strain tensor [22]
given by

$\mathbf{D}_{\Gamma}(\mathbf{U}):=\frac{1}{2}\mathbf{P}(\nabla\mathbf{U}+(\nabla\mathbf{U})^{T})\mathbf{P}=\frac{1}{2}(\nabla_{\Gamma}\mathbf{U}+(\nabla_{\Gamma}\mathbf{U})^{T}).$
(2.4)

The surface divergence for a vector $\mathbf{g}:\Gamma\to\mathbb{R}^{3}$ and a
tensor $\mathbf{A}:\Gamma\to\mathbb{R}^{3\times 3}$ are defined as:
${\mathop{\,\rm div}}_{\Gamma}\mathbf{g}:={\rm
tr}(\nabla_{\Gamma}\mathbf{g})$, ${\mathop{\,\rm
div}}_{\Gamma}\mathbf{A}:=\mbox{\small$\left({\mathop{\,\rm
div}}_{\Gamma}(\mathbf{e}_{1}^{T}\mathbf{A}),\,{\mathop{\,\rm
div}}_{\Gamma}(\mathbf{e}_{2}^{T}\mathbf{A}),\,{\mathop{\,\rm
div}}_{\Gamma}(\mathbf{e}_{3}^{T}\mathbf{A})\right)^{T}$ },$where
$\mathbf{e}_{i}$ is the $i$th standard basis vector. With this notation, the
surface Navier–Stokes equations take the form:

$\displaystyle\rho_{\Gamma}\accentset{\mbox{\large.}}{\partial}_{t}\mathbf{U}$
$\displaystyle=-\nabla_{\Gamma}\pi+2\mu_{\Gamma}{\mathop{\,\rm
div}}_{\Gamma}\mathbf{D}_{\Gamma}(\mathbf{U})+\mathbf{f}_{\Gamma}+\mathbf{b}^{e}+\pi\kappa\mathbf{n}$
$\displaystyle\text{ on }\Gamma,$ (2.5) $\displaystyle{\mathop{\,\rm
div}}_{\Gamma}\mathbf{U}$ $\displaystyle=0$ $\displaystyle\text{ on }\Gamma,$
(2.6)

where $\rho_{\Gamma}$ is the surface fluid density, $\mu_{\Gamma}$ is the
surface fluid dynamic viscosity, $\kappa$ denotes point-wise doubled mean
curvature on $\Gamma$, and $\pi$ is the surface fluid pressure. The material
derivative in (2.5) is taken with respect to surface fluid trajectories, i.e.
$\accentset{\mbox{\large.}}{\partial}_{t}\mathbf{U}=\frac{\partial\mathbf{U}}{\partial
t}+(\mathbf{U}\cdot\nabla)\mathbf{U}$. Note that
$\accentset{\mbox{\large.}}{\partial}_{t}\mathbf{U}$ is an intrinsic surface
quantity, although both terms $\frac{\partial\mathbf{U}}{\partial t}$ and
$(\mathbf{U}\cdot\nabla)\mathbf{U}$ depend on extension of $\mathbf{U}$ in the
bulk. On the right hand side of (2.5), $\mathbf{f}_{\Gamma}$ denotes the
external area force acting on the surface as a result of the interaction with
the bulk fluids (specified below), while $\mathbf{b}^{e}$ denotes other
possible area force (such as elastic bending forces) and not specified further
in the present paper.

Next, we turn to the coupling conditions between equations (2.1)–(2.2) posed
in the bulk and equations (2.5)–(2.6) posed on $\Gamma$. First, the
immiscibility condition means that the bulk fluid does not penetrate through
$\Gamma$, which implies that

$\mathbf{u}^{+}\cdot\mathbf{n}=U_{N}=\mathbf{u}^{-}\cdot\mathbf{n}\qquad\qquad\text{
on }\Gamma.$ (2.7)

Normal velocity $U_{N}$ determines radial deformations of $\Gamma(t)$ and so
it governs the geometric evolution of the interface, which can be defined
through the Lagrangian mapping $\Psi(t,\cdot)$ from $\Gamma(0)$ to
$\Gamma(t)$: for $\mathbf{x}\in\Gamma(0)$, $\Psi(t,\mathbf{x})$ solves the ODE
system

$\Psi(0,\mathbf{x})=\mathbf{x},\quad\frac{\partial\Psi(t,\mathbf{x})}{\partial
t}=U_{N}(t,\Psi(t,\mathbf{x})),\quad t\in[0,T].$ (2.8)

In fluid vesicles and cells, typically a viscous and dense lipid membrane,
represented by $\Gamma$, is surrounded by a less viscous and less dense
liquid. We are interested in modelling slip with friction between the bulk
fluid and the viscous membrane. Thus, we consider Navier-type conditions

$\displaystyle\mathbf{P}{\boldsymbol{\sigma}^{+}\mathbf{n}}$
$\displaystyle=f^{+}(\mathbf{P}\mathbf{u}^{+}-\mathbf{U}_{T})$
$\displaystyle\text{ on }\Gamma,$ (2.9)
$\displaystyle\mathbf{P}{\boldsymbol{\sigma}^{-}\mathbf{n}}$
$\displaystyle=-f^{-}(\mathbf{P}\mathbf{u}^{-}-\mathbf{U}_{T})$
$\displaystyle\text{ on }\Gamma,$ (2.10)

where $f^{-}$ and $f^{+}$ are friction coefficients at $\Gamma$ on the
$\Omega_{-}$ and $\Omega_{+}$ side, respectively. Conditions (2.9)–(2.10)
model an incomplete adhesion of a bulk fluid to the material surface with
$1/f^{\pm}$ often referred to as a “slip length” [navier1823memoire]; see,
e.g., [6, 31] for the modern description of experimental and theoretical
validations. In particular, the acceptance of non-zero slip length resolves
the well-known “no-collision paradox” [15, 27, 1], thus suggesting
(2.9)–(2.10) to be an important modeling assumption in the simulation of a
lipid vesicle – cell membrane contact (and fusion). We finally note that the
Navier conditions are not an uncommon choice in numerical models, if the flow
in boundary region is under-resolved [30].

The area force in (2.5) coming from the bulk fluid is defined by the jump of
the normal stress on $\Gamma$:

$\mathbf{f}_{\Gamma}=[\boldsymbol{\sigma}\mathbf{n}]^{+}_{-}=\boldsymbol{\sigma}^{+}\mathbf{n}-\boldsymbol{\sigma}^{-}\mathbf{n}~{}\text{
on }\Gamma.$ (2.11)

We summarize the complete system of equations and coupling conditions below

$\left\\{\begin{aligned}
\rho^{\pm}\accentset{\mbox{\large.}}{\partial}_{t}\mathbf{u}^{\pm}-\mu^{\pm}\Delta\mathbf{u}^{\pm}+\nabla
p^{\pm}&=\mathbf{f}^{\pm}&\text{ in }\Omega_{\pm},\\\ \mathop{\rm
div}{}\mathbf{u}^{\pm}&=0&\text{ in }\Omega_{\pm},\\\
\rho_{\Gamma}\accentset{\mbox{\large.}}{\partial}_{t}\mathbf{U}-2\mu_{\Gamma}{\mathop{\,\rm
div}}_{\Gamma}\mathbf{D}_{\Gamma}(\mathbf{U})+\nabla_{\Gamma}\pi-\pi\kappa\mathbf{n}&=[\boldsymbol{\sigma}\mathbf{n}]^{+}_{-}+\mathbf{b}^{e}&\text{
on }\Gamma,\\\ {\mathop{\,\rm div}}_{\Gamma}\mathbf{U}&=0&\text{ on
}\Gamma,\\\
\mathbf{u}^{+}\cdot\mathbf{n}=\mathbf{u}^{-}\cdot\mathbf{n}&=U_{N}&\text{ on
}\Gamma,\\\ \mathbf{P}{\boldsymbol{\sigma}^{\pm}\mathbf{n}}&=\pm
f^{\pm}(\mathbf{P}\mathbf{u}^{\pm}-\mathbf{U}_{T})&\text{ on
}\Gamma,\end{aligned}\right.$ (2.12)

with the evolution of $\Omega_{\pm}$ and $\Gamma$ defined by the velocity
solving the system through (2.8). On $\partial\Omega$ the system is endowed
with boundary conditions either for the bulk velocity or for the bulk normal
stress:

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

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. At
$t=0$, initial velocity is given $\mathbf{u}^{\pm}=\mathbf{u}_{0}^{\pm}$ in
$\Omega_{\pm}(0)$ and $\mathbf{U}=\mathbf{U}_{0}$ on $\Gamma(0)$.

###  2.1 Energy balance of the continuous coupled problem

We look for the energy balance of the coupled system (2.12)–(2.14). We make
use of the following identities for time-dependent domains $\Omega_{\pm}(t)$,
whose only moving part of the boundary is $\Gamma(t)$:

$\displaystyle\frac{1}{2}\frac{d}{dt}\int_{\Omega_{\pm}(t)}|\mathbf{u}^{\pm}|^{2}\,dV$
$\displaystyle=\int_{\Omega_{\pm}(t)}\mathbf{u}^{\pm}\cdot{}\frac{\partial\mathbf{u}^{\pm}}{\partial
t}\,dV+\frac{1}{2}\int_{\Gamma(t)}|\mathbf{u}^{\pm}|^{2}\mathbf{u}^{\pm}\cdot\mathbf{n}^{\pm}\,dS,$
(2.15)
$\displaystyle\int_{\Omega_{\pm}(t)}(\mathbf{u}^{\pm}\cdot\nabla)\mathbf{u}^{\pm}\cdot\mathbf{u}^{\pm}\,dV$
$\displaystyle=\frac{1}{2}\int_{\partial\Omega_{\pm}(t)}|\mathbf{u}^{\pm}|^{2}\mathbf{u}^{\pm}\cdot\mathbf{n}^{\pm}\,dS.$
(2.16)

Identity (2.15) is the Reynolds transport theorem, while identity (2.16) is
obtained from integration by parts.

Let us start from the kinetic energy of the fluid in $\Omega_{-}$:

$\displaystyle\frac{d}{dt}E^{-}$
$\displaystyle=\frac{1}{2}\frac{d}{dt}\int_{\Omega_{-}}\\!\\!\rho^{-}|\mathbf{u}^{-}|^{2}\,dV=\int_{\Omega_{-}}\rho^{-}\mathbf{u}^{-}\cdot{}\frac{\partial\mathbf{u}^{-}}{\partial
t}\,dV+\frac{1}{2}\int_{\Gamma}\rho^{-}|\mathbf{u}^{-}|^{2}\mathbf{u}^{-}\cdot\mathbf{n}\,dS$
$\displaystyle=\int_{\Omega_{-}}\mathbf{u}^{-}\cdot(\mathop{\rm
div}\boldsymbol{\sigma}^{-}+\mathbf{f}^{-}-\rho^{-}(\mathbf{u}^{-}\cdot\nabla)\mathbf{u}^{-})\,dV+\frac{1}{2}\int_{\Gamma}\rho^{-}|\mathbf{u}^{-}|^{2}\mathbf{u}^{-}\cdot\mathbf{n}\,dS$
$\displaystyle=\int_{\Omega_{-}}\nabla\cdot(\boldsymbol{\sigma}^{-}\mathbf{u}^{-})\,dV-\int_{\Omega_{-}}\boldsymbol{\sigma}^{-}:\nabla\mathbf{u}^{-}\,dV+\int_{\Omega_{-}}\mathbf{u}^{-}\cdot\mathbf{f}^{-}\,dV$
$\displaystyle=\int_{\Gamma}\mathbf{u}^{-}\cdot(\boldsymbol{\sigma}{}^{-}\mathbf{n})\,dS-\int_{\Omega_{-}}\boldsymbol{\sigma}^{-}:\mathbf{D}(\mathbf{u}^{-})dV+\int_{\Omega_{-}}\mathbf{u}^{-}\cdot\mathbf{f}^{-}\,dV$
$\displaystyle=\int_{\Gamma}\mathbf{u}^{-}\cdot(\boldsymbol{\sigma}{}^{-}\mathbf{n})\,dS-2\mu^{-}\int_{\Omega_{-}}\|\mathbf{D}(\mathbf{u}^{-})\|^{2}\,dV+\int_{\Omega_{-}}\mathbf{u}^{-}\cdot\mathbf{f}^{-}\,dV$
(2.17)

Above, we have used (2.15), (2.1), (2.2), (2.16), and integration by parts.

We repeat similar steps for the kinetic energy of the fluid in $\Omega_{+}$,
the main difference being that
$\overline{\partial\Omega_{+}}=\overline{\Gamma}\cup\overline{\partial\Omega_{D}}\cup\overline{\partial\Omega_{N}}$
while $\overline{\partial\Omega_{-}}=\overline{\Gamma}$. We obtain:

$\displaystyle\frac{d}{dt}E^{+}=-\int_{\Gamma}\mathbf{u}^{+}\cdot{}(\boldsymbol{\sigma}^{+}\mathbf{n})\,dS-2\mu^{+}\int_{\Omega_{+}}\left\lVert\mathbf{D}(\mathbf{u}^{+})\right\rVert^{2}dV+\int_{\Omega_{+}}\mathbf{u}^{+}\cdot\mathbf{f}^{+}\,dV+B$
(2.18)

where

$\displaystyle
B=\int_{\partial\Omega_{D}}\left(\mathbf{g}\cdot(\boldsymbol{\sigma}^{+}\mathbf{n})-\frac{1}{2}\rho^{+}|\mathbf{g}|^{2}\mathbf{g}\cdot\mathbf{n}^{+}\right)\,dS+\int_{\partial\Omega_{N}}\left(\mathbf{u}^{+}\cdot{}\mathbf{f}{}-\frac{1}{2}\rho^{+}|\mathbf{u}^{+}|^{2}\mathbf{u}^{+}\cdot\mathbf{n}^{+}\right)\,dS.$
(2.19)

Putting together (2.17) and (2.18), we obtain the total kinetic energy for the
bulk flow:

$\displaystyle\frac{dE}{dt}=$
$\displaystyle\frac{d}{dt}(E^{+}+E^{-})=-2\mu^{-}\int_{\Omega_{-}}\|\mathbf{D}(\mathbf{u}^{-})\|^{2}\,dV-2\mu^{+}\int_{\Omega_{+}}\left\lVert\mathbf{D}(\mathbf{u}^{+})\right\rVert^{2}dV-\int_{\Gamma}U_{N}[\mathbf{n}^{T}\sigma\mathbf{n}]^{+}_{-}dS$
$\displaystyle-\int_{\Gamma}f^{-}\mathbf{P}\mathbf{u}^{-}\cdot(\mathbf{P}\mathbf{u}^{-}-\mathbf{U}_{T})dS-\int_{\Gamma}f^{+}\mathbf{P}\mathbf{u}^{+}\cdot(\mathbf{P}\mathbf{u}^{+}-\mathbf{U}_{T})dS$
$\displaystyle+\int_{\Omega^{-}}\mathbf{u}^{-}\cdot\mathbf{f}^{-}\,dV+\int_{\Omega^{+}}\mathbf{u}^{+}\cdot\mathbf{f}^{+}\,dV+B$
(2.20)

where we have used (2.7), (2.9), and (2.10). For the kinetic energy of the
surface flow, we will use the surface analogue of the Reynolds transport
theorem (see, e.g. [16, Lemma 2.1]):

$\displaystyle\frac{{\rm d}}{{\rm
d}t}\int_{\Gamma}f\,dS=\int_{\Gamma}\left(\accentset{\mbox{\large.}}{\partial}_{t}f+f{\mathop{\,\rm
div}}_{\Gamma}\mathbf{U}\right)\,dS.$ (2.21)

We also need the surface integration by part identity:

$\int_{\Gamma}(g{\mathop{\,\rm
div}}_{\Gamma}\mathbf{f}+\mathbf{f}\cdot\nabla_{\Gamma}g)\,dS=\int_{\Gamma}\kappa
g(\mathbf{f}\cdot\mathbf{n})\,dS,$

which is valid for a smooth closed surface with smooth scalar field $g$ and
vector field $\mathbf{f}$. A componentwise application of this equality yields
the identity

$\int_{\Gamma}\mathbf{g}\cdot{\mathop{\,\rm
div}}_{\Gamma}\mathbf{G}\,dS=-\int_{\Gamma}\mathbf{G}:\nabla_{\Gamma}\mathbf{g}\,dS$
(2.22)

for a vector field $\mathbf{g}\in$ $({C}^{1}(\Gamma))^{3}$ and a matrix
function $\mathbf{G}\in({C}^{1}(\Gamma))^{3\times{}3}$ such that
$\mathbf{G}=\mathbf{P}\mathbf{G}\mathbf{P}$.

The energy balance on $\Gamma$ is given by:

$\displaystyle\frac{dE_{\Gamma}}{dt}$
$\displaystyle=\frac{1}{2}\frac{d}{dt}\int_{\Gamma}\\!\\!\rho_{\Gamma}|\mathbf{U}|^{2}\,dS=\int_{\Gamma}\left(\rho_{\Gamma}\mathbf{U}\cdot\accentset{\mbox{\large.}}{\partial}_{t}\mathbf{U}+\frac{1}{2}\rho_{\Gamma}|\mathbf{U}|^{2}{\mathop{\,\rm
div}}_{\Gamma}\mathbf{U}\right)\,dS$
$\displaystyle=\int_{\Gamma}\mathbf{U}\cdot\left(-\nabla_{\Gamma}\pi+2\mu_{\Gamma}{\mathop{\,\rm
div}}_{\Gamma}\mathbf{D}_{\Gamma}(\mathbf{U})+\mathbf{f}_{\Gamma}+\mathbf{b}^{e}\right)\,dS$
$\displaystyle=-2\mu_{\Gamma}\int_{\Gamma}\mathbf{D}_{\Gamma}(\mathbf{U}):\nabla_{\Gamma}\mathbf{U}\,dS+\int_{\Gamma}\mathbf{U}\cdot[\boldsymbol{\sigma}\mathbf{n}]^{+}_{-}\,dS+\int_{\Gamma}\mathbf{U}_{T}\cdot\mathbf{b}^{e}\,dS$
$\displaystyle=-2\mu_{\Gamma}\int_{\Gamma}\|\mathbf{D}_{\Gamma}(\mathbf{U})\|^{2}\,dS+\int_{\Gamma}U_{N}[\mathbf{n}^{T}\boldsymbol{\sigma}\mathbf{n}]^{+}_{-}\,dS$
$\displaystyle\quad+\int_{\Gamma}\mathbf{U}_{T}\cdot\left(f^{+}(\mathbf{P}\mathbf{u}^{+}-\mathbf{U}_{T})+f^{-}(\mathbf{P}\mathbf{u}^{-}-\mathbf{U}_{T})\right)\,dS+\int_{\Gamma}\mathbf{U}\cdot\mathbf{b}^{e}\,dS,$
(2.23)

where we have applied (2.22) and used (2.21), (2.5)-(2.7), (2.9)-(2.11), and
integration by parts.

Combining (2.20) and (2.23), we get the kinetic energy for the bulk and
surface flows:

$\displaystyle\frac{dE}{dt}+\frac{dE_{\Gamma}}{dt}$
$\displaystyle=\underbrace{-2\mu^{-}\int_{\Omega_{-}}\|\mathbf{D}(\mathbf{u}^{-})\|^{2}\,dV-2\mu^{+}\int_{\Omega_{+}}\left\lVert\mathbf{D}(\mathbf{u}^{+})\right\rVert^{2}dV}_{\text{Bulk
fluid viscous
dissipation}}\underbrace{-2\mu_{\Gamma}\int_{\Gamma}\|\mathbf{D}_{\Gamma}(\mathbf{U}_{T})\|^{2}\,dS}_{\text{Surface
fluid viscous dissipation}}$
$\displaystyle\quad\underbrace{-\int_{\Gamma}f^{-}\|\mathbf{P}\mathbf{u}^{-}-\mathbf{U}_{T}\|^{2}dS-\int_{\Gamma}f^{+}\|\mathbf{P}\mathbf{u}^{+}-\mathbf{U}_{T}\|^{2}d}_{\text{Frictional
energy dissipation}}$
$\displaystyle\quad\underbrace{+\int_{\Gamma}\mathbf{u}^{-}\cdot\mathbf{f}^{-}\,dV+\int_{\Gamma}\mathbf{u}^{+}\cdot\mathbf{f}^{+}\,dV+\int_{\Gamma}\mathbf{U}\cdot\mathbf{b}^{e}dS}_{\text{work
of external forces }}\underbrace{\phantom{\int_{\Gamma}}+B}_{\text{work of
b.c.}}$ (2.24)

where $B$, i.e. the work of the boundary conditions, is defined in (2.19).

From (2.24), we see that in the absence of external forces and with no energy
inflow through the boundary the system is dissipative, i.e. thermodynamically
consistent.

##  3 A simplified steady problem

In this section, we consider a (strongly) simplified version of the problem
presented in Sec. 2. Our main assumption is that the coupled bulk and surface
fluid system has reached a steady state and inertia terms can be neglected.
Since the steady state implies $\Gamma(t)=\Gamma(0)$, we have $U_{N}=0$ and
hence $\mathbf{U}=\mathbf{U}_{T}$. This simplified surface–bulk Stokes problem
models a viscosity dominated two-phase flow with the viscous interface in a
dynamical equilibrium; see also Remark 1. It is an interesting model problem
for the purpose of numerical analysis. With these simplifications, the
equations (2.1)–(2.2) become:

$\displaystyle-\mu^{\pm}\Delta\mathbf{u}^{\pm}+\nabla p$
$\displaystyle=\mathbf{f}^{\pm}\qquad\text{in }\Omega_{\pm},$ (3.1)
$\displaystyle\mathop{\rm div}{}\mathbf{u}^{\pm}$
$\displaystyle=0\qquad~{}~{}\text{in }\Omega_{\pm}.$ (3.2)

We impose a non-homogeneous Dirichlet condition on the entire outer boundary
of $\Omega$, i.e. problem (3.1)–(3.2) is supplemented with boundary condition
(2.13) with $\mathbf{g}\in[H^{1/2}(\partial\Omega_{D})]^{3}$ and
$\partial\Omega_{D}=\partial\Omega$. Under our assumption, the momentum
equation for the surface fluid simplifies to $-2\mu_{\Gamma}{\mathop{\,\rm
div}}_{\Gamma}\mathbf{D}_{\Gamma}(\mathbf{U}_{T})+\nabla_{\Gamma}\pi-\pi\kappa\mathbf{n}=[\boldsymbol{\sigma}\mathbf{n}]^{+}_{-}+\mathbf{b}^{e}.$
The tangential part of the above momentum equation together with the
inextensibility condition (2.6) leads to the surface Stokes problem

$\displaystyle-2\mu_{\Gamma}\mathbf{P}{\mathop{\,\rm
div}}_{\Gamma}\mathbf{D}_{\Gamma}(\mathbf{U}_{T})+\nabla_{\Gamma}\pi$
$\displaystyle=[\mathbf{P}\boldsymbol{\sigma}\mathbf{n}]^{+}_{-}+\mathbf{P}\mathbf{b}^{e}\qquad$
$\displaystyle\text{ on }\Gamma,$ (3.3) $\displaystyle{\mathop{\,\rm
div}}_{\Gamma}\mathbf{U}_{T}$ $\displaystyle=0\qquad$ $\displaystyle\text{ on
}\Gamma,$ (3.4)

while the normal part simplifies to

$[\mathbf{n}^{T}\boldsymbol{\sigma}\mathbf{n}]^{-}_{+}=\pi\kappa\qquad\text{
on }\Gamma.$ (3.5)

The interface condition above is standard in many models of two-phase flows,
where $\pi$ has the meaning of the surface tension coefficient.

Coupling condition (2.7) is replaced by:

$\mathbf{u}^{+}\cdot\mathbf{n}=\mathbf{u}^{-}\cdot\mathbf{n}\qquad\text{ on
}\Gamma,$ (3.6)

while conditions (2.9) and (2.10) still hold:

$\mathbf{P}{\boldsymbol{\sigma}^{\pm}\mathbf{n}}=\pm
f^{\pm}(\mathbf{P}\mathbf{u}^{\pm}-\mathbf{U}_{T})\qquad\text{ on }\Gamma.$
(3.7)

Finally, we will see that the (weak formulation of the) problem is well-posed
under two mean conditions for the bulk pressure:

$\int_{\Omega^{\pm}}p^{\pm}\,dx=0.$

######  Remark 1.

Since $\mathbf{U}\cdot\mathbf{n}=0$, condition (3.6) allows the flow through
the steady interface $\Gamma$. This is inconsistent with (2.7), which assumes
immiscibility of fluids. For a physically consistent formulation that
describes the true equilibrium one has to set
$\mathbf{u}^{+}\cdot\mathbf{n}=\mathbf{u}^{-}\cdot\mathbf{n}=0$ on $\Gamma$,
but allow the shape of $\Gamma$ to be the unknown, i.e. to be determined as a
part of the problem. For such equilibrium to exist, external forces and
boundary conditions may have to satisfy additional constraints. Finding such
constraints and solving the resulting non-linear problem is outside the scope
of this paper. We rather follow a common convention in the analysis of models
for steady two-phase problems and allow (3.6) for the steady interface; see,
e.g. [18, 42, 21, 8].

###  3.1 Variational formulation

The purpose of this section is to derive the variational formulation of
coupled problem (3.1)–(3.7). Let us introduce some standard notation. The
space of functions whose square is integrable in a domain $\omega$ is denoted
by $L^{2}(\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)^{3}$. $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_{-})^{3},~{}V^{+}=\\{\mathbf{u}\in
H^{1}(\Omega_{+})^{3},\mathbf{u}{\big{|}_{\partial\Omega_{D}}}=\mathbf{g}\\},~{}V^{+}_{0}=\\{\mathbf{u}\in
H^{1}(\Omega_{+})^{3},\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 L^{2}_{\pm}$ $\displaystyle=\\{p=(p^{-},p^{+})\in
L^{2}(\Omega_{-})\times
L^{2}(\Omega_{+}),~{}s.t.~{}\int_{\Omega^{\pm}}p^{\pm}\,dx=0\\},$
$\displaystyle V_{\Gamma}$ $\displaystyle=\\{\mathbf{U}\in
H^{1}(\Gamma)^{3}:\mathbf{U}\cdot\mathbf{n}=0\\}.$

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

Multiplying (3.1) by $\mathbf{v}\in V^{\pm}_{0}$ and (3.2) by $q\in
L^{2}_{0}(\Omega)$ and integrating over each subdomain, we see that smooth
bulk velocity and pressure satisfy integral identity:

$\displaystyle-(p^{-},\mathop{\rm
div}\mathbf{v}^{-})_{\Omega_{-}}-(p^{+},\mathop{\rm
div}\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\pi\kappa,\mathbf{v}^{-}\cdot\mathbf{n}\rangle_{\Gamma}+\langle
f^{-}(\mathbf{P}\mathbf{u}^{-}-\mathbf{U}),\mathbf{P}\mathbf{v}^{-}\rangle_{\Gamma}+\langle
f^{+}(\mathbf{P}\mathbf{u}^{+}-\mathbf{U}),\mathbf{P}\mathbf{v}^{+}\rangle_{\Gamma}$
$\displaystyle\quad\quad=(\mathbf{f}^{-},\mathbf{v}^{-})_{\Omega_{-}}+(\mathbf{f}^{+},\mathbf{v}^{+})_{\Omega_{+}}$
(3.8) $\displaystyle(\mathop{\rm
div}\mathbf{u}^{-},q^{-})_{\Omega_{-}}+(\mathop{\rm
div}\mathbf{u}^{+},q^{+})_{\Omega_{+}}=0$ (3.9)

for all $(\mathbf{v},q)\in V^{\pm}_{0}\times L^{2}_{0}(\Omega)$. The interface
terms in (3.8) have been obtained using coupling conditions (3.7) and (3.5) 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{U}),\mathbf{P}\mathbf{v}^{-}\rangle_{\Gamma}+\langle
f^{+}(\mathbf{P}\mathbf{u}^{+}-\mathbf{U}),\mathbf{P}\mathbf{v}^{+}\rangle_{\Gamma}$
$\displaystyle\quad-\langle\pi\kappa,\mathbf{v}^{-}\cdot\mathbf{n}\rangle_{\Gamma}.$

Likewise, we find that the surface velocity and pressure satisfy the following
integral identities:

$\displaystyle-(\pi,{\mathop{\,\rm
div}}_{\Gamma}\mathbf{V})_{\Gamma}+2(\mu_{\Gamma}\mathbf{D}_{\Gamma}(\mathbf{U}),\mathbf{D}_{\Gamma}(\mathbf{V}))_{\Gamma}-\langle
f^{-}(\mathbf{P}\mathbf{u}^{-}-\mathbf{U}),\mathbf{V}\rangle_{\Gamma}$
$\displaystyle\quad\quad-\langle
f^{+}(\mathbf{P}\mathbf{u}^{+}-\mathbf{U}),\mathbf{V}\rangle_{\Gamma}=(\mathbf{P}\mathbf{b}^{e},\mathbf{V})_{\Gamma}$
(3.10) $\displaystyle({\mathop{\,\rm
div}}_{\Gamma}\mathbf{U},\tau)_{\Gamma}=0$ (3.11)

for all $(\mathbf{V},\tau)\in V_{\Gamma}\times L^{2}_{0}(\Gamma)$.

The weak formulation of the coupled problem (3.1)–(3.7) follows by combining
(3.8)–(3.9) and (3.10)–(3.11). In order to write it, we introduce the
following forms for all $\mathbf{u}\in V^{\pm}$,
$\mathbf{v}\in\mathbf{V}^{\pm}_{0}$, $\mathbf{U},\mathbf{V}\in V_{\Gamma}$,
$p\in L^{2}(\Omega)$, $\pi\in L^{2}(\Gamma)$:

$\displaystyle a(\\{\mathbf{u},\mathbf{U}\\},\\{\mathbf{v},\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_{+}}$
$\displaystyle+2(\mu_{\Gamma}\mathbf{D}_{\Gamma}(\mathbf{U}),\mathbf{D}_{\Gamma}(\mathbf{V}))_{\Gamma}+\langle
f^{-}(\mathbf{P}\mathbf{u}^{-}-\mathbf{U}),\mathbf{P}\mathbf{v}^{-}-\mathbf{V}\rangle_{\Gamma}$
$\displaystyle+\langle
f^{+}(\mathbf{P}\mathbf{u}^{+}-\mathbf{U}),\mathbf{P}\mathbf{v}^{+}-\mathbf{V}\rangle_{\Gamma},$
$\displaystyle b(\\{\mathbf{v},\mathbf{V}\\},\\{p,\pi\\})=$
$\displaystyle-(p^{-},\mathop{\rm
div}\mathbf{v}^{-})_{\Omega_{-}}-(p^{+},\mathop{\rm
div}\mathbf{v}^{+})_{\Omega_{+}}-(\pi,{\mathop{\,\rm
div}}_{\Gamma}\mathbf{V})_{\Gamma},$ $\displaystyle s(\mathbf{v},\pi)=$
$\displaystyle-\langle\pi\kappa,\mathbf{v}\cdot\mathbf{n}\rangle_{\Gamma},$
$\displaystyle r(\mathbf{v},\mathbf{V})=$
$\displaystyle(\mathbf{f}^{-},\mathbf{v}^{-})_{\Omega_{-}}+(\mathbf{f}^{+},\mathbf{v}^{+})_{\Omega_{+}}+(\mathbf{P}\mathbf{b}^{e},\mathbf{V})_{\Gamma}.$

Then, the weak formulation reads: Find $(\mathbf{u},p)\in V^{\pm}\times
L^{2}_{\pm}$, and $(\mathbf{U},\pi)\in V_{\Gamma}\times L^{2}(\Gamma)$ such
that

$\left\\{\begin{split}a(\\{\mathbf{u},\mathbf{U}\\},\\{\mathbf{v},\mathbf{V}\\})+b(\\{\mathbf{v},\mathbf{V}\\},\\{p,\pi\\})+s(\mathbf{v},\pi)&=r(\mathbf{v},\mathbf{V})\\\
b(\\{\mathbf{u},\mathbf{U}\\},\\{q,\tau\\})&=0\end{split}\right.$ (3.12)

for all $(\mathbf{v},q)\in V_{0}^{\pm}\times L^{2}_{0}(\Omega)$ and
$(\mathbf{V},\tau)\in V_{\Gamma}\times L^{2}_{0}(\Gamma)$. Note that test and
trial pressure spaces both involve two (different) gauge conditions.

###  3.2 Well-Posedness

With the goal of proving the well-posedness of the stationary problem, we
start by showing that $a(\\{\cdot,\cdot\\},\\{\cdot,\cdot\\})$ is coercive.
Let
$\left\lVert\mathbf{u}\right\rVert^{2}_{H^{1}(\Omega_{\pm})}=\left\lVert\mathbf{u}^{+}\right\rVert^{2}_{H^{1}(\Omega_{+})}+\left\lVert\mathbf{u}^{-}\right\rVert^{2}_{H^{1}(\Omega_{-})}$.
Let $\left\lVert p\right\rVert_{L^{2}(\Omega_{\pm})}^{2}=\left\lVert
p\right\rVert_{\Omega_{+}}^{2}+\left\lVert p\right\rVert_{\Omega_{-}}^{2}$. We
define the following additional norms:

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{\mathbf{v},\mathbf{V}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}=\left\lVert\mathbf{v}\right\rVert_{H^{1}(\Omega_{\pm})}^{2}+\left\lVert\mathbf{V}\right\rVert_{H^{1}(\Gamma)}^{2},\quad{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}=\left\lVert
p\right\rVert^{2}_{\Omega_{\pm}}+\left\lVert\pi\right\rVert^{2}_{\Gamma}$

The coercivity result is formulated in the form of a lemma.

######  Lemma 2.

For any $\mathbf{v}\in V^{\pm}_{0}$ and $\mathbf{V}\in V_{\Gamma}$ it holds

$a(\\{\mathbf{v},\mathbf{V}\\},\\{\mathbf{v},\mathbf{V}\\})\geq
C{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{\mathbf{v},\mathbf{V}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}$
(3.13)

with a positive constant $C$, which may depend on the viscosity values and
$\Omega_{\pm}$.

###### Proof.

One readily computes that

$a(\\{\mathbf{v},\mathbf{V}\\},\\{\mathbf{v},\mathbf{V}\\})=2(\mu^{-}\mathbf{D}(\mathbf{v}^{-}),\mathbf{D}(\mathbf{v}^{-}))_{\Omega_{-}}+2(\mu^{+}\mathbf{D}(\mathbf{v}^{+}),\mathbf{D}(\mathbf{v}^{+}))_{\Omega_{+}}\\\
+2(\mu_{\Gamma}\mathbf{D}_{\Gamma}(\mathbf{V}),\mathbf{D}_{\Gamma}(\mathbf{V}))_{\Gamma}+f^{-}\|\mathbf{P}\mathbf{v}^{-}-\mathbf{V}\|_{\Gamma}^{2}+f^{+}\|\mathbf{P}\mathbf{v}^{+}-\mathbf{V}\|_{\Gamma}^{2}.$
(3.14)

Since function $\mathbf{v}^{+}$ satisfies homogeneous Dirichlet boundary
condition on $\partial\Omega_{+}\setminus\Gamma$, we apply the following
Korn’s inequality in $\Omega_{+}$:

$\left\lVert\mathbf{v}^{+}\right\rVert_{H^{1}(\Omega_{+})}\leq
C\|\mathbf{D}(\mathbf{v}^{+})\|_{\Omega_{+}}$ (3.15)

By the triangle and trace inequalities in $\Gamma$, we get

$\left\lVert\mathbf{V}\right\rVert_{\Gamma}\leq\|\mathbf{P}\mathbf{v}^{+}-\mathbf{V}\|_{\Gamma}+\|\mathbf{P}\mathbf{v}^{+}\|_{\Gamma}\leq\|\mathbf{P}\mathbf{v}^{+}-\mathbf{V}\|_{\Gamma}+C\left\lVert\mathbf{v}^{+}\right\rVert_{H^{1}(\Omega_{+})}.$
(3.16)

We further apply Korn’s inequality on $\Gamma$ [28]:

$\left\lVert\mathbf{V}\right\rVert_{H^{1}(\Gamma)}\leq
C\left(\left\lVert\mathbf{V}\right\rVert_{\Gamma}+\|\mathbf{D}_{\Gamma}(\mathbf{V})\|_{\Gamma}\right).$
(3.17)

Next, we can estimate the trace of $\mathbf{v}^{-}$ on $\Gamma$ through the
triangle inequality:

$\left\lVert\mathbf{v}^{-}\right\rVert_{\Gamma}\leq\|\mathbf{P}\mathbf{v}^{-}-\mathbf{V}\|_{\Gamma}+\|\mathbf{V}\|_{\Gamma}\leq\|\mathbf{P}\mathbf{v}^{-}-\mathbf{V}\|_{\Gamma}+\left\lVert\mathbf{V}\right\rVert_{H^{1}(\Gamma)}.$
(3.18)

We finally apply the following Korn’s inequality in $\Omega_{-}$:

$\left\lVert\mathbf{v}\right\rVert_{H^{1}(\Omega_{-})}\leq
C\left(\|\mathbf{D}(\mathbf{v}^{-})\|_{\Omega_{-}}+\left\lVert\mathbf{v}^{-}\right\rVert_{\Gamma}\right).$
(3.19)

Identity (3.14) and inequalities (3.15)–(3.19) lead to (3.13) after easy
computations.
∎

The continuity of the bilinear forms $a(\\{\cdot,\cdot\\},\\{\cdot,\cdot\\})$,
$b(\\{\cdot,\cdot\\},\\{\cdot,\cdot\\})$ and $s(\cdot,\cdot)$ follows from
standard arguments based on the Cauchy–Schwarz and triangle inequalities:

$\begin{split}a(\\{\mathbf{u},\mathbf{U}\\},\\{\mathbf{v},\mathbf{V}\\})\leq
C{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{\mathbf{u},\mathbf{U}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{\mathbf{v},\mathbf{V}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\quad\text{for
all}~{}\mathbf{u},\mathbf{v}\in V^{\pm}_{0},~{}\mathbf{U},\mathbf{V}\in
V_{\Gamma},\\\ b(\\{\mathbf{u},\mathbf{U}\\},\\{p,\pi\\})\leq
C{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{\mathbf{u},\mathbf{U}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\quad\text{for
all}~{}\mathbf{u}\in V^{\pm}_{0},~{}\mathbf{U}\in V_{\Gamma},~{}p\in
L^{2}(\Omega),~{}\pi\in L^{2}(\Gamma),\\\ s(\mathbf{v},\pi)\leq
C{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{\mathbf{v},0\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{0,\pi\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\quad\text{for
all}~{}\mathbf{v}\in V^{\pm},~{}\pi\in L^{2}(\Gamma).\end{split}$ (3.20)

Problem (3.12) falls into the class of so-called generalized saddle point
problems. An abstract well-posedness result for such problems can be found,
e.g. in [5, 37], which extend the Babuşka–Brezzi theory. Applied to (3.12),
this well-posedness result requires coercivity (3.13), continuity (3.20) and
two inf-sup conditions formulated in the following lemma.

######  Lemma 3.

The following inf-sup conditions hold with positive constants $\gamma_{1}$ and
$\gamma_{2}$:

$\displaystyle\sup_{\mathbf{v}\in V^{\pm}_{0},\mathbf{V}\in
V_{\Gamma}}\frac{b(\\{\mathbf{v},\mathbf{V}\\},\\{p,\pi\\})+s(\mathbf{v},\pi)}{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{\mathbf{v},\mathbf{V}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}}\geq\gamma_{1}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|},\quad\forall~{}p\in
L^{2}_{\pm},~{}\pi\in L^{2}(\Gamma),$ (3.21) $\displaystyle\sup_{\mathbf{v}\in
V^{\pm}_{0},\mathbf{V}\in
V_{\Gamma}}\frac{b(\\{\mathbf{v},\mathbf{V}\\},\\{p,\pi\\})}{{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{\mathbf{v},\mathbf{V}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}}\geq\gamma_{2}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|},\quad\forall~{}p\in
L^{2}_{0}(\Omega),~{}\pi\in L^{2}_{0}(\Gamma).$ (3.22)

###### Proof.

The proof follows by combining well-known results about the existence of a
continuous right inverse of the divergence operator in $H^{1}_{0}(\Omega)^{3}$
[7] and $V_{\Gamma}$ [28]: For arbitrary $p\in L^{2}_{0}(\Omega)$ and $\pi\in
L^{2}_{0}(\Gamma)$ there exist $\mathbf{v}\in H^{1}_{0}(\Omega)^{3}$ and
$\mathbf{V}\in V_{\Gamma}$ such that

$\begin{split}p=\mathop{\rm
div}\mathbf{v}~{}\text{~{}in}~{}\Omega,\quad\text{and}\quad\|\mathbf{v}\|_{H^{1}(\Omega)}\leq
c_{\Omega}\|p\|_{L^{2}(\Omega)},\\\ \pi={\mathop{\,\rm
div}}_{\Gamma}\mathbf{V}~{}\text{~{}on}~{}\Gamma,\quad\text{and}\quad\|\mathbf{V}\|_{H^{1}(\Gamma)}\leq
c_{\Gamma}\|\pi\|_{\Gamma}.\end{split}$ (3.23)

Letting $\mathbf{v}^{\pm}=\mathbf{v}|_{\Omega_{\pm}}$,
$(\mathbf{v}^{-},\mathbf{v}^{+})^{T}\in V^{\pm}_{0}$, and adding estimates in
(3.23) we get

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}$
$\displaystyle\leq b(\\{\mathbf{v},\mathbf{V}\\},\\{p,\pi\\}),$ (3.24)
$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{\mathbf{v},\mathbf{V}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}$
$\displaystyle\leq
c_{\Omega}\|p\|_{\Omega}+c_{\Gamma}\|\pi\|_{\Gamma}\leq(c_{\Omega}+c_{\Gamma}){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}.$
(3.25)

This proves (3.22) with $\gamma_{2}=1/(c_{\Omega}+c_{\Gamma})$.

To show (3.21), we split $\pi=\pi_{0}+\pi^{\perp}$ with $\pi_{0}\in
L^{2}_{0}(\Gamma)$ and $\pi^{\perp}=|\Gamma|^{-1}\int_{\Gamma}\pi\,ds$. For
the $\pi_{0}$ part of $\pi$, we use again (3.23) as above, while for
$p^{\pm}\in L^{2}_{0}(\Omega^{\pm})$ we use the existence of a continuous
right inverse of $\mathop{\rm div}$ in $H^{1}_{0}(\Omega^{\pm})^{3}$ to claim
the existence of $\mathbf{v}\in H^{1}_{0}(\Omega^{-})^{3}\times
H^{1}_{0}(\Omega^{+})^{3}\subset V^{\pm}_{0}$ and $\mathbf{V}\in V_{\Gamma}$
such that

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi_{0}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq
b(\\{\mathbf{v},\mathbf{V}\\},\\{p,\pi_{0}\\})+s(\mathbf{v},\pi),\qquad{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{\mathbf{v},\mathbf{V}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq(c_{\Omega}+c_{\Gamma}){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi_{0}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|},$
(3.26)

with some positive $c_{\Omega},c_{\Gamma}$ depending only on $\Gamma$ and
$\Omega$. We also used that $\mathbf{v}=\boldsymbol{0}$ on $\Gamma$ implies
$s(\mathbf{v},\pi)=0$.

Let $C^{\pm}=\pm|\Omega^{\pm}|^{-1}|\Gamma|\int_{\Gamma}\kappa\,ds$. To
control $\|\pi^{\perp}\|_{\Gamma}$, we need $\mathbf{v}_{1}\in V^{\pm}_{0}$
such that

$\mathop{\rm
div}\mathbf{v}_{1}=-C^{\pm}~{}~{}\text{in}~{}\Omega^{\pm},~{}~{}\mathbf{v}_{1}\cdot\mathbf{n}=\kappa~{}~{}\text{on}~{}\Gamma~{}~{}\text{and}~{}~{}\|\mathbf{v}_{1}\|_{H^{1}(\Omega^{\pm})}\leq
C.$ (3.27)

Such $\mathbf{v}_{1}$ can be built, for example, as follows: Let
$\mathbf{v}_{1}^{-}=\nabla\psi$, where $\psi\in H^{2}(\Omega^{-})$ solves the
Neumann problem $-\Delta\psi=C^{-}$ in $\Omega^{-}$,
$\mathbf{n}\cdot\nabla\psi=\kappa$ on $\Gamma$. Since
$\Gamma=\partial\Omega^{-}$ is smooth, by the $H^{2}$-regularity of the
Neumann problem we have that
$\|\mathbf{v}_{1}^{-}\|_{H^{1}(\Omega^{-})}\leq\|\psi\|_{H^{2}(\Omega^{-})}\leq
C$. The boundary $\partial\Omega$ is only Lipschits and so the Neumann problem
in $\Omega^{+}$ is not necessarily $H^{2}$-regular. To handle this, we first
extend $\mathbf{v}_{1}^{-}$ from $\Omega^{-}$ to a function
$\tilde{\mathbf{v}}_{1}$ in $H^{1}_{0}(\Omega)^{3}$ such that
$\|\tilde{\mathbf{v}}_{1}\|_{H^{1}(\Omega^{+})}\leq
c\|\mathbf{v}_{1}^{-}\|_{H^{1}(\Omega^{-})}$ [50]. Next, we consider
$\mathbf{w}\in H^{1}_{0}(\Omega^{+})^{3}$ such that $\mathop{\rm
div}\mathbf{w}=C^{+}-\mathop{\rm div}\tilde{\mathbf{v}}_{1}\in
L^{2}_{0}(\Omega^{+})$, and $\|\mathbf{w}\|_{H^{1}(\Omega^{+})}\leq
c_{\Omega^{+}}\|\mathop{\rm div}\mathbf{w}\|_{L^{2}(\Omega^{+})}\leq C$ [7].
The desired $\mathbf{v}_{1}^{+}$ is given in $\Omega^{+}$ by
$\mathbf{v}_{1}^{+}=\tilde{\mathbf{v}}_{1}+\mathbf{w}$. Since $\mathop{\rm
div}\mathbf{v}_{1}=-C^{\pm}$, for $p\in L^{2}_{\pm}(\Omega)$ and $\pi\in
L^{2}(\Gamma)$ we have identities

$b(\\{\mathbf{v}_{1},0\\},\\{p,\pi\\})=0=({\mathop{\,\rm
div}}_{\Gamma}\mathbf{V},\pi^{\perp})_{\Gamma}.$ (3.28)

We also note the equality
$\|\pi^{\perp}\|_{\Gamma}^{2}=\hat{c}\,s(\mathbf{v}_{1},\pi^{\perp})$, with

$\hat{c}=\pi^{\perp}|\Gamma|/\int_{\Gamma}\kappa^{2}ds.$

The denominator above is positive, since $\Gamma$ is closed and so $\kappa$
cannot be zero everywhere on $\Gamma$. We use (3.24)–(3.28) to estimate for
some $\beta>0$:

$\begin{split}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}&={\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi_{0}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}+\beta\|\pi^{\perp}\|_{\Gamma}^{2}\\\
&\leq
b(\\{\mathbf{v},\mathbf{V}\\},\\{p,\pi_{0}\\})+s(\mathbf{v},\pi)+s(\beta\hat{c}\,\mathbf{v}_{1},\pi^{\perp})\\\
&=b(\\{\mathbf{v}+\beta\hat{c}\mathbf{v}_{1},\mathbf{V}\\},\\{p,\pi\\})+s(\mathbf{v}+\beta\hat{c}\mathbf{v}_{1},\pi)-(\beta\hat{c}\mathbf{v}_{1},\pi_{0})_{\Gamma}\\\
&\leq
b(\\{\mathbf{v}+\beta\hat{c}\mathbf{v}_{1},\mathbf{V}\\},\\{p,\pi\\})+s(\mathbf{v}+\beta\hat{c}\mathbf{v}_{1},\pi)+\frac{\beta^{2}\hat{c}^{2}}{2}\|\mathbf{v}_{1}\|^{2}_{\Gamma}+\frac{1}{2}\|\pi_{0}\|^{2}_{\Gamma}\\\
&\leq
b(\\{\mathbf{v}+\beta\hat{c}\mathbf{v}_{1},\mathbf{V}\\},\\{p,\pi\\})+s(\mathbf{v}+\beta\hat{c}\mathbf{v}_{1},\pi)+c_{3}\,\beta^{2}\|\pi^{\perp}\|_{\Gamma}^{2}+\frac{1}{2}\|\pi_{0}\|^{2}_{\Gamma}.\end{split}$

with some $c_{3}>0$ depending only on $\Gamma$ and $\Omega$. For $\beta>0$
sufficiently small such that $\frac{\beta}{2}-c_{3}\,\beta^{2}\geq 0$, we get

$c{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}\leq
b(\\{\mathbf{v}+\beta\hat{c}\mathbf{v}_{1},\mathbf{V}\\},\\{p,\pi\\})+s(\mathbf{v}+\beta\hat{c}\mathbf{v}_{1},\pi),$
(3.29)

with $c>0$ depending only on $\Gamma$ and $\Omega$. Thanks to the triangle
inequality, the second estimate in (3.26) and the definition of $\hat{c}$ and
$\mathbf{v}_{1}$, we find the bound

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v}+\beta\hat{c}\mathbf{v}_{1},\mathbf{V}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\mathbf{v},\mathbf{V}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+\beta\hat{c}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\beta\hat{c}\mathbf{v}_{1},0\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}\leq(c_{\Omega}+c_{\Gamma}){\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi_{0}\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}+C\,\|\pi^{\perp}\|_{\Gamma}\leq
C\,{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|\\{p,\pi\\}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|},$

with $C>0$ depending only on $\Gamma$ and $\Omega$. The combination of the
above bound and (3.29) completes the proof of the lemma. ∎

###  3.3 Finite element discretization

Let $\Omega\subset\mathbb{R}^{3}$ be a fixed polygonal domain that strictly
contains $\Gamma$. We consider a family of shape regular tetrahedral
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}=\mbox{diam}(T)$ for
$T\in\mathcal{T}_{h}$. 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\\},$
(3.30)

respectively. We assume $\\{\mathcal{T}_{h}^{\Gamma}\\}_{h>0}$ to be quasi-
uniform. 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)$
(3.31)

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\\}.$
(3.32)

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

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

$\displaystyle L^{2}_{\pm}(\Omega)_{h}=\\{p=(p^{-},p^{+})\in Q_{h}^{-}\times
Q_{h}^{+}\,:\,\int_{\Omega^{-}}p^{-}=\int_{\Omega^{+}}p^{+}=0\\},$

and the test space by $Q_{h}^{\pm}=Q_{h}^{-}\times Q_{h}^{+}\cap
L^{2}_{0}(\Omega)$. Let

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

with the analogous definition for $V_{h}^{+}$. Our FE velocity space is given
by:

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

Functions in $L^{2}_{\pm}(\Omega)$ 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]^{-}_{+}$. Moreover, we define the following averages:

$\\{\mathbf{u}\\}=\alpha\mathbf{u}^{+}+\beta\mathbf{u}^{-},\quad\langle\mathbf{u}\rangle=\beta\mathbf{u}^{+}+\alpha\mathbf{u}^{-},$
(3.33)

where $\alpha$ and $\beta$ are weights to be chosen such that
$\alpha+\beta=1$, $0\leq\alpha,\beta\leq 1$. For example, in [14] 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 (3.37) and
(3.42) we will use relationship:

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

For the discretization of the surface Stokes problem, we first consider the
generalized Taylor–Hood bulk spaces in the strip $\Omega_{h}^{\Gamma}$:

$\displaystyle V_{\Gamma,h}$ $\displaystyle=\\{\mathbf{U}\in
C(\Omega_{h}^{\Gamma})^{3}:\mathbf{U}|_{T}\in\mathbf{P}_{k+1}(T)~{}\forall
T\in\mathcal{T}^{\Gamma}_{h}\\},$ $\displaystyle Q_{\Gamma,h}$
$\displaystyle=\\{\pi\in C(\Omega_{h}^{\Gamma}):\pi|_{T}\in P_{k}(T)~{}\forall
T\in\mathcal{T}^{\Gamma}_{h}\\},$

$Q_{\Gamma,h}^{0}=Q_{\Gamma,h}\cap L^{2}_{0}(\Gamma)$. In the trace finite
element method, we use the traces of functions from $V_{\Gamma,h}$ and
$Q_{\Gamma,h}$ on $\Gamma$. The inf-sup stability of the resulting trace FEM
was analyzed in [41] for $k=1$ and extended to higher order isoparametric
trace elements in [29].

In the treatment of the surface Stokes problem, one has to enforce the
tangentiality condition $\mathbf{U}\cdot\mathbf{n}=0$ on $\Gamma$. In order to
enforce it while avoiding locking, we follow [23, 24, 28, 47, 39] and add a
penalty term to the weak formulation.

A discrete variational analogue of problem (3.12) reads: Find
$(\mathbf{u}_{h},p_{h})\in V^{\pm}_{h}\times L^{2}_{\pm}(\Omega)_{h}$, and
$(\mathbf{U}_{h},\pi_{h})\in V_{\Gamma,h}\times Q_{\Gamma,h}$ such that

$\left\\{\begin{split}a_{h}(\\{\mathbf{u}_{h},\mathbf{U}_{h}\\},\\{\mathbf{v}_{h},\mathbf{V}_{h}\\})+b_{h}(\\{\mathbf{v}_{h},\mathbf{V}_{h}\\},\\{p_{h},\pi_{h}\\})+s_{h}(\mathbf{v}_{h},\pi_{h})&=r_{h}(\mathbf{v}_{h},\mathbf{V}_{h})\\\
b_{h}(\\{\mathbf{u}_{h},\mathbf{U}_{h}\\},\\{q_{h},\tau_{h}\\})-b_{p}(p_{h},q_{h})-b_{s}(\pi_{h},\tau_{h})&=0\end{split}\right.$
(3.35)

for all $(\mathbf{v}_{h},q_{h})\in V_{0,h}^{\pm}\times Q_{h}^{\pm}$ and
$(\mathbf{V}_{h},\tau_{h})\in V_{\Gamma,h}\times Q_{\Gamma,h}^{0}$. We define
all the bilinear forms in (3.35) for all $\mathbf{u}_{h}\in V_{h}^{\pm}$,
$\mathbf{v}_{h}\in V^{\pm}_{0,h}$, $\mathbf{U},\mathbf{V}\in V_{\Gamma,h}$,
$p\in L^{2}(\Omega)$, $\pi\in L^{2}(\Gamma)$. Let us start from form
$a_{h}(\\{\cdot,\cdot\\},\\{\cdot,\cdot\\})$:

$\displaystyle
a_{h}(\\{\mathbf{u}_{h},\mathbf{U}_{h}\\},\\{\mathbf{v}_{h},\mathbf{V}_{h}\\})=$
$\displaystyle
a_{i}(\\{\mathbf{u}_{h},\mathbf{U}_{h}\\},\\{\mathbf{v}_{h},\mathbf{V}_{h}\\})+a_{n}(\mathbf{u}_{h},\mathbf{v}_{h})$
$\displaystyle+a_{p}(\\{\mathbf{u}_{h},\mathbf{U}_{h}\\},\\{\mathbf{v}_{h},\mathbf{V}_{h}\\})+a_{s}(\mathbf{U}_{h},\mathbf{V}_{h}),$
(3.36)

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{U}_{h}\\},\\{\mathbf{v}_{h},$
$\displaystyle\mathbf{V}_{h}\\})=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_{+}}$
$\displaystyle+\langle
f^{-}(\mathbf{P}\mathbf{u}_{h}^{-}-\mathbf{U}_{h}),\mathbf{P}\mathbf{v}_{h}^{-}-\mathbf{V}_{h}\rangle_{\Gamma}+\langle
f^{+}(\mathbf{P}\mathbf{u}_{h}^{+}-\mathbf{U}_{h}),\mathbf{P}\mathbf{v}_{h}^{+}-\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(\mu_{\Gamma}\mathbf{D}_{\Gamma}(\mathbf{U}_{h}),\mathbf{D}_{\Gamma}(\mathbf{V}_{h}))_{\Gamma}$
(3.37)

the terms that enforce condition (3.6) 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\\}([\mathbf{u}_{h}\cdot\mathbf{n}],[\mathbf{v}_{h}\cdot\mathbf{n}])_{\Gamma}-2\langle\\{\mu\mathbf{n}^{T}\mathbf{D}(\mathbf{v})\mathbf{n}\\},[\mathbf{u}_{h}\cdot\mathbf{n}]\rangle_{\Gamma},$
(3.38)

and the stabilization and penalty terms:

$\displaystyle
a_{p}(\\{\mathbf{u}_{h},\mathbf{U}_{h}\\},\\{\mathbf{v}_{h},\mathbf{V}_{h}\\})$
$\displaystyle=\mathbf{J}_{h}^{-}(\mathbf{u}_{h},\mathbf{v}_{h})+\mathbf{J}_{h}^{+}(\mathbf{u}_{h},\mathbf{v}_{h})+\tau_{s}(\mathbf{U}_{h}\cdot\mathbf{n},\mathbf{V}_{h}\cdot\mathbf{n})_{\Gamma},$
(3.39) $\displaystyle\mathbf{J}_{h}^{\pm}(\mathbf{u}_{h},\mathbf{v}_{h})$
$\displaystyle=\sum_{\ell=1}^{k+1}|e|^{2\ell-1}\sum_{e\in\mathcal{E}_{h}^{\Gamma,\pm}}\gamma_{\mathbf{u}}^{\pm}\mu^{\pm}([\partial_{n}^{\ell}\mathbf{u}_{h}^{\pm}],[\partial_{n}^{\ell}\mathbf{u}_{h}^{\pm}])_{e}.$
(3.40)

In (3.40), $\partial_{n}^{\ell}\mathbf{u}_{h}^{-}$ denotes the derivative of
order $\ell$ of $\mathbf{u}_{h}^{-}$ in the direction of $\mathbf{n}$. The
$\mathbf{J}_{h}$ terms in (3.39) are so called ghost-penalty stabilization
[11, 12] included to avoid poorly conditioned algebraic systems due to
possible small cuts of tetrahedra from $\mathcal{T}^{\Gamma}_{h}$ by the
interface. The terms in (3.41) and (3.43) have the same role for the surface
bilinear forms.

The last form in (3.36) is related to the algebraic stability of the surface
Stokes problem:

$\displaystyle
a_{s}(\mathbf{U}_{h},\mathbf{V}_{h})=\rho_{\mathbf{u}}(\nabla\mathbf{u}_{h}\mathbf{n},\nabla\mathbf{v}_{h}\mathbf{n})_{\Omega_{h}^{\Gamma}}.$
(3.41)

Similarly, the terms coming from the integration by parts of the divergence of
the stress tensors are contained in

$\displaystyle
b_{h}(\\{\mathbf{v}_{h},\mathbf{V}_{h}\\},\\{p_{h},\pi_{h}\\})=$
$\displaystyle-(p_{h}^{-},\mathop{\rm
div}\mathbf{v}_{h}^{-})_{\Omega_{-}}-(p_{h}^{+},\mathop{\rm
div}\mathbf{v}_{h}^{+})_{\Omega_{+}}$
$\displaystyle+\left\langle\\{p_{h}\\},[\mathbf{v}_{h}\cdot\mathbf{n}]\right\rangle_{\Gamma}+(\nabla_{\Gamma}\pi_{h},\mathbf{V}_{h})_{\Gamma},$
(3.42)

the penalty terms are grouped together in

$b_{p}(p_{h},q_{h})=J_{h}^{-}(p_{h},q_{h})+J_{h}^{+}(p_{h},q_{h}),\quad
J_{h}^{\pm}(p_{h},q_{h})=\frac{\gamma_{p}^{\pm}}{\mu^{\pm}}\sum_{e\in\mathcal{E}_{h}^{\Gamma,\pm}}\sum_{\ell=1}^{k}|e|^{2\ell+1}([\partial_{n}^{\ell}p_{h}^{\pm}],[\partial_{n}^{\ell}q_{h}^{\pm}])_{e},$

and we have a term related to algebraic stability of the surface Stokes
problem in:

$\displaystyle b_{s}(\pi_{h},\tau_{h})=\rho_{p}(\nabla
p_{h}\cdot\mathbf{n},\nabla p_{h}\cdot\mathbf{n})_{\Omega_{h}^{\Gamma}}.$
(3.43)

Finally,

$\displaystyle s_{h}(\mathbf{v}_{h},\pi_{h})=$
$\displaystyle-\langle\pi_{h}\kappa,\left\langle\mathbf{v}_{h}\cdot\mathbf{n}\right\rangle\rangle_{\Gamma},$
$\displaystyle r_{h}(\mathbf{v}_{h},\mathbf{V}_{h})=$
$\displaystyle(\mathbf{f}_{h}^{-},\mathbf{v}_{h}^{-})_{\Omega_{-}}+(\mathbf{f}_{h}^{+},\mathbf{v}_{h}^{+})_{\Omega_{+}}+(\mathbf{P}\mathbf{b}_{h}^{e},\mathbf{V}_{h})_{\Gamma}.$

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

Parameters $\gamma_{\mathbf{u}}^{\pm}$, $\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 (3.38) needs to be
large enough to provide the bilinear form
$a_{h}(\\{\cdot,\cdot\\},\\{\cdot,\cdot\\})$ with coercivity. Parameters
$\gamma_{\mathbf{u}}^{\pm}$ and $\gamma_{p}^{\pm}$ can be tuned to improve the
numerical performance of the method. As for the parameters required by the
discretization of the surface Stokes problem, we allow:

$\displaystyle\tau_{s}=c_{\tau}h^{-2},\quad\rho_{p}=c_{p}h,\quad\rho_{\mathbf{u}}\in[c_{\mathbf{u}}h,C_{\mathbf{u}}h^{-1}],$
(3.44)

where $c_{\tau}$, $c_{p}$, $c_{\mathbf{u}}$, and $C_{\mathbf{u}}$ are positive
constants independent of $h$ and how $\Gamma$ cuts the bulk mesh.

The definition of bilinear forms requires integration over $\Gamma\cap T$ and
$T\cap\Omega^{\pm}$ for $T$ from $\Omega_{h}^{\Gamma}$. In general, there are
no exact quadrature formulas to accomplish this task [43]. In practice,
approximations should be made which introduce geometric errors. To keep these
geometric errors of the order consistent with the approximation properties of
the finite element spaces, we use isoparametric variants of the above spaces
introduced in [32]; see also [34, 20].

We expect that the stability of the finite element formulation can be analyzed
largely following the same steps of the well-posedness analysis for the weak
formulation in Sec. 3.2, with a special treatment of cut elements, Nitsche
terms and surface elements as available in the literature for bulk Stokes
interface and surface Stokes problems. For the sake of brevity, we do not work
out these details here, but will present them in a follow-up paper.

##  4 A partitioned method for the coupled bulk-surface flow

For the solution of the coupled problem described in Sec. 3 we intend to use a
partitioned strategy, i.e. each sub-problem is solved separately and the
coupling conditions are enforced in an iterative fashion. Partitioned method
are appealing for solving coupled problems because they allow to reuse
existing solvers with minimal modifications. In order to devise such a method
for the simplified problem in Sec. 3, let us take a step back and look at the
original problem (2.12).

Discretize problem (2.12) in time with, e.g., the Backward Euler method and
consider the coupled problem at a particular time $t=t^{n+1}$. Let $S_{b}$ be
the map that associates the jump in the normal stress across the interface to
any given surface flow velocity $\mathbf{U}=\mathbf{U}_{T}+U_{N}\mathbf{n}$:

$S_{b}(\mathbf{U})=[\boldsymbol{\sigma}\mathbf{n}]^{+}_{-}=\boldsymbol{\sigma}^{+}(\mathbf{u}^{+},p^{+})\mathbf{n}-\boldsymbol{\sigma}^{-}(\mathbf{u}^{-},p^{-})\mathbf{n}\quad\text{on}~{}\Gamma,$

where $(\mathbf{u}^{+},p^{+})$ and $(\mathbf{u}^{-},p^{-})$ represent the
solution of the two-phase time-discrete Navier-Stokes problem at time $t$
associated to (2.1)–(2.2) endowed with interface conditions (2.7), (2.9), and
(2.10). Moreover, let $S_{s}$ be the operator associated to the surface flow
such that to any given surface flow velocity $\mathbf{U}$ it associates the
load $\mathbf{f}_{\Gamma}$:

$S_{s}(\mathbf{U})=\mathbf{f}_{\Gamma}\quad\text{on}~{}\Gamma,$

through the time-discrete surface Navier-Stokes problem at time $t$ associated
to (2.5)–(2.6). Note that $S_{b}$ and $S_{s}$ are nonlinear and their
definitions can involve also forcing terms and, in the case of the bulk fluid
problem, terms due to the boundary conditions. For the surface operator, we
can define $S^{-1}_{s}$ as the map that associates the surface flow velocity
$\mathbf{U}$ to any given load $\mathbf{f}_{\Gamma}$ on $\Gamma$.

With the above definitions, we can express the time discrete version of
coupled problem (2.12) in terms of the solution $\mathbf{U}$ of a nonlinear
equation defined only on $\Gamma$. This interface equation is usually
presented in one of three formulations that are equivalent from the
mathematical point of view, but give rise to different iterative algorithms.
The first and perhaps most used formulation is the fixed-point one: Find
$\mathbf{U}$ such that

$S_{s}^{-1}(S_{b}(\mathbf{U}))=\mathbf{U}\quad\text{on}~{}\Gamma.$ (4.1)

The second formulation is a slight modification of (4.1), which lends itself
to a Newton iterative method: Find $\mathbf{U}$ such that

$S_{s}^{-1}(S_{b}(\mathbf{U}))-\mathbf{U}=\boldsymbol{0}\quad\text{on}~{}\Gamma.$

The third approach is given by the Steklov-Poincaré equation: Find
$\mathbf{U}$ such that

$S_{b}(\mathbf{U})-S_{s}(\mathbf{U})=\boldsymbol{0}\quad\text{on}~{}\Gamma.$

See, e.g., [44] for more details on these three formulations.

A standard algorithm for eq. (4.1) uses fixed-point iterations: Given
$\mathbf{U}^{k}$, compute

$\mathbf{U}^{k+1}=\mathbf{U}^{k}+\omega^{k}(\overline{\mathbf{U}}^{k}-\mathbf{U}^{k})\quad\text{with}~{}\overline{\mathbf{U}}^{k}=S_{s}^{-1}(S_{b}(\mathbf{U}^{k})).$
(4.2)

The choice of the relaxation parameter $\omega^{k}$ determines the efficiency
of the algorithm or it might be crucial for convergence in certain ranges of
the physical parameters. An effective strategy for setting $\omega^{k}$ is the
Aitken’s acceleration method.

For simplicity, we present algorithm (4.2) applied to the time discrete
version of coupled problem (2.12) with $\omega^{k}=1$ for all $k$ (i.e., no
relaxation). At time $t=t^{n+1}$, assuming that $\mathbf{U}^{k}$ is known,
perform the following steps:

  * -

Step 1: solve the two-phase time-discrete Navier-Stokes problem at time $t$
associated to (2.1)–(2.2) for the bulk flow variables
$(\mathbf{u}^{-}_{k+1},p^{-}_{k+1})$ and $(\mathbf{u}^{+}_{k+1},p^{+}_{k+1})$
with interface conditions

$\displaystyle\mathbf{u}^{+}_{k+1}\cdot\mathbf{n}$
$\displaystyle=U_{N}^{k}=\mathbf{u}^{-}_{k+1}\cdot\mathbf{n}$
$\displaystyle\text{on }\Gamma$
$\displaystyle\mathbf{P}{\boldsymbol{\sigma}^{+}_{k+1}\mathbf{n}}$
$\displaystyle=f^{+}(\mathbf{P}\mathbf{u}^{+}_{k+1}-\mathbf{U}_{T}^{k})$
$\displaystyle\text{on }\Gamma,$
$\displaystyle\mathbf{P}{\boldsymbol{\sigma}^{-}_{k+1}\mathbf{n}}$
$\displaystyle=-f^{-}(\mathbf{P}\mathbf{u}^{-}_{k+1}-\mathbf{U}_{T}^{k})$
$\displaystyle\text{on }\Gamma.$

  * -

Step 2: solve the time-discrete surface Navier-Stokes problem at time $t$
associated to (2.5)–(2.6) for variables $(\mathbf{U}^{k+1},\pi^{k+1})$ with
interface condition

$\displaystyle\mathbf{f}_{\Gamma}^{k+1}$
$\displaystyle=[\boldsymbol{\sigma}_{k+1}\mathbf{n}]^{+}_{-}$
$\displaystyle\text{ on }\Gamma.$

  * -

Step 3: Check the stopping criterion

$\displaystyle||\mathbf{U}^{k+1}-\mathbf{U}^{k}||_{\Gamma}<\epsilon||\mathbf{U}^{k}||_{\Gamma},$

where $\epsilon$ is a given stopping tolerance.

Notice that the bulk and surface flow problems are solved separately and
sequentially. In general, this algorithm is easy to implement but convergence
could be slow in certain ranges of the physical parameters and require
relaxation for speed-up.

The above algorithm adapted to simplified problem (3.1)–(3.7), reads as
follows. At iteration $k+1$, assuming that $(\mathbf{U}_{T}^{k},\pi^{k})$ are
known, perform the following steps:

  * -

Step 1: solve two-phase problem (3.1)–(3.2) for the bulk flow variables
$(\mathbf{u}^{-}_{k+1},p^{-}_{k+1})$ and $(\mathbf{u}^{+}_{k+1},p^{+}_{k+1})$
with interface conditions

$\displaystyle\mathbf{u}^{+}_{k+1}\cdot\mathbf{n}$
$\displaystyle=\mathbf{u}^{-}_{k+1}\cdot\mathbf{n}$ $\displaystyle\text{on
}\Gamma$ (4.3)
$\displaystyle\mathbf{P}{\boldsymbol{\sigma}^{+}_{k+1}\mathbf{n}}$
$\displaystyle=f^{+}(\mathbf{P}\mathbf{u}^{+}_{k+1}-\mathbf{U}_{T}^{k})$
$\displaystyle\text{on }\Gamma,$ (4.4)
$\displaystyle\mathbf{P}{\boldsymbol{\sigma}^{-}_{k+1}\mathbf{n}}$
$\displaystyle=-f^{-}(\mathbf{P}\mathbf{u}^{-}_{k+1}-\mathbf{U}_{T}^{k})$
$\displaystyle\text{on }\Gamma,$ (4.5)
$\displaystyle\left[\mathbf{n}^{T}\boldsymbol{\sigma}_{k+1}\mathbf{n}\right]^{-}_{+}$
$\displaystyle=\pi^{k}\kappa$ $\displaystyle\text{on }\Gamma.$ (4.6)

  * -

Step 2: solve surface flow problem (3.3)–(3.4) for variables
$(\mathbf{U}_{T}^{k+1},\pi^{k+1})$ with interface condition

$\displaystyle\mathbf{P}\mathbf{f}_{\Gamma}^{k+1}$
$\displaystyle=[\mathbf{P}\boldsymbol{\sigma}_{k+1}\mathbf{n}]^{+}_{-}$
$\displaystyle\text{ on }\Gamma.$ (4.7)

  * -

Step 3: Check the stopping criterion

$\displaystyle||\mathbf{U}_{T}^{k+1}-\mathbf{U}_{T}^{k}||_{\Gamma}<\epsilon||\mathbf{U}_{T}^{k}||_{\Gamma}.$
(4.8)

Notice that only interface condition (4.4)–(4.7) are coupling conditions for
bulk and surface flows. If one was to compute the load exerted on the surface
fluid in (4.7) directly from the solution of the problem at Step 1, the
overall accuracy of the method would be spoiled. Instead, one can compute
$\mathbf{P}\mathbf{f}_{\Gamma}^{k+1}$ by plugging (4.4)–(4.5) into (4.7):

$\displaystyle\mathbf{P}\mathbf{f}_{\Gamma}^{k+1}$
$\displaystyle=f^{+}\mathbf{P}\mathbf{u}^{+}_{k+1}+f^{-}\mathbf{P}\mathbf{u}^{-}_{k+1}-(f^{+}+f^{-})\mathbf{U}_{T}^{k}$
$\displaystyle\text{ on }\Gamma.$

However, we prefer to use a more implicit version of the above condition:

$\displaystyle\mathbf{P}\mathbf{f}_{\Gamma}^{k+1}$
$\displaystyle=f^{+}\mathbf{P}\mathbf{u}^{+}_{k+1}+f^{-}\mathbf{P}\mathbf{u}^{-}_{k+1}-(f^{+}+f^{-})\mathbf{U}_{T}^{k+1}$
$\displaystyle\text{ on }\Gamma.$

since it could help have a better control of approximate rigid rotations
(Killing vector fields).

##  5 Numerical results

The aim of the numerical results collected in this section is to provide
evidence of the robustness of the proposed finite element approach with
respect to the contrast in viscosity in the bulk fluid, surface fluid
viscosity, value of the slip coefficients, and position of the interface
relative to the fixed computational mesh.

For the averages in (3.33), we set $\alpha=0$ and $\beta=1$ for all the
numerical experiments since we have $\mu_{-}\leq\mu_{+}$. In addition, we set
$\gamma^{\pm}_{\mathbf{u}}=0.05$, $\gamma_{p}^{\pm}=0.05$, and $\gamma=80$.
The value of all other parameters will depend on the specific test. The
stopping tolerance for criterion (4.8) is set to $\epsilon=10^{-6}$. For all
the simulations, we choose to use finite element pair $\mathbf{P}_{2}-P_{1}$
for both the bulk and surface fluid problems.

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

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

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

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

Such weighted norm naturally arise in the error analysis of the Stokes
interface problem [42]. In addition, we will report the $L^{2}$ and $H^{1}$
errors for the surface velocity and $L^{2}$ error for the surface pressure.

###  5.1 Sphere embedded in a cube

We perform a series of tests where domain $\Omega$ is the cube
$[-1.5,1.5]^{3}$ and interface $\Gamma$ is the unit sphere centered at the
origin. Let $\mathbf{x}=(x,y,z)\in\Omega$. Surface $\Gamma$ is characterized
as the zero level set of function $\phi(\mathbf{x})=||\mathbf{x}||^{2}_{2}-1$.
We consider the following solution for the bulk flow:

$\displaystyle p^{-}$
$\displaystyle=3x\sqrt{x^{2}+y^{2}+z^{2}}-2x(x^{2}+y^{2}+z^{2}),\quad\mathbf{u}^{-}=\frac{2f^{-}}{f^{-}-\mu^{-}}\mathbf{a}(x,y,z),$
(5.3) $\displaystyle p^{+}$
$\displaystyle=6x\sqrt{x^{2}+y^{2}+z^{2}}-4x(x^{2}+y^{2}+z^{2}),\quad\mathbf{u}^{+}=\frac{2f^{+}}{f^{+}+\mu^{+}}\mathbf{a}(x,y,z),$
(5.4)

where

$\displaystyle\mathbf{a}(x,y,z)=\left(\frac{3}{2}-\sqrt{x^{2}+y^{2}+z^{2}}\right)\left[\begin{array}[]{c}(-y-z)x+y^{2}+z^{2}\\\
(-x-z)y+x^{2}+z^{2}\\\ (-x-y)z+y^{2}+x^{2}\end{array}\right],$ (5.8)

coupled to the following exact solution for the surface flow:

$\displaystyle\pi=x,\quad\mathbf{U}=\left[\begin{array}[]{c}(-y-z)x+y^{2}+z^{2}\\\
(-x-z)y+x^{2}+z^{2}\\\ (-x-y)z+y^{2}+x^{2}\end{array}\right],$ (5.12)

The forcing terms $\mathbf{f}^{-}$ and $\mathbf{f}^{+}$ are found by plugging
the solution (5.3)–(5.4) in (3.1). We impose a Dirichlet condition (2.13) on
the faces $x=1.5$, $y=-1.5$, $z=-1.5$, where function $\mathbf{g}$ is found
from $\mathbf{u}^{+}$ in (5.4). On the remaining part of the boundary, we
impose a Neumann condition (2.14) where $\mathbf{f}_{N}$ is found from $p^{+}$
in (5.3) and $\mathbf{u}^{+}$ in (5.4).

The value of the physical parameters will be specified for each test.

Spatial convergence. To check the spatial accuracy of the finite element
method described in Sec. 3.3, we consider exact solution (5.3)–(5.12) with
viscosities $\mu^{-}=1$, $\mu^{+}=10$ and $\mu_{\Gamma}=1$, and friction
coefficients $f^{-}=2$ and $f^{+}=10$. Notice that the fluid outside the
sphere has a larger viscosity than the fluid inside the sphere, which has the
same viscosity as the surface fluid. We consider structured meshes of
tetrahedra with five levels of refinement, the coarsest mesh having mesh size
$h=0.5$ while the finest mesh has $h=0.05$. All the meshes feature a local
one-level refinement near the corners of $\Omega$. Table 1 reports the number
of DOFs for each mesh. Fig. 2 (left) shows the $L^{2}$ error and weighted
$H^{1}$ error (5.1) for the bulk velocity, weighted $L^{2}$ error (5.2) for
the bulk pressure, $L^{2}$ and $H^{1}$ errors for the surface velocity and
$L^{2}$ error for the surface pressure against the mesh size $h$. We observe
optimal convergence rates for all the norms under consideration. Fig. 2
(right) shows the number of bulk-surface iterations to satisfy stopping
criterion (4.8) as $h$ varies. As we can see, the number of iterations is
fairly insensitive to a mesh refinement or coarsening.

$h$ | 0.5 | 0.25 | 0.125 | 0.0625 | 0.05
---|---|---|---|---|---
# bulk velocity DOFs | $1.1e4$ | $7.4e4$ | $5.2e5$ | $3.6e6$ | $6.4e6$
# bulk pressure DOFs | $6.2e2$ | $3.7e3$ | $2.3e4$ | $1.6e5$ | $2.8e5$
# surface velocity DOFs | $2.4e3$ | $1.0e4$ | $4.0e4$ | $1.5e5$ | $2.2e5$
# surface pressure DOFs | $1.4e2$ | $5.9e2$ | $2.3e3$ | $8.5e3$ | $1.3e4$
Table 1: Sphere: DOFs for bulk and surface variables for all the meshes under
consideration in the spatial convergence test.

Fig. 2: Sphere: (left) Bulk and surface FE errors against the mesh size $h$.
(right) Number of bulk-surface iterations of the partitioned method as $h$
varies.

Robustness with respect to the viscosity contrast. It is known that 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 with respect to the viscosity contrast in the bulk, we consider exact
solution (5.3)–(5.12) and fix $\mu^{-}=1$, while we let $\mu^{+}$ vary from 1
to 256. We set $\mu_{\Gamma}=1$ and friction coefficients $f^{-}=2$ and
$f^{+}=10$.

We consider one of the meshes adopted for the previous sets of simulations
(with $h=0.125$). Fig. 3 (left) shows the $L^{2}$ error and weighted $H^{1}$
error (5.1) for the bulk velocity, weighted $L^{2}$ error (5.2) for the bulk
pressure, $L^{2}$ and $H^{1}$ errors for the surface velocity and $L^{2}$
error for the surface pressure against the value of $\mu^{+}$. We see that the
errors remain mostly unchanged as $\mu^{+}$ varies, with the exception of the
weighted $L^{2}$ error for the bulk pressure, which decreases as $\mu^{+}$
increases. In [40], which focuses only on two-phase bulk flow, we found that
such error reaches a plateau as $\mu^{+}$ is further increased. Fig. 3 (left)
shows that our approach is substantially robust with respect to the viscosity
contrast $\mu^{+}/\mu^{-}$.

Fig. 3: Sphere: (left) Bulk and surface FE errors against the value of
$\mu^{+}$. (right) Number of bulk-surface iterations of the partitioned method
as $\mu^{+}$ varies.

Fig. 3 (right) reports the number of bulk-surface iterations to satisfy
stopping criterion (4.8) as $\mu^{+}$ varies. We observe that the number of
iterations increases as the $\mu^{+}/\mu^{-}$ ratio decreases, indicating that
the coupled bulk-surface problem becomes more stiff as $\mu^{+}$ decreases to
match $\mu^{-}$ and $\mu_{\Gamma}$.

Robustness with respect to the value of the surface viscosity. We now let
$\mu_{\Gamma}$ vary from 1 to 256 and keep all the other physical parameters
fixed to the following values: $\mu^{-}=1$, $\mu^{+}=10$, $f^{-}=2$ and
$f^{+}=10$. Again, we consider exact solution (5.3)–(5.12) and the mesh with
mesh size $h=0.125$. Fig. 4 (left) shows all the errors we have considered so
far against the value of $\mu_{\Gamma}$. We notice that all the bulk errors
stay constant as $\mu_{\Gamma}$ varies. The $L^{2}$ errors for the surface
velocity and pressure increase as $\mu_{\Gamma}$ increases, while the $H^{1}$
error for the surface velocity slightly decreases as $\mu_{\Gamma}$ increases.
This experiment suggests that more viscous embedded layer is less controlled
by the bulk fluid which effects the numerical stability of the complete
system. In a water – lipid membrane system, the ratio of lateral dynamic
viscosities of the embedded bi-layer and bulk water is typically 1–10 $\mu m$
(depending on the temperature and composition) with the size of a vesicle
being generally between 0.1 and 10 $\mu m$. Hence the observed increase of the
numerical error does not look critical for this application.

Fig. 4: Sphere: (left) Bulk and surface FE errors against the value of
$\mu_{\Gamma}$. (right) Number of bulk-surface iterations of the partitioned
method as $\mu_{\Gamma}$ varies.

Fig. 4 (right) shows the number of bulk-surface iterations to satisfy stopping
criterion (4.8) as $\mu_{\Gamma}$ varies. Our partitioned method seems to be
insensitive to a variation in the value of $\mu_{\Gamma}$. In particular, for
the range of $\mu_{\Gamma}$ under consideration the number of iterations stays
constant and equal to 12.

Robustness with respect to the slip coefficients. To check the sensitivity of
the errors and partitioned method to the value of the slip coefficients, we
run two sets of experiments, both involving exact solution (5.3)–(5.12). In
the first set we fix $f^{+}=2$ and vary $f^{-}$ from 1 to 256, while in the
second set we take $f^{+}=f^{-}$ and let them both vary from 1 to 256 . The
viscosities are set as follows: $\mu^{-}=1$, $\mu^{+}=10$, and
$\mu_{\Gamma}=1$. We consider again the mesh with mesh size $h=0.125$. Fig. 5
(left) and 6 (left) show all the errors under consideration against the value
of the slip coefficient(s) for both sets of tests. The only error that shows a
substantial variation is the weighted $H^{1}$ error the bulk velocity, which
increases as the slip coefficient(s) increase. However, such error seems to
reach a plateau in both cases.

Fig. 5: Sphere: (left) Bulk and surface FE errors against the value of
$f^{+}$. (right) Number of bulk-surface iterations of the partitioned method
as $f^{+}$ varies.

Fig. 5 (right) and 6 (right) report the number of bulk-surface iterations to
satisfy stopping criterion (4.8) as the value of the coefficient(s) varies for
both sets of tests. In Fig. 5 (right), we see a rather sharp increase in the
number of iterations as $f^{-}$ increases. This is even more true when both
slip coefficients are increased together, as we can see from Fig. 6 (right).
Fig. 7 reports the relative difference of the surface velocity between
subsequent iterations in $L^{2}$ norm until stopping criterion (4.8) is met
for $f^{+}=f^{-}=2^{2}$ and $f^{+}=f^{-}=2^{8}$. We see that such relative
difference decreases regularly for $f^{+}=f^{-}=2^{2}$, while for
$f^{+}=f^{-}=2^{8}$ it decreases quickly for the first few iterations and then
it slows down. A heuristic explanation we have for this is that as the two
friction coefficients increase interface conditions (2.9)–(2.10) become close
to Dirichlet conditions, making the surface flow more “passive”. Thus,
separating the surface flow from the bulk flow as in the partitioned algorithm
might not make much sense.

Fig. 6: Sphere: (left) Bulk and surface FE errors against the value of
$f^{+}=f^{-}$. (right) Number of bulk-surface iterations of the partitioned
method as the value of $f^{+}$ and $f^{-}$ (with $f^{+}=f^{-}$) varies. Fig.
7: Sphere: relative difference of the surface velocity between subsequent
iterations in $L^{2}$ norm until stopping criterion (4.8) is met.

###  5.2 Torus embedded in a cube

The domain $\Omega$ is cube $[-2,2]^{2}$ and surface $\Gamma$ is a torus
centered at $\mathbf{c}=(c_{1},c_{2},c_{3})$. Let
$(x,y)=(\tilde{x}-c_{1},\tilde{y}-c_{2},\tilde{z}-c_{3})$,
$(\tilde{x},\tilde{y},\tilde{z})\in\Omega$. We can characterize $\Gamma$ as
the zero level set of function
$\phi(\mathbf{x})=\sqrt{z^{2}+(\sqrt{x^{2}+y^{2}}-1)^{2}}-\frac{1}{2}$.
Finding an exact solution problem (3.1)–(3.5), (2.9), and (2.10) with this
more complicated surface is highly non-trivial. To simplify the task, we relax
interface conditions (2.9), (2.10), and (3.5) as follows:

$\displaystyle\mathbf{P}{\boldsymbol{\sigma}^{\pm}\mathbf{n}}$
$\displaystyle=\pm f^{\pm}(\mathbf{P}\mathbf{u}^{\pm}-\mathbf{U})+g^{\pm}$
$\displaystyle\text{ on }\Gamma,$
$\displaystyle\left[\mathbf{n}^{T}\boldsymbol{\sigma}\mathbf{n}\right]^{-}_{+}$
$\displaystyle=\pi\kappa+g^{n}$ $\displaystyle\text{ on }\Gamma.$

where $g^{+}$, $g^{-}$, and $g^{n}$ are computed such that exact solution
given below satisfy these relaxed interface conditions. The solution is given
by

$\displaystyle p^{-}$
$\displaystyle=\left(\frac{1}{2}-\frac{2-4\sqrt{x^{2}+y^{2}}}{\sqrt{x^{2}+y^{2}}}\right)(x^{3}+x),~{}p^{+}=\frac{1}{2}(x^{3}+x),~{}\mathbf{u}^{-}=\mathbf{u}^{+}=\left[\begin{array}[]{c}x^{2}y\\\
5-xy^{2}+z^{2}\\\ -xy\end{array}\right],$ (5.16)

for the bulk and:

$\displaystyle\pi=x^{3}+x,\quad\mathbf{U}=\left[\displaystyle{\frac{-zx}{\sqrt{x^{2}+y^{2}}}},~{}\displaystyle{\frac{-zy}{\sqrt{x^{2}+y^{2}}}},~{}\sqrt{x^{2}+y^{2}}-1\right]^{T},$
(5.17)

for the surface. The forcing terms $\mathbf{f}^{-}$ and $\mathbf{f}^{+}$ are
found by plugging the solution (5.16)–(5.17) in (3.1). We impose a Dirichlet
condition (2.13) on the faces $x=2$, $y=-2$, $z=-2$, where function
$\mathbf{g}$ is found from $\mathbf{u}^{+}$ in (5.16). On the remaining part
of the boundary, we impose a Neumann condition (2.14) where $\mathbf{f}_{n}$
is found from $p^{+}$ and $\mathbf{u}^{+}$ in (5.16).

Spatial convergence. Once again, we start by checking spatial accuracy. To
this end, we consider exact solution (5.16)–(5.17) with $\mathbf{c}=(0,0,0)$,
viscosities $\mu^{-}=1$, $\mu^{+}=10$, $\mu_{\Gamma}=1$, and friction
coefficients $f^{-}=2$ and $f^{+}=10$. Just like in the case of the sphere, we
consider structured meshes of tetrahedra that feature a local one-level
refinement near the corners of $\Omega$. The details of the meshes under
consideration are reported in Table 2. Fig. 8 shows the $L^{2}$ error and
weighted $H^{1}$ error (5.1) for the bulk velocity, weighted $L^{2}$ error
(5.2) for the bulk pressure, $L^{2}$ and $H^{1}$ errors for the surface
velocity and $L^{2}$ error for the surface pressure against the mesh size $h$.
Also for this second convergence test, we observe optimal convergence rates
for all the norms.

$h$ | 0.25 | 0.125 | 0.0625 | 0.05
---|---|---|---|---
# bulk velocity DOFs | $1.6e5$ | $1.2e6$ | $8.5e6$ | $1.5e7$
# bulk pressure DOFs | $7.6e3$ | $5.4e4$ | $3.7e5$ | $6.7e5$
# surface velocity DOFs | $1.6e4$ | $6.0e4$ | $2.3e5$ | $3.4e5$
# surface pressure DOFs | $9.0e2$ | $3.4e3$ | $1.3e4$ | $2.0e4$
Table 2: Torus: DOFs for bulk and surface variables for all the meshes under
consideration in the spatial convergence test. Fig. 8: Torus: Bulk and surface
FE errors against the mesh size $h$.

Robustness with respect to the position of the interface. We conclude our
series of numerical results 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. We vary the center $\mathbf{c}=(c_{1},c_{2},c_{3})$ of
the torus that represents $\Gamma$:

$\displaystyle
c_{1}=h\frac{k}{20}\sin\left(\frac{k\pi}{10}\right),~{}c_{2}=h\frac{k\sqrt{2}}{40}\cos\left(\frac{k\pi}{10}\right),~{}c_{3}=h\frac{k\sqrt{2}}{40}\cos\left(\frac{k\pi}{10}\right),$
(5.18)

where $h$ is the mesh size. The physical parameters are set like in the
convergence test. We consider the mesh in Table 2 with $h=0.125$. Fig. 9 shows
all the errors against the value of $k$ in (5.18). We see that all the errors
are fairly insensitive to the position of $\Gamma$ with respect to the
background mesh, indicating robustness.

Fig. 9: Torus: Bulk and surface FE errors against the value of $k$ in (5.18).

## Acknowlegement

We are grateful to Dr. Christoph Lehrenfeld for providing us with an ngsxfem
implementation of isoparametric unfitted finite elements.

## References

  * [1] The influence of boundary conditions on the contact problem in a 3d.
  * [2] Netgen/NGSolve.  https://ngsolve.org/.
  * [3] ngsxfem.  https://github.com/ngsxfem/ngsxfem/tree/49205a1ae637771a0ed56d4993ce99008f3a00e0.
  * [4] J. W. Barrett, H. Garcke, and R. Nürnberg, A stable numerical method for the dynamics of fluidic membranes, Numerische Mathematik, (2014), pp. 1–40.
  * [5] C. Bernardi, C. Canuto, and Y. Maday, Generalized inf-sup conditions for Chebyshev spectral approximation of the Stokes problem, SIAM Journal on Numerical Analysis, 25 (1988), pp. 1237–1271.
  * [6] L. Bocquet and J.-L. Barrat, Flow boundary conditions from nano-to micro-scales, Soft matter, 3 (2007), pp. 685–693.
  * [7] M. E. Bogovskii, Solution of the first boundary value problem for the equation of continuity of an incompressible medium, in Doklady Akademii Nauk, vol. 248, Russian Academy of Sciences, 1979, pp. 1037–1040.
  * [8] N. Bootland, A. Bentley, C. Kees, and A. Wathen, Preconditioners for two-phase incompressible Navier –Stokes flow, SIAM Journal on Scientific Computing, 41 (2019), pp. B843–B869.
  * [9] S. Bordas, E. Burman, M. Larson, and e. M. A. Olshanskii, Geometrically Unfitted Finite Element Methods and Applications, vol. Lecture Notes in Computational Science and Engineering 121, Springer, Berlin, 2018.
  * [10] D. Bothe and J. Prüss, On the two-phase Navier –Stokes equations with Boussinesq–Scriven surface fluid, Journal of Mathematical Fluid Mechanics, 12 (2010), pp. 133–150.
  * [11] E. Burman, Ghost penalty, C. R. Math. Acad. Sci. Paris, 348 (2010), pp. 1217–1220.
  * [12] E. Burman, S. Claus, P. Hansbo, M. G. Larson, and A. Massing, Cutfem: Discretizing geometry and partial differential equations, International Journal for Numerical Methods in Engineering, 104 (2015), pp. 472–501.
  * [13] E. Cáceres, J. Guzmán, and M. Olshanskii, New stability estimates for an unfitted finite element method for two-phase Stokes problem, SIAM Journal on Numerical Analysis, 58 (2020), pp. 2165–2192.
  * [14] S. Claus and P. Kerfriden, A cutfem method for two-phase flow problems, Computer Methods in Applied Mechanics and Engineering, 348 (2019), pp. 185 – 206.
  * [15] M. Cooley and M. O’neill, On the slow motion generated in a viscous fluid by the approach of a sphere to a plane wall or stationary sphere, Mathematika, 16 (1969), pp. 37–49.
  * [16] G. Dziuk and C. M. Elliott, $L^{2}$-estimates for the evolving surface finite element method, Mathematics of Computation, 82 (2013), pp. 1–24.
  * [17] T. Frachon and S. Zahedi, A cut finite element method for incompressible two-phase Navier –Stokes flows, Journal of Computational Physics, 384 (2019), pp. 77 – 98.
  * [18] S. Ganesan, G. Matthies, and L. Tobiska, On spurious velocities in incompressible flow problems with interfaces, Preprint 05-35, Department of Mathematics, University of Magdeburg, 2005.
  * [19] P. Gangl, K. Sturm, M. Neunteufel, and J. Schöberl, Fully and semi-automated shape differentiation in NGSolve, arXiv:2004.06783, (2020).
  * [20] J. Grande, C. Lehrenfeld, and A. Reusken, Analysis of a high-order trace finite element method for pdes on level set surfaces, SIAM Journal on Numerical Analysis, 56 (2018), pp. 228–255.
  * [21] S. Gross and A. Reusken, Numerical methods for two-phase incompressible flows, vol. 40, Springer Science & Business Media, 2011.
  * [22] M. E. Gurtin and A. I. Murdoch, A continuum theory of elastic material surfaces, Archive for Rational Mechanics and Analysis, 57 (1975), pp. 291–323.
  * [23] P. Hansbo, M. Larson, and A. Massing, A stabilized cut finite element method for the darcy problem on surfaces, Computer Methods in Applied Mechanics and Engineering, 326 (2017), pp. 298 – 318.
  * [24] P. Hansbo, M. G. Larson, and K. Larsson, Analysis of finite element methods for vector laplacians on surfaces, IMA Journal of Numerical Analysis, 40 (2020), pp. 1652–1701.
  * [25] P. Hansbo, M. G. Larson, and S. Zahedi, A cut finite element method for a Stokes interface problem, Applied Numerical Mathematics, 85 (2014), pp. 90 – 114.
  * [26] 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).
  * [27] L. M. Hocking, The effect of slip on the motion of a sphere close to a wall and of two adjacent spheres, Journal of Engineering Mathematics, 7 (1973), pp. 207–221.
  * [28] T. Jankuhn, M. A. Olshanskii, and A. Reusken, Incompressible fluid problems on embedded surfaces: Modeling and variational formulations, Interfaces and Free Boundaries, 20 (2018), pp. 353–377.
  * [29] T. Jankuhn, M. A. Olshanskii, A. Reusken, and A. Zhiliakov, Error analysis of higher order trace finite element methods for the surface Stokes equation, Journal of Numerical Mathematics, (2021).
  * [30] V. John, Slip with friction and penetration with resistance boundary conditions for the Navier –Stokes equations—numerical tests and aspects of the implementation, Journal of Computational and Applied Mathematics, 147 (2002), pp. 287–300.
  * [31] E. Lauga, M. Brenner, and H. Stone, Microfluidics: The no-slip boundary condition, in Springer Handbooks, Springer, 2007, pp. 1219–1240.
  * [32] C. Lehrenfeld, High order unfitted finite element methods on level set domains using isoparametric mappings, Computer Methods in Applied Mechanics and Engineering, 300 (2016), pp. 716–733.
  * [33] , A higher order isoparametric fictitious domain method for level set domains, in Geometrically Unfitted Finite Element Methods and Applications, S. P. A. Bordas, E. Burman, M. G. Larson, and M. A. Olshanskii, eds., Cham, 2017, Springer International Publishing, pp. 65–92.
  * [34] C. Lehrenfeld and A. Reusken, Analysis of a high-order unfitted finite element method for elliptic interface problems, IMA Journal of Numerical Analysis, 38 (2018), pp. 1351–1387.
  * [35] A. Massing, M. Larson, A. Logg, and M. Rognes, A stabilized Nitsche overlapping mesh method for the Stokes problem, Numer. Math., 128 (2014), pp. 73 – 101.
  * [36] T.-H. Miura, On singular limit equations for incompressible fluids in moving thin domains, Quarterly of Applied Mathematics, 76 (2018), pp. 215–251.
  * [37] R. A. Nicolaides, Existence, uniqueness and approximation for generalized saddle point problems, SIAM Journal on Numerical Analysis, 19 (1982), pp. 349–357.
  * [38] I. Nitschke, A. Voigt, and J. Wensch, A finite element approach to incompressible two-phase flow on manifolds, Journal of Fluid Mechanics, 708 (2012), pp. 418–7438.
  * [39] M. Olshanskii, A. Quaini, A. Reusken, and V. Yushutin, A finite element method for the surface Stokes problem, SIAM Journal on Scientific Computing, 40 (2018), pp. A2492–A2518.
  * [40] M. Olshanskii, A. Quaini, and Q. Sun, An unfitted finite element method for two-phase Stokes problems with slip between phases, arXiv preprint arXiv:2101.09627, (2021).
  * [41] M. Olshanskii, A. Reusken, and A. Zhiliakov, Inf-sup stability of the trace $\mathbf{P}_{2}$–${P}_{1}$ taylor–hood elements for surface pdes, Mathematics of Computation, (2021).
  * [42] M. A. Olshanskii and A. Reusken, Analysis of a Stokes interface problem, Numerische Mathematik, 103 (2006), pp. 129–149.
  * [43] M. A. Olshanskii and D. Safin, Numerical integration over implicitly defined domains for higher order unfitted finite element methods, Lobachevskii Journal of Mathematics, 37 (2016), pp. 582–596.
  * [44] A. Quarteroni and A. Valli, Domain Decomposition Methods for Partial Differential Equations, Oxford Science Publications, 1999.
  * [45] A. Reusken and Y. Zhang, Numerical simulation of incompressible two-phase flows with a Boussinesq-Scriven surface stress tensor, Int. J. Numerical Methods in Fluids, 73 (2013), pp. 1042–1058.
  * [46] S. Reuther and A. Voigt, The interplay of curvature and vortices in flow on curved surfaces, Multiscale Modeling & Simulation, 13 (2015), pp. 632–643.
  * [47] S. Reuther and A. Voigt, Solving the incompressible surface Navier -Stokes equation by surface finite elements, Physics of Fluids, 30 (2018), p. 012107.
  * [48] D. S. Rodrigues, R. F. Ausas, F. Mut, and G. C. Buscaglia, A semi-implicit finite element method for viscous lipid membranes, Journal of Computational Physics, 298 (2015), pp. 565–584.
  * [49] G. Salbreux and F. Jülicher, Mechanics of active surfaces, Physical Review E, 96 (2017), p. 032404.
  * [50] E. M. Stein, Singular integrals and differentiability properties of functions, vol. 2, Princeton university press, 1970.
  * [51] A. Torres-Sánchez, D. Millán, and M. Arroyo, Modelling fluid deformable surfaces with an emphasis on biological interfaces, Journal of fluid mechanics, 872 (2019), pp. 218–271.
  * [52] N. Wang and J. Chen, A nonconforming Nitsche’s extended finite element method for Stokes interface problems., J Sci Comput, 81 (2019), pp. 342–374.
  * [53] A. Yavari, A. Ozakin, and S. Sadik, Nonlinear elasticity in a deforming ambient space, Journal of Nonlinear Science, 26 (2016), pp. 1651–1692.