Daniele C Corti
email: daniele.corti@inria.fr

Guillaume Delay
email: guillaume.delay@sorbonne-universite.fr

Miguel A Fernández
email: miguel.fernandez@inria.fr

Fabien Vergnet
email: fabien.vergnet@sorbonne-universite.fr

Marina Vidrascu
email: marina.vidrascu@inria.fr

M C A Corti

Fernández

Mathematical Modelling and Numerical Analysis Will be set by the publisher Modélisation Mathématique et Analyse Numérique LOW-ORDER FICTITIOUS DOMAIN METHOD WITH ENHANCED MASS CONSERVATION FOR AN INTERFACE STOKES PROBLEM *

Keywords: Stokes problem, fictitious domain method, Lagrange multiplier, interfacial mass conservation

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.

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., [START_REF] Liu | Rheology of red blood cell aggregation by computer simulation[END_REF][START_REF] Nakata | A fluid-structure interaction model of insect flight with flexible wings[END_REF][START_REF] Takizawa | Computational methods for parachute fluid-structure interactions[END_REF][START_REF] Van Loon | A three-dimensional fluid-structure interaction method for heart valve modelling[END_REF][START_REF] Weymouth | Advancements in Cartesian-grid methods for computational ship hydrodynamics[END_REF]). 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.

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., [START_REF] Alauzet | A changing-topology moving mesh technique for large displacements[END_REF][START_REF] Takizawa | Space-time interface-tracking with topology change (ST-TC)[END_REF]) 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., [START_REF] Bertrand | A three-dimensional fictitious domain method for incompressible fluid flow problems[END_REF][START_REF] Boffi | A fictitious domain approach with distributed Lagrange multiplier for fluid-structure interactions[END_REF][START_REF] Glowinski | A distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF][START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Singh | A distributed Lagrange multiplier/fictitious domain method for viscoelastic particulate flows[END_REF]) 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 [START_REF] Atamian | Control approach to fictitious-domain methods. application to fluid dynamics and electro-magnetics[END_REF][START_REF] Fabrèges | A smooth extension method for the simulation of fluid/particles flows[END_REF][START_REF] Vergnet | Active structures in viscous flows: models, mathematical analysis and numerical simulations[END_REF], 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., [START_REF] Burman | Fictitious domain methods using cut elements: III. A stabilized Nitsche method for Stokes' problem[END_REF][START_REF] Groß | An extended pressure finite element space for two-phase incompressible flows with surface tension[END_REF][START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF]), 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 ϕ-FEM method (see, e.g., [START_REF] Duprez | ϕ-FEM: an optimally convergent and easily implementable immersed boundary method for particulate flows and Stokes equations[END_REF][START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by level-sets[END_REF][START_REF] Li | The shifted interface method: A flexible approach to embedded interface computations[END_REF]) 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., [START_REF] Boilevin-Kayl | Numerical methods for immersed FSI with thin-walled structures[END_REF][START_REF] Casquero | NURBS-based numerical proxies for red blood cells and circulating tumor cells in microscale blood flow[END_REF][START_REF] Galvin | Stabilizing poor mass conservation in incompressible flow problems with large irrotational forcing and application to thermal convection[END_REF][START_REF] Kamensky | An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves[END_REF][START_REF] Liu | A cutFEM divergence-free discretization for the Stokes problem[END_REF]), 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., [START_REF] Alauzet | Nitsche-XFEM for the coupling of an incompressible fluid with immersed thin-walled structures[END_REF][START_REF] Becker | A Nitsche extended finite element method for incompressible elasticity with discontinuous modulus of elasticity[END_REF][START_REF] Kirchhart | Analysis of an XFEM discretization for Stokes interface problems[END_REF][START_REF] Sawada | LLM and X-FEM based interface modeling of fluid-thin structure interactions on a non-interfacefitted mesh[END_REF][START_REF] Zilian | The enriched space-time finite element method (EST) for simultaneous solution of fluid-structure interaction[END_REF]), 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., [START_REF] Baaijens | A fictitious domain/mortar element method for fluid-structure interaction[END_REF][START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF][START_REF] Boffi | Stabilized Stokes elements and local mass conservation[END_REF]). Similar observations can be made on the combination of unfitted mesh methods and divergence free approximations (see, e.g., [START_REF] Casquero | The divergence-conforming immersed boundary method: Application to vesicle and capsule dynamics[END_REF][START_REF] Liu | A cutFEM divergence-free discretization for the Stokes problem[END_REF]), with the exception of the cut-FEM method reported in [START_REF] Burman | Cut finite element method for divergence free approximation of incompressible flow: optimal error estimates and pressure independence[END_REF] based on the minimal divergence free element introduced in [START_REF] Christiansen | Generalized finite element systems for smooth differential forms and Stokes' problem[END_REF].

Numerical evidence (see, e.g., [START_REF] Boilevin-Kayl | Numerical methods for immersed FSI with thin-walled structures[END_REF]) 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 [START_REF] Hisada | Mathematical considerations for fluid-structure interaction simulation of heart valves[END_REF], 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 [START_REF] Formaggia | Numerical treatment of defective boundary conditions for the Navier-Stokes equations[END_REF][START_REF] Ohmori | Flux-free finite element method with Lagrange multipliers for two-fluid flows[END_REF], 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

• • • • D N D N ⌦ 1 ⌦ 2 ⌃ • • • • • • • • • • • • D N D N ⌦ 1 ⌦ 2 ⌃ Figure 1
. Left: admissible geometric configuration. Right: non-admissible configuration (see Assumption 2.1 and Remark 2.2).

will be extensively used in the paper. Let be a scalar or tensorial field defined in ⌦ and possibly discontinuous across ⌃, we define its sided restrictions, denoted by 1 and 2 , as

1 (x) def = lim !0 (x + n 1 ) , 2 (x) def = lim !0 (x + n 2 ) , 8x 2 ⌃,
and the following jump operators across the interface ⌃:

J K def = 1 2 , J nK def = 1 n 1 + 2 n 2 = J Kn ⌃ .
TITLE WILL BE SET BY THE PUBLISHER 3 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.

Problem setting

