Michel Duprez
email: michel.duprez@inria.fr

Vanessa Lleras
email: vanessa.lleras@umontpellier.fr

Alexei Lozinski
email: alexei.lozinski@univ-fcomte.fr

ϕ-FEM: an optimally convergent and easily implementable immersed boundary method for particulate flows and Stokes equations *

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

Introduction

The main goal of the present article is to demonstrate that the recently proposed ϕ-FEM methodology [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF][START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF][START_REF] Cotin | ϕ-fem: an efficient simulation tool using simple meshes for problems in structure mechanics and heat transfer[END_REF] is suitable for numerical simulation of incompressible viscous fluid flow past moving rigid bodies. This approach allows us to use simple (Cartesian) computational meshes, not evolving in time and not fitted to the moving rigid bodies, while achieving the optimal accuracy with classical finite element (FE) spaces of any order and performing the usual numerical integration on the whole mesh cells and facets, allowing for the use of standard FEM libraries for the implementation. We consider here only the creeping motion regime (zero Reynolds number), neglecting all the inertial terms in the equation governing both the fluid and the rigid bodies.

Numerical simulations of flows around moving rigid or elastic structures using immobile simple grids is a popular approach in, for instance, biomechanics, starting from the work of Peskin [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF]. Different approaches have emerged since then, such as the Immersed Boundary method [START_REF] Lai | An immersed boundary method with formal second-order accuracy and reduced numerical viscosity[END_REF][START_REF] Mittal | Immersed boundary methods[END_REF], the Fictitious Domain method [START_REF] Glowinski | A distributed Lagrange multiplier/fictitious domain method for particulate flows[END_REF][START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF], the penalty approximation [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF], etc. All these classical methods suffer from poor accuracy due to the necessity to approximate the singularities near the fluid-solid interfaces which arise as the artifact of extending the fluid velocity field inside the solid domain. More recently, several optimally convergent fictitious domain-type methods have been proposed for the Stokes equations, which can also be used to simulate the fluid-solid motions.

• Extend the governing equations from Ω to Ω h and write down a formal variational formulation on Ω h without taking into account the boundary conditions on Γ (the relevant part of the boundary of Ω).

• Impose the boundary conditions on Γ using an appropriate ansatz or some additional variables, explicitly involving the level set ϕ which provides the link to the actual boundary.

For instance, the homogeneous Dirichlet boundary conditions (u = 0 on Γ) can be imposed by the ansatz u = ϕw thus reformulating the problem in terms of the new unknown w.

• Add appropriate stabilization, typically combining the ghost penalty [START_REF] Burman | Ghost penalty[END_REF] with a least square imposition of the governing equation on the mesh cells intersected by Γ, to guarantee coerciveness/stability on the discrete level.

This program has been successfully carried out for elliptic scalar PDEs with Dirichlet boundary conditions in [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF] and for Neumann boundary conditions in [START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF]. Its feasibility is also demonstrated in [START_REF] Cotin | ϕ-fem: an efficient simulation tool using simple meshes for problems in structure mechanics and heat transfer[END_REF] for the linear elasticity with mixed boundary conditions including the cases of internal interfaces between different materials or cracks, and for the heat equation. However, the adaptation to the equations governing the fluid flow around the moving particles is not straightforward. In particular, the following challenges are dealt with in the present article:

• The discrete inf-sup stability theory should be adapted to the case of a non-standard variational formulation of the Stokes equations posed on Ω h rather than on Ω, and lacking the saddle-point structure. We shall show that this is possible by adapting the ghost penalty, which should be taken slightly more complicated than in the case of scalar elliptic equations [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF]. We shall do it here for Taylor-Hood finite elements of any order, but similar ideas should be also applicable to other classical inf-sup stable FE spaces.

• The motion equations for the solid particles involve the forces exerted on them by the surrounding fluid. These are defined through the integrals of some functions of fluid velocity and pressure on the particle boundary. However, the whole point of ϕ-FEM is to avoid such integrals. Indeed, the particle boundary is not resolved by the mesh, and our goal is to provide a method that necessitates the integration on the whole mesh cells or facets only. The way out of this paradox, pursued in the present paper, lies in providing a weak formulation of the governing equations, extended to Ω h , that incorporates in an appropriate way the force balance equations, without stating them directly. This formulation is similar in spirit but different from that in [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF].

We note that the method of this article shares some similarities to the shifted boundary method (SBM) proposed in [START_REF] Main | The shifted boundary method for embedded domain computations. Part I: Poisson and Stokes problems[END_REF] and analysed in [START_REF] Atallah | Analysis of the shifted boundary method for the Stokes problem[END_REF] in the case of Stokes equations. In particular, SBM also gives an optimal accurate solution (at least with the lowest order finite elements) without introducing integrals on the cut cells. It is however not evident how one can deal with the computation of the forces on the particles in the SBM framework. The paper is structured as follows. In the next section, we properly introduce the governing equations, develop an appropriate weak formulation, discretize it (thus introducing our ϕ-FEM scheme), and announce the main theorem about the convergence of the scheme. Section 3 is devoted to the proof of this theorem. As a by-product, we also introduce a ϕ-FEM approach to discretize the Stokes equations alone (on a fixed geometry) on a non-fitted mesh. The details about this (comparatively simple) particular case are given in Appendix A. In section 4, we illustrate our theoretical results with numerical examples both for the Stokes equations and for the fluid/rigid particle motion problem. ϕ-FEM is also compared there with a standard (non isoparametric) FEM on fitted meshes, demonstrating the superiority of ϕ-FEM in terms of the accuracy achieved on comparable meshes. We end up with the last section giving some conclusions and perspectives.

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

2 Construction of the ϕ-FEM scheme and main results

Governing equations

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

-2ν div D(u) + ∇p = ρ f g, in Ω (1a) div u = 0, in Ω (1b) u = U + ψ × r, on Γ (1c) u = 0, on Γ w (1d) Γ (2νD(u) -pI)n = mg (1e) Γ (2νD(u) -pI)n × r = 0 (1f) Ω p = 0 (1g)
Here, the unknowns are the fluid velocity u : Ω → R d and the pressure p : Ω → R, the velocity of the particle barycenter U ∈ R d , and the angular velocity of the particle

ψ ∈ R d ′ (d ′ = 1 if d = 2 and d ′ = 3 if d = 3
). In these equations, D(u) =1 2 (∇u + ∇u T ) denotes the strain tensor, r denotes the vector from the barycenter of the solid S, and n denotes the unit normal on Γ looking into the solid. Equations (1e)-(1f) come from the balance of forces exerted on the particle (the force exerted by the fluid and the gravitational force).

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

A formal derivation of the appropriate weak formulation

Let T O h be a regular simplicial mesh on O (the background mesh). Assume that the solid and fluid domains are given by the level-set function ϕ: S = {ϕ > 0} and Ω = O ∩{ϕ < 0}. Introduce the active computational mesh T h as a submesh of T O h covering Ω, i.e. excluding the cells of T O h lying completely inside S. Let Ω h ⊃ Ω be the domain of T h and G h be the component of ∂Ω h , other than Γ w , and thus lying inside S, cf. Fig. 1. 1 Assume (on a formal level, just to derive the scheme) that u and p can be extended from Ω

to Ω h as solution to the Stokes equations so that -2ν div D(u) + ∇p = ρ f g and div u = 0 in Ω h .

Taking any sufficiently smooth test functions v and q on Ω h such that v = 0 on Γ w , an integration by parts gives 2ν

Ω h D(u) : D(v) - Ω h p div v - Ω h q div u - G h (2νD(u) -pI)n • v = Ω h ρ f g • v. (2)
Assuming u = 0 on Γ w , this imposes already the boundary condition (1d) on Γ w , which we suppose to fit to the mesh T h . On the contrary, this formulation does not take into account any boundary conditions on Γ. In order to incorporate boundary conditions (1c) we make the ansatz

u = ϕw + χ(U + ψ × r) ( 3
)
where ϕ is the level-set for Ω so that ϕ = 0 on Γ, and χ is a sufficiently smooth function on O such that χ = 1 on the solid S and χ = 0 on Γ w . This introduces the new vector valued unknown w on Ω h that should vanish on Γ w (indeed w = 0 on Γ w implies u = 0 on Γ w thanks to the choice of χ; in fact, the reason for χ is to decouple the boundary conditions on Γ and Γ w from one another).

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

v = ϕs + χ(V + ω × r) (4)
for all vector-valued functions s on Ω h vanishing on Γ w , and

V ∈ R d , ω ∈ R d ′ .
In particular, the test functions of the form χ(V + ω × r) can be used to take into account the force balance (1e-1f). To this end, we introduce B h = Ω h \ Ω, i.e. the strip between Γ and G h , and use the divergence theorem on B h to transfer the boundary term in [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] from G h to Γ where it can be evaluated by (1e-1f):

G h (2νD(u) -pI)n • χ(V + ω × r) = G h (2νD(u) -pI)n • (V + ω × r) = Γ (2νD(u) -pI)n • (V + ω × r) + B h div(2νD(u) -pI)n • (V + ω × r) = mg • V - B h ρ f g • χ(V + ω × r) = • • • (5)
(the unit normal n on G h in the first line is exterior with respect to domain Ω h , whereas n on Γ in the second line is the exterior unit normal with respect to domain Ω, so that the exterior normals with respect to B h are n on G h and -n on Γ). We now remark B h = Ω h \ (O \ S) to rewrite the above as

• • • = - Ω h ρ f g • χ(V + ω × r) + O ρ f g • χ(V + ω × r) - S ρ f g • (V + ω × r) + mg • V = - Ω h ρ f g • χ(V + ω × r) + O ρ f g • χ(V + ω × r) + 1 - ρ f ρ s mg • V. (6)
The last line is justified by observing S ρ s g • (V + ω × r) = mg • V with ρ s being the constant density of the solid. Indeed, S ρ s = m and S ρ s r = 0 since r = x -x b is the vector pointing from the barycenter of the solid x b = 1 m S ρ s x to the current position x.

Substituting the ansatzes (3)-( 4) for u and v into (2) and rewriting the boundary term using ( 5)-( 6) we arrive at the following formal variational formulation of our problem in terms of the new unknowns w, U , ψ:

find w : Ω h → R d vanishing on Γ w , U ∈ R d , ψ ∈ R d ′ , and p : Ω h → R such that 2ν Ω h D(ϕw + χ(U + ψ × r)) : D(ϕs + χ(V + ω × r)) - G h (2νD(ϕw + χ(U + ψ × r)) -pI)n • ϕs - Ω h p div(ϕs + χ(V + ω × r)) - Ω h q div(ϕw + χ(U + ψ × r)) = Ω h ρ f g • ϕs + O ρ f g • χ(V + ω × r) + 1 - ρ f ρ s mg • V (7) for all s : Ω h → R d vanishing on Γ w , V ∈ R d , ω ∈ R d ′ ,
and q : Ω h → R. In addition, the pressure p should satisfy the constraint (1g). Note that the formulation above contains only the integrals on Ω h , G h , O which can be easily approximated by quadrature rules on meshes T h and T O h . We can thus discretize using the usual finite elements for the trial and test functions.

The ϕ-FEM scheme: discretization with Taylor-Hood finite elements

We fix an integer k ⩾ 2 and introduce the approximations ϕ h and χ h to the levelset ϕ and to the cut-off χ, given by the standard nodal interpolation to the continuous FE spaces of degree k on the mesh T O h . The active computational mesh T h , its domain Ω h and the internal boundary component G h are actually defined as follows, cf. Fig. 1,

