#  $\phi$-FEM: an optimally convergent and easily implementable immersed
boundary method for particulate flows and Stokes equations111This work was
supported by the Agence Nationale de la Recherche, Project PhiFEM, under grant
ANR-22-CE46-0003-01.

Michel Duprez222MIMESIS team, Inria Nancy - Grand Est, MLMS team, Université
de Strasbourg, France. michel.duprez@inria.fr,Vanessa Lleras333IMAG, Univ
Montpellier, CNRS, Montpellier, France<EMAIL_ADDRESS>Alexei
Lozinski444Université de Franche-Comté, Laboratoire de Mathématiques de
Besançon, UMR CNRS 6623, 16, route de Gray 25000 Besançon, France.
<EMAIL_ADDRESS>
(August 27, 2024)

###### Abstract

We present an immersed boundary method to simulate the creeping motion of a
rigid particle in a fluid described by the Stokes equations discretized thanks
to a finite element strategy on unfitted meshes, called $\phi$-FEM, that uses
the description of the solid with a level-set function. One of the advantages
of our method is the use of standard finite element spaces and classical
integration tools, while maintaining the optimal convergence (theoretically in
the $H^{1}$ norm for the velocity and $L^{2}$ for pressure; numerically also
in the $L^{2}$ norm for the velocity).

##  1 Introduction

The main goal of the present article is to demonstrate that the recently
proposed $\phi$-FEM methodology [12, 11, 10] is suitable for numerical
simulation of incompressible viscous fluid flow past moving rigid bodies. This
approach allows us to use simple (Cartesian) computational meshes, not
evolving in time and not fitted to the moving rigid bodies, while achieving
the optimal accuracy with classical finite element (FE) spaces of any order
and performing the usual numerical integration on the whole mesh cells and
facets, allowing for the use of standard FEM libraries for the implementation.
We consider here only the creeping motion regime (zero Reynolds number),
neglecting all the inertial terms in the equation governing both the fluid and
the rigid bodies.

Numerical simulations of flows around moving rigid or elastic structures using
immobile simple grids is a popular approach in, for instance, biomechanics,
starting from the work of Peskin [22]. Different approaches have emerged since
then, such as the Immersed Boundary method [18, 21], the Fictitious Domain
method [15, 16], the penalty approximation [2], etc. All these classical
methods suffer from poor accuracy due to the necessity to approximate the
singularities near the fluid-solid interfaces which arise as the artifact of
extending the fluid velocity field inside the solid domain. More recently,
several optimally convergent fictitious domain-type methods have been proposed
for the Stokes equations, which can also be used to simulate the fluid-solid
motions. We cite in particular [8, 20, 17] following the CutFEM paradigm, and
[1, 14] following the X-FEM paradigm. The common feature of all these methods
is that they discretize the variational formulation of the Stokes equation on
the physical fluid domain $\Omega$ using the FE spaces defined on the
background mesh occupying a domain $\Omega_{h}$, slightly larger than
$\Omega$. On the one hand, this permits to avoid a non-smooth extension of the
solution outside its natural domain and to retrieve the optimal accuracy of
the employed finite elements. On the other hand, this introduces integrals on
the cut cells into the FE scheme, i.e. the numerical integration should be
performed on the portions of mesh cells, cut by the fluid-solid interface,
making the methods difficult to implement.

The $\phi$-FEM approach, which is the subject of the present paper, aims at
combining the advantages of both classical Immersed Boundary/Fictitious Domain
methods, and more recent CutFEM/X-FEM. Similarly to the former, $\phi$-FEM
does not need non-standard numerical integration on the cut cells; similarly
to the latter, $\phi$-FEM achieves the optimal accuracy of the finite elements
employed. The general procedure of $\phi$-FEM can be summarized as follows:

  * •

Supposing that the physical domain $\Omega$ is given by a level set function
$\Omega=\\{\phi<0\\}$ and that it is embedded into a simple background mesh,
introduce the active computational mesh $\mathcal{T}_{h}$ by getting rid of
the mesh cells lying completely outside $\Omega$. The active mesh thus
occupies a domain $\Omega_{h}\supset\Omega$, cf. Fig. 1, as in CutFEM/X-FEM.

  * •

Extend the governing equations from $\Omega$ to $\Omega_{h}$ and write down a
formal variational formulation on $\Omega_{h}$ without taking into account the
boundary conditions on $\Gamma$ (the relevant part of the boundary of
$\Omega$).

  * •

Impose the boundary conditions on $\Gamma$ using an appropriate ansatz or some
additional variables, explicitly involving the level set $\phi$ which provides
the link to the actual boundary. For instance, the homogeneous Dirichlet
boundary conditions ($u=0$ on $\Gamma$) can be imposed by the ansatz
${u}=\phi{w}$ thus reformulating the problem in terms of the new unknown
${w}$.

  * •

Add appropriate stabilization, typically combining the ghost penalty [7] with
a least square imposition of the governing equation on the mesh cells
intersected by $\Gamma$, to guarantee coerciveness/stability on the discrete
level.

This program has been successfully carried out for elliptic scalar PDEs with
Dirichlet boundary conditions in [12] and for Neumann boundary conditions in
[11]. Its feasibility is also demonstrated in [10] for the linear elasticity
with mixed boundary conditions including the cases of internal interfaces
between different materials or cracks, and for the heat equation. However, the
adaptation to the equations governing the fluid flow around the moving
particles is not straightforward. In particular, the following challenges are
dealt with in the present article:

  * •

The discrete inf-sup stability theory should be adapted to the case of a non-
standard variational formulation of the Stokes equations posed on $\Omega_{h}$
rather than on $\Omega$, and lacking the saddle-point structure. We shall show
that this is possible by adapting the ghost penalty, which should be taken
slightly more complicated than in the case of scalar elliptic equations [12].
We shall do it here for Taylor-Hood finite elements of any order, but similar
ideas should be also applicable to other classical inf-sup stable FE spaces.

  * •

The motion equations for the solid particles involve the forces exerted on
them by the surrounding fluid. These are defined through the integrals of some
functions of fluid velocity and pressure on the particle boundary. However,
the whole point of $\phi$-FEM is to avoid such integrals. Indeed, the particle
boundary is not resolved by the mesh, and our goal is to provide a method that
necessitates the integration on the whole mesh cells or facets only. The way
out of this paradox, pursued in the present paper, lies in providing a weak
formulation of the governing equations, extended to $\Omega_{h}$, that
incorporates in an appropriate way the force balance equations, without
stating them directly. This formulation is similar in spirit but different
from that in [16].

We note that the method of this article shares some similarities to the
shifted boundary method (SBM) proposed in [19] and analysed in [3] in the case
of Stokes equations. In particular, SBM also gives an optimal accurate
solution (at least with the lowest order finite elements) without introducing
integrals on the cut cells. It is however not evident how one can deal with
the computation of the forces on the particles in the SBM framework.

Figure 1: Left: an example of geometry for the fluid $\Omega$ with a solid
$\mathcal{S}$ inside; Right: the non-conforming active mesh $\mathcal{T}_{h}$
on $\Omega_{h}$ with its internal boundary $G_{h}$.

The paper is structured as follows. In the next section, we properly introduce
the governing equations, develop an appropriate weak formulation, discretize
it (thus introducing our $\phi$-FEM scheme), and announce the main theorem
about the convergence of the scheme. Section 3 is devoted to the proof of this
theorem. As a by-product, we also introduce a $\phi$-FEM approach to
discretize the Stokes equations alone (on a fixed geometry) on a non-fitted
mesh. The details about this (comparatively simple) particular case are given
in Appendix A. In section 4, we illustrate our theoretical results with
numerical examples both for the Stokes equations and for the fluid/rigid
particle motion problem. $\phi$-FEM is also compared there with a standard
(non isoparametric) FEM on fitted meshes, demonstrating the superiority of
$\phi$-FEM in terms of the accuracy achieved on comparable meshes. We end up
with the last section giving some conclusions and perspectives.

Various notations for different parts of geometry and triangulation appear
throughout the article. For readers’ convenience, they are gathered in
Appendix B.

##  2 Construction of the $\phi$-FEM scheme and main results

###  2.1 Governing equations

We consider the motion of a viscous incompressible fluid around a solid
(rigid) particle in the regime of creeping motion, i.e. neglecting all the
inertial terms (for simplicity, we restrict ourselves here to the case of one
particle, the extension to multiple particles being straight-forward). The
particle is mobile and it moves under the action of the forces exerted by the
surrounding fluid and the external forces (gravity). Let the fluid occupy (at
a given time $t$) the domain $\Omega\subset\mathbb{R}^{d}$ ($d=2$ or $3$), the
particle occupy the domain $\mathcal{S}\subset\mathbb{R}^{d}$, and denote
$\mathcal{O}=\Omega\cup\bar{\mathcal{S}}$. Let
$\Gamma_{w}=\partial\mathcal{O}$ be the external boundary of the fluid domain
(the immobile wall) where the fluid velocity is assumed to vanish,
$\Gamma=\partial\mathcal{S}$ be the fluid/solid interface, and assume that
$\Gamma$ does not touch $\Gamma_{w}$, so that $\partial\Omega$ contains two
disjoint components $\Gamma_{w}$ and $\Gamma$. For simplicity, we assume that
the only external body force is gravitation with the constant acceleration
$g$. Hence, the body force density in the fluid is $\rho_{f}g$ where
$\rho_{f}$ is the constant fluid density. Let $\rho_{s}$ be the constant
density of the solid. Then, the resultant external force on the particle is
$mg$ where $m=\rho_{s}|\mathcal{S}|$ is the mass of the particle, and the
resultant moment of the external force with respect to the barycenter of the
particle is 0. Denoting the constant fluid viscosity by $\nu$, the equations
governing the motion of the fluid/particle system can be now given as:

$\displaystyle-2\nu\operatorname{div}D(u)+\nabla p$ $\displaystyle=\rho_{f}g,$
$\displaystyle\text{ in }\Omega\text{}$ (1a)
$\displaystyle\operatorname{div}u$ $\displaystyle=0,$ $\displaystyle\text{ in
}\Omega$ (1b) $\displaystyle u$ $\displaystyle=U+\psi\times r,$
$\displaystyle\text{ on }\Gamma\quad$ (1c) $\displaystyle u$
$\displaystyle=0,$ $\displaystyle\text{ on }\Gamma_{w}$ (1d)
$\displaystyle\int_{\Gamma}(2\nu D(u)-pI)n$ $\displaystyle=mg$ (1e)
$\displaystyle\int_{\Gamma}(2\nu D(u)-pI)n\times r$ $\displaystyle=0$ (1f)
$\displaystyle\int_{\Omega}p$ $\displaystyle=0$ (1g)

Here, the unknowns are the fluid velocity $u:\Omega\to\mathbb{R}^{d}$ and the
pressure $p:\Omega\to\mathbb{R}$, the velocity of the particle barycenter
$U\in\mathbb{R}^{d}$, and the angular velocity of the particle
$\psi\in\mathbb{R}^{d^{\prime}}$ ($d^{\prime}=1$ if $d=2$ and $d^{\prime}=3$
if $d=3$). In these equations, $D(u)=\frac{1}{2}(\nabla u+\nabla u^{T})$
denotes the strain tensor, $r$ denotes the vector from the barycenter of the
solid $\mathcal{S}$, and $n$ denotes the unit normal on $\Gamma$ looking into
the solid. Equations (1e)–(1f) come from the balance of forces exerted on the
particle (the force exerted by the fluid and the gravitational force). 1b

From a numerical simulation perspective, it is natural to introduce an
immobile computational mesh on the immobile box $\mathcal{O}$ containing both
the fluid and the particle. On the other hand, the solid $\mathcal{S}$ will be
moving with velocities $U=U(t)$, $\psi=\psi(t)$ at all-time $t$, thus
permanently changing the shape of the fluid domain $\Omega$. It is therefore
interesting to design numerical methods for the system (1a)–(1g) that
discretize $u$ and $p$ on a mesh non fitted to $\Omega$.

###  2.2 A formal derivation of the appropriate weak formulation

Let $\mathcal{T}_{h}^{\mathcal{O}}$ be a regular simplicial mesh on
$\mathcal{O}$ (the background mesh). Assume that the solid and fluid domains
are given by the level-set function $\phi$: $\mathcal{S}=\\{\phi>0\\}$ and
$\Omega=\mathcal{O}\cap\\{\phi<0\\}$. Introduce the active computational mesh
$\mathcal{T}_{h}$ as a submesh of $\mathcal{T}_{h}^{\mathcal{O}}$ covering
$\Omega$, i.e. excluding the cells of $\mathcal{T}_{h}^{\mathcal{O}}$ lying
completely inside $\mathcal{S}$. Let $\Omega_{h}\supset\Omega$ be the domain
of $\mathcal{T}_{h}$ and $G_{h}$ be the component of $\partial\Omega_{h}$,
other than $\Gamma_{w}$, and thus lying inside $\mathcal{S}$, cf. Fig. 1.555In
practice, the geometrical setting may be slightly more complicated. The
rigorous theoretical definitions of $\Omega_{h}$ and $G_{h}$ will be given in
(8) and (9) and will be based on an approximation $\phi_{h}$ to the levelset
$\phi$, rather than on $\phi$ itself. This may occasionally result in
situations where some tiny portions of $\Omega$ lie outside $\Omega_{h}$ so
that $G_{h}$ slightly penetrates $\Omega$. These technical details are not
important for the forthcoming formal derivation of the FE scheme, while the
rigorous proofs will be done assuming definitions (8) and (9). The actual
implementation may introduce yet more geometrical approximations, as mentioned
in Remark 1, which are not covered by our theory.

Assume (on a formal level, just to derive the scheme) that $u$ and $p$ can be
extended from $\Omega$ to $\Omega_{h}$ as solution to the Stokes equations so
that

$-2\nu\operatorname{div}D(u)+\nabla p=\rho_{f}g\quad\text{ and
}\quad\operatorname{div}u=0\text{ in }\Omega_{h}.$

Taking any sufficiently smooth test functions $v$ and $q$ on $\Omega_{h}$ such
that $v=0$ on $\Gamma_{w}$, an integration by parts gives

$2\nu\int_{\Omega_{h}}D(u):D(v)-\int_{\Omega_{h}}p\operatorname{div}v-\int_{\Omega_{h}}q\operatorname{div}u-\int_{G_{h}}(2\nu
D(u)-pI)n\cdot v=\int_{\Omega_{h}}\rho_{f}g\cdot v.$ (2)

Assuming $u=0$ on $\Gamma_{w}$, this imposes already the boundary condition
(1d) on $\Gamma_{w}$, which we suppose to fit to the mesh $\mathcal{T}_{h}$.
On the contrary, this formulation does not take into account any boundary
conditions on $\Gamma$. In order to incorporate boundary conditions (1c) we
make the ansatz

$u=\phi w+\chi(U+\psi\times r)$ (3)

where $\phi$ is the level-set for $\Omega$ so that $\phi=0$ on $\Gamma$, and
$\chi$ is a sufficiently smooth function on $\mathcal{O}$ such that $\chi=1$
on the solid $\mathcal{S}$ and $\chi=0$ on $\Gamma_{w}$. This introduces the
new vector valued unknown $w$ on $\Omega_{h}$ that should vanish on
$\Gamma_{w}$ (indeed $w=0$ on $\Gamma_{w}$ implies $u=0$ on $\Gamma_{w}$
thanks to the choice of $\chi$; in fact, the reason for $\chi$ is to decouple
the boundary conditions on $\Gamma$ and $\Gamma_{w}$ from one another).

The test functions $v$ in (2) can be represented in the same way as the
solution (3):

$v=\phi s+\chi(V+\omega\times r)$ (4)

for all vector-valued functions $s$ on $\Omega_{h}$ vanishing on $\Gamma_{w}$,
and $V\in\mathbb{R}^{d}$, $\omega\in\mathbb{R}^{d^{\prime}}$. In particular,
the test functions of the form $\chi(V+\omega\times r)$ can be used to take
into account the force balance (1e–1f). To this end, we introduce
$B_{h}=\Omega_{h}\setminus\Omega$, i.e. the strip between $\Gamma$ and
$G_{h}$, and use the divergence theorem on $B_{h}$ to transfer the boundary
term in (2) from $G_{h}$ to $\Gamma$ where it can be evaluated by (1e–1f):

$\int_{G_{h}}(2\nu D(u)-pI)n\cdot\chi(V+\omega\times r)=\int_{G_{h}}(2\nu
D(u)-pI)n\cdot(V+\omega\times r)\\\ =\int_{\Gamma}(2\nu
D(u)-pI)n\cdot(V+\omega\times r)+\int_{B_{h}}\text{div}(2\nu
D(u)-pI)n\cdot(V+\omega\times r)\\\ =mg\cdot
V-\int_{B_{h}}\rho_{f}g\cdot\chi(V+\omega\times r)=\cdots$ (5)

(the unit normal $n$ on $G_{h}$ in the first line is exterior with respect to
domain $\Omega_{h}$, whereas $n$ on $\Gamma$ in the second line is the
exterior unit normal with respect to domain $\Omega$, so that the exterior
normals with respect to $B_{h}$ are $n$ on $G_{h}$ and $-n$ on $\Gamma$). We
now remark $B_{h}=\Omega_{h}\setminus(\mathcal{O}\setminus\mathcal{S})$ to
rewrite the above as

$\cdots=-\int_{\Omega_{h}}\rho_{f}g\cdot\chi(V+\omega\times
r)+\int_{\mathcal{O}}\rho_{f}g\cdot\chi(V+\omega\times
r)-\int_{\mathcal{S}}\rho_{f}g\cdot(V+\omega\times r)+mg\cdot V\\\
=-\int_{\Omega_{h}}\rho_{f}g\cdot\chi(V+\omega\times
r)+\int_{\mathcal{O}}\rho_{f}g\cdot\chi(V+\omega\times
r)+\left(1-\frac{\rho_{f}}{\rho_{s}}\right)mg\cdot V.$ (6)

The last line is justified by observing
$\int_{\mathcal{S}}\rho_{s}g\cdot(V+\omega\times r)=mg\cdot V$ with $\rho_{s}$
being the constant density of the solid. Indeed,
$\int_{\mathcal{S}}\rho_{s}=m$ and $\int_{\mathcal{S}}\rho_{s}r=0$ since
$r=x-x_{b}$ is the vector pointing from the barycenter of the solid
$x_{b}=\frac{1}{m}\int_{\mathcal{S}}\rho_{s}x$ to the current position $x$.

Substituting the ansatzes (3)-(4) for $u$ and $v$ into (2) and rewriting the
boundary term using (5)-(6) we arrive at the following formal variational
formulation of our problem in terms of the new unknowns
$w,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}U},\psi$:
find $w:\Omega_{h}\to\mathbb{R}^{d}$ vanishing on $\Gamma_{w}$,
$U\in\mathbb{R}^{d}$, $\psi\in\mathbb{R}^{d^{\prime}}$, and
$p:\Omega_{h}\to\mathbb{R}$ such that

$2\nu\int_{\Omega_{h}}D(\phi w+\chi(U+\psi\times r)):D(\phi
s+\chi(V+\omega\times r))-\int_{G_{h}}(2\nu D(\phi w+\chi(U+\psi\times
r))-pI)n\cdot\phi s\\\ -\int_{\Omega_{h}}p\operatorname{div}(\phi
s+\chi(V+\omega\times r))-\int_{\Omega_{h}}q\operatorname{div}(\phi
w+\chi(U+\psi\times r))\\\ =\int_{\Omega_{h}}\rho_{f}g\cdot\phi
s+\int_{\mathcal{O}}\rho_{f}g\cdot\chi(V+\omega\times
r)+\left(1-\frac{\rho_{f}}{\rho_{s}}\right)mg\cdot V$ (7)