Let ⌦ ⇢ R d , d = 2 or 3, be a bounded polygonal/polyhedral domain and ⌃ be a polygonal oriented manifold of co-dimension 1 fully immersed in ⌦ and with unit normal vector n ⌃ . We assume that ⌃ divides ⌦ into two sub-domains ⌦ 1 and ⌦ 2 such that ⌦ 1 \ ⌦ 2 = ;, @⌦ 1 \ @⌦ 2 = ⌃ and @⌦ = (@⌦ 1 [ @⌦ 2 ) \ ⌃, where @⌦ i and @⌦ stand for the boundaries of ⌦ i and ⌦, respectively. Moreover, @⌦ is partitioned as @⌦ = D [ N with D (respectively 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 ⌦ i . Note that ⌦ 1 and ⌦ 2 are arbitrarily chosen so that n 1 = n ⌃ and n 2 = n ⌃ on ⌃. We also denote by n the outward unit normal vector to @⌦.

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

Assumption 2.1. The Neumann boundary N is away from the interface ⌃; i.e. N \ ⌃ = ;. Moreover, both sub-domains have a Neumann boundary; i.e. 8i 2 {1, 2}, @⌦ i \ N 6 = ;.

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

• • • • D N D N ⌦ 1 ⌦ 2 ⌃ • • • • • • • • • • • • D N D N ⌦ 1 ⌦ 2 ⌃ Figure 1
. Left: admissible geometric configuration. Right: non-admissible configuration (see Assumption 2.1 and Remark 2.2).

will be extensively used in the paper. Let be a scalar or tensorial field defined in ⌦ and possibly discontinuous across ⌃, we define its sided restrictions, denoted by 1 and 2 , as

1 (x) def = lim !0 (x + n 1 ) , 2 (x) def = lim !0 (x + n 2 ) , 8x 2 ⌃,
and the following jump operators across the interface ⌃: will be extensively used in the paper. Let ϕ be a scalar or tensorial field defined in Ω and possibly discontinuous across Σ, we define its sided restrictions, denoted by ϕ 1 and ϕ 2 , as

J K def = 1 2 , J nK def = 1 n 1 + 2 n 2 = J Kn ⌃ .
ϕ 1 (x) def = lim δ→0 - ϕ (x + δn 1 ) , ϕ 2 (x) def = lim δ→0 - ϕ (x + δn 2 ) , ∀x ∈ Σ,
and the following jump operators across the interface Σ:

ϕ def = ϕ 1 -ϕ 2 , ϕn def = ϕ 1 n 1 + ϕ 2 n 2 = ϕ n Σ .
For any subset ω ⊂ Ω, we will use the standard Sobolev spaces H m (ω) (m ≥ 0), equipped with the standard norm ∥•∥ H m (ω) , while (•, •) ω denotes the usual L 2 (ω)-scalar product. We also denote by |ω| the usual Lebesgue measure of ω. For S ⊂ ∂Ω ∪ Σ, the closed subspace H 1 S represents the space of H 1 functions with zero value on S.

Let

f ∈ [L 2 (Ω)] d , g ∈ [L 2 (Γ N )] d and u Σ ∈ [H 1 0 (Σ)
] d be given fields. We consider the following Stokes problem in Ω with an immersed Dirichlet condition: Find u : Ω → R d and p : Ω → R, fluid velocity and pressure respectively, such that:

               -∇ • σ(u, p) = f in Ω i , ∇ • u = 0 in Ω i , u = 0 on Γ D , σ(u, p)n = g on Γ N , u = u Σ on Σ, (1)
for i = 1, 2. Here, the symbol σ(u, p) stands for the Cauchy stress tensor of the fluid, given by

σ(u, p) def = 2µε(u) -pI, ε(u) def = 1 2 ∇u + ∇u T ,
where µ > 0 denotes the fluid dynamic viscosity.

Note that system (1) corresponds to two independent Stokes problems with a shared Dirichlet condition on the interface Σ. This kind of problems arises, for instance, in fluid-structure interaction problems with immersed thin-walled solids (see, e.g., [START_REF] Boilevin-Kayl | Numerical methods for immersed FSI with thin-walled structures[END_REF]). 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 Σ. To this purpose, we first introduce a global variational formulation of (1) in which Lagrange multipliers are used to enforce the Dirichlet constraint on Σ.

We consider the following velocity and pressure functional spaces

V def = [H 1 ΓD (Ω)] d , Q def = L 2 (Ω),
respectively. We make use of the standard Stokes bilinear form a :

(V × Q) × (V × Q) → R and the linear form ℓ : V → R , given by a((u, p), (v, q)) def = 2µ (ε (u) , ε (v)) Ω -(p, ∇ • v) Ω + (q, ∇ • u) Ω , ℓ (v) def = (f , v) Ω + (g, v) ΓN .
For the weak treatment of the Dirichlet interface condition on Σ, we consider the Lagrange multiplier

Λ def = [H 1 2 00 (Σ)] d ′ ,
(u, p, λ) ∈ V × Q × Λ such that A( (u, p, λ) , (v, q, ξ) ) = F (v, ξ) , (2)
for all (v, q, ξ) ∈ V × Q × Λ, where the operators A and F are defined as follows:

A( (u, p, λ) , (v, q, ξ)

) def = a( (u, p) , (v, q) ) -b (λ, v) + b (ξ, u) , F (v, ξ) def = ℓ (v) + b(ξ, u Σ ). (3)
Using a standard argument of integration by parts in Ω i , it can be shown that problem (1) is equivalent to problem (2) with λ = σ(u, p)n on Σ.

In other words, the unknown Lagrange multiplier represents the jump of the fluid stress across the interface Σ resulting from the enforcement of the Dirichlet condition on Σ. 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 [START_REF] Alauzet | A changing-topology moving mesh technique for large displacements[END_REF], the potential presence of discontinuities on the interface Σ for the solution (resp. its gradient) reduces the regularity of p (resp. u). More precisely, instead of the standard spaces H 1 (Ω) and [H 2 (Ω)] d for the pressure and velocity, we use the spaces H 1/2-ε (Ω) and [H 3/2-ε (Ω)] d , for some ε ∈ (0, 1 2 ). Furthermore, the regularity of p (resp. u) in each subdomain

Ω i , with i = 1, 2, is limited to H 1 2 +ε (Ω i ) (resp. [H 3 2 +ε (Ω i )] d ) due
to the potential presence of reentrant corners on the interface and corners between Γ D and Γ N (see Figure 1). For further details on these aspects the reader is referred to [34,[START_REF] Nguyen | Boundary stabilization of the Navier-Stokes equations in the case of mixed boundary conditions[END_REF]. Accordingly, we introduce the following spaces:

V ε = v ∈ [H 3 2 -ε ΓD (Ω)] d | v| Ω i ∈ [H 3 2 +ε (Ω i )] d , i = 1, 2 , Q ε = q ∈ H 1 2 -ε (Ω) | q| Ω i ∈ H 1 2 +ε (Ω i ) , i = 1, 2 .
(5) Note that the traces of elements of these spaces are well defined on each sub-domain boundary ∂Ω i .

Discrete problem

Let us introduce {T Ω

h } 0<h<1 and {T Σ H } 0<H<1 two families of quasi uniform simplicial meshes (see, e.g., [START_REF] Ern | Finite elements I-Approximation and interpolation[END_REF]) of Ω and Σ, respectively. The symbol h (respectively H) denotes the mesh parameter for Ω (respectively Σ). We assume that T Ω h is fitted to ∂Ω but, in general, not to Σ. We consider the two following standard spaces of continuous piecewise affine functions:

X Ω h def = χ h ∈ C 0 (Ω)| χ h | K ∈ P 1 (K), ∀K ∈ T Ω h , X Σ H def = χ H ∈ C 0 (Σ)| χ H | K ∈ P 1 (K), ∀K ∈ T Σ H .
Finally, we introduce the discrete space Q h def = X Ω h ⊂ Q and the discrete spaces V h , Q h , Λ H for the approximation of the velocity, pressure and Lagrange multiplier as follows:

V h def = [X Ω h ] d ∩ V , Q h def = Q h ⊕ ⟨1 Ω1 ⟩ ⊂ Q, Λ H def = [X Σ H ] d ⊂ Λ, (6)
where the symbol 1 Ω1 denotes the characteristic function of Ω 1 and ⟨1 Ω1 ⟩ the vector space spanned by 1 Ω1 . We recall that the choice of Ω 1 and Ω 2 is arbitrary, so that the characteristic function could have been taken over Ω 2 instead of over Ω 1 . The Lagrange multiplier space Λ H could alternatively be made of discontinuous functions (as in, e.g., [START_REF] Barrenechea | A stabilized finite element method for a fictitious domain problem allowing small inclusions[END_REF][START_REF] Girault | Error analysis of a fictitious domain method applied to a Dirichlet problem[END_REF]).

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 Ω, but does not prevent spurious leakage across the interface Σ (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 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., [START_REF] Groß | An extended pressure finite element space for two-phase incompressible flows with surface tension[END_REF]) 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 ∈ Q h can be decomposed into the sum of two contributions, namely,

q h = q h + q h , (7)
where q h ∈ Q h is a continuous piecewise affine function over the whole computational domain Ω and q h ∈ ⟨1 Ω1 ⟩ is a Heaviside function that allows strong discontinuities across the interface Σ. Moreover, the following relations hold

q h = q h 1 Ω1 , q h = q h -q h 1 Ω1 , (8)
where q h is the constant jump of q h across Σ.

In order to overcome the instability of the saddle-point problem due to the choice of the Lagrange multiplier finite element space Λ H in [START_REF] Barbosa | The finite element method with Lagrange multipliers on the boundary: circumventing the Babuška-Brezzi condition[END_REF] and to avoid any constraint on the ratio between the fluid mesh size h and the solid mesh size H (as can, for instance, be seen in [START_REF] Boffi | A fictitious domain approach with distributed Lagrange multiplier for fluid-structure interactions[END_REF]), we introduce the following stabilization term (in the spirit of the Barbosa-Hughes stabilization [START_REF] Barbosa | The finite element method with Lagrange multipliers on the boundary: circumventing the Babuška-Brezzi condition[END_REF]):

s BH h ( (p h , λ H ) , (q h , ξ H ) ) def = h γ λ µ (λ H + p h n Σ , ξ H + θ q h n Σ ) Σ , (9)
with γ λ > 0 a user-defined dimensionless parameter and θ ∈ {0, 1}. The choice θ = 1 yields a symmetric stabilization term, while for θ = 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 [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF] is considered to overcome the lack of inf-sup compatibility between the velocity and pressure spaces V h and Q h :

s BP h ( p h , q h ) def = γ p h 2 µ (∇ p h , ∇ q h ) Ω ,
with γ p > 0 a user-defined dimensionless parameter. Note that the pressure stabilization operator s BP h only acts on the continuous part, p h , of the discrete pressure,

p h ∈ Q h .
The proposed discrete approximation of (2) reads as follows:

Find (u h , p h , λ H ) ∈ V h × Q h × Λ H such that A h ( (u h , p h , λ H ) , (v h , q h , ξ H ) ) = F (v h , ξ H ) (10)
for all

(v h , q h , ξ H ) ∈ V h × Q h × Λ H , with the notation A h ( (u h , p h , λ H ) , (v h , q h , ξ H ) ) def = A( (u h , p h , λ H ) , (v h , q h , ξ H ) ) + s BP h ( p h , q h ) + s BH h ( (p h , λ H ) , (q h , ξ H ) ). ( 11
)
Remark 3.1. It is worth noticing that the stabilization operator ( 9) is only weakly consistent since, according to (4), the Lagrange multiplier λ is not necessarily equal top n Σ . In the spirit of [START_REF] Barbosa | The finite element method with Lagrange multipliers on the boundary: circumventing the Babuška-Brezzi condition[END_REF], a strongly consistent method would be

h γ λ µ (λ H -σ (u h , p h ) n , ξ H -θ σ (v h , q h ) n ) Σ .
The reason of excluding the normal component of the strain tensor ε (v h ) n Σ 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 Q h for the pressure, the stabilization operator (9) reduces to the one investigated in [START_REF] Barrenechea | A stabilized finite element method for a fictitious domain problem allowing small inclusions[END_REF] 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

(∇ • v h , q h ) Ω = (∇ • v h , q h ) Ω + q h (v h • n 1 , 1) ∂Ω1 = (∇ • v h , q h ) Ω + q h ∂Ω1 v h • n 1 for all (v h , q h ) ∈ V h × Q h .
Therefore, owing to ( 7)-(8), the system matrix associated to [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] 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,

(∇ • u h , q h ) Ω + s BP h ( p h , q h ) = 0 ∀ q h ∈ Q h , (12)
which yields global mass conservation in Ω, while the second enhances interfacial mass conservation by enforcing

∂Ω1 u h • n 1 + θ h γ λ µ Σ (λ H + p h n Σ ) • n Σ = 0. (13)
For θ = 0, this relation guaranties global mass conservation in Ω 1 (and also in Ω 2 , thanks to [START_REF] Boffi | A fictitious domain approach with distributed Lagrange multiplier for fluid-structure interactions[END_REF]). On the contrary, for θ = 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 [START_REF] Boilevin-Kayl | Numerical methods for immersed FSI with thin-walled structures[END_REF] with θ = 0 for enhancing interfacial mass conservation in a fictitious domain context was reported for the first time in [START_REF] Hisada | Mathematical considerations for fluid-structure interaction simulation of heart valves[END_REF], for a different method and without error analysis. Similar ideas have also been advocated in [START_REF] Formaggia | Numerical treatment of defective boundary conditions for the Navier-Stokes equations[END_REF][START_REF] Ohmori | Flux-free finite element method with Lagrange multipliers for two-fluid flows[END_REF] in a different physical framework.

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 ≲ B (resp. A ≳ B ) to abbreviate the inequality A ≤ CB (resp. CA ≥ B) for positive real numbers A and B, where the constant C > 0 does not depend on h, H, µ, neither on the way the interface mesh T Σ H cuts the fluid mesh T Ω h nor on the solutions of the continuous and discrete problems. We also use the notation H s (ω) def = [H s (ω)] d for all s > 0 and all ω ⊂ Ω. We now introduce some geometric definitions that are widely used in the current section. Let K be an element of T Ω h and x K,i its i-th node. We denote by T c

h def = {K ∈ T Ω h : K ∩ Σ ̸ = ∅} the set of elements of T Ω h cut by the interface Σ, by N h def = {x K,i : K ∈ T Ω h } the set of nodes, by N c h def = {x K,i : K ∈ T c h }
q def = ✓ 1 |⌃| Z ⌃ JqK ◆ 1 ⌦1 , (14)
the mean value of the pressure jump across the interface ⌃, and

e ⇧ h q def = ⇧h (q b ⇧ h q),
with ⇧h : H . The elements of T c h are marked in , the elements of T e h are marked in and the nodes in N c h are marked by .

Q ✏ ! e Q h a given quasi-interpolation operator onto e Q h . (3) ⇧ H : [L 2 (⌃)] d ! ⇤ H , a given quasi-interpolation operator onto ⇤ H .
TITLE WILL BE SET BY THE PUBLISHER Assumption 4.1. We assume that the interface is well resolved by T Ω h , i.e. h |Σ ∩ K| ≲ |K| ∀K ∈ T Ω h . Assumption 4.2. We assume that every cut cell K ∈ T c h shares at least one vertex with an uncut cell K 1 ⊂ Ω 1 and at least another one with an uncut cell K 2 ⊂ Ω 2 , with K 1 , K 2 ∈ T Ω h \ T c h . These assumptions are fulfilled when the mesh is fine enough.

Approximability

In this section, quasi-interpolation operators are considered to approximate non-smooth functions (see, e.g., [START_REF] Ern | Finite elements I-Approximation and interpolation[END_REF]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)

Π h : V → V h is a given quasi-interpolation operator onto V h ; (2) Π h : Q ε → Q h , is defined as Π h q def = Π h q + Π h q ∀q ∈ Q ε , with Π h q def = 1 |Σ| Σ q 1 Ω1 , (14)
the mean value of the pressure jump across the interface Σ, and

Π h q def = Πh (q -Π h q),
with Πh :

Q ϵ → Q h a given quasi-interpolation operator onto Q h . (3) Π H : [L 2 (Σ)] d → Λ H , a given quasi-interpolation operator onto Λ H .
Owing to the properties of the quasi-interpolation operators (see, e.g., [START_REF] Ern | Finite elements I-Approximation and interpolation[END_REF]Chapter 22]), we have the following result.

Lemma 4.3. Let ε ∈ 0, 1 2 . Assume that p ∈ Q ε and λ ∈ [L 2 (Σ)] d .
For the above defined operators the following estimates hold:

∥ Π h p ∥ L 2 (Σ) ≤ ∥ p ∥ L 2 (Σ) , (15)
∥Π H λ∥ L 2 (Σ) ≲ ∥λ∥ L 2 (Σ) , (16)
h ∇ Π h p L 2 (Ω) ≲ h 1 2 -ε p -Π h p H 1 2 -ε (Ω) , (17)
Π H λ + Π h p n Σ L 2 (Σ) ≲ λ + Π h p n Σ L 2 (Σ) . (18)
Proof. The estimates [START_REF] Burman | Cut finite element method for divergence free approximation of incompressible flow: optimal error estimates and pressure independence[END_REF] and [START_REF] Casquero | The divergence-conforming immersed boundary method: Application to vesicle and capsule dynamics[END_REF] are standard (see, e.g., [START_REF] Ern | Finite elements I-Approximation and interpolation[END_REF]Corollary 22.8]). We focus on estimate [START_REF] Casquero | NURBS-based numerical proxies for red blood cells and circulating tumor cells in microscale blood flow[END_REF].

Let K ∈ T Ω h and p K be the mean value of p -Π h p over S(K), the patch of elements around K. Using the definition of Π h , we have

h ∇ Π h p L 2 (K) = h ∇( Π h p -p K ) L 2 (K) ≲ Π h p -p K L 2 (K) ≲ Πh p -Π h p -p K L 2 (K) ≲ Πh (p -Π h p -p K ) L 2 (K) ≲ p -Π h p -p K L 2 (S(K)) ≲ h 1 2 -ε p -Π h p H 1 2 -ε (S(K))
, where we used an inverse inequality, the local stability of Πh (see, e.g., [START_REF] Ern | Finite elements I-Approximation and interpolation[END_REF]Corollary 22.8]) and approximation properties of the mean value. Summing the square of this relation over all K ∈ T Ω h and using the quasi-uniformity of the mesh gives [START_REF] Casquero | NURBS-based numerical proxies for red blood cells and circulating tumor cells in microscale blood flow[END_REF]. For [START_REF] Burman | Fictitious domain methods using cut elements: III. A stabilized Nitsche method for Stokes' problem[END_REF], we use the fact that Π h p is a constant and that Π h p 1 Ω1 = Π h p, so that we have

∥ Π h p ∥ 2 L 2 (Σ) = |Σ| | Π h p | 2 = |Σ| -1 Σ p 2 ≤ ∥ p ∥ 2 L 2 (Σ)
, where for the last inequality we used Cauchy-Schwarz on the squared term. This completes the proof.

□

For u ∈ V , p ∈ Q ε and λ ∈ [L 2 (Σ)] d
, we introduce the following approximation errors

e u def = u -Π h u, e p def = p -Π h p, e λ def = λ -Π H λ. ( 19
)
We have the following approximation result.

Lemma 4.4. Let ε ∈ 0, 1 2 . Assume that u ∈ V ∩ H 3 2 -ε (Ω), p ∈ Q ε and λ ∈ H 1 2 -ε (Σ).
The following estimates hold:

∥∇ • e u ∥ L 2 (Ω) ≲ ∥ε (e u )∥ L 2 (Ω) ≲ h 1 2 -ε ∥u∥ H 3 2 -ε (Ω) , ( 20
)
∥e u ∥ L 2 (Σ) ≲ h 1-ε ∥u∥ H 3 2 -ε (Ω) , ( 21
)
∥e λ ∥ L 2 (Σ) ≲ H 1 2 -ε ∥λ∥ H 1 2 -ε (Σ) , (22)
∥e p ∥ L 2 (Ω) ≲ h 1 2 -ε p -Π h p H 1 2 -ε (Ω) , (23)
λ + Π h p n Σ -Π H λ + Π h p n Σ L 2 (Σ) ≲ H 1 2 -ε λ + Π h p n Σ H 1 2 -ε (Σ) . (24)
Proof. Estimates ( 20)-( 22) and ( 24) can be easily proven using the approximation properties of quasi-interpolators (see, e.g., [START_REF] Ern | Finite elements I-Approximation and interpolation[END_REF]Theorem 22.6]). For [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by level-sets[END_REF], we use the definition of Π h to get

∥e p ∥ L 2 (Ω) = ∥p -Π h p∥ L 2 (Ω) = p -Π h p -Π h p L 2 (Ω) = p -Π h p -Πh p -Π h p L 2 (Ω) ≲ h 1 2 -ε p -Π h p H 1 2 -ε (Ω)
,

where we used the approximation properties of the quasi-interpolation operator Πh , so that the proof is complete. □

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 , ∇ • v h ) Ω , 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 Σ. This interpolation is used to establish the inf-sup stability of the bilinear form (p h , ∇ • v h ) Ω , 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 [START_REF] Ern | Finite elements I-Approximation and interpolation[END_REF]Lemma 12.8] to our unfitted setting. Similar propositions are also employed in [START_REF] Fournié | Stability and optimal convergence of unfitted extended finite element methods with Lagrange multipliers for the Stokes equations[END_REF][START_REF] Guzmán | Inf-sup stability of geometrically unfitted Stokes finite elements[END_REF].

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 Π 0 h : V → V h such that for every v ∈ V , we have:

• Π 0 h v = 0 on Σ; • Π 0 h v = 0 on ∂Ω if v = 0 on ∂Ω. Moreover, for each v ∈ V such that v = 0 on Σ, it holds v -Π 0 h v L 2 (Ω) ≲ h |v| H 1 (Ω) , (25)
Π 0 h v H 1 (Ω) ≲ ∥v∥ H 1 (Ω) . (26)
Proof. We start by considering the Scott-Zhang interpolation operator Π Z h : V → V h (see, e.g., [START_REF] Scott | Finite element interpolation of nonsmooth functions satisfying boundary conditions[END_REF]). For all v ∈ V , Π 0 h v is defined as the sole element of V h satisfying

Π 0 h v(x) = 0 ∀x ∈ N c h , (27)
Π 0 h v(x) = Π Z h v(x) ∀x ∈ N h \ N c h . (28)
Note that, by construction, we have Π 0 h v = 0 on Σ and

Π Z h v| K = Π 0 h v| K ∀K ∈ T Ω h \ T e h . (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

Π 0 h v = 0 on ∂Ω for all v ∈ V with v = 0 on ∂Ω. Let us now consider v ∈ V such that v = 0 on Σ.
We start with the proof of [START_REF] Fabrèges | A smooth extension method for the simulation of fluid/particles flows[END_REF]. By using the triangle inequality and the approximation properties of the Scott-Zhang operator (see, e.g., [START_REF] Scott | Finite element interpolation of nonsmooth functions satisfying boundary conditions[END_REF]Theorem 4

.1]), it follows that v -Π 0 h v L 2 (Ω) ≤ v -Π Z h v L 2 (Ω) + Π Z h v -Π 0 h v L 2 (Ω) ≲ h |v| H 1 (Ω) + Π Z h v -Π 0 h v L 2 (Ω) . (30)
It only remains to control the last term in the previous bound. Owing to [START_REF] Franca | Error analysis of some Galerkin least squares methods for the elasticity equations[END_REF], we can focus only on the mesh elements in T e h . Let K ∈ T e h , from [24, Proposition 12.5], ( 27) and ( 28), we have

Π Z h v -Π 0 h v L 2 (K) ≲ h d 2 max i Π Z h v(x K,i ) -Π 0 h v(x K,i ) ≲ h d 2 max i: x K,i ∈N c h Π Z h v(x K,i ) -Π 0 h v(x K,i ) ≲ h d 2 max i: x K,i ∈N c h Π Z h v(x K,i ) ≲ h d 2 max i Π Z h v(x K,i ) ≲ Π Z h v L 2 (K) . (31)
In addition, from the stability of the Scott-Zhang interpolation operator (see, e.g., [49, Theorem 3.1]), we have

Π Z h v L 2 (K) ≤ C ∥v∥ L 2 (S(K)) + h |v| H 1 (S(K)) ,
where S(K) denotes the patch of elements around K. Finally, since v vanishes on Σ ∩ S(K) ̸ = ∅, the Poincaré's inequality ∥v∥ L 2 (S(K)) ≲ h |v| H 1 (S(K)) holds (see [START_REF] Ern | Finite elements I-Approximation and interpolation[END_REF]Lemma 3.30]), so that from (31) we get

Π Z h v -Π 0 h v L 2 (K) ≲ h |v| H 1 (S(K)) ∀K ∈ T e h . (32)
Since from [START_REF] Franca | Error analysis of some Galerkin least squares methods for the elasticity equations[END_REF] the estimate (32) holds for every K ∈ T Ω h , we can sum (32) over all K ∈ T Ω h and, by using the quasi-uniformity of the mesh, we finally get

Π Z h v -Π 0 h v L 2 (Ω) ≲ h |v| H 1 (Ω) . (33)
The estimate [START_REF] Fabrèges | A smooth extension method for the simulation of fluid/particles flows[END_REF] simply follows by inserting ( 33) into [START_REF] Galvin | Stabilizing poor mass conservation in incompressible flow problems with large irrotational forcing and application to thermal convection[END_REF].

In order to prove [START_REF] Fabrèges | Approximation of single layer distributions by dirac masses in finite element computations[END_REF], we first notice that owing to [START_REF] Franca | Error analysis of some Galerkin least squares methods for the elasticity equations[END_REF], it suffices to consider the elements K in T e h . By combining [START_REF] Girault | Error analysis of a fictitious domain method applied to a Dirichlet problem[END_REF] with the triangle inequality, we have

Π 0 h v L 2 (K) ≲ Π Z h v L 2 (K) ∀K ∈ T e h . (34)
Similarly, using inverse inequalities and ( 32), for the gradient we have

∇Π 0 h v L 2 (K) ≲ h -1 Π Z h v -Π 0 h v L 2 (K) + ∇Π Z h v L 2 (K) ≲ |v| H 1 (S(K)) + ∇Π Z h v L 2 (K) ∀K ∈ T e h . (35)
Owing to [START_REF] Franca | Error analysis of some Galerkin least squares methods for the elasticity equations[END_REF], the relations (34) and ( 35) hold for every K ∈ T Ω h , so that by summing over K and using the H 1 -stability of the Scott-Zhang interpolator (see, e.g., [START_REF] Scott | Finite element interpolation of nonsmooth functions satisfying boundary conditions[END_REF]Corollary 4.1]), we finally get

Π 0 h v 2 H 1 (Ω) ≲ Π Z h v 2 H 1 (Ω) + |v| 2 H 1 (Ω) ≲ ∥v∥ 2 H 1 (Ω) ,
which completes the proof. □

In the sequel, for i ∈ {1, 2} and for any element q h of Q h , we denote by q h,i def = q h | Ω i the restriction of q h to Ω i . In a similar way, we denote by

Q h,i def = q h | Ω i | q h ∈ Q h the set of such restrictions. Then, for every q h ∈ Q h we have q h = 2 i=1 q h,i 1 Ω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̸ =v h ∈V h ∩[H 1 Σ (Ω)] d (q h , ∇ • v h ) Ω ∥v h ∥ H 1 (Ω) + h ∥∇ q h ∥ L 2 (Ω) ≳ ∥q h ∥ L 2 (Ω) (36)
for all q h ∈ Q h .

Proof. For i ∈ {1, 2}, we proceed in a classical way by decomposing q h,i as

q h,i = qh,i + * q h,i , qh,i def = 1 |Ω i | Ωi q h,i , * q h,i def = q h,i -qh,i ∈ L 2 0 (Ω i ) ∩ Q h,i . (37)
Note that, by restricting the decomposition [START_REF] Barrenechea | A stabilized finite element method for a fictitious domain problem allowing small inclusions[END_REF] to each sub-domain Ω i , we also have q h,i = q h,i + q h,i with q h,i

def = q h | Ω i and q h,i def = q h | Ω i . Furthermore, the fact that * q h,i + qh,i = q h,i = q h,i + q h,i yields ∇ * q h,i = ∇ q h,i . (38)
We start by treating the zero-mean contributions in [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF], namely, *

q h def = i * q h,i 1 Ωi . Since * q h,i ∈ L 2 0 (Ω i ), there exists 0 ̸ = * v i ∈ H 1 0 (Ω i
) such that (see, e.g., [START_REF] Girault | Finite Element Methods for Navier-Stokes Equations[END_REF])

( * q h,i , ∇ • * v i ) Ωi = ∥ * q h,i ∥ 2 L 2 (Ωi) , ∥ * v i ∥ H 1 (Ωi) ≲ ∥ * q h,i ∥ L 2 (Ωi) . (39)
By setting

* v def = i * v i 1 Ωi , we get * v ∈ H 1 ∂Ω∪Σ (Ω) , ( * q h , ∇ • * v) Ω = ∥ * q h ∥ 2 L 2 (Ω) , ∥ * v∥ H 1 (Ω) ≲ ∥ * q h ∥ L 2 (Ω) . (40)
By using the properties of Π 0 h (Proposition 4.5), the fact that (q h,i , ∇ • Π 0 h * v i ) Ωi = 0 and (38), we have

∥ * q h,i ∥ 2 L 2 (Ωi) = ( * q h,i , ∇ • * v i ) Ωi = ( * q h,i , ∇ • ( * v i -Π 0 h * v i )) Ωi + ( * q h,i , ∇ • Π 0 h * v i ) Ωi = -(∇ * q h,i , * v i -Π 0 h * v i ) Ωi + ( * q h,i , ∇ • Π 0 h * v i ) Ωi ≲ h ∥∇ * q h,i ∥ L 2 (Ωi) ∥ * v i ∥ H 1 (Ωi) + (q h,i , ∇ • Π 0 h * v i ) Ωi = h ∥∇ q h,i ∥ L 2 (Ωi) ∥ * v i ∥ H 1 (Ωi) + (q h,i , ∇ • Π 0 h * v i ) Ωi .
Notice that here we have made a little abuse of notation in Π 0 h * v i , since it has to be understood as Π 0 h applied to the extension of * v i by zero out of Ω i . By summing the above relation over i ∈ {1, 2} and since Π 0

h * v = Π 0 h * v 1 + Π 0 h * v 2 (
which follows from the linearity of Π 0 h and Proposition 4.5), we get

∥ * q h ∥ 2 L 2 (Ω) ≲ h ∥∇ q h ∥ L 2 (Ω) ∥ * v∥ H 1 (Ω) + (q h , ∇ • Π 0 h * v) Ω ,
which, with (40), yields,

∥ * q h ∥ L 2 (Ω) ≲ h ∥∇ q h ∥ L 2 (Ω) + (q h , ∇ • Π 0 h * v) Ω ∥ * v∥ H 1 (Ω) . (41)
Let us now handle the constant part of the pressure qh,i . Let Ω i be a subdomain of Ω i such that Γ N ∩ Ω i ⊂ Ω i and such that the mesh is fitted to Ω i . From [21, Lemma 3.3] (see also [START_REF] Franca | Error analysis of some Galerkin least squares methods for the elasticity equations[END_REF]Lemma 3.3]), there exists a constant

C i > 0 and 0 ̸ = v h,i ∈ {w h | Ωi | w h ∈ V h } ∩ [H 1 ∂ Ωi\ΓN ( Ω i )] d , such that C i ∥q h,i ∥ L 2 ( Ωi) ∥ v h,i ∥ H 1 ( Ωi) ≤ (q h,i , ∇ • v h,i ) Ωi . (42)
Notice that, since Ω i is fitted to the mesh, the constant C i is independent from the way the unfitted interface cuts the mesh. Moreover, denoting by vh,i

∈ {w h | Ωi | w h ∈ V h } ∩ [H 1 ∂Ωi\ΓN (Ω i )]
d the extension by zero of v h,i to the whole domain Ω i , the estimate (42) holds with Ω i = Ω i , v h,i = vh,i and with the constant C i rescaled by | Ω i |/|Ω i |. According to Assumption 2.1, Ω i can be chosen in a way that the ratio | Ω i |/|Ω i | does not vanish when h tends to zero. We established

∥q h,i ∥ L 2 (Ωi) ≲ (q h,i , ∇ • vh,i ) Ωi ∥v h,i ∥ H 1 (Ωi) . (43)
We now consider

zh def = i 1 √ 2 ∥v h,i ∥ H 1 (Ωi) vh,i 1 Ωi , so that ∥z h ∥ H 1 (Ω) = 1
, and summing [START_REF] Liu | Rheology of red blood cell aggregation by computer simulation[END_REF] for i ∈ {1, 2} we have

∥q h ∥ L 2 (Ω) ≲ (q h , ∇ • zh ) Ω = (q h , ∇ • zh ) Ω ∥z h ∥ H 1 (Ω) ≲ (q h , ∇ • zh ) Ω ∥z h ∥ H 1 (Ω) + ∥ * q h ∥ L 2 (Ω) . (44)
Setting [START_REF] Guzmán | Inf-sup stability of geometrically unfitted Stokes finite elements[END_REF] follows from ( 41) and [START_REF] Nakata | A fluid-structure interaction model of insect flight with flexible wings[END_REF], which completes the proof. □ Proposition 4.7. Under Assumptions 4.1 and 4.2, we have