T h = {T ∈ T O h : T ∩ {ϕ h < 0} ̸ = ∅}, Ω h := (∪ T ∈T h T ) o , (8)
G h = ∂Ω h \ Γ w = {E (boundary facets of T h ) such that ϕ h ≥ 0 on E} . (9)
Moreover, we shall need the collections of the mesh cells T Γ h and facets F Γ h near the boundary Γ, as illustrated in Fig. 2, to include the appropriate stabilization into the FE scheme. More specifically, we introduce the submesh T Γ h ⊂ T h and the corresponding subdomain Ω Γ h ⊂ Ω containing the mesh elements intersected by the approximate interface

Γ h = {ϕ h = 0}, i.e. T Γ h = {T ∈ T h : T ∩ Γ h ̸ = ∅}, Ω Γ h := ∪ T ∈T Γ h T o . (10)
Finally, we set F Γ h as the collection of the interior facets of the mesh T h either cut by Γ h or belonging to a cut mesh element

F Γ h = {E (an internal facet of T h ) such that ∃T ∈ T Γ h and E ∈ ∂T }.
Remark 1. The definitions of T h and T Γ h above assume an idealized setting where one can check the sign of ϕ h at any point of any given mesh cell T . In practice, one would rather check this sign only at the vertices of the mesh or, eventually, at some other well chosen points. In our current implementation, we attribute the cells to T h or T Γ h based on the sign of ϕ (equivalently, the sign of ϕ h ) at the vertices only: in 2D, a triangle T ∈ T O h is selected to be in T h if ϕ ⩽ 0 on at least one vertex of T ; a triangle T ∈ T h is then selected to be in T Γ h if ϕ ⩾ 0 on at least one vertex of T . This deviation from definitions (8)-( 9) is not covered by our theory.

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

Introduce the FE spaces for velocity and pressure on the mesh T h :

V h = v h ∈ C( Ω) d : v h | T ∈ P k (T ) d ∀T ∈ T h , v h = 0 on Γ w and M h = q h ∈ C( Ω) : q h | T ∈ P k-1 (T ) ∀T ∈ T h , Ω q h = 0 .
Remark 2. Note that the definition of the pressure space involves an integral on Ω, which is incompatible with our ϕ-FEM framework since its whole point is to avoid integrals on Ω and Γ. In practice, we shall rather impose Ω h q h = 0, introducing a mismatch in the additive pressure constant (which, anyway, has no physical meaning) with respect to the exact solution satisfying (1g). We prefer however to keep the unimplementable constraint in the definition above to avoid some technical difficulties in theory. In practice, a special care will have to be taken in the interpretation of the error in pressure. We shall return to this technical point in the numerical results section.

The stabilized scheme inspired by [START_REF] Burman | Ghost penalty[END_REF] can be now written as: find

w h ∈ V h , U h ∈ R d , ψ h ∈ R d ′ , p h ∈ M h such that 2ν Ω h D(ϕ h w h + χ h (U h + ψ h × r)) : D(ϕ h s h + χ h (V h + ω h × r)) - G h (2νD(χ h (U h + ψ h × r) + ϕ h w h ) -p h I)n • ϕ h s h - Ω h p h div(ϕ h s h + χ h (V h + ω h × r)) - Ω h q h div(ϕ h w h + χ h (U h + ψ h × r)) + σ u J u (χ h (U h + ψ h × r) + ϕ h w h , χ h (V h + ω h × r) + ϕ h s h ) + σh 2 T ∈T Γ h T (-ν∆(ϕ h w h + χ h (U h + ψ h × r)) + ∇p h ) • (-ν∆(ϕ h s h + χ h (V h + ω h × r)) -∇q h ) + σ T ∈T Γ h T div(ϕ h w h + χ h (U h + ψ h × r)) div(ϕ h s h + χ h (V h + ω h × r)) = Ω h ρ f g • ϕ h s h + O ρ f g • χ h (V h + ω h × r) + 1 - ρ f ρ s mg • V h + σh 2 T ∈T Γ h T ρ f g • (-ν∆(ϕ h s h + χ h (V h + ω h × r)) -∇q h ) (11) for all s h ∈ V h , V h ∈ R d , ω h ∈ R d ′ , q h ∈ M h .
Here J u is the ghost penalties for the velocity, cf. [START_REF] Burman | Ghost penalty[END_REF]:

J u (u, v) = h E∈F Γ h E [∂ n u] • [∂ n v] + h 3 E∈F Γ h E ∂ 2 n u • ∂ 2 n v .
Note that, unlike [START_REF] Burman | Ghost penalty[END_REF][START_REF] Burman | Fictitious domain methods using cut elements: III. A stabilized Nitsche method for Stokes' problem[END_REF], we do not penalize the jumps of all the derivatives of the velocity; only the derivatives of order up to 2 are included in J u . There is no penalization on the pressure either. This alleviation of the ghost penalty is possible thanks to the additional least-squarestype stabilization (the terms multiplied by σ), cf. Lemmas 1 and 2. These least-squares terms are also necessary in themselves to control the fictitious extension of the solution outside Ω, cf. the proof of Lemma 9. Note that this extension is not present in CutFEM (this is indeed the principal difference between CutFEM and ϕ-FEM). We also mention that the version of the ghost penalty in ϕ-FEM for Poisson problem in [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF] is even more reduced: only the jumps of the first order derivatives are penalized there. The inclusion of the second order derivatives in J u in the present case of Stokes equations allows us to control both velocity and pressure in the forthcoming proofs, cf. Lemma 1.

Assumptions on the mesh and main results

Prior to stating our main results on the numerical convergence of our method, we begin with some geometrical assumptions on Ω and the functions ϕ and χ. • Each Π k is composed of a mesh element T k lying inside Ω and some elements cut by Γ, more precisely andΠ Γ k contains at most N mesh elements;

Π k = T k ∪ Π Γ k where T k ∈ T h , T k ⊂ Ω, Π Γ k ⊂ T Γ h ,
• Each mesh element in a patch Π k shares at least a facet with another mesh element in the same patch. In particular, T k shares a facet

F k with an element in Π Γ k ; • T Γ h = ∪ NΠ k=1 Π Γ k ; • Π k and Π l are disjoint if k ̸ = l.
Assumption 5. Any mesh cell T ∈ T h has at least d facets not lying on Γ w . Remark 3. Assumptions 1, 3, 4 are similar to those made in the previous ϕ-FEM publications [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF][START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF], which contain a more detailed discussion about them and some illustrations. In what concerns the mesh, these assumptions are satisfied if the mesh is sufficiently refined and Γ is sufficiently smooth. In what concerning the level-set function ϕ, we require essentially that it behaves like the signed distance to Γ near Γ and it is bounded away from zero far from Γ, while remaining globally smooth. In general, one cannot thus take the signed distance to Γ as ϕ everywhere since it is guaranteed to be smooth only in a vicinity of Γ.

Remark 4. The last assumption 5 is usually required in the theoretical analysis of Taylor-Hood elements for the Stokes equation in the geometrically conforming setting [START_REF] Ern | Theory and practice of finite elements[END_REF], although it can be significantly relaxed, at least in the 2D setting [START_REF] Boffi | Finite elements for the Stokes problem[END_REF]. Note that this assumption only affects the mesh near the outer wall Γ w , more particularly in the corners of O, which we treat in the standard geometrically conforming manner anyway. It does not impose any further restriction on the active mesh T h near the interface Γ, where ϕ-FEM is effectively employed.

Let us now state our main results: Theorem 1. Suppose that Assumptions 1-5 hold true and the mesh [START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF]. Denoting

T h is quasi-uniform. Let (u, U, ψ, p) ∈ H k+1 (Ω) d × R d × R d ′ × H k (Ω) be the solution to (1a)-(1g) and (w h , U h , ψ h , p h ) ∈ V h × R d × R d ′ × M h be the solution to
u h := χ h (U h + ψ h × r) + ϕ h w h it holds for h ≤ h 0 |u -u h | 1,Ω∩Ω h + 1 ν ∥p -p h ∥ 0,Ω∩Ω h ≤ Ch k (∥u∥ k+1,Ω + 1 ν ∥p∥ k,Ω ) (12)
and

|U -U h | + |ψ -ψ h | ≤ Ch k (∥u∥ k+1,Ω + 1 ν ∥p∥ k,Ω ) (13)
with some C > 0 and h 0 > 0 depending on the parameters C 0 , m, M , N in Assumptions 1-5, on the maximum of the derivatives of ϕ and χ of order up to k + 1, on the mesh regularity, and on the polynomial degree k, but independent of h, f , and u.

Moreover, supposing Ω ⊂ Ω h ∥u -u h ∥ 0,Ω ≤ Ch k+1/2 (∥u∥ k+1,Ω + 1 ν ∥p∥ k,Ω ) (14)
with a constant C > 0 of the same type as above.

Remark 5. The numerical results in Section 4 suggest that the convergence order for the particle velocity error and that for the L 2 -error of the fluid velocity is k + 1. This suggests that both estimates ( 13) and ( 14) are not sharp. This is similar to our previous studies [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF] and [START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF], in what concerns the L 2 -error.

Proofs

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

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

Some technical lemmas.

Lemmas 1, 2, and 4 are adaptions of, respectively, Lemmas 3.2, 3.3, and 3.4-3.5 from [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF]. Lemma 3 is a version of the well known Korn inequality stating the uniformity of the constant in this inequality for a family of h-dependent domains Ω h . All these results, most notably Lemma 2, are necessary to prove the inf-sup stability of our scheme.

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

v = ∂v ∂n = ∂ 2 v ∂n 2 = 0 on E (15)
and -∆v + ∇q = 0 and div v = 0 on T.

Then v = 0 and q = const on T .

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

v k = i,j≥0 v k ij x i y j and q = i,j≥0
q ij x i y j .

We will prove by strong induction on m that

v k ij = 0 ∀k ∈ {1, 2}, i ≥ 0, j ∈ {0, ..., m}. (S m ) Using (15), it holds i≥0 v k i0 x i = i≥0 v k i1 x i = i≥0 v k i2 x i = 0 for all x and k ∈ {1, 2}. Hence v k i0 = v k i1 = v k i2 = 0 for all i ≥ 0, k ∈ {1, 2}
and we obtain (S 2 ). Assume that for a given m ≥ 2, (S m ) holds true. Thanks to ( 16), one has

-(i + 2)(i + 1)v 1 (i+2)j -(j + 2)(j + 1)v 1 i(j+2) + (i + 1)q (i+1)j = 0, (17)
-(i + 2)(i + 1)v 2 (i+2)j -(j + 2)(j + 1)v 2 i(j+2) + (j + 1)q i(j+1) = 0 (18)
and

(i + 1)v 1 (i+1)j + (j + 1)v 2 i(j+1) = 0. ( 19
)
From ( 17) and ( 18), for each i, j

(i + 2)v 1 (i+2)(j+1) + (j + 3)(j + 2) i + 1 v 1 i(j+3) = (i + 3)(i + 2) j + 1 v 2 (i+3)j + (j + 2)v 2 (i+1)(j+2) .
The last equality for j = m -2 combined with (S m ) implies that

v 1 i(m+1) = 0 ∀i ≥ 0.
Relation [START_REF] Main | The shifted boundary method for embedded domain computations. Part I: Poisson and Stokes problems[END_REF] for j = m and (S m ) gives

v 2 i(m+1) = 0 ∀i ≥ 0,
which leads to (S m+1 ). Thus v = 0. This also implies ∇q = 0 on T thanks to [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF].