for all $s:\Omega_{h}\to\mathbb{R}^{d}$ vanishing on $\Gamma_{w}$,
$V\in\mathbb{R}^{d}$, $\omega\in\mathbb{R}^{d^{\prime}}$, and
$q:\Omega_{h}\to\mathbb{R}$. In addition, the pressure $p$ should satisfy the
constraint (1g).

Note that the formulation above contains only the integrals on $\Omega_{h}$,
$G_{h}$, $\mathcal{O}$ which can be easily approximated by quadrature rules on
meshes $\mathcal{T}_{h}$ and $\mathcal{T}_{h}^{\mathcal{O}}$. We can thus
discretize using the usual finite elements for the trial and test functions.

###  2.3 The $\phi$-FEM scheme: discretization with Taylor-Hood finite
elements

We fix an integer $k\geqslant 2$ and introduce the approximations $\phi_{h}$
and $\chi_{h}$ to the levelset $\phi$ and to the cut-off $\chi$, given by the
standard nodal interpolation to the continuous FE spaces of degree $k$ on the
mesh $\mathcal{T}_{h}^{\mathcal{O}}$. The active computational mesh
$\mathcal{T}_{h}$, its domain $\Omega_{h}$ and the internal boundary component
$G_{h}$ are actually defined as follows, cf. Fig. 1,

$\mathcal{T}_{h}=\\{T\in\mathcal{T}_{h}^{\mathcal{O}}:T\cap\\{\phi_{h}<0\\}\neq\varnothing\\},\quad\Omega_{h}:=\left(\cup_{T\in\mathcal{T}_{h}}T\right)^{o}\,,$
(8)
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}G_{h}=\partial\Omega_{h}\setminus\Gamma_{w}=\\{E\text{
(boundary facets of }\mathcal{T}_{h})\text{ such that
}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi_{h}\geq
0\text{ on }E}\\}\,.}$ (9)

Moreover, we shall need the collections of the mesh cells
$\mathcal{T}_{h}^{\Gamma}$ and facets $\mathcal{F}_{h}^{\Gamma}$ near the
boundary $\Gamma$, as illustrated in Fig. 2, to include the appropriate
stabilization into the FE scheme. More specifically, we introduce the submesh
$\mathcal{T}_{h}^{\Gamma}\subset\mathcal{T}_{h}$ and the corresponding
subdomain $\Omega_{h}^{\Gamma}\subset\Omega$ containing the mesh elements
intersected by the approximate interface

$\Gamma_{h}=\\{\phi_{h}=0\\},$

i.e.

$\mathcal{T}_{h}^{\Gamma}=\\{T\in\mathcal{T}_{h}:T\cap\Gamma_{h}\neq\varnothing\\},\quad\Omega_{h}^{\Gamma}:=\left(\cup_{T\in\mathcal{T}_{h}^{\Gamma}}T\right)^{o}\,.$
(10)

Finally, we set $\mathcal{F}_{h}^{\Gamma}$ as the collection of the interior
facets of the mesh $\mathcal{T}_{h}$ either cut by $\Gamma_{h}$ or belonging
to a cut mesh element

$\mathcal{F}_{h}^{\Gamma}=\\{E\text{ (an internal facet of
}\mathcal{T}_{h})\text{ such that }\exists T\in\mathcal{T}_{h}^{\Gamma}\text{
and }E\in\partial T\\}.$

Figure 2: Example of $\mathcal{T}_{h}^{\Gamma}$ and $\mathcal{F}_{h}^{\Gamma}$
with the geometry given in Fig. 1 (left).

######  Remark 1.

The definitions of $\mathcal{T}_{h}$ and $\mathcal{T}_{h}^{\Gamma}$ above
assume an idealized setting where one can check the sign of $\phi_{h}$ at any
point of any given mesh cell $T$. In practice, one would rather check this
sign only at the vertices of the mesh or, eventually, at some other well
chosen points. In our current implementation, we attribute the cells to
$\mathcal{T}_{h}$ or $\mathcal{T}_{h}^{\Gamma}$ based on the sign of $\phi$
(equivalently, the sign of $\phi_{h}$) at the vertices only: in 2D, a triangle
$T\in\mathcal{T}_{h}^{\mathcal{O}}$ is selected to be in $\mathcal{T}_{h}$ if
$\phi\leqslant 0$ on at least one vertex of $T$; a triangle
$T\in\mathcal{T}_{h}$ is then selected to be in $\mathcal{T}_{h}^{\Gamma}$ if
$\phi\geqslant 0$ on at least one vertex of $T$. This deviation from
definitions (8)-(9) is not covered by our theory.

We also note that in more advanced applications of $\phi$-FEM, $\phi_{h}$ and
$\chi_{h}$ may be given directly on the discrete level, for instance by a
discrete level-set equation. This possibility is however outside of the scope
of the present article. We suppose here that the shape of the particle is
sufficiently simple so that $\phi$ and $\chi$ are known analytically.

Introduce the FE spaces for velocity and pressure on the mesh
$\mathcal{T}_{h}$:

$\mathcal{V}_{h}=\left\\{v_{h}\in
C(\bar{\Omega})^{d}:v_{h}|_{T}\in\mathbb{P}^{k}(T)^{d}\hskip 10.00002pt\forall
T\in\mathcal{T}_{h},\hskip 10.00002ptv_{h}=0\text{ on }\Gamma_{w}\right\\}$

and

$\mathcal{M}_{h}=\left\\{q_{h}\in
C(\bar{\Omega}):q_{h}|_{T}\in\mathbb{P}^{k-1}(T)\hskip 10.00002pt\forall
T\in\mathcal{T}_{h},\hskip 10.00002pt\int_{\Omega}q_{h}=0\right\\}.$

######  Remark 2.

Note that the definition of the pressure space involves an integral on
$\Omega$, which is incompatible with our $\phi$-FEM framework since its whole
point is to avoid integrals on $\Omega$ and $\Gamma$. In practice, we shall
rather impose $\int_{\Omega_{h}}q_{h}=0$, introducing a mismatch in the
additive pressure constant (which, anyway, has no physical meaning) with
respect to the exact solution satisfying (1g). We prefer however to keep the
unimplementable constraint in the definition above to avoid some technical
difficulties in theory. In practice, a special care will have to be taken in
the interpretation of the error in pressure. We shall return to this technical
point in the numerical results section.

The stabilized scheme inspired by (7) can be now written as: find
$w_{h}\in\mathcal{V}_{h}$,
$U_{h}\in\mathbb{R}^{d},\psi_{h}\in\mathbb{R}^{d^{\prime}}$,
$p_{h}\in\mathcal{M}_{h}$ such that

$2\nu\int_{\Omega_{h}}D(\phi_{h}w_{h}+\chi_{h}(U_{h}+\psi_{h}\times
r)):D(\phi_{h}s_{h}+\chi_{h}(V_{h}+\omega_{h}\times r))\\\ -\int_{G_{h}}(2\nu
D(\chi_{h}(U_{h}+\psi_{h}\times
r)+\phi_{h}w_{h})-p_{h}I)n\cdot\phi_{h}s_{h}-\int_{\Omega_{h}}p_{h}\operatorname{div}(\phi_{h}s_{h}+\chi_{h}(V_{h}+\omega_{h}\times
r))\\\
-\int_{\Omega_{h}}q_{h}\operatorname{div}(\phi_{h}w_{h}+\chi_{h}(U_{h}+\psi_{h}\times
r))\\\ +\sigma_{u}J_{u}(\chi_{h}(U_{h}+\psi_{h}\times
r)+\phi_{h}w_{h},\chi_{h}(V_{h}+\omega_{h}\times r)+\phi_{h}s_{h})\\\ +\sigma
h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}(-\nu\Delta(\phi_{h}w_{h}+\chi_{h}(U_{h}+\psi_{h}\times
r))+\nabla
p_{h})\cdot(-\nu\Delta(\phi_{h}s_{h}+\chi_{h}(V_{h}+\omega_{h}\times
r))-\nabla q_{h})\\\
+\sigma\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}\operatorname{div}(\phi_{h}w_{h}+\chi_{h}(U_{h}+\psi_{h}\times
r))\operatorname{div}(\phi_{h}s_{h}+\chi_{h}(V_{h}+\omega_{h}\times r))\\\
=\int_{\Omega_{h}}\rho_{f}g\cdot\phi_{h}s_{h}+\int_{\mathcal{O}}\rho_{f}g\cdot\chi_{h}(V_{h}+\omega_{h}\times
r)+\left(1-\frac{\rho_{f}}{\rho_{s}}\right)mg\cdot V_{h}\\\ +\sigma
h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}\rho_{f}g\cdot(-\nu\Delta(\phi_{h}s_{h}+\chi_{h}(V_{h}+\omega_{h}\times
r))-\nabla q_{h})$ (11)

for all $s_{h}\in\mathcal{V}_{h}$,
$V_{h}\in\mathbb{R}^{d},\omega_{h}\in\mathbb{R}^{d^{\prime}}$,
$q_{h}\in\mathcal{M}_{h}$.

Here $J_{u}$ is the ghost penalties for the velocity, cf. [7]:

$\displaystyle J_{u}(u,v)$
$\displaystyle=h\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}\left[{\partial_{n}}u\right]\cdot\left[{\partial_{n}}v\right]+h^{3}\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}\left[{\partial_{n}^{2}}u\right]\cdot\left[{\partial_{n}^{2}}v\right]\,.$

Note that, unlike [7, 8], we do not penalize the jumps of all the derivatives
of the velocity; only the derivatives of order up to $2$ are included in
$J_{u}$. There is no penalization on the pressure either. This alleviation of
the ghost penalty is possible thanks to the additional least-squares-type
stabilization (the terms multiplied by $\sigma$), cf. Lemmas 1 and 2. These
least-squares terms are also necessary in themselves to control the fictitious
extension of the solution outside $\Omega$, cf. the proof of Lemma 9. Note
that this extension is not present in CutFEM (this is indeed the principal
difference between CutFEM and $\phi$-FEM). We also mention that the version of
the ghost penalty in $\phi$-FEM for Poisson problem in [12] is even more
reduced: only the jumps of the first order derivatives are penalized there.
The inclusion of the second order derivatives in $J_{u}$ in the present case
of Stokes equations allows us to control both velocity and pressure in the
forthcoming proofs, cf. Lemma 1.

###  2.4 Assumptions on the mesh and main results

Prior to stating our main results on the numerical convergence of our method,
we begin with some geometrical assumptions on $\Omega$ and the functions
$\phi$ and $\chi$.

######  Assumption 1.

The boundary $\Gamma$ can be covered by open sets $\mathcal{O}_{i}$,
$i=1,\ldots,I$ and one can introduce on every $\mathcal{O}_{i}$ local
coordinates $\xi_{1},\ldots,\xi_{d}$ with $\xi_{d}=\phi$ such that all the
partial derivatives $\partial^{\alpha}\xi/\partial x^{\alpha}$ and
$\partial^{\alpha}x/\partial\xi^{\alpha}$ up to order $k+1$ are bounded by
some $C_{0}>0$. Thus, $\phi$ is of class $C^{k+1}$ on $\mathcal{O}$.

######  Assumption 2.

$\chi\in H^{k+1}(\mathcal{O})$, $\chi=1$ on $\mathcal{S}$, $\chi=0$ on
$\Gamma_{w}$.

We continue with assumptions on the mesh. To this end, we introduce an
extended band of mesh elements near the boundary $\Gamma$, namely the submesh
$\mathcal{T}_{h}^{\Gamma,ext}$ with
$\mathcal{T}_{h}^{\Gamma}\subset\mathcal{T}_{h}^{\Gamma,ext}\subset\mathcal{T}_{h}$
by adding to $\mathcal{T}_{h}^{\Gamma}$ the cells which are neighbors and
neighbors of neighbors of cells in $\mathcal{T}_{h}^{\Gamma}$.

######  Assumption 3.

$|\nabla\phi|\geq m$, $|\nabla\phi_{h}|\geq\frac{m}{2}$ on all the mesh cells
in $\mathcal{T}_{h}^{\Gamma,ext}$, $|\phi|\geq mh$ on
$\mathcal{T}_{h}\setminus\mathcal{T}_{h}^{\Gamma,ext}$, and
$|\nabla\phi_{h}|\leq M$ on $\Omega_{h}$ with some $m,M>0$.

######  Assumption 4.

The approximate interface $G_{h}$ can be covered by element patches
$\\{\Pi_{k}\\}_{k=1,\ldots,N_{\Pi}}$ having the following properties:

  * •

Each $\Pi_{k}$ is composed of a mesh element $T_{k}$ lying inside $\Omega$ and
some elements cut by $\Gamma$, more precisely
$\Pi_{k}=T_{k}\cup\Pi_{k}^{\Gamma}$ where $T_{k}\in\mathcal{T}_{h}$,
$T_{k}\subset\bar{\Omega}$, $\Pi_{k}^{\Gamma}\subset\mathcal{T}_{h}^{\Gamma}$,
and $\Pi_{k}^{\Gamma}$ contains at most $N$ mesh elements;

  * •

Each mesh element in a patch $\Pi_{k}$ shares at least a facet with another
mesh element in the same patch. In particular, $T_{k}$ shares a facet $F_{k}$
with an element in $\Pi_{k}^{\Gamma}$;

  * •

$\mathcal{T}_{h}^{\Gamma}=\cup_{k=1}^{N_{\Pi}}\Pi_{k}^{\Gamma}$;

  * •

$\Pi_{k}$ and $\Pi_{l}$ are disjoint if $k\neq l$.

######  Assumption 5.

Any mesh cell $T\in\mathcal{T}_{h}$ has at least $d$ facets not lying on
$\Gamma_{w}$.

######  Remark 3.

Assumptions 1, 3, 4 are similar to those made in the previous $\phi$-FEM
publications [12, 11], which contain a more detailed discussion about them and
some illustrations. In what concerns the mesh, these assumptions are satisfied
if the mesh is sufficiently refined and $\Gamma$ is sufficiently smooth. In
what concerning the level-set function $\phi$, we require essentially that it
behaves like the signed distance to $\Gamma$ near $\Gamma$ and it is bounded
away from zero far from $\Gamma$, while remaining globally smooth. In general,
one cannot thus take the signed distance to $\Gamma$ as $\phi$ everywhere
since it is guaranteed to be smooth only in a vicinity of $\Gamma$.

######  Remark 4.

The last assumption 5 is usually required in the theoretical analysis of
Taylor-Hood elements for the Stokes equation in the geometrically conforming
setting [13], although it can be significantly relaxed, at least in the 2D
setting [5]. Note that this assumption only affects the mesh near the outer
wall $\Gamma_{w}$, more particularly in the corners of $\mathcal{O}$, which we
treat in the standard geometrically conforming manner anyway. It does not
impose any further restriction on the active mesh $\mathcal{T}_{h}$ near the
interface $\Gamma$, where $\phi$-FEM is effectively employed.

Let us now state our main results:

######  Theorem 1.

Suppose that Assumptions 1–5 hold true and the mesh $\mathcal{T}_{h}$ is
quasi-uniform. Let $(u,U,\psi,p)\in
H^{k+1}(\Omega)^{d}\times\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}\times
H^{k}(\Omega)$ be the solution to (1a)-(1g) and
$(w_{h},U_{h},\psi_{h},p_{h})\in\mathcal{V}_{h}\times\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}\times\mathcal{M}_{h}$
be the solution to (11). Denoting

$u_{h}:=\chi_{h}(U_{h}+\psi_{h}\times r)+\phi_{h}w_{h}$

it holds for $h\leq h_{0}$

$|u-u_{h}|_{1,\Omega\cap\Omega_{h}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{\nu}}\|p-p_{h}\|_{0,\Omega\cap\Omega_{h}}\leq
Ch^{k}(\|u\|_{k+1,\Omega}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{\nu}}\|p\|_{k,\Omega})$
(12)

and

$|U-U_{h}|+|\psi-\psi_{h}|\leq
Ch^{k}(\|u\|_{k+1,\Omega}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{\nu}}\|p\|_{k,\Omega})$
(13)

with some $C>0$ and $h_{0}>0$ depending on the parameters $C_{0}$, $m$, $M$,
$N$ in Assumptions 1–5, on the maximum of the derivatives of $\phi$ and $\chi$
of order up to $k+1$, on the mesh regularity, and on the polynomial degree
$k$, but independent of $h$, $f$, and $u$.

Moreover, supposing $\Omega\subset\Omega_{h}$

$\|u-u_{h}\|_{0,\Omega}\leq
Ch^{k+1/2}(\|u\|_{k+1,\Omega}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{\nu}}\|p\|_{k,\Omega})$
(14)

with a constant $C>0$ of the same type as above.

######  Remark 5.

The numerical results in Section 4 suggest that the convergence order for the
particle velocity error and that for the $L^{2}$-error of the fluid velocity
is $k+1$. This suggests that both estimates (13) and (14) are not sharp. This
is similar to our previous studies [12] and [11], in what concerns the
$L^{2}$-error.

##  3 Proofs

From now on, we put the viscosity of the fluid to $\nu=1$ to simplify the
formulas. The general case can be easily recovered by dividing the governing
equations by $\nu$ and redefining $\frac{p}{\nu}$ as $p$. We shall also
tacitly suppose that Assumptions 1–5 holds true.

This section is organized as follows: we start with some technical lemmas in
Sections 3.1 and 3.2, essentially adapting the corresponding results from
[12]. Note however that the interpolation error bound in Section 3.2 is
sharper than its counterpart in [12]; it is now optimal with respect to the
Sobolev norm expected from the interpolated function. The proofs of Taylor-
Hood inf-sup stability, the generalized coercivity of the bilinear form and
finally the a priori error estimates are then given, respectively, in Sections
3.3, 3.4 and 3.5, thus establishing Theorem 1.

###  3.1 Some technical lemmas.

Lemmas 1, 2, and 4 are adaptions of, respectively, Lemmas 3.2, 3.3, and
3.4-3.5 from [12]. Lemma 3 is a version of the well known Korn inequality
stating the uniformity of the constant in this inequality for a family of
$h$-dependent domains $\Omega_{h}$. All these results, most notably Lemma 2,
are necessary to prove the inf-sup stability of our scheme.

######  Lemma 1.

Let $T$ be a triangle/tetrahedron, $E$ one of its sides, $v$ a vector-valued
polynomial function on $T$, $q$ a scalar-valued polynomial function on $T$
such that

$v=\frac{\partial v}{\partial n}=\frac{\partial^{2}v}{\partial n^{2}}=0\text{
on }E$ (15)

and

$-\Delta v+\nabla q=0\mbox{ and }\operatorname{div}v=0\text{ on }T.$ (16)

Then $v=0$ and $q=const$ on $T$.

###### Proof.

We shall give the proof only in the two dimensional setting, the
generalization to the case $d=3$ being straightforward. Without loss of
generality, we can choose the Cartesian coordinates $(x,y)$ such that the edge
$E$ lies on the $x$-axis. We shall denote the components of the vector-valued
function $v$ by $(v_{1},v_{2})$.

Let us write, for $k\in\\{1,2\\}$,

$v_{k}=\sum\limits_{\scriptstyle i,j\geq 0}v_{ij}^{k}x^{i}y^{j}\mbox{ and
}q=\sum\limits_{\scriptstyle i,j\geq 0}q_{ij}x^{i}y^{j}.$

We will prove by strong induction on $m$ that

$v^{k}_{ij}=0\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\forall k\in\\{1,2\\},\leavevmode\nobreak\ i\geq 0,\leavevmode\nobreak\
j\in\\{0,...,m\\}.$ ($S_{m}$)

