# Low-order fictitious domain method with enhanced mass conservation for an interface Stokes problem

Daniele Corti, Guillaume Delay, Miguel Angel Fernández, Fabien Vergnet, Marina Vidrascu

► To cite this version:

Daniele Corti, Guillaume Delay, Miguel Angel Fernández, Fabien Vergnet, Marina Vidrascu. Low-order fictitious domain method with enhanced mass conservation for an interface Stokes problem. ESAIM: Mathematical Modelling and Numerical Analysis, 2024, 58 (1), pp.303-333. 10.1051/m2an/2023103. hal-04084162v2

HAL Id: hal-04084162

https://inria.hal.science/hal-04084162v2

Submitted on 18 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# LOW-ORDER FICTITIOUS DOMAIN METHOD WITH ENHANCED MASS CONSERVATION FOR AN INTERFACE STOKES PROBLEM*

DANIELE C. CORTI<fnref n="1" />, GUILLAUME DELAY<fnref n="1" />, MIGUEL A. FERNÁNDEZ<fnref n="1" />, FABIEN VERGNET<fnref n="1" />
AND MARINA VIDRASCU<fnref n="1" />

**Abstract.** One of the main difficulties that has to be faced with fictitious domain approximation of incompressible flows with immersed interfaces is related to the potential lack of mass conservation across the interfaces. In this paper, we propose and analyze a low order fictitious domain stabilized finite element method which mitigates this issue with the addition of a single velocity constraint. We provide a complete a priori numerical analysis of the method under minimal regularity assumptions. A comprehensive numerical study illustrates the capabilities of the proposed method, including comparisons with alternative fitted and unfitted mesh methods.

**Résumé.** Il est connu que les méthodes de type domaine fictif pour l'approximation d'écoulements de fluides incompressibles en présence d'interfaces immergées présentent un défaut de conservation de masse à travers les interfaces. Dans cet article, nous proposons et analysons une telle méthode basée sur des éléments finis d'ordre bas et stabilisée qui atténue ce problème en ajoutant une seule contrainte scolaire sur la vitesse. Nous réalisons l'analyse numérique a priori de la méthode sous des hypothèses de régularité minimales. Une étude numérique complète illustre les capacités de la méthode présentée, en proposant notamment des comparaisons avec d'autres méthodes de la littérature.

**2020 Mathematics Subject Classification.** 65N15, 65N30, 65N85.

December 18, 2023.

## 1. INTRODUCTION

The simulation of incompressible flows with immersed moving interfaces plays a fundamental role in a wide variety of engineering fields: from the biomechanics of heart valves to the aeroelasticity of parachutes (see, e.g., [43, 44, 51, 53, 55]). The spatial approximation of this type of problems typically falls into one of the following two categories: fitted and unfitted mesh numerical methods.

In fitted mesh methods, the computational mesh has an explicit representation of the immersed interface. This facilitates the enforcement of the interface conditions and, furthermore, (weak and strong) discontinuities of the solution can be straightforwardly incorporated at the discrete level, yielding an optimally accurate method.

<footnote>*Keywords and phrases: Stokes problem, fictitious domain method, Lagrange multiplier, interfacial mass conservation.</footnote>

<footnote n="1">¹D.C. Corti, M.A. Fernández and M. Vidrascu were partially supported by the French National Research Agency (ANR), through the SIMR project (ANR-19-CE45-0020).</footnote>

<footnote n="1">¹ Sorbonne Université, CNRS, Laboratoire Jacques-Louis Lions (LJLL), Inria, Paris, France
e-mail: daniele.corti@inria.fr, guillaume.delay@sorbonne-universite.fr, miguel.fernandez@inria.fr, fabien.vergnet@sorbonne-universite.fr, marina.vidrascu@inria.fr</footnote>
Nevertheless, the body-fitted nature of the mesh can become cumbersome for large interface deflections, due to highly distorted meshes which may require remeshing or topological modifications of the mesh (see, e.g., [1,52]) and thus increase the computational cost.

Unfitted mesh methods are a widespread approach to avoid these issues. In these approaches, the immersed interface is, by construction, free to move independently of the background computational mesh. This flexibility comes however at a price: the lack of interfacial representation within the computational mesh can lead to accuracy issues. Among these types of methods, the immersed boundary and fictitious domain methods (see, e.g., [9,12,33,47,50]) are ones of the most widespread, due to their simplicity of implementation. Basically, the idea consists in enforcing the interface Dirichlet condition as an additional constraint, using Lagrange multipliers, or penalization. In general, these methods have the reputation of being sub-optimally accurate and of suffering from interfacial mass conservation issues. This is a consequence of the fact that the discrete velocity and pressure do not allow for weak and strong discontinuities, respectively, across the interface. Different approaches have been proposed in the literature to circumvent these two different issues.

In the fictitious domain control approach proposed in [4,25,54], optimal accuracy is obtained by constructing a smooth extension which removes the weak discontinuity in the velocity. The method introduces, however, additional unknowns and lacks complete numerical analysis. Cut-FEM approaches achieve optimal accuracy by integrating the equations only in the physical region and by adding suitable stabilization terms for robustness (see, e.g., [15,35,37]), but they require a more involved computer implementation due to the specific tracking of the interface intersections and quadrature over arbitrary polygons. Recent approaches, such as the shifted boundary method or the $\phi$-FEM method (see, e.g., [22,23,41]) avoid these computational geometry difficulties, but their interface formulation is not straightforward (due to the introduction of additional unknowns or the need of a level-set description of the interface).

Penalized grad-div interfacial stabilization is known to enhance mass conservation across the interface (see, e.g., [13,17,30,39,42]), but at the price of degradation of the system matrix conditioning, which drastically limits the applicability of the method. The extended finite element method (XFEM) is an elegant way of guaranteeing interfacial mass conservation by introducing strong pressure discontinuities across the interface at the discrete level (see, e.g., [2,8,40,48,56]), but this requires the integration over cut-elements and the number of degrees of freedom (and hence the size of the system matrix) may change between different interface locations. Enhancing mass conservation in fictitious domain methods by using globally discontinuous pressure is an alternative to XFEM, but inf-sup stability prevents the use of low order elements for the velocity, unless stabilization terms are introduced that would compromise local mass conservation (see, e.g., [5,10,11]). Similar observations can be made on the combination of unfitted mesh methods and divergence free approximations (see, e.g., [18,42]), with the exception of the cut-FEM method reported in [16] based on the minimal divergence free element introduced in [19].

Numerical evidence (see, e.g., [13]) suggests that, in practice, the inaccuracies of traditional fictitious domains methods for incompressible flows are mainly driven by the artificial interfacial mass loss, rather than by the sub-optimality of the convergence rate. In this paper, we introduce a new low order fictitious domain finite element method based on continuous piecewise affine approximations. We consider as model problem the Stokes equations with an immersed Dirichlet interface condition. At the discrete level, the interface condition is enforced via Lagrange multipliers and mass conservation across the interface is enhanced by enforcing an additional velocity constraint. This constraint can be interpreted as an enrichment of the pressure approximation space with a single Heaviside function across the interface. Similar ideas have been advocated in [38], but for a different fictitious domain method (involving discrete Dirac delta functions) and without error analysis. It is also worth mentioning the work reported in [27,46], where interface flux constraints are introduced in the context of bi-fluid simulations, but therein the kinematic continuity can directly be embedded within the discrete space. Suitable stabilization terms are introduced to guarantee inf-sup stability (Lemma 4.8) and a priori error estimates for the velocity and pressure are also derived (Theorem 4.12). These error estimates are robust with respect to the magnitude of the pressure jump across the interface, due to the additional basis function. Finally, a comprehensive numerical study provides numerical evidence on the error analysis and illustrates the
capabilities of the proposed method. This method is compared with different fitted and unfitted mesh methods in two benchmarks. The numerical results show, in particular, that the present approach provides a simple and robust low order method that mitigates the accuracy issues of traditional fictitious domain methods for incompressible flows with immersed interfaces.

The rest of the paper is organized as follows. The continuous setting is introduced in Section 2. Section 3 presents the proposed numerical method. A priori stability and error estimates of the method are derived in Section 4. Section 5 is devoted to the numerical investigations. Finally, some conclusions and lines of future work are drawn in Section 6.

## 2. PROBLEM SETTING

Let $\Omega \subset \mathbb{R}^d$, $d=2$ or $3$, be a bounded polygonal/polyhedral domain and $\Sigma$ be a polygonal oriented manifold of co-dimension 1 fully immersed in $\Omega$ and with unit normal vector $n_{\Sigma}$. We assume that $\Sigma$ divides $\Omega$ into two sub-domains $\Omega_1$ and $\Omega_2$ such that $\Omega_1 \cap \Omega_2 = \emptyset$, $\partial\Omega_1 \cap \partial\Omega_2 = \Sigma$ and $\partial\Omega = (\partial\Omega_1 \cup \partial\Omega_2) \setminus \Sigma$, where $\partial\Omega_i$ and $\partial\Omega$ stand for the boundaries of $\Omega_i$ and $\Omega$, respectively. Moreover, $\partial\Omega$ is partitioned as $\partial\Omega = \Gamma_D \cup \Gamma_N$ with $\Gamma_D$ (respectively $\Gamma_N$) being the part of the boundary on which Dirichlet (respectively Neumann) conditions are enforced. The symbol $n_i$ denotes the unit outward normal to $\Omega_i$. Note that $\Omega_1$ and $\Omega_2$ are arbitrarily chosen so that $n_1 = n_{\Sigma}$ and $n_2 = -n_{\Sigma}$ on $\Sigma$. We also denote by $n$ the outward unit normal vector to $\partial\Omega$.

In order to ease the presentation, we shall make use of the following assumption.

**Assumption 2.1.** The Neumann boundary $\Gamma_N$ is away from the interface $\Sigma$; i.e. $\Gamma_N \cap \Sigma = \emptyset$. Moreover, both sub-domains have a Neumann boundary; i.e. $\forall i \in \{1, 2\}$, $\partial\Omega_i \cap \Gamma_N \neq \emptyset$.

**Remark 2.2.** This assumption is chosen to deal with a configuration of two heart cavities separated by a closed valve. The case of one or both sub-domains not having a Neumann boundary can be easily treated by changing the spaces representing the pressure at continuous and discrete levels in order to take into account a zero-mean constraint.

An illustration of admissible and non-admissible configurations is given in Figure 1. The following notation

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Left: admissible geometric configuration. Right: non-admissible configuration (see Assumption 2.1 and Remark 2.2).</figcaption></figure>

will be extensively used in the paper. Let $\phi$ be a scalar or tensorial field defined in $\Omega$ and possibly discontinuous across $\Sigma$, we define its sided restrictions, denoted by $\phi_1$ and $\phi_2$, as

$$ \phi_1(\mathbf{x}) \stackrel{\text{def}}{=} \lim_{\delta \to 0^-} \phi(\mathbf{x} + \delta \mathbf{n}_1), \quad \phi_2(\mathbf{x}) \stackrel{\text{def}}{=} \lim_{\delta \to 0^-} \phi(\mathbf{x} + \delta \mathbf{n}_2), \quad \forall \mathbf{x} \in \Sigma, $$

and the following jump operators across the interface $\Sigma$:

$$ [[\phi]] \stackrel{\text{def}}{=} \phi_1 - \phi_2, \quad [[\phi n]] \stackrel{\text{def}}{=} \phi_1 n_1 + \phi_2 n_2 = [[\phi]] n_{\Sigma}. $$
For any subset $\omega \subset \Omega$, we will use the standard Sobolev spaces $H^m(\omega)$ ($m \ge 0$), equipped with the standard norm $\|\cdot\|_{H^m(\omega)}$, while $(\cdot, \cdot)_\omega$ denotes the usual $L^2(\omega)$-scalar product. We also denote by $|\omega|$ the usual Lebesgue measure of $\omega$. For $S \subset \partial\Omega \cup \Sigma$, the closed subspace $H_S^1$ represents the space of $H^1$ functions with zero value on $S$.

Let $f \in [L^2(\Omega)]^d$, $g \in [L^2(\Gamma_N)]^d$ and $u_\Sigma \in [H_0^1(\Sigma)]^d$ be given fields. We consider the following Stokes problem in $\Omega$ with an immersed Dirichlet condition: Find $u: \Omega \to \mathbb{R}^d$ and $p: \Omega \to \mathbb{R}$, fluid velocity and pressure respectively, such that:

$$ egin{cases} -\nabla \cdot \sigma(\mathbf{u}, p) = \mathbf{f} & \text{in } \Omega_i, \\ \nabla \cdot \mathbf{u} = 0 & \text{in } \Omega_i, \\ \mathbf{u} = \mathbf{0} & \text{on } \Gamma_D, \\ \sigma(\mathbf{u}, p)\mathbf{n} = \mathbf{g} & \text{on } \Gamma_N, \\ \mathbf{u} = \mathbf{u}_\Sigma & \text{on } \Sigma, \end{cases} \quad (1) $$

for $i=1,2$. Here, the symbol $\sigma(\mathbf{u}, p)$ stands for the Cauchy stress tensor of the fluid, given by

$$ \sigma(\mathbf{u}, p) \stackrel{\text{def}}{=} 2\mu\varepsilon(\mathbf{u}) - p\mathbf{I}, \quad \varepsilon(\mathbf{u}) \stackrel{\text{def}}{=} \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T), $$

where $\mu > 0$ denotes the fluid dynamic viscosity.

Note that system (1) corresponds to two independent Stokes problems with a shared Dirichlet condition on the interface $\Sigma$. This kind of problems arises, for instance, in fluid-structure interaction problems with immersed thin-walled solids (see, e.g., [13]). The goal of the present paper is to propose and analyze a fictitious domain approximation of (1) in which the underlying computational mesh is not fitted to $\Sigma$. To this purpose, we first introduce a global variational formulation of (1) in which Lagrange multipliers are used to enforce the Dirichlet constraint on $\Sigma$.

We consider the following velocity and pressure functional spaces

$$ V \stackrel{\text{def}}{=} [H_{\Gamma_D}^1(\Omega)]^d, \quad Q \stackrel{\text{def}}{=} L^2(\Omega), $$

respectively. We make use of the standard Stokes bilinear form $a: (V \times Q) \times (V \times Q) \to \mathbb{R}$ and the linear form $\ell: V \to \mathbb{R}$, given by

$$ a((\mathbf{u}, p), (\mathbf{v}, q)) \stackrel{\text{def}}{=} 2\mu(\varepsilon(\mathbf{u}), \varepsilon(\mathbf{v}))_{\Omega} - (p, \nabla \cdot \mathbf{v})_{\Omega} + (q, \nabla \cdot \mathbf{u})_{\Omega}, \quad \ell(\mathbf{v}) \stackrel{\text{def}}{=} (\mathbf{f}, \mathbf{v})_{\Omega} + (\mathbf{g}, \mathbf{v})_{\Gamma_N}. $$

For the weak treatment of the Dirichlet interface condition on $\Sigma$, we consider the Lagrange multiplier $\mathbf{\Lambda} \stackrel{\text{def}}{=} ([H_{00}^{\frac{1}{2}}(\Sigma)]^d)'$, defined as the dual space of $[H_{00}^{\frac{1}{2}}(\Sigma)]^d$, and the bilinear form $b: \mathbf{\Lambda} \times [H_{00}^{\frac{1}{2}}(\Sigma)]^d \to \mathbb{R}$, given by

$$ b(\xi, z) \stackrel{\text{def}}{=} \langle \xi, z \rangle, $$

where $\langle \cdot, \cdot \rangle$ denotes the duality pairing between $\mathbf{\Lambda}$ and $[H_{00}^{\frac{1}{2}}(\Sigma)]^d$.

The weak formulation of problem (1) then reads as follows: Find $(\mathbf{u}, p, \mathbf{\lambda}) \in V \times Q \times \mathbf{\Lambda}$ such that

$$ \mathcal{A}((\mathbf{u}, p, \mathbf{\lambda}), (\mathbf{v}, q, \xi)) = \mathcal{F}(\mathbf{v}, \xi), \quad (2) $$

for all $(\mathbf{v}, q, \xi) \in V \times Q \times \mathbf{\Lambda}$, where the operators $\mathcal{A}$ and $\mathcal{F}$ are defined as follows:

$$ \mathcal{A}((\mathbf{u}, p, \mathbf{\lambda}), (\mathbf{v}, q, \xi)) \stackrel{\text{def}}{=} a((\mathbf{u}, p), (\mathbf{v}, q)) - b(\mathbf{\lambda}, \mathbf{v}) + b(\xi, \mathbf{u}), \quad \mathcal{F}(\mathbf{v}, \xi) \stackrel{\text{def}}{=} \ell(\mathbf{v}) + b(\xi, \mathbf{u}_\Sigma). \quad (3) $$

Using a standard argument of integration by parts in $\Omega_i$, it can be shown that problem (1) is equivalent to problem (2) with