Lemma 2. For any β > 0 and any integers s, r ≥ 1 there exists 0 < α < 1 depending only on the mesh regularity and s, r such that for any continuous vector-valued P s FE function v h on T h and any continuous scalar P r FE function q h it holds

∥D(v h )∥ 2 0,Ω Γ h + (1 -α)h 2 |q h | 2 1,Ω Γ h ≤ α∥D(v h )∥ 2 0,Ω h + β   h 2 ∥ -∆v h + ∇q h ∥ 2 0,Ω Γ h + ∥ div v h ∥ 2 0,Ω Γ h + E∈F Γ h (h∥[∂ n v h ]∥ 2 0,E + h 3 ∥[∂ 2 n v h ]∥ 2 0,E )   . ( 20
)
Proof. Thanks to Assumption 4, the boundary Γ can be covered by patches {Π k } k=1,...,NΠ . Take β > 0 and set α := max

Π k ,v h ,q h F (Π k , v h , q h ), (21)
where

F (Π k , v h , q h ) = ∥D(v h )∥ 2 0,Π Γ k +h 2 |q h | 2 1,Π Γ k -βG(Π k , v h , q h ) ∥D(v h )∥ 2 0,Π k +h 2 |q h | 2 1,Π Γ k with G(Π k , v h , q h ) = h 2 ∥ -∆v h + ∇q h ∥ 2 0,Π Γ k + ∥ div v h ∥ 2 0,Π Γ k + h ∥[∂ n v h ]∥ 2 F k + h 3 ∂ 2 n v h 2 F k .
The maximum in ( 21) is taken over all the continuous vector-valued P s FE functions v h on Π k , all the continuous scalar P r FE functions q h on Π Γ k , such that the denominator in the expression for F does not vanish, and over all the possible configurations of patches Π k satisfying Assumption 4. The notation F k stands for the set of mesh facets inside the patch Π k which includes thus

F k separating T k from Π Γ k and the other facets inside Π Γ k . The norm ∥ • ∥ F k should be understood as ( F ∈F k ∥ • ∥ 2 F ) 1/2 . Since the maximized function F is invariant with respect to the transformation x → 1 h x, v h → 1 h v h , q h → q h , we can assume that h = 1 in (21). Furthermore F (Π k , v h , q h ) = F (Π k , λv h , λq h )
for any λ ̸ = 0. Hence the maximum ( 21) is attained since it can be taken over all admissible patches with h = 1 and all v h , q h such that

∥D(v h )∥ 2 0,Π k +h 2 |q h | 2 1,Π Γ k = 1
, forming the unit sphere in the finite dimensional space of all (v h , q h ) factored by rigid body motions on Π k and constants on Π Γ k . Clearly α ⩽ 1. Let us prove by contradiction that α < 1. Assume that α = 1. Consider the patch Π k (with h = 1) and v h ,

q h with ∥D(v h )∥ 2 0,Π k +|q h | 2 1,Π Γ k = 1 on which the maximum (21) is attained. Then, ∥D(v h )∥ 2 0,T k + βG(Π k , v h , q h ) = 0, since Π k = T k ∪ Π Γ k . We deduce that v h = ∂v h ∂n = ∂ 2 v h ∂n 2 = 0 on all the facets in F k and -∆v h + ∇q h = 0, div(v h ) = 0 on all T ∈ Π Γ k . Moreover, v h is a rigid body motion on T k . Let v rbm h be the rigid body motion velocity on R d coinciding with v h on T k . Thanks to Lemma 1 applied to v h -v rbm h and to q h on the cells in Π Γ k starting from the cell adjacent in Π Γ k to T k , we have v h = v rbm h on Π k and q h = const on Π Γ k (recall that q h is continuous). We have thus reached a contradiction with the assumptions ∥D(v h )∥ 2 0,Π k +|q h | 2 1,Π Γ k = 1 and α = 1. This proves that there exists α < 1 such that ∥D(v h )∥ 2 0,Π Γ k + (1 -α)h 2 |q h | 2 1,Π Γ k ≤ α∥D(v h )∥ 2 0,Π k + β h 2 ∥ -∆v h + ∇q h ∥ 2 0,Π Γ k + ∥ div v h ∥ 2 0,Π Γ k + h ∥[∂ n v h ]∥ 2 F k + h 3 ∂ 2 n v h 2 F k
on all the patches Π k and for all v h , q h . Summing this over all Π k gives [START_REF] Massing | A stabilized Nitsche fictitious domain method for the Stokes problem[END_REF].

Lemma 3. For any v ∈ H 1 (Ω) d vanishing on Γ w |v| 1,Ω h ≤ C∥D(v)∥ 0,Ω h . (22)
Proof. Since v = 0 on Γ ω , we have the following Korn inequality

|v| 1,Ω ≤ C∥D(v)∥ 0,Ω (23)
with a constant C > 0 depending only on the shape of Ω, cf. [START_REF] Ciarlet | Mathematical Elasticity: Volume I: three-dimensional elasticity[END_REF]. This implies

|v| 1,Ω i h ≤ C∥D(v)∥ 0,Ω h , (24)
where Ω i h denotes the mesh cells inside Ω. Now, for any pair of mesh cells T, T ′ sharing a facet E, we can prove

|v| 1,T ≤ C∥D(v)∥ 0,T + C|v| 1,T ′ (25)
with a constant C independent of h. Indeed, combining the Korn inequalities (23) and the trace theorem on the reference element leads to |v| 1,T ≤ C∥D(v)∥ 0,T + C|v| 1/2,E . Employing again the trace inequality |v| 1/2,E ≤ C|v| 1,T ′ (see [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]Lemma 7.5.26 ]) leads to (25).

Let Ω i,1 h be Ω i h plus the cells which are not in Ω i h but have a neighbor in Ω i h . For any such cell T , we take T ′ as its neighbor in Ω i h , apply the estimate above and sum which gives

|v| 1,Ω i,1 h \Ω i h ≤ C∥D(v)∥ 0,Ω i,1 h \Ω i h + C|v| 1,Ω i h , hence, using (24), |v| 1,Ω i,1 h ≤ C∥D(v)∥ 0,Ω h .
Let Ω i,2 h be Ω i,1 h plus the cells which are not in Ω i,1 h but have a neighbor in Ω i,1 h . We have similar to above

|v| 1,Ω i,2 h ≤ C∥D(v)∥ 0,Ω h
and so on. After a finite number of steps, say k, we arrive at Ω i,k h = Ω h . And

|v| 1,Ω h = |v| 1,Ω i,k h ≤ C∥D(v)∥ 0,Ω h . Lemma 4. For any s h ∈ V h and any V h ∈ R d , ω h ∈ R d ′ , ∥ϕ h s h ∥ 0,Ω Γ h +   F ∈F Γ h h∥ϕ h s h ∥ 2 0,F   1/2 + √ h∥ϕ h s h ∥ 0,G h ≤ Ch∥D(ϕ h s h + χ h (V h + ω h × r))∥ 0,Ω h and |V h | + |ω h | ⩽ C∥ϕ h s h + χ h (V h + ω h × r)∥ 1,Ω h . (26)
Proof. Take any

s h ∈ V h , V h ∈ R d , ω h ∈ R d ′ and denote v h = ϕ h s h + χ h (V h + ω h × r). By [12, Lemma 3.4] ∥ϕ h s h ∥ 0,Ω Γ h ⩽ Ch|ϕ h s h | 1,Ω Γ h ⩽ Ch(|v h | 1,Ω Γ h + |χ h (V h + ω h × r)| 1,Ω Γ h ) ⩽ Ch(|v h | 1,Ω Γ h + ∥χ h ∥ 1,Ω Γ h (|V h | + |ω h |)). (27) By equivalence of norms |V h | + |ω h | ⩽ C∥V h + ω h × r∥ 0,Γ . (28)
Denote by B Γ h the band between Γ and Γ h . Applying the divergence theorem to the vector field |V h + ω h × r| 2 ∇ϕ and noting that the normal on Γ (resp. Γ h ) is given by ± ∇ϕ |∇ϕ| (resp.

± ∇ϕ h |∇ϕ h | ) gives Γ |V h + ω h × r| 2 |∇ϕ| ⩽ Γ h |V h + ω h × r| 2 |∇ϕ • ∇ϕ h | |∇ϕ h | + B Γ h div |V h + ω h × r| 2 ∇ϕ
We now note that |∇ϕ| (resp. |∇ϕ h | and χ h ) are both positive and bounded away from 0 on Γ (resp. on Γ h ) uniformly in h for h small enough, and the measure of B Γ h is of order h k+1 . The inequality above implies thus

∥V h + ω h × r∥ 2 0,Γ ≤ C(∥χ h (V h + ω h × r)∥ 2 0,Γ h + h k+1 (|V h | + |ω h |) 2
). Combining this with (28) gives, for h small enough

|V h | + |ω h | ⩽ C∥χ h (V h + ω h × r)∥ 0,Γ h = C∥v h ∥ 0,Γ h

and, by the trace inequality, |V