Using (15), it holds

$\sum\limits_{i\geq 0}v_{i0}^{k}x^{i}=\sum\limits_{i\geq
0}v_{i1}^{k}x^{i}=\sum\limits_{i\geq 0}v_{i2}^{k}x^{i}=0$

for all $x$ and $k\in\\{1,2\\}$. Hence $v_{i0}^{k}=v_{i1}^{k}=v_{i2}^{k}=0$
for all $i\geq 0$, $k\in\\{1,2\\}$ and we obtain ($S_{2}$). Assume that for a
given $m\geq 2$, ($S_{m}$) holds true. Thanks to (16), one has

$-(i+2)(i+1)v_{(i+2)j}^{1}-(j+2)(j+1)v_{i(j+2)}^{1}+(i+1)q_{(i+1)j}=0,$ (17)
$-(i+2)(i+1)v_{(i+2)j}^{2}-(j+2)(j+1)v_{i(j+2)}^{2}+(j+1)q_{i(j+1)}=0$ (18)

and

$(i+1)v_{(i+1)j}^{1}+(j+1)v_{i(j+1)}^{2}=0.$ (19)

From (17) and (18), for each $i,j$

$(i+2)v_{(i+2)(j+1)}^{1}+\frac{(j+3)(j+2)}{i+1}v_{i(j+3)}^{1}=\frac{(i+3)(i+2)}{j+1}v_{(i+3)j}^{2}+(j+2)v_{(i+1)(j+2)}^{2}.$

The last equality for $j=m-2$ combined with ($S_{m}$) implies that

$v^{1}_{i(m+1)}=0\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \forall i\geq 0.$

Relation (19) for $j=m$ and ($S_{m}$) gives

$v^{2}_{i(m+1)}=0\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \forall i\geq 0,$

which leads to ($S_{m+1}$). Thus $v=0$. This also implies $\nabla q=0$ on $T$
thanks to (16). ∎

######  Lemma 2.

For any $\beta>0$ and any integers $s,r\geq 1$ there exists $0<\alpha<1$
depending only on the mesh regularity and $s,r$ such that for any continuous
vector-valued $\mathbb{P}^{s}$ FE function $v_{h}$ on $\mathcal{T}_{h}$ and
any continuous scalar $\mathbb{P}^{r}$ FE function $q_{h}$ it holds

$\|D(v_{h})\|_{0,\Omega_{h}^{\Gamma}}^{2}+(1-\alpha)h^{2}|q_{h}|_{1,\Omega_{h}^{\Gamma}}^{2}\leq\alpha\|D(v_{h})\|_{0,\Omega_{h}}^{2}\\\
+\beta\left(h^{2}\|-\Delta v_{h}+\nabla
q_{h}\|_{0,\Omega_{h}^{\Gamma}}^{2}+\|\operatorname{div}v_{h}\|_{0,\Omega_{h}^{\Gamma}}^{2}+\sum_{E\in\mathcal{F}_{h}^{\Gamma}}(h\|[\partial_{n}v_{h}]\|_{0,E}^{2}+h^{3}\|[\partial_{n}^{2}v_{h}]\|_{0,E}^{2})\right).$
(20)

###### Proof.

Thanks to Assumption 4, the boundary $\Gamma$ can be covered by patches
$\\{\Pi_{k}\\}_{k=1,...,N_{\Pi}}$. Take $\beta>0$ and set

$\alpha:=\max\limits_{\Pi_{k},v_{h},q_{h}}F(\Pi_{k},v_{h},q_{h}),$ (21)

where

$F(\Pi_{k},v_{h},q_{h})=\dfrac{\|D(v_{h})\|_{0,\Pi_{k}^{\Gamma}}^{2}{+h^{2}|q_{h}|_{1,\Pi_{k}^{\Gamma}}^{2}}-\beta
G(\Pi_{k},v_{h},q_{h})}{\|D(v_{h})\|_{0,\Pi_{k}}^{2}{+h^{2}|q_{h}|_{1,\Pi_{k}^{\Gamma}}^{2}}}$

with

$G(\Pi_{k},v_{h},q_{h})=h^{2}\|-\Delta v_{h}+\nabla
q_{h}\|_{0,\Pi_{k}^{\Gamma}}^{2}+\|\operatorname{div}v_{h}\|_{0,\Pi_{k}^{\Gamma}}^{2}+h\left\|\left[\partial_{n}v_{h}\right]\right\|_{\mathcal{F}_{k}}^{2}+h^{3}\left\|\left[\partial_{n}^{2}v_{h}\right]\right\|_{\mathcal{F}_{k}}^{2}.$

The maximum in (21) is taken over all the continuous vector-valued
$\mathbb{P}^{s}$ FE functions $v_{h}$ on $\Pi_{k}$, all the continuous scalar
$\mathbb{P}^{r}$ FE functions $q_{h}$ on $\Pi_{k}^{\Gamma}$, such that the
denominator in the expression for $F$ does not vanish, and over all the
possible configurations of patches $\Pi_{k}$ satisfying Assumption 4. The
notation $\mathcal{F}_{k}$ stands for the set of mesh facets inside the patch
$\Pi_{k}$ which includes thus $F_{k}$ separating $T_{k}$ from
$\Pi_{k}^{\Gamma}$ and the other facets inside $\Pi_{k}^{\Gamma}$. The norm
$\|\cdot\|_{\mathcal{F}_{k}}$ should be understood as
$(\sum_{F\in\mathcal{F}_{k}}\|\cdot\|_{F}^{2})^{1/2}$.

Since the maximized function $F$ is invariant with respect to the
transformation $x\mapsto\frac{1}{h}x$, $v_{h}\mapsto\frac{1}{h}v_{h}$,
$q_{h}\mapsto q_{h}$, we can assume that $h=1$ in (21). Furthermore
$F(\Pi_{k},v_{h},q_{h})=F(\Pi_{k},\lambda v_{h},\lambda q_{h})$ for any
$\lambda\neq 0$. Hence the maximum (21) is attained since it can be taken over
all admissible patches with $h=1$ and all $v_{h},q_{h}$ such that
$\|D(v_{h})\|_{0,\Pi_{k}}^{2}{+h^{2}|q_{h}|_{1,\Pi_{k}^{\Gamma}}^{2}}=1$,
forming the unit sphere in the finite dimensional space of all $(v_{h},q_{h})$
factored by rigid body motions on $\Pi_{k}$ and constants on
$\Pi_{k}^{\Gamma}$.

Clearly $\alpha\leqslant 1$. Let us prove by contradiction that $\alpha<1$.
Assume that $\alpha=1$. Consider the patch $\Pi_{k}$ (with $h=1$) and $v_{h}$,
$q_{h}$ with
$\|D(v_{h})\|_{0,\Pi_{k}}^{2}{+|q_{h}|_{1,\Pi_{k}^{\Gamma}}^{2}}=1$ on which
the maximum (21) is attained. Then,

$\|D(v_{h})\|_{0,T_{k}}^{2}+\beta G(\Pi_{k},v_{h},q_{h})=0,$

since $\Pi_{k}=T_{k}\cup\Pi_{k}^{\Gamma}$. We deduce that
$v_{h}=\frac{\partial v_{h}}{\partial n}=\frac{\partial^{2}v_{h}}{\partial
n^{2}}=0$ on all the facets in $\mathcal{F}_{k}$ and $-\Delta v_{h}+\nabla
q_{h}=0$, $\operatorname{div}(v_{h})=0$ on all $T\in\Pi_{k}^{\Gamma}$.
Moreover, $v_{h}$ is a rigid body motion on $T_{k}$. Let $v_{h}^{rbm}$ be the
rigid body motion velocity on $\mathbb{R}^{d}$ coinciding with $v_{h}$ on
$T_{k}$. Thanks to Lemma 1 applied to $v_{h}-v_{h}^{rbm}$ and to $q_{h}$ on
the cells in $\Pi_{k}^{\Gamma}$ starting from the cell adjacent in
$\Pi_{k}^{\Gamma}$ to $T_{k}$, we have $v_{h}=v_{h}^{rbm}$ on $\Pi_{k}$ and
$q_{h}=const$ on $\Pi_{k}^{\Gamma}$ (recall that $q_{h}$ is continuous). We
have thus reached a contradiction with the assumptions
$\|D(v_{h})\|_{0,\Pi_{k}}^{2}{+|q_{h}|_{1,\Pi_{k}^{\Gamma}}^{2}}=1$ and
$\alpha=1$.

This proves that there exists $\alpha<1$ such that

$\|D(v_{h})\|_{0,\Pi_{k}^{\Gamma}}^{2}+(1-\alpha)h^{2}|q_{h}|_{1,\Pi_{k}^{\Gamma}}^{2}\leq\alpha\|D(v_{h})\|_{0,\Pi_{k}}^{2}\\\
+\beta\left(h^{2}\|-\Delta v_{h}+\nabla
q_{h}\|_{0,\Pi_{k}^{\Gamma}}^{2}+\|\operatorname{div}v_{h}\|_{0,\Pi_{k}^{\Gamma}}^{2}+h\left\|\left[\partial_{n}v_{h}\right]\right\|_{\mathcal{F}_{k}}^{2}+h^{3}\left\|\left[\partial_{n}^{2}v_{h}\right]\right\|_{\mathcal{F}_{k}}^{2}\right)$

on all the patches $\Pi_{k}$ and for all $v_{h},q_{h}$. Summing this over all
$\Pi_{k}$ gives (20). ∎

######  Lemma 3.

For any $v\in H^{1}(\Omega)^{d}$ vanishing on $\Gamma_{w}$

$|v|_{1,\Omega_{h}}\leq C\|D(v)\|_{0,\Omega_{h}}\,.$ (22)

###### Proof.

Since $v=0$ on $\Gamma_{\omega}$, we have the following Korn inequality

$|v|_{1,\Omega}\leq C\|D(v)\|_{0,\Omega}$ (23)

with a constant $C>0$ depending only on the shape of $\Omega$, cf. [9, Theorem
6.3-4]. This implies

$|v|_{1,\Omega_{h}^{i}}\leq C\|D(v)\|_{0,\Omega_{h}},$ (24)

where $\Omega_{h}^{i}$ denotes the mesh cells inside $\Omega$. Now, for any
pair of mesh cells $T,T^{\prime}$ sharing a facet $E$, we can prove

$|v|_{1,T}\leq C\|D(v)\|_{0,T}+C|v|_{1,T^{\prime}}$ (25)

with a constant $C$ independent of $h$. Indeed, combining the Korn
inequalities (23) and the trace theorem on the reference element leads to
$|v|_{1,T}\leq C\|D(v)\|_{0,T}+C|v|_{1/2,E}$. Employing again the trace
inequality $|v|_{1/2,E}\leq C|v|_{1,T^{\prime}}$ (see [6, Lemma 7.5.26 ])
leads to (25).

Let $\Omega_{h}^{i,1}$ be $\Omega_{h}^{i}$ plus the cells which are not in
$\Omega_{h}^{i}$ but have a neighbor in $\Omega_{h}^{i}$. For any such cell
$T$, we take $T^{\prime}$ as its neighbor in $\Omega_{h}^{i}$, apply the
estimate above and sum which gives

$|v|_{1,\Omega_{h}^{i,1}\backslash\Omega_{h}^{i}}\leq
C\|D(v)\|_{0,\Omega_{h}^{i,1}\backslash\Omega_{h}^{i}}+C|v|_{1,\Omega_{h}^{i}},$

hence, using (24),

$|v|_{1,\Omega_{h}^{i,1}}\leq C\|D(v)\|_{0,\Omega_{h}}.$

Let $\Omega_{h}^{i,2}$ be $\Omega_{h}^{i,1}$ plus the cells which are not in
$\Omega_{h}^{i,1}$ but have a neighbor in $\Omega_{h}^{i,1}$. We have similar
to above

$|v|_{1,\Omega_{h}^{i,2}}\leq C\|D(v)\|_{0,\Omega_{h}}$

and so on. After a finite number of steps, say $k$, we arrive at
$\Omega_{h}^{i,k}=\Omega_{h}$. And

$|v|_{1,\Omega_{h}}=|v|_{1,\Omega_{h}^{i,k}}\leq C\|D(v)\|_{0,\Omega_{h}}.$

∎

######  Lemma 4.

For any $s_{h}\in\mathcal{V}_{h}$ and any
$V_{h}\in\mathbb{R}^{d},\omega_{h}\in\mathbb{R}^{d^{\prime}}$,

$\|\phi_{h}s_{h}\|_{0,\Omega_{h}^{\Gamma}}+\left(\sum_{F\in\mathcal{F}_{h}^{\Gamma}}h\|\phi_{h}s_{h}\|_{0,F}^{2}\right)^{1/2}+\sqrt{h}\|\phi_{h}s_{h}\|_{0,G_{h}}\leq
Ch\|D(\phi_{h}s_{h}+\chi_{h}(V_{h}+\omega_{h}\times r))\|_{0,\Omega_{h}}$

and

$|V_{h}|+|\omega_{h}|\leqslant
C\|\phi_{h}s_{h}+\chi_{h}(V_{h}+\omega_{h}\times r)\|_{1,\Omega_{h}}\,.$ (26)

###### Proof.

Take any $s_{h}\in\mathcal{V}_{h},$
$V_{h}\in\mathbb{R}^{d},\leavevmode\nobreak\
\omega_{h}\in\mathbb{R}^{d^{\prime}}$ and denote
$v_{h}=\phi_{h}s_{h}+\chi_{h}(V_{h}+\omega_{h}\times r)$. By [12, Lemma 3.4]

$\|\phi_{h}s_{h}\|_{0,\Omega_{h}^{\Gamma}}\leqslant
Ch|\phi_{h}s_{h}|_{1,\Omega_{h}^{\Gamma}}\leqslant
Ch(|v_{h}|_{1,\Omega_{h}^{\Gamma}}+|\chi_{h}(V_{h}+\omega_{h}\times
r)|_{1,\Omega_{h}^{\Gamma}})\\\ \leqslant
Ch(|v_{h}|_{1,\Omega_{h}^{\Gamma}}+\|\chi_{h}\|_{1,\Omega_{h}^{\Gamma}}(|V_{h}|+|\omega_{h}|)).$
(27)

By equivalence of norms

$|V_{h}|+|\omega_{h}|\leqslant C\|V_{h}+\omega_{h}\times r\|_{0,\Gamma}.$ (28)

Denote by $B_{h}^{\Gamma}$ the band between $\Gamma$ and $\Gamma_{h}$.
Applying the divergence theorem to the vector field $|V_{h}+\omega_{h}\times
r|^{2}{\nabla\phi}$ and noting that the normal on $\Gamma$ (resp.
$\Gamma_{h}$) is given by $\pm\frac{\nabla\phi}{|\nabla\phi|}$ (resp.
$\pm\frac{\nabla\phi_{h}}{|\nabla\phi_{h}|}$) gives

$\int_{\Gamma}|V_{h}+\omega_{h}\times
r|^{2}|\nabla\phi|\leqslant\int_{\Gamma_{h}}|V_{h}+\omega_{h}\times
r|^{2}\frac{|{\nabla\phi}\cdot\nabla\phi_{h}|}{|\nabla\phi_{h}|}+\left|\int_{B_{h}^{\Gamma}}\operatorname{div}\left(|V_{h}+\omega_{h}\times
r|^{2}{\nabla\phi}\right)\right|$

We now note that $|\nabla\phi|$ (resp. $|\nabla\phi_{h}|$ and $\chi_{h}$) are
both positive and bounded away from 0 on $\Gamma$ (resp. on $\Gamma_{h}$)
uniformly in $h$ for $h$ small enough, and the measure of $B_{h}^{\Gamma}$ is
of order $h^{k+1}$. The inequality above implies thus

$\|V_{h}+\omega_{h}\times r\|_{0,\Gamma}^{2}\leq
C(\|\chi_{h}(V_{h}+\omega_{h}\times
r)\|^{2}_{0,\Gamma_{h}}+h^{k+1}(|V_{h}|+|\omega_{h}|)^{2}).$

Combining this with (28) gives, for $h$ small enough

$|V_{h}|+|\omega_{h}|\leqslant C\|\chi_{h}(V_{h}+\omega_{h}\times
r)\|_{0,\Gamma_{h}}=C\|v_{h}\|_{0,\Gamma_{h}}$

and, by the trace inequality, $|V_{h}|+|\omega_{h}|\leqslant
C\|v_{h}\|_{1,\Omega_{h}}$, i.e. (26). Substituting this into (27) and
combining with the Korn inequality (22) yields the announced estimate for
$\|\phi_{h}s_{h}\|_{0,\Omega_{h}^{\Gamma}}$ since
$\|\chi_{h}\|_{1,\Omega_{h}^{\Gamma}}$ is bounded uniformly in $h$. The
remaining part of the estimate follows by trace inverse inequalities as in
[12, Lemma 3.5]. ∎

###  3.2 Interpolation by finite elements multiplied with the level set

We recall first a Hardy-type inequality, cf. [12].

######  Lemma 5.

For any integer $s\in[0,k]$ and any $u\in H^{s+1}(\Omega_{h})$ vanishing on
$\Gamma$, it holds
$\left\|\displaystyle{\frac{u}{\phi}}\right\|_{s,\Omega_{h}}\leq
C\|u\|_{s+1,\Omega_{h}}$ with $C>0$ depending only on the constants in
Assumption 1 and on $s$.

This allows us to prove the following bound for interpolation by the products
of finite elements with $\phi_{h}$.

######  Lemma 6.

Let $t$ be an integer $1\leqslant t\leqslant k+1$. For any $v\in
H^{t}(\Omega_{h})^{d}\cap H^{1}_{0}(\Omega)$ there exists
$w_{h}\in\mathcal{V}_{h}$ s.t.

$\|v-\phi_{h}w_{h}\|_{s,\Omega_{h}}\leqslant
Ch^{t-s}\|v\|_{t,\Omega_{h}},\quad s=0,1$ (29)

with $C>0$ depending only on $t,s$, the constants in Assumptions 1–3, and the
mesh regularity.

###### Proof.

Let $v\in H^{t}(\Omega_{h})^{d}$, $v=0$ on $\Gamma$, and set $w=v/\phi$.
Thanks to Lemma 5, $w\in H^{t-1}(\Omega_{h})^{d}$ and
$\left\|w\right\|_{t-1,\Omega_{h}}\leq C\|v\|_{t,\Omega_{h}}$. Consider
$w_{h}=I_{h}^{c}w$, where $I_{h}^{c}$ is a Scott-Zhang interpolation operator.
For any $T\in\mathcal{T}_{h}$, let $\omega_{T}$ denote the patch of mesh cells
adjacent to $T$ (not necessarily all the adjacent cells) regrouping the cells
affected by the construction of $I_{h}^{c}w$ on $T$, so that $I_{h}^{c}w$ on
$T$ depends on $w$ only through its restriction to $\omega_{T}$. The Scott-
Zhang interpolation operator can be constructed so that
$\omega_{T}\subset\mathcal{T}_{h}^{\Gamma,ext}$ for all
$T\in\mathcal{T}_{h}^{\Gamma,ext}$, and
$\omega_{T}\subset\mathcal{T}_{h}\setminus\mathcal{T}_{h}^{\Gamma,ext}$ for
all $T\in\mathcal{T}_{h}\setminus\mathcal{T}_{h}^{\Gamma,ext}$. In what
follows, we assume that the operator $I_{h}^{c}$ enjoys this property together
with the usual interpolation error estimates see for instance [6].