w h = ∥ * v∥ -1 H 1 (Ω) Π 0 h * v + δ ∥z h ∥ -1 H 1 (Ω) zh , with δ > 0 sufficiently small, results in w h ∈ V h ∩ [H 1 Σ (Ω)] d and ∥w h ∥ H 1 (Ω) ≲ 1 + δ. The relation
h 1 2 ∥ q h ∥ L 2 (Σ) ≲ ∥q h ∥ L 2 (Ω) + h ∥∇ q h ∥ L 2 (Ω) (45)
for all q h ∈ Q h .

Proof. Let us introduce the operators

E 1 h : Q h → Q h and E 2 h : Q h → Q h defined as E 1 h (q h ) def = q h + q h , E 2 h (q h ) def = q h
for all q h ∈ Q h . We recall that by construction of the pressure space Q h in (6), the pressure jump q h is constant. In particular, we have

E i h (q h ) = q h,i in Ω i . ( 46
) Let K ∈ T c
h be given. Since, q h is constant on Σ, using Assumption 4.1 and the triangle inequality, we have

h 1 2 ∥ q h ∥ L 2 (Σ∩K) = h 1 2 |Σ ∩ K| 1 2 | q h | ≲ ∥ q h ∥ L 2 (K) = E 1 h (q h ) -E 2 h (q h ) L 2 (K) ≤ i E i h (q h ) L 2 (K) . ( 47
)
In order to bound E i h (q h ) L 2 (K) for every i ∈ {1, 2}, we first note that since E i h (q h ) is affine in K, we have

E i h (q h )(x) = E i h (q h )(x K,j ) + ∇E i h (q h ) • (x -x K,j ),
where x K,j denotes one of the nodes of the element K. Thus, using an inverse inequality we get

E i h (q h ) L 2 (K) ≲ h d 2 E i h (q h ) L ∞ (K) ≲ h d 2 E i h (q h )(x K,j ) + h 1+ d 2 ∇E i h (q h ) L ∞ (K) . (48)
Furthermore, owing to Assumption 4.2, relation [START_REF] Ohmori | Flux-free finite element method with Lagrange multipliers for two-fluid flows[END_REF] and [24, Proposition 12.5], we have

h d 2 E i h (q h )(x K,j ) ≤ h d 2 max l q h,i (x K i ,l ) ≲ ∥q h,i ∥ L 2 ( Ki) . (49)
In addition, since

∇E i h (q h ) is constant in K and ∇E 1 h (q h ) = ∇E 2 h (q h ) = ∇ q h , we have h d 2 ∇E i h (q h ) L ∞ (K) ≲ ∥∇ q h ∥ L 2 (K) . (50)
By inserting the bounds ( 48)-( 50) into (47), we have

h 1 2 ∥ q h ∥ L 2 (K∩Σ) ≲ i E i h (q h ) L 2 (K) ≲ ∥q h ∥ L 2 ( K1) + ∥q h ∥ L 2 ( K2) + h ∥∇ q h ∥ L 2 (K) .
The estimate [START_REF] Nguyen | Boundary stabilization of the Navier-Stokes equations in the case of mixed boundary conditions[END_REF] then follows by exploiting the quasi-uniformity of the mesh and summing the above relation over all K ∈ T c h , which completes the proof. □

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 , ξ H )∥ def = µ ∥ε (v h )∥ 2 L 2 (Ω) + h 2 µ -1 ∥∇ q h ∥ 2 L 2 (Ω) + µ -1 ∥q h ∥ 2 L 2 (Ω) + hµ -1 ∥ξ H ∥ 2 L 2 (Σ) 1 2 for all (v h , q h , ξ H ) ∈ V h × Q h × Λ H .
The following lemma provides the inf-sup stability for the full bilinear form A h . Lemma 4.8. Let Assumptions 2.1, 4.1 and 4.2 be fulfilled. For θ = 0, we furthermore assume that γ λ is large enough. Then, there exists a constant β > 0 independent of µ, the mesh and the interface position such that inf

(u h ,p h ,λ H ) ∈ V h ×Q h ×Λ H \(0,0,0) sup (v h ,q h ,ξ H ) ∈ V h ×Q h ×Λ H \(0,0,0) A h ( (u h , p h , λ H ) , (v h , q h , ξ H ) ) ∥(u h , p h , λ H )∥ ∥(v h , q h , ξ H )∥ ≥ β. ( 51
)
TITLE WILL BE SET BY THE PUBLISHER Proof. The inf-sup condition ( 51) is a consequence of the following estimate

∥u h , p h , λ H ∥ ≲ sup (v h ,q h ,ξ H ) ∈ V h ×Q h ×Λ H \(0,0,0) A h ( (u h , p h , λ H ) , (v h , q h , ξ H ) ) ∥(v h , q h , ξ H )∥ def = S (52)
for all

(u h , p h , λ H ) ∈ V h × Q h × Λ H \ (0, 0, 0)
. Now, let us prove [START_REF] Takizawa | Space-time interface-tracking with topology change (ST-TC)[END_REF]. Using Assumption 2.1, we can invoke Proposition 4.6. Let w h ∈ V h ∩ H 1 Σ (Ω) be a function for which the supremum in [START_REF] Guzmán | Inf-sup stability of geometrically unfitted Stokes finite elements[END_REF] is attained (for q h = p h ). We have

∥p h ∥ L 2 (Ω) ≲ (p h , ∇ • w h ) Ω ∥w h ∥ H 1 (Ω) + h ∥∇ p h ∥ L 2 (Ω) . (53)
In addition, from [START_REF] Boffi | Stabilized Stokes elements and local mass conservation[END_REF], we have

(p h , ∇ • w h ) Ω = A h ( (u h , p h , λ H ) , (-w h , 0, 0) ) + 2µ(ε (u h ) , ε (w h )) Ω ≲ S ∥(w h , 0, 0)∥ + 2µ ∥ε (u h )∥ L 2 (Ω) ∥ε (w h )∥ L 2 (Ω) ≲ (µ 1 2 S + 2µ ∥ε (u h )∥ L 2 (Ω) ) ∥w h ∥ H 1 (Ω) ,
so that (53) yields

∥p h ∥ L 2 (Ω) ≲ µ 1 2 S + 2µ ∥ε (u h )∥ L 2 (Ω) + h ∥∇ p h ∥ L 2 (Ω) . (54)
Case 11), we have

θ = 1. Let (u h , p h , λ H ) ∈ V h × Q h × Λ H \ (0, 0, 0). From (
2µ ∥ε (u h )∥ 2 L 2 (Ω) + γ p h 2 µ ∥∇ p h ∥ 2 L 2 (Ω) + h γ λ µ ∥λ H + p h n Σ ∥ 2 L 2 (Σ) = A h ( (u h , p h , λ H ) , (u h , p h , λ H ) ) (55)
≤ S ∥(u h , p h , λ H )∥ .

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

h ∥λ H ∥ 2 L 2 (Σ) ≤ h ∥λ H + p h n Σ ∥ 2 L 2 (Σ) + h ∥ p h ∥ 2 L 2 (Σ) ≲ h ∥λ H + p h n Σ ∥ 2 L 2 (Σ) + ∥p h ∥ 2 L 2 (Ω) + h 2 ∥∇ p h ∥ 2 L 2 (Ω) . (56)
We now conclude by inserting ( 54) and ( 56) (multiplied by α

1 2
p µ -1 2 and α λ µ -1 γ -1 λ , respectively for α p , α p two positive coefficients) into [START_REF] Weymouth | Advancements in Cartesian-grid methods for computational ship hydrodynamics[END_REF], which yields

C 1 (γ p , γ λ , α p , α λ ) ∥(u h , p h , λ H )∥ 2 ≲ α p S 2 + S ∥(u h , p h , λ H )∥ , with C 1 (γ p , γ λ , α p , α λ ) def = min 2 -4α p , γ p -α p - α λ γ λ , α p - α λ γ λ , α λ γ λ ,
which is positive for an adequate choice of α p and α λ . Then (52) simply follows from Young's inequality, which completes the proof. Case θ = 0. Similarly to the previous case, let

(u h , p h , λ H ) ∈ V h × Q h × Λ H \ (0, 0, 0). From (11), we have 2µ ∥ε (u h )∥ 2 L 2 (Ω) + γ p h 2 µ ∥∇ p h ∥ 2 L 2 (Ω) + h γ λ µ (λ H + p h n Σ , λ H ) Σ = A h ( (u h , p h , λ H ) , (u h , p h , λ H ) ) (57) ≤ S ∥(u h , p h , λ H )∥ .
Moreover, using Cauchy-Schwarz and Young inequalities and Proposition 4.7 (under Assumptions 4.1 and 4.2), we get

h 2 ∥λ H ∥ 2 L 2 (Σ) ≤ h(λ H + p h n Σ , λ H ) L 2 (Σ) + C λ 2 ∥p h ∥ 2 L 2 (Ω) + C λ h 2 2 ∥∇ p h ∥ 2 L 2 (Ω) . (58)
We now conclude by inserting ( 54) and (58) (multiplied respectively by α

1 2
p µ -1 2 and γ -1 λ µ -1 for α p a positive coefficient) into (57), which yields

C 0 (γ p , γ λ , α p ) ∥(u h , p h , λ H )∥ 2 ≲ α p S 2 + S ∥(u h , p h , λ H )∥ , with C 0 (γ p , γ λ , α p ) def = min 2 -4α p , γ p -α p - C λ 2γ λ , α p - C λ 2γ λ , 1 2γ λ ,
where C λ > 0 is the hidden constant in the statement of Proposition 4.7. The constant C 0 (γ p , γ λ , α p ) is positive for an adequate choice of α p , under the assumption that γ λ is large enough. Then (52) simply follows from Young's inequality, which completes the proof. □ Corollary 4.9. The discrete problem (10) admits a unique solution.

Remark 4.10. In the case where θ = 0, the requirement for γ λ to be large enough stems from the definition of the constant C 0 (γ p , γ λ , α p ) (as detailed in the proof). This constitutes a limitation of the θ = 0 variant of the proposed method. In practical terms, we aim for a γ λ 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).

A priori error estimates

We first introduce the following discrete errors

e u,h = Π h u -u h , e p,h = Π h p -p h , e λ,H = Π H λ -λ H ,
where (u, p, λ) is the solution to the continuous problem (2), (u h , p h , λ H ) is the solution to the discrete problem [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] and the interpolation operators Π h , Π h and Π H are defined in Section 4.1. The following result provides an estimate for the discrete error.

Lemma 4.11. Assume that (u, p, λ) ∈ V ε × Q ε × H 1 2 -ε (Σ) with ε ∈ (0, 1/2).
The following estimate holds

|A h ( (e u,h , e p,h , e λ,H ) , (v h , q h , ξ H ) )| ≲ max{h, H} 1 2 -ε µ 1 2 ∥u∥ H 3 2 -ε (Ω) + µ -1 2 p -Π h p H 1 2 -ε (Ω) +µ -1 2 ∥ p ∥ L 2 (Σ) + µ -1 2 ∥λ∥ H 1 2 -ε (Σ) ∥v h , q h , ξ H ∥ (59) for all (v h , q h , ξ H ) ∈ V h × Q h × Λ H .
Moreover, this estimate can be improved to

|A h ( (e u,h , e p,h , e λ,h ) , (v h , q h , ξ H ) )| ≲ max{h, H} 1 2 -ε µ 1 2 ∥u∥ H 3 2 -ε (Ω) + µ -1 2 p -Π h p H 1 2 -ε (Ω) + µ -1 2 λ + Π h p n Σ H 1 2 -ε (Σ) +µ -1 2 ∥ p ∥ L 2 (Σ) ∥n Σ -Π H n Σ ∥ L 2 (Σ) ∥v h , q h , ξ H ∥ . ( 60
)

16

TITLE WILL BE SET BY THE PUBLISHER Proof. Using the relations ( 2) and ( 10) and the definitions (3), ( 11) and ( 19), we have

A h ( (e u,h , e p,h , e λ,H ) , (v h , q h , ξ H ) ) = A h ( (Π h u, Π h p, Π H λ) , (v h , q h , ξ H ) ) -F (v h , ξ H ) = A h ( (Π h u, Π h p, Π H λ) , (v h , q h , ξ H ) ) -A( (u, p, λ) , (v h , q h , ξ H ) ) = -2µ (ε (e u ) , ε (v h )) Ω + (e p , ∇ • v h ) Ω -(q h , ∇ • e u ) Ω + (e λ , v h ) Σ -(ξ H , e u ) Σ + s BP h Π h p, q h + s BH h ( (Π h p, Π H λ) , (q h , ξ H ) ). (61)
We now proceed by estimating term by term. Using the Cauchy-Schwarz inequality and the approximation error bounds ( 20), ( 21) and ( 23), we have