$$ \mathbf{\lambda} = [[\sigma(\mathbf{u}, p)\mathbf{n}]]]_{\Sigma}. \quad (4) $$
In other words, the unknown Lagrange multiplier represents the jump of the fluid stress across the interface $\Sigma$ resulting from the enforcement of the Dirichlet condition on $\Sigma$. Moreover, it is classical to show that problem (2) is well-posed.

As regards the regularity of the solution $(u, p)$, it should be noted that the Dirichlet constraint in (1) does not ensure continuity of the pressure or of the velocity gradient. Concerning the functional framework of problem (1), the potential presence of discontinuities on the interface $\Sigma$ for the solution (resp. its gradient) reduces the regularity of $p$ (resp. $u$). More precisely, instead of the standard spaces $H^1(\Omega)$ and $[H^2(\Omega)]^d$ for the pressure and velocity, we use the spaces $H^{1/2-\epsilon}(\Omega)$ and $[H^{3/2-\epsilon}(\Omega)]^d$, for some $\epsilon \in (0, \frac{1}{2})$. Furthermore, the regularity of $p$ (resp. $u$) in each subdomain $\Omega_i$, with $i = 1, 2$, is limited to $H^{\frac{1}{2}+\epsilon}(\Omega_i)$ (resp. $[H^{\frac{3}{2}+\epsilon}(\Omega_i)]^d$) due to the potential presence of reentrant corners on the interface and corners between $\Gamma_D$ and $\Gamma_N$ (see Figure 1). For further details on these aspects the reader is referred to [34, 45]. Accordingly, we introduce the following spaces:

$$ 
\begin{aligned}
V_{\epsilon} &= \{v \in [H_{\Gamma_D}^{\frac{3}{2}-\epsilon}(\Omega)]^d | v|_{\Omega_i} \in [H^{\frac{3}{2}+\epsilon}(\Omega_i)]^d, i = 1, 2\}, \\
Q_{\epsilon} &= \{q \in H^{\frac{1}{2}-\epsilon}(\Omega) | q|_{\Omega_i} \in H^{\frac{1}{2}+\epsilon}(\Omega_i), i = 1, 2\}.
\end{aligned} 
\quad (5) 
$$

Note that the traces of elements of these spaces are well defined on each sub-domain boundary $\partial\Omega_i$.

### 3. DISCRETE PROBLEM

Let us introduce $\{\mathcal{T}_h^\Omega\}_{0<h<1}$ and $\{\mathcal{T}_h^\Sigma\}_{0<h<1}$ two families of quasi uniform simplicial meshes (see, e.g., [24]) of $\Omega$ and $\Sigma$, respectively. The symbol $h$ (respectively $\mathcal{H}$) denotes the mesh parameter for $\Omega$ (respectively $\Sigma$). We assume that $\mathcal{T}_h^\Omega$ is fitted to $\partial\Omega$ but, in general, not to $\Sigma$. We consider the two following standard spaces of continuous piecewise affine functions:

$$ X_h^\Omega \stackrel{\text{def}}{=} \{\chi_h \in C^0(\bar{\Omega}) | \chi_h|_K \in \mathbb{P}_1(K), \forall K \in \mathcal{T}_h^\Omega\}, \quad X_h^\Sigma \stackrel{\text{def}}{=} \{\chi_\mathcal{H} \in C^0(\bar{\Sigma}) | \chi_\mathcal{H}|_K \in \mathbb{P}_1(K), \forall K \in \mathcal{T}_h^\Sigma\}. $$

Finally, we introduce the discrete space $\tilde{Q}_h \stackrel{\text{def}}{=} X_h^\Omega \subset Q$ and the discrete spaces $V_h, Q_h, \Lambda_\mathcal{H}$ for the approximation of the velocity, pressure and Lagrange multiplier as follows:

$$ V_h \stackrel{\text{def}}{=} [X_h^\Omega]^d \cap V, \quad Q_h \stackrel{\text{def}}{=} \tilde{Q}_h \oplus \langle 1_{\Omega_1} \rangle \subset Q, \quad \Lambda_\mathcal{H} \stackrel{\text{def}}{=} [X_\mathcal{H}^\Sigma]^d \subset \Lambda, \quad (6) $$

where the symbol $1_{\Omega_1}$ denotes the characteristic function of $\Omega_1$ and $\langle 1_{\Omega_1} \rangle$ the vector space spanned by $1_{\Omega_1}$. We recall that the choice of $\Omega_1$ and $\Omega_2$ is arbitrary, so that the characteristic function could have been taken over $\Omega_2$ instead of over $\Omega_1$. The Lagrange multiplier space $\Lambda_\mathcal{H}$ could alternatively be made of discontinuous functions (as in, e.g., [7, 31]).

As mentioned in the introduction, it is well known that the approximation of the pressure in (2) with continuous finite elements guarantees global mass conservation in the computational fluid domain $\Omega$, but does not prevent spurious leakage across the interface $\Sigma$ (see, e.g., Section 5.2, and particularly Figure 9a). The choice of the non-standard pressure space $Q_h$ can be considered as the simplest enrichment of the usual discrete continuous space $\tilde{Q}_h$ in order to overcome the artificial interfacial mass loss induced by the continuous nature of the pressure approximation. Further advanced enrichment techniques have been suggested. Notably, the extended finite element method (XFEM) (see, e.g., [35]) involves the addition of supplementary basis functions for each cut cell. Notice that, in the present method, by using a single supplementary basis function for the pressure space, the size of the linear system is independent of the interface position.

Owing to the definition of the discrete pressure space (6), every $q_h \in Q_h$ can be decomposed into the sum of two contributions, namely,

$$ q_h = \tilde{q}_h + \hat{q}_h, \quad (7) $$
where $\tilde{q}_h \in \tilde{Q}_h$ is a continuous piecewise affine function over the whole computational domain $\Omega$ and $\hat{q}_h \in \langle 1_{\Omega_1} \rangle$ is a Heaviside function that allows strong discontinuities across the interface $\Sigma$. Moreover, the following relations hold

$$ \hat{q}_h = [[q_h]]1_{\Omega_1}, \quad \tilde{q}_h = q_h - [[q_h]]1_{\Omega_1}, \qquad (8) $$

where $[[q_h]]$ is the constant jump of $q_h$ across $\Sigma$.

In order to overcome the instability of the saddle-point problem due to the choice of the Lagrange multiplier finite element space $\Lambda_H$ in (6) and to avoid any constraint on the ratio between the fluid mesh size $h$ and the solid mesh size $\mathcal{H}$ (as can, for instance, be seen in [12]), we introduce the following stabilization term (in the spirit of the Barbosa-Hughes stabilization [6]):

$$ s_h^{\text{BH}}((p_h, \lambda_H), (q_h, \xi_H)) \stackrel{\text{def}}{=} \frac{h}{\gamma_\lambda \mu} (\lambda_H + [[p_h]]\mathbf{n}_\Sigma, \xi_H + \theta[[q_h]]\mathbf{n}_\Sigma)_\Sigma, \qquad (9) $$

with $\gamma_\lambda > 0$ a user-defined dimensionless parameter and $\theta \in \{0, 1\}$. The choice $\theta = 1$ yields a symmetric stabilization term, while for $\theta = 0$ we get a non-symmetric method. Both variants will be considered in the sequel and, as discussed below in Remark 3.2, they provide different accuracy in terms of interfacial mass conservation. Furthermore, the following Brezzi-Pitkäranta stabilization [14] is considered to overcome the lack of inf-sup compatibility between the velocity and pressure spaces $V_h$ and $Q_h$:

$$ s_h^{\text{BP}}(\tilde{p}_h, \tilde{q}_h) \stackrel{\text{def}}{=} \frac{\gamma_p h^2}{\mu} (\nabla \tilde{p}_h, \nabla \tilde{q}_h)_\Omega, $$

with $\gamma_p > 0$ a user-defined dimensionless parameter. Note that the pressure stabilization operator $s_h^{\text{BP}}$ only acts on the continuous part, $\tilde{p}_h$, of the discrete pressure, $p_h \in Q_h$.

The proposed discrete approximation of (2) reads as follows: Find $(u_h, p_h, \lambda_H) \in V_h \times Q_h \times \Lambda_H$ such that

$$ \mathcal{A}_h((u_h, p_h, \lambda_H), (v_h, q_h, \xi_H)) = \mathcal{F}(v_h, \xi_H) \qquad (10) $$

for all $(v_h, q_h, \xi_H) \in V_h \times Q_h \times \Lambda_H$, with the notation

$$ \mathcal{A}_h((u_h, p_h, \lambda_H), (v_h, q_h, \xi_H)) \stackrel{\text{def}}{=} \mathcal{A}((u_h, p_h, \lambda_H), (v_h, q_h, \xi_H)) + s_h^{\text{BP}}(\tilde{p}_h, \tilde{q}_h) + s_h^{\text{BH}}((p_h, \lambda_H), (q_h, \xi_H)). \quad (11) $$

**Remark 3.1.** It is worth noticing that the stabilization operator (9) is only weakly consistent since, according to (4), the Lagrange multiplier $\lambda$ is not necessarily equal to $-[[p]]\mathbf{n}_\Sigma$. In the spirit of [6], a strongly consistent method would be

$$ \frac{h}{\gamma_\lambda \mu} (\lambda_H - [[\boldsymbol{\sigma}(u_h, p_h)\mathbf{n}]], \xi_H - \theta[[\boldsymbol{\sigma}(v_h, q_h)\mathbf{n}]])_\Sigma. $$

The reason of excluding the normal component of the strain tensor $[[\boldsymbol{\varepsilon}(v_h)]]\mathbf{n}_\Sigma$ is related to the fact that for piecewise affine velocity approximations, this contribution vanishes in any cut element whose boundary is not conformal to the interface. If, instead of $Q_h$, we consider the standard piecewise affine continuous space $\tilde{Q}_h$ for the pressure, the stabilization operator (9) reduces to the one investigated in [7] with piecewise constant Lagrange multipliers. The resulting unconstrained formulation will be considered for comparison purposes in the numerical experiments of Section 5.

**Remark 3.2.** Due to the choice of the finite element space $Q_h$ in (6), we have

$$ (\nabla \cdot v_h, q_h)_\Omega = (\nabla \cdot v_h, \tilde{q}_h)_\Omega + [[q_h]](v_h \cdot \mathbf{n}_1, 1)_{\partial\Omega_1} = (\nabla \cdot v_h, \tilde{q}_h)_\Omega + [[q_h]] \int_{\partial\Omega_1} v_h \cdot \mathbf{n}_1 $$

for all $(v_h, q_h) \in V_h \times Q_h$. Therefore, owing to (7)-(8), the system matrix associated to (10) can be assembled as a combination of two mass conservation constraints and two stabilization terms. The first mass conservation
constraint corresponds to the traditional divergence free constraint imposed via continuous piecewise affine pressures, namely,

$$ (\nabla \cdot \mathbf{u}_h, \tilde{q}_h)_\Omega + s_h^{\text{BP}} (\tilde{p}_h, \tilde{q}_h) = 0 \quad \forall \tilde{q}_h \in \tilde{Q}_h, \qquad (12) $$

which yields global mass conservation in $\Omega$, while the second enhances interfacial mass conservation by enforcing

$$ \int_{\partial \Omega_1} \mathbf{u}_h \cdot \mathbf{n}_1 + \theta \frac{h}{\gamma \lambda \mu} \int_{\Sigma} (\lambda_H + [\![p_h]\!] \mathbf{n}_\Sigma) \cdot \mathbf{n}_\Sigma = 0. \qquad (13) $$

For $\theta = 0$, this relation guarantees global mass conservation in $\Omega_1$ (and also in $\Omega_2$, thanks to (12)). On the contrary, for $\theta = 1$, interfacial mass conservation is only generally enforced in a weakly consistent fashion. The constraint (13) is imposed via the constant Lagrange multiplier $[\![p_h]\!]$, which also represents the value of the discrete pressure jump across the interface. To the best of our knowledge, the introduction of the additional constraint (13) with $\theta = 0$ for enhancing interfacial mass conservation in a fictitious domain context was reported for the first time in [38], for a different method and without error analysis. Similar ideas have also been advocated in [27, 46] in a different physical framework.

## 4. NUMERICAL ANALYSIS

This section is devoted to the a priori stability and error analysis of the numerical method provided by (10) and is organized as follows. We first introduce several interpolation operators, based on quasi-interpolation. Then, we prove inf-sup stability of the bilinear form $A_h$. We then bound the consistency error of the scheme. Finally, we combine all these results to establish the desired a priori error bounds fulfilled by the numerical solution.

In what follows, we write $A \lesssim B$ (resp. $A \gtrsim B$) to abbreviate the inequality $A \le CB$ (resp. $CA \ge B$) for positive real numbers $A$ and $B$, where the constant $C > 0$ does not depend on $h$, $\mathcal{H}$, $\mu$, neither on the way the interface mesh $\mathcal{T}_h^\Omega$ cuts the fluid mesh $\mathcal{T}_h^\Omega$ nor on the solutions of the continuous and discrete problems. We also use the notation $\mathbf{H}^s(\omega) \stackrel{\text{def}}{=} [\![H^s(\omega)]\!]^d$ for all $s > 0$ and all $\omega \subset \Omega$.

We now introduce some geometric definitions that are widely used in the current section. Let $K$ be an element of $\mathcal{T}_h^\Omega$ and $\mathbf{x}_{K,i}$ its $i$-th node. We denote by $\mathcal{T}_h^c$ the set of elements of $\mathcal{T}_h^\Omega$ cut by the interface $\Sigma$, by $\mathcal{N}_h$ the set of nodes, by $\mathcal{N}_h^c$ the set of nodes belonging to at least one cut element and by $\mathcal{T}_h^e$ the set of elements of $\mathcal{T}_h^\Omega$ sharing at least one node with a cut element. See Figure 2. Let us now introduce the following assumptions.

<figure><img src="image_3.png" /><figcaption>FIGURE 2. In gray the fluid domain triangulation $\mathcal{T}_h^\Omega$ and in black the interface domain triangulation $\mathcal{T}_h^c$. The elements of $\mathcal{T}_h^c$ are marked in $\square$, the elements of $\mathcal{T}_h^e$ are marked in $\square$ and the nodes in $\mathcal{N}_h^c$ are marked by $\bullet$.</figcaption></figure>
**Assumption 4.1.** We assume that the interface is well resolved by $T_h^\Omega$, i.e. $h|\Sigma \cap K| \lesssim |K| \quad \forall K \in T_h^\Omega$.

**Assumption 4.2.** We assume that every cut cell $K \in T_h^c$ shares at least one vertex with an uncut cell $\tilde{K}_1 \subset \Omega_1$ and at least another one with an uncut cell $\tilde{K}_2 \subset \Omega_2$, with $\tilde{K}_1, \tilde{K}_2 \in T_h^\Omega \setminus T_h^c$.

These assumptions are fulfilled when the mesh is fine enough.

## 4.1. Approximability

In this section, quasi-interpolation operators are considered to approximate non-smooth functions (see, e.g., [24, Chapter 22]). We recall that the exact solution of (2) has a low regularity due to the presence of weak and strong discontinuities across the interface (velocity and pressure, respectively). This is why we propose the following adapted approximation operators, based on quasi-interpolation, that are defined even for low-regularity functions.

(1) $\Pi_h : V \to V_h$ is a given quasi-interpolation operator onto $V_h$;

(2) $\Pi_h : Q_\epsilon \to Q_h$, is defined as

$$ \Pi_h q \stackrel{\text{def}}{=} \tilde{\Pi}_h q + \hat{\Pi}_h q \quad \forall q \in Q_\epsilon, $$

with

$$ \hat{\Pi}_h q \stackrel{\text{def}}{=} \left( \frac{1}{|\Sigma|} \int_{\Sigma} \|q\| \right) 1_{\Omega_1}, \quad (14) $$

the mean value of the pressure jump across the interface $\Sigma$, and

$$ \tilde{\Pi}_h q \stackrel{\text{def}}{=} \bar{\Pi}_h (q - \hat{\Pi}_h q), $$

with $\bar{\Pi}_h : Q_\epsilon \to \tilde{Q}_h$ a given quasi-interpolation operator onto $\tilde{Q}_h$.

(3) $\Pi_H : [L^2(\Sigma)]^d \to \Lambda_H$, a given quasi-interpolation operator onto $\Lambda_H$.

Owing to the properties of the quasi-interpolation operators (see, e.g., [24, Chapter 22]), we have the following result.

**Lemma 4.3.** Let $\epsilon \in (0, \frac{1}{2})$. Assume that $p \in Q_\epsilon$ and $\lambda \in [L^2(\Sigma)]^d$. For the above defined operators the following estimates hold:

$$ \| \llbracket \Pi_h p \rrbracket \|_{L^2(\Sigma)} \le \| \llbracket p \rrbracket \|_{L^2(\Sigma)}, \quad (15) $$

$$ \| \Pi_H \lambda \|_{L^2(\Sigma)} \lesssim \| \lambda \|_{L^2(\Sigma)}, \quad (16) $$