Our first goal is to prove (29) for $s=0$. Taking any $T\in\mathcal{T}_{h}$.
Recall that $\phi$ is supposed to be of class (at least) $C^{t}$ so that
$\|\phi-\phi_{h}\|_{\infty,T}\leqslant Ch^{t}$. Hence,

$\displaystyle\|v-\phi_{h}w_{h}\|_{0,T}$ $\displaystyle=\|\phi
w-\phi_{h}w_{h}\|_{0,T}$ (30)
$\displaystyle\leqslant\|\phi\|_{\infty,T}\|w-w_{h}\|_{0,T}+\|\phi-\phi_{h}\|_{\infty,T}\|w_{h}\|_{0,T}$
$\displaystyle\leqslant\|\phi\|_{\infty,T}\|w-w_{h}\|_{0,T}+Ch^{t}\|w\|_{0,\omega_{T}}.$

To continue this proof, we distinguish two cases: the cells
$T\in\mathcal{T}_{h}^{\Gamma,ext}$ close to $\Gamma_{h}$ and the remaining
cells, which are at the distance of at least order $h$ from $\Gamma_{h}$.

  * (i)

Consider $T\in\mathcal{T}_{h}^{\Gamma,ext}$. We have
$\|\phi\|_{\infty,T}\leqslant Ch$ on these cells since they are at the
distance $\sim h$ from $\Gamma$. Noting that $\|w-w_{h}\|_{0,T}\leqslant
Ch^{t-1}|w|_{t-1,\omega_{T}}$ by the usual interpolation estimate, we derive
from (30)

$\|v-\phi_{h}w_{h}\|_{0,T}\leqslant
Ch^{t}\left(|w|_{t-1,\omega_{T}}+\|w\|_{0,\omega_{T}}\right)\,.$ (31)

  * (ii)

Now consider $T\in\mathcal{T}_{h}\setminus\mathcal{T}_{h}^{\Gamma,ext}$. We
note that $\phi$ does not vanish on $\omega_{T}$ for such $T$ (recall that
$\omega_{T}\subset\mathcal{T}_{h}\setminus\mathcal{T}_{h}^{\Gamma,ext}$), so
that $w\in H^{t}(\omega_{T})$ and, by (30) and the usual approximation
estimates,

$\|v-\phi_{h}w_{h}\|_{0,T}\leqslant
Ch^{t}\left(\|\phi\|_{\infty,T}|w|_{t,\omega_{T}}+\|w\|_{0,\omega_{T}}\right)\,.$
(32)

In order to bound $|w|_{t,\omega_{T}}$ here, we recall the Leibniz rule valid
for any multi-index $\alpha\in\mathbb{N}^{d}$

$\partial^{\alpha}v=\partial^{\alpha}(\phi
w)=\sum_{\scriptsize{\begin{array}[]{c}\beta\in\mathbb{N}^{d}\\\
\beta\leqslant\alpha\end{array}}}C_{\alpha}^{\beta}\,(\partial^{\beta}\phi)\,(\partial^{\alpha-\beta}w)$

with binomial coefficients $C_{\alpha}^{\beta}$ depending only on the multi-
indices $\alpha$ and $\beta$ (this formula can be easily proven by induction
on the length of $\alpha$) $\beta\leqslant\alpha$ means $\forall
i=1,\dots,d,\beta_{i}\leqslant\alpha_{i}$. If $\alpha\not=0$, this can be
rewritten, by separating the term with $\beta=0$ (note that
$C^{0}_{\alpha}=1$) and dividing by $\phi$, as

$\partial^{\alpha}w=\frac{1}{\phi}\partial^{\alpha}v\quad-\sum_{\scriptsize{\begin{array}[]{c}\beta\in\mathbb{N}^{d}\\\
\beta\leqslant\alpha,\beta\neq
0\end{array}}}C_{\alpha}^{\beta}\frac{\partial^{\beta}\phi}{\phi}\partial^{\alpha-\beta}w$
(33)

Applying (33) to $w$ on $\omega_{T}$ gives

$|w|_{t,\omega_{T}}\leqslant\frac{1}{\min_{\omega_{T}}|\phi|}(|v|_{t,\omega_{T}}+C\|w\|_{t-1,\omega_{T}}).$

Hence, by (32),

$\|v-\phi_{h}w_{h}\|_{0,T}\leqslant
Ch^{t}\left(\frac{\|\phi\|_{\infty,T}}{\min_{\omega_{T}}|\phi|}(|v|_{t,\omega_{T}}+\|w\|_{t-1,\omega_{T}})+\|w\|_{0,\omega_{T}}\right).$
(34)

Recall that $\min_{\omega_{T}}|\phi|\geqslant mh$ by Assumption 3. The
distance between any point on $T$ and any point on $\omega_{T}$ is at most
$2h$ so that

$\frac{\|\phi\|_{\infty,T}}{\min_{\omega_{T}}|\phi|}=1+\frac{\max_{T}|\phi|-\min_{\omega_{T}}|\phi|}{\min_{\omega_{T}}|\phi|}\leqslant
1+\frac{2Mh}{mh}\leqslant 1+\frac{2M}{m}$

with $M$ denoting an upper bound on $|\nabla\phi|$. Substituting this into
(34) gives

$\|v-\phi_{h}w_{h}\|_{0,T}\leqslant
Ch^{t}(|v|_{t,\omega_{T}}+\|w\|_{t-1,\omega_{T}}).$ (35)

Summing (31) over all the cells $T\in\mathcal{T}_{h}^{\Gamma,ext}$ and (35)
over all the remaining cells of mesh $\mathcal{T}_{h}$ gives

$\|v-\phi_{h}w_{h}\|_{0,\Omega_{h}}\leqslant
Ch^{t}(|v|_{t,\Omega_{h}}+\|w\|_{t-1,\Omega_{h}})\,.$

This yields (29) with $s=0$ thanks to the estimate
$\|w\|_{t-1,\Omega_{h}}\leqslant\|v\|_{t,\Omega_{h}}$ given by Lemma 5.

Let us now prove (29) for $s=1$. Introduce $v_{h}=I_{h}^{c}v$ so that

$|v-v_{h}|_{1,\Omega_{h}}\leqslant
Ch^{t-1}|v|_{t,\Omega_{h}}\quad\text{and}\quad\|v-v_{h}\|_{0,\Omega_{h}}\leqslant
Ch^{t}|v|_{t,\Omega_{h}}.$

Then, combining the already proven estimate (29) for $s=0$, the inverse
inequality, and the interpolation estimates above, we obtain

$\displaystyle|v-\phi_{h}w_{h}|_{1,\Omega_{h}}$
$\displaystyle\leqslant|v_{h}-\phi_{h}{w_{h}}|_{1,\Omega_{h}}+|v-v_{h}|_{1,\Omega_{h}}$
$\displaystyle\leqslant\frac{C}{h}\|v_{h}-\phi_{h}w_{h}\|_{0,\Omega_{h}}+|v-v_{h}|_{1,\Omega_{h}}$
$\displaystyle\leqslant\frac{C}{h}\|v-\phi_{h}w_{h}\|_{0,\Omega_{h}}+\frac{C}{h}\|v-v_{h}\|_{0,\Omega_{h}}+|v-v_{h}|_{1,\Omega_{h}}$
$\displaystyle\leqslant Ch^{t-1}\|v\|_{t,\Omega_{h}}.$

∎

###  3.3 An adaptation of the Taylor-Hood inf-sup stability to $\phi$-FEM

In this section, we prove some velocity-pressure inf-sup conditions that will
be used to establish the generalized coercivity (inf-sup) for the full
bilinear form in the next section. The proofs are inspired by [17] and start
from an auxiliary inf-sup condition of Lemma 7 with respect to an
$h$-dependent norm for the pressure. The final result in Lemma 8 is weaker
than the usual inf-sup since it does not provide the control of the $L^{2}$
norm of the pressure over the whole domain $\Omega_{h}$. It will be however
sufficient for our purposes since the encumbering term
$-Ch^{2}|p_{h}|^{2}_{1,\Omega_{h}^{\Gamma}}$ will be controlled by the
stabilization present in the scheme, cf. Lemma 2 for this matter.

######  Lemma 7.

There exists an $h$-independent constant $C>0$ such that $\forall
p_{h}\in\mathcal{M}_{h}$, $\exists w_{h}\in\mathcal{V}_{h}$ satisfying

$h^{2}|p_{h}|^{2}_{1,\Omega_{h}}-Ch^{2}|p_{h}|^{2}_{1,\Omega_{h}^{\Gamma}}\leqslant\int_{\Omega_{h}}\nabla
p_{h}\cdot(\phi_{h}w_{h})\quad\text{and}\quad|\phi_{h}w_{h}|_{1,\Omega_{h}}\leqslant
Ch|p_{h}|_{1,\Omega_{h}}.$ (36)

###### Proof.

Let us introduce the submesh
$\mathcal{T}_{h}^{i}=\mathcal{T}_{h}\setminus\mathcal{T}_{h}^{\Gamma}$ so that
$\phi_{h}\leqslant 0$ on $T\in\mathcal{T}_{h}^{i}$. Denote by
$\mathcal{E}_{h}^{i}$ the set of the edges of the mesh $\mathcal{T}_{h}^{i}$
including those shared with $\mathcal{T}_{h}^{\Gamma}$, but excluding those
lying on $\Gamma_{w}$. For any edge $E\in\mathcal{E}_{h}^{i}$, let $t_{E}$ be
the unit tangent vector to $E$ (any of two, but fixed), $x_{E}$ be the
midpoint of $E$, $\omega(E)$ be the set of the mesh cells sharing $E$, and
$\psi_{E}$ be the piecewise quadratic function such that $\psi_{E}(x_{E})=1$
and $\psi_{E}$ vanishes at all the other edge midpoints and at all the nodes
of $\mathcal{T}_{h}$. Moreover, define for all $E\in\mathcal{E}_{h}^{i}$,

$\bar{\phi}_{E}=\left\\{\begin{array}[]{ll}-h,&\text{if }E\text{ belongs to a
cell from }\mathcal{T}_{h}^{\Gamma,ext},\\\
\displaystyle{\frac{1}{|E|}\int_{E}\phi_{h}},&\text{otherwise.}\end{array}\right.$

Take any $p_{h}\in\mathcal{M}_{h}$ and set $w_{h}\in\mathcal{V}_{h}$ as

$w_{h}=\sum_{E\in\mathcal{E}^{i}_{h}}\frac{h^{2}}{\bar{\phi}_{E}}\psi_{E}(t_{E}\cdot\nabla
p_{h})t_{E}.$ (37)

We have indeed $w_{h}\in\mathcal{V}_{h}$, since the pressure tangential
derivative $t_{E}\cdot\nabla p_{h}$ is a continuous piecewise polynomial of
degree
$\leqslant{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}k-2}$
on $\omega(E)$ and $\psi_{E}$ is a continuous piecewise polynomial of degree
$2$, vanishing outside $\omega(E)$. Note also that $w_{h}=0$ on $\Gamma_{w}$
since $\mathcal{E}^{i}_{h}$ does not contain the edges lying on $\Gamma_{w}$.

Take any $E\in\mathcal{E}_{h}^{i}$ and any cell
$T\in\omega(E)\cap\mathcal{T}_{h}^{i}$. We shall see that

$\int_{T}\frac{\phi_{h}}{\bar{\phi}_{E}}\psi_{E}|t_{E}\cdot\nabla
p_{h}|^{2}\geq c\int_{T}|t_{E}\cdot\nabla p_{h}|^{2}\,.$ (38)

Here and elsewhere, the constants $c>0$ depend only on the polynomial degree
$k$, the shape regularity, and the parameters of Assumption 3. To prove (38),
we set $\tilde{\phi}_{h}=\frac{\phi_{h}}{\bar{\phi}_{E}}$ and note that
$\tilde{\phi}_{h}\geqslant 0$ on $T$ since $\phi_{h}\leqslant 0$ on
$T\in\mathcal{T}_{h}^{i}$. To derive further properties of $\tilde{\phi}_{h}$
from Assumption 3, we consider 3 following cases with respect to the placement
of $T$ and $E$ in the mesh (we recall that $T\in\mathcal{T}_{h}^{i}$ in any
case and $E$ is an edge belonging to $T$).

Case 1

$T\in\mathcal{T}_{h}^{\Gamma,ext}$. We have then $\bar{\phi}_{E}=-h$ so that

$\tilde{\phi}_{h}\geqslant 0\text{ and
}h\left|\nabla{\tilde{\phi}_{h}}\right|\geq\frac{m}{2}\text{ on }T\,.$ (39)

Case 2

$T\not\in\mathcal{T}_{h}^{\Gamma,ext}$, $E$ is not shared with any cell from
$\mathcal{T}_{h}^{\Gamma,ext}$. Then $\bar{\phi}_{h}=\phi_{h}(x_{E})$ for some
point $x_{E}\in E$ and we have for any $x\in T$

$\frac{\bar{\phi}_{E}}{\phi_{h}(x)}=\frac{\phi_{h}(x_{E})}{\phi_{h}(x)}=1+\frac{\phi_{h}(x_{E})-\phi_{h}(x)}{\phi_{h}(x)}=1+\frac{\nabla\phi_{h}(c)\cdot(x_{E}-x)}{\phi_{h}(x)}\leq
1+\frac{Mh}{mh}=\frac{m+M}{m}$

where $c\in T$ is a point on the segment connecting $x_{E}$ with $x$. This
implies

$\tilde{\phi}_{h}\geq\frac{m}{m+M}\text{ on }T\,.$ (40)

Case 3

$T\not\in\mathcal{T}_{h}^{\Gamma,ext}$, $E$ is shared with a cell from
$\mathcal{T}_{h}^{\Gamma,ext}$. Then $\bar{\phi}_{E}=-h$ and
$\frac{1}{|E|}\int_{E}\phi_{h}\leq-c_{1}h$ with some $c_{1}>0$ depending only
on the constants in Assumption 3 (since the distance between $E$ and
$\Gamma_{h}$, where $\phi_{h}$ vanishes, is of order $h$ and
$|\nabla\phi_{h}|$ is bounded away from 0 on $\mathcal{T}_{h}^{\Gamma,ext}$).
Combining this with the arguments of the previous case, we arrive at

$\tilde{\phi}_{h}\geq\frac{c_{1}m}{m+M}\text{ on }T\,.$ (41)

Moreover, in all of the 3 cases above,

$|\tilde{\phi}_{h}|\leqslant C\text{ on }T\,.$ (42)

with some $C>0$ depending only on the constants in Assumption 3. In case 1,
this follows from the bound $|\nabla\phi_{h}|\leq M$ and the fact that the
maximal distance between $T$ and $\Gamma_{h}$ is of order $h$. In cases 2 and
3, (42) can be proven in the same way as (40) and (41).

Inequality (38) can be now proven setting

$c=\min_{T,\tilde{\phi}_{h},q_{h}}\frac{\int_{T}{\tilde{\phi}_{h}}\psi_{E}q_{h}^{2}}{\int_{T}q_{h}^{2}}$
(43)

where the minimum is taken over all the simplexes $T$ permitted by the mesh
regularity, all the polynomials
$\tilde{\phi}_{h}:=\frac{\phi_{h}}{\bar{\phi}_{E}}$ of degree $k$ satisfying
(42) and either of (39)–(40)–(41), and all the polynomials
$q_{h}:=t_{E}\cdot\nabla p_{h}\not=0$ of degree $k-2$. By homogeneity and
rescaling, one can safely assume that $h=1$ and $\|q_{h}\|_{0,T}=1$. The sets
of possible $T$, $\tilde{\phi}_{h}$, $q_{h}$ are bounded and closed, so that
the minimum in (43) is indeed attained and $c>0$. Indeed, either of
(39)–(40)–(41) excludes the possibility of $\tilde{\phi}_{h}$ vanishing
everywhere on $T$. This concludes the proof of (38).

Thanks to (38), we have, setting
$\Omega_{h}^{i}=\Omega_{h}\setminus\Omega_{h}^{\Gamma}$ and denoting by
$\mathcal{E}(T)$ the set of edges of a cell $T$ excluding the edges on
$\Gamma_{w}$,

$\int_{\Omega_{h}^{i}}\phi_{h}w_{h}\cdot\nabla
p_{h}=\sum_{E\in\mathcal{E}_{h}^{i}}h^{2}\int_{\omega(E)\cap\Omega_{h}^{i}}\frac{\phi_{h}}{\bar{\phi}_{E}}\psi_{E}|t_{E}\cdot\nabla
p_{h}|^{2}\geq
c\sum_{T\in\mathcal{T}_{h}^{i}}{\sum_{E\in\mathcal{E}(T)}}h^{2}\int_{T}|t_{E}\cdot\nabla
p_{h}|^{2}\,.$

Taking into account Assumption 5, we have by scaling and the equivalence of
norms on all $T\in\mathcal{T}_{h}^{i}$

$\sum_{E\in\mathcal{E}(T)}\int_{T}|t_{E}\cdot\nabla p_{h}|^{2}\geq
c|p_{h}|^{2}_{1,T}\,.$

Hence,

$\int_{\Omega_{h}^{i}}\phi_{h}w_{h}\cdot\nabla p_{h}\geqslant
ch^{2}|p_{h}|^{2}_{1,\Omega^{i}_{h}}.$ (44)

Now, on a mesh cell $T\in\mathcal{T}_{h}^{\Gamma}$ having common edges with
cells of $\mathcal{T}_{h}^{i}$, definition (37) clearly gives

$\left|\int_{T}\phi_{h}w_{h}\cdot\nabla
p_{h}\right|=\left|\sum_{E\in\mathcal{E}(T)\cap\mathcal{E}_{h}^{i}}h\int_{T}\phi_{h}\psi_{E}|t_{E}\cdot\nabla
p_{h}|^{2}\right|\leqslant Ch^{2}|p_{h}|^{2}_{1,T}$

since $\phi_{h}$ is of order $h$ on such cells. Combining this with (44) gives

$\displaystyle ch^{2}|p_{h}|^{2}_{1,\Omega_{h}}$
$\displaystyle\leqslant\int_{\Omega_{h}}\phi_{h}w_{h}\cdot\nabla
p_{h}-\int_{\Omega_{h}^{\Gamma}}\phi_{h}w_{h}\cdot\nabla
p_{h}+ch^{2}|p_{h}|^{2}_{1,\Omega^{\Gamma}_{h}}$
$\displaystyle\leqslant\int_{\Omega_{h}}\phi_{h}w_{h}\cdot\nabla
p_{h}+Ch^{2}|p_{h}|^{2}_{1,\Omega^{\Gamma}_{h}}.$ (45)

We also have

$|\phi_{h}w_{h}|_{1,\Omega_{h}}^{2}=\int_{\Omega_{h}}|(\nabla\phi_{h})w_{h}+\phi_{h}(\nabla
w_{h})|^{2}\\\ \leq
C\sum_{T\in\mathcal{T}_{h}}\sum_{E\in\mathcal{E}(T)\cap\mathcal{E}_{h}^{i}}h^{4}\int_{T}\left(\frac{|\nabla\phi_{h}|^{2}}{\bar{\phi}_{E}^{2}}|\nabla
p_{h}|^{2}+\frac{\phi_{h}^{2}}{\bar{\phi}_{E}^{2}}|\nabla\psi_{E}|^{2}|\nabla
p_{h}|^{2}+\frac{\phi_{h}^{2}}{\bar{\phi}_{E}^{2}}|\nabla^{2}p_{h}|^{2}\right)\\\
\leq C\sum_{T\in\mathcal{T}_{h}}h^{2}\int_{T}|\nabla
p_{h}|^{2}=Ch^{2}|p_{h}|_{1,\Omega_{h}}^{2}.$ (46)