h | + |ω h | ⩽ C∥v h ∥ 1,Ω h , i.e. (

Interpolation by finite elements multiplied with the level set

We recall first a Hardy-type inequality, cf. [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF]. Lemma 5. For any integer s ∈ [0, k] and any u ∈ H s+1 (Ω h ) vanishing on Γ, it holds u ϕ s,Ω h ≤ C∥u∥ s+1,Ω h with C > 0 depending only on the constants in Assumption 1 and on s.

This allows us to prove the following bound for interpolation by the products of finite elements with ϕ h . Lemma 6. Let t be an integer

1 ⩽ t ⩽ k + 1. For any v ∈ H t (Ω h ) d ∩ H 1 0 (Ω) there exists w h ∈ V h s.t. ∥v -ϕ h w h ∥ s,Ω h ⩽ Ch t-s ∥v∥ t,Ω h , s = 0, 1 (29)
with C > 0 depending only on t, s, the constants in Assumptions 1-3, and the mesh regularity.

Proof. Let v ∈ H t (Ω h ) d , v = 0 on Γ, and set w = v/ϕ. Thanks to Lemma 5, w ∈ H t-1 (Ω h ) d and ∥w∥ t-1,Ω h ≤ C∥v∥ t,Ω h . Consider w h = I c h w
, where I c h is a Scott-Zhang interpolation operator. For any T ∈ T h , let ω T denote the patch of mesh cells adjacent to T (not necessarily all the adjacent cells) regrouping the cells affected by the construction of I c h w on T , so that I c h w on T depends on w only through its restriction to ω T . The Scott-Zhang interpolation operator can be constructed so that ω T ⊂ T Γ,ext h for all T ∈ T Γ,ext h , and

ω T ⊂ T h \ T Γ,ext h for all T ∈ T h \ T Γ,ext h .
In what follows, we assume that the operator I c h enjoys this property together with the usual interpolation error estimates see for instance [START_REF] Brenner | The mathematical theory of finite element methods[END_REF].

Our first goal is to prove (29) for s = 0. Taking any T ∈ T h . Recall that ϕ is supposed to be of class (at least) C t so that ∥ϕ -

ϕ h ∥ ∞,T ⩽ Ch t . Hence, ∥v -ϕ h w h ∥ 0,T = ∥ϕw -ϕ h w h ∥ 0,T (30)
⩽ ∥ϕ∥ ∞,T ∥w -w h ∥ 0,T + ∥ϕ -ϕ h ∥ ∞,T ∥w h ∥ 0,T ⩽ ∥ϕ∥ ∞,T ∥w -w h ∥ 0,T + Ch t ∥w∥ 0,ω T .
To continue this proof, we distinguish two cases: the cells T ∈ T Γ,ext h close to Γ h and the remaining cells, which are at the distance of at least order h from Γ h .

(i) Consider T ∈ T Γ,ext h . We have ∥ϕ∥ ∞,T ⩽ Ch on these cells since they are at the distance ∼ h from Γ. Noting that ∥w-w h ∥ 0,T ⩽ Ch t-1 |w| t-1,ω T by the usual interpolation estimate, we derive from (30)

∥v -ϕ h w h ∥ 0,T ⩽ Ch t (|w| t-1,ω T + ∥w∥ 0,ω T ) . ( 31
) (ii) Now consider T ∈ T h \ T Γ,ext h
. We note that ϕ does not vanish on ω T for such T (recall that ω T ⊂ T h \ T Γ,ext h ), so that w ∈ H t (ω T ) and, by (30) and the usual approximation estimates, ∥v -

ϕ h w h ∥ 0,T ⩽ Ch t (∥ϕ∥ ∞,T |w| t,ω T + ∥w∥ 0,ω T ) . (32)
In order to bound |w| t,ω T here, we recall the Leibniz rule valid for any multi-index α ∈ N d

∂ α v = ∂ α (ϕw) = β ∈ N d β ⩽ α C β α (∂ β ϕ) (∂ α-β w)
with binomial coefficients C β α depending only on the multi-indices α and β (this formula can be easily proven by induction on the length of α) β ⩽ α means ∀i = 1, . . . , d, β i ⩽ α i . If α ̸ = 0, this can be rewritten, by separating the term with β = 0 (note that C 0 α = 1) and dividing by ϕ, as

∂ α w = 1 ϕ ∂ α v - β ∈ N d β ⩽ α, β ̸ = 0 C β α ∂ β ϕ ϕ ∂ α-β w (33)
Applying (33) to w on ω T gives

|w| t,ω T ⩽ 1 min ω T |ϕ| (|v| t,ω T + C∥w∥ t-1,ω T ).
Hence, by (32),

∥v -ϕ h w h ∥ 0,T ⩽ Ch t ∥ϕ∥ ∞,T min ω T |ϕ| (|v| t,ω T + ∥w∥ t-1,ω T ) + ∥w∥ 0,ω T . ( 34
)
Recall that min ω T |ϕ| ⩾ mh by Assumption 3. The distance between any point on T and any point on ω T is at most 2h so that

∥ϕ∥ ∞,T min ω T |ϕ| = 1 + max T |ϕ| -min ω T |ϕ| min ω T |ϕ| ⩽ 1 + 2M h mh ⩽ 1 + 2M m
with M denoting an upper bound on |∇ϕ|. Substituting this into (34) gives

∥v -ϕ h w h ∥ 0,T ⩽ Ch t (|v| t,ω T + ∥w∥ t-1,ω T ). (35)
Summing (31) over all the cells T ∈ T Γ,ext h and (35) over all the remaining cells of mesh T h gives

∥v -ϕ h w h ∥ 0,Ω h ⩽ Ch t (|v| t,Ω h + ∥w∥ t-1,Ω h ) .
This yields (29) with s = 0 thanks to the estimate ∥w∥ t-1,Ω h ⩽ ∥v∥ t,Ω h given by Lemma 5.

Let us now prove (29

) for s = 1. Introduce v h = I c h v so that |v -v h | 1,Ω h ⩽ Ch t-1 |v| t,Ω h and ∥v -v h ∥ 0,Ω h ⩽ Ch t |v| t,Ω h .
Then, combining the already proven estimate (29) for s = 0, the inverse inequality, and the interpolation estimates above, we obtain

|v -ϕ h w h | 1,Ω h ⩽ |v h -ϕ h w h | 1,Ω h + |v -v h | 1,Ω h ⩽ C h ∥v h -ϕ h w h ∥ 0,Ω h + |v -v h | 1,Ω h ⩽ C h ∥v -ϕ h w h ∥ 0,Ω h + C h ∥v -v h ∥ 0,Ω h + |v -v h | 1,Ω h ⩽ Ch t-1 ∥v∥ t,Ω h .

An adaptation of the Taylor-Hood inf-sup stability to ϕ-FEM

In this section, we prove some velocity-pressure inf-sup conditions that will be used to establish the generalized coercivity (inf-sup) for the full bilinear form in the next section. The proofs are inspired by [START_REF] Guzmán | Inf-sup stability of geometrically unfitted Stokes finite elements[END_REF] and start from an auxiliary inf-sup condition of Lemma 7 with respect to an h-dependent norm for the pressure. The final result in Lemma 8 is weaker than the usual inf-sup since it does not provide the control of the L 2 norm of the pressure over the whole domain Ω h . It will be however sufficient for our purposes since the encumbering term -Ch

2 |p h | 2 1,Ω Γ h
will be controlled by the stabilization present in the scheme, cf. Lemma 2 for this matter.

Lemma 7. There exists an

h-independent constant C > 0 such that ∀p h ∈ M h , ∃w h ∈ V h satisfying h 2 |p h | 2 1,Ω h -Ch 2 |p h | 2 1,Ω Γ h ⩽ Ω h ∇p h • (ϕ h w h ) and |ϕ h w h | 1,Ω h ⩽ Ch|p h | 1,Ω h . ( 36
)
Proof. Let us introduce the submesh

T i h = T h \ T Γ h so that ϕ h ⩽ 0 on T ∈ T i h .
Denote by E i h the set of the edges of the mesh T i h including those shared with T Γ h , but excluding those lying on Γ w . For any edge E ∈ E i h , let t E be the unit tangent vector to E (any of two, but fixed), x E be the midpoint of E, ω(E) be the set of the mesh cells sharing E, and ψ E be the piecewise quadratic function such that ψ E (x E ) = 1 and ψ E vanishes at all the other edge midpoints and at all the nodes of T h . Moreover, define for all

E ∈ E i h , φE =    -h, if E belongs to a cell from T Γ,ext h , 1 |E| E ϕ h , otherwise.
Take any p h ∈ M h and set w h ∈ V h as

w h = E∈E i h h 2 φE ψ E (t E • ∇p h )t E . (37)
We have indeed w h ∈ V h , since the pressure tangential derivative t E •∇p h is a continuous piecewise polynomial of degree ⩽ k -2 on ω(E) and ψ E is a continuous piecewise polynomial of degree 2, vanishing outside ω(E). Note also that w h = 0 on Γ w since E i h does not contain the edges lying on Γ w .

Take any E ∈ E i h and any cell T ∈ ω(E) ∩ T i h . We shall see that

T ϕ h φE ψ E |t E • ∇p h | 2 ≥ c T |t E • ∇p h | 2 . ( 38
)
Here and elsewhere, the constants c > 0 depend only on the polynomial degree k, the shape regularity, and the parameters of Assumption 3. To prove (38), we set φh = ϕ h φE and note that φh ⩾ 0 on T since ϕ h ⩽ 0 on T ∈ T i h . To derive further properties of φh from Assumption 3, we consider 3 following cases with respect to the placement of T and E in the mesh (we recall that T ∈ T i h in any case and E is an edge belonging to T ). . Then φh = ϕ h (x E ) for some point x E ∈ E and we have for any x ∈ T φE

ϕ h (x) = ϕ h (x E ) ϕ h (x) = 1 + ϕ h (x E ) -ϕ h (x) ϕ h (x) = 1 + ∇ϕ h (c) • (x E -x) ϕ h (x) ≤ 1 + M h mh = m + M m
where c ∈ T is a point on the segment connecting x E with x. This implies

φh ≥ m m + M on T . ( 40
) Case 3 T ̸ ∈ T Γ,ext
h , E is shared with a cell from T Γ,ext h . Then φE = -h and 1 |E| E ϕ h ≤ -c 1 h with some c 1 > 0 depending only on the constants in Assumption 3 (since the distance between E and Γ h , where ϕ h vanishes, is of order h and |∇ϕ h | is bounded away from 0 on T Γ,ext h ). Combining this with the arguments of the previous case, we arrive at

φh ≥ c 1 m m + M on T . (41)
Moreover, in all of the 3 cases above,

| φh | ⩽ C on T . ( 42
)
with some C > 0 depending only on the constants in Assumption 3. In case 1, this follows from the bound |∇ϕ h | ≤ M and the fact that the maximal distance between T and Γ h is of order h.

In cases 2 and 3, (42) can be proven in the same way as ( 40) and (41). Inequality (38) can be now proven setting

c = min T, φh ,q h T φh ψ E q 2 h T q 2 h ( 43
)
where the minimum is taken over all the simplexes T permitted by the mesh regularity, all the polynomials φh := ϕ h φE of degree k satisfying (42) and either of (39)-( 40)-( 41), and all the polynomials q h := t E • ∇p h ̸ = 0 of degree k -2. By homogeneity and rescaling, one can safely assume that h = 1 and ∥q h ∥ 0,T = 1. The sets of possible T , φh , q h are bounded and closed, so that the minimum in (43) is indeed attained and c > 0. Indeed, either of (39)-( 40)-(41) excludes the possibility of φh vanishing everywhere on T . This concludes the proof of (38).

Thanks to (38), we have, setting Ω i h = Ω h \ Ω Γ h and denoting by E(T ) the set of edges of a cell T excluding the edges on Γ w ,

Ω i h ϕ h w h • ∇p h = E∈E i h h 2 ω(E)∩Ω i h ϕ h φE ψ E |t E • ∇p h | 2 ≥ c T ∈T i h E∈E(T ) h 2 T |t E • ∇p h | 2 .
Taking into account Assumption 5, we have by scaling and the equivalence of norms on all

T ∈ T i h E∈E(T ) T |t E • ∇p h | 2 ≥ c|p h | 2 1,T .
Hence,

Ω i h ϕ h w h • ∇p h ⩾ ch 2 |p h | 2 1,Ω i h . (44)
Now, on a mesh cell T ∈ T Γ h having common edges with cells of T i h , definition (37) clearly gives

T ϕ h w h • ∇p h = E∈E(T )∩E i h h T ϕ h ψ E |t E • ∇p h | 2 ⩽ Ch 2 |p h | 2 1,T
since ϕ h is of order h on such cells. Combining this with (44) gives

ch 2 |p h | 2 1,Ω h ⩽ Ω h ϕ h w h • ∇p h - Ω Γ h ϕ h w h • ∇p h + ch 2 |p h | 2 1,Ω Γ h ⩽ Ω h ϕ h w h • ∇p h + Ch 2 |p h | 2 1,Ω Γ h . ( 45
)
We also have

|ϕ h w h | 2 1,Ω h = Ω h |(∇ϕ h )w h + ϕ h (∇w h )| 2 ≤ C T ∈T h E∈E(T )∩E i h h 4 T |∇ϕ h | 2 φ2 E |∇p h | 2 + ϕ 2 h φ2 E |∇ψ E | 2 |∇p h | 2 + ϕ 2 h φ2 E |∇ 2 p h | 2 ≤ C T ∈T h h 2 T |∇p h | 2 = Ch 2 |p h | 2 1,Ω h . ( 46
)
We have used here the finite element inverse estimates on ϕ h , ψ E , p h , and the uniform upper bound (42

) on | φh | = |ϕ h | | φE | .
Redefining w h as w h /c with the constant c from (45) shows (36) as a combination of (45) and (46).

Lemma 8. There exists an

h-independent constant C > 0 such that ∀p h ∈ M h ∃s h ∈ V h ∥p h ∥ 2 0,Ω h -Ch 2 |p h | 2 1,Ω Γ h ≤ Ω h ∇p h • (ϕ h s h ) and |ϕ h s h | 1,Ω h ⩽ C∥p h ∥ 0,Ω h . (47)
Proof. Taking p h ∈ M h . By continuous velocity-pressure inf-sup (recall that Ω p h = 0), there exists v ∈ H 1 0 (Ω) s.t. div v = -p h on Ω, and ∥v∥ 1,Ω ⩽ C∥p h ∥ 0,Ω .

Let ṽ ∈ H 1 (Ω h ) be the extension of v by 0 outside Ω. Lemma 6 with t = 1 implies ∃v h ∈ V h s.t.

∥ṽ -ϕ h v h ∥ 0,Ω h ⩽ Ch∥ṽ∥ 1,Ω h ⩽ Ch∥p h ∥ 0,Ω h and |ϕ h v h | 1,Ω h ⩽ C∥ṽ∥ 1,Ω h ⩽ C∥p h ∥ 0,Ω h .
Thus,

∥p h ∥ 2 0,Ω = - Ω p h div v = Ω ∇p h • ṽ = Ω h ∇p h • (ϕ h v h ) + Ω h ∇p h • (ṽ -ϕ h v h ) ⩽ Ω h ∇p h • (ϕ h v h ) + Ch|p h | 1,Ω h ∥p h ∥ 0,Ω h . ( 48
)
We have

∥p h ∥ 2 0,Ω h ⩽ C 1 (∥p h ∥ 2 0,Ω + h 2 |p h | 2 1,Ω Γ h ).
This can be proven by an argument similar to that in Lemma 2: one can consider the maximum of

∥p h ∥ 2 0,Π k ∥p h ∥ 2 0,T k +h 2 |p h | 2 1,Π Γ k
taken over all the admissible patches Π k , as in Assumption 4, and piecewise polynomials p h , observe that this maximum is attained and takes the value C 1 > 0, and sum up over all the patches covering Ω Γ h . We can thus pass from the norm on Ω to that on Ω h in (48):

∥p h ∥ 2 0,Ω h ⩽ C 1 Ω h ∇p h • (ϕ h v h ) + Ch|p h | 1,Ω h ∥p h ∥ 0,Ω h + C 1 h 2 |p h | 2 1,Ω Γ h .
Hence, by Young inequality,

1 2 ∥p h ∥ 2 0,Ω h ⩽ C 1 Ω h ∇p h • (ϕ h v h ) + C 2 h 2 |p h | 2 1,Ω h + C 1 h 2 |p h | 2 1,Ω Γ h ⩽ C 1 Ω h ∇p h • (ϕ h v h ) + C 2 Ω h ∇p h • (ϕ h w h ) + Ch 2 |p h | 2 1,Ω Γ h
with w h given by Lemma 7. Thus,

1 2 ∥p h ∥ 2 0,Ω h -Ch 2 |p h | 2 1,Ω Γ h ⩽ Ω h ∇p h • ϕ h (C 1 v h + C 2 w h ) and |ϕ h (C 1 v h + C 2 w h )| 1,Ω h ⩽ C∥p h ∥ 0,Ω h since |ϕ h w h | 1,Ω h ⩽ Ch|p h | 1,Ω h ⩽ C∥p h ∥ 0,Ω h by inverse finite element estimates. Setting s h = 2(C 1 v h + C 2 w h ) proves (47).

3.4

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

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

V rbm h = {χ h (V h + ω h × r) + ϕ h s h with s h ∈ V h , V h ∈ R d , ω h ∈ R d ′ }. (49)
In the scheme [START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF], we shall now combine the test functions s h , V h , ω h into v h ∈ V rbm h as in the definition above. Similarly, we shall combine the trial functions w [START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF] can be then rewritten in the compact form:

h , U h , ψ h into u h ∈ V rbm h setting u h = χ h (U h + ψ h × r) + ϕ h w h ). Scheme
find u h ∈ V rbm h and p h ∈ M h such that c h (u h , p h ; v h , q h ) = L h (v h , q h ), ∀v h ∈ V rbm h , q h ∈ M h , (50)
where the bilinear form c h is given by