$$ h \| \nabla \tilde{\Pi}_h p \|_{L^2(\Omega)} \lesssim h^{\frac{1}{2}-\epsilon} \| p - \hat{\Pi}_h p \|_{H^{\frac{1}{2}-\epsilon}(\Omega)}, \quad (17) $$

$$ \| \Pi_H (\lambda + \llbracket \hat{\Pi}_h p \rrbracket n_\Sigma) \|_{L^2(\Sigma)} \lesssim \| \lambda + \llbracket \hat{\Pi}_h p \rrbracket n_\Sigma \|_{L^2(\Sigma)}. \quad (18) $$

*Proof.* The estimates (16) and (18) are standard (see, e.g., [24, Corollary 22.8]). We focus on estimate (17). Let $K \in T_h^\Omega$ and $\bar{p}_K$ be the mean value of $p - \hat{\Pi}_h p$ over $S(K)$, the patch of elements around $K$. Using the definition of $\tilde{\Pi}_h$, we have

$$ \begin{aligned} h \| \nabla \tilde{\Pi}_h p \|_{L^2(K)} &= h \| \nabla (\tilde{\Pi}_h p - \bar{p}_K) \|_{L^2(K)} \lesssim \| \tilde{\Pi}_h p - \bar{p}_K \|_{L^2(K)} \lesssim \| \bar{\Pi}_h (p - \hat{\Pi}_h p) - \bar{p}_K \|_{L^2(K)} \\ &\lesssim \| \bar{\Pi}_h (p - \hat{\Pi}_h p - \bar{p}_K) \|_{L^2(K)} \lesssim \| p - \hat{\Pi}_h p - \bar{p}_K \|_{L^2(S(K))} \lesssim h^{\frac{1}{2}-\epsilon} \| p - \hat{\Pi}_h p \|_{H^{\frac{1}{2}-\epsilon}(S(K))}, \end{aligned} $$

where we used an inverse inequality, the local stability of $\bar{\Pi}_h$ (see, e.g., [24, Corollary 22.8]) and approximation properties of the mean value. Summing the square of this relation over all $K \in T_h^\Omega$ and using the quasi-uniformity
of the mesh gives (17). For (15), we use the fact that $\llbracket \Pi_h p \rrbracket$ is a constant and that $\llbracket \Pi_h p \rrbracket_{1_{\Omega_1}} = \hat{\Pi}_h p$, so that we have

$$ \|\llbracket \Pi_h p \rrbracket\|_{L^2(\Sigma)}^2 = |\Sigma| \|\llbracket \Pi_h p \rrbracket\|^2 = |\Sigma|^{-1} \left( \int_{\Sigma} \|p\| \right)^2 \le \|\llbracket p \rrbracket\|_{L^2(\Sigma)}^2, $$

where for the last inequality we used Cauchy-Schwarz on the squared term. This completes the proof. □

For $u \in V$, $p \in Q_{\varepsilon}$ and $\lambda \in [L^2(\Sigma)]^d$, we introduce the following approximation errors

$$ e_u \stackrel{\text{def}}{=} u - \Pi_h u, \quad e_p \stackrel{\text{def}}{=} p - \Pi_h p, \quad e_{\lambda} \stackrel{\text{def}}{=} \lambda - \Pi_{\mathcal{H}} \lambda. \qquad (19) $$

We have the following approximation result.

**Lemma 4.4.** Let $\varepsilon \in (0, \frac{1}{2})$. Assume that $u \in V \cap H^{\frac{3}{2}-\varepsilon}(\Omega)$, $p \in Q_{\varepsilon}$ and $\lambda \in H^{\frac{1}{2}-\varepsilon}(\Sigma)$. The following estimates hold:

$$ \|\nabla \cdot e_u\|_{L^2(\Omega)} \lesssim \|\varepsilon(e_u)\|_{L^2(\Omega)} \lesssim h^{\frac{1}{2}-\varepsilon} \|u\|_{H^{\frac{3}{2}-\varepsilon}(\Omega)}, \qquad (20) $$

$$ \|e_u\|_{L^2(\Sigma)} \lesssim h^{1-\varepsilon} \|u\|_{H^{\frac{3}{2}-\varepsilon}(\Omega)}, \qquad (21) $$

$$ \|e_{\lambda}\|_{L^2(\Sigma)} \lesssim \mathcal{H}^{\frac{1}{2}-\varepsilon} \|\lambda\|_{H^{\frac{1}{2}-\varepsilon}(\Sigma)}, \qquad (22) $$

$$ \|e_p\|_{L^2(\Omega)} \lesssim h^{\frac{1}{2}-\varepsilon} \|p - \hat{\Pi}_h p\|_{H^{\frac{1}{2}-\varepsilon}(\Omega)}, \qquad (23) $$

$$ \|\lambda + \llbracket \hat{\Pi}_h p \rrbracket n_{\Sigma} - \Pi_{\mathcal{H}} (\lambda + \llbracket \hat{\Pi}_h p \rrbracket n_{\Sigma})\|_{L^2(\Sigma)} \lesssim \mathcal{H}^{\frac{1}{2}-\varepsilon} \|\lambda + \llbracket \hat{\Pi}_h p \rrbracket n_{\Sigma}\|_{H^{\frac{1}{2}-\varepsilon}(\Sigma)}. \qquad (24) $$

*Proof.* Estimates (20)–(22) and (24) can be easily proven using the approximation properties of quasi-interpolators (see, e.g., [24, Theorem 22.6]). For (23), we use the definition of $\tilde{\Pi}_h$ to get

$$ \begin{aligned} \|e_p\|_{L^2(\Omega)} &= \|p - \Pi_h p\|_{L^2(\Omega)} = \|p - \hat{\Pi}_h p - \tilde{\Pi}_h p\|_{L^2(\Omega)} \\ &\lesssim h^{\frac{1}{2}-\varepsilon} \|p - \hat{\Pi}_h p\|_{H^{\frac{1}{2}-\varepsilon}(\Omega)}, \end{aligned} $$

where we used the approximation properties of the quasi-interpolation operator $\tilde{\Pi}_h$, so that the proof is complete. □

<figure><img src="image_4.png" /></figure>

## 4.2. Stability

In this section we establish inf-sup stability of the bilinear form $A_h$ (see Lemma 4.8 below). The associated proof is based on two main steps: first, the proof of an inf-sup condition for the bilinear form $(p_h, \nabla \cdot v_h)_{\Omega}$, see Proposition 4.6; second, the extension of inf-sup stability to the full bilinear form $A_h$. Even if these steps are standard in the context of fitted meshes, some intermediate lemmas have to be proven in the present context of unfitted meshes. In Proposition 4.5, we introduce an *ad hoc* interpolation operator vanishing on the elements intersected by the interface $\Sigma$. This interpolation is used to establish the inf-sup stability of the bilinear form $(p_h, \nabla \cdot v_h)_{\Omega}$, where the pressure is decomposed into a constant and a zero-mean valued function. More particularly, this interpolation operator is used to construct an adequate test function for the zero-mean part of the pressure. In Proposition 4.7, we extend the well-known discrete trace inequality on fitted meshes [24, Lemma 12.8] to our unfitted setting. Similar propositions are also employed in [28, 36].

### 4.2.1. Preliminary results

We start by constructing an adequate interpolation operator that will be used to establish the inf-sup condition.

**Proposition 4.5.** There exists a linear operator $\Pi_h^0 : V \to V_h$ such that for every $v \in V$, we have:
• $\Pi_h^0 v = 0$ on $\Sigma$;

• $\Pi_h^0 v = 0$ on $\partial\Omega$ if $v = 0$ on $\partial\Omega$.

Moreover, for each $v \in V$ such that $v = 0$ on $\Sigma$, it holds

$$ \|v - \Pi_h^0 v\|_{L^2(\Omega)} \lesssim h |v|_{H^1(\Omega)}, \qquad (25) $$

$$ \|\Pi_h^0 v\|_{H^1(\Omega)} \lesssim \|v\|_{H^1(\Omega)}. \qquad (26) $$

*Proof.* We start by considering the Scott-Zhang interpolation operator $\Pi_h^Z : V \to V_h$ (see, e.g., [49]). For all $v \in V$, $\Pi_h^0 v$ is defined as the sole element of $V_h$ satisfying

$$ \Pi_h^0 v(x) = 0 \quad \forall x \in \mathcal{N}_h^c, \qquad (27) $$

$$ \Pi_h^0 v(x) = \Pi_h^Z v(x) \quad \forall x \in \mathcal{N}_h \setminus \mathcal{N}_h^c. \qquad (28) $$

Note that, by construction, we have $\Pi_h^0 v = 0$ on $\Sigma$ and

$$ \Pi_h^Z v|_K = \Pi_h^0 v|_K \quad \forall K \in \mathcal{T}_h^\Omega \setminus \mathcal{T}_h^e. \qquad (29) $$

Furthermore, owing to the preservation of homogeneous Dirichlet boundary conditions by the Scott-Zhang interpolation (see, e.g., [49, Theorem 2.1]), we have that $\Pi_h^0 v = 0$ on $\partial\Omega$ for all $v \in V$ with $v = 0$ on $\partial\Omega$.

Let us now consider $v \in V$ such that $v = 0$ on $\Sigma$. We start with the proof of (25). By using the triangle inequality and the approximation properties of the Scott-Zhang operator (see, e.g., [49, Theorem 4.1]), it follows that

$$ \begin{aligned} \|v - \Pi_h^0 v\|_{L^2(\Omega)} &\le \|v - \Pi_h^Z v\|_{L^2(\Omega)} + \|\Pi_h^Z v - \Pi_h^0 v\|_{L^2(\Omega)} \\ &\lesssim h |v|_{H^1(\Omega)} + \|\Pi_h^Z v - \Pi_h^0 v\|_{L^2(\Omega)}. \end{aligned} \qquad (30) $$

It only remains to control the last term in the previous bound. Owing to (29), we can focus only on the mesh elements in $\mathcal{T}_h^e$. Let $K \in \mathcal{T}_h^e$, from [24, Proposition 12.5], (27) and (28), we have

$$ \begin{aligned} \|\Pi_h^Z v - \Pi_h^0 v\|_{L^2(K)} &\lesssim h^{\frac{d}{2}} \max_i \|\Pi_h^Z v(x_{K,i}) - \Pi_h^0 v(x_{K,i})\| \\ &\lesssim h^{\frac{d}{2}} \max_{i: x_{K,i} \in \mathcal{N}_h^c} \|\Pi_h^Z v(x_{K,i}) - \Pi_h^0 v(x_{K,i})\| \\ &\lesssim h^{\frac{d}{2}} \max_{i: x_{K,i} \in \mathcal{N}_h^c} \|\Pi_h^Z v(x_{K,i})\| \lesssim h^{\frac{d}{2}} \max_i \|\Pi_h^Z v(x_{K,i})\| \\ &\lesssim \|\Pi_h^Z v\|_{L^2(K)}. \end{aligned} \qquad (31) $$

In addition, from the stability of the Scott-Zhang interpolation operator (see, e.g., [49, Theorem 3.1]), we have

$$ \|\Pi_h^Z v\|_{L^2(K)} \le C \left( \|v\|_{L^2(S(K))} + h |v|_{H^1(S(K))} \right), $$

where $S(K)$ denotes the patch of elements around $K$. Finally, since $v$ vanishes on $\Sigma \cap S(K) \ne \emptyset$, the Poincaré's inequality $\|v\|_{L^2(S(K))} \lesssim h |v|_{H^1(S(K))}$ holds (see [24, Lemma 3.30]), so that from (31) we get

$$ \|\Pi_h^Z v - \Pi_h^0 v\|_{L^2(K)} \lesssim h |v|_{H^1(S(K))} \quad \forall K \in \mathcal{T}_h^e. \qquad (32) $$

Since from (29) the estimate (32) holds for every $K \in \mathcal{T}_h^\Omega$, we can sum (32) over all $K \in \mathcal{T}_h^\Omega$ and, by using the quasi-uniformity of the mesh, we finally get

$$ \|\Pi_h^Z v - \Pi_h^0 v\|_{L^2(\Omega)} \lesssim h |v|_{H^1(\Omega)}. \qquad (33) $$
The estimate (25) simply follows by inserting (33) into (30).

In order to prove (26), we first notice that owing to (29), it suffices to consider the elements $K$ in $\mathcal{T}_h^e$. By combining (31) with the triangle inequality, we have

$$ \|\Pi_h^0 v\|_{L^2(K)} \lesssim \|\Pi_h^Z v\|_{L^2(K)} \quad \forall K \in \mathcal{T}_h^e. \qquad (34) $$

Similarly, using inverse inequalities and (32), for the gradient we have

$$ \begin{aligned} \|\nabla \Pi_h^0 v\|_{L^2(K)} & \lesssim h^{-1} \|\Pi_h^Z v - \Pi_h^0 v\|_{L^2(K)} + \|\nabla \Pi_h^Z v\|_{L^2(K)} \\ & \lesssim |v|_{H^1(S(K))} + \|\nabla \Pi_h^Z v\|_{L^2(K)} \quad \forall K \in \mathcal{T}_h^e. \end{aligned} \qquad (35) $$

Owing to (29), the relations (34) and (35) hold for every $K \in \mathcal{T}_h^\Omega$, so that by summing over $K$ and using the $H^1$-stability of the Scott-Zhang interpolator (see, e.g., [49, Corollary 4.1]), we finally get

$$ \|\Pi_h^0 v\|_{H^1(\Omega)}^2 \lesssim \|\Pi_h^Z v\|_{H^1(\Omega)}^2 + |v|_{H^1(\Omega)}^2 \lesssim \|v\|_{H^1(\Omega)}^2, $$

which completes the proof.

<figure><img src="image_5.png" /></figure>

In the sequel, for $i \in \{1, 2\}$ and for any element $q_h$ of $Q_h$, we denote by $q_{h,i} \stackrel{\text{def}}{=} q_h|_{\Omega_i}$ the restriction of $q_h$ to $\Omega_i$. In a similar way, we denote by $Q_{h,i} \stackrel{\text{def}}{=} \{q_h|_{\Omega_i} | q_h \in Q_h\}$ the set of such restrictions. Then, for every $q_h \in Q_h$ we have

$$ q_h = \sum_{i=1}^{2} q_{h,i} 1_{\Omega_i}. $$

We can now prove the following intermediate result, which is a first inf-sup condition.

**Proposition 4.6.** *Under Assumption 2.1, there holds*

$$ \sup_{0 \neq v_h \in V_h \cap [H_h^1(\Omega)]^d} \frac{(q_h, \nabla \cdot v_h)_\Omega}{\|v_h\|_{H^1(\Omega)}} + h \|\nabla \tilde{q}_h\|_{L^2(\Omega)} \gtrsim \|q_h\|_{L^2(\Omega)} \qquad (36) $$

for all $q_h \in Q_h$.

*Proof.* For $i \in \{1, 2\}$, we proceed in a classical way by decomposing $q_{h,i}$ as

$$ q_{h,i} = \bar{q}_{h,i} + \tilde{q}_{h,i}, \quad \bar{q}_{h,i} \stackrel{\text{def}}{=} \frac{1}{|\Omega_i|} \int_{\Omega_i} q_{h,i}, \quad \tilde{q}_{h,i} \stackrel{\text{def}}{=} q_{h,i} - \bar{q}_{h,i} \in L_0^2(\Omega_i) \cap Q_{h,i}. \qquad (37) $$

Note that, by restricting the decomposition (7) to each sub-domain $\Omega_i$, we also have $q_{h,i} = \tilde{q}_{h,i} + \hat{q}_{h,i}$ with $\tilde{q}_{h,i} \stackrel{\text{def}}{=} \tilde{q}_h|_{\Omega_i}$ and $\hat{q}_{h,i} \stackrel{\text{def}}{=} \hat{q}_h|_{\Omega_i}$. Furthermore, the fact that $\tilde{q}_{h,i} + \hat{q}_{h,i} = q_{h,i} = \tilde{q}_{h,i} + \hat{q}_{h,i}$ yields

$$ \nabla \tilde{q}_{h,i} = \nabla \tilde{q}_{h,i}. \qquad (38) $$

We start by treating the zero-mean contributions in (37), namely, $\tilde{q}_h \stackrel{\text{def}}{=} \sum_i \tilde{q}_{h,i} 1_{\Omega_i}$. Since $\tilde{q}_{h,i} \in L_0^2(\Omega_i)$, there exists $0 \neq \tilde{v}_i \in H_0^1(\Omega_i)$ such that (see, e.g., [32])

$$ (\tilde{q}_{h,i}, \nabla \cdot \tilde{v}_i)_{\Omega_i} = \|\tilde{q}_{h,i}\|_{L^2(\Omega_i)}^2, \quad \|\tilde{v}_i\|_{H^1(\Omega_i)} \lesssim \|\tilde{q}_{h,i}\|_{L^2(\Omega_i)}. \qquad (39) $$

By setting $\tilde{v} \stackrel{\text{def}}{=} \sum_i \tilde{v}_i 1_{\Omega_i}$, we get