µ| (ε (e u ) , ε (v h )) Ω | ≤ µ ∥ε(e u )∥ L 2 (Ω) ∥ε (v h )∥ L 2 (Ω) ≲ µ 1 2 h 1 2 -ε ∥u∥ H 3 2 -ε (Ω) ∥(v h , q h , ξ H )∥ , | (q h , ∇ • e u ) Ω | ≤ ∥q h ∥ L 2 (Ω) ∥∇ • e u ∥ L 2 (Ω) ≲ µ 1 2 h 1 2 -ε ∥u∥ H 3 2 -ε (Ω) ∥(v h , q h , ξ H )∥ , | (ξ H , e u ) Σ | ≤ ∥ξ H ∥ L 2 (Σ) ∥e u ∥ L 2 (Σ) ≲ µ 1 2 h 1 2 -ε ∥u∥ H 3 2 -ε (Ω) ∥(v h , q h , ξ H )∥ , | (e p , ∇ • v h ) Ω | ≤ ∥e p ∥ L 2 (Ω) ∥∇ • v h ∥ L 2 (Ω) ≲ µ -1 2 h 1 2 -ε p -Π h p H 1 2 -ε (Ω) ∥(v h , q h , ξ H )∥ . (62)
For the pressure stabilization term, we use Cauchy-Schwarz inequality and ( 17) to get

|s BP h Π h p, q h | ≲ µ -1 h 2 | ∇ Π h p, ∇ q h Ω | ≲ µ -1 h ∇ Π h p L 2 (Ω) h ∥∇ q h ∥ L 2 (Ω) ≲ µ -1 2 h 1 2 -ε p -Π h p H 1 2 -ε (Ω) ∥(v h , q h , ξ H )∥ . ( 63
)
It only remains to estimate the terms (e λ , v h ) Σ and s BH h ( (Π h p, Π H λ) , (q h , ξ H ) ). Proof of (59). For the term (e λ , v h ) Σ , using the Cauchy-Schwarz, trace and Korn's inequalities and [START_REF] Duprez | ϕ-FEM: an optimally convergent and easily implementable immersed boundary method for particulate flows and Stokes equations[END_REF], we have

| (Π H λ -λ, v h ) Σ | ≲ H 1 2 -ε ∥λ∥ H 1 2 -ε (Σ) ∥ε (v h )∥ L 2 (Ω) ≲ µ -1 2 H 1 2 -ε ∥λ∥ H 1 2 -ε (Σ) ∥(v h , q h , ξ H )∥ . (64)
Finally, for the term s BH h ( (Π h p, Π H λ) , (q h , ξ H ) ), we use once more the Cauchy-Schwarz inequality, Proposition 4.7 (for θ = 1), [START_REF] Burman | Fictitious domain methods using cut elements: III. A stabilized Nitsche method for Stokes' problem[END_REF] and ( 16), to get

µ -1 h | (Π H λ + Π h p n Σ , ξ H + θ q h n Σ ) Σ | ≲ µ -1 h 1 2 ∥Π H λ∥ L 2 (Σ) + ∥ Π h p ∥ L 2 (Σ) h 1 2 ∥ξ H ∥ L 2 (Σ) + θ ∥ q h ∥ L 2 (Σ) ≲ µ -1 2 h 1 2 ∥λ∥ L 2 (Σ) + ∥ p ∥ L 2 (Σ) ∥(v h , q h , ξ H )∥ . ( 65
)
The estimate (59) then follows from (61) and the estimates (62)-(65).

Proof of (60). We first recall that, from [START_REF] Brezzi | On the stabilization of finite element approximations of the Stokes equations[END_REF], Π 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

| (Π H λ -λ, v h ) Σ | ≤ Π H λ + Π h p n Σ -λ + Π h p n Σ + Π h p (n Σ -Π H n Σ ) L 2 (Σ) ∥v h ∥ L 2 (Σ) ≲ µ -1 2 (H 1 2 -ε λ + Π h p n Σ H 1 2 -ε (Σ) + | Π h p | ∥(n Σ -Π H n Σ )∥ L 2 (Σ) ) ∥(v h , q h , ξ H )∥ . (66)
Similarly, for the last term we use the Cauchy-Schwarz inequality, Proposition 4.7 (for θ = 1), ( 18) and ( 15), which yields