c h (u h , p h ; v h , q h ) = 2 Ω h D(u h ) : D(v h ) - G h (2D(u h ) -p h I)n • ϕ h s h - Ω h q h div u h - Ω h p h div v h + σh 2 T ∈T Γ h T (-∆u h + ∇p h ) • (-∆v h -∇q h ) + σ T ∈T Γ h T (div u h )(div v h ) + σ u h E∈F Γ h E [∂ n u h ] • [∂ n v h ] + σ u h 3 E∈F Γ h E ∂ 2 n u h • ∂ 2 n v h
and the linear form L h is given by

L h (v h , q h ) = Ω h ρ f g • ϕ h s h + O ρ f g • χ h (V h + ω h × r) + 1 - ρ f ρ s mg • V h + σh 2 T ∈T Γ h T ρ f g • (-∆v h -∇q h ). ( 51
)
In both expressions above, s h , V h , and ω h are related to v h ∈ V rbm h as in (49).

Lemma 9. Introduce the norm on

V rbm h × M h v h , q h h :=   |v h | 2 1,Ω h + ∥q h ∥ 2 0,Ω h + h 2 T ∈T Γ h ∥ -∆v h + ∇q h ∥ 2 0,T + J u (v h , v h )   1/2
.

The following inf-sup condition holds provided σ and σ u are sufficiently large:

∀(u h , p h ) ∈ V rbm h × M h ∃(v h , q h ) ∈ V rbm h × M h such that c h (u h , p h ; v h , q h ) v h , q h h ⩾ θ u h , p h h ( 52
)
with a constant θ > 0 depending only on the mesh regularity.

Proof. Let us take p h ∈ M h and u

h = χ h (U h + ψ h × r) + ϕ h w h ∈ V rbm h with w h ∈ V h , U h ∈ R d and ψ h ∈ R d ′ .
Step 1: controlling the velocity. By choosing (u h , p h ) as the trial function and (u h , -p h ) as the test function in the bilinear form c h , we obtained:

c h (u h , p h ; u h , -p h ) = 2 Ω h |D(u h )| 2 - G h (2D(u h ) -p h I)n • ϕ h w h + σ u h E∈F Γ h E |[∂ n u h ]| 2 + σ u h 3 E∈F Γ h E ∂ 2 n u h 2 + σh 2 T ∈T Γ h T | -∆u h + ∇p h | 2 + σ T ∈T Γ h T | div u h | 2 . ( 53
)
Let Bh be the strip between Γ h = {ϕ h = 0} and G h , i.e.

Bh = {ϕ h > 0} ∩ Ω h . Since ϕ h w h = 0 on Γ h , G h D(u h )n • ϕ h w h = ∂ Bh D(u h )n • ϕ h w h = T ∈T Γ h ∂( Bh ∩T ) D(u h )n T • ϕ h w h - T ∈T Γ h E∈F cut h (T ) Bh ∩E D(u h )n T • ϕ h w h ,
where T Γ h is defined in [START_REF] Cotin | ϕ-fem: an efficient simulation tool using simple meshes for problems in structure mechanics and heat transfer[END_REF], F cut h (T ) regroups the facets of a mesh element T cut by Γ h , and n T is the unit normal pointing outside of T on the boundary of a mesh cell T . Applying the divergence theorem to the integrals on ∂( Bh ∩ T ) and regrouping the integrals on the facets gives

G h D(u h )n • ϕ h w h = Bh D(u h ) : D(ϕ h w h ) + T ∈T Γ h Bh ∩T div D(u h ) • ϕ h w h - E∈F Γ h E∩ Bh [D(u h )n] • ϕ h w h = Bh |D(u h )| 2 + T ∈T Γ h Bh ∩T div D(u h ) • ϕ h w h - E∈F Γ h E∩ Bh [D(u h )n] • ϕ h w h , since u h -ϕ h w h is
the velocity of a rigid motion on Bh .

Similarly (and simpler)

G h p h n • ϕ h w h = ∂ Bh p h n • ϕ h w h = Bh p h div u h + Bh ∇p h • ϕ h w h .
Substituting this into (53) and rewriting 2 div

D(u h ) -∇p h = ∆u h -∇p h + ∇ div u h on the cells T ∈ T Γ h yields c h (u h , p h ; u h , -p h ) = 2 Ω h |D(u h )| 2 -2 Bh |D(u h )| 2 - T ∈T Γ h Bh ∩T (∆u h -∇p h ) • ϕ h w h Young with ε1 + 2 F ∈F Γ h F ∩ Bh [D(u h )n] • ϕ h w h Young with ε2 + Bh p h div u h Young with ε3 - T ∈T Γ h Bh ∩T (∇ div u h ) • ϕ h w h Young with ε4 + σ u h E∈F Γ h E | [∂ n u h ] | 2 + σ u h 3 E∈F Γ h E ∂ 2 n u h 2 + σh 2 T ∈T Γ h T |∆u h -∇p h | 2 + σ T ∈T Γ h T | div u h | 2 .
Several terms above are marked with "Young with ε i " meaning that we are going to apply the Young inequality with some weights ε 1 , . . . , ε 4 > 0 (multiplied by the appropriate powers of h) to these terms. We recall that Lemma 4 implies

T ∈T Γ h 1 h 2 ∥ϕ h w h ∥ 2 Bh ∩T ≤ C∥D(u h )∥ 2 0,Ω h , F ∈F Γ h 1 h ∥ϕ h w h ∥ 2 F ∩ Bh ≤ C∥D(u h )∥ 2 0,Ω h ,
which allows us to absorb the norms of ϕ h w h into the first term with ∥D(u h )∥ 0,Ω h . We also use the inverse inequality h∥∇ div

u h ∥ 0,Ω Γ h ≤ C∥ div u h ∥ 0,Ω Γ h . This yields c h (u h , p h ; u h , -p h ) ⩾ 2 -C ε 1 + ε 2 + ε 4 2 ∥D(u h )∥ 2 0,Ω h -2∥D(u h )∥ 2 0,Ω Γ h - ε 3 2 ∥p h ∥ 0,Ω Γ h + h 2 σ - 1 2ε 1 T ∈T Γ h ∥∆u h -∇p h ∥ 2 0,T + σ - 1 2ε 3 - C ε 4 ∥ div u h ∥ 2 0,Ω Γ h + h σ u - 1 2ε 2 E∈F Γ h ∥ [∂ n u h ] ∥ 2 0,E + σ u h 3 E∈F Γ h ∂ 2 n u h 2 0,E
Thanks to Lemma 2, this can be further bounded as

c h (u h , p h ; u h , -p h ) ⩾ 2(1 -α) -C ε 1 + ε 2 + ε 4 2 ∥D(u h )∥ 2 0,Ω h + 2(1 -α)h 2 |p h | 1,Ω Γ h - ε 3 2 ∥p h ∥ 0,Ω Γ h + h 2 σ - 1 2ε 1 -2β T ∈T Γ h ∥∆u h -∇p h ∥ 2 0,T + σ - 1 2ε 3 - C ε 4 -2β ∥ div u h ∥ 2 0,Ω Γ h + h σ u - 1 2ε 2 -2β E∈F Γ h ∥ [∂ n u h ] ∥ 2 0,E + h 3 (σ u -2β) E∈F Γ h ∂ 2 n u h 2 0,E (54)
with some β > 0 and α ∈ (0, 1).

Step 2: controlling the pressure. Let now s h ∈ V h be the function given by Lemma 8 and set v p h = ϕ h s h . Noting that