$$ \tilde{v} \in H_{\partial\Omega \cup \Sigma}^1(\Omega), \quad (\tilde{q}_h, \nabla \cdot \tilde{v})_{\Omega} = \|\tilde{q}_h\|_{L^2(\Omega)}^2, \quad \|\tilde{v}\|_{H^1(\Omega)} \lesssim \|\tilde{q}_h\|_{L^2(\Omega)}. \qquad (40) $$
By using the properties of $\Pi_h^0$ (Proposition 4.5), the fact that $(\bar{q}_{h,i}, \nabla \cdot \Pi_h^0 \hat{v}_i)_{\Omega_i} = 0$ and (38), we have

$$
\begin{align*}
\| \hat{q}_{h,i} \|_{L^2(\Omega_i)}^2 &= (\hat{q}_{h,i}, \nabla \cdot \hat{v}_i)_{\Omega_i} = (\hat{q}_{h,i}, \nabla \cdot (\hat{v}_i - \Pi_h^0 \hat{v}_i))_{\Omega_i} + (\hat{q}_{h,i}, \nabla \cdot \Pi_h^0 \hat{v}_i)_{\Omega_i} \\
&= -(\nabla \hat{q}_{h,i}, \hat{v}_i - \Pi_h^0 \hat{v}_i)_{\Omega_i} + (\hat{q}_{h,i}, \nabla \cdot \Pi_h^0 \hat{v}_i)_{\Omega_i} \\
&\lesssim h \| \nabla \hat{q}_{h,i} \|_{L^2(\Omega_i)} \| \hat{v}_i \|_{H^1(\Omega_i)} + (q_{h,i}, \nabla \cdot \Pi_h^0 \hat{v}_i)_{\Omega_i} \\
&= h \| \nabla \tilde{q}_{h,i} \|_{L^2(\Omega_i)} \| \hat{v}_i \|_{H^1(\Omega_i)} + (q_{h,i}, \nabla \cdot \Pi_h^0 \hat{v}_i)_{\Omega_i}.
\end{align*}
$$

Notice that here we have made a little abuse of notation in $\Pi_h^0 \hat{v}_i$, since it has to be understood as $\Pi_h^0$ applied to the extension of $\hat{v}_i$ by zero out of $\Omega_i$.

By summing the above relation over $i \in \{1, 2\}$ and since $\Pi_h^0 \hat{v} = \Pi_h^0 \hat{v}_1 + \Pi_h^0 \hat{v}_2$ (which follows from the linearity of $\Pi_h^0$ and Proposition 4.5), we get

$$
\| \hat{q}_h \|_{L^2(\Omega)}^2 \lesssim h \| \nabla \tilde{q}_h \|_{L^2(\Omega)} \| \hat{v} \|_{H^1(\Omega)} + (q_h, \nabla \cdot \Pi_h^0 \hat{v})_{\Omega},
$$

which, with (40), yields,

$$
\| \hat{q}_h \|_{L^2(\Omega)} \lesssim h \| \nabla \tilde{q}_h \|_{L^2(\Omega)} + \frac{(q_h, \nabla \cdot \Pi_h^0 \hat{v})_{\Omega}}{\| \hat{v} \|_{H^1(\Omega)}}. \quad (41)
$$

Let us now handle the constant part of the pressure $\tilde{q}_{h,i}$. Let $\tilde{\Omega}_i$ be a subdomain of $\Omega_i$ such that $\Gamma_N \cap \bar{\Omega}_i \subset \tilde{\Omega}_i$ and such that the mesh is fitted to $\tilde{\Omega}_i$. From [21, Lemma 3.3] (see also [29, Lemma 3.3]), there exists a constant $C_i > 0$ and $\mathbf{0} \neq \hat{v}_{h,i} \in \{w_h|_{\tilde{\Omega}_i} | w_h \in V_h\} \cap [H^1_{\partial\tilde{\Omega}_i \setminus \Gamma_N}(\tilde{\Omega}_i)]^d$, such that

$$
C_i \| \bar{q}_{h,i} \|_{L^2(\bar{\Omega}_i)} \| \hat{v}_{h,i} \|_{H^1(\bar{\Omega}_i)} \le (\bar{q}_{h,i}, \nabla \cdot \hat{v}_{h,i})_{\bar{\Omega}_i}. \quad (42)
$$

Notice that, since $\tilde{\Omega}_i$ is fitted to the mesh, the constant $C_i$ is independent from the way the unfitted interface cuts the mesh. Moreover, denoting by $\bar{v}_{h,i} \in \{w_h|_{\tilde{\Omega}_i} | w_h \in V_h\} \cap [H^1_{\partial\tilde{\Omega}_i \setminus \Gamma_N}(\tilde{\Omega}_i)]^d$ the extension by zero of $\hat{v}_{h,i}$ to the whole domain $\Omega_i$, the estimate (42) holds with $\tilde{\Omega}_i = \Omega_i$, $\hat{v}_{h,i} = \bar{v}_{h,i}$ and with the constant $C_i$ rescaled by $|\tilde{\Omega}_i|/|\Omega_i|$. According to Assumption 2.1, $\tilde{\Omega}_i$ can be chosen in a way that the ratio $|\tilde{\Omega}_i|/|\Omega_i|$ does not vanish when $h$ tends to zero. We established

$$
\| \bar{q}_{h,i} \|_{L^2(\Omega_i)} \lesssim \frac{(\bar{q}_{h,i}, \nabla \cdot \bar{v}_{h,i})_{\Omega_i}}{\| \bar{v}_{h,i} \|_{H^1(\Omega_i)}}. \quad (43)
$$

We now consider

$$
\bar{z}_h \stackrel{\text{def}}{=} \sum_i \frac{1}{\sqrt{2} \| \bar{v}_{h,i} \|_{H^1(\Omega_i)}} \bar{v}_{h,i} 1_{\Omega_i},
$$

so that $\|\bar{z}_h\|_{H^1(\Omega)} = 1$, and summing (43) for $i \in \{1, 2\}$ we have

$$
\| \bar{q}_h \|_{L^2(\Omega)} \lesssim (\bar{q}_h, \nabla \cdot \bar{z}_h)_{\Omega} = \frac{(\bar{q}_h, \nabla \cdot \bar{z}_h)_{\Omega}}{\| \bar{z}_h \|_{H^1(\Omega)}} \lesssim \frac{(q_h, \nabla \cdot \bar{z}_h)_{\Omega}}{\| \bar{z}_h \|_{H^1(\Omega)}} + \| \hat{q}_h \|_{L^2(\Omega)}. \quad (44)
$$

Setting $w_h = \|\hat{v}\|_{H^1(\Omega)}^{-1} \Pi_h^0 \hat{v} + \delta \|\bar{z}_h\|_{H^1(\Omega)}^{-1} \bar{z}_h$, with $\delta > 0$ sufficiently small, results in $w_h \in V_h \cap [H^1(\Omega)]^d$ and $\|w_h\|_{H^1(\Omega)} \lesssim 1 + \delta$. The relation (36) follows from (41) and (44), which completes the proof. $\square$

**Proposition 4.7.** *Under Assumptions 4.1 and 4.2, we have*

$$
h^{\frac{1}{2}} \| [q_h] \|_{L^2(\Sigma)} \lesssim \|q_h\|_{L^2(\Omega)} + h \| \nabla \tilde{q}_h \|_{L^2(\Omega)} \quad (45)
$$

for all $q_h \in Q_h$.
Proof. Let us introduce the operators $E_h^1: Q_h \to \tilde{Q}_h$ and $E_h^2: Q_h \to \tilde{Q}_h$ defined as

$$E_h^1(q_h) \stackrel{\text{def}}{=} \tilde{q}_h + \llbracket q_h \rrbracket, \quad E_h^2(q_h) \stackrel{\text{def}}{=} \tilde{q}_h$$

for all $q_h \in Q_h$. We recall that by construction of the pressure space $Q_h$ in (6), the pressure jump $\llbracket q_h \rrbracket$ is constant. In particular, we have

$$E_h^i(q_h) = q_{h,i} \quad \text{in } \Omega_i. \qquad (46)$$

Let $K \in \mathcal{T}_h^c$ be given. Since, $\llbracket q_h \rrbracket$ is constant on $\Sigma$, using Assumption 4.1 and the triangle inequality, we have

$$h^{\frac{1}{2}} \|\llbracket q_h \rrbracket\|_{L^2(\Sigma \cap K)} = h^{\frac{1}{2}} |\Sigma \cap K|^{\frac{1}{2}} \|\llbracket q_h \rrbracket\|_{L^2(K)} = \|E_h^1(q_h) - E_h^2(q_h)\|_{L^2(K)} \le \sum_i \|E_h^i(q_h)\|_{L^2(K)}. \quad (47)$$

In order to bound $\|E_h^i(q_h)\|_{L^2(K)}$ for every $i \in \{1, 2\}$, we first note that since $E_h^i(q_h)$ is affine in $K$, we have

$$E_h^i(q_h)(\boldsymbol{x}) = E_h^i(q_h)(\boldsymbol{x}_{\kappa,j}) + \nabla E_h^i(q_h) \cdot (\boldsymbol{x} - \boldsymbol{x}_{\kappa,j}),$$

where $\boldsymbol{x}_{\kappa,j}$ denotes one of the nodes of the element $K$. Thus, using an inverse inequality we get

$$\|E_h^i(q_h)\|_{L^2(K)} \lesssim h^{\frac{d}{2}} \|E_h^i(q_h)\|_{L^\infty(K)} \lesssim h^{\frac{d}{2}} |E_h^i(q_h)(\boldsymbol{x}_{\bar{\kappa}_i, l})| + h^{1+\frac{d}{2}} \|\nabla E_h^i(q_h)\|_{L^\infty(K)}. \quad (48)$$

Furthermore, owing to Assumption 4.2, relation (46) and [24, Proposition 12.5], we have

$$h^{\frac{d}{2}} |E_h^i(q_h)(\boldsymbol{x}_{\kappa,j})| \le h^{\frac{d}{2}} \max_l \|q_{h,i}(\boldsymbol{x}_{\bar{\kappa}_i,l})\| \lesssim \|q_{h,i}\|_{L^2(\bar{\kappa}_i)}. \quad (49)$$

In addition, since $\nabla E_h^i(q_h)$ is constant in $K$ and $\nabla E_h^1(q_h) = \nabla E_h^2(q_h) = \nabla \tilde{q}_h$, we have

$$h^{\frac{d}{2}} \|\nabla E_h^i(q_h)\|_{L^\infty(K)} \lesssim \|\nabla \tilde{q}_h\|_{L^2(K)}. \quad (50)$$

By inserting the bounds (48)-(50) into (47), we have

$$h^{\frac{1}{2}} \|\llbracket q_h \rrbracket\|_{L^2(K \cap \Sigma)} \lesssim \sum_i \|E_h^i(q_h)\|_{L^2(K)} \lesssim (\|q_h\|_{L^2(\bar{\kappa}_1)} + \|q_h\|_{L^2(\bar{\kappa}_2)} + h \|\nabla \tilde{q}_h\|_{L^2(K)})^{\frac{1}{2}}.$$

The estimate (45) then follows by exploiting the quasi-uniformity of the mesh and summing the above relation over all $K \in \mathcal{T}_h^c$, which completes the proof. $\square$

### 4.2.2. Main stability result

We now have all the ingredients to state the main result of this section. To this purpose, we introduce the following norm

$$\|(v_h, q_h, \boldsymbol{\xi}_H)\| \stackrel{\text{def}}{=} \left( \mu \|\boldsymbol{\varepsilon}(v_h)\|_{L^2(\Omega)}^2 + h^2 \mu^{-1} \|\nabla \tilde{q}_h\|_{L^2(\Omega)}^2 + \mu^{-1} \|q_h\|_{L^2(\Omega)}^2 + h \mu^{-1} \|\boldsymbol{\xi}_H\|_{L^2(\Sigma)}^2 \right)^{\frac{1}{2}}$$

for all $(v_h, q_h, \boldsymbol{\xi}_H) \in V_h \times Q_h \times \Lambda_H$. The following lemma provides the inf-sup stability for the full bilinear form $\mathcal{A}_h$.

**Lemma 4.8.** Let Assumptions 2.1, 4.1 and 4.2 be fulfilled. For $\theta = 0$, we furthermore assume that $\gamma_\lambda$ is large enough. Then, there exists a constant $\beta > 0$ independent of $\mu$, the mesh and the interface position such that

$$(u_h, p_h, \boldsymbol{\lambda}_H) \in \mathcal{V}_h \times Q_h \times \Lambda_H \setminus (0,0,0) \quad (v_h, q_h, \boldsymbol{\xi}_H) \in \mathcal{V}_h \times Q_h \times \Lambda_H \setminus (0,0,0) \quad \frac{\mathcal{A}_h((u_h, p_h, \boldsymbol{\lambda}_H), (v_h, q_h, \boldsymbol{\xi}_H))}{\|(u_h, p_h, \boldsymbol{\lambda}_H)\| \|(v_h, q_h, \boldsymbol{\xi}_H)\|} \ge \beta. \quad (51)$$
*Proof.* The inf-sup condition (51) is a consequence of the following estimate

$$ \|\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}}\| \lesssim \underbrace{\sup_{(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}}) \in \mathbf{V}_h \times Q_h \times \mathbf{\Lambda}_{\mathcal{H}} \setminus (\mathbf{0}, \mathbf{0}, \mathbf{0})}}_{\text{def}=S} \frac{\mathcal{A}_h((\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}}), (v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}}))}{\|(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}})\|} \quad (52) $$

for all $(\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}}) \in \mathbf{V}_h \times Q_h \times \mathbf{\Lambda}_{\mathcal{H}} \setminus (\mathbf{0}, \mathbf{0}, \mathbf{0})$. Now, let us prove (52). Using Assumption 2.1, we can invoke Proposition 4.6. Let $\mathbf{w}_h \in \mathbf{V}_h \cap \mathbf{H}_{\Sigma}^{\mathcal{H}}(\Omega)$ be a function for which the supremum in (36) is attained (for $q_h = p_h$). We have

$$ \|p_h\|_{L^2(\Omega)} \lesssim \frac{(p_h, \nabla \cdot \mathbf{w}_h)_{\Omega}}{\|\mathbf{w}_h\|_{H^1(\Omega)}} + h \|\nabla \tilde{p}_h\|_{L^2(\Omega)}. \quad (53) $$

In addition, from (11), we have

$$ \begin{aligned} (p_h, \nabla \cdot \mathbf{w}_h)_{\Omega} &= \mathcal{A}_h((\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}}), (-\mathbf{w}_h, \mathbf{0}, \mathbf{0})) + 2\mu(\boldsymbol{\varepsilon}(\mathbf{u}_h), \boldsymbol{\varepsilon}(\mathbf{w}_h))_{\Omega} \\ &\lesssim \mathcal{S} \|(\mathbf{w}_h, \mathbf{0}, \mathbf{0})\| + 2\mu \|\boldsymbol{\varepsilon}(\mathbf{u}_h)\|_{L^2(\Omega)} \|\boldsymbol{\varepsilon}(\mathbf{w}_h)\|_{L^2(\Omega)} \\ &\lesssim (\mu^{\frac{1}{2}} \mathcal{S} + 2\mu \|\boldsymbol{\varepsilon}(\mathbf{u}_h)\|_{L^2(\Omega)}) \|\mathbf{w}_h\|_{H^1(\Omega)}, \end{aligned} $$

so that (53) yields

$$ \|p_h\|_{L^2(\Omega)} \lesssim \mu^{\frac{1}{2}} \mathcal{S} + 2\mu \|\boldsymbol{\varepsilon}(\mathbf{u}_h)\|_{L^2(\Omega)} + h \|\nabla \tilde{p}_h\|_{L^2(\Omega)}. \quad (54) $$

Case $\theta = 1$. Let $(\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}}) \in \mathbf{V}_h \times Q_h \times \mathbf{\Lambda}_{\mathcal{H}} \setminus (\mathbf{0}, \mathbf{0}, \mathbf{0})$. From (11), we have

$$ \begin{aligned} 2\mu \|\boldsymbol{\varepsilon}(\mathbf{u}_h)\|_{L^2(\Omega)}^2 &+ \frac{\gamma_p h^2}{\mu} \|\nabla \tilde{p}_h\|_{L^2(\Omega)}^2 + \frac{h}{\gamma_\lambda \mu} \|\boldsymbol{\lambda}_{\mathcal{H}} + \llbracket p_h \rrbracket \mathbf{n}_\Sigma\|_{L^2(\Sigma)}^2 \\ & = \mathcal{A}_h((\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}}), (\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}})) \quad (55) \\ & \leq \mathcal{S} \|(\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}})\|. \end{aligned} $$

Moreover, using the triangle inequality and Proposition 4.7 (under Assumptions 4.1 and 4.2), we get