We have used here the finite element inverse estimates on $\phi_{h},\psi_{E}$,
$p_{h}$, and the uniform upper bound (42) on
$|\tilde{\phi}_{h}|=\frac{|\phi_{h}|}{|\bar{\phi}_{E}|}$.

Redefining $w_{h}$ as $w_{h}/c$ with the constant $c$ from (45) shows (36) as
a combination of (45) and (46). ∎

######  Lemma 8.

There exists an $h$-independent constant $C>0$ such that $\forall
p_{h}\in\mathcal{M}_{h}$ $\exists s_{h}\in\mathcal{V}_{h}$

$\|p_{h}\|_{0,\Omega_{h}}^{2}-Ch^{2}|p_{h}|^{2}_{1,\Omega_{h}^{\Gamma}}\leq\int_{\Omega_{h}}\nabla
p_{h}\cdot(\phi_{h}s_{h})\quad\text{ and
}\quad|\phi_{h}s_{h}|_{1,\Omega_{h}}\leqslant C\|p_{h}\|_{0,\Omega_{h}}.$ (47)

###### Proof.

Taking $p_{h}\in\mathcal{M}_{h}$. By continuous velocity-pressure inf-sup
(recall that $\int_{\Omega}p_{h}=0$), there exists $v\in H_{0}^{1}(\Omega)$
s.t.

$\operatorname{div}v=-p_{h}\text{ on }\Omega,\text{ and
}\|v\|_{1,\Omega}\leqslant C\|p_{h}\|_{0,\Omega}\,.$

Let $\tilde{v}\in H^{1}(\Omega_{h})$ be the extension of $v$ by 0 outside
$\Omega$. Lemma 6 with $t=1$ implies $\exists v_{h}\in\mathcal{V}_{h}$ s.t.

$\|\tilde{v}-\phi_{h}v_{h}\|_{0,\Omega_{h}}\leqslant
Ch\|\tilde{v}\|_{1,\Omega_{h}}\leqslant
Ch\|p_{h}\|_{0,\Omega_{{}_{h}}}\quad\text{ and
}\quad|\phi_{h}v_{h}|_{1,\Omega_{h}}\leqslant
C\|\tilde{v}\|_{1,\Omega_{h}}\leqslant C\|p_{h}\|_{0,\Omega_{{}_{h}}}.$

Thus,

$\|p_{h}\|_{0,\Omega}^{2}=-\int_{\Omega}p_{h}\operatorname{div}v=\int_{\Omega}\nabla
p_{h}\cdot\tilde{v}=\int_{\Omega_{h}}\nabla
p_{h}\cdot(\phi_{h}v_{h})+\int_{\Omega_{h}}\nabla
p_{h}\cdot(\tilde{v}-\phi_{h}v_{h})\\\ \leqslant\int_{\Omega_{h}}\nabla
p_{h}\cdot(\phi_{h}v_{h})+Ch|p_{h}|_{1,\Omega_{h}}\|p_{h}\|_{0,\Omega_{{}_{h}}}.$
(48)

We have

$\|p_{h}\|_{0,\Omega_{h}}^{2}\leqslant
C_{1}(\|p_{h}\|_{0,\Omega}^{2}+h^{2}|p_{h}|_{1,\Omega_{h}^{\Gamma}}^{2}).$

This can be proven by an argument similar to that in Lemma 2: one can consider
the maximum of
$\frac{\|p_{h}\|_{0,\Pi_{k}}^{2}}{\|p_{h}\|_{0,T_{k}}^{2}+h^{2}|p_{h}|_{1,\Pi_{k}^{\Gamma}}^{2}}$
taken over all the admissible patches $\Pi_{k}$, as in Assumption 4, and
piecewise polynomials $p_{h}$, observe that this maximum is attained and takes
the value $C_{1}>0$, and sum up over all the patches covering
$\Omega_{h}^{\Gamma}$.
We can thus pass from the norm on $\Omega$ to that on $\Omega_{h}$ in (48):

$\|p_{h}\|_{0,\Omega_{h}}^{2}\leqslant C_{1}\int_{\Omega_{h}}\nabla
p_{h}\cdot(\phi_{h}v_{h})+Ch|p_{h}|_{1,\Omega_{h}}\|p_{h}\|_{0,\Omega_{h}}+C_{1}h^{2}|p_{h}|_{1,\Omega_{h}^{\Gamma}}^{2}.$

Hence, by Young inequality,

$\frac{1}{2}\|p_{h}\|_{0,\Omega_{h}}^{2}\leqslant C_{1}\int_{\Omega_{h}}\nabla
p_{h}\cdot(\phi_{h}v_{h})+C_{2}h^{2}|p_{h}|_{1,\Omega_{h}}^{2}+C_{1}h^{2}|p_{h}|_{1,\Omega_{h}^{\Gamma}}^{2}\\\
\leqslant C_{1}\int_{\Omega_{h}}\nabla
p_{h}\cdot(\phi_{h}v_{h})+C_{2}\int_{\Omega_{h}}\nabla
p_{h}\cdot(\phi_{h}w_{h})+Ch^{2}|p_{h}|_{1,\Omega_{h}^{\Gamma}}^{2}$

with $w_{h}$ given by Lemma 7. Thus,

$\frac{1}{2}\|p_{h}\|_{0,\Omega_{h}}^{2}-Ch^{2}|p_{h}|_{1,\Omega_{h}^{\Gamma}}^{2}\leqslant\int_{\Omega_{h}}\nabla
p_{h}\cdot\phi_{h}(C_{1}v_{h}+C_{2}w_{h})$

and

$|\phi_{h}(C_{1}v_{h}+C_{2}w_{h})|_{1,\Omega_{h}}\leqslant
C\|p_{h}\|_{0,\Omega_{{}_{h}}}$

since $|\phi_{h}w_{h}|_{1,\Omega_{h}}\leqslant
Ch|p_{h}|_{1,\Omega_{h}}\leqslant C\|p_{h}\|_{0,\Omega_{h}}$ by inverse finite
element estimates.

Setting $s_{h}=2(C_{1}v_{h}+C_{2}w_{h})$ proves (47). ∎

###  3.4 The generalized coercivity (the inf-sup condition) for the bilinear
form.

To ease the forthcoming calculations, let us introduce the finite element
space of velocities combining the rigid body motion on the approximate
boundary and the contributions involving the level set:

$\mathcal{V}_{h}^{rbm}=\\{\chi_{h}(V_{h}+\omega_{h}\times
r)+\phi_{h}s_{h}\text{ with
}s_{h}\in\mathcal{V}_{h},V_{h}\in\mathbb{R}^{d},\omega_{h}\in\mathbb{R}^{d^{\prime}}\\}.$
(49)

In the scheme (11), we shall now combine the test functions
$s_{h},V_{h},\omega_{h}$ into $v_{h}\in\mathcal{V}_{h}^{rbm}$ as in the
definition above. Similarly, we shall combine the trial functions
$w_{h},U_{h},\psi_{h}$ into $u_{h}\in\mathcal{V}_{h}^{rbm}$ setting
$u_{h}=\chi_{h}(U_{h}+\psi_{h}\times r)+\phi_{h}w_{h})$. Scheme (11) can be
then rewritten in the compact form: find $u_{h}\in\mathcal{V}_{h}^{rbm}$ and
$p_{h}\in\mathcal{M}_{h}$ such that

$c_{h}(u_{h},p_{h};v_{h},q_{h})=L_{h}(v_{h},q_{h}),\quad\forall
v_{h}\in\mathcal{V}_{h}^{rbm},\ q_{h}\in\mathcal{M}_{h}\,,$ (50)

where the bilinear form $c_{h}$ is given by

$c_{h}(u_{h},p_{h};v_{h},q_{h})=2\int_{\Omega_{h}}D(u_{h}):D(v_{h})-\int_{G_{h}}(2D(u_{h})-p_{h}I)n\cdot\phi_{h}s_{h}\\\
-\int_{\Omega_{h}}q_{h}\operatorname{div}u_{h}-\int_{\Omega_{h}}p_{h}\operatorname{div}v_{h}\\\
+\sigma h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}(-\Delta u_{h}+\nabla
p_{h})\cdot(-\Delta v_{h}-\nabla
q_{h})+\sigma\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}(\operatorname{div}u_{h})(\operatorname{div}v_{h})\\\
+\sigma_{u}h\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}\left[\partial_{n}u_{h}\right]\cdot\left[\partial_{n}v_{h}\right]+\sigma_{u}h^{3}\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}\left[\partial_{n}^{2}u_{h}\right]\cdot\left[\partial_{n}^{2}v_{h}\right]$

and the linear form $L_{h}$ is given by

$L_{h}(v_{h},q_{h})=\int_{\Omega_{h}}\rho_{f}g\cdot\phi_{h}s_{h}+\int_{\mathcal{O}}\rho_{f}g\cdot\chi_{h}(V_{h}+\omega_{h}\times
r)+\left(1-\frac{\rho_{f}}{\rho_{s}}\right)mg\cdot V_{h}\\\ +\sigma
h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}\rho_{f}g\cdot(-\Delta
v_{h}-\nabla q_{h}).$ (51)

In both expressions above, $s_{h}$, $V_{h}$, and $\omega_{h}$ are related to
$v_{h}\in\mathcal{V}_{h}^{rbm}$ as in (49).

######  Lemma 9.

Introduce the norm on $\mathcal{V}_{h}^{rbm}\times\mathcal{M}_{h}$

$\interleave
v_{h},q_{h}\interleave_{h}:=\left(|v_{h}|_{1,\Omega_{h}}^{2}+\|q_{h}\|_{0,\Omega_{h}}^{2}+h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\|-\Delta
v_{h}+\nabla q_{h}\|_{0,T}^{2}+J_{u}(v_{h},v_{h})\right)^{1/2}.$

The following inf-sup condition holds provided $\sigma$ and $\sigma_{u}$ are
sufficiently large:

$\forall(u_{h},p_{h})\in\mathcal{V}_{h}^{rbm}\times\mathcal{M}_{h}\quad\exists(v_{h},q_{h})\in\mathcal{V}_{h}^{rbm}\times\mathcal{M}_{h}$

such that

$\frac{c_{h}(u_{h},p_{h};v_{h},q_{h})}{\interleave
v_{h},q_{h}\interleave_{h}}\geqslant\theta\interleave
u_{h},p_{h}\interleave_{h}$ (52)

with a constant $\theta>0$ depending only on the mesh regularity.

###### Proof.

Let us take $p_{h}\in\mathcal{M}_{h}$ and $u_{h}=\chi_{h}(U_{h}+\psi_{h}\times
r)+\phi_{h}w_{h}\in\mathcal{V}_{h}^{rbm}$ with $w_{h}\in\mathcal{V}_{h}$,
$U_{h}\in\mathbb{R}^{d}$ and $\psi_{h}\in\mathbb{R}^{d^{\prime}}$.

Step 1: controlling the velocity. By choosing $(u_{h},p_{h})$ as the trial
function and $(u_{h},-p_{h})$ as the test function in the bilinear form
$c_{h}$, we obtained:

$c_{h}(u_{h},p_{h};u_{h},-p_{h})=2\int_{\Omega_{h}}|D(u_{h})|^{2}-\int_{G_{h}}(2D(u_{h})-p_{h}I)n\cdot\phi_{h}w_{h}\\\
+\sigma_{u}h\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}\left|\left[\partial_{n}u_{h}\right]\right|^{2}+\sigma_{u}h^{3}\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}\left|\left[\partial_{n}^{2}u_{h}\right]\right|^{2}\\\
+\sigma h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}|-\Delta u_{h}+\nabla
p_{h}|^{2}+\sigma\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}|\operatorname{div}u_{h}|^{2}.$
(53)

Let $\tilde{B}_{h}$ be the strip between $\Gamma_{h}=\\{\phi_{h}=0\\}$ and
$G_{h}$, i.e. $\tilde{B}_{h}=\\{\phi_{h}>0\\}\cap\Omega_{h}$. Since
$\phi_{h}w_{h}=0$ on $\Gamma_{h}$,

$\int_{G_{h}}D(u_{h})n\cdot\phi_{h}w_{h}=\int_{\partial\tilde{B}_{h}}D(u_{h})n\cdot\phi_{h}w_{h}\\\
=\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{\partial(\tilde{B}_{h}\cap
T)}D(u_{h})n_{T}\cdot\phi_{h}w_{h}-\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\sum_{E\in\mathcal{F}_{h}^{cut}(T)}\int_{\tilde{B}_{h}\cap
E}D(u_{h})n_{T}\cdot\phi_{h}w_{h},$

where $\mathcal{T}_{h}^{\Gamma}$ is defined in (10),
$\mathcal{F}_{h}^{cut}(T)$ regroups the facets of a mesh element $T$ cut by
$\Gamma_{h}$, and $n_{T}$ is the unit normal pointing outside of $T$ on the
boundary of a mesh cell $T$. Applying the divergence theorem to the integrals
on ${\partial(\tilde{B}_{h}\cap T)}$ and regrouping the integrals on the
facets gives

$\int_{G_{h}}D(u_{h})n\cdot\phi_{h}w_{h}=\int_{\tilde{B}_{h}}D(u_{h}):D(\phi_{h}w_{h})+\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{\tilde{B}_{h}\cap
T}\operatorname{div}D(u_{h})\cdot\phi_{h}w_{h}\\\
-\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E\cap\tilde{B}_{h}}[D(u_{h})n]\cdot\phi_{h}w_{h}\\\
=\int_{\tilde{B}_{h}}|D(u_{h})|^{2}+\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{\tilde{B}_{h}\cap
T}\operatorname{div}D(u_{h})\cdot\phi_{h}w_{h}-\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E\cap\tilde{B}_{h}}[D(u_{h})n]\cdot\phi_{h}w_{h},$

since $u_{h}-\phi_{h}w_{h}$ is the velocity of a rigid motion on
$\tilde{B}_{h}$.

Similarly (and simpler)

$\int_{G_{h}}p_{h}n\cdot\phi_{h}w_{h}=\int_{\partial\tilde{B}_{h}}p_{h}n\cdot\phi_{h}w_{h}=\int_{\tilde{B}_{h}}p_{h}\operatorname{div}u_{h}+\int_{\tilde{B}_{h}}\nabla
p_{h}\cdot\phi_{h}w_{h}.$

Substituting this into (53) and rewriting $2\operatorname{div}D(u_{h})-\nabla
p_{h}=\Delta u_{h}-\nabla p_{h}+\nabla\operatorname{div}u_{h}$ on the cells
$T\in\mathcal{T}_{h}^{\Gamma}$ yields

$c_{h}(u_{h},p_{h};u_{h},-p_{h})=2\int_{\Omega_{h}}|D(u_{h})|^{2}-2\int_{\tilde{B}_{h}}|D(u_{h})|^{2}-\underbrace{\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{\tilde{B}_{h}\cap
T}(\Delta u_{h}-\nabla p_{h})\cdot\phi_{h}w_{h}}_{\text{Young with
$\varepsilon_{1}$}}\\\
+\underbrace{2\sum_{F\in\mathcal{F}_{h}^{\Gamma}}\int_{F\cap\tilde{B}_{h}}[D(u_{h})n]\cdot\phi_{h}w_{h}}_{\text{Young
with
$\varepsilon_{2}$}}+\underbrace{\int_{\tilde{B}_{h}}p_{h}\operatorname{div}u_{h}}_{\text{Young
with
$\varepsilon_{3}$}}-\underbrace{\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{\tilde{B}_{h}\cap
T}(\nabla\operatorname{div}u_{h})\cdot\phi_{h}w_{h}}_{\text{Young with
$\varepsilon_{4}$}}\\\
+\sigma_{u}h\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}|\left[\partial_{n}u_{h}\right]|^{2}+\sigma_{u}h^{3}\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}\left|\left[\partial_{n}^{2}u_{h}\right]\right|^{2}\\\
+\sigma h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}|\Delta u_{h}-\nabla
p_{h}|^{2}+\sigma\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}|\operatorname{div}u_{h}|^{2}.$

Several terms above are marked with “Young with $\varepsilon_{i}$” meaning
that we are going to apply the Young inequality with some weights
$\varepsilon_{1},\ldots,\varepsilon_{4}>0$ (multiplied by the appropriate
powers of $h$) to these terms. We recall that Lemma 4 implies

$\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\frac{1}{h^{2}}\|\phi_{h}w_{h}\|_{\tilde{B}_{h}\cap
T}^{2}\leq
C\|D(u_{h})\|_{0,\Omega_{h}}^{2}\,,\quad\sum_{F\in\mathcal{F}_{h}^{\Gamma}}\frac{1}{h}\|\phi_{h}w_{h}\|_{F\cap\tilde{B}_{h}}^{2}\leq
C\|D(u_{h})\|_{0,\Omega_{h}}^{2}\,,$

which allows us to absorb the norms of $\phi_{h}w_{h}$ into the first term
with $\|D(u_{h})\|_{0,\Omega_{h}}$. We also use the inverse inequality
$h\|\nabla\operatorname{div}u_{h}\|_{0,\Omega_{h}^{\Gamma}}\leq
C\|\operatorname{div}u_{h}\|_{0,\Omega_{h}^{\Gamma}}$. This yields

$c_{h}(u_{h},p_{h};u_{h},-p_{h})\geqslant\left(2-C{\frac{\varepsilon_{1}+\varepsilon_{2}+\varepsilon_{4}}{2}}\right)\|D(u_{h})\|_{0,\Omega_{h}}^{2}-2\|D(u_{h})\|_{0,\Omega^{\Gamma}_{h}}^{2}-\frac{\varepsilon_{3}}{2}\|p_{h}\|_{0,\Omega_{h}^{\Gamma}}\\\
+h^{2}\left(\sigma-\frac{1}{2\varepsilon_{1}}\right)\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\|\Delta
u_{h}-\nabla
p_{h}\|_{0,T}^{2}+\left(\sigma-\frac{1}{2\varepsilon_{3}}-\frac{C}{\varepsilon_{4}}\right)\|\operatorname{div}u_{h}\|_{0,\Omega_{h}^{\Gamma}}^{2}\\\
+h\left(\sigma_{u}-\frac{1}{2\varepsilon_{2}}\right)\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\|\left[\partial_{n}u_{h}\right]\|_{0,E}^{2}+\sigma_{u}h^{3}\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\left\|\left[\partial_{n}^{2}u_{h}\right]\right\|_{0,E}^{2}$

Thanks to Lemma 2, this can be further bounded as

$c_{h}(u_{h},p_{h};u_{h},-p_{h})\geqslant\left(2(1-\alpha)-C{\frac{\varepsilon_{1}+\varepsilon_{2}+\varepsilon_{4}}{2}}\right)\|D(u_{h})\|_{0,\Omega_{h}}^{2}+2(1-\alpha)h^{2}|p_{h}|_{1,\Omega_{h}^{\Gamma}}\\\
-\frac{\varepsilon_{3}}{2}\|p_{h}\|_{0,\Omega_{h}^{\Gamma}}+h^{2}\left(\sigma-\frac{1}{2\varepsilon_{1}}-2\beta\right)\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\|\Delta
u_{h}-\nabla p_{h}\|_{0,T}^{2}\\\
+\left(\sigma-\frac{1}{2\varepsilon_{3}}-\frac{C}{\varepsilon_{4}}-2\beta\right)\|\operatorname{div}u_{h}\|_{0,\Omega_{h}^{\Gamma}}^{2}+h\left(\sigma_{u}-\frac{1}{2\varepsilon_{2}}-2\beta\right)\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\|\left[\partial_{n}u_{h}\right]\|_{0,E}^{2}\\\
+h^{3}(\sigma_{u}-2\beta)\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\left\|\left[\partial_{n}^{2}u_{h}\right]\right\|_{0,E}^{2}$
(54)