- Ω h p h div v p h + G h p h n • v p h = Ω h ∇p h • v p h ⩾ ∥p h ∥ 2 0,Ω h -Ch 2 |p h | 2 1,Ω Γ h we get c h (u h , p h ; v p h , 0) ⩾ ∥p h ∥ 2 0,Ω h -Ch 2 |p h | 2 1,Ω Γ h + 2 Ω h D(u h ) : D(v p h ) -2 G h D(u h )n • v p h + σh 2 T ∈T Γ T (∆u h -∇p h ) • ∆v p h + σ T ∈T Γ h T (div u h )(div v p h ) + σ u h E∈F Γ h E [∂ n u h ] • [∂ n v p h ] + σ u h 3 E∈F Γ h E ∂ 2 n u h • ∂ 2 n v p h .

Recalling that |v

p h | 1,Ω h ⩽ C∥p h ∥ 0,Ω h , remarking that ∥v p h ∥ 0,G h ≤ C √ h |v p h | 1,Ω Γ h ,

and applying Young and inverse inequalities allows us to conclude

c h (u h , p h ; v p h , 0) ⩾ 1 2 ∥p h ∥ 2 0,Ω h -Ch 2 |p h | 2 1,Ω Γ h -C∥D(u h )∥ 2 0,Ω h -C   σ 2 h 2 T ∈T Γ h ∥ -∆u h + ∇p h ∥ 2 0,T + σ 2 T ∈T Γ h ∥ div u h ∥ 2 0,T σ 2 u h E∈F Γ h ∥[∂ n u h ]∥ 2 0,E + σ 2 u h 3 E∈F Γ h ∂ 2 n u h 2 0,E   . ( 55
)
Step 3: combining the estimates. Multiply (55) by λ > 0 and add it to (54). This gives

c h (u h , p h ; u h + λv p h , -p h ) ⩾ 2(1 -α) -C ε 1 + ε 2 + ε 4 2 -Cλ ∥D(u h )∥ 2 0,Ω h + λ -ε 3 2 ∥p h ∥ 2 0,Ω h + (2(1 -α) -Cλ)h 2 |p h | 2 1,Ω Γ h + h 2 σ - 1 2ε 1 -2β -Cλσ 2 T ∈T Γ h ∥∆u h -∇p h ∥ 2 0,T + σ - 1 2ε 3 - C ε 4 -2β -Cλσ 2 ∥ div u h ∥ 2 0,Ω Γ h + h σ u - 1 2ε 2 -2β -Cλσ 2 u E∈F Γ h ∥ [∂ n u h ] ∥ 2 0,E + (σ u -2β -Cλσ 2 u )h 3 E∈F Γ h ∂ 2 n u h 2 0,E .
Taking λ, ε 1 , ε 2 , ε 4 small enough, ε 3 = λ/2 and σ, σ u big enough, and recalling Korn inequality [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF], this amounts to

c h (u h , p h ; w h + λv p h , -p h ) ⩾ c u h , p h 2 h
with some c > 0. We also have easily

w h + λv p h , -p h h ⩽ C w h , p h h ,
hence the inf-sup estimate (52) with v h = u h + λv p h , q h = -p h .

A priori error estimates.

In this section, we will prove Theorem 1 following the argumentation of [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF], which is ameliorated since we require only the optimal regularity H k+1 (Ω) d × H k (Ω) for the velocity-pressure pair (u, p) given by (1a)-(1g).

Proof of the H 1 a priori error estimate (12): Let (u, p) ∈ H k+1 (Ω) d × H k (Ω) with u = U +ψ ×r on Γ be the solution to the continuous problem (1a)-(1g) and (u h , p h ) ∈ V rbm h ×M h with u h = ϕ h w h +χ h (U h +ψ h ×r) be the solution to the discrete problem [START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF]. Choose sufficiently smooth extension ũ and p of u and p on Ω h such that ũ = u, p = p on Ω, and

∥ũ∥ k+1,Ω h ⩽ C∥u∥ k+1,Ω , ∥p∥ k,Ω h ⩽ C∥p∥ k,Ω .
Applying Lemma 6 to ũ -χ(U + ψ × r), which vanishes on Γ and on Γ w , we see that there exists wh ∈ V h such that

∥ũ -χ(U + ψ × r) -ϕ h wh ∥ 1,Ω h ⩽ Ch k ∥ũ -χ(U + ψ × r)∥ k+1,Ω h .

This allows us to introduce ũh

= ϕ h wh + χ h (U + ψ × r) ∈ V rbm h satisfying ∥ũ -ũh ∥ 1,Ω h ⩽ ∥ũ -χ(U + ψ × r) -ϕ h wh ∥ 1,Ω h + ∥(χ -χ h )(U + ψ × r)∥ 1,Ω h ⩽ Ch k (∥ũ∥ k+1,Ω h + ∥χ∥ k+1,Ω h (|U | + |ψ|)) ⩽ Ch k ∥ũ∥ k+1,Ω h ⩽ Ch k ∥u∥ k+1,Ω , (56)
thanks to the standard interpolation of χ ∈ H k+1 (Ω h ) and to the bounds |U |, |ψ| ≤ C∥u∥ 1,Ω valid by the trace inequality (recall that u = U + ψ × r on Γ).

Similarly, ∥ũ -ũh ∥ 0,Ω h ⩽ Ch k+1 ∥u∥ k+1,Ω . We define moreover ph ∈ M h by the standard FE nodal interpolation ph = I h p such that

∥p -ph ∥ 0,Ω h ⩽ Ch k ∥p∥ k,Ω . (57)
Thanks to Lemma 9,

∃(v h , q h ) ∈ V rbm h × M h such that ũh -u h , ph -p h h ≤ C c h (ũ h -u h , ph -p h ; v h , q h ) v h , q h h . (58)
We should now substitute (ũ, p) into the form c h. To this end, we introduce the fictitious righthand sides F and Q on Ω h so that

-2 div D(ũ) + ∇p = F and div ũ = Q in Ω h .
We observe then, taking any

v h = ϕ h s h + χ h (V h + ω h × r) ∈ V rbm h , q h ∈ M h , 2 Ω h D(ũ) : D(v h ) - Ω h p div v h - G h (2D(ũ) -pI)n • v h = Ω h F • v h and, recalling B h = Ω h \ Ω, G h (2D(ũ) -pI)n • (V h + ω h × r) = Γ (2D(ũ) -pI)n • (V h + ω h × r) + B h div(2D(ũ) -pI)n • (V h + ω h × r) = mg • V h - B h F • (V h + ω h × r).
Hence,

c h (ũ, p; v h , q h ) = Ω h F • v h - B h F • (V h + ω h × r) + mg • V h - Ω h q h div ũ + σh 2 T ∈T Γ h T F • (-∆v h -∇q h ) + σ Ω Γ h (div ũ)(div v h ).
Also note that the RHS (51) of the scheme (50) can be rewritten as

L h (v h , q h ) = Ω h ρ f g • v h - B h ρ f g • (V h + ω h × r) + mg • V h + σh 2 T ∈T Γ h T ρ f g • (-∆v h -∇q h )
This allows us to establish the following Galerkin orthogonality relation, valid for all

v h ∈ V rbm h , q h ∈ M h , c h (ũ -u h , p -p h ; v h , q h ) = R h (v h , q h ), (59)
where

R h (v h , q h ) = B h ( F -ρ f g) • ϕ h s h - B h q h div ũ + σh 2 T ∈T Γ h T ( F -ρ f g) • (-∆v h -∇q h ) + σ Ω Γ h (div ũ)(div v h ).
The integrals of F -ρ f g and div ũ on Ω h have been rewritten as integrals on B h since both F -ρ f g and div ũ vanish on Ω. Combination of (58) and (59) entails

ũh -u h , ph -p h h ≤ C c h (ũ h -ũ, ph -p; v h , q h ) + R h (v h , q h ) v h , q h h .
We can now use interpolation inequalities as in [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF]Section 3.4]. In particular, the term with c h in the nominator of the fraction above is bounded by Ch k (∥u∥ k+1,Ω + ∥p∥ k,Ω ) v h , q h h thanks to (56)-( 57) and to the estimates of Lemma 4. To bound R h (v h , q h ) we recall that F -ρ f g and div ũ vanish on Ω. Thus, thanks to [12, Lemma 3.6]

∥ F -ρ f g∥ 0,Ω Γ h ≤ Ch k-1 ∥ F -ρ f g∥ k-1,Ω Γ h ≤ Ch k-1 (∥u∥ k+1,Ω + ∥p∥ k,Ω ) and ∥ div ũ∥ 0,Ω Γ h ≤ Ch k ∥ div ũ∥ k,Ω Γ h ≤ Ch k ∥u∥ k+1,Ω . (60)
This, combined with the estimates of Lemma 4, in particular

∥ϕ h s h ∥ 0,Ω Γ h ≤ Ch|v h | 1,Ω Γ h , leads to |R h (v h , q h )| ≤ Ch k (∥u∥ k+1,Ω + ∥p∥ k,Ω ) v h , q h h and ũh -u h , ph -p h h ≤ Ch k (∥u∥ k+1,Ω + ∥p∥ k,Ω ) . (61)
Recalling again the interpolation estimates ( 56)-(57), we obtain the error estimates for u and p, announced by [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF].

Proof of the a priori error estimate (13) on the velocity of the solid: We have by the construction of the interpolant ũh = ϕ h wh + χ h (U + ψ × r) and thanks to (26) [START_REF] Ern | Theory and practice of finite elements[END_REF] thanks to (61).

|U -U h | + |ψ -ψ h | ⩽ C∥ϕ h ( wh -w h ) + χ h (U -U h + (ψ -ψ h ) × r)∥ 1,Ω h = C∥ũ h -u h ∥ 1,Ω h which proves
Proof of the L 2 a priori error estimate (14

): Let (v, q, V, ω) ∈ H 2 (Ω) d ×H 1 (Ω)×R d ×R d ′ the solution to                      Let (ṽ, q) ∈ H 2 (Ω h ) d × H 1 (Ω h
) coincide with (v, q) on Ω. They can be constructed by a bounded extension operator in H 2 × H 1 so that

∥ṽ∥ 2,Ω h + ∥q∥ 1,Ω h ⩽ C(∥v∥ 2,Ω + ∥q∥ 1,Ω ) ⩽ C∥u -u h ∥ 0,Ω . (63)
We now further rewrite (62) using Galerkin orthogonality (59) with the test functions

v h = ϕ h s h + χ h (V h + ω h × r) ∈ V rbm h and q h ∈ M h and recalling B h = Ω h \Ω, ∥u -u h ∥ 2 0,Ω = Ω h 2D(ũ -u h ) : D(ṽ -v h ) - Ω h (q -q h ) div(ũ -u h ) - Ω h (p -p h ) div(ṽ -v h ) I - B h 2D(ũ -u h ) : D(ṽ) - B h q div(ũ -u h ) - B h (p -p h ) div ṽ II + G h (2D(ũ -u h ) -(p -p h )I)n • ϕ h s h III (64) -σ u h E∈F Γ h E ∂ ∂n (ũ -u h ) • ∂v h ∂n -σ u h 3 E∈F Γ h E ∂ 2 ∂n 2 -u h ) • ∂ 2 v h ∂n 2 IV -σh 2 T ∈T Γ h T -∆(ũ -u h ) + ∇(p -p h )) • (-∆v h -∇q h ) V -σ T ∈T Γ h T div(ũ -u h ) div v h V I + R h (v h , q h ) V II + Γ ((ϕ h -ϕ)w h + (χ h -χ)(U h + ϕ h × r)) • (2D(v) -qI)n V III
We now take V h = V, ω h = ω and set s h ∈ V h so that ϕ h s h is an optimal interpolant of v -χ h (V + ω × r), as guaranteed by Lemma 6. We also set q h = Ĩh q using an appropriate Clément interpolation Ĩh . We can now estimate all the terms of (64) using the already proven estimate (61) and the interpolation estimates for ṽ -v h and p -p h . This gives