$$ \begin{aligned} h \|\boldsymbol{\lambda}_{\mathcal{H}}\|_{L^2(\Sigma)}^2 &\leq h \|\boldsymbol{\lambda}_{\mathcal{H}} + \llbracket p_h \rrbracket \mathbf{n}_\Sigma\|_{L^2(\Sigma)}^2 + h \|\llbracket p_h \rrbracket\|_{L^2(\Sigma)}^2 \\ &\lesssim h \|\boldsymbol{\lambda}_{\mathcal{H}} + \llbracket p_h \rrbracket \mathbf{n}_\Sigma\|_{L^2(\Sigma)}^2 + \|p_h\|_{L^2(\Omega)}^2 + h^2 \|\nabla \tilde{p}_h\|_{L^2(\Omega)}^2. \end{aligned} \quad (56) $$

We now conclude by inserting (54) and (56) (multiplied by $\alpha_p^{\frac{1}{2}}\mu^{-\frac{1}{2}}$ and $\alpha_\lambda\mu^{-1}\gamma_\lambda^{-1}$, respectively for $\alpha_p, \alpha_\lambda$ two positive coefficients) into (55), which yields

$$ C_1(\gamma_p, \gamma_\lambda, \alpha_p, \alpha_\lambda) \|(\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}})\|^2 \lesssim \alpha_p \mathcal{S}^2 + \mathcal{S} \|(\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}})\|, $$

with

$$ C_1(\gamma_p, \gamma_\lambda, \alpha_p, \alpha_\lambda) \stackrel{\text{def}}{=} \min \left\{ 2 - 4\alpha_p, \gamma_p - \alpha_p - \frac{\alpha_\lambda}{\gamma_\lambda}, \alpha_p - \frac{\alpha_\lambda}{\gamma_\lambda}, \frac{\alpha_\lambda}{\gamma_\lambda} \right\}, $$

which is positive for an adequate choice of $\alpha_p$ and $\alpha_\lambda$. Then (52) simply follows from Young's inequality, which completes the proof.

Case $\theta = 0$. Similarly to the previous case, let $(\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}}) \in \mathbf{V}_h \times Q_h \times \mathbf{\Lambda}_{\mathcal{H}} \setminus (\mathbf{0}, \mathbf{0}, \mathbf{0})$. From (11), we have

$$ \begin{aligned} 2\mu \|\boldsymbol{\varepsilon}(\mathbf{u}_h)\|_{L^2(\Omega)}^2 &+ \frac{\gamma_p h^2}{\mu} \|\nabla \tilde{p}_h\|_{L^2(\Omega)}^2 + \frac{h}{\gamma_\lambda \mu} (\boldsymbol{\lambda}_{\mathcal{H}} + \llbracket p_h \rrbracket \mathbf{n}_\Sigma, \boldsymbol{\lambda}_{\mathcal{H}})_\Sigma \\ & = \mathcal{A}_h((\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}}), (\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}})) \quad (57) \\ & \leq \mathcal{S} \|(\mathbf{u}_h, p_h, \boldsymbol{\lambda}_{\mathcal{H}})\|. \end{aligned} $$
Moreover, using Cauchy-Schwarz and Young inequalities and Proposition 4.7 (under Assumptions 4.1 and 4.2), we get

$$ \frac{h}{2} \|\lambda_H\|_{L^2(\Sigma)}^2 \le h(\lambda_H + \|p_h\|_H \mathbf{n}_\Sigma, \lambda_H)_{L^2(\Sigma)} + \frac{C_\lambda}{2} \|p_h\|_{L^2(\Omega)}^2 + \frac{C_\lambda h^2}{2} \|\nabla \tilde{p}_h\|_{L^2(\Omega)}^2. \quad (58) $$

We now conclude by inserting (54) and (58) (multiplied respectively by $\alpha_p^{-\frac{1}{2}}\mu^{-\frac{1}{2}}$ and $\gamma_\lambda^{-1}\mu^{-1}$ for $\alpha_p$ a positive coefficient) into (57), which yields

$$ C_0(\gamma_p, \gamma_\lambda, \alpha_p) \|(u_h, p_h, \lambda_H)\|^2 \lesssim \alpha_p S^2 + S \|(u_h, p_h, \lambda_H)\|, $$

with

$$ C_0(\gamma_p, \gamma_\lambda, \alpha_p) \stackrel{\text{def}}{=} \min \left\{ 2 - 4\alpha_p, \gamma_p - \alpha_p - \frac{C_\lambda}{2\gamma_\lambda}, \alpha_p - \frac{C_\lambda}{2\gamma_\lambda}, \frac{1}{2\gamma_\lambda} \right\}, $$

where $C_\lambda > 0$ is the hidden constant in the statement of Proposition 4.7. The constant $C_0(\gamma_p, \gamma_\lambda, \alpha_p)$ is positive for an adequate choice of $\alpha_p$, under the assumption that $\gamma_\lambda$ is large enough. Then (52) simply follows from Young's inequality, which completes the proof. $\square$

**Corollary 4.9.** The discrete problem (10) admits a unique solution.

**Remark 4.10.** In the case where $\theta = 0$, the requirement for $\gamma_\lambda$ to be large enough stems from the definition of the constant $C_0(\gamma_p, \gamma_\lambda, \alpha_p)$ (as detailed in the proof). This constitutes a limitation of the $\theta = 0$ variant of the proposed method. In practical terms, we aim for a $\gamma_\lambda$ that is small enough to effectively ensure the stabilization of the Lagrange multiplier, while simultaneously being large enough to prevent excessive perturbation of the mass conservation (see Remark 3.2).

### 4.3. A priori error estimates

We first introduce the following discrete errors

$$ e_{u,h} = \Pi_h u - u_h, \quad e_{p,h} = \Pi_h p - p_h, \quad e_{\lambda,\mathcal{H}} = \Pi_{\mathcal{H}} \lambda - \lambda_{\mathcal{H}}, $$

where $(u, p, \lambda)$ is the solution to the continuous problem (2), $(u_h, p_h, \lambda_H)$ is the solution to the discrete problem (10) and the interpolation operators $\Pi_h$, $\Pi_h$ and $\Pi_{\mathcal{H}}$ are defined in Section 4.1. The following result provides an estimate for the discrete error.

**Lemma 4.11.** Assume that $(u, p, \lambda) \in V_\varepsilon \times Q_\varepsilon \times H^{\frac{1}{2}-\varepsilon}(\Sigma)$ with $\varepsilon \in (0, 1/2)$. The following estimate holds

$$ |\mathcal{A}_h((e_{u,h}, e_{p,h}, e_{\lambda,\mathcal{H}}), (v_h, q_h, \xi_{\mathcal{H}}))| \lesssim \max\{h, \mathcal{H}\}^{\frac{1}{2}-\varepsilon} \left( \mu^{\frac{1}{2}} \|u\|_{H^{\frac{3}{2}-\varepsilon}(\Omega)} + \mu^{-\frac{1}{2}} \|p - \hat{\Pi}_h p\|_{H^{\frac{1}{2}-\varepsilon}(\Omega)} \right. \\ \left. + \mu^{-\frac{1}{2}} \|\|p\|\|_{L^2(\Sigma)} + \mu^{-\frac{1}{2}} \|\lambda\|_{H^{\frac{1}{2}-\varepsilon}(\Sigma)} \right) \|v_h, q_h, \xi_{\mathcal{H}}\| \quad (59) $$

for all $(v_h, q_h, \xi_{\mathcal{H}}) \in V_h \times Q_h \times \Lambda_{\mathcal{H}}$. Moreover, this estimate can be improved to

$$ |\mathcal{A}_h((e_{u,h}, e_{p,h}, e_{\lambda,h}), (v_h, q_h, \xi_{\mathcal{H}}))| \\ \lesssim \left[ \max\{h, \mathcal{H}\}^{\frac{1}{2}-\varepsilon} \left( \mu^{\frac{1}{2}} \|u\|_{H^{\frac{3}{2}-\varepsilon}(\Omega)} + \mu^{-\frac{1}{2}} \|p - \hat{\Pi}_h p\|_{H^{\frac{1}{2}-\varepsilon}(\Omega)} + \mu^{-\frac{1}{2}} \|\lambda + \|\hat{\Pi}_h p\|_H \mathbf{n}_\Sigma\|_{H^{\frac{1}{2}-\varepsilon}(\Sigma)} \right) \right. \\ \left. + \mu^{-\frac{1}{2}} \|\|p\|\|_{L^2(\Sigma)} \|n_\Sigma - \Pi_{\mathcal{H}} n_\Sigma\|_{L^2(\Sigma)} \right] \|v_h, q_h, \xi_{\mathcal{H}}\|. \quad (60) $$
*Proof.* Using the relations (2) and (10) and the definitions (3), (11) and (19), we have

$$
\begin{align}
\mathcal{A}_h((e_{u,h}, e_{p,h}, e_{\lambda,\mathcal{H}}), (v_h, q_h, \xi_{\mathcal{H}})) & = \mathcal{A}_h((\Pi_h u, \Pi_h p, \Pi_{\mathcal{H}} \lambda), (v_h, q_h, \xi_{\mathcal{H}})) - \mathcal{F}(v_h, \xi_{\mathcal{H}}) \nonumber \\
& = \mathcal{A}_h((\Pi_h u, \Pi_h p, \Pi_{\mathcal{H}} \lambda), (v_h, q_h, \xi_{\mathcal{H}})) - \mathcal{A}((u, p, \lambda), (v_h, q_h, \xi_{\mathcal{H}})) \tag{61} \\
& = -2\mu (\varepsilon(e_u), \varepsilon(v_h))_{\Omega} + (e_p, \nabla \cdot v_h)_{\Omega} - (q_h, \nabla \cdot e_u)_{\Omega} + (e_{\lambda}, v_h)_{\Sigma} - (\xi_{\mathcal{H}}, e_u)_{\Sigma} \nonumber \\
& \quad + s_h^{\text{BP}}(\tilde{\Pi}_h p, \tilde{q}_h) + s_h^{\text{BH}}((\Pi_h p, \Pi_{\mathcal{H}} \lambda), (q_h, \xi_{\mathcal{H}})). \nonumber
\end{align}
$$

We now proceed by estimating term by term. Using the Cauchy-Schwarz inequality and the approximation error bounds (20), (21) and (23), we have

$$
\begin{align}
\mu | (\boldsymbol{\varepsilon}(e_u), \boldsymbol{\varepsilon}(v_h))_{\Omega} | &\le \mu \|\boldsymbol{\varepsilon}(e_u)\|_{L^2(\Omega)} \|\boldsymbol{\varepsilon}(v_h)\|_{L^2(\Omega)} \lesssim \mu^{\frac{1}{2}} h^{\frac{1}{2}-\varepsilon} \|\boldsymbol{u}\|_{H^{\frac{3}{2}-\varepsilon}(\Omega)} \|(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}})\|, \nonumber \\
| (q_h, \nabla \cdot e_u)_{\Omega} | &\le \|q_h\|_{L^2(\Omega)} \|\nabla \cdot e_u\|_{L^2(\Omega)} \lesssim \mu^{\frac{1}{2}} h^{\frac{1}{2}-\varepsilon} \|\boldsymbol{u}\|_{H^{\frac{3}{2}-\varepsilon}(\Omega)} \|(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}})\|, \tag{62} \\
| (\boldsymbol{\xi}_{\mathcal{H}}, e_u)_{\Sigma} | &\le \|\boldsymbol{\xi}_{\mathcal{H}}\|_{L^2(\Sigma)} \|e_u\|_{L^2(\Sigma)} \lesssim \mu^{\frac{1}{2}} h^{\frac{1}{2}-\varepsilon} \|\boldsymbol{u}\|_{H^{\frac{3}{2}-\varepsilon}(\Omega)} \|(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}})\|, \nonumber \\
| (e_p, \nabla \cdot v_h)_{\Omega} | &\le \|e_p\|_{L^2(\Omega)} \|\nabla \cdot v_h\|_{L^2(\Omega)} \lesssim \mu^{-\frac{1}{2}} h^{\frac{1}{2}-\varepsilon} \|p - \hat{\Pi}_h p\|_{H^{\frac{1}{2}-\varepsilon}(\Omega)} \|(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}})\|. \nonumber
\end{align}
$$

For the pressure stabilization term, we use Cauchy-Schwarz inequality and (17) to get

$$
\begin{align}
|s_h^{\text{BP}}(\tilde{\Pi}_h p, \tilde{q}_h)| & \lesssim \mu^{-1} h^2 |(\nabla \tilde{\Pi}_h p, \nabla \tilde{q}_h)_\Omega| \lesssim \mu^{-1} h \|\nabla \tilde{\Pi}_h p\|_{L^2(\Omega)} h \|\nabla \tilde{q}_h\|_{L^2(\Omega)} \nonumber \\
& \lesssim \mu^{-\frac{1}{2}} h^{\frac{1}{2}-\varepsilon} \|p - \hat{\Pi}_h p\|_{H^{\frac{1}{2}-\varepsilon}(\Omega)} \|(v_h, q_h, \xi_{\mathcal{H}})\|. \tag{63}
\end{align}
$$

It only remains to estimate the terms $(e_\lambda, v_h)_\Sigma$ and $s_h^{\text{BH}}((\Pi_h p, \Pi_\lambda \lambda), (q_h, \xi_\lambda))$.

*Proof of (59).* For the term $(e_\lambda, v_h)_\Sigma$, using the Cauchy-Schwarz, trace and Korn's inequalities and (22), we have

$$
| (\Pi_{\mathcal{H}} \lambda - \lambda, v_h)_{\Sigma} | \lesssim \mathcal{H}^{\frac{1}{2}-\epsilon} \|\lambda\|_{H^{\frac{1}{2}-\epsilon}(\Sigma)} \|\boldsymbol{\varepsilon}(v_h)\|_{L^2(\Omega)} \lesssim \mu^{-\frac{1}{2}} \mathcal{H}^{\frac{1}{2}-\epsilon} \|\lambda\|_{H^{\frac{1}{2}-\epsilon}(\Sigma)} \|(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}})\|. \quad (64)
$$

Finally, for the term $s_h^{\text{BH}}((\Pi_h p, \Pi_\lambda \lambda), (q_h, \xi_\lambda))$, we use once more the Cauchy-Schwarz inequality, Proposition 4.7 (for $\theta = 1$), (15) and (16), to get

$$
\begin{align}
\mu^{-1} h \left| \left( \Pi_{\mathcal{H}} \lambda + \left[ \Pi_h p \right]_{\Sigma} \mathbf{n}_{\Sigma} + \theta \left[ q_h \right]_{\Sigma} \mathbf{n}_{\Sigma} \right)_{\Sigma} \right| & \lesssim \mu^{-1} h^{\frac{1}{2}} \left( \| \Pi_{\mathcal{H}} \lambda \|_{L^2(\Sigma)} + \| \left[ \Pi_h p \right]_{\Sigma} \| \right) h^{\frac{1}{2}} \left( \| \boldsymbol{\xi}_{\mathcal{H}} \|_{L^2(\Sigma)} + \theta \| \left[ q_h \right]_{\Sigma} \| \right) \tag{65} \\
& \lesssim \mu^{-\frac{1}{2}} h^{\frac{1}{2}} \left( \| \boldsymbol{\lambda} \|_{L^2(\Sigma)} + \| \left[ p \right]_{\Sigma} \| \right) \|(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}})\|. \nonumber
\end{align}
$$

The estimate (59) then follows from (61) and the estimates (62)–(65).

*Proof of (60).* We first recall that, from (14), $\|\hat{\Pi}_h p\|$ is constant. By adding and subtracting suitable terms and by using the Cauchy-Schwarz, triangle, trace and Korn's inequalities and (24), we have

$$
\begin{align}
| (\Pi_{\mathcal{H}} \lambda - \lambda, v_h)_{\Sigma} | & \le \| \Pi_{\mathcal{H}} (\lambda + \|\hat{\Pi}_h p\| \mathbf{n}_{\Sigma}) - (\lambda + \|\hat{\Pi}_h p\| \mathbf{n}_{\Sigma}) + \|\hat{\Pi}_h p\| (\mathbf{n}_{\Sigma} - \Pi_{\mathcal{H}} \mathbf{n}_{\Sigma}) \|_{L^2(\Sigma)} \|v_h\|_{L^2(\Sigma)} \nonumber \\
& \lesssim \mu^{-\frac{1}{2}} (\mathcal{H}^{\frac{1}{2}-\epsilon} \|\lambda + \|\hat{\Pi}_h p\| \mathbf{n}_{\Sigma} \|_{H^{\frac{1}{2}-\epsilon}(\Sigma)} + \|\hat{\Pi}_h p\| \|(\mathbf{n}_{\Sigma} - \Pi_{\mathcal{H}} \mathbf{n}_{\Sigma})\|_{L^2(\Sigma)} ) \|(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}})\|. \tag{66}
\end{align}
$$
Similarly, for the last term we use the Cauchy-Schwarz inequality, Proposition 4.7 (for $\theta = 1$), (18) and (15), which yields