with some $\beta>0$ and $\alpha\in(0,1)$.

Step 2: controlling the pressure. Let now $s_{h}\in\mathcal{V}_{h}$ be the
function given by Lemma 8 and set $v_{h}^{p}=\phi_{h}s_{h}$. Noting that

$-\int_{\Omega_{h}}p_{h}\operatorname{div}v_{h}^{p}+\int_{G_{h}}p_{h}n\cdot
v_{h}^{p}=\int_{\Omega_{h}}\nabla p_{h}\cdot
v_{h}^{p}\geqslant\|p_{h}\|_{0,\Omega_{h}}^{2}-Ch^{2}|p_{h}|^{2}_{1,\Omega_{h}^{\Gamma}}$

we get

$c_{h}(u_{h},p_{h};v_{h}^{p},0)\geqslant\|p_{h}\|_{0,\Omega_{h}}^{2}-Ch^{2}|p_{h}|^{2}_{1,\Omega_{h}^{\Gamma}}+2\int_{\Omega_{h}}D(u_{h}):D(v_{h}^{p})-2\int_{G_{h}}D(u_{h})n\cdot
v_{h}^{p}\\\ +\sigma h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}(\Delta
u_{h}-\nabla p_{h})\cdot\Delta
v_{h}^{p}+\sigma\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}(\operatorname{div}u_{h})(\operatorname{div}v_{h}^{p})\\\
+\sigma_{u}h\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}\left[\partial_{n}u_{h}\right]\cdot\left[\partial_{n}v_{h}^{p}\right]+\sigma_{u}h^{3}\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}\left[\partial_{n}^{2}u_{h}\right]\cdot\left[\partial_{n}^{2}v_{h}^{p}\right].$

Recalling that $|v_{h}^{p}|_{1,\Omega_{h}}\leqslant
C\|p_{h}\|_{0,\Omega_{h}}$, remarking that
$\|v_{h}^{p}\|_{0,G_{h}}\leq\frac{C}{\sqrt{h}}|v_{h}^{p}|_{1,\Omega_{h}^{\Gamma}}$,
and applying Young and inverse inequalities allows us to conclude

$c_{h}(u_{h},p_{h};v_{h}^{p},0)\geqslant\frac{1}{2}\|p_{h}\|_{0,\Omega_{h}}^{2}-Ch^{2}|p_{h}|^{2}_{1,\Omega_{h}^{\Gamma}}-C\|D(u_{h})\|_{0,\Omega_{h}}^{2}\\\
-C\left(\sigma^{2}h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\|-\Delta
u_{h}+\nabla
p_{h}\|_{0,T}^{2}+\sigma^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\|\operatorname{div}u_{h}\|_{0,T}^{2}\right.\\\
\left.\sigma_{u}^{2}h\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\left\|\left[\partial_{n}u_{h}\right]\right\|_{0,E}^{2}+\sigma_{u}^{2}h^{3}\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\left\|\left[\partial_{n}^{2}u_{h}\right]\right\|_{0,E}^{2}\right).$
(55)

Step 3: combining the estimates. Multiply (55) by $\lambda>0$ and add it to
(54). This gives

$c_{h}(u_{h},p_{h};u_{h}+\lambda
v_{h}^{p},-p_{h})\geqslant\left(2(1-\alpha)-C\frac{\varepsilon_{1}+\varepsilon_{2}+\varepsilon_{4}}{2}-C\lambda\right)\|D(u_{h})\|_{0,\Omega_{h}}^{2}\\\
+\frac{\lambda-\varepsilon_{3}}{2}\|p_{h}\|_{0,\Omega_{h}}^{2}+(2(1-\alpha)-C\lambda)h^{2}|p_{h}|^{2}_{1,\Omega_{h}^{\Gamma}}\\\
+h^{2}\left(\sigma-\frac{1}{2\varepsilon_{1}}-2\beta-C\lambda\sigma^{2}\right)\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\|\Delta
u_{h}-\nabla p_{h}\|_{0,T}^{2}\\\
+\left(\sigma-\frac{1}{2\varepsilon_{3}}-\frac{C}{\varepsilon_{4}}-2\beta-C\lambda\sigma^{2}\right)\|\operatorname{div}u_{h}\|_{0,\Omega_{h}^{\Gamma}}^{2}\\\
+h\left(\sigma_{u}-\frac{1}{2\varepsilon_{2}}-2\beta-C\lambda\sigma_{u}^{2}\right)\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\|\left[\partial_{n}u_{h}\right]\|_{0,E}^{2}+(\sigma_{u}-2\beta-C\lambda\sigma_{u}^{2})h^{3}\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\left\|\left[\partial_{n}^{2}u_{h}\right]\right\|_{0,E}^{2}.$

Taking $\lambda,\varepsilon_{1},\varepsilon_{2},\varepsilon_{4}$ small enough,
$\varepsilon_{3}=\lambda/2$ and $\sigma,\sigma_{u}$ big enough, and recalling
Korn inequality (22), this amounts to

$c_{h}(u_{h},p_{h};w_{h}+\lambda v_{h}^{p},-p_{h})\geqslant c\interleave
u_{h},p_{h}\interleave_{h}^{2}$

with some $c>0$. We also have easily

$\interleave w_{h}+\lambda v_{h}^{p},-p_{h}\interleave_{h}\leqslant
C\interleave w_{h},p_{h}\interleave_{h},$

hence the inf-sup estimate (52) with $v_{h}=u_{h}+\lambda v_{h}^{p}$,
$q_{h}=-p_{h}$. ∎

###  3.5 A priori error estimates.

In this section, we will prove Theorem 1 following the argumentation of [12],
which is ameliorated since we require only the optimal regularity
$H^{k+1}(\Omega)^{d}\times H^{k}(\Omega)$ for the velocity-pressure pair
$(u,p)$ given by (1a)-(1g).

Proof of the $H^{1}$ a priori error estimate (12): Let $(u,p)\in
H^{k+1}(\Omega)^{d}\times H^{k}(\Omega)$ with $u=U+\psi\times r$ on $\Gamma$
be the solution to the continuous problem (1a)-(1g) and
$(u_{h},p_{h})\in\mathcal{V}_{h}^{rbm}\times\mathcal{M}_{h}$ with
$u_{h}=\phi_{h}w_{h}+\chi_{h}(U_{h}+\psi_{h}\times r)$ be the solution to the
discrete problem (11). Choose sufficiently smooth extension $\tilde{u}$ and
$\tilde{p}$ of $u$ and $p$ on $\Omega_{h}$ such that $\tilde{u}=u$,
$\tilde{p}=p$ on $\Omega$, and

$\|\tilde{u}\|_{k+1,\Omega_{h}}\leqslant
C\|u\|_{k+1,\Omega},\quad\|\tilde{p}\|_{k,\Omega_{h}}\leqslant
C\|p\|_{k,\Omega}.$

Applying Lemma 6 to $\tilde{u}-\chi(U+\psi\times r)$, which vanishes on
$\Gamma$ and on $\Gamma_{w}$, we see that there exists
$\tilde{w}_{h}\in\mathcal{V}_{h}$ such that

$\|\tilde{u}-\chi(U+\psi\times
r)-\phi_{h}\tilde{w}_{h}\|_{1,\Omega_{h}}\leqslant
Ch^{k}\|\tilde{u}-\chi(U+\psi\times r)\|_{k+1,\Omega_{h}}.$

This allows us to introduce
$\tilde{u}_{h}=\phi_{h}\tilde{w}_{h}+\chi_{h}(U+\psi\times
r)\in\mathcal{V}_{h}^{rbm}$ satisfying

$\displaystyle\|\tilde{u}-\tilde{u}_{h}\|_{1,\Omega_{h}}$
$\displaystyle\leqslant$ $\displaystyle\|\tilde{u}-\chi(U+\psi\times
r)-\phi_{h}\tilde{w}_{h}\|_{1,\Omega_{h}}+\|(\chi-\chi_{h})(U+\psi\times
r)\|_{1,\Omega_{h}}$ (56) $\displaystyle\leqslant$ $\displaystyle
Ch^{k}(\|\tilde{u}\|_{k+1,\Omega_{h}}+\|\chi\|_{k+1,\Omega_{h}}(|U|+|\psi|))$
$\displaystyle\leqslant$ $\displaystyle
Ch^{k}\|\tilde{u}\|_{k+1,\Omega_{h}}\leqslant Ch^{k}\|u\|_{k+1,\Omega},$

thanks to the standard interpolation of $\chi\in H^{k+1}(\Omega_{h})$ and to
the bounds $|U|,|\psi|\leq C\|u\|_{1,\Omega}$ valid by the trace inequality
(recall that $u=U+\psi\times r$ on $\Gamma$).

Similarly, $\|\tilde{u}-\tilde{u}_{h}\|_{0,\Omega_{h}}\leqslant
Ch^{k+1}\|u\|_{k+1,\Omega}$. We define moreover
$\tilde{p}_{h}\in\mathcal{M}_{h}$ by the standard FE nodal interpolation
$\tilde{p}_{h}=I_{h}\tilde{p}$ such that

$\|\tilde{p}-\tilde{p}_{h}\|_{0,\Omega_{h}}\leqslant
Ch^{k}\|p\|_{k,\Omega}\,.$ (57)

Thanks to Lemma 9,
$\exists(v_{h},q_{h})\in\mathcal{V}_{h}^{rbm}\times\mathcal{M}_{h}$ such that

$\interleave\tilde{u}_{h}-u_{h},\tilde{p}_{h}-p_{h}\interleave_{h}\leq
C\dfrac{c_{h}(\tilde{u}_{h}-u_{h},\tilde{p}_{h}-p_{h};v_{h},q_{h})}{\interleave
v_{h},q_{h}\interleave_{h}}.$ (58)

We should now substitute $(\tilde{u},\tilde{p})$ into the form $c_{h.}$ To
this end, we introduce the fictitious right-hand sides $\tilde{F}$ and
$\tilde{Q}$ on $\Omega_{h}$ so that

$-2\operatorname{div}D(\tilde{u})+\nabla\tilde{p}=\tilde{F}\text{ and
}\operatorname{div}\tilde{u}=\tilde{Q}\text{ in }\Omega_{h}.$

We observe then, taking any
$v_{h}=\phi_{h}s_{h}+\chi_{h}(V_{h}+\omega_{h}\times
r)\in\mathcal{V}_{h}^{rbm}$, $q_{h}\in\mathcal{M}_{h}$,

$2\int_{\Omega_{h}}D(\tilde{u}):D(v_{h})-\int_{\Omega_{h}}\tilde{p}\operatorname{div}v_{h}-\int_{G_{h}}(2D(\tilde{u})-\tilde{p}I)n\cdot
v_{h}=\int_{\Omega_{h}}\tilde{F}\cdot v_{h}$

and, recalling $B_{h}=\Omega_{h}\setminus\Omega$,

$\int_{G_{h}}(2D(\tilde{u})-\tilde{p}I)n\cdot(V_{h}+\omega_{h}\times
r)=\int_{\Gamma}(2D(\tilde{u})-pI)n\cdot(V_{h}+\omega_{h}\times r)\\\
+\int_{B_{h}}\text{div}(2D(\tilde{u})-\tilde{p}I)n\cdot(V_{h}+\omega_{h}\times
r)=mg\cdot V_{h}-\int_{B_{h}}\tilde{F}\cdot(V_{h}+\omega_{h}\times r).$

Hence,

$c_{h}(\tilde{u},\tilde{p};v_{h},q_{h})=\int_{\Omega_{h}}\tilde{F}\cdot
v_{h}-\int_{B_{h}}\tilde{F}\cdot(V_{h}+\omega_{h}\times r)+mg\cdot
V_{h}-\int_{\Omega_{h}}q_{h}\operatorname{div}\tilde{u}\\\ +\sigma
h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}\tilde{F}\cdot(-\Delta
v_{h}-\nabla
q_{h})+\sigma\int_{\Omega_{h}^{\Gamma}}(\operatorname{div}\tilde{u})(\operatorname{div}v_{h}).$

Also note that the RHS (51) of the scheme (50) can be rewritten as

$L_{h}(v_{h},q_{h})=\int_{\Omega_{h}}\rho_{f}g\cdot
v_{h}-\int_{{B_{h}}}\rho_{f}g\cdot(V_{h}+\omega_{h}\times r)+mg\cdot V_{h}\\\
+\sigma h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}\rho_{f}g\cdot(-\Delta
v_{h}-\nabla q_{h})$

This allows us to establish the following Galerkin orthogonality relation,
valid for all $v_{h}\in\mathcal{V}_{h}^{rbm}$, $q_{h}\in\mathcal{M}_{h}$,

$c_{h}(\tilde{u}-u_{h},\tilde{p}-p_{h};v_{h},q_{h})=R_{h}(v_{h},q_{h}),$ (59)

where

$R_{h}(v_{h},q_{h})=\int_{B_{h}}(\tilde{F}-\rho_{f}g)\cdot\phi_{h}s_{h}-\int_{B_{h}}q_{h}\operatorname{div}\tilde{u}\\\
+\sigma
h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}(\tilde{F}-\rho_{f}g)\cdot(-\Delta
v_{h}-\nabla
q_{h})+\sigma\int_{\Omega_{h}^{\Gamma}}(\operatorname{div}\tilde{u})(\operatorname{div}v_{h}).$

The integrals of $\tilde{F}-\rho_{f}g$ and $\operatorname{div}\tilde{u}$ on
$\Omega_{h}$ have been rewritten as integrals on $B_{h}$ since both
$\tilde{F}-\rho_{f}g$ and $\operatorname{div}\tilde{u}$ vanish on $\Omega$.

Combination of (58) and (59) entails

$\interleave\tilde{u}_{h}-u_{h},\tilde{p}_{h}-p_{h}\interleave_{h}\leq
C\dfrac{c_{h}(\tilde{u}_{h}-\tilde{u},\tilde{p}_{h}-\tilde{p};v_{h},q_{h})+R_{h}(v_{h},q_{h})}{\interleave
v_{h},q_{h}\interleave_{h}}.$

We can now use interpolation inequalities as in [12, Section 3.4]. In
particular, the term with $c_{h}$ in the nominator of the fraction above is
bounded by $Ch^{k}(\|u\|_{k+1,\Omega}+\|p\|_{k,\Omega}){\interleave
v_{h},q_{h}\interleave_{h}}$ thanks to (56)–(57) and to the estimates of Lemma
4. To bound $R_{h}(v_{h},q_{h})$ we recall that $\tilde{F}-\rho_{f}g$ and
$\operatorname{div}\tilde{u}$ vanish on $\Omega$. Thus, thanks to [12, Lemma
3.6]

$\|\tilde{F}-\rho_{f}g\|_{0,\Omega_{h}^{\Gamma}}\leq
Ch^{k-1}\|\tilde{F}-\rho_{f}g\|_{k-1,\Omega_{h}^{\Gamma}}\leq
Ch^{k-1}(\|u\|_{k+1,\Omega}+\|p\|_{k,\Omega})$

and

$\|\operatorname{div}\tilde{u}\|_{0,\Omega_{h}^{\Gamma}}\leq
Ch^{k}\|\operatorname{div}\tilde{u}\|_{k,\Omega_{h}^{\Gamma}}\leq
Ch^{k}\|u\|_{k+1,\Omega}.$ (60)

This, combined with the estimates of Lemma 4, in particular
$\|\phi_{h}s_{h}\|_{0,\Omega_{h}^{\Gamma}}\leq
Ch|v_{h}|_{1,\Omega_{h}^{\Gamma}}$, leads to $|R_{h}(v_{h},q_{h})|\leq
Ch^{k}(\|u\|_{k+1,\Omega}+\|p\|_{k,\Omega}){\interleave
v_{h},q_{h}\interleave_{h}}$ and

$\interleave\tilde{u}_{h}-u_{h},\tilde{p}_{h}-p_{h}\interleave_{h}\leq
Ch^{k}(\|u\|_{k+1,\Omega}+\|p\|_{k,\Omega})\,.$ (61)

Recalling again the interpolation estimates (56)–(57), we obtain the error
estimates for $u$ and $p$, announced by (12).

Proof of the a priori error estimate (13) on the velocity of the solid: We
have by the construction of the interpolant
$\tilde{u}_{h}=\phi_{h}\tilde{w}_{h}+\chi_{h}(U+\psi\times r)$ and thanks to
(26)

$|U-U_{h}|+|\psi-\psi_{h}|\leqslant
C\|\phi_{h}(\tilde{w}_{h}-w_{h})+\chi_{h}(U-U_{h}+(\psi-\psi_{h})\times
r)\|_{1,\Omega_{h}}=C\|\tilde{u}_{h}-u_{h}\|_{1,\Omega_{h}}$

which proves (13) thanks to (61).

Proof of the $L^{2}$ a priori error estimate (14): Let $(v,q,V,\omega)\in
H^{2}(\Omega)^{d}\times
H^{1}(\Omega)\times\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}$ the solution
to

$\left\\{\begin{array}[]{ll}-2\operatorname{div}D(v)+\nabla q=u-u_{h},&\text{
in }\Omega,\\\ \operatorname{div}v=0,&\text{ in }\Omega,\\\ v=V+\omega\times
r,&\text{ on }\Gamma,\\\ v=0,&\text{ on }\Gamma_{w},\\\\[1.0pt]
\int_{\Gamma}(2D(v)-qI)n=0,&\\\\[3.0pt] \int_{\Gamma}(2D(v)-qI)n\times
r=0,&\\\\[3.0pt] \int_{\Omega}q=0.&\end{array}\right.$

An integration by parts gives

$\displaystyle\|u-u_{h}\|_{0,\Omega}^{2}$
$\displaystyle=\int_{\Omega}(u-u_{h})(-2\operatorname{div}D(v)+\nabla q)$ (62)
$\displaystyle=2\int_{\Omega}D(u-u_{h}):D(v)-\int_{\Omega}q\operatorname{div}(u-u_{h})-\int_{\Omega}(p-p_{h})\operatorname{div}v.$
$\displaystyle\qquad+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{\Gamma}((\phi_{h}-\phi)w_{h}+(\chi_{h}-\chi)(U_{h}+\phi_{h}\times
r))\cdot(2D(v)-qI)n}$