∥u -u h ∥ 2 0,Ω ≤ Ch k+1/2 (∥u∥ k+1,Ω + ∥p∥ k,Ω )(∥ṽ∥ 2,Ω h + ∥q∥ 1,Ω h ). (65)
In particular, term I is completely standard and gives in fact a contribution of the optimal order h k+1 . Rather than go to the details of the tedious calculations leading to the bounds of the remaining terms, we prefer here to refer to the similar arguments used in [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF] to estimate the terms in eq. (3.24). Indeed, the terms II -III in (64) can be treated as the terms II -III in eq. (3.24) of [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF]. Terms IV -V in (64) can be treated as term IV in eq. (3.24) of [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF]. Terms V I in (64) is also similar to the latter (note, in particular, ∥ div ṽ∥ 0,Ω Γ h ≤ Ch∥v∥ 2,Ω similarly to (60)). Finally, term V II in (64) can be treated as term V in eq. (3.24) of [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF]. As in [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF], all these terms result in the sub-optimal estimate of order O(h k+1/2 ). The origin of this sub-optimality lies in the lack of adjoint consistency in formulation (50): the adjoint discrete problem cannot be interpreted as a consistent discretization of a meaningful continuous problem.

The only term in (64), which does not have a direct analogue in [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF], is term V III. To bound it, we apply Cauchy-Schwarz inequality together with the interpolation estimates on ϕ -ϕ h and χ -χ h , recalling the hypotheses ϕ ∈ C k+1 (Ω h ) and χ ∈ H k+1 (Ω h ):

|V III| ⩽ Ch k+1 (∥w h ∥ 0,Γ + |U h | + |ϕ h |)(∥D(v)∥ 0,Γ + ∥q∥ 0,Γ ) .
Then, to bound w h in L 2 (Γ), we start by a trace inverse inequality ∥w h ∥ 0,Γ ⩽ C √ h ∥w h ∥ 0,Ω h and apply Hardy inequality of Lemma 5 to w h = ϕw h ϕ :

∥w h ∥ 0,Ω h ⩽ C∥ϕw h ∥ 1,Ω h ⩽ C(∥(ϕ -ϕ h )w h ∥ 1,Ω h + ∥ϕ h w h ∥ 1,Ω h ) .
Noting that, by interpolation and inverse inequalities,

∥(ϕ -ϕ h )w h ∥ 1,Ω h ⩽ C(∥ϕ -ϕ h ∥ L ∞ (Ω h ) |w h | 1,Ω h + ∥∇ϕ -∇ϕ h ∥ L ∞ (Ω h ) ∥w h ∥ 0,Ω h ) ⩽ Ch k ∥w h ∥ 0,Ω h we conclude (1 -Ch k )∥w h ∥ 0,Ω h ⩽ C∥ϕ h w h ∥ 1,Ω h
Hence, for h small enough,

|V III| ⩽ Ch k+1 1 √ h ∥ϕ h w h ∥ 1,Ω h + |U h | + |ϕ h | (∥D(v)∥ 0,Γ + ∥q∥ 0,Γ ) . Recalling that u h = ϕ h w h + χ h (U h + ϕ h × r), we conclude by Lemma 4 that ∥ϕ h w h ∥ 1,Ω h , |U h |, |ϕ h
| can be all bounded by ∥u h ∥ 1,Ω h . Applying the trace inequalities to v and q, we arrive at

|V III| ⩽ Ch k+1/2 ∥u h ∥ 1,Ω h (∥v∥ 2,Ω + ∥q∥ 1,Ω ) .
Since we know that ∥u h ∥ 1,Ω h is bounded by the norms of u and p thanks to the already proven error estimates for the velocity in H 1 norm, we conclude that term V III contributes to (65) in the same manner as all the other terms. Combining (65) with (63) proves [START_REF] Fournié | Stability and optimal convergence of unfitted extended finite element methods with Lagrange multipliers for the Stokes equations[END_REF].

Numerical tests

In this section, we present numerical results, first in the particular case of a fixed particle, i.e. for the Stokes equations aone in a fixed domain (cf. Appendix A and the ϕ-FEM scheme (72)), and second in the case of the particulate flows (equations (1a)-(1g) and the ϕ-FEM scheme [START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF]). These schemes will be compared with standard FEM on fitted triangular meshes as on Fig. 3 left (we do not introduce higher order approximations of the curvilinear boundary of the domain, as would be the case in the isoparametric FEM for example). In the case of Stokes equations, the error is measured with respect to a manufactured solution, while a reference solution obtained by standard FEM on a fitted fine mesh is used in the case of particulate flows. We have implemented ϕ-FEM in multiphenics [START_REF] Ballarin | multiphenics[END_REF]. The implementation scripts can be consulted on GitHub. 2 The fluid/solid domain in both our test cases is O = (0, 1) 2 ⊂ R 2 and we take the particle S as a disk of radius R = 0.21 centered at a point (0.5, 0.5). Then Ω = [0, 1] 2 \S. The geometry is presented in Fig. 1 (left). In ϕ-FEM, we use the following level-set function, well defined and smooth for all (x, y) ∈ R 2 , ϕ(x, y) = R 2 -(x -0.5) 2 -(y -0.5) 2 .

(66)

We present only the results with the lowest order Taylor-Hood elements, i.e. setting k = 2 and thus using P 2 elements for w h , ϕ h , χ h (the approximation ϕ h for ϕ is exact in this case). The tests with elements of higher order would lead to essentially the same observations. The stabilization parameters are set to σ = σ u = 20 (as in [START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF]). We shall test the ϕ-FEM scheme given by (72) in Appendix A and compare it with a standard Taylor-Hood FEM on a fitted mesh. To this end, we introduce a quasi-uniform triangular mesh T f it h fitted to Ω, in the sense that the boundary nodes of the mesh lie on Γ ∪ Γ w . The domain occupied by this mesh, denoted by Ω f it h , is a polygonal approximation of Ω, since the interface Γ is curvilinear and cannot be represented exactly by the straight edges. We introduce then the FE spaces

V f it,u D h = v h ∈ C( Ωfit h ) d : v h | T ∈ P 2 (T ) d ∀T ∈ T f it h , v h = I h u D on Γ w ∩ Γ f it h (67) M f it h = q h ∈ C( Ωfit h ) : q h | T ∈ P 1 (T ) ∀T ∈ T f it h , Ω f it h q h = 0 , (68)
where Γ f it h is the part of the boundary of Ω f it h approximating Γ. A standard fitted Taylor-Hood FEM can be written as: find

(u h , p h ) ∈ V f it,u D h × M f it h such that Ω f it h 2D(u h ) : D(v h ) - Ω f it h p h div v h - Ω f it h q h div u h = Ω f v h , (69)
for all (v h , q h ) ∈ V f it,0 h × M f it h . Examples of meshes for the standard FEM formulation and the ϕ-FEM scheme are given in Fig. 3 (left) and (right), respectively. In Fig. 4, we report convergence results for the standard fitted Taylor-Hood FEM (69) and for ϕ-FEM (72) in the case of Stokes equations. We recover the theoretical rates of convergence of ϕ-FEM stated in Theorem 2: 2nd order in h for both the H 1 -error in velocity and the L 2 -error in pressure. The L 2 -error in velocity is actually better than theoretically predicted: it is of order 3 instead of theoretically expected 2.5. We observe thus that ϕ-FEM is fully optimal in practice: it demonstrates the same convergence rates in all the aforementioned norms as the standard FEM would demonstrate in the ideal situation of a fitted mesh on a convex polygonal domain. Actually, our setting is not ideal, Ω is neither convex, not polygonal. It is thus not surprising that the standard Taylor-Hood FEM underperforms (we recall that the mesh T f it h is triangular with straight edges and no higher order geometrical approximation is introduced there). Experimentally observed convergence rates for this scheme are ≈ 2 for the L 2 -error in velocity, ≈ 1.5 for the H 1 -error in velocity, and slightly smaller than 2 for the L 2 -error in pressure (the error in pressure is thus the only one for which the optimal convergence order seems to be retained in practice). Moreover, all the errors on all the considered meshes are systematically smaller for ϕ-FEM than for the fitted FEM. Remark 6. As already noted in Remark 2, in ϕ-FEM, it is impossible to impose Ω p h = 0. In our implementation, we rather impose Ω h p h = 0 with the help of a Lagrange multiplier, i.e. we add λ h Ω h p h + µ h Ω h q h (with λ h , µ h ∈ R) to the formulation. To compute the relative L 2error for the pressure, we should compare p h with the exact pressure whose integral vanishes over Ω h (recall that the pressure is physically defined up to an additive constant any way). We thus introduce p = p -c Ω h with c Ω h = 1 |Ω h | Ω h p and compute the errors with respect to p. Similarly, in the case of standard fitted FEM, we impose Ω f it h p h = 0 by a Lagrange multiplier and compute the relative error against p

= p -c Ω f it h with c Ω f it h = 1 |Ω f it h | Ω f it h p. 10 -2 10 -1 10 -7 10 -6 10 -5 10 -4 10 -3 1 2 1 3 h ∥u -uh∥ 0,Ω f it h /∥u∥ 0,Ω f it h st.-FEM ∥u -uh∥0,Ω h /∥u∥0,Ω h ϕ-FEM 10 -2 10 -1 10 -5 10 -4 10 -3 10 -2 1 2 h ∥u -uh∥ 1,Ω f it h /∥u∥ 1,Ω f it h st.-FEM ∥u -uh∥1,Ω h /∥u∥1,Ω h ϕ-FEM 10 -2 10 -1 10 -4 10 -3 10 -2 10 -1 1 2 h ∥p -ph∥ 0,Ω f it h /∥ f it p∥ 0,Ω f it h st.-FEM ∥p -ph∥0,Ω h /∥p∥0,Ω h ϕ-FEM

Particulate flows

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

χ(r) =      1, for r < r 0 1 + f (r0,r1) (r1-r0) 5 , for r 0 < r < r 1 0, for r > r 1
where f (r 0 , r 1 ) = (-6r 5 + 15(r 0 + r 1 )r 4 -10(r 2 0 + 4r 0 r 1 + r 2 1 )r 3 + 30r 0 r 1 (r 0 + r 1 )r 2 -30r 2 0 r 2 1 r + r 3 0 (r 2 0 -5r 1 r 0 + 10r 2 1 )).

Let us introduce a Taylor-Hood scheme which will be compared with our ϕ-FEM scheme. We introduce first the fitting mesh T f it h on domain Ω f it h as in the preceding section, cf. Fig. 3 (left), and adapt the Taylor-Hood FE space velocity space, cf. (67), as

V f it h = v h ∈ C( Ωfit h ) d : v h | T ∈ P 2 (T ) d ∀T ∈ T f it h , v h = 0 on Γ w ,
while keeping the pressure space (68) as before. Note that the velocity FE space does no longer contain any restrictions on the boundary part Γ f it h approximating Γ. We shall impose the boundary conditions there with the help of Lagrange multipliers, introducing the space

Λ f it h = µ h ∈ C( Ωfit h ) : µ h | F ∈ P 2 (F ) ∀F ∈ F f it h ,
where F f it h is the set of the boundary facets on Γ f it h . A fitted Taylor-Hood FE formulation is then written as: find