$$ 
\begin{aligned}
& \mu^{-1}h | ( \mathbf{\Pi}_{\mathcal{H}} \boldsymbol{\lambda} + \|\Pi_h p\| \mathbf{n}_{\Sigma}, \boldsymbol{\xi}_{\mathcal{H}} + \theta \|q_h\| \mathbf{n}_{\Sigma})_{\Sigma} | \\
& \lesssim \mu^{-1}h^{\frac{1}{2}} \| \mathbf{\Pi}_{\mathcal{H}} (\boldsymbol{\lambda} + \|\hat{\Pi}_h p\| \mathbf{n}_{\Sigma}) + \|\hat{\Pi}_h p\| (\mathbf{n}_{\Sigma} - \mathbf{\Pi}_{\mathcal{H}} \mathbf{n}_{\Sigma}) \|_{L^2(\Sigma)} h^{\frac{1}{2}} ( \|\boldsymbol{\xi}_{\mathcal{H}}\|_{L^2(\Sigma)} + \theta \|q_h\|_{L^2(\Sigma)} ) \\
& \lesssim \mu^{-\frac{1}{2}} h^{\frac{1}{2}} ( \|\boldsymbol{\lambda} + \|\hat{\Pi}_h p\| \mathbf{n}_{\Sigma}\|_{L^2(\Sigma)} + \|\|\hat{\Pi}_h p\|\| \mathbf{n}_{\Sigma} - \mathbf{\Pi}_{\mathcal{H}} \mathbf{n}_{\Sigma}\|_{L^2(\Sigma)} ) \|(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}})\| \\
& \lesssim \mu^{-\frac{1}{2}} h^{\frac{1}{2}} ( \|\boldsymbol{\lambda} + \|\hat{\Pi}_h p\| \mathbf{n}_{\Sigma}\|_{L^2(\Sigma)} + \|\|p\|\|_{L^2(\Sigma)} \|\mathbf{n}_{\Sigma} - \mathbf{\Pi}_{\mathcal{H}} \mathbf{n}_{\Sigma}\|_{L^2(\Sigma)} ) \|(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}})\|.
\end{aligned} 
\quad (67) $$

The bound (60) then follows from (61) and the estimates (62), (63), (66) and (67). This completes the proof. □

We now have all the ingredients to derive a priori error estimates. This is stated in the next theorem.

**Theorem 4.12.** Let Assumptions 2.1, 4.1 and 4.2 be fulfilled. Let $\varepsilon \in (0, 1/2)$ and assume that $(u, p, \lambda) \in V_{\varepsilon} \times Q_{\varepsilon} \times H^{\frac{1}{2}-\varepsilon}(\Sigma)$ is the solution to the continuous problem (2) and that $(u_h, p_h, \lambda_{\mathcal{H}})$ is the solution to the discrete problem (10). Furthermore, in the case $\theta = 0$, assume that $\gamma_{\lambda}$ is large enough. The following a priori error bounds hold:

$$ 
\begin{aligned}
& \mu^{\frac{1}{2}} \|\boldsymbol{\varepsilon}(u - u_h)\|_{L^2(\Omega)} + \mu^{-\frac{1}{2}} \|p - p_h\|_{L^2(\Omega)} + \mu^{-\frac{1}{2}} h^{\frac{1}{2}} \|\boldsymbol{\lambda} - \boldsymbol{\lambda}_{\mathcal{H}}\|_{L^2(\Sigma)} \\
& \lesssim \max\{h, \mathcal{H}\}^{\frac{1}{2}-\varepsilon} \left( \mu^{\frac{1}{2}} \|\boldsymbol{u}\|_{H^{\frac{3}{2}-\varepsilon}(\Omega)} + \mu^{-\frac{1}{2}} \|p - \hat{\Pi}_h p\|_{H^{\frac{1}{2}-\varepsilon}(\Omega)} + \mu^{-\frac{1}{2}} \|\|p\|\|_{L^2(\Sigma)} + \mu^{-\frac{1}{2}} \|\boldsymbol{\lambda}\|_{H^{\frac{1}{2}-\varepsilon}(\Sigma)} \right), \quad (68)
\end{aligned} $$

$$ 
\begin{aligned}
& \mu^{\frac{1}{2}} \|\boldsymbol{\varepsilon}(u - u_h)\|_{L^2(\Omega)} + \mu^{-\frac{1}{2}} \|p - p_h\|_{L^2(\Omega)} + \mu^{-\frac{1}{2}} h^{\frac{1}{2}} \|\boldsymbol{\lambda} - \boldsymbol{\lambda}_{\mathcal{H}}\|_{L^2(\Sigma)} \\
& \lesssim \max\{h, \mathcal{H}\}^{\frac{1}{2}-\varepsilon} \left( \mu^{\frac{1}{2}} \|\boldsymbol{u}\|_{H^{\frac{3}{2}-\varepsilon}(\Omega)} + \mu^{-\frac{1}{2}} \|p - \hat{\Pi}_h p\|_{H^{\frac{1}{2}-\varepsilon}(\Omega)} + \mu^{-\frac{1}{2}} \|\boldsymbol{\lambda} + \|\hat{\Pi}_h p\| \mathbf{n}_{\Sigma}\|_{H^{\frac{1}{2}-\varepsilon}(\Sigma)} \right) \\
& \qquad + \mu^{-\frac{1}{2}} \|\|p\|\|_{L^2(\Sigma)} \|\mathbf{n}_{\Sigma} - \mathbf{\Pi}_{\mathcal{H}} \mathbf{n}_{\Sigma}\|_{L^2(\Sigma)}.
\end{aligned} 
\quad (69) $$

*Proof.* The estimate (68) follows immediately using the classical error decomposition

$$ 
\begin{aligned}
& \mu^{\frac{1}{2}} \|\boldsymbol{\varepsilon}(u - u_h)\|_{L^2(\Omega)} + \mu^{-\frac{1}{2}} \|p - p_h\|_{L^2(\Omega)} + \mu^{-\frac{1}{2}} h^{\frac{1}{2}} \|\boldsymbol{\lambda} - \boldsymbol{\lambda}_{\mathcal{H}}\|_{L^2(\Sigma)} \\
& \leq \mu^{\frac{1}{2}} \|\boldsymbol{\varepsilon}(e_u)\|_{L^2(\Omega)} + \mu^{-\frac{1}{2}} \|e_p\|_{L^2(\Omega)} + \mu^{-\frac{1}{2}} h^{\frac{1}{2}} \|e_{\lambda}\|_{L^2(\Sigma)} + \|(e_{u,h}, e_{p,h}, e_{\lambda,\mathcal{H}})\|.
\end{aligned} $$

Indeed, using Lemmas 4.8 and 4.11, we have

$$ 
\begin{aligned}
\|e_{u,h}, e_{p,h}, e_{\lambda,\mathcal{H}}\| & \lesssim \sup_{(v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}}) \in V_h \times Q_h \times \mathbf{\Lambda}_{\mathcal{H}} \setminus (0,0,0)} \frac{A_h((e_{u,h}, e_{p,h}, e_{\lambda,\mathcal{H}}); (v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}}))}{\|v_h, q_h, \boldsymbol{\xi}_{\mathcal{H}}\|} \\
& \lesssim \max\{h, \mathcal{H}\}^{\frac{1}{2}-\varepsilon} (\mu^{\frac{1}{2}} \|\boldsymbol{u}\|_{H^{\frac{3}{2}-\varepsilon}(\Omega)} + \mu^{-\frac{1}{2}} \|p - \hat{\Pi}_h p\|_{H^{\frac{1}{2}-\varepsilon}(\Omega)} \\
& \qquad + \mu^{-\frac{1}{2}} \|\|p\|\|_{L^2(\Sigma)} + \mu^{-\frac{1}{2}} \|\boldsymbol{\lambda}\|_{H^{\frac{1}{2}-\varepsilon}(\Sigma)}).
\end{aligned} $$

The estimate (68) then follows from (20), (22) and (23). Finally, the error estimate (69) follows similarly, by using (60) instead of (59). This completes the proof. □

Estimate (68) corresponds to the expected sub-optimal a priori error bound. Since the degrees of freedom are not doubled in the cut cells (like in XFEM), we cannot hope for a better convergence rate which exploits the local regularity (5) in each of the two sub-domains. Numerical evidence of this is provided in Section 5.1. Estimate (69) is an improved version of (68) in the presence of high pressure drops through the interface. Indeed,
the mean value of the pressure jump $\|\hat{\Pi}_h p\|_{\Sigma}$ is subtracted to the Lagrange multiplier, which is expected to improve the constant of the error bound. Moreover, the term $\|\|p\|\|_{L^2(\Sigma)} \|n_{\Sigma} - \Pi_{\mathcal{H}} n_{\Sigma}\|_{L^2(\Sigma)}$ that carries the pressure jump has an improved convergence rate in this estimate, assuming that $n_{\Sigma}$ has a better regularity than $H^{\frac{1}{2}-\varepsilon}(\Sigma)$. The numerical results of Section 5.2.1 illustrate this improved behavior.

**Remark 4.13.** The result stated in Theorem 4.12 does not provide convergence for the Lagrange multiplier in $L^2$-norm. We could however prove convergence of the multiplier in its natural norm, i.e. $H^{-\frac{1}{2}}(\Sigma)$. Elements of such a proof are given in [3, Propositions A.1, A.2].

## 5. NUMERICAL EXPERIMENTS

In this section, we illustrate the behavior of the numerical method (10) in three different two-dimensional examples. In all these test cases, we consider a rectangular domain $\Omega = (-1, 1) \times (0, 1)$. Dirichlet boundary conditions are imposed at the top and bottom boundaries $\Gamma_D = (-1, 1) \times \{0, 1\}$, and Neumann boundary conditions are enforced on the remaining boundary $\Gamma_N = \{-1, 1\} \times (0, 1)$. A different interface $\Sigma$ is considered for each test case, as shown in Figure 3. Without loss of generality, the interface is oriented so that the sub-domain $\Omega_1$ (resp. $\Omega_2$) lies on the left (resp. right) side of the interface.

<figure><img src="image_6.png" /><figcaption>FIGURE 3. Geometric configurations considered in the numerical tests.</figcaption></figure>

In the subsequent sections, all physical quantities and parameters are expressed in the centimeter-gram-second (CGS) unit system, the stabilization parameters are $\gamma_p = 10^{-2}$ and $\gamma_\lambda = 10$.

In order to highlight the impact of the additional velocity constraint in the accuracy of the proposed method (or of the additional degree of freedom for the pressure), we compare the results with those of its corresponding unenriched counterpart (i.e., without the additional discontinuous pressure basis function). This method, referred to as FD-LM in the sequel, simply consists in the discrete formulation (10) with the pressure space $\tilde{Q}_h$ instead of $Q_h$, namely: Find $(u_h, p_h, \lambda_{\mathcal{H}}) \in V_h \times \tilde{Q}_h \times \Lambda_{\mathcal{H}}$ such that

$$ \mathcal{A}_h((u_h, p_h, \lambda_{\mathcal{H}}), (v_h, q_h, \xi_{\mathcal{H}})) = \mathcal{F}(v_h, \xi_{\mathcal{H}}) $$

for all $(v_h, q_h, \xi_{\mathcal{H}}) \in V_h \times \tilde{Q}_h \times \Lambda_{\mathcal{H}}$. It should be noted that for this method the Lagrange multiplier stabilization term (9) reduces to the one considered in [7], namely,

$$ s_h^{\text{BH}}(\lambda_{\mathcal{H}}, \xi_{\mathcal{H}}) \stackrel{\text{def}}{=} \frac{h}{\gamma_{\lambda}\mu}(\lambda_{\mathcal{H}}, \xi_{\mathcal{H}})_{\Sigma}. $$

### 5.1. Convergence study for a non-trivial solution

This first example aims at providing numerical evidence of the convergence rate stated in Theorem 4.12. For the sake of simplicity, a straight interface $\Sigma = \{0\} \times (0, 1)$ is considered (see Figure 3a). A non-trivial analytical solution is considered and numerical errors are computed separately for the velocity, the pressure and the Lagrange multiplier.
<figure><img src="image_7.png" /><figcaption>(a) Velocity relative energy error.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Pressure relative approximation error.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Lagrange multiplier relative approximation error.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) Mass loss in $\Omega_1$. The results for $\theta = 0$ are omitted since this variant guarantees interfacial mass conservation.</figcaption></figure>

<caption>FIGURE 4. Convergence histories of the velocity, pressure, Lagrange multiplier and mass loss in $\Omega_1$ obtained with the present method ($\theta = 0, \theta = 1$) and the FD-LM method.</caption>

For $i \in \{1, 2\}$, the velocity field $\mathbf{u}_i \stackrel{\text{def}}{=} \mathbf{u}|_{\Omega_i}$ in (1) is chosen as

$$ \mathbf{u}_1 = \left( \frac{\partial}{\partial y} \phi_1 (x, y), -\frac{\partial}{\partial x} \phi_1 (x, y) \right), \quad \mathbf{u}_2 = \left( \frac{\partial}{\partial y} \phi_2 (x, y), -\frac{\partial}{\partial x} \phi_2 (x, y) \right), $$

with the function $\phi_i(x, y)$ defined as

$$ \phi_i (x, y) \stackrel{\text{def}}{=} (x - x_{l,i})^2 (x - x_{r,i})^2 (y - y_{b,i})^2 (y - y_{t,i})^2 \varphi_i (x, y), $$

with $\varphi_1(x, y) \stackrel{\text{def}}{=} x$, $\varphi_2(x, y) \stackrel{\text{def}}{=} y$ and the parameters ($x_{l,i}$, $x_{r,i}$, $y_{b,i}$, $y_{t,i}$) are chosen as $(-1, 0, 0, 1)$ and $(0, 1, 0, 1)$ for $i=1$ and $i=2$, respectively. The resulting analytical solution to problem (1) reads:

$$ \mathbf{u}_1(x, y) = \begin{pmatrix} 20(x+1)^2 x^3 y(y-1)(2y-1) \\ 10(x+1)^2 x^2 y^2 (y-1)^2 (5x+3) \end{pmatrix}, \quad \mathbf{u}_2(x, y) = \begin{pmatrix} 10(y-1)^2 x^2 y^2 (x-1)^2 (5y-3) \\ 20(y-1)^2 y^3 x(x-1)(1-2x) \end{pmatrix}, \quad (70) $$

$$ p_1(x, y) = x + 1, $$

$$ p_2(x, y) = y, $$
<figure><img src="image_11.png" /><figcaption>FIGURE 5. Convergence histories of the velocity, pressure, Lagrange multiplier and mass loss in $\Omega_1$ obtained with the present method for $\theta = 1$ and different ratios between $h$ and $\mathcal{H}$.</figcaption></figure>

where the pressure profile has been chosen to provide a non-constant jump through the interface. The right-hand side functions $f$, $g$ and $u_\Sigma$ are chosen so that (1) holds with (70). In particular, zero velocity $u_\Sigma = 0$ is prescribed on the interface. The expression for the Lagrange multiplier $\lambda$ is obtained by plugging the exact solution (70) into (4). The fluid dynamic viscosity is set to $\mu = 0.035$.

In order to numerically assess the accuracy of (10), we consider a family of triangulations $T_h^\Omega$ resulting from the uniform subdivision of $\Gamma_D$ into $2N+1$ segments and of $\Gamma_N$ into $N+2$ segments, with $N \in \{5, 10, 20, 40, 80, 160\}$. The unfitted interfacial mesh $T_h^\Sigma$ is obtained by uniformly subdividing the domain $\Sigma$ into $3N$ segments. We hence have $\mathcal{H} \approx \frac{1}{3}h$. The associated convergence histories for the $L^2$-error of the velocity strain tensor, pressure and Lagrange multiplier are reported in Figure 4. As predicted by Theorem 4.12, a convergence rate of $1/2$ is observed for the approximation errors of both the velocity strain tensor and the pressure, as illustrated in Figures 4a and 4b, respectively. Although Theorem 4.12 does not guarantee convergence of the Lagrange multiplier in $L^2$-norm (see Remark 4.13), a convergence rate of $1/2$ is observed in Figure 4c. We can also notice the superior accuracy of the two variants of (10) with respect the unenriched FD-LM method. In particular, we can observe in Figure 4 the effect of the smaller constant in the error bound, as discussed in Section 4. The
enhanced mass conservation provided by the present method with respect to the FD-LM method (see Remark 3.2) is clearly visible in Figure 4d.

Finally, Figure 5 presents the relative $L^2$-norm errors for the velocity strain tensor, pressure, and Lagrange multiplier for the $\theta = 1$ method across various $h/\mathcal{H}$ ratios. The superior regularity of $\lambda$, with respect to the one assumed in Theorem 4.12, delivers a better convergence rate in terms of $\mathcal{H}$. Indeed, in this case, the estimate (60) scales as $\mathcal{H} + h^{\frac{1}{2}-\epsilon}$, so that we can take a coarser Lagrange multiplier mesh, $\mathcal{H} = h^{\frac{1}{2}}$, while preserving the overall convergence rate.