µ -1 h | ( Π H λ + Π h p n Σ , ξ H + θ q h n Σ ) Σ | ≲ µ -1 h 1 2 Π H λ + Π h p n Σ + Π h p (n Σ -Π H n Σ ) L 2 (Σ) h 1 2 ∥ξ H ∥ L 2 (Σ) + θ ∥ q h ∥ L 2 (Σ) ≲ µ -1 2 h 1 2 λ + Π h p n Σ L 2 (Σ) + Π h p ∥n Σ -Π H n Σ ∥ L 2 (Σ) ∥(v h , q h , ξ H )∥ ≲ µ -1 2 h 1 2 λ + Π h p n Σ L 2 (Σ) + ∥ p ∥ L 2 (Σ) ∥n Σ -Π H n Σ ∥ L 2 (Σ) ∥(v h , q h , ξ H )∥ . (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 ε ∈ (0, 1/2) and assume that (u, p, λ)

∈ V ε × Q ε × H 1 2 -ε (Σ)
is the solution to the continuous problem (2) and that (u h , p h , λ H ) is the solution to the discrete problem [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF]. Furthermore, in the case θ = 0, assume that γ λ is large enough. The following a priori error bounds hold:

µ 1 2 ∥ε (u -u h )∥ L 2 (Ω) + µ -1 2 ∥p -p h ∥ L 2 (Ω) + µ -1 2 h 1 2 ∥λ -λ H ∥ L 2 (Σ) ≲ max{h, H} 1 2 -ε µ 1 2 ∥u∥ H 3 2 -ε (Ω) + µ -1 2 p -Π h p H 1 2 -ε (Ω) + µ -1 2 ∥ p ∥ L 2 (Σ) + µ -1 2 ∥λ∥ H 1 2 -ε (Σ) , (68)
µ 1 2 ∥ε (u -u h )∥ L 2 (Ω) + µ -1 2 ∥p -p h ∥ L 2 (Ω) + µ -1 2 h 1 2 ∥λ -λ H ∥ L 2 (Σ) ≲ max{h, H} 1 2 -ε µ 1 2 ∥u∥ H 3 2 -ε (Ω) + µ -1 2 p -Π h p H 1 2 -ε (Ω) + µ -1 2 λ + Π h p n Σ H 1 2 -ε (Σ) + µ -1 2 ∥ p ∥ L 2 (Σ) ∥n Σ -Π H n Σ ∥ L 2 (Σ) . ( 69
)
Proof. The estimate (68) follows immediately using the classical error decomposition

µ 1 2 ∥ε (u -u h )∥ L 2 (Ω) +µ -1 2 ∥p -p h ∥ L 2 (Ω) + µ -1 2 h 1 2 ∥λ -λ H ∥ L 2 (Σ) ≤µ 1 2 ∥ε (e u )∥ L 2 (Ω) + µ -1 2 ∥e p ∥ L 2 (Ω) + µ -1 2 h 1 2 ∥e λ ∥ L 2 (Σ) + ∥(e u,
(v h ,q h ,ξ H ) ∈ V h ×Q h ×Λ H \(0,0,0) A h ((e u,h , e p,h , e λ,H ) ; (v h , q h , ξ H )) ∥v h , q h , ξ H ∥ ≲max{h, H} 1 2 -ε (µ 1 2 ∥u∥ H 3 2 -ε (Ω) + µ -1 2 p -Π h p H 1 2 -ε (Ω) + µ -1 2 ∥ p ∥ L 2 (Σ) + µ -1 2 ∥λ∥ H 1 2 -ε (Σ)
).

The estimate (68) then follows from (20), ( 22) and [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by level-sets[END_REF]. 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 [START_REF] Baaijens | A fictitious domain/mortar element method for fluid-structure interaction[END_REF] 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, TITLE WILL BE SET BY THE PUBLISHER the mean value of the pressure jump Π h p Σ is subtracted to the Lagrange multiplier, which is expected to improve the constant of the error bound. Moreover, the term ∥ p ∥ L 2 (Σ) ∥n Σ -Π H n Σ ∥ L 2 (Σ) that carries the pressure jump has an improved convergence rate in this estimate, assuming that n Σ has a better regularity than H 1 2 -ε (Σ). 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 -1 2 (Σ). Elements of such a proof are given in [3, Propositions A.1, A.2].

Numerical experiments

In this section, we illustrate the behavior of the numerical method [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] in three different two-dimensional examples. In all these test cases, we consider a rectangular domain Ω = (-1, 1) × (0, 1). Dirichlet boundary conditions are imposed at the top and bottom boundaries Γ D = (-1, 1) × {0, 1}, and Neumann boundary conditions are enforced on the remaining boundary Γ N = {-1, 1} × (0, 1). A different interface Σ is considered for each test case, as shown in Figure 3. Without loss of generality, the interface is oriented so that the sub-domain Ω 1 (resp. Ω 2 ) lies on the left (resp. right) side of the interface.

18

TITLE WILL BE SET BY THE PUBLISHER the mean value of the pressure jump Jb ⇧ h pK ⌃ is subtracted to the Lagrange multiplier, which is expected to improve the constant of the error bound. Moreover, the term kJpKk L 2 (⌃) kn ⌃ ⇧ H n ⌃ k L 2 (⌃) that carries the pressure jump has an improved convergence rate in this estimate, assuming that n ⌃ has a better regularity than H 1 2 " (⌃). 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

Numerical experiments

In this section, we illustrate the behavior of the numerical method [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] in three different two-dimensional examples. In all these test cases, we consider a rectangular domain ⌦ = ( 1, 1) ⇥ (0, 1). Dirichlet boundary conditions are imposed at the top and bottom boundaries D = ( 1, 1) ⇥ {0, 1}, and Neumann boundary conditions are enforced on the remaining boundary N = { 1, 1} ⇥ (0, 1). A different interface ⌃ is considered for each test case, as shown in Figure 3. Without loss of generality, the interface is oriented so that the sub-domain ⌦ 1 (resp. ⌦ 2 ) lies on the left (resp. right) side of the interface. In the subsequent sections, all physical quantities and parameters are expressed in the centimeter-gram-second (CGS) unit system, the stabilization parameters are p = 10 2 and = 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 [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] with the pressure space e

• • D N D N ⌦ 1 ⌦ 2 ⌃ (a) Straight interface. • • D N D N ⌦ 1 ⌦ 2 ⌃ (b) Curved interface.
Q h instead of Q h , namely: Find (u h , p h , H ) 2 V h ⇥ e Q h ⇥ ⇤ H such that A h ( (u h , p h , H ) , (v h , q h , ⇠ H ) ) = F (v h , ⇠ H ) for all (v h , q h , ⇠ H ) 2 V h ⇥ e Q h ⇥⇤ H .
It should be noted that for this method the Lagrange multiplier stabilization term (9) reduces to the one considered in [START_REF] Barrenechea | A stabilized finite element method for a fictitious domain problem allowing small inclusions[END_REF], namely,

s BH h ( H , ⇠ H ) def = h µ ( H , ⇠ H ) ⌃ .

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 ⌃ = {0} ⇥ (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.

18

TITLE WILL BE SET BY THE PUBLISHER the mean value of the pressure jump Jb ⇧ h pK ⌃ is subtracted to the Lagrange multiplier, which is expected to improve the constant of the error bound. Moreover, the term kJpKk L 2 (⌃) kn ⌃ ⇧ H n ⌃ k L 2 (⌃) that carries the pressure jump has an improved convergence rate in this estimate, assuming that n ⌃ has a better regularity than H 1 2 " (⌃). 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

Numerical experiments

In this section, we illustrate the behavior of the numerical method [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] in three different two-dimensional examples. In all these test cases, we consider a rectangular domain ⌦ = ( 1, 1) ⇥ (0, 1). Dirichlet boundary conditions are imposed at the top and bottom boundaries D = ( 1, 1) ⇥ {0, 1}, and Neumann boundary conditions are enforced on the remaining boundary N = { 1, 1} ⇥ (0, 1). A different interface ⌃ is considered for each test case, as shown in Figure 3. Without loss of generality, the interface is oriented so that the sub-domain ⌦ 1 (resp. ⌦ 2 ) lies on the left (resp. right) side of the interface. In the subsequent sections, all physical quantities and parameters are expressed in the centimeter-gram-second (CGS) unit system, the stabilization parameters are p = 10 2 and = 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 [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] with the pressure space e

• • D N D N ⌦ 1 ⌦ 2 ⌃ (a) Straight interface. • • D N D N ⌦ 1 ⌦ 2 ⌃ (b) Curved interface.
Q h instead of Q h , namely: Find (u h , p h , H ) 2 V h ⇥ e Q h ⇥ ⇤ H such that A h ( (u h , p h , H ) , (v h , q h , ⇠ H ) ) = F (v h , ⇠ H ) for all (v h , q h , ⇠ H ) 2 V h ⇥ e Q h ⇥⇤ H .
It should be noted that for this method the Lagrange multiplier stabilization term (9) reduces to the one considered in [START_REF] Barrenechea | A stabilized finite element method for a fictitious domain problem allowing small inclusions[END_REF], namely,

s BH h ( H , ⇠ H ) def = h µ ( H , ⇠ H ) ⌃ .

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 ⌃ = {0} ⇥ (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.

18

TITLE WILL BE SET BY THE PUBLISHER the mean value of the pressure jump Jb ⇧ h pK ⌃ is subtracted to the Lagrange multiplier, which is expected to improve the constant of the error bound. Moreover, the term kJpKk L 2 (⌃) kn ⌃ ⇧ H n ⌃ k L 2 (⌃) that carries the pressure jump has an improved convergence rate in this estimate, assuming that n ⌃ has a better regularity than H 1 2 " (⌃). 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

Numerical experiments

In this section, we illustrate the behavior of the numerical method [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] in three different two-dimensional examples. In all these test cases, we consider a rectangular domain ⌦ = ( 1, 1) ⇥ (0, 1). Dirichlet boundary conditions are imposed at the top and bottom boundaries D = ( 1, 1) ⇥ {0, 1}, and Neumann boundary conditions are enforced on the remaining boundary N = { 1, 1} ⇥ (0, 1). A different interface ⌃ is considered for each test case, as shown in Figure 3. Without loss of generality, the interface is oriented so that the sub-domain ⌦ 1 (resp. ⌦ 2 ) lies on the left (resp. right) side of the interface. In the subsequent sections, all physical quantities and parameters are expressed in the centimeter-gram-second (CGS) unit system, the stabilization parameters are p = 10 2 and = 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 [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] with the pressure space e

• • D N D N ⌦ 1 ⌦ 2 ⌃ (a) Straight interface. • • D N D N ⌦ 1 ⌦ 2 ⌃ (b) Curved interface.
Q h instead of Q h , namely: Find (u h , p h , H ) 2 V h ⇥ e Q h ⇥ ⇤ H such that A h ( (u h , p h , H ) , (v h , q h , ⇠ H ) ) = F (v h , ⇠ H ) for all (v h , q h , ⇠ H ) 2 V h ⇥ e Q h ⇥⇤ H .
It should be noted that for this method the Lagrange multiplier stabilization term (9) reduces to the one considered in [START_REF] Barrenechea | A stabilized finite element method for a fictitious domain problem allowing small inclusions[END_REF], namely,

s BH h ( H , ⇠ H ) def = h µ ( H , ⇠ H ) ⌃ .

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 ⌃ = {0} ⇥ (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. In the subsequent sections, all physical quantities and parameters are expressed in the centimeter-gram-second (CGS) unit system, the stabilization parameters γ p = 10 -2 and γ λ = 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 [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] with the pressure space

Q h instead of Q h , namely: Find (u h , p h , λ H ) ∈ V h × Q h × Λ H such that A h ( (u h , p h , λ H ) , (v h , q h , ξ H ) ) = F (v h , ξ H ) for all (v h , q h , ξ H ) ∈ V h × Q h ×Λ H .
It should be noted that for this method the Lagrange multiplier stabilization term (9) reduces to the one considered in [START_REF] Barrenechea | A stabilized finite element method for a fictitious domain problem allowing small inclusions[END_REF], namely,

s BH h (λ H , ξ H ) def = h γ λ µ (λ H , ξ H ) Σ .

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 Σ = {0} × (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. For i 2 {1, 2}, the velocity field u i def = u| ⌦i in ( 1) is chosen as

10 0 h k"(u u h )k L 2 (⌦) k"(u)k L 2 (⌦) h 1 2 ✓ = 1 ✓ = 0 FD-LM (
10 0 h k H k L 2 (⌃) k k L 2 (⌃) h 1 2 ✓ = 1 ✓ = 0 FD-LM (
u 1 = ✓ @ @y 1 (x, y) , @ @x 1 (x, y) ◆ , u 2 = ✓ @ @y 2 (x, y) , @ @x 2 (x, y) ◆ ,
with the function i (x, y) defined as For i 2 {1, 2}, the velocity field u i def = u| ⌦i in ( 1) is chosen as

i (x, y) def = (x x l,i ) 2 (x x r,i ) 2 (y y b,i ) 2 (y y t,i ) 2 ' i (x, y) ,
10 0 h k"(u u h )k L 2 (⌦) k"(u)k L 2 (⌦) h 1 2 ✓ = 1 ✓ = 0 FD-LM (
10 0 h k H k L 2 (⌃) k k L 2 (⌃) h 1 2 ✓ = 1 ✓ = 0 FD-LM (
u 1 = ✓ @ @y 1 (x, y) , @ @x 1 (x, y) ◆ , u 2 = ✓ @ @y 2 (x, y) , @ @x 2 (x, y) ◆ ,
with the function i (x, y) defined as For i 2 {1, 2}, the velocity field u i def = u| ⌦i in ( 1) is chosen as

i (x, y) def = (x x l,i ) 2 (x x r,i ) 2 (y y b,i ) 2 (y y t,i ) 2 ' i (x, y) , (
10 0 h k"(u u h )k L 2 (⌦) k"(u)k L 2 (⌦) h 1 2 ✓ = 1 ✓ = 0 FD-LM (
10 0 h k H k L 2 (⌃) k k L 2 (⌃) h 1 2 ✓ = 1 ✓ = 0 FD-LM (
u 1 = ✓ @ @y 1 (x, y) , @ @x 1 (x, y) ◆ , u 2 = ✓ @ @y 2 (x, y) , @ @x 2 (x, y) ◆ ,
with the function i (x, y) defined as For i 2 {1, 2}, the velocity field u i def = u| ⌦i in ( 1) is chosen as

i (x, y) def = (x x l,i ) 2 (x x r,i ) 2 (y y b,i ) 2 (y y t,i ) 2 ' i (x,
10 0 h k"(u u h )k L 2 (⌦) k"(u)k L 2 (⌦) h 1 2 ✓ = 1 ✓ = 0 FD-LM (
10 0 h k H k L 2 (⌃) k k L 2 (⌃) h 1 2 ✓ = 1 ✓ = 0 FD-LM (
u 1 = ✓ @ @y 1 (x, y) , @ @x 1 (x, y) ◆ , u 2 = ✓ @ @y 2 (x, y) , @ @x 2 (x, y) ◆ ,
with the function i (x, y) defined as For i ∈ {1, 2}, the velocity field u i def = u| Ωi in ( 1) is chosen as

i (x, y) def = (x x l,i ) 2 (x x r,i ) 2 (y y b,i ) 2 (y y t,i ) 2 ' i (x, y) , (
u 1 = ∂ ∂y ϕ 1 (x, y) , - ∂ ∂x ϕ 1 (x, y) , u 2 = ∂ ∂y ϕ 2 (x, y) , - ∂ ∂x ϕ 2 (x, y) ,
with the function ϕ i (x, y) defined as

ϕ i (x, y) def = (x -x l,i ) 2 (x -x r,i ) 2 (y -y b,i ) 2 (y -y t,i ) 2 φ i (x, y) ,
with φ 1 (x, y) def = x, φ 2 (x, y) 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:

u 1 (x, y) = 20 (x + 1) 2 x 3 y (y -1) (2y -1) 10 (x + 1) x 2 y 2 (y -1) 2 (5x + 3) , p 1 (x, y) = x + 1, u 2 (x, y) = 10 (y -1) x 2 y 2 (x -1) 2 (5y -3) 20 (y -1) 2 y 3 x (x -1) (1 -2x) , p 2 (x, y) = y, (70)
10 0 h k"(u u h )k L 2 (⌦) k"(u)k L 2 (⌦) H = h/3 H = h H = h 1 2 h 1 2
(a) Velocity relative energy error.

10 2 10 1 10 2 10 1 h kp p h k L 2 (⌦) kpk L 2 (⌦) H = h/3 H = h H = h 1 2 h 1 2
(b) Pressure relative approximation error.

10 2 10 1 10 1 10 0 h k H k L 2 (⌃) k k L 2 (⌃) H = h/3 H = h H = h 1 2 h 1 2
(c) Lagrange multiplier relative approximation error. with ' 1 (x, y) def = x, ' 2 (x, y) 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:

2 h | R @⌦1 u h • n| H = h/3 H = h H = h
u 1 (x, y) = ✓ 20 (x + 1) 2 x 3 y (y 1) (2y 1) 10 (x + 1) x 2 y 2 (y 1) 2 (5x + 3) ◆ , p 1 (x, y) = x + 1, u 2 (x, y) = ✓ 10 (y 1) x 2 y 2 (x 1) 2 (5y 3) 20 (y 1) 2 y 3 x (x 1) (1 2x) ◆ , p 2 (x, y) = (70)
where the pressure profile has been chosen to provide a non-constant jump through the interface. The righthand side functions f , g and u ⌃ are chosen so that (1) holds with (70). In particular, zero velocity u ⌃ = 0 is prescribed on the interface. The expression for the Lagrange multiplier is obtained by plugging the exact solution (70) into (4). The fluid dynamic viscosity is set to µ = 0.035.

In order to numerically asses the accuracy of (10), we consider a family of triangulations T ⌦ h resulting from the uniform subdivision of D into 2N +1 segments and of N into N +2 segments, with N 2 {5, 10, 20, 40, 80, 160}.

The unfitted interfacial mesh T ⌃

H is obtained by uniformly subdividing the domain ⌃ into 3N segments. We

20

TITLE WILL BE SET BY THE PUBLISHER

10 2 10 1 10 1 10 0 h k"(u u h )k L 2 (⌦) k"(u)k L 2 (⌦) H = h/3 H = h H = h 1 2 h 1 2
(a) Velocity relative energy error.

10 2 10 1 10 2 10 1 h kp p h k L 2 (⌦) kpk L 2 (⌦) H = h/3 H = h H = h 1 2 h 1 2
(b) Pressure relative approximation error.

10 2 10 1 10 1 10 0 h k H k L 2 (⌃) k k L 2 (⌃) H = h/3 H = h H = h 1 2 h 1 2
(c) Lagrange multiplier relative approximation error. with ' 1 (x, y) def = x, ' 2 (x, y) 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:

2 h | R @⌦1 u h • n| H = h/3 H = h H = h
u 1 (x, y) = ✓ 20 (x + 1) 2 x 3 y (y 1) (2y 1) 10 (x + 1) x 2 y 2 (y 1) 2 (5x + 3) ◆ , p 1 (x, y) = x + 1, u 2 (x, y) = ✓ 10 (y 1) x 2 y 2 (x 1) 2 (5y 3) 20 (y 1) 2 y 3 x (x 1) (1 2x) ◆ , p 2 (x, y) = y, (70)
where the pressure profile has been chosen to provide a non-constant jump through the interface. The righthand side functions f , g and u ⌃ are chosen so that (1) holds with (70). In particular, zero velocity u ⌃ = 0 is prescribed on the interface. The expression for the Lagrange multiplier is obtained by plugging the exact solution (70) into (4). The fluid dynamic viscosity is set to µ = 0.035.

In order to numerically asses the accuracy of (10), we consider a family of triangulations T ⌦ h resulting from the uniform subdivision of D into 2N +1 segments and of N into N +2 segments, with N 2 {5, 10, 20, 40, 80, 160}.

The unfitted interfacial mesh T ⌃

H is obtained by uniformly subdividing the domain ⌃ into 3N segments. We (b) Pressure relative approximation error.

BE SET BY THE PUBLISHER

h = h 10 2 10 1 10 2 10 1 h kp p h k L 2 (⌦) kpk L 2 (⌦) H = h/3 H = h H = h 1 2 h 1 2
(b) Pressure relative approximation error.

2 h | R @⌦1 u h • n| H = h/3 H = h H = h 1 2 h 2 (d) Mass loss in ⌦1.
e velocity, pressure, Lagrange multiplier and mass loss d for ✓ = 1 and different ratios between h and H. eters (x l,i , x r,i , y b,i , y t,i ) are chosen as ( 1, 0, 0, 1) and (0, 1, 0, 1) analytical solution to problem (1) reads:

) 3) ◆ , u 2 (x, y) = ✓ 10 (y 1) x 2 y 2 (x 1) 2 (5y 3) 20 (y 1) 2 y 3 x (x 1) (1 2x) ◆ , p 2 (x, y) = y, (70)
rovide a non-constant jump through the interface. The rightthat (1) holds with (70). In particular, zero velocity u ⌃ = 0 or the Lagrange multiplier is obtained by plugging the exact ity is set to µ = 0.035. 10), we consider a family of triangulations T ⌦ h resulting from the and of N into N +2 segments, with N 2 {5, 10, 20, 40, 80, 160}. y uniformly subdividing the domain ⌃ into 3N segments. We

10 0 h k"(u u h )k L 2 (⌦) k"(u)k L 2 (⌦) H = h/3 H = h H = h 1 2 h 1 2
(a) Velocity relative energy error.

10 2 10 2 10 1 kp p h k L 2 (⌦) kpk L 2 (⌦) (b) Pressure r 10 2 10 1 10 1 10 0 h k H k L 2 (⌃) k k L 2 (⌃) H = h/3 H = h H = h 1 2 h 1 2
(c) Lagrange multiplier relative approximation error. = y and the parameters (x l,i , x r,i , y b,i , y t,i ) are for i = 1 and i = 2, respectively. The resulting analytical solution to proble

u 1 (x, y) = ✓ 20 (x + 1) 2 x 3 y (y 1) (2y 1) 10 (x + 1) x 2 y 2 (y 1) 2 (5x + 3) ◆ , p 1 (x, y) = x + 1, u 2 (x, y) = ✓ 10 (y 1 20 (y p 2 (x, y) = y,
where the pressure profile has been chosen to provide a non-constant jump hand side functions f , g and u ⌃ are chosen so that (1) holds with (70). is prescribed on the interface. The expression for the Lagrange multiplier solution (70) into (4). The fluid dynamic viscosity is set to µ = 0.035.

In order to numerically asses the accuracy of (10), we consider a family of uniform subdivision of D into 2N +1 segments and of N into N +2 segmen The unfitted interfacial mesh T ⌃ H is obtained by uniformly subdividing the where the pressure profile has been chosen to provide a non-constant jump through the interface. The righthand side functions f , g and u Σ are chosen so that (1) holds with (70). In particular, zero velocity u Σ = 0 is prescribed on the interface. The expression for the Lagrange multiplier λ is obtained by plugging the exact solution (70) into (4). The fluid dynamic viscosity is set to µ = 0.035.

In order to numerically asses the accuracy of (10), we consider a family of triangulations T Ω h resulting from the uniform subdivision of Γ D into 2N + 1 segments and of Γ N into N + 2 segments, with N ∈ {5, 10, 20, 40, 80, 160}.

The unfitted interfacial mesh T Σ

H is obtained by uniformly subdividing the domain Σ into 3N segments. We hence have H ≈ 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 and4b, 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 [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] 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 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 BH h , as well as the incompressibility constraint for the additional basis function (see Remark 3.2). The interface is discretized using a mesh T Σ h , 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 T Σ H . This piecewise constant normal will be denoted by n 0 Σ . We also consider a P 1 -approximation, denoted n 1 Σ , as a weighted average of n 0 Σ with respect to the length of the surrounding elements. (⌦)

h 3 ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (a)
h 3 ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (b)
Pressure relative approximation error. For the present method (✓ = 0 and ✓ = 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 ✓ = 0 or ✓ = 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 0 ⌃ of the normal is used, (⌦)

h k H k L 2 (⌃) k k L 2 (⌃) h ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 1 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (c)
3 h | R @⌦1 u h • n i | h 2 ✓ = 1, n 1 ⌃ FD-LM, n 1 ⌃ h ✓ = 1, n 0 ⌃ FD-LM, n 0 ⌃ (d)
h 3 ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (a)
h 3 ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (b)
Pressure relative approximation error. For the present method (✓ = 0 and ✓ = 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 ✓ = 0 or ✓ = 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 0 ⌃ of the normal is used, (⌦)

h k H k L 2 (⌃) k k L 2 (⌃) h ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 1 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (c)
3 h | R @⌦1 u h • n i | h 2 ✓ = 1, n 1 ⌃ FD-LM, n 1 ⌃ h ✓ = 1, n 0 ⌃ FD-LM, n 0 ⌃ (d)
h 3 ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (a)
h 3 ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (b)
Pressure relative approximation error. For the present method (✓ = 0 and ✓ = 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 ✓ = 0 or ✓ = 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 0 ⌃ of the normal is used, (c) Lagrange multiplier relative approximation error.

10 0 h k H k L 2 (⌃) k k L 2 (⌃) h ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 1 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (c)
3 h | R @⌦1 u h • n i | h 2 ✓ = 1, n 1 ⌃ FD-LM, n 1 ⌃ h ✓ = 1, n 0 ⌃ FD-LM, n 0 ⌃ (d)
h 3 ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (a)
h 3 ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (b)
Pressure relative approximation error. For the present method (✓ = 0 and ✓ = 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 ✓ = 0 or ✓ = 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 0 ⌃ of the normal is used, For the present method (θ = 0 and θ = 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 θ = 0 or θ = 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 0 Σ of the normal is used, while we reach 3 if n 1 Σ if used instead. This is a consequence of the fact that, assuming the geometry is exactly represented, the exact solution (u, p, λ) belongs to the discrete space V h × Q h × Λ H . The error observed is then due to geometry approximation only. More precisely, taking into account the discrete normal representation n i Σ TITLE WILL BE SET BY THE PUBLISHER

10 0 h k H k L 2 (⌃) k k L 2 (⌃) h ✓ = 1, n 1 ⌃ ✓ = 0, n 1 ⌃ FD-LM, n 1 ⌃ h 1 2 ✓ = 1, n 0 ⌃ ✓ = 0, n 0 ⌃ FD-LM, n 0 ⌃ (c)
3 h | R @⌦1 u h • n i | h 2 ✓ = 1, n 1 ⌃ FD-LM, n 1 ⌃ h ✓ = 1, n 0 ⌃ FD-LM, n 0 ⌃ (d)
h n ⌃ n i ⌃ L 2 (⌃) h n 0 ⌃ h 2 n 1 ⌃ Figure 8.
Convergence histories for the P 0 and P 1 interface normals.

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

µ 1 2 k" (u u h )k L 2 (⌦) + µ 1 2 kp p h k L 2 (⌦) + µ 1 2 h 1 2 k H k L 2 (⌃) . µ 1 2 kJpKk L 2 (⌃) (kn ⌃ ⇧ H n ⌃ k L 2 (⌃) + n i ⌃ ⇧ H n ⌃ L 2 (⌃) ).
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 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 1 2 convergence rate. In addition, we have also considered test case (71) with a straight interface ⌃ = {0} ⇥ (0, 1) (see Figure 3a). In this geometrical setting, n ⌃ = n 1 ⌃ = n 0 ⌃ , thus the geometry of the interface is exactly represented by the mesh T ⌃ H and n ⌃ belongs to ⇤ H . As expected by Theorem 4.12 the numerical solution corresponds to the analytical solution up to machine error.

Comparison with other fictitious domain methods

We now compare the behavior of the proposed unfitted mesh method [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] 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 e Q h of Section 3) and the interface Dirichlet constraint is treated in a collocated-penalized fashion (see, e.g., [?, Section 3.2.1]). This amounts to consider in (2) the following non-conforming approximation of the Lagrange multiplier space ⇤:

⇤ C H = ( ⇠ H = N H X i=1 ⇠ i x s i | ⇠ i 2 R 2 , i = 1, . . . , N H ) , where {x s i } N H
i=1 are the points of the interface mesh T ⌃ H and x s i denotes the Dirac's measure at point x s i (see, e.g., [?, ?, ?]). The discrete treatment of the Dirichlet constraint takes the penalized-collocated form

u h (x s i ) = u ⌃ (x s i ) + " i 8i = 1, . . . , N H , (72)
where " > 0 is a given (non-dimensionless) user-defined parameter. Note that the relation (72) enables the elimination of the Lagrange multiplier , so that the first alternative fictitious domain method reads: Find (i ∈ {0, 1}) and the actual value of the exact solution, estimate (69) can be modified as

(u h , e p h ) 2 V h ⇥ e Q h such that a FD h," ( (u h , e p h ) , (v h , e q h ) ) = `" (v h ) (73)
µ 1 2 ∥ε (u -u h )∥ L 2 (Ω) + µ -1 2 ∥p -p h ∥ L 2 (Ω) + µ -1 2 h 1 2 ∥λ -λ H ∥ L 2 (Σ) ≲ µ -1 2 ∥ p ∥ L 2 (Σ) (∥n Σ -Π H n Σ ∥ L 2 (Σ) + n i Σ -Π H n Σ L 2 (Σ) ).
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 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 1 2 convergence rate. In addition, we have also considered test case (71) with a straight interface Σ = {0} × (0, 1) (see Figure 3a). In this geometrical setting, n Σ = n 1 Σ = n 0 Σ , thus the geometry of the interface is exactly represented by the mesh T Σ H and n Σ belongs to Λ H . As expected by Theorem 4.12 the numerical solution corresponds to the analytical solution up to machine error.

Comparison with other fictitious domain methods

We now compare the behavior of the proposed unfitted mesh method [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] 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 Q h of Section 3) and the interface Dirichlet constraint is treated in a collocated-penalized fashion (see, e.g., [START_REF] Boilevin-Kayl | Numerical methods for immersed FSI with thin-walled structures[END_REF]Section 3.2.1]). This amounts to consider in (2) the following non-conforming approximation of the Lagrange multiplier space Λ:

Λ C H = ξ H = N H i=1 ξ i δ x s i | ξ i ∈ R 2 , i = 1, . . . , N H , where {x s i } N H
i=1 are the points of the interface mesh T Σ H and δ x s i denotes the Dirac's measure at point x s i (see, e.g., [START_REF] Santos | A partitioned fluid-structure algorithm for elastic thin valves with contact[END_REF][START_REF] Fabrèges | Approximation of single layer distributions by dirac masses in finite element computations[END_REF][START_REF] Glowinski | A distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF]). The discrete treatment of the Dirichlet constraint takes the penalized-collocated form

u h (x s i ) = u Σ (x s i ) + ελ i ∀i = 1, . . . , N H ,
where ε > 0 is a given (non-dimensionless) user-defined parameter. Note that the relation (72) enables the elimination of the Lagrange multiplier λ, so that the first alternative fictitious domain method reads: Find

(u h , p h ) ∈ V h × Q h such that a FD h,ε ( (u h , p h ) , (v h , q h ) ) = ℓ ε (v h ) (73) for all (v h , q h ) ∈ V h × Q h , with the notations a FD h,ε ( (u h , p h ) , (v h , q h ) ) def = a( (u h , p h ) , (v h , q h ) ) + s BP h ( p h , q h ) + d ε (u h , v h ) , ℓ ε (v h ) def = ℓ (v h ) + d ε (u Σ , v h ), d ε (u h , v h ) def = 1 ε N H i=1 u h (x s i ) • v h (x s i ) .
In the above definitions, the bilinear form a and the linear functional ℓ are those defining the weak formulation of the Stokes problem (1), in Section 2. The results of the method (73) have been obtained with ε = 10 -5 and γ 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 T e h of the interface Σ as follows (see, e.g., [START_REF] Casquero | NURBS-based numerical proxies for red blood cells and circulating tumor cells in microscale blood flow[END_REF][START_REF] Galvin | Stabilizing poor mass conservation in incompressible flow problems with large irrotational forcing and application to thermal convection[END_REF][START_REF] Kamensky | An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves[END_REF]): where γ d > 0 and 0 < γ gd ≪ 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).

s BP-GD ( (u h , p h ) , (v h , q h ) ) def = K∈T Ω h K δ h γ p h 2 µ ∇ p h • ∇ q h + K∈T Ω h K γ d µ δ h ∇ • u h ∇ • v h , δ h def = 1 K ∈ T Ω h \ T e h , γ gd K ∈ T e h , (
The resulting discrete problem reads as follows:

Find (u h , p h ) ∈ V h × Q h such that a FD-GD h,ε ( (u h , p h ) , (v h , q h ) ) = ℓ ε (v h ) (74) for all (v h , q h ) ∈ V h × Q h , where a FD-GD h,ε ( (u h , p h ) , (v h , q h ) ) def = a( (u h , p h ) , (v h , q h ) ) + s BP-GD h ( (u h , p h ) , (v h , q h ) ) + d ε (u h , v h ) .
The results of the method (74) have been obtained with γ d = 1, γ gd = 10 -4 , ε = 10 -5 and γ p = 10 -2 . We will refer to this method as FD-GD in the sequel. 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).

FD ε FD -GD ε FD-LM Present method θ = 1 θ = 0 13
Figures 9a and9c 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).

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. The Neumann boundary is given by Γ N = {1} × (0, 1) and a Poiseuille velocity profile u = -5y(1y)n is prescribed on {-1} × (0, 1). Additionally, a traction-free boundary condition g = 0 is imposed on {1} × (0, 1). We also consider a null right-hand side function f = 0 for the Stokes equation (1) and a velocity profile u Σ = 0 on the immersed interface Σ. The fluid dynamic viscosity is set to µ = 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 Σ fic that connects the two endpoints of Σ immersed in Ω. The sub-domain Ω 1 (resp. Ω 2 ) is then defined as the left (resp. right) side of the closed interface Σ ∪ Σ fic . It should be noted that, since no unknowns are attached to Σ 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 Ω 1 . Nevertheless, from a computational perspective, when replacing the mass conservation contraint on Ω 1 with its equivalent on ∂Ω 1 , as highlighted in Remark 3.2, it is essential to choose an appropriate discretization of Σ 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 ≈ 6.25 × 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 Σ is discretized with 120 segments (i.e., 60 on each side). The fictitious interface Σ fic is discretized with 25 segments.

The value of the user defined parameter ε is set to 10 -5 and γ p = 10 -2 . Furthermore, for the FD-GD method, we have γ d = 1 and γ gd = 10 -4 . For the FD-LM method, the stabilization parameter γ λ is set to 10 2 . In the present method, the normal is represented by n 1 Σ defined in Section 5.2.1 at the discrete level. Similar results are obtained when considering n 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 and12b-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 θ = 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 Ω 1 is slightly lower than the reference solution (compare Figures 13a and13e). For the non-symmetric variant (θ = 0), despite overall mass conservation in Ω 1 (see Table 2), Figures 11f and12f reveal fluid leakage through the interface and a significantly lower pressure in Ω 1 , as shown in Figure 13f. Taking γ λ larger in the non-symmetric variant (θ = 0) slightly improves the velocity results, but at the price of degrading the approximation of λ. A slight pressure discontinuity is noticeable across the fictitious interface Σ 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 Ω 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 2. Mass loss in Ω 1 .

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 fluidsolid 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 1 2ε 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 (θ = 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 (θ = 0) show that the exact satisfaction of

Figure 1 .

 1 Figure 1. Left: admissible geometric configuration. Right: non-admissible configuration (see Assumption 2.1 and Remark 2.2).

defined as the dual space of [H 1 2 00 2 00 2 00

 222 (Σ)] d , and the bilinear form b : Λ × [H 1 (Σ)] d → R, given by b(ξ, z) def = ⟨ξ, z⟩, where ⟨•, •⟩ denotes the duality pairing between Λ and [H 1 (Σ)] d . The weak formulation of problem (1) then reads as follows: Find

Figure 2 .

 2 Figure 2. In gray the fluid domain triangulation T Ω h and in black the interface domain triangulation T Σ H . The elements of T c h are marked in , the elements of T e h are marked in and the nodes in N c h are marked by .

1 2 (

 2 ⌃). Elements of such a proof are given in [3, Propositions A.1, A.2].

Figure 3 .

 3 Figure 3. Geometric configurations considered in the numerical tests.

1 2 (

 2 ⌃). Elements of such a proof are given in [3, Propositions A.1, A.2].

Figure 3 .

 3 Figure 3. Geometric configurations considered in the numerical tests.

1 2 (

 2 ⌃). Elements of such a proof are given in [3, Propositions A.1, A.2].

Figure 3 .

 3 Figure 3. Geometric configurations considered in the numerical tests.

Figure 3 .

 3 Figure 3. Geometric configurations considered in the numerical tests.

  a) Velocity relative energy error.

  Pressure relative approximation error.

  c) Lagrange multiplier relative approximation error.

  Mass loss in ⌦1. The results for ✓ = 0 are omitted since this variant guarantees interfacial mass conservation.

Figure 4 .

 4 Figure 4. Convergence histories of the velocity, pressure, Lagrange multiplier and mass loss in ⌦ 1 obtained with the present method (✓ = 0, ✓ = 1) and the FD-LM method.

  a) Velocity relative energy error.

  Pressure relative approximation error.

  c) Lagrange multiplier relative approximation error.

  Mass loss in ⌦1. The results for ✓ = 0 are omitted since this variant guarantees interfacial mass conservation.

Figure 4 .

 4 Figure 4. Convergence histories of the velocity, pressure, Lagrange multiplier and mass loss in ⌦ 1 obtained with the present method (✓ = 0, ✓ = 1) and the FD-LM method.

  a) Velocity relative energy error.

  Pressure relative approximation error.

  c) Lagrange multiplier relative approximation error.

  Mass loss in ⌦1. The results for ✓ = 0 are omitted since this variant guarantees interfacial mass conservation.

Figure 4 .

 4 Figure 4. Convergence histories of the velocity, pressure, Lagrange multiplier and mass loss in ⌦ 1 obtained with the present method (✓ = 0, ✓ = 1) and the FD-LM method.

  a) Velocity relative energy error.

  Pressure relative approximation error.

  c) Lagrange multiplier relative approximation error.

  Mass loss in ⌦1. The results for ✓ = 0 are omitted since this variant guarantees interfacial mass conservation.

Figure 4 .

 4 Figure 4. Convergence histories of the velocity, pressure, Lagrange multiplier and mass loss in ⌦ 1 obtained with the present method (✓ = 0, ✓ = 1) and the FD-LM method.

  d) Mass loss in Ω1. The results for θ = 0 are omitted since this variant guarantees interfacial mass vation.

Figure 4 .

 4 Figure 4. Convergence histories of the velocity, pressure, Lagrange multiplier mass loss in Ω 1 obtained with the present method (θ = 0, θ = 1) and the FD-LM method.

20 TITLE

 20

10

 10

2 (

 2 d) Mass loss in ⌦1.

Figure 5 .

 5 Figure 5. Convergence histories of the velocity, pressure, Lagrange multiplier and mass loss in ⌦ 1 obtained with the present method for ✓ = 1 and different ratios between h and H.

  (a) Velocity relative energy error.

10

 10

2 (

 2 d) Mass loss in ⌦1.

Figure 5 .

 5 Figure 5. Convergence histories of the velocity, pressure, Lagrange multiplier and mass loss in ⌦ 1 obtained with the present method for ✓ = 1 and different ratios between h and H.

10

 10

  (c) Lagrange multiplier relative approximation error.

20 TITLE

 20

Figure 5 .

 5 Figure 5. Convergence histories of the velocity, pressure, Lagran in ⌦ 1 obtained with the present method for ✓ = 1 and different rat

  (d) Mass loss in Ω1.

Figure 5 .

 5 Figure 5. Convergence histories of the pressure, Lagrange multiplier and mass loss in Ω 1 obtained with the present method for θ = 1 and different ratios between h and H.

Figure 6 .(c) n 1 Σ : P1 representationFigure 6 .

 616 Figure 6. Geometric approximation of the normal vectors to the interface ⌃. In gray the interface ⌃ and in black the mesh T ⌃ H . (c) n 1 Σ : P1 representation Figure 6. Geometric approximation of the normal vectors to the interface Σ. In gray the interface Σ and in black the mesh T Σ H .

22 TITLE

 22

  Velocity energy error.

  Lagrange multiplier relative approximation error.

10

 10

  Mass loss in ⌦1. The results for ✓ = 0 are omitted since this variant guarantees interfacial mass conservation.

Figure 7 .

 7 Figure 7. Convergence histories with the P 0 and P 1 approximation of the interface normal are compared among the present method in both the ✓ = 1 and ✓ = 0 variants and the un-enriched FD-LM method.

  Velocity energy error.

  Lagrange multiplier relative approximation error.

10

 10

  Mass loss in ⌦1. The results for ✓ = 0 are omitted since this variant guarantees interfacial mass conservation.

Figure 7 .

 7 Figure 7. Convergence histories with the P 0 and P 1 approximation of the interface normal are compared among the present method in both the ✓ = 1 and ✓ = 0 variants and the un-enriched FD-LM method.

  Velocity energy error.

  Lagrange multiplier relative approximation error.

10

 10

  Mass loss in ⌦1. The results for ✓ = 0 are omitted since this variant guarantees interfacial mass conservation.

Figure 7 .

 7 Figure 7. Convergence histories with the P 0 and P 1 approximation of the interface normal are compared among the present method in both the ✓ = 1 and ✓ = 0 variants and the un-enriched FD-LM method.

  Velocity energy error.

  Lagrange multiplier relative approximation error.

10

 10

  Mass loss in ⌦1. The results for ✓ = 0 are omitted since this variant guarantees interfacial mass conservation.

Figure 7 .

 7 Figure 7. Convergence histories with the P 0 and P 1 approximation of the interface normal are compared among the present method in both the ✓ = 1 and ✓ = 0 variants and the un-enriched FD-LM method.

  (d) Mass loss in Ω1. The results for θ = 0 are omitted since this variant guarantees interfacial mass conservation.

Figure 7 .

 7 Figure 7. Convergence histories with the P 0 and P 1 approximation of the interface normal are compared among the present method in both the θ = 1 and θ = 0 variants and the un-enriched FD-LM method.

Figure 8 .

 8 Figure 8. Convergence histories for the P 0 and P 1 interface normals.

  (a) FD method. (b) FD-GD method. (c) FD-LM method. (d) Present method: θ = 1. (e) Present method: θ = 0.

Figure 9 .

 9 Figure 9. Snapshots of the velocity magnitude.

  a) FD method. (b) FD-GD method.(c) FD-LM method.(d) Present method: θ = 1.(e) Present method: θ = 0.

Figure 10 .

 10 Figure 10. Snapshots of the pressure.

5 . 4 • 10 -4 17 3 . 4 •

 541034 10 -10 2.0 • 10 -14

  (a) Fitted mesh -Reference (b) FD method. (c) FD-GD method. (d) FD-LM method.(e) Present method: θ = 1.(f) Present method: θ = 0.

Figure 11 .

 11 Figure 11. Snapshots of the velocity magnitude.

( a )

 a Fitted mesh -Reference (b) FD method.(c) FD-GD method.(d) FD-LM method.(e) Present method: θ = 1.(f) Present method: θ = 0.

Figure 12 .

 12 Figure 12. Snapshots of the velocity.

  TITLE WILL BE SET BY THE PUBLISHER (a) Fitted mesh -Reference (b) FD method. (c) FD-GD method. (d) FD-LM method.(e) Present method: θ = 1.(f) Present method: θ = 0.

Figure 13 .

 13 Figure 13. Snapshots of the pressure.

the set of nodes belonging to at least one cut element and by T e

	h one node with a cut element. See Figure 2. Let us now introduce the following assumptions. def = {K ∈ T Ω h : {x K,i } ∩ N c h ̸ = ∅} the set of elements of T Ω h sharing at least 8 TITLE WILL BE SET BY THE PUBLISHER
	Figure 2. In gray the fluid domain triangulation T ⌦ h and in black the interface domain trian-gulation T ⌃ H . The elements of T c h are marked in , the elements of T e h are marked in and the nodes in N c h are marked by .
	with
	b ⇧ h

  b) Pressure relative approximation error.

			TITLE WILL BE SET BY THE PUBLISHER	19
	10 1	10 2	10 1

Table 1 .

 1 Mass loss in Ω 1 .

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 scalaire 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 presenté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.

* Dand M. Vidrascu were partially supported by the French National Research Agency (ANR), through the SIMR project (ANR-19-CE45-0020).

enhanced mass conservation provided by the pressent method with respect to the FD-LM method (see Remark 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 θ = 1 method across various h/H ratios. The superior regularity of λ, with respect to the one assumed in Theorem 4.12, delivers a better convergence rate in terms of H. Indeed, in this case, the estimate (60) scales as H + h 1 2 -ε , so that we can take a coarser Lagrange multiplier mesh, H = h 1 2 , while preserving the overall convergence rate.

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., [START_REF] Groß | An extended pressure finite element space for two-phase incompressible flows with surface tension[END_REF][START_REF] Kamensky | An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves[END_REF]). 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 Σ = {(0.2 sin(πt), t)| t ∈ [0, 1]}. A traction force g = -3 • 10 5 n is imposed at the channel inlet {-1} × (0, 1), while no traction g = 0 is prescribed at the outlet {1} × (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 Σ is set to be u Σ = 0. The fluid dynamic viscosity is set to µ = 10. The solution to problem (1) reads

Convergence study

hence have H ⇡ 3 h. The associated convergence histories for the L -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 and4b, 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 [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] 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 pressent method with respect to the FD-LM method (see Remark 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 ✓ = 1 method across various h/H ratios. The superior regularity of , with respect to the one assumed in Theorem 4.12, delivers a better convergence rate in terms of H. Indeed, in this case, the estimate (60) scales as H + h 1 2 " , so that we can take a coarser Lagrange multiplier mesh, H = h 1 2 , while preserving the overall convergence rate.

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., [START_REF] Groß | An extended pressure finite element space for two-phase incompressible flows with surface tension[END_REF][START_REF] Kamensky | An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves[END_REF]). 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 ⌃ = {(0.2 sin(⇡t), t)| t 2 [0, 1]}. A traction force g = 3 • 10 5 n is imposed at the channel inlet { 1} ⇥ (0, 1), while no traction g = 0 is prescribed at the outlet {1} ⇥ (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 ⌃ is set to be u ⌃ = 0. The fluid dynamic viscosity is set to µ = 10. The solution to problem (1) reads (a) nΣ : normal vector to Σ.

hence have H ⇡ 3 h. The associated convergence histories for the L -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 and4b, 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 [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] 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 pressent method with respect to the FD-LM method (see Remark 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 ✓ = 1 method across various h/H ratios. The superior regularity of , with respect to the one assumed in Theorem 4.12, delivers a better convergence rate in terms of H. Indeed, in this case, the estimate (60) scales as H + h 1 2 " , so that we can take a coarser Lagrange multiplier mesh, H = h 1 2 , while preserving the overall convergence rate.

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., [START_REF] Groß | An extended pressure finite element space for two-phase incompressible flows with surface tension[END_REF][START_REF] Kamensky | An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves[END_REF]). 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 ⌃ = {(0.2 sin(⇡t), t)| t 2 [0, 1]}. A traction force g = 3 • 10 5 n is imposed at the channel inlet { 1} ⇥ (0, 1), while no traction g = 0 is prescribed at the outlet {1} ⇥ (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 ⌃ is set to be u ⌃ = 0. The fluid dynamic viscosity is set to µ = 10. The solution to problem (1) reads hence have H ⇡ 3 h. The associated convergence histories for the L -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 and4b, 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 [START_REF] Boffi | Local mass conservation of Stokes finite elements[END_REF] 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 pressent method with respect to the FD-LM method (see Remark 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 ✓ = 1 method across various h/H ratios. The superior regularity of , with respect to the one assumed in Theorem 4.12, delivers a better convergence rate in terms of H. Indeed, in this case, the estimate (60) scales as H + h 1 2 " , so that we can take a coarser Lagrange multiplier mesh, H = h 1 2 , while preserving the overall convergence rate.

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., [START_REF] Groß | An extended pressure finite element space for two-phase incompressible flows with surface tension[END_REF][START_REF] Kamensky | An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves[END_REF]). 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 ⌃ = {(0.2 sin(⇡t), t)| t 2 [0, 1]}. A traction force g = 3 • 10 5 n is imposed at the channel inlet { 1} ⇥ (0, 1), while no traction g = 0 is prescribed at the outlet {1} ⇥ (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 ⌃ is set to be u ⌃ = 0. The fluid dynamic viscosity is set to µ = 10.