(u h , p h , λ h , U h , ψ h ) ∈ V f it h × M f it h × Λ f it h × R × R d ′ such that Ω f it h 2D(u h ) : D(v h ) - Ω f it h p h div v h - Ω f it h q h div u h + Γ f it h λ h • (v h -V h -ω h × r) + Γ f it h µ h • (u h -U h -ψ h × r) = Ω f it h ρ f gv h + mg • V h , (70) for all (v h , q h , µ h , V h , ω h ) ∈ V f it h × M f it h × Λ f it h × R d × R d ′ .
We present in Fig. 5 the velocity obtained with the standard Taylor-Hood FEM scheme (70). Such a velocity and the accompanying pressure, computed on a very fine fitted grid, will be used as the reference solution in the subsequent numerical experiments and will be denoted as u, p in what follows.

A comparison between the standard Taylor-Hood FEM (70) and ϕ-FEM (11) is presented in Figs. 6 and7 (we do not report the error in the particle rotation velocity whose exact value is 0 and which is accurately predicted by all the schemes up to machine precision; this can be attributed to the symmetry of our test case). Since the error is computed with respect to a fine reference solution solution living on a fine fitted mesh, the numerical solution computed by either (70) or [START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF] should be projected to this fine mesh in order to compute the errors. This is reflected in the legends of the convergence curves: Ω f ine stands for the fine fitted approximation there. Similarly to Remark 6, we impose the pressure mean by a Lagrange multiplier in both schemes, and adjust the additive constants properly when computing the errors in pressure.

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

Conclusions and perspectives

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

• ϕ-FEM is readily available for finite elements of any order (without the need of any additional higher order approximation of the geometry). • (corollary of the previous point) ϕ-FEM outperforms the standard fitted FEM on polygonal meshes if the order of piecewise polynomials is > 1 (unless one implements more complicated versions of FEM in order to treat properly the curvilinear boundary, such as the isoparametric FEM).

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

To counter-balance the last point, we should mention that the variational formulation at the base of ϕ-FEM can be rather complicated. Typically, it contains more terms than a fitted FEM or a CutFEM scheme for the same problem. The implementation, although straightforward, may thus require some extra programming efforts. As a reward, one easily gets a good description of the geometry which may result in smaller computing times than those for standard FEM, as highlighted in [START_REF] Cotin | ϕ-fem: an efficient simulation tool using simple meshes for problems in structure mechanics and heat transfer[END_REF].

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

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

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

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

• A practical construction of the levelset function, which should satisfy some assumptions summarized in Remark 3. In the present article, ϕ was given analytically, but in more realistic applications one will have to construct an appropriate ϕ h on the discrete level directly. A good candidate, in the vicinity of Γ would be the signed distance to Γ, for which there exist efficient libraries, but it will remain to extend it in a smooth way (again directly on the discrete level) far from Γ. As an alternative, we note that in other versions of ϕ-FEM, ϕ h may be required only on mesh cells near Γ, cf. [START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF][START_REF] Cotin | ϕ-fem: an efficient simulation tool using simple meshes for problems in structure mechanics and heat transfer[END_REF].

• The current construction of stabilization in ϕ-FEM relies heavily on the linearity of the governing equations. Indeed, the terms with prefactor σ in [START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF] reproduce the governing equations for both trial and test functions. If the equations are non-linear, one cannot do this since the formulation should remain linear in the test functions. Thus, going from Stokes to Navier-Stokes, for example, is not straightforward. Various options of linearization of the stabilization terms should be yet tested numerically and theoretically.

A ϕ-FEM for Stokes equations in a fixed domain

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

   -2 div(D(u)) + ∇p = f, in Ω, div u = 0, in Ω, u = u D , on Γ ∪ Γ w . (71)
Assume that u D and f are defined in the whole discrete domain Ω h . Inspired by the ϕ-FEM scheme for particulate flow given in [START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF], we can derive the following ϕ-FEM scheme for the non-homogeneous Stokes equations (71): find w h ∈ V h , p h ∈ M h satisfying

Ω h 2D(u D + ϕ h w h ) : D(ϕ h s h ) - G h (2D(u D + ϕ h w h ) -p h I)n • ϕ h s h (72)
-

Ω h p h div(ϕ h s h ) - Ω h q h div(u D + ϕ h w h ) + σ u J u (u D + ϕ h w h , ϕ h s h ) +σh 2 T ∈T Γ h T (-∆(u D + ϕ h w h ) + ∇p h ) • (-∆(ϕ h s h ) -∇q h ) +σ T ∈T Γ h T div(u D + ϕ h w h ) div(ϕ h v h ) = Ω h f ϕ h s h + σh 2 T ∈T Γ h T f (-∆(u D + ϕ h s h ) -∇q h ),
for all s h ∈ V h , q h ∈ M h . We now state our second main result for the Stokes equations:

Theorem 2. Suppose that Assumptions 1, 3, 4 and 5 hold true, the mesh T h is quasi-uniform. Let (u, p) ∈ H k+1 (Ω) d × H k (Ω) be the solution to (71) and (w h , p h ) ∈ V h × M h be the solution to (72). Denoting u h := ϕ h w h , it holds

|u -u h | 1,Ω∩Ω h + 1 ν |p -p h | 0,Ω∩Ω h ≤ Ch k (∥u∥ k+1,Ω + 1 ν ∥p∥ k,Ω )
with a constant C > 0 depending on the C 0 , m, M in Assumptions 1, 4, on the maximum of the derivatives of ϕ, on the mesh regularity, and on the polynomial degree k, but independent of h, f , and u. Moreover, supposing Ω ⊂ Ω h ∥u -u h ∥ 0,Ω ≤ Ch k+1/2 (∥u∥ k+1,Ω + 1 ν ∥p∥ k,Ω )

with a constant C > 0 of the same type.

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

B A glossary of geometrical notations.

Γ h the approximate interface: Γ h = {ϕ h = 0} T h the active mesh:

T h = {T ∈ T O h : T ∩ {ϕ h < 0} ̸ = ∅} T Γ h intersection of T h with Γ h T Γ h = {T ∈ T h : T ∩ Γ h ̸ = ∅} T Γ,

Figure 1 :

 1 Figure 1: Left: an example of geometry for the fluid Ω with a solid S inside; Right: the nonconforming active mesh T h on Ω h with its internal boundary G h .

Figure 2 :

 2 Figure 2: Example of T Γ h and F Γ h with the geometry given in Fig. 1 (left).

Case 1 TCase 2 T

 12 ∈ T Γ,ext h . We have then φE = -h so that φh ⩾ 0 and h ∇ φh ≥ ̸ ∈ T Γ,ext h , E is not shared with any cell from T Γ,ext h

4. 1

 1 Particular case of a fixed particle: Stokes equations We start by Stokes equations (71) in the domain Ω, as above, with the right-hand side such that the exact solution is as follows, cf. [14], u(x, y) = (cos(πx) sin(πy), -sin(πx) cos(πy)), p(x, y) = (y -0.5) cos(2πx) + (x -0.5) sin(2πy), taking ν = 1.

Figure 3 :

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

Figure 4 :

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

Figure 5 :

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

10 - 1 . 8 hFigure 6 :

 186 Figure6: Rates of convergence for the standard Taylor-Hood FEM scheme (70) and the ϕ-FEM scheme[START_REF] Duprez | A new ϕ-FEM approach for problems with natural boundary conditions[END_REF] in the case of particulate flows. The L 2 relative error of the velocity (left) and the H 1 relative error of the velocity (right).

10 - 1 . 8 Figure 7 :

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

  Assumption 1. The boundary Γ can be covered by open sets O i , i = 1, . . . , I and one can introduce on every O i local coordinates ξ 1 , . . . , ξ d with ξ d = ϕ such that all the partial derivatives ∂ α ξ/∂x α and ∂ α x/∂ξ α up to order k + 1 are bounded by some C 0 > 0. Thus, ϕ is of class C k+1 on O.We continue with assumptions on the mesh. To this end, we introduce an extended band of mesh elements near the boundary Γ, namely the submesh T Γ,ext The approximate interface G h can be covered by element patches {Π k } k=1,...,NΠ having the following properties:

	h h the cells which are neighbors and neighbors of neighbors of cells in T Γ with T Γ h ⊂ T Γ,ext h adding to T Γ h .	⊂ T h by
	Assumption 3. |∇ϕ| ≥ m, |∇ϕ h | ≥ m 2 on all the mesh cells in T Γ,ext h	, |ϕ| ≥ mh on T h \ T Γ,ext h	,
	and |∇ϕ h | ≤ M on Ω h with some m, M > 0.
	Assumption 4.
	Assumption 2. χ ∈ H k+1 (O), χ = 1 on S, χ = 0 on Γ w .

  26). Substituting this into (27) and combining with the Korn inequality (22) yields the announced estimate for ∥ϕh s h ∥ 0,Ω Γ h since ∥χ h ∥ 1,Ω Γh is bounded uniformly in h. The remaining part of the estimate follows by trace inverse inequalities as in[START_REF] Duprez | ϕ-FEM: a finite element method on domains defined by levelsets[END_REF] Lemma 3.5].

  {E (an internal facet of Th ) such that ∃T ∈ T Γ h and E ∈ ∂T }. G h the internal component of ∂Ω h , corresponding to the interface Γ: G h = ∂Ω h \ Γ w B hthe strip between Γ and G h :B h = Ω h \ Ω B Γ h the strip between Γ and Γ h Bh the strip between Γ h = {ϕ h = 0} and G h : Bh = {ϕ h > 0} ∩ Ω hWe also recall that the domain occupied by the active mesh T h is denoted by Ω h , i.e. Ω h := (∪ T ∈T h T )o . The same convention is applied to the submeshes T Γ h and T i h , giving respectively Ω Γ

	h T i h F Γ h	ext	T Γ h and the cells which are neighbors and neighbors of neighbors of cells of T Γ h in T h T h \ T Γ h F Γ h =

h

and Ω i h .

In practice, the geometrical setting may be slightly more complicated. The rigorous theoretical definitions of Ω h and G h will be given in[START_REF] Burman | Fictitious domain methods using cut elements: III. A stabilized Nitsche method for Stokes' problem[END_REF] and[START_REF] Ciarlet | Mathematical Elasticity: Volume I: three-dimensional elasticity[END_REF] and will be based on an approximation ϕ h to the levelset ϕ, rather than on ϕ itself. This may occasionally result in situations where some tiny portions of Ω lie outside Ω h so that G h slightly penetrates Ω. These technical details are not important for the forthcoming formal derivation of the FE scheme, while the rigorous proofs will be done assuming definitions[START_REF] Burman | Fictitious domain methods using cut elements: III. A stabilized Nitsche method for Stokes' problem[END_REF] and[START_REF] Ciarlet | Mathematical Elasticity: Volume I: three-dimensional elasticity[END_REF]. The actual implementation may introduce yet more geometrical approximations, as mentioned in Remark 1, which are not covered by our theory.

-2 div D(v) + ∇q = u -u h , in Ω, div v = 0, in Ω, v = V + ω × r, on Γ, v = 0,

https://github.com/michelduprez/phi-FEM-particulate-flows-Stokes.git or https://doi.org/10.5281/zenodo.6817135

This work was supported by the Agence Nationale de la Recherche, Project PhiFEM, under grant ANR-22-CE46-0003-01.

Note that the boundary term Γ u • (2D(v) -qI)n vanishes since u is a rigid body motion on Γ.

For the same reason, Γ (ϕw h + χ(U h + ϕ h × r) • (2D(v) -qI)n vanishes.