## 5.2. A constant pressure jump through a curved interface

As second example we consider the classical benchmark of two fluids separated by a curved interface (see Figure 3b) under a constant pressure drop (see, e.g., [35, 39]). This yields a very simple analytical solution: the two fluids are at rest with a constant pressure (in each sub-domain). This problem is critical for an unfitted mesh numerical method in terms of mass conservation, since continuous pressure approximations are known to generate spurious (non-zero) velocity approximations (which are proportional to the amount of pressure drop across the interface). We first study the convergence of the present method and then compare its behavior on a fixed mesh with respect to alternative methods.

In this example, the interface is given by $\Sigma = \{(0.2 \sin(\pi t), t) | t \in [0, 1]\}$. A traction force $g = -3 \cdot 10^5 \mathbf{n}$ is imposed at the channel inlet $\{-1\} \times (0, 1)$, while no traction $g = 0$ is prescribed at the outlet $\{1\} \times (0, 1)$. Additionally, we consider a null right-hand side function $f = 0$ for the Stokes equation (1) and the velocity profile on the immersed interface $\Sigma$ is set to be $u_\Sigma = 0$. The fluid dynamic viscosity is set to $\mu = 10$. The solution to problem (1) reads

$$ \mathbf{u} = 0, \quad p|_{\Omega_1} = 3 \cdot 10^5, \quad p|_{\Omega_2} = 0, \quad \lambda = -[\![p]\!]n_\Sigma. \qquad (71) $$

### 5.2.1. Convergence study

<figure><img src="image_12.png" /><figcaption>FIGURE 6. Geometric approximation of the normal vectors to the interface $\Sigma$. In gray the interface $\Sigma$ and in black the mesh $\mathcal{T}_h^\Sigma$.</figcaption></figure>

Since the interface is curved, we consider different discrete representations for the normal, see Figure 6. In practice, the choice of the normal approximation has an impact on the stabilization term $s_h^{BH}$, as well as the incompressibility constraint for the additional basis function (see Remark 3.2). The interface is discretized using a mesh $\mathcal{T}_h^\Sigma$, which is constituted of points linked by segments. The simplest way of representing the normal is then to define it by a piecewise constant function that takes the value of the exact normal to $\mathcal{T}_h^\Sigma$. This piecewise constant normal will be denoted by $n_\Sigma^0$. We also consider a $\mathbb{P}^1$-approximation, denoted $n_\Sigma^1$, as a weighted average of $n_\Sigma^0$ with respect to the length of the surrounding elements.
<figure><img src="image_13.png" /><figcaption>(a) Velocity energy error.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) Lagrange multiplier relative approximation error.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Pressure relative approximation error.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) Mass loss in $\Omega_1$. The results for $\theta = 0$ are omitted since this variant guarantees interfacial mass conservation.</figcaption></figure>

<caption>FIGURE 7. Convergence histories with the $\mathbb{P}^0$ and $\mathbb{P}^1$ approximation of the interface normal are compared among the present method in both the $\theta = 1$ and $\theta = 0$ variants and the un-enriched FD-LM method.</caption>

For the present method ($\theta = 0$ and $\theta = 1$) and the FD-LM method, the errors obtained with the two normal representations for the velocity, the pressure and the Lagrange multiplier are reported in Figure 7, using the same triangulations as in Section 5.1. The error of the normal discrete representation is given in Figure 8. For $\theta = 0$ or $\theta = 1$, we observe improved convergence rates for all the variables compared to what was obtained in Section 5.1. We also note the superior accuracy with respect to the stands FD-LM method. For the velocity and the pressure, we reach a convergence rate of 2 if a piecewise constant representation $n_{\Sigma}^0$ of the normal is used, while we reach 3 if $n_{\Sigma}^1$ is used instead. This is a consequence of the fact that, assuming the geometry is exactly represented, the exact solution $(u, p, \lambda)$ belongs to the discrete space $V_h \times Q_h \times \Lambda_H$. The error observed is then due to geometry approximation only. More precisely, taking into account the discrete normal representation $n_{\Sigma}^i$
<figure><img src="image_17.png" /><figcaption>FIGURE 8. Convergence histories for the $\mathbb{P}^0$ and $\mathbb{P}^1$ interface normals.</figcaption></figure>

$(i \in \{0, 1\})$ and the actual value of the exact solution, estimate (69) can be modified as

$$ \mu^{\frac{1}{2}} \|\varepsilon(u - u_h)\|_{L^2(\Omega)} + \mu^{-\frac{1}{2}} \|p - p_h\|_{L^2(\Omega)} + \mu^{-\frac{1}{2}} h^{\frac{1}{2}} \|\lambda - \lambda_H\|_{L^2(\Sigma)} \\ \lesssim \mu^{-\frac{1}{2}} \|[p]\|_{L^2(\Sigma)} (\|n_{\Sigma} - \Pi_H n_{\Sigma}\|_{L^2(\Sigma)} + \|n_{\Sigma}^i - \Pi_H n_{\Sigma}\|_{L^2(\Sigma)}). $$

This shows that the approximation errors are bounded by the geometric error on the normals. We observe this effect in Figure 7 where the convergence rate of the solution is linked to the convergence rate of the normal representation geometric error, see Figure 8. We also observe that a $\mathbb{P}_1$-representation for the normal is more efficient. On the contrary, for the FD-LM method the exact solution does not belong to the approximation space, so that we retrieve the classical $h^{\frac{1}{2}}$ convergence rate.

In addition, we have also considered test case (71) with a straight interface $\Sigma = \{0\} \times (0, 1)$ (see Figure 3a). In this geometrical setting, $n_{\Sigma} = n_{\Sigma}^1 = n_{\Sigma}^0$, thus the geometry of the interface is exactly represented by the mesh $\mathcal{T}_{\mathcal{H}}^{\Sigma}$ and $n_{\Sigma}$ belongs to $\Lambda_{\mathcal{H}}$. As expected by Theorem 4.12 the numerical solution corresponds to the analytical solution up to machine error.

### 5.2.2. Comparison with other fictitious domain methods

We now compare the behavior of the proposed unfitted mesh method (10) with two alternative fictitious domain approaches. In these two methods the discrete velocity and pressure are based on continuous piecewise affine functions (i.e., the spaces $V_h$ and $\tilde{Q}_h$ of Section 3) and the interface Dirichlet constraint is treated in a collocated-penalized fashion (see, e.g., [13, Section 3.2.1]). This amounts to consider in (2) the following non-conforming approximation of the Lagrange multiplier space $\Lambda$:

$$ \Lambda_{\mathcal{H}}^{\mathrm{C}} = \left\{ \xi_{\mathcal{H}} = \sum_{i=1}^{N_{\mathcal{H}}} \xi_i \delta_{x_i^s} \mid \xi_i \in \mathbb{R}^2, \quad i = 1, \dots, N_{\mathcal{H}} \right\}, $$