Note that the boundary term $\int_{\Gamma}u\cdot(2D(v)-qI)n$ vanishes since
$u$ is a rigid body motion on $\Gamma$. For the same reason,
$\int_{\Gamma}(\phi w_{h}+\chi(U_{h}+\phi_{h}\times r)\cdot(2D(v)-qI)n$
vanishes.

Let $(\tilde{v},\tilde{q})\in H^{2}(\Omega_{h})^{d}\times H^{1}(\Omega_{h})$
coincide with $(v,q)$ on $\Omega$. They can be constructed by a bounded
extension operator in $H^{2}\times H^{1}$ so that

$\|\tilde{v}\|_{2,\Omega_{h}}+\|\tilde{q}\|_{1,\Omega_{h}}\leqslant
C(\|v\|_{2,\Omega}+\|q\|_{1,\Omega})\leqslant C\|u-u_{h}\|_{0,\Omega}.$ (63)

We now further rewrite (62) using Galerkin orthogonality (59) with the test
functions $v_{h}=\phi_{h}s_{h}+\chi_{h}(V_{h}+\omega_{h}\times
r)\in\mathcal{V}_{h}^{rbm}$ and $q_{h}\in\mathcal{M}_{h}$ and recalling
$B_{h}=\Omega_{h}\backslash\Omega$,

$\|u-u_{h}\|_{0,\Omega}^{2}\\\
=\underbrace{\int_{\Omega_{h}}2D(\tilde{u}-u_{h}):D(\tilde{v}-v_{h})-\int_{\Omega_{h}}(\tilde{q}-q_{h})\operatorname{div}(\tilde{u}-u_{h})-\int_{\Omega_{h}}(\tilde{p}-p_{h})\operatorname{div}(\tilde{v}-v_{h})}\limits_{I}\\\
-\underbrace{\left(\int_{B_{h}}2D(\tilde{u}-u_{h}):D(\tilde{v})-\int_{B_{h}}\tilde{q}\operatorname{div}(\tilde{u}-u_{h})-\int_{B_{h}}(\tilde{p}-p_{h})\operatorname{div}\tilde{v}\right)}\limits_{II}\\\
+\underbrace{\int_{G_{h}}(2D(\tilde{u}-u_{h})-(\tilde{p}-p_{h})I)n\cdot\phi_{h}s_{h}}\limits_{III}$
(64)
$\underbrace{-\sigma_{u}h\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}\left[\frac{\partial}{\partial
n}(\tilde{u}-u_{h})\right]\cdot\left[\frac{\partial v_{h}}{\partial
n}\right]-\sigma_{u}h^{3}\sum_{E\in\mathcal{F}_{h}^{\Gamma}}\int_{E}\left[\frac{\partial^{2}}{\partial
n^{2}}(\tilde{u}-u_{h})\right]\cdot\left[\frac{\partial^{2}v_{h}}{\partial
n^{2}}\right]}\limits_{IV}\\\ -\underbrace{\sigma
h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}-\Delta(\tilde{u}-u_{h})+\nabla(\tilde{p}-p_{h}))\cdot(-\Delta
v_{h}-\nabla q_{h})}\limits_{V}\\\
-\underbrace{\sigma\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}\operatorname{div}(\tilde{u}-u_{h})\operatorname{div}v_{h}}\limits_{VI}+\underbrace{R_{h}(v_{h},q_{h})}\limits_{VII}\\\
+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\underbrace{\int_{\Gamma}((\phi_{h}-\phi)w_{h}+(\chi_{h}-\chi)(U_{h}+\phi_{h}\times
r))\cdot(2D(v)-qI)n}\limits_{VIII}}$

We now take $V_{h}=V,\omega_{h}=\omega$ and set $s_{h}\in\mathcal{V}_{h}$ so
that $\phi_{h}s_{h}$ is an optimal interpolant of $v-\chi_{h}(V+\omega\times
r)$, as guaranteed by Lemma 6. We also set $q_{h}=\tilde{I}_{h}\tilde{q}$
using an appropriate Clément interpolation $\tilde{I}_{h}$. We can now
estimate all the terms of (64) using the already proven estimate (61) and the
interpolation estimates for $\tilde{v}-v_{h}$ and $\tilde{p}-p_{h}$ . This
gives

$\|u-u_{h}\|_{0,\Omega}^{2}\leq
Ch^{k+1/2}(\|u\|_{k+1,\Omega}+\|p\|_{k,\Omega})(\|\tilde{v}\|_{2,\Omega_{h}}+\|\tilde{q}\|_{1,\Omega_{h}}).$
(65)

In particular, term $I$ is completely standard and gives in fact a
contribution of the optimal order $h^{k+1}$. Rather than go to the details of
the tedious calculations leading to the bounds of the remaining terms, we
prefer here to refer to the similar arguments used in [12] to estimate the
terms in eq. (3.24). Indeed, the terms $II-III$ in (64) can be treated as the
terms $II-III$ in eq. (3.24) of [12]. Terms $IV-V$ in (64) can be treated as
term $IV$ in eq. (3.24) of [12]. Terms $VI$ in (64) is also similar to the
latter (note, in particular,
$\|\operatorname{div}\tilde{v}\|_{0,\Omega_{h}^{\Gamma}}\leq
Ch\|v\|_{2,\Omega}$ similarly to (60)). Finally, term $VII$ in (64) can be
treated as term $V$ in eq. (3.24) of [12]. As in [12], all these terms result
in the sub-optimal estimate of order $O(h^{k+1/2})$. The origin of this sub-
optimality lies in the lack of adjoint consistency in formulation (50): the
adjoint discrete problem cannot be interpreted as a consistent discretization
of a meaningful continuous problem.

The only term in (64), which does not have a direct analogue in [12], is term
$VIII$. To bound it, we apply Cauchy-Schwarz inequality together with the
interpolation estimates on $\phi-\phi_{h}$ and $\chi-\chi_{h}$, recalling the
hypotheses $\phi\in C^{k+1}(\Omega_{h})$ and $\chi\in H^{k+1}(\Omega_{h})$:

$|VIII|\leqslant
Ch^{k+1}(\|w_{h}\|_{0,\Gamma}+|U_{h}|+|\phi_{h}|)(\|D(v)\|_{0,\Gamma}+\|q\|_{0,\Gamma})\,.$

Then, to bound $w_{h}$ in $L^{2}(\Gamma)$, we start by a trace inverse
inequality
$\|w_{h}\|_{0,\Gamma}\leqslant\frac{C}{\sqrt{h}}\|w_{h}\|_{0,\Omega_{h}}$ and
apply Hardy inequality of Lemma 5 to $w_{h}=\frac{\phi w_{h}}{\phi}$:

$\|w_{h}\|_{0,\Omega_{h}}\leqslant C\|\phi w_{h}\|_{1,\Omega_{h}}\leqslant
C(\|(\phi-\phi_{h})w_{h}\|_{1,\Omega_{h}}+\|\phi_{h}w_{h}\|_{1,\Omega_{h}})\,.$

Noting that, by interpolation and inverse inequalities,

$\|(\phi-\phi_{h})w_{h}\|_{1,\Omega_{h}}\leqslant
C(\|\phi-\phi_{h}\|_{L^{\infty}(\Omega_{h})}|w_{h}|_{1,\Omega_{h}}+\|\nabla\phi-\nabla\phi_{h}\|_{L^{\infty}(\Omega_{h})}\|w_{h}\|_{0,\Omega_{h}})\leqslant
Ch^{k}\|w_{h}\|_{0,\Omega_{h}}$

we conclude

$(1-Ch^{k})\|w_{h}\|_{0,\Omega_{h}}\leqslant
C\|\phi_{h}w_{h}\|_{1,\Omega_{h}}$

Hence, for $h$ small enough,

$|VIII|\leqslant
Ch^{k+1}\left(\frac{1}{\sqrt{h}}\|\phi_{h}w_{h}\|_{1,\Omega_{h}}+|U_{h}|+|\phi_{h}|\right)(\|D(v)\|_{0,\Gamma}+\|q\|_{0,\Gamma})\,.$

Recalling that $u_{h}=\phi_{h}w_{h}+\chi_{h}(U_{h}+\phi_{h}\times r)$, we
conclude by Lemma 4 that $\|\phi_{h}w_{h}\|_{1,\Omega_{h}}$, $|U_{h}|$,
$|\phi_{h}|$ can be all bounded by $\|u_{h}\|_{1,\Omega_{h}}$. Applying the
trace inequalities to $v$ and $q$, we arrive at

$|VIII|\leqslant
Ch^{k+1/2}\|u_{h}\|_{1,\Omega_{h}}(\|v\|_{2,\Omega}+\|q\|_{1,\Omega})\,.$

Since we know that $\|u_{h}\|_{1,\Omega_{h}}$ is bounded by the norms of $u$
and $p$ thanks to the already proven error estimates for the velocity in
$H^{1}$ norm, we conclude that term $VIII$ contributes to (65) in the same
manner as all the other terms.

Combining (65) with (63) proves (14).

##  4 Numerical tests

In this section, we present numerical results, first in the particular case of
a fixed particle, i.e. for the Stokes equations aone in a fixed domain (cf.
Appendix A and the $\phi$-FEM scheme (72)), and second in the case of the
particulate flows (equations (1a)-(1g) and the $\phi$-FEM scheme (11)). These
schemes will be compared with standard FEM on fitted triangular meshes as on
Fig. 3 left (we do not introduce higher order approximations of the
curvilinear boundary of the domain, as would be the case in the isoparametric
FEM for example). In the case of Stokes equations, the error is measured with
respect to a manufactured solution, while a reference solution obtained by
standard FEM on a fitted fine mesh is used in the case of particulate flows.
We have implemented $\phi$-FEM in multiphenics [4]. The implementation scripts
can be consulted on GitHub.666https://github.com/michelduprez/phi-FEM-
particulate-flows-Stokes.git
or https://doi.org/10.5281/zenodo.6817135

The fluid/solid domain in both our test cases is
$\mathcal{O}=(0,1)^{2}\subset\mathbb{R}^{2}$ and we take the particle
$\mathcal{S}$ as a disk of radius $R=0.21$ centered at a point $(0.5,0.5)$.
Then $\Omega=[0,1]^{2}\backslash\mathcal{S}$. The geometry is presented in
Fig. 1 (left). In $\phi$-FEM, we use the following level-set function, well
defined and smooth for all $(x,y)\in\mathbb{R}^{2}$,

$\phi(x,y)=R^{2}-(x-0.5)^{2}-(y-0.5)^{2}.$ (66)

We present only the results with the lowest order Taylor-Hood elements, i.e.
setting $k=2$ and thus using $P_{2}$ elements for $w_{h}$, $\phi_{h}$,
$\chi_{h}$ (the approximation $\phi_{h}$ for $\phi$ is exact in this case).
The tests with elements of higher order would lead to essentially the same
observations. The stabilization parameters are set to $\sigma=\sigma_{u}=20$
(as in [12]).

###  4.1 Particular case of a fixed particle: Stokes equations

We start by Stokes equations (71) in the domain $\Omega$, as above, with the
right-hand side such that the exact solution is as follows, cf. [14],

$\displaystyle u(x,y)$ $\displaystyle=(\cos(\pi x)\sin(\pi y),-\sin(\pi
x)\cos(\pi y)),$ $\displaystyle p(x,y)$ $\displaystyle=(y-0.5)\cos(2\pi
x)+(x-0.5)\sin(2\pi y),$

taking $\nu=1$.

Figure 3: Mesh used for the standard FEM formulation (left) and mesh used in
the $\phi$-FEM schemes (right).

We shall test the $\phi$-FEM scheme given by (72) in Appendix A and compare it
with a standard Taylor-Hood FEM on a fitted mesh. To this end, we introduce a
quasi-uniform triangular mesh $\mathcal{T}_{h}^{fit}$ fitted to $\Omega$, in
the sense that the boundary nodes of the mesh lie on $\Gamma\cup\Gamma_{w}$.
The domain occupied by this mesh, denoted by $\Omega_{h}^{fit}$, is a
polygonal approximation of $\Omega$, since the interface $\Gamma$ is
curvilinear and cannot be represented exactly by the straight edges. We
introduce then the FE spaces

$\mathcal{V}_{h}^{fit,u_{D}}=\left\\{v_{h}\in
C(\bar{\Omega}_{h}^{fit})^{d}:v_{h}|_{T}\in\mathbb{P}^{2}(T)^{d}\hskip
10.00002pt\forall T\in\mathcal{T}_{h}^{fit},\hskip
10.00002ptv_{h}=I_{h}u_{D}\text{ on
}\Gamma_{w}\cap\Gamma_{h}^{fit}\right\\}\,$ (67)
$\mathcal{M}_{h}^{fit}=\left\\{q_{h}\in
C(\bar{\Omega}_{h}^{fit}):q_{h}|_{T}\in\mathbb{P}^{1}(T)\hskip
10.00002pt\forall T\in\mathcal{T}_{h}^{fit},\hskip
10.00002pt\int_{\Omega_{h}^{fit}}q_{h}=0\right\\},$ (68)

where $\Gamma_{h}^{fit}$ is the part of the boundary of $\Omega_{h}^{fit}$
approximating $\Gamma$. A standard fitted Taylor-Hood FEM can be written as:
find $(u_{h},p_{h})\in\mathcal{V}_{h}^{fit,u_{D}}\times\mathcal{M}_{h}^{fit}$
such that

$\int_{\Omega_{h}^{fit}}2D(u_{h}):D(v_{h})-\int_{\Omega_{h}^{fit}}p_{h}\operatorname{div}v_{h}-\int_{\Omega_{h}^{fit}}q_{h}\operatorname{div}u_{h}=\int_{\Omega}fv_{h},$
(69)

for all $(v_{h},q_{h})\in\mathcal{V}_{h}^{fit,0}\times\mathcal{M}_{h}^{fit}$.

Examples of meshes for the standard FEM formulation and the $\phi$-FEM scheme
are given in Fig. 3 (left) and (right), respectively. In Fig. 4, we report
convergence results for the standard fitted Taylor-Hood FEM (69) and for
$\phi$-FEM (72) in the case of Stokes equations. We recover the theoretical
rates of convergence of $\phi$-FEM stated in Theorem 2: 2nd order in $h$ for
both the $H^{1}$-error in velocity and the $L^{2}$-error in pressure. The
$L^{2}$-error in velocity is actually better than theoretically predicted: it
is of order 3 instead of theoretically expected 2.5. We observe thus that
$\phi$-FEM is fully optimal in practice: it demonstrates the same convergence
rates in all the aforementioned norms as the standard FEM would demonstrate in
the ideal situation of a fitted mesh on a convex polygonal domain. Actually,
our setting is not ideal, $\Omega$ is neither convex, not polygonal. It is
thus not surprising that the standard Taylor-Hood FEM underperforms (we recall
that the mesh $\mathcal{T}_{h}^{fit}$ is triangular with straight edges and no
higher order geometrical approximation is introduced there). Experimentally
observed convergence rates for this scheme are $\approx 2$ for the
$L^{2}$-error in velocity, $\approx 1.5$ for the $H^{1}$-error in velocity,
and slightly smaller than 2 for the $L^{2}$-error in pressure (the error in
pressure is thus the only one for which the optimal convergence order seems to
be retained in practice). Moreover, all the errors on all the considered
meshes are systematically smaller for $\phi$-FEM than for the fitted FEM.

######  Remark 6.

As already noted in Remark 2, in $\phi$-FEM, it is impossible to impose
$\int_{\Omega}p_{h}=0$. In our implementation, we rather impose
$\int_{\Omega_{h}}p_{h}=0$ with the help of a Lagrange multiplier, i.e. we add
$\lambda_{h}\int_{\Omega_{h}}p_{h}+\mu_{h}\int_{\Omega_{h}}q_{h}$ (with
$\lambda_{h},\mu_{h}\in\mathbb{R}$) to the formulation. To compute the
relative $L^{2}$-error for the pressure, we should compare $p_{h}$ with the
exact pressure whose integral vanishes over $\Omega_{h}$ (recall that the
pressure is physically defined up to an additive constant any way). We thus
introduce $\tilde{p}=p-c_{\Omega_{h}}$ with
$c_{\Omega_{h}}=\frac{1}{|\Omega_{h}|}\int_{\Omega_{h}}p$ and compute the
errors with respect to $\tilde{p}$. Similarly, in the case of standard fitted
FEM, we impose $\int_{\Omega_{h}^{fit}}p_{h}=0$ by a Lagrange multiplier and
compute the relative error against $\tilde{p}=p-c_{\Omega_{h}^{fit}}$ with
$c_{\Omega_{h}^{fit}}=\frac{1}{|\Omega_{h}^{fit}|}\int_{\Omega_{h}^{fit}}p$.

$10^{-2}$$10^{-1}$$10^{-7}$$10^{-6}$$10^{-5}$$10^{-4}$$10^{-3}$1213$h$$\|u-u_{h}\|_{0,\Omega_{h}^{fit}}/\|u\|_{0,\Omega_{h}^{fit}}$
st.-FEM$\|u-u_{h}\|_{0,\Omega_{h}}/\|u\|_{0,\Omega_{h}}$ $\phi$-FEM

$10^{-2}$$10^{-1}$$10^{-5}$$10^{-4}$$10^{-3}$$10^{-2}$12$h$$\|u-u_{h}\|_{1,\Omega_{h}^{fit}}/\|u\|_{1,\Omega_{h}^{fit}}$
st.-FEM$\|u-u_{h}\|_{1,\Omega_{h}}/\|u\|_{1,\Omega_{h}}$ $\phi$-FEM

$10^{-2}$$10^{-1}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$12$h$$\|\tilde{p}-p_{h}\|_{0,\Omega_{h}^{fit}}/\|^{fit}p\|_{0,\Omega_{h}^{fit}}$
st.-FEM$\|\tilde{p}-p_{h}\|_{0,\Omega_{h}}/\|\tilde{p}\|_{0,\Omega_{h}}$
$\phi$-FEM

Figure 4: Rates of convergence for the standard Taylor-Hood FEM scheme (69)
and the $\phi$-FEM scheme (72) in the case of Stokes equations. The $L^{2}$
relative error of the velocity (left), the $H^{1}$ relative error of the
velocity (middle) and the $L^{2}$ relative error of the pressure (right).

###  4.2 Particulate flows

We now turn to the creeping particulate flow equations (1a)-(1g) in the same
geometry as before. The level-set $\phi$ is again defined by (66). The
vertical gravity is assumed to be equal to $10$. The density of the fluid and
the solid are equal to $\rho_{f}=1$ and $\rho_{s}=2$, respectively, and the
viscosity $\nu=1$. We deduce that the mass of the solid is equal to
$m=\rho_{s}\pi^{2}R^{2}$. For the cut-off $\chi$, we consider the radial
polynomial of degree $5$ on the interval $(r_{0},r_{1})$ with $r_{0}=R$ and
$r_{1}=0.45$ such that $\chi(r_{0})=1$ and
$\chi^{\prime}(r_{0})=\chi^{\prime\prime}(r_{0})=\chi(r_{1})=\chi^{\prime}(r_{1})=\chi^{\prime\prime}(r_{1})=0$
so that, setting $\chi=1$ for $r<R$ and $\chi=0$ for $r>R$, the resulting
$\chi$ is two times differentiable and thus $\chi\in H^{3}(\mathcal{O})$ as
required by Assumption 2:

$\chi(r)=\begin{cases}1,&\text{ for }r<r_{0}\\\
1+\frac{f(r_{0},r_{1})}{(r_{1}-r_{0})^{5}},&\text{ for }r_{0}<r<r_{1}\\\
0,&\text{ for }r>r_{1}\\\ \end{cases}$

where

$f(r_{0},r_{1})=(-6r^{5}+15(r_{0}+r_{1})r^{4}-10(r_{0}^{2}+4r_{0}r_{1}+r_{1}^{2})r^{3}+30r_{0}r_{1}(r_{0}+r_{1})r^{2}\\\
-30r_{0}^{2}r_{1}^{2}r+r_{0}^{3}(r_{0}^{2}-5r_{1}r_{0}+10r_{1}^{2})).$

Let us introduce a Taylor-Hood scheme which will be compared with our
$\phi$-FEM scheme. We introduce first the fitting mesh $\mathcal{T}_{h}^{fit}$
on domain $\Omega_{h}^{fit}$ as in the preceding section, cf. Fig. 3 (left),
and adapt the Taylor-Hood FE space velocity space, cf. (67), as

$\mathcal{V}_{h}^{fit}=\left\\{v_{h}\in
C(\bar{\Omega}_{h}^{fit})^{d}:v_{h}|_{T}\in\mathbb{P}^{2}(T)^{d}\hskip
10.00002pt\forall T\in\mathcal{T}_{h}^{fit},\hskip 10.00002ptv_{h}=0\text{ on
}\Gamma_{w}\right\\},$

while keeping the pressure space (68) as before. Note that the velocity FE
space does no longer contain any restrictions on the boundary part
$\Gamma_{h}^{fit}$ approximating $\Gamma$. We shall impose the boundary
conditions there with the help of Lagrange multipliers, introducing the space