where $\{x_i^s\}_{i=1}^{N_{\mathcal{H}}} \text{ are the points of the interface mesh } \mathcal{T}_{\mathcal{H}}^{\Sigma} \text{ and } \delta_{x_i^s} \text{ denotes the Dirac's measure at point } x_i^s \text{ (see, e.g., [20, 26, 33])}$. The discrete treatment of the Dirichlet constraint takes the penalized-collocated form

$$ u_h(x_i^s) = u_{\Sigma}(x_i^s) + \varepsilon \lambda_i, \quad \forall i = 1, \dots, N_{\mathcal{H}}, \qquad (72) $$

where $\varepsilon > 0$ is a given (non-dimensionless) user-defined parameter. Note that the relation (72) enables the elimination of the Lagrange multiplier $\lambda$, so that the first alternative fictitious domain method reads: Find
<figure><img src="image_18.png" /><figcaption>FIGURE 9. Snapshots of the velocity magnitude.</figcaption></figure>

($\mathbf{u}_h, \tilde{\mathbf{p}}_h) \in \mathbf{V}_h \times \tilde{\mathbf{Q}}_h$ such that

$$a_{h,\varepsilon}^{\text{FD}}((\mathbf{u}_h, \tilde{\mathbf{p}}_h), (\mathbf{v}_h, \tilde{\mathbf{q}}_h)) = \ell_{\varepsilon}(\mathbf{v}_h) \quad (73)$$

for all $(\mathbf{v}_h, \tilde{\mathbf{q}}_h) \in \mathbf{V}_h \times \tilde{\mathbf{Q}}_h$, with the notations

$$a_{h,\varepsilon}^{\text{FD}}((\mathbf{u}_h, \tilde{\mathbf{p}}_h), (\mathbf{v}_h, \tilde{\mathbf{q}}_h)) \stackrel{\text{def}}{=} a((\mathbf{u}_h, \tilde{\mathbf{p}}_h), (\mathbf{v}_h, \tilde{\mathbf{q}}_h)) + s_h^{\text{BP}}(\tilde{\mathbf{p}}_h, \tilde{\mathbf{q}}_h) + d_{\varepsilon}(\mathbf{u}_h, \mathbf{v}_h),$$

$$\ell_{\varepsilon}(\mathbf{v}_h) \stackrel{\text{def}}{=} \ell(\mathbf{v}_h) + d_{\varepsilon}(\mathbf{u}_{\Sigma}, \mathbf{v}_h),$$

In the above definitions, the bilinear form $a$ and the linear functional $\ell$ are those defining the weak formulation of the Stokes problem (1), in Section 2. The results of the method (73) have been obtained with $\varepsilon = 10^{-5}$ and $\gamma_p = 10^{-2}$. This method will be denoted by FD in the sequel.

The second considered approach is a variant of (73) with an enhanced interfacial mass conservation. This is achieved by modifying the pressure and grad-div stabilizations in the neighborhood $\mathcal{T}_h^e$ of the interface $\Sigma$ as follows (see, e.g., [17, 30, 39]):

$$s^{\text{BP-GD}}((\mathbf{u}_h, \tilde{\mathbf{p}}_h), (\mathbf{v}_h, \tilde{\mathbf{q}}_h)) \stackrel{\text{def}}{=} \sum_{K \in \mathcal{T}_h^{\Omega}} \int_K \frac{\gamma_p h^2}{\mu} \nabla \tilde{\mathbf{p}}_h \cdot \nabla \tilde{\mathbf{q}}_h + \sum_{K \in \mathcal{T}_h^{\Omega}} \int_K \frac{\gamma_d \mu}{\delta_h} \nabla \cdot \mathbf{u}_h \nabla \cdot \mathbf{v}_h,$$

$$\delta_h \stackrel{\text{def}}{=} \begin{cases} 1 & K \in \mathcal{T}_h^{\Omega} \setminus \mathcal{T}_h^e, \\ \gamma_{\text{gd}} & K \in \mathcal{T}_h^e, \end{cases}$$
<figure><img src="image_19.png" /><figcaption>FIGURE 10. Snapshots of the pressure.</figcaption></figure>

where $\gamma_d > 0$ and $0 < \gamma_{gd} \ll 1$ are user-defined parameters. The fundamental idea of the method hence consists in boosting the grad-div stabilization (which enhances mass conservation) while reducing the pressure stabilization near the interface (which limits the impact of the interfacial pressure gradient inconsistencies). The resulting discrete problem reads as follows: Find $(\mathbf{u}_h, \tilde{\mathbf{p}}_h) \in \mathbf{V}_h \times \tilde{\mathbf{Q}}_h$ such that

$$a_{h,\varepsilon}^{\text{FD-GD}}((\mathbf{u}_h, \tilde{\mathbf{p}}_h), (\mathbf{v}_h, \tilde{\mathbf{q}}_h)) = \ell_{\varepsilon}(\mathbf{v}_h) \qquad (74)$$

for all $(\mathbf{v}_h, \tilde{\mathbf{q}}_h) \in \mathbf{V}_h \times \tilde{\mathbf{Q}}_h$, where

$$a_{h,\varepsilon}^{\text{FD-GD}}((\mathbf{u}_h, \tilde{\mathbf{p}}_h), (\mathbf{v}_h, \tilde{\mathbf{q}}_h)) \stackrel{\text{def}}{=} a((\mathbf{u}_h, \tilde{\mathbf{p}}_h), (\mathbf{v}_h, \tilde{\mathbf{q}}_h)) + s_h^{\text{BP-GD}}((\mathbf{u}_h, \tilde{\mathbf{p}}_h), (\mathbf{v}_h, \tilde{\mathbf{q}}_h)) + d_{\varepsilon}(\mathbf{u}_h, \mathbf{v}_h).$$

The results of the method (74) have been obtained with $\gamma_d = 1$, $\gamma_{gd} = 10^{-4}$, $\varepsilon = 10^{-5}$ and $\gamma_p = 10^{-2}$. We will refer to this method as FD-GD in the sequel.

<table><caption>TABLE 1. Mass loss in $\Omega_1$.</caption><thead><tr><th rowspan="2">FD<sub>ε</sub></th><th rowspan="2">FD - GD<sub>ε</sub></th><th rowspan="2">FD-LM</th><th colspan="2">Present method</th></tr><tr><th>θ = 1</th><th>θ = 0</th></tr></thead><tbody><tr><td>13</td><td>5.4 · 10<sup>-4</sup></td><td>17</td><td>3.4 · 10<sup>-10</sup></td><td>2.0 · 10<sup>-14</sup></td></tr></tbody></table>

Figures 9 and 10 report snapshots of the velocity magnitude and pressure obtained with the three considered fictitious domain methods, using the meshes of the previous convergence tests for $N = 40$. We recall that the exact solution is given by (71).
Figures 9a and 9c show the poor robustness that characterizes the FD method and the FD-LM method, respectively. Though the interface Dirichlet condition $u_Σ = 0$ is correctly imposed, the lack of additional mass conservation constraints results in a significant mass loss (see Table 1) which appears as a spurious velocity field in the whole computational domain $Ω$.

This issue is overcome in both the FD-GD method (Figure 9b) and the method proposed in this article (Figures 9d-e). These last two methods provide an adequate solution for the velocity and the pressure for this test case. As shown in Figure 10, due to the enrichment of the finite element space with an Heaviside function, the present method is, among the four, the only one capable of capturing the sharp pressure jump at the interface (see Figures 10d-e). Nevertheless, a good approximation of the pressure profile is also obtained for the FD, FD-GD and FD-LM methods (see Figures 10a-b-c).

### 5.3. Open interface

The last example is motivated by applications related to fluid-structure interaction problems involving topological changes, such as in cardiac valve simulations which involve both closed and open interfacial configurations. In the case of an open configuration, the separation between the two sub-domains $Ω_1$ and $Ω_2$ is not clear anymore. We propose to cope with this difficulty by adding a fictitious interface $Σ_{fic}$ closing $Σ$, as shown in Figure 3c.

<figure><img src="image_20.png" /><figcaption>FIGURE 11. Snapshots of the velocity magnitude.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>FIGURE 12. Snapshots of the velocity.</figcaption></figure>

In this example, the interface is made of two distinct parts, which are attached to the channel boundary but are disconnected in the middle, as illustrated in Figure 3c. More precisely, we set

$$ \Sigma = \{(0.2(1 - \cos(2\pi t)) - 0.1, t) | t \in [0, 0.4] \cup [0.6, 1]\}. $$

The Neumann boundary is given by $\Gamma_N = \{1\} \times (0, 1)$ and a Poiseuille velocity profile $u = -5y(1-y)n$ is prescribed on $\{-1\} \times (0, 1)$. Additionally, a traction-free boundary condition $g = 0$ is imposed on $\{1\} \times (0, 1)$. We also consider a null right-hand side function $f = 0$ for the Stokes equation (1) and a velocity profile $u_\Sigma = 0$ on the immersed interface $\Sigma$. The fluid dynamic viscosity is set to $\mu = 0.035$.

As mentioned above, in order to facilitate the definition of the discrete pressure space (6) in the present geometrical setting, we introduce a fictitious interface $\Sigma_{fic}$ that connects the two endpoints of $\Sigma$ immersed in $\Omega$. The sub-domain $\Omega_1$ (resp. $\Omega_2$) is then defined as the left (resp. right) side of the closed interface $\Sigma \cup \Sigma_{fic}$.
<figure><img src="image_22.png" /><figcaption>FIGURE 13. Snapshots of the pressure.</figcaption></figure>

In the present test case, we consider the straight line

$$ \Sigma_{\text{fic}} = \{(0.2(1 - \cos(0.8\pi)) - 0.1, y) | y \in [0.4, 0.6]\}, $$

as shown in Figure 3c.

It should be noted that, since no unknowns are attached to $\Sigma_{\text{fic}}$, this fictitious interface can be arbitrarily chosen, depending on the specific geometry of the problem, and its purpose is solely to facilitate the definition of $\Omega_1$. Nevertheless, from a computational perspective, when replacing the mass conservation constraint on $\Omega_1$ with its equivalent on $\partial\Omega_1$, as highlighted in Remark 3.2, it is essential to choose an appropriate discretization of $\Sigma_{\text{fic}}$ that enables the constraint to be accurately imposed. In this regard, a straight fictitious interface is a good choice.

In the sequel, we compare the results obtained using the method introduced in Section 3, with the two fictitious domain methods, FD and FD-GD, discussed in Section 5.2.2 and the standard FD-LM method without additional velocity constraint. As a reference solution, we consider the numerical approximation obtained through a fitted mesh approach using a highly refined mesh ($h \approx 6.25 \times 10^{-3}$). The results for the fictitious domain methods, are obtained using meshes as defined in Section 5.1 for $N = 40$. Specifically, the interface $\Sigma$ is discretized with 120 segments (i.e., 60 on each side). The fictitious interface $\Sigma_{\text{fic}}$ is discretized with 25 segments.
The value of the user defined parameter $\epsilon$ is set to $10^{-5}$ and $\gamma_p = 10^{-2}$. Furthermore, for the FD-GD method, we have $\gamma_d = 1$ and $\gamma_{gd} = 10^{-4}$. For the FD-LM method, the stabilization parameter $\gamma_\lambda$ is set to $10^2$. In the present method, the normal is represented by $n_\Sigma^1$ defined in Section 5.2.1 at the discrete level. Similar results are obtained when considering $n_\Sigma^0$.

Figures 11–13 display the velocity and pressure fields obtained with the four methods. For comparison purposes, the reference solution obtained with the fitted mesh method is also reported. Figures 11b-d show that the FD and FD-LM methods exhibit significant lack of mass conservation across the interface. This is particularly noticeable in the velocity profile near the interface, by comparing Figures 12a and 12b-d. In the FD and FD-LM methods the fluid passes through the interface, resulting in a lower velocity of the jet and the absence of a vortex behind the valve. On the contrary, the FD-GD method and the present method with $\theta = 1$ are able to capture the dynamics of the reference solution, as it can be inferred from Figures 12c-e. Figure 13c shows however that the grad-div penalty term induces an important perturbation of the pressure field in the FD-GD method. The pressure provided by the present method in $\Omega_1$ is slightly lower than the reference solution (compare Figures 13a and 13e). For the non-symmetric variant ($\theta = 0$), despite overall mass conservation in $\Omega_1$ (see Table 2), Figures 11f and 12f reveal fluid leakage through the interface and a significantly lower pressure in $\Omega_1$, as shown in Figure 13f. Taking $\gamma_\lambda$ larger in the non-symmetric variant ($\theta = 0$) slightly improves the velocity results, but at the price of degrading the approximation of $\lambda$. A slight pressure discontinuity is noticeable across the fictitious interface $\Sigma_{fic}$, which is compensated by the continuous part of the pressure. Indeed, as discussed in Remark 3.2, the amplitude $\|p_h\|$ of the additional Heaviside function of $Q_h$ can be interpreted as the Lagrange multiplier enforcing the mass conservation constraint in $\Omega_1$ and also as the mean of the pressure jump through the interface, which must not be zero in order to prevent mass leakage at the interface.

<table><caption>TABLE 2. Mass loss in $\Omega_1$.</caption><thead><tr><th rowspan="2">FD<sub>&epsilon;</sub></th><th rowspan="2">FD - GD<sub>&epsilon;</sub></th><th rowspan="2">FD-LM</th><th colspan="2">Present method</th></tr><tr><th>&theta; = 1</th><th>&theta; = 0</th></tr></thead><tbody><tr><td>9.4 &middot; 10<sup>-2</sup></td><td>1.8 &middot; 10<sup>-3</sup></td><td>9.8 &middot; 10<sup>-2</sup></td><td>2.2 &middot; 10<sup>-2</sup></td><td>2.2 &middot; 10<sup>-16</sup></td></tr></tbody></table>

## 6. CONCLUSION

In this paper, we have proposed and analyzed a new fictitious domain method for a Stokes problem with a Dirichlet constraint on an immersed interface. Although simplified, this setting appears in more complex fluid-solid interaction problems with immersed thin-walled solids. This work was motivated by numerical evidence showing that inaccuracy issues in standard fictitious domain methods for incompressible flows are mainly driven by the lack of mass conservation across the interface. We have shown that this difficulty can be mitigated by adding a single global velocity constraint enforced on one side of the interface via a scalar Lagrange multiplier. Another salient feature of the method is that the size of the system matrix is constant irrespectively of the location of the interface with respect to the background computational mesh, which is particularly appealing for the case of moving interfaces. A complete a priori numerical analysis of the method has been provided. The additional basis function leads to a priori error bounds that are robust with respect to the magnitude of the pressure jump. As expected, sub-optimal convergence rate of order $\frac{1}{2} - \epsilon$ is obtained for both the velocity and pressure, under minimal regularity assumptions. The numerical results of Section 5 show that the symmetric variant of the proposed method ($\theta = 1$) provides similar or superior accuracy to alternative fictitious domain methods, without the need of resorting to penalty terms which jeopardize the conditioning of the resulting system matrix. Another salient conclusion that can be drawn from the numerics is that there is an intriguing interplay between perturbation of the velocity constraint and control/consistency provided by the Lagrange multiplier stabilization. In particular, the results of the asymmetric variant ($\theta = 0$) show that the exact satisfaction of
the velocity constraint does not guarantee robustness. Current investigations are devoted to the application
of the present method to the simulation of fluid-structure interaction of heart valves.

REFERENCES

[1] F. Alauzet. A changing-topology moving mesh technique for large displacements. *Engineering with Computers*, 30(2):175–200, 2014.

[2] F. Alauzet, B. Fabrèges, M. A. Fernández, and M. Landajuela. Nitsche-XFEM for the coupling of an incompressible fluid with immersed thin-walled structures. *Computer Methods in Applied Mechanics and Engineering*, 301:300–335, 2016.

[3] M. Annese, M. A. Fernández, and L. Gastaldi. Splitting schemes for a Lagrange multiplier formulation of FSI with immersed thin-walled structure: stability and convergence analysis. *IMA Journal of Numerical Analysis*, March 2022.

[4] C. Atamian, Q. V. Dinh, R. Glowinski, J. He, and J. Périaux. Control approach to fictitious-domain methods. application to fluid dynamics and electro-magnetics. In *Domain decomposition methods for partial differential equations*, 1991.

[5] F. Baaijens. A fictitious domain/mortar element method for fluid-structure interaction. *International Journal for Numerical Methods in Fluids*, 35:743–761, 04 2001.

[6] H. J. C. Barbosa and T. J. R. Hughes. The finite element method with Lagrange multipliers on the boundary: circumventing the Babuška-Brezzi condition. *Computer Methods in Applied Mechanics and Engineering*, 85(1):109–128, 1991.

[7] G.R. Barrenechea and C. González. A stabilized finite element method for a fictitious domain problem allowing small inclusions. *Numer. Methods Partial Differential Equations*, 34(1):167–183, 2018.

[8] R. Becker, E. Burman, and P. Hansbo. A Nitsche extended finite element method for incompressible elasticity with discontinuous modulus of elasticity. *Computer Methods in Applied Mechanics and Engineering*, 198(41-44):3352–3360, 2009.

[9] F. Bertrand, P. A. Tanguy, and F. Thibault. A three-dimensional fictitious domain method for incompressible fluid flow problems. *International Journal for Numerical Methods in Fluids*, 25(6):719–736, 1997.

[10] D. Boffi, N. Cavallini, F. Gardini, and L. Gastaldi. Local mass conservation of Stokes finite elements. *Journal of Scientific Computing*, 52:383–400, 2012.

[11] D. Boffi, N. Cavallini, F. Gardini, and L. Gastaldi. Stabilized Stokes elements and local mass conservation. *Bollettino dell'Unione Matematica Italiana*, 5(3):543–573, 10 2012.

[12] D. Boffi and L. Gastaldi. A fictitious domain approach with distributed Lagrange multiplier for fluid-structure interactions. *Numerische Mathematik*, 135, 03 2017.

[13] L. Boilevin-Kayl, M. A. Fernández, and J.-F. Gerbeau. Numerical methods for immersed FSI with thin-walled structures. *Computers & Fluids*, 179:744–763, 2019.

[14] F. Brezzi and J. Pitkäranta. On the stabilization of finite element approximations of the Stokes equations. In W. Hackbush, editor, *Efficient Solutions of Elliptic Systems*, volume 10 of *Notes on Numerical Fluid Mechanics*. Vieweg, Braunschweig, 1984.

[15] E. Burman and P. Hansbo. Fictitious domain methods using cut elements: III. A stabilized Nitsche method for Stokes' problem. *European Series in Applied and Industrial Mathematics (ESAIM): Mathematical Modelling and Numerical Analysis*, 48, 05 2014.

[16] E. Burman, P. Hansbo, and M. G. Larson. Cut finite element method for divergence free approximation of incompressible flow: optimal error estimates and pressure independence, 2022. arXiv:2207.04734.

[17] H. Casquero, C. Bona-Casas, and H. Gomez. NURBS-based numerical proxies for red blood cells and circulating tumor cells in microscale blood flow. *Computer Methods in Applied Mechanics and Engineering*, 316:646–667, 2017.

[18] H. Casquero, C. Bona-Casas, D. Toshniwal, T.J.R. Hughes, H. Gomez, and Y. Jessica Zhang. The divergence-conforming immersed boundary method: Application to vesicle and capsule dynamics. *Journal of Computational Physics*, 425:109872, 2021.

[19] S. H. Christiansen and K. Hu. Generalized finite element systems for smooth differential forms and Stokes' problem. *Numerische Mathematik*, 140:327–371, 2016.

[20] N. D. Dos Santos, J.-F. Gerbeau, and J.-F. Bourgat. A partitioned fluid-structure algorithm for elastic thin valves with contact. *Computer Methods in Applied Mechanics and Engineering*, 197(19):1750–1761, 2008.

[21] Q. Du, M. D. Gunzburger, L. S. H. Hou, and J. Lee. Semidiscrete finite element approximations of a linear fluid-structure interaction problem. *SIAM Journal on Numerical Analysis*, 42:1–29, 2004.

[22] M. Duprez, V. Lleras, and A. Lozinski. φ-FEM: an optimally convergent and easily implementable immersed boundary method for particulate flows and Stokes equations, 2023. arXiv:2211.07012.

[23] M. Duprez and A. Lozinski. φ-FEM: a finite element method on domains defined by level-sets. *SIAM Journal on Numerical Analysis*, 58(2):1008–1028, 2020.

[24] Alexandre Ern and Jean-Luc Guermond. *Finite elements I—Approximation and interpolation*, volume 72 of *Texts in Applied Mathematics*. Springer, Cham, [2021] ©2021.

[25] B. Fabrèges. *A smooth extension method for the simulation of fluid/particles flows*. PhD thesis, Université Paris Sud - Paris XI, 2012. (In French).
[26] B. Fabrèges and B. Maury. Approximation of single layer distributions by dirac masses in finite element computations. *Journal of Scientific Computing*, 58(1):25-40, 2014.

[27] L. Formaggia, J.-F. Gerbeau, F. Nobile, and A. Quarteroni. Numerical treatment of defective boundary conditions for the Navier-Stokes equations. *SIAM J. Numer. Anal.*, 40(1):376-401, 2002.

[28] M. Fournié and A. Lozinski. Stability and optimal convergence of unfitted extended finite element methods with Lagrange multipliers for the Stokes equations. In Stéphane P. A. Bordas, Erik Burman, Mats G. Larson, and Maxim A. Olshanskii, editors, *Geometrically Unfitted Finite Element Methods and Applications*, pages 143-182, Cham, 2017. Springer International Publishing.

[29] L. P. Franca and R. Stenberg. Error analysis of some Galerkin least squares methods for the elasticity equations. *SIAM Journal on Numerical Analysis*, 28(6):1680-1697, 1991.

[30] K.J. Galvin, A. Linke, L.G. Rebholz, and N.E. Wilson. Stabilizing poor mass conservation in incompressible flow problems with large irrotational forcing and application to thermal convection. *Computer Methods in Applied Mechanics and Engineering*, 237/240:166-176, 2012.

[31] V. Girault and R. Glowinski. Error analysis of a fictitious domain method applied to a Dirichlet problem. *Japan Journal of Industrial and Applied Mathematics*, 12(3):487-514, 1995.

[32] V. Girault and P.A. Raviart. *Finite Element Methods for Navier-Stokes Equations*. Springer Series in Computational Mathematics, Vol. 5. Springer-Verlag, Berlin, 1986.

[33] 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.

[34] P. Grisvard. *Elliptic problems in nonsmooth domains*, volume 69 of *Classics in Applied Mathematics*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2011. Reprint of the 1985 original [MR0775683], With a foreword by Susanne C. Brenner.

[35] S. Groß and A. Reusken. An extended pressure finite element space for two-phase incompressible flows with surface tension. *Journal of Computational Physics*, 224(1):40-58, 2007.

[36] J. Guzmán and M. A. Olshanskii. Inf-sup stability of geometrically unfitted Stokes finite elements. *Mathematics of Computation*, 87:2091-2112, 2018.

[37] J. Haslinger and Y. Renard. A new fictitious domain approach inspired by the extended finite element method. *SIAM Journal on Numerical Analysis*, 47, 01 2009.

[38] T. Hisada and T. Washio. Mathematical considerations for fluid-structure interaction simulation of heart valves. *Bulletin of the Japan Society for Industrial and Applied Mathematics*, 16(2):142-156, 2006. (In Japanese).

[39] D. Kamensky, M.-C. Hsu, D. Schillinger, J. A. Evans, A. Aggarwal, Y. Bazilevs, M. S. Sacks, and T. J. R. Hughes. An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves. *Computer methods in applied mechanics and engineering*, 284:1005-1053, 2015.

[40] M. Kirchhart, S. Groß, and A. Reusken. Analysis of an XFEM discretization for Stokes interface problems. *SIAM Journal on Scientific Computing*, 38(2):A1019-A1043, 2016.

[41] K. Li, N. M. Atallah, G. A. Main, and G. Scovazzi. The shifted interface method: A flexible approach to embedded interface computations. *International Journal for Numerical Methods in Engineering*, 121(3):492-518, 2020.

[42] H. Liu, M. Neilan, and M. Olshanskii. A cutFEM divergence-free discretization for the Stokes problem. *ESAIM: M2AN*, 57:143-165, 2023.

[43] Y. Liu and W. K. Liu. Rheology of red blood cell aggregation by computer simulation. *Journal of Computational Physics*, 220(1):139-154, 2006.

[44] T. Nakata and H. Liu. A fluid-structure interaction model of insect flight with flexible wings. *Journal of Computational Physics*, 231(4):1822-1847, 2012.

[45] P. A. Nguyen and J.-P. Raymond. Boundary stabilization of the Navier-Stokes equations in the case of mixed boundary conditions. *SIAM Journal on Control and Optimization*, 53(5):3006-3039, 2015.

[46] K. Ohmori and N. Saito. Flux-free finite element method with Lagrange multipliers for two-fluid flows. *J. Sci. Comput.*, 32(2):147-173, 2007.

[47] C. S. Peskin. The immersed boundary method. *Acta Numerica*, 11:479-517, 2002.

[48] T. Sawada and A. Tezuka. LLM and X-FEM based interface modeling of fluid-thin structure interactions on a non-interface-fitted mesh. *Computational Mechanics*, 48(3):319-332, 2011.

[49] R. L. Scott and S. Zhang. Finite element interpolation of nonsmooth functions satisfying boundary conditions. *Mathematics of Computation*, 54(190):483-493, 1990.

[50] P. Singh, D.D. Joseph, T.I. Hesla, R. Glowinski, and T.-W. Pan. A distributed Lagrange multiplier/fictitious domain method for viscoelastic particulate flows. *Journal of Non-Newtonian Fluid Mechanics*, 91(2):165-188, 2000.

[51] K. Takizawa and T. E. Tezduyar. Computational methods for parachute fluid-structure interactions. *Archives of Computational Methods in Engineering*, 19(1):125-169, 2012.

[52] K. Takizawa, T.E. Tezduyar, A. Buscher, and S. Asada. Space-time interface-tracking with topology change (ST-TC). *Computational Mechanics*, 54(4):955-971, 2014.
[53] R. Van Loon, P. D. Anderson, F. P. T. Baaijens, and F. N. Van de Vosse. A three-dimensional fluid-structure interaction method for heart valve modelling. *Comptes Rendus Mécanique*, 333(12):856-866, 2005.

[54] F. Vergnet. *Active structures in viscous flows: models, mathematical analysis and numerical simulations*. PhD thesis, Université Paris-Saclay, 2019.

[55] G. D. Weymouth, D. G. Dommermuth, K. Hendrickson, and D.K.P. Yue. Advancements in Cartesian-grid methods for computational ship hydrodynamics. In *26th Symposium on Naval Hydrodynamics (16/09/06 - 21/09/06)*, Rome, Italy, September 2006.

[56] A. Zilian and A. Legay. The enriched space-time finite element method (EST) for simultaneous solution of fluid-structure interaction. *International Journal for Numerical Methods in Engineering*, 75(3):305-334, 2008.