$\Lambda_{h}^{fit}=\left\\{\mu_{h}\in
C(\bar{\Omega}_{h}^{fit}):\mu_{h}|_{F}\in\mathbb{P}^{2}(F)\hskip
10.00002pt\forall F\in\mathcal{F}_{h}^{fit}\right\\},$

where $\mathcal{F}_{h}^{fit}$ is the set of the boundary facets on
$\Gamma_{h}^{fit}$. A fitted Taylor-Hood FE formulation is then written as:
find
$(u_{h},p_{h},\lambda_{h},U_{h},\psi_{h})\in\mathcal{V}_{h}^{fit}\times\mathcal{M}_{h}^{fit}\times\Lambda_{h}^{fit}\times\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}$
such that

$\int_{\Omega_{h}^{fit}}2D(u_{h}):D(v_{h})-\int_{\Omega_{h}^{fit}}p_{h}\operatorname{div}v_{h}-\int_{\Omega_{h}^{fit}}q_{h}\operatorname{div}u_{h}\\\
+\int_{\Gamma_{h}^{fit}}\lambda_{h}\cdot(v_{h}-V_{h}-\omega_{h}\times
r)+\int_{\Gamma_{h}^{fit}}\mu_{h}\cdot(u_{h}-U_{h}-\psi_{h}\times
r)={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\int_{\Omega_{h}^{fit}}\rho_{f}gv_{h}}+mg\cdot
V_{h},$ (70)

for all
$(v_{h},q_{h},\mu_{h},V_{h},\omega_{h})\in\mathcal{V}_{h}^{fit}\times\mathcal{M}_{h}^{fit}\times\Lambda_{h}^{fit}\times\mathbb{R}^{d}\times\mathbb{R}^{d^{\prime}}$.
We present in Fig. 5 the velocity obtained with the standard Taylor-Hood FEM
scheme (70). Such a velocity and the accompanying pressure, computed on a very
fine fitted grid, will be used as the reference solution in the subsequent
numerical experiments and will be denoted as $u,p$ in what follows.

Figure 5: Velocity obtained with the standard Taylor-Hood FEM scheme (70).

A comparison between the standard Taylor-Hood FEM (70) and $\phi$-FEM (11) is
presented in Figs. 6 and 7 (we do not report the error in the particle
rotation velocity whose exact value is 0 and which is accurately predicted by
all the schemes up to machine precision; this can be attributed to the
symmetry of our test case). Since the error is computed with respect to a fine
reference solution solution living on a fine fitted mesh, the numerical
solution computed by either (70) or (11) should be projected to this fine mesh
in order to compute the errors. This is reflected in the legends of the
convergence curves: $\Omega^{fine}$ stands for the fine fitted approximation
there. Similarly to Remark 6, we impose the pressure mean by a Lagrange
multiplier in both schemes, and adjust the additive constants properly when
computing the errors in pressure.

The conclusions are essentially the same as in the previous test case (Stokes
equations alone): $\phi$-FEM exhibits optimal convergence rates, while the
fitted standard FEM is suboptimal (with the exception of the $L^{2}$ error in
pressure). It seems again that our theoretical estimates for the $L^{2}$-error
of the fluid velocity is not sharp: the experimental convergence rate is $k+1$
rather than $k+\frac{1}{2}$. The same observation can be made about the
particle velocity: the experimental convergence rate is $k+1$ rather than
theoretically predicted $k$.

$10^{-1.8}$$10^{-1.6}$$10^{-1.4}$$10^{-1.2}$$10^{-1}$$10^{-0.8}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$$10^{0}$1213$h$$\|u-u_{h}\|_{0,\Omega^{fine}}/\|u\|_{0,\Omega^{fine}}$
st.-FEM$\|u-u_{h}\|_{0,\Omega^{fine}}/\|u\|_{0,\Omega^{fine}}$ $\phi$-FEM

$10^{-1.8}$$10^{-1.6}$$10^{-1.4}$$10^{-1.2}$$10^{-1}$$10^{-0.8}$$10^{-2}$$10^{-1}$$10^{0}$1211.5$h$$\|u-u_{h}\|_{1,\Omega^{fine}}/\|u\|_{1,\Omega^{fine}}$
st.-FEM$\|u-u_{h}\|_{1,\Omega^{fine}}/\|u\|_{1,\Omega^{fine}}$ $\phi$-FEM

Figure 6: Rates of convergence for the standard Taylor-Hood FEM scheme (70)
and the $\phi$-FEM scheme (11) in the case of particulate flows. The $L^{2}$
relative error of the velocity (left) and the $H^{1}$ relative error of the
velocity (right).

$10^{-1.8}$$10^{-1.6}$$10^{-1.4}$$10^{-1.2}$$10^{-1}$$10^{-0.8}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$12$h$$\|p-p_{h}\|_{0,\Omega^{fine}}/\|p\|_{0,\Omega^{fine}}$
st.-FEM$\|p-p_{h}\|_{0,\Omega^{fine}}/\|p\|_{0,\Omega^{fine}}$ $\phi$-FEM

$10^{-1.8}$$10^{-1.6}$$10^{-1.4}$$10^{-1.2}$$10^{-1}$$10^{-0.8}$$10^{-4}$$10^{-3}$$10^{-2}$$10^{-1}$1213$h$$|U-U_{h}|/|U|$
st.-FEM$|U-U_{h}|/|U|$ $\phi$-FEM

Figure 7: Rates of convergence for the standard Taylor-Hood FEM scheme (69)
and the $\phi$-FEM scheme (72) in the case of particulate flows. The $L^{2}$
relative error of the pressure (left) and relative error of the displacement
of the solid (right).

##  5 Conclusions and perspectives

The main goal of the $\phi$-FEM approach is to avoid the non standard
quadrature on the cut mesh cells inherent to CutFEM. Some of the attractive
features of $\phi$-FEM are:

  * •

$\phi$-FEM is readily available for finite elements of any order (without the
need of any additional higher order approximation of the geometry).

  * •

(corollary of the previous point) $\phi$-FEM outperforms the standard fitted
FEM on polygonal meshes if the order of piecewise polynomials is $>1$ (unless
one implements more complicated versions of FEM in order to treat properly the
curvilinear boundary, such as the isoparametric FEM).

  * •

$\phi$-FEM uses standard FE spaces and is based on a variational formulation
of the problem, so that it can be easily implemented in existing general-
purpose FEM libraries, provided they allow one to compute the jumps on
selected facets and the second derivatives on selected cells.

To counter-balance the last point, we should mention that the variational
formulation at the base of $\phi$-FEM can be rather complicated. Typically, it
contains more terms than a fitted FEM or a CutFEM scheme for the same problem.
The implementation, although straightforward, may thus require some extra
programming efforts. As a reward, one easily gets a good description of the
geometry which may result in smaller computing times than those for standard
FEM, as highlighted in [10].

In the present article we have adapted $\phi$-FEM to the Stokes equations and
to the combination of the Stokes equations with the motion of a rigid particle
inside the fluid, providing a first brick in future applications of this
technology to fluid structure interaction problems.

Of course, there remains a lot of open questions. To conclude, we list some of
those (not necessarily related to the particular case of Stokes equations or
flows with particles) and envisage potential solutions:

  * •

Theoretical suboptimality of $\phi$-FEM in the $L^{2}$-norm. However, all the
numerical experiments show the optimal convergence in this norm, which gives
us hope that a sharper theoretical estimate could be found.

  * •

A mismatch between the theoretical construction of the active mesh
$\mathcal{T}_{h}$ and its practical implementation, cf. Remark 1. The
theoretical construction of $\mathcal{T}_{h}$ is crucial for the current proof
of coercivity, cf. Lemma 9. An alternative proof should be found.

  * •

A practical construction of the levelset function, which should satisfy some
assumptions summarized in Remark 3. In the present article, $\phi$ was given
analytically, but in more realistic applications one will have to construct an
appropriate $\phi_{h}$ on the discrete level directly. A good candidate, in
the vicinity of $\Gamma$ would be the signed distance to $\Gamma$, for which
there exist efficient libraries, but it will remain to extend it in a smooth
way (again directly on the discrete level) far from $\Gamma$. As an
alternative, we note that in other versions of $\phi$-FEM, $\phi_{h}$ may be
required only on mesh cells near $\Gamma$, cf. [11, 10].

  * •

The current construction of stabilization in $\phi$-FEM relies heavily on the
linearity of the governing equations. Indeed, the terms with prefactor
$\sigma$ in (11) reproduce the governing equations for both trial and test
functions. If the equations are non-linear, one cannot do this since the
formulation should remain linear in the test functions. Thus, going from
Stokes to Navier-Stokes, for example, is not straightforward. Various options
of linearization of the stabilization terms should be yet tested numerically
and theoretically.

##  Appendix A $\phi$-FEM for Stokes equations in a fixed domain

In this section, we propose a $\phi$-FEM scheme for the simpler case of a
fixed solid in the fluid. The governing equations are the non-homogeneous
Stokes equations given by

$\left\\{\begin{array}[]{ll}-2\operatorname{div}(D(u))+\nabla p=f,&\text{in
}\Omega,\\\ \operatorname{div}u=0,&\text{in }\Omega,\\\ u=u_{D},&\text{on
}\Gamma\cup\Gamma_{w}.\end{array}\right.$ (71)

Assume that $u_{D}$ and $f$ are defined in the whole discrete domain
$\Omega_{h}$. Inspired by the $\phi$-FEM scheme for particulate flow given in
(11), we can derive the following $\phi$-FEM scheme for the non-homogeneous
Stokes equations (71): find $w_{h}\in\mathcal{V}_{h}$,
$p_{h}\in\mathcal{M}_{h}$ satisfying

$\int_{\Omega_{h}}2D(u_{D}+\phi_{h}w_{h}):D(\phi_{h}s_{h})-\int_{G_{h}}(2D(u_{D}+\phi_{h}w_{h})-p_{h}I)n\cdot\phi_{h}s_{h}\text{}$
(72)
$-\int_{\Omega_{h}}p_{h}\operatorname{div}(\phi_{h}s_{h})-\int_{\Omega_{h}}q_{h}\operatorname{div}(u_{D}+\phi_{h}w_{h})+\sigma_{u}J_{u}(u_{D}+\phi_{h}w_{h},\phi_{h}s_{h})$
$+\sigma
h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}(-\Delta(u_{D}+\phi_{h}w_{h})+\nabla
p_{h})\cdot(-\Delta(\phi_{h}s_{h})-\nabla q_{h})$
$+\sigma\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}\operatorname{div}(u_{D}+\phi_{h}w_{h})\operatorname{div}(\phi_{h}v_{h})$
$=\int_{\Omega_{h}}f\phi_{h}s_{h}+\sigma
h^{2}\sum_{T\in\mathcal{T}_{h}^{\Gamma}}\int_{T}f(-\Delta(u_{D}+\phi_{h}s_{h})-\nabla
q_{h}),$

for all $s_{h}\in\mathcal{V}_{h}$, $q_{h}\in\mathcal{M}_{h}$.

We now state our second main result for the Stokes equations:

######  Theorem 2.

Suppose that Assumptions 1, 3, 4 and 5 hold true, the mesh $\mathcal{T}_{h}$
is quasi-uniform. Let $(u,p)\in H^{k+1}(\Omega)^{d}\times H^{k}(\Omega)$ be
the solution to (71) and
$(w_{h},p_{h})\in\mathcal{V}_{h}\times\mathcal{M}_{h}$ be the solution to
(72). Denoting $u_{h}:=\phi_{h}w_{h}$, it holds

$|u-u_{h}|_{1,\Omega\cap\Omega_{h}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{\nu}}|p-p_{h}|_{0,\Omega\cap\Omega_{h}}\leq
Ch^{k}(\|u\|_{k+1,\Omega}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{1}{\nu}}\|p\|_{k,\Omega})$

with a constant $C>0$ depending on the $C_{0}$, $m$, $M$ in Assumptions 1, 4,
on the maximum of the derivatives of $\phi$, on the mesh regularity, and on
the polynomial degree $k$, but independent of $h$, $f$, and $u$.
Moreover, supposing $\Omega\subset\Omega_{h}$

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\|u-u_{h}\|_{0,\Omega}\leq
Ch^{k+1/2}(\|u\|_{k+1,\Omega}+\frac{1}{\nu}\|p\|_{k,\Omega})}$

with a constant $C>0$ of the same type.

The proof of Theorem 2 can be adapted from the proof of Theorem 1. It is even
more simple.

##  Appendix B A glossary of geometrical notations.

$\Gamma_{h}$ |  the approximate interface: $\Gamma_{h}=\\{\phi_{h}=0\\}$
---|---
$\mathcal{T}_{h}$ |  the active mesh: $\mathcal{T}_{h}=\\{T\in\mathcal{T}_{h}^{\mathcal{O}}:T\cap\\{\phi_{h}<0\\}\neq\varnothing\\}$
$\mathcal{T}_{h}^{\Gamma}$ |  intersection of $\mathcal{T}_{h}$ with $\Gamma_{h}$ $\mathcal{T}_{h}^{\Gamma}=\\{T\in\mathcal{T}_{h}:T\cap\Gamma_{h}\neq\varnothing\\}$
$\mathcal{T}_{h}^{\Gamma,ext}$ |  $\mathcal{T}_{h}^{\Gamma}$ and the cells which are neighbors and neighbors of neighbors of cells
|  of $\mathcal{T}_{h}^{\Gamma}$ in $\mathcal{T}_{h}$
$\mathcal{T}_{h}^{i}$ | $\mathcal{T}_{h}\setminus\mathcal{T}_{h}^{\Gamma}$
$\mathcal{F}_{h}^{\Gamma}$ | $\mathcal{F}_{h}^{\Gamma}=\\{E\text{ (an internal facet of }\mathcal{T}_{h})\text{ such that }\exists T\in\mathcal{T}_{h}^{\Gamma}\text{ and }E\in\partial T\\}.$
$G_{h}$ |  the internal component of $\partial\Omega_{h}$, corresponding to the interface $\Gamma$: $G_{h}=\partial\Omega_{h}\setminus\Gamma_{w}$
$B_{h}$ |  the strip between $\Gamma$ and $G_{h}$: $B_{h}=\Omega_{h}\setminus\Omega$
$B_{h}^{\Gamma}$ |  the strip between $\Gamma$ and $\Gamma_{h}$
$\tilde{B}_{h}$ |  the strip between $\Gamma_{h}=\\{\phi_{h}=0\\}$ and $G_{h}$: $\tilde{B}_{h}=\\{\phi_{h}>0\\}\cap\Omega_{h}$

We also recall that the domain occupied by the active mesh $\mathcal{T}_{h}$
is denoted by $\Omega_{h}$, i.e.
$\Omega_{h}:=\left(\cup_{T\in\mathcal{T}_{h}}T\right)^{o}$. The same
convention is applied to the submeshes $\mathcal{T}_{h}^{\Gamma}$ and
$\mathcal{T}_{h}^{i}$, giving respectively $\Omega_{h}^{\Gamma}$ and
$\Omega_{h}^{i}$.

## References

  * [1] S. Amdouni, K. Mansouri, Y. Renard, M. Arfaoui, and M. Moakher.  Numerical convergence and stability of mixed formulation with x-fem cut-off.  European Journal of Computational Mechanics, 21(3-6):160–173, 2012\.
  * [2] P. Angot, C.-H. Bruneau, and P. Fabrie.  A penalization method to take into account obstacles in incompressible viscous flows.  Numerische Mathematik, 81(4):497–520, 1999.
  * [3] N. M. Atallah, C. Canuto, and G. Scovazzi.  Analysis of the shifted boundary method for the Stokes problem.  Comput. Methods Appl. Mech. Engrg., 358:112609, 33, 2020.
  * [4] F. Ballarin and G. Rozza.  multiphenics.  https://mathlab.sissa.it/multiphenics, 2020.
  * [5] D. Boffi, F. Brezzi, and M. Fortin.  Finite elements for the Stokes problem., volume 1939 of Lecture Notes in Mathematics.  Springer-Verlag, Berlin; Fondazione C.I.M.E., Florence, 2008.  Mixed finite elements, compatibility conditions, and applications. Lectures given at the C.I.M.E. Summer School held in Cetraro, June 26–July 1, 2006, Edited by Daniele Boffi and Lucia Gastaldi.
  * [6] S. C. Brenner and L. R. Scott.  The mathematical theory of finite element methods, volume 3.  Springer, 2008.
  * [7] E. Burman.  Ghost penalty.  C. R. Math. Acad. Sci. Paris, 348(21-22):1217–1220, 2010.
  * [8] E. Burman and P. Hansbo.  Fictitious domain methods using cut elements: III. A stabilized Nitsche method for Stokes’ problem.  ESAIM Math. Model. Numer. Anal., 48(3):859–874, 2014.
  * [9] P. G. Ciarlet.  Mathematical Elasticity: Volume I: three-dimensional elasticity.  North-Holland, 1988.
  * [10] S. Cotin, M. Duprez, V. Lleras, A. Lozinski, and K. Vuillemot.  $\phi$-fem: an efficient simulation tool using simple meshes for problems in structure mechanics and heat transfer.  In S. Bordas and A. Menk, editors, Partition of Unity Methods (Wiley Series in Computational Mechanics) 1st Edition. Wiley; 1st edition, 2022\.
  * [11] M. Duprez, V. Lleras, and A. Lozinski.  A new $\phi$-FEM approach for problems with natural boundary conditions.  NMPDE, accepted, 2021.
  * [12] M. Duprez and A. Lozinski.  $\phi$-FEM: a finite element method on domains defined by level-sets.  SIAM J. Numer. Anal., 58(2):1008–1028, 2020.
  * [13] A. Ern and J.-L. Guermond.  Theory and practice of finite elements, volume 159.  Springer Science & Business Media, 2013.
  * [14] M. Fournié and A. Lozinski.  Stability and optimal convergence of unfitted extended finite element methods with Lagrange multipliers for the Stokes equations.  In Geometrically unfitted finite element methods and applications, volume 121 of Lect. Notes Comput. Sci. Eng., pages 143–182. Springer, Cham, 2017.
  * [15] R. Glowinski, T.-W. Pan, T. I. Hesla, and D. D. Joseph.  A distributed Lagrange multiplier/fictitious domain method for particulate flows.  International Journal of Multiphase Flow, 25(5):755–794, 1999.
  * [16] R. Glowinski, T.-W. Pan, T. I. Hesla, D. D. Joseph, and J. Periaux.  A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow.  Journal of computational physics, 169(2):363–426, 2001.
  * [17] J. Guzmán and M. Olshanskii.  Inf-sup stability of geometrically unfitted Stokes finite elements.  Math. Comp., 87(313):2091–2112, 2018.
  * [18] M.-C. Lai and C. S. Peskin.  An immersed boundary method with formal second-order accuracy and reduced numerical viscosity.  Journal of computational Physics, 160(2):705–719, 2000.
  * [19] A. Main and G. Scovazzi.  The shifted boundary method for embedded domain computations. Part I: Poisson and Stokes problems.  J. Comput. Phys., 372:972–995, 2018.
  * [20] A. Massing, M. G. Larson, A. Logg, and M. E. Rognes.  A stabilized Nitsche fictitious domain method for the Stokes problem.  J. Sci. Comput., 61(3):604–628, 2014.
  * [21] R. Mittal and G. Iaccarino.  Immersed boundary methods.  Annu. Rev. Fluid Mech., 37:239–261, 2005.
  * [22] C. S. Peskin.  Numerical analysis of blood flow in the heart.  Journal of computational physics, 25(3):220–252, 1977.