Mirco Ciallella
email: mirco.ciallella@ensam.eu

Thomas Milcent
email: thomas.milcent@ensam.eu

SEMI-IMPLICIT EULERIAN METHOD FOR THE FLUID STRUCTURE INTERACTION OF ELASTIC MEMBRANES *

Keywords: Eulerian elasticity, Fluid-structure interaction, Semi-implicit, Immersed boundary, Incompressible Navier-Stokes MSC codes. 65M08, 65M12, 65M85, 76D05, 74B20

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1. Introduction. Fluid-structure interaction problems arise in a wide range of applications and their numerical simulation is extremely challenging. In biomechanics, the correct modeling and computation of membranes immersed in an incompressible fluid is a key point to study biological capsules and cells, like red blood cells (RBCs). Capsules are in general modeled as a liquid drop protected by a thin elastic membrane and their use is widely spread in food, cosmetics and pharmaceutical industries. Simulations of such fluid-structure problems are expensive and require the correct coupling between the flow behavior and the structure behavior, which often comes with the imposition of some kind of interface conditions. One of the most classical approaches is the so-called Arbitrary Lagrangian-Eulerian (ALE) method [START_REF] Liu | A numerical study of insect flight[END_REF][START_REF] Legay | An eulerian-lagrangian method for fluid-structure interaction based on level sets[END_REF][START_REF] Sahin | An arbitrary lagrangian-eulerian formulation for the numerical simulation of flow patterns generated by the hydromedusa aequorea victoria[END_REF][START_REF] Donea | An arbitrary lagrangian-eulerian finite element method for transient dynamic fluid-structure interactions[END_REF][START_REF] Fanion | Deriving adequate formulations for fluidstructure interaction problems: from ale to transpiration[END_REF]. This approach consists in solving the fluid-structure problem on a moving mesh that follows the displacement of the interface, and enforce the coupling conditions right on the grid points that describe the interface. Although this philosophy may bring accurate results for some problems, its limitations are clear when dealing with large deformations or topology changes.

In [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF], Peskin introduced a new approach that simplifies the fluid-structure coupling for applications with structures modeled as lower dimensional manifolds of the problem. For immersed boundaries [START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Lee | An immersed interface method for incompressible navier-stokes equations[END_REF][START_REF] Mittal | Immersed boundary methods[END_REF], the coupling interface conditions reduce to a forcing term in the Navier-Stokes equations that modifies the flow behavior by taking into account the deformation of the immersed membrane. In this case, the structure is tracked using Lagrangian markers and the source term is spread on the fluid with a discretized Dirac delta. When large deformations are under consideration the need to insert or delete markers for stretching and shrinking of the membrane remains a delicate aspect, which requires the introduction of additional parameters to tune.

Both ALE and immersed boundary (IB) methods belong to the class of fronttracking methods that need an explicit definition and displacement of the interface. Fully Eulerian models for fluid-structure interaction have been recently introduced [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF] for computational purposes as they allow to simulate quite easily multi-dimensional problems with large deformations thanks to an implicit treatment of interfaces and deformations through level-set fields [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF][START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF]. In the last decade, Eulerian elasticity has drawn a lot of attention and its implementation turned out to be very successful for both incompressible [START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF][START_REF] Maitre | Applications of level set methods in computational biophysics[END_REF][START_REF] Richter | A fully eulerian formulation for fluid-structure-interaction problems[END_REF][START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF][START_REF] Deborde | Numerical simulations of the interaction of solitary waves and elastic structures with a fully eulerian method[END_REF][START_REF] Bergmann | An eulerian finite-volume approach of fluidstructure interaction problems on quadtree meshes[END_REF] and compressible [START_REF] Gorsse | A simple cartesian scheme for compressible multimaterials[END_REF][START_REF] Brauer | A cartesian scheme for compressible multimaterial models in 3d[END_REF][START_REF] Brauer | A cartesian scheme for compressible multimaterial hyperelastic models with plasticity[END_REF] applications. These models rely on the so-called backward characteristics which are defined as the inverse map of the deformation, and correlates the location of a material point to its initial position. As mentioned above, this vector is also defined with Eulerian variables that are advected by the Eulerian fluid velocity, and its spatial derivatives are used to compute the elastic stress tensor. Rather than following the backward characteristics, it was proposed in [START_REF] Sugiyama | A full eulerian finite difference approach for solving fluid-structure coupling problems[END_REF] an Eulerian approach based on following the evolution of the symmetrical left Cauchy-Green tensor through an inhomogeneous advection equation with extra terms depending on the velocity gradient. In this case, six components in three dimensions (three in two dimensions) must be advected, rather than three in three dimensions (two in two dimensions) for the backward characteristics.

Whether one deals with ALE, IB or fully Eulerian models for incompressible fluidstructure interaction problems, stability is a well-known issue of monolithic frameworks in which both the fluid and the solid are solved at the same time. For Eulerian models this was proven to be related to the velocity of elastic waves that (for stiff materials) are order of magnitude higher than that of the fluid, which leads to dramatic time step restrictions [START_REF] Brauer | A cartesian scheme for compressible multimaterial models in 3d[END_REF]. In order to overcome such limitations related to monolithic frameworks, several works have been dedicated to the derivation of implicit schemes [START_REF] Tu | Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods[END_REF][START_REF] Mayo | An implicit numerical method for fluid dynamics problems with immersed elastic boundaries[END_REF][START_REF] Stockie | Analysis of stiffness in the immersed boundary method and implications for time-stepping schemes[END_REF][START_REF] Newren | Enhancing the immersed boundary method: stability, volume conservation, and implicit solvers[END_REF]. However, due to the strong nonlinearity of the problem, these schemes require computationally expensive and cumbersome iterative approaches. With the goal of moving towards methods that are easier to introduce in different computational frameworks, semi-implicit or approximate implicit schemes are a very interesting and widespread research topic since their goal is to provide a less expensive and efficient alternative to standard implicit methods [START_REF] Tu | Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods[END_REF][START_REF] Stockie | Analysis of stiffness in the immersed boundary method and implications for time-stepping schemes[END_REF][START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF][START_REF] Fernández | Explicit robin-neumann schemes for the coupling of incompressible fluids with thin-walled structures[END_REF]. Similar interest in such schemes was also shown in the context of multi-phase flows for the computation of surface tension [START_REF] Hysing | A new implicit surface tension implementation for interfacial flows[END_REF][START_REF] Sussman | A stable and efficient method for treating surface tension in incompressible two-phase flow[END_REF] to avoid the time step restriction proportional to ∆x 3/2 .

In the context of fully Eulerian models for fluid-structure interaction, the attempts to relax the time step restrictions have been proposed in [START_REF] Ii | An implicit full eulerian method for the fluid-structure interaction problem[END_REF][START_REF] Cottet | A level set method for multiphase flows and fluidstructure interaction problems[END_REF] with methods that seem to be tailored to the application under consideration. In particular, the first reference [START_REF] Ii | An implicit full eulerian method for the fluid-structure interaction problem[END_REF] addresses the simulation of solid bulks immersed in an incompressible fluid by introducing a linear expansion for the nonlinear elastic stress with a Jacobian tensor, to improve the performances of the scheme initially proposed in [START_REF] Sugiyama | A full eulerian finite difference approach for solving fluid-structure coupling problems[END_REF]. Whereas in [START_REF] Cottet | A level set method for multiphase flows and fluidstructure interaction problems[END_REF] a semi-implicit approach for thin membranes subjected to area variation is developed by solving a diffusion equation to predict the interface position, in the context of the model proposed in [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF]. In this work, we aim at proposing a general way of deriving semi-implicit schemes for Eulerian models for fluid-structure interaction that can be directly applicable to general (nonlinear) elastic models without any need of linearization. This approach consists in discretizing the evolution equation of the stress tensor in a semi-implicit way, and writing the value of the stress at the new time level as a function of that at the previous time plus some additional terms. In order to show the generality and potential of this approach, we applied it to the twodimensional version of the full membrane model introduced in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF], which is capable of dealing with both area and shear variations. While the level-set approach proposed in [START_REF] Cottet | A level set method for multiphase flows and fluidstructure interaction problems[END_REF] only works for the elastic force that models the area variation introduced in [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF], our idea can be ideally applied to more general and nonlinear elastic models, because it depends on the deformation vector. Indeed, the model treated in this paper is the starting point to develop semi-implicit schemes for full membrane elasticity that includes both area and shear variations [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF]. The complete three-dimensional analysis will be tackled in a future work.

The paper is organized as follows. In Section 2, we recall the basic definitions to derive the fully Eulerian model for fluid-structure interaction focusing on backward characteristics, level-set and the membrane elastic model treated herein. In Section 3 we recall and summarize the main stability results obtained for a simplified onedimensional model, and derive our semi-implicit scheme in this framework. Here it is shown that the new scheme consists in adding a consistent dissipation to the system which allows us to prove its unconditional stability. In Section 4, we begin by describing the classical explicit coupling of the elastic source term and the general algorithm to solve the full model. Here, we also derive an evolution equation on the stress tensor that will be then used to design the semi-implicit approach. After that, we focus on the derivation of the semi-implicit coupling which is performed by introducing the discretized evolution equation of the stress in the momentum equation. Here we show that also the approach for the full model consists in adding a diffusion that depends on stiffness and area variation. In Section 5, we test the convergence properties of the discretized semi-implicit operators to show that they are second-order accurate in space. In Section 6, we provide numerical illustrations on the ability of the newly developed semi-implicit scheme to introduce an important relaxation on the time step restriction to simulate a membrane immersed in a simple shear flow. In particular, we focus on the influence of membrane stiffness and mesh refinement on the computational gain and quality of the solution. In Section 7, we conclude with a summary of the work and give some perspectives on aspects that need further investigation. Finally, in Appendix A we give a proof on the parameter chosen to measure the area variation, in Appendix B we give the proof of the evolution equation of the stress tensor, and in Appendix C we provide examples on how to discretize the semi-implicit operators on Cartesian staggered meshes.

Fully Eulerian Model.

2.1. Forward and backward characteristics. Let Ω 0 ⊂ R 2 be the reference configuration of a continuous medium and assume that this medium is deformed by a smooth map X : R + × R 2 -→ R 2 (the forward characteristics) to Ω t = X(t, Ω 0 ). A velocity field u : R 2 × R + -→ R 2 is naturally associated with X:

(2.1) ∂ t X(t, ξ) = u(X(t, ξ), t), X(0, ξ) = ξ, ξ ∈ Ω 0 .
We introduce the backward characteristics Y : R 2 × R + -→ R 2 by the formula Y (X(t, ξ), t) = ξ. The physical interpretation of Y (x, t) is the position at time 0 of a material particle lying in x at time t and moving at speed u. The derivative of this relation with respect to t and ξ in turn gives with (2.1)

Y (x, t) = ξ x = X(t, ξ) n(x, t) Initial configuration Ω 0 Deformed configuration Ω t
(2.2) ∂ t Y + (u • ∇ x )Y = 0, Y (x, 0) = x, x ∈ Ω t , and
(2.3) [∇ ξ X(t, ξ)] = [∇ x Y (x, t)] -1 , for x = X(t, ξ).
The relation (2.2) is the Eulerian equivalent of the characteristic equation (2.1). In addition, equation (2.3) allows to compute the gradient of the deformation in the Eulerian frame via Y . The next sections are devoted to the description of the mathematical model used to treat the membrane elastic deformations with the backward characteristics Y . To simplify the notation, we drop the subscript of the operator ∇ x .

2.2. Hyperelastic membrane models. The notations and results summarized in this section are detailed in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF] for the general three-dimensional model. We consider a surface Γ t = {x ∈ R 2 / ϕ(x, t) = 0} captured by a level set function ϕ : R 2 × R + -→ R and advected by the Eulerian velocity field u:

(2.4) ∂ t ϕ + u • ∇ϕ = 0.
The normal n(x, t) for x ∈ Γ t is then expressed in terms of the normalized gradient of the level set:

(2.5) n(x, t) = ∇ϕ(x, t) |∇ϕ(x, t)| .
To measure the deformations on the surface Γ t we introduce the tensor

(2.6) A = B - (Bn) ⊗ (Bn) (Bn) • n .
where B = [∇Y ] -1 [∇Y ] -T is the left Cauchy-Green tensor. Tensor A measures the surface deformations by projecting the deformations (measured by B) on the surface Γ t (represented locally by n). The vector n is an eigenvector of A associated to the eigenvalue 0 so det(A) = 0. In two dimensions, the other invariant (the trace) is used to define (2.7) Z = Tr(A).

Following the same approach in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF], it is possible to prove that Z measures the local area variation for two-dimensional problems (see Appendix A for the proof).

Remark 2.1. Using the backward characteristics Y , rather than only the level-set ϕ as in [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF], allows us to have the deformations that can also be used to measure the shear variation for tree-dimensional membranes. Indeed, the transport equation for the level-set function only records information on the surface area variations, and ignores any tangential component of u. This is due to the fact that ∇ϕ is normal to the interface.

In order to compute the elastic force, we then introduce the regularized membrane energy,

(2.8) E = Q E(Z) 1 ε ζ ϕ ε dx.
Here Q is a box containing the membrane, E(Z) is the constitutive law, ε the width of the interface and ζ is a cut-off function used to spread the interface near {ϕ = 0}. The associated force, computed with the principle of virtual power, is given by (2.9)

F = div (σ) δ ε (ϕ), where σ = E ′ (Z)Z(I -n ⊗ n)
where

δ ε (ϕ) = 1 ε ζ ϕ ε
is introduced to simplify the notation. The Evan-Skalak constitutive law is used in this article:

(2.10)

E ′ (Z) = K(Z -1),
where K is the elastic modulus of the membrane.

2.3.

Fluid-structure interaction model. The elastic membrane is immersed in a incompressible fluid modeled by the Navier-Stokes equations. The overall fully Eulerian model is given by (2.11)

         ρ(∂ t u + (u • ∇)u) + ∇p -div(2µD(u)) = F (ϕ, Y ), div(u) = 0, ∂ t ϕ + u • ∇ϕ = 0, ∂ t Y + (u • ∇)Y = 0,
where u is the fluid velocity, D(u) = ([∇u] + [∇u] T )/2 the strain rate tensor, ρ the density, p the pressure and µ the viscosity. Different values of density and viscosity, outside (1) and inside (2) the membrane, are taken into account by defining

ρ = H ϕ ε ρ 1 + 1 -H ϕ ε ρ 2 , µ = H ϕ ε µ 1 + 1 -H ϕ ε µ 2 ,
where H represents a smooth Heaviside function. The membrane force F is given by (2.9). These equations are completed with appropriate initial and boundary conditions that will be detailed in the numerical validation section.

The fluid-structure problem is governed by two parameters, the Reynolds number Re and the capillary number Ca, that can be recovered by rewriting the momentum equation in a non-dimensional form. Therefore, we introduce the reference variables φ * and the non-dimensional variables φ, which give the following operators and physical variables:

ρ = ρ * ρ, u = u * ū, p = p * p, µ = µ * μ, F = K * (L * ) 2 F , ∇ = ∇ L * .
By replacing all variables with their counterpart, and taking p * = µ * u * /L * , the nondimensional momentum equation reads,

(2.12) Re ρ ∂ t ū + (ū • ∇)ū + ∇p -div(μ([ ∇ū] + [ ∇ū] T )) = 1 Ca F ,
where Re = ρ * u * L * /µ * , and Ca = µ * u * /K * .

Reinitialization and extrapolation.

It is a well-known problem that the level-set field ϕ stops being a signed distance function after its time evolution with the associated advection equation. In this work, for its reinitialization, we follow the idea proposed in [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] and solve for a fictitious time τ

(2.13) ∂ τ ϕ + sgn(ϕ 0 )(|∇ϕ| -1) = 0.
Moreover, as discussed in [START_REF] Deborde | Numerical simulations of the interaction of solitary waves and elastic structures with a fully eulerian method[END_REF][START_REF] Bergmann | An eulerian finite-volume approach of fluidstructure interaction problems on quadtree meshes[END_REF] also the backward characteristics Y can get extremely distorted with time. Being the elastic force computed by means of Y and ϕ, if not treated, these distortions are likely to cause spurious disturbances close to the fluidstructure interface.

To summarize previous works on this topic, we perform the linear extrapolation proposed in [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF] for the outer characteristics, i.e. ϕ > 0,

∂ τ Y n + H ϕ ε (n • Y n ) = 0, with Y n = n • ∇Y, (2.14) ∂ τ Y + H ϕ ε (n • Y -Y n ) = 0. (2.15)
3. Stability analysis for a simplified one-dimensional model. In this section, we recall recent stability results in the context of explicit and implicit coupling for a simplified model, and we introduce the new semi-implicit approach in this framework. We consider here the same linearization and assumptions taken for similar models [START_REF] Bost | Linear stability analysis of a level set model of immersed elastic membrane[END_REF][START_REF] Fondaneche | Interaction fluide-structure dans un dispositif vasculaire actif[END_REF]. To summarize we set ρ = 1 and we decompose the solution as the small perturbation of a stationary solution (u 0 , Y 0 ) = (0, x). By avoiding the treatment of nonlinear terms we end up with the following linearized one-dimensional model:

(3.1)        ∂u ∂t -µ ∂ 2 u ∂x 2 = - K ε ∂ 2 Y ∂x 2 , ∂Y ∂t + u = 0.
Here we recall the stability conditions for such models when using a classical explicit coupling similarly to what was introduced in [START_REF] Bost | Linear stability analysis of a level set model of immersed elastic membrane[END_REF]. For simplicity we consider a classical centered discretization of second order derivatives, with grid size ∆x, with an implicit treatment for the viscous term:

(3.2)          u n+1 j -u n j ∆t -µ u n+1 j+1 -2u n+1 j + u n+1 j-1 (∆x) 2 = - K ε Y n j+1 -2Y n j + Y n j-1 (∆x) 2 , Y n+1 j -Y n j ∆t + u n+1 j = 0.
For this system, the stability condition was proven in [START_REF] Bost | Linear stability analysis of a level set model of immersed elastic membrane[END_REF] to be

(3.3) ∆t < µε + max(µε, √ Kε∆x) K .
Instead, when treating implicitly the elastic term, the discretized system reads:

(3.4)          u n+1 j -u n j ∆t -µ u n+1 j+1 -2u n+1 j + u n+1 j-1 (∆x) 2 = - K ε Y n+1 j+1 -2Y n+1 j + Y n+1 j-1 (∆x) 2 , Y n+1 j -Y n j ∆t + u n+1 j = 0.
Once again in [START_REF] Bost | Linear stability analysis of a level set model of immersed elastic membrane[END_REF], this system was proven to be unconditionally stable. However, for real applications, this scheme is very cumbersome to design and implement due to the full fluid-structure Navier-Stokes system with nonlinear source term and advection equations.

For this reason, following the same motivation of [START_REF] Cottet | A level set method for multiphase flows and fluidstructure interaction problems[END_REF], we introduce and study the stability analysis of the new semi-implicit scheme in this simplified framework. The new approach developed for the full nonlinear model is described in Section 4.2. We start by introducing an auxiliary variable Ỹ for the implicit elastic term,

           u n+1 j -u n j ∆t -µ u n+1 j+1 -2u n+1 j + u n+1 j-1 (∆x) 2 = - K ε Ỹ n+1 j+1 -2 Ỹ n+1 j + Ỹ n+1 j-1 (∆x) 2 , Y n+1 j -Y n j ∆t + u n+1 j = 0. (3.5)
For the simplified model, writing an evolution equation on the stress tensor, discretized in a semi-implicit manner, entails that Ỹ n+1 j = Y n j -∆tu n+1 j . This gives us an explicit part and an implicit one that depends on the second order derivatives of u,

Ỹ n+1 j+1 -2 Ỹ n+1 j + Ỹ n+1 j-1 (∆x) 2 = Y n j+1 -2Y n j + Y n j-1 (∆x) 2 -∆t u n+1 j+1 -2u n+1 j + u n+1 j-1 (∆x) 2 .
When introducing the previous relation in System (3.5), the overall system reads

(3.6)            u n+1 j -u n j ∆t -µ + ∆t K ε u n+1 j+1 -2u n+1 j + u n+1 j-1 (∆x) 2 = - K ε Y n j+1 -2Y n j + Y n j-1 (∆x) 2 , Y n+1 j -Y n j ∆t + u n+1 j = 0.
It should be noticed that this approach reduces to solving a PDE similar to the one arising from the explicit discretization, but with an additional consistent viscosity term ∆t K ε that depends on the time step and the elastic parameter of the structure. Proposition 3.1. The semi-implicit scheme (3.6) is unconditionally stable. Proof. With the assumption of periodic boundary conditions we perform a Von-Neumann stability analysis. Therefore, we decompose the solution (u, Y ) as Fourier expansions of the grid values,

u n j = ûn (k)e 2iπkj∆x , Y n j = Ŷ n (k)e 2iπkj∆x θ:=2πk∆x ------→ u n j = ûn (θ)e ijθ , Y n j = Ŷ n (θ)e ijθ .
We then obtain from (3.6)

       ûn+1 (θ) 1 + 4 ∆t (∆x) 2 µ + ∆t K ε sin 2 θ 2 = ûn (θ) + Ŷ n (θ) 4 ∆t (∆x) 2 K ε sin 2 θ 2 , ∆tû n+1 (θ) + Ŷ n+1 (θ) = Ŷ n (θ).
By setting

α θ = 1 + 4 ∆t (∆x) 2 µ + ∆t K ε sin 2 θ 2 , β θ = 4 ∆t (∆x) 2 K ε sin 2 θ 2 ,
we can recast the system in matrix form

α θ 0 ∆t 1 A θ u n+1 (θ) Ŷ n+1 (θ) = 1 β θ 0 1 B θ u n (θ) Ŷ n (θ) .
The eigenvalues λ 1 and λ 2 are the solution of the characteristic polynomial det(A -1 θ B θ -λI) = 0, which is det

1 α θ -λ β θ α θ -∆t α θ 1 -∆t β θ α θ -λ = α θ λ 2 -λ(1 + α θ -∆tβ θ ) + 1 = 0,
whose invariants are

λ 1 λ 2 = 1 α θ < 1, λ 1 + λ 2 = 1 + α θ -∆tβ θ α θ > 0.
If (λ 1 + λ 2 ) 2 < 4λ 1 λ 2 , the eigenvalues are complex conjugates, and their modulus is always less than 1. If (λ 1 + λ 2 ) 2 > 4λ 1 λ 2 , the eigenvalues are real and always positive.

In particular,

λ 1 = 1 + α θ -∆tβ θ + (1 + α θ -∆tβ θ ) 2 -4α θ 2α θ ≤ 1 + α θ + (1 + α θ ) 2 -4α θ 2α θ = 1,
which proves that the spectral radius of A -1 θ B θ is always less than one. Remark 3.2 (Stability of the fully nonlinear semi-implicit schemes). Observe that the above stability condition has been proven for a linear simplified one-dimensional model and it cannot be generally applied to the fully nonlinear fluid-structure problem. In particular, for realistic applications, a relaxation of the time step restriction is experienced thanks to the additional dissipation added depending on the local deformation of the media.

■ ■ ■ ■ ■ ■ ■ ■ ■ ▲ ▶ ▲ ▶ ▲ ▶ ▲ ▶ ▲ ▶ ▲ ▶ ▲ ▶ ▲ ▶ ▲ ▶ ▶ ▶ ▶ ▲ ▲ ▲ i -1 i i + 1 i -3 2 i -1 2 i + 1 2 i + 3 2 j -1 j j + 1 j -3 2 j -1 2 j + 1 2 j + 3 2
Fig. 2: Staggered grid with position of unknowns. ■: scalar field (p, ρ, Y , ϕ). ▶: vector field -component along x. ▲: vector field -component along y.

4. Numerical schemes. In this section we present how to develop both explicit and semi-implicit couplings for the full nonlinear model.

Explicit method.

In general, The equations (2.11) are discretized with finite volume schemes on a staggered grid (see Figure 4.1 for a two-dimensional configuration). It should be noticed that the force term F (ϕ, Y ) is here written as the divergence of a stress tensor multiplied by a cut-off function.

Let ∆t be the time step and u n , p n , ϕ n , Y n , ρ n , µ n the time discretization of the variables at t n = n∆t. The discretization is given by the projection method (see [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] for an overview on projection methods) for the Navier Stokes equations and an explicit scheme for the advection equations for the interface and backward characteristics:

Step 1 :

ρ n u ⋆ -u n ∆t + div(u n ⊗ u ⋆ ) + ∇p n (4.1) -div(2µ n D(u ⋆ )) = δ ε (ϕ n ) div(σ n ),
Step 2 : div ∆t

ρ n ∇ψ n+1 = div(u ⋆ ), (4.2)
Step 3 :

u n+1 = u ⋆ - ∆t ρ n ∇ψ n+1 , p n+1 = p n + ψ n+1 , (4.3)
Step 4 :

ϕ n+1 -ϕ n ∆t + u n+1 • ∇ϕ n = 0, Y n+1 -Y n ∆t + (u n+1 • ∇)Y n = 0. (4.4)
In Step 1 (4.1) a prediction of the velocity u ⋆ is computed with an Euler implicit scheme for the viscous term, and an explicit treatment of the elastic source term σ n . The convection term is treated following a semi-implicit approach, with a standard second-order discretization. In Step 2 (4.2), the Poisson equation for the pressure increment ψ n+1 is solved with appropriate boundary conditions [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. The resulting linear systems are solved with the GMRES algorithm of the HYPRE library [START_REF] Falgout | Hypre: a library of high performance preconditioners[END_REF][START_REF] Falgout | The design and implementation of hypre, a library of parallel high performance preconditioners, chapter in Numerical Solution of Partial Differential Equations on Parallel Computers[END_REF] with preconditioning. In Step 3 (4.3) the velocity is corrected to enforce the incompressibility condition and the pressure is updated. In Step 4 (4.4) the transport equations are discretized with an explicit Runge-Kutta3 scheme in time and a WENO5 [START_REF] Jiang | Efficient implementation of weighted eno schemes[END_REF] scheme in space for ϕ and Y .

For the cut-off function, we considered the following expression:

(4.5) ζ(r) = 1 2 (1 + cos(πr)), on [-1, 1] 0, elsewhere.
In our simulations ε is fixed at 2∆x which is the standard value used in the literature to spread the interface.

4.2. Semi-implicit method. In the context of fluid-structure interaction problems, a classical explicit treatment of the elastic source term introduces very strict time step restrictions, especially when the structure is stiff. This restriction leads to long and cumbersome simulations, particularly when working with incompressible flows for which long times are needed to simulate realistic configurations, and the speed of elastic waves is much faster than the fluid speed.

For such problems, a good approximation of the time step is given by the surface tension condition introduced in [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF], which limits the time step to ∼ ρ∆x 3 /K. Although capturing properly all elastic waves becomes crucial for compressible materials in which transients are under study [START_REF] Brauer | A cartesian scheme for compressible multimaterial models in 3d[END_REF], in incompressible simulations one is more interested in the long-time dynamics of deformable structures. For instance when dealing with very stiff structures, the admissible time step will be extremely low even if the structure will not deform that much.

For this reason here we present the new semi-implicit method for thin elastic membranes, following the approach introduced in Section 3 for which it was possible to prove unconditional stability in a simplified framework. The general idea consists in writing the evolution equation of the stress tensor, and then discretizing it in a semi-implicit way with respect to the Eulerian velocity u. Therefore, we introduce here an evolution equation for the nonlinear membrane stress tensor σ defined in (2.9). Proposition 4.1. Under the smoothness assumption made on u, the membrane stress tensor σ verifies

∂ t σ+u • ∇σ = (f ′ (Z)Z[∇u] : C) C+ f (Z) [∇u] T (n ⊗ n) + (n ⊗ n)[∇u] -2([∇u]n • n)(n ⊗ n) . (4.6)
where we defined f (Z) = E ′ (Z)Z, and C = I -n ⊗ n.

Proof. See Appendix B.

In order to develop the new method for the full model, we start by replacing the explicit term σ n by σ ⋆ , within the prediction step (4.1), (4.7)

ρ n u ⋆ -u n ∆t + div(u n ⊗ u ⋆ ) + ∇p n -div(2µ n D(u ⋆ )) = δ ε (ϕ n ) div(σ ⋆ ).
Once one has obtained the evolution equation for σ (4.6), we discretize it in a semiimplicit way for all the terms that depend on u as,

σ ⋆ -σ n ∆t + (u ⋆ • ∇)σ n = T (Z n , ∇u ⋆ , n) + f (Z n ) [∇u ⋆ ] T (n ⊗ n) + (n ⊗ n)[∇u ⋆ ] (4.8)
where

T (Z n , ∇u ⋆ , n) = (f ′ (Z n )Z n [∇u ⋆ ] : C) C -2f (Z n )([∇u ⋆ ]n • n)(n ⊗ n).
When replacing Equation (4.8) into (4.7) the fluid-structure momentum equation becomes

ρ n u ⋆ -u n ∆t + div(u n ⊗ u ⋆ ) + ∇p n + δ ε (ϕ n )∆t div (u ⋆ • ∇)σ n -T (Z n , ∇u ⋆ , n) + (4.9) div µ n I + δ ε (ϕ n )∆tf (Z n )(n ⊗ n) [∇u ⋆ ] + [∇u ⋆ ] T µ n I + δ ε (ϕ n )∆tf (Z n )(n ⊗ n) = δ ε (ϕ n ) div (σ n ) .
It should be noticed that the consistent terms coming from the evolution equation of σ add some viscosity to the classical one. Moreover, the considered elastic model also gives additional viscosity terms appearing in the T term of (4.9). All these additional discretized terms depend on the time step ∆t, the geometric features of the membrane described by n, and its area variation measured by Z. These new viscosity terms differ significantly from 2µD(u) because µ is scalar, while the new terms, like f (Z n )(n ⊗ n), are matrices. For this reason they are going to be referred to as tensorial viscosity terms.

Remark 4.2 (Computational costs of the semi-implicit method). It should be noticed that the resolution of the momentum equation for the explicit scheme already deals with a linear system to treat implicitly and semi-implicitly the strain rate tensor and the advection term, respectively. Therefore, the additional computational cost needed to perform the semi-implicit computations is given by adding the tensorial viscosity terms to the linear system. As shown in Appendix C a wider stencil is needed with respect to the classical viscous terms. These differences do not impact the costs that much also because only those elements around the zero level-set require the tensorial viscosity, meaning when δ ε (ϕ n ) ̸ = 0.

Manufactured solution.

For validation purposes in this section we analyze the convergence properties of the semi-implicit terms, discussed in Section 4.2, discretized with standard second-order accurate formulas. In Appendix C we give further details on how to discretize these operators.

The following equation is solved

10u + div (u • ∇)M -M [∇u] -[∇u] T M -([∇u] : M )M = S,
which reduces to solving a linear system M u h = S in a square domain

[0, 1] × [0, 1],
where S is computed by following the method of manufactured solutions and considering the exact solution (5.1) u = sin (πx) sin (πy) sin (πx) sin (πy) with M = 1 + sin (πxy) 1 + sin (πxy) 1 + sin (πxy) 1 + sin (πxy) .

Remark 5.1 (Manufactured solution). It should be noticed that an additional term 10u is added to the diagonal because solving only the semi-implicit operators gives rise to an ill-conditioned linear system which is never considered in practical 10 -2.6 10 -2.4 10 -2.2 10 -2 10 -1.8 10 -1.6 10 -1.4 10 -1.2 10 -5

10 -4 10 -3 10 -2 ∆x ∥u h -u ex ∥ L2 Numerical results
first order second order Fig. 3: Validation: convergence analysis with the discretized semi-implicit operators. u h and u ex represent the numerical and exact solution, respectively. applications. As a matter of fact, when solving the Navier-Stokes equation, we always have some additional terms to the diagonal given by the advection and diffusion parts of the system. This allows us to study the convergence of the discretized operators without having conditioning issues.

In Figure 5.1 the discretization error is computed with a L 2 norm on a set of nested meshes: 20 × 20, 40 × 40, 80 × 80, 160 × 160, 320 × 320. The convergence analysis shows that the discretized operators needed to design the semi-implicit scheme are indeed second-order accurate.

Numerical experiments.

In this section, we study the effectiveness of the proposed scheme by comparing it to the original explicit one for a thin membrane immersed in a simple shear flow. The simulation is setup by considering a membrane of radius a = 0.5 immersed within a computational domain [-4, 4]×[-2, 2] discretized with uniform Cartesian meshes. From (2.12) the dimensionless parameters are modified by taking the reference speed equal to a γ, and read

(6.1) Re = ρa 2 γ µ = 0.1, Ca = µa γ K ,
where the capillary number Ca, is modified according to the configuration chosen. For all cases, ρ = ρ 1 = ρ 2 and the angular speed γ are both set equal to 1. From the non-dimensional momentum Equation (2.12) it is possible to notice that when Ca tends to zero, the elastic source term becomes stiff. This also relates to what happens when dealing with low Mach flows or materials with the compressible Euler equations [START_REF] Abbate | An all-speed relaxation scheme for gases and compressible materials[END_REF][START_REF] Abbate | An asymptotic-preserving all-speed scheme for fluid dynamics and nonlinear elasticity[END_REF][START_REF] Thomann | Implicit relaxed all mach number schemes for gases and compressible materials[END_REF]. For these experiments, four capillary numbers are considered: 0.02, 0.01, 0.008 and 0.001. It should be noticed that when dealing with Ca = 0.001, the explicit scheme no longer provide any physical result (when taking µ 2 /µ 1 = 1) due to the stiff nature of the system, and therefore no comparison can be done with the new semi-implicit method. Therefore, we take the viscosity ratio equal to 10 for the stiffest membrane to perform an adequate comparison, and 1 for the others. This also shows the ability of the method to deal easily with a viscosity ratio that may differ from one, as it may be essential for practical applications [START_REF] Foessel | Influence of internal viscosity on the large deformation and buckling of a spherical capsule in a simple shear flow[END_REF]. The initial conditions are given by (6.2)

u 0 (x, y) = γy 0 , ϕ 0 (x, y) = x 2 + y 2 -a, Y 0 (x, y) = x y ,
to impose the linear shear velocity field, a circular membrane, and no pre-deformation on the initial configuration. We set Neumann boundary conditions on the left and on the right of the domain, and moving wall at the top and at the bottom. When a thin membrane is immersed into a shear flow, it starts elongating in the direction of the velocity vectors until it reaches a steady state position in which its shape no longer changes but its dynamics is still evolving. Indeed, in this configuration, the velocity vectors are perfectly tangent to the membrane and there is a constant stress applied on it, which makes it rotating around itself: this phenomenon is called tank-treading [START_REF] Keller | Motion of a tank-treading ellipsoidal particle in a shear flow[END_REF]. An overview of the initial and final position of the membrane for this experiment is shown in Figure 6.1. To compare properly the final shapes of the membrane, the reference time to stop the simulation is set for all cases equal to 1.5, which was observed to be enough reach the configuration described above.

As mentioned before, we are interested in studying the dynamics of the membrane when the new semi-implicit scheme is used, and specifically how the additional viscosity influences the computational cost and the quality of the solution when compared to the original explicit scheme. In particular, we focus on the simulation of the tank-treading phenomenon by changing the stiffness coefficient of the structure, ranging from a soft membrane to a very stiff one, and performing a grid refinement study to analyze its influence.

In Table 6.1 we present the ∆t needed to perform several experiments with varying capillary number, from the softest to the stiffest, on a set of three nested uniform meshes 128 × 64, 256 × 128, and 512 × 256. It can be noticed that in all cases the semi-implicit method provides higher maximum time steps, proving also numerically speaking that the additional dissipation introduces consistent viscosity where it is Table 1: Membrane in simple shear flow: comparisons between explicit (EX) and semi-implicit (SI) schemes in terms of maximum ∆t to reach a final time t f = 1.5. Ca = 0.001 Ca = 0.008 Ca = 0.01 Ca = 0.02

N x × N y EX SI EX SI EX SI EX SI
128 × 64 6.0 10 -3 6.0 10 -2 2.5 10 -2 1.0 10 -1 3.0 10 -2 1.0 10 -1 5.0 10 -2 1.5 10 -1

256 × 128 1.0 10 -3 2.5 10 -2 1.0 10 -2 8.5 10 -2 1.5 10 -2 9.5 10 -2 2.5 10 -2 1.0 10 -1 512 × 256 2.0 10 -4 5.0 10 -3 6.0 10 -3 6.0 10 -2 5.0 10 -3 5.0 10 -2 1.0 10 -2 7.5 10 -2 needed allowing us to use larger time steps. Moreover, the new method provides larger time steps when stiffer structures are under study, which is extremely promising especially for real applications. We believe that this numerical behavior can be explained by considering that when dealing with stiffer structure the elastic waves seem to become less relevant to capture the general dynamics of the problem. This entails larger time steps for stiff structures that do not deform that much, for which the CFL condition given by the elastic waves is very strict. Instead, since the deformations are more important for softer membranes, smaller time step ratios seem to be crucial to have good predictions.

The mesh refinement analysis also provides interesting results regarding the effect of the new method when refining the mesh. In particular, it is observed that a larger gain in terms of speedup is experienced on finer grids with respect to coarser grids, which is very promising also in view of the future development of the method for three-dimensional models. Table 6.2 shows the ratio between the semi-implicit time step ∆t SI and the explicit time step ∆t EX , showing its observed increasing trend when using finer grids and stiffer membranes.

The increased dissipation experienced for softer structures can also be observed when studying the membrane final shape once the final time is reached. In the following, we discuss the effect of the numerical scheme on the physics of the problem, which is extremely important to understand if the maximum admissible time steps given in Table 6.1 are enough to provide a good approximated solution. For this reason, we perform a set of numerical simulations and compare the final shape obtained through the original explicit scheme with those simulated with the semi-implicit scheme. Moreover, to investigate further the impact of the time step on the diffusion we compare the explicit result with a set of three simulations, for each problem, for which the time step has been adapted ranging from the explicit upper bound to the semi-implicit one. For all cases, the semi-implicit simulations are run with larger time steps with respect to the explicit one. For instance, in Figure 6.2 where we show the numerical results Fig. 5: Membrane in simple shear flow: final configuration on the finest mesh for Ca = 0.02 with both explicit (EX) and semi-implicit (SI) methods. For the solution computed with the explicit method, the maximum ∆t shown in Table 6.1 is taken.

For the semi-implicit method, a set of three simulations have been performed with ∆t ranging from the explicit bound to semi-implicit one shown in Table 6.1.

for the softest membrane Ca = 0.02 on the finest grid, we compare the zero level-set of the membrane obtained with the explicit method with ∆t = 1.0 10 -2 to the same curve given by the semi-implicit scheme run with ∆t = 2.5 10 -2 , 5.0 10 -2 , 7.5 10 -2 .

For all experiments we perform the same study and show the results on the finest grid.

As expected since the proposed scheme is consistent we obtain that all shapes follow the trace of the explicit one. For the very soft membrane shown in Figure 6.2 we experience a greater difference in the local error with respect to the explicit solution which increases for larger time steps. It is observed that, when deformations are important, damping the elastic waves with the tensorial viscosity terms comes with less accurate solutions. Moreover, we conclude that for soft structures CFL restrictions are not that important, and therefore the use of more sophisticated methods is not needed. Indeed, it is well-known that stricter restrictions occur for stiffer structures, and it is for those applications that an implicit or semi-implicit approach is crucial. In Figure 6.3, two intermediate test cases are considered with lower capillary numbers to study the effect of the new scheme in these configurations. For simplicity, the experiments performed with Ca = 0.01 and Ca = 0.008 will be discussed together due to similarities in the final shapes, but also for the maximum ∆t and gains shown in Tables 6.1 and 6.2. For both cases, the semi-implicit scheme introduces much less dissipation giving a final shape of the membrane that is much closer to the reference one with a higher speedup ∆t SI /∆t EX , especially when compared to the results Fig. 6: Membrane in simple shear flow: final configuration on the fine mesh for Ca = 0.01, 0.008 with both explicit (EX) and semi-implicit (SI) methods. For the solution computed with the explicit method, the maximum ∆t shown in Table 6.1 is taken. For the semi-implicit method, a set of three simulations have been performed with ∆t ranging from the explicit bound to semi-implicit one shown in Table 6.1.

observed in Figure 6.2. In particular, in this case a speedup of 10 is observed on the finest mesh. As expected from previous results and assumptions, the best results are obtained for the stiffest membrane simulated with Ca = 0.001. Figure 6.4 shows that all semi-implicit results with higher time step are almost superimposed on the reference shape computed with ∆t = 2.0 10 -4 . In particular, for this experiment, we observe a speedup of 10 on the coarse mesh and 25 on the fine one that translates into a dramatic reduction of the computational time which also comes with a very precise prediction of the membrane shape. Although elastic waves are extremely fast for this experiment, the additional viscosity seem to be able to dampen them and relax the time step restriction, still providing a good dynamics of the membrane with almost no additional costs. We conclude by pointing out that the proposed approach behaves the best for the stiffest structure on the finest mesh. This is critical for both real applications that often concern stiff materials and computational efficiency, which one would like to improve especially when refining the mesh.

7. Summary and future perspectives. We presented a novel approach for the construction of semi-implicit schemes for Eulerian models for fluid-structure interaction. This approach is very interesting for the simulation of incompressible flows in which the speed of elastic waves is much higher than that of the fluid causing the need of using extremely low time steps when compared to the simulation time. The focus of this work was on hyperelastic membranes and on the two-dimensional version of the model for full membrane elasticity introduced in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF]. The tensorial viscosity terms provide additional but consistent dissipation that activates only on the few cells where the membrane interface is spread. This was also demonstrated for a simplified Fig. 7: Membrane in simple shear flow: final configuration on the fine mesh for Ca = 0.001 with both explicit (EX) and semi-implicit (SI) methods. For the solution computed with the explicit method, the maximum ∆t shown in Table 6.1 is taken.

For the semi-implicit method, a set of three simulations have been performed with ∆t ranging from the explicit bound to semi-implicit one shown in Table 6.1.

one-dimensional model in Section 3, where this additional term allowed us to prove that the scheme is unconditionally stable. As shown in Section 4, in the full nonlinear fluid-structure system these consistent terms appear in the momentum equation (4.9) and depend on the time step, the stiffness parameter, and the local deformation of the membrane. In Section 6 a set of numerical experiments is performed for a membrane immersed in a simple shear flow. It was shown that the dependences mentioned before influence the performances of the new scheme, which is observed to be more efficient and accurate for stiffer structures and finer grids. In particular, the numerical method proposed herein performed best on the experiment with the lowest capillary number showing a remarkable speedup of 10 on the coarse mesh and 25 on the fine one while preserving a very good agreement with the reference solution.

The perspectives of this work are several and range from the development of the method for other hyperelastic models to its application to advance simulations of biological systems.

Regarding the investigations on the method, in this work we considered the twodimensional version of the full membrane model introduced in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF], and future studies will also focus on the development of the method for the full three-dimensional model. Following the same reasoning, new evolution equations will be written to take into account the modeling of shear variation. Therefore, this will introduce viscosity terms that depend on both the area and shear variation. Other models will also be considered to study the effect of this approach on surface tension (special case of the membrane We first focus on the equation related to f (Z). Considering the time derivative of f (Z) gives us

(B.2) ∂ t f (Z) + u • ∇f (Z) = f ′ (Z)Z [∇u] : C,
where the last step is obtained with the equation on Z taken from [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF], recalled here

(B.3) ∂ t Z + u • ∇Z = Z [∇u] : C.
To compute the equation on C, an equation on n is needed. The evolution of n is computed by deriving Equation (2.5) as follows Appendix C. Discretization of semi-implicit operators.

∂ t n = 1 |∇ϕ| ∂ t (
In this appendix, we briefly describe how we discretize the semi-implicit operators for the staggered mesh framework used by the incompressible Navier-Stokes solver. In particular, for a two-dimensional problem, up to 16 terms need to be discretized: ∂ i (m ∂ j u k ), with i, j, k = x, y, for both u x and u y located on vertical and horizontal faces, respectively (see Figure 4.1). Indeed, since the position of the x-and y-components of a vector field are located in different points the discretization of the same operator on one direction is going to be different to that on the other direction. For compactness, here we do not present the discretization of all terms but, we mainly focus on few of them to give an idea on how to discretize these operators. There are of course several ways to discretize such terms but as long as they are consistent no remarkable difference is observed. For simplicity a uniform discretization is considered, meaning that each cell has a ∆x × ∆y area. To simplify the notation, only in this appendix we take the following abuse of notation: u = u x and v = u y .

For both components, we discretize two trivial operators and a third trickier one. For the x-direction (on vertical faces) we consider the discretization of ∂ x (m ∂ x u), ∂ y (m ∂ y u), and ∂ y (m ∂ y v): ∂x(m ∂xu)| i-1/2,j = 1 ∆x mi,j u i+1/2,j -u i-1/2,j ∆x -mi-1,j u i-1/2,j -u i- (m i+1,j-1 + m i+1,j + m i,j + m i,j-1 ), and the velocity component u as u * i+1,j-1/2 = 1 4 (u i+3/2,j-1 + u i+3/2,j + u i+1/2,j + u i+1/2,j-1 ), u * i,j-1/2 = 1 4 (u i+1/2,j-1 + u i+1/2,j + u i-1/2,j + u i-1/2,j-1 ),

u * i-1,j-1/2 = 1 4
(u i-1/2,j-1 + u i-1/2,j + u i-3/2,j + u i-3/2,j-1 ).

Fig. 1 :

 1 Fig. 1: Forward and backward characteristics.

Fig. 4 :

 4 Fig. 4: Membrane in simple shear flow: initial and final configurations for soft membranes.

  ∇ϕ) -∂ t (∇ϕ) • ∇ϕ |∇ϕ| ∇ϕ |∇ϕ| .This relation together with the evolution equation on the gradient of ϕ (2.6)∂ t (∇ϕ) + u • ∇(∇ϕ) = -[∇u] T ∇ϕ, gives us (B.4) ∂ t n + u • ∇n = -[∇u] T n + (([∇u]n) • n)n.With the identities (Ab) ⊗ c = A(b ⊗ c) and b ⊗ (Ac) = (b ⊗ c)A T , it is possible to write the equation on C, which reads(B.5) ∂ t C + u • ∇C = [∇u] T (n ⊗ n) + (n ⊗ n)[∇u] -2(([∇u]n) • n)(n ⊗ n).Using (B.5) and (B.2) in (B.1) we get the evolution equation (4.6) stated in the proposition.

Table 2 :

 2 Membrane in simple shear flow: time step ratio ∆t SI /∆t EX with data taken from Table6.1.

		Ca = 0.001 Ca = 0.008 Ca = 0.01 Ca = 0.02
	Nx × Ny		∆t SI /∆t EX
	128 × 64	10	4	3.33	3
	256 × 128	25	8.5	6.33	4
	512 × 256	25	10	10	7.5

  identifies an interpolated variable. In this case, we can approximate the multiplicative factor m as + m i,j+1+ m i-1,j+1 + m i-1,j ), + m i,j + m i-1,j + m i-1,j-1 ), + v i,j+3/2 + v i-1,j+3/2 + v i-1,j+1/2 ), + v i,j+1/2 + v i-1,j+1/2 + v i-1,j-1/2 ), + v i,j-1/2 + v i-1,j-1/2 + v i-1,j-3/2 ).For the y-direction (on horizontal faces) we consider the discretization of∂ y (m ∂ y v), ∂ x (m ∂ x v), and ∂ x (m ∂ x u): ∂ y (m ∂ y v)| i,j-1/2 = 1 ∆y m i,j v i,j+1/2 -v i,j-1/2∆y -m i,j-1 v i,j-1/2 -v i,j-3/2 ∆y , ∂

	3/2,j i-1/2,j-1/2 ∆x -m * u i-1/2,j -u i-1/2,j-1 , ∆y i-1/2,j-1/2 -m * i,j-1/2 u i-1/2,j+1 -u i-1/2,j ∆y v * i-1/2,j+1 -v * i-1/2,j ∆y (m i,j m * ∂y(m ∂yu)| i-1/2,j = 1 ∆y m * i-1/2,j+1/2 ∂y(m ∂yv)| i-1/2,j = 1 ∆y m * i-1/2,j+1/2 i-1/2,j+1/2 = 1 4 i-1/2,j-1/2 = 1 4 (m i,j-1 and the velocity component v as v * i-1/2,j+1 = 1 4 (v i,j+1/2 v * i-1/2,j = 1 4 (v i,j-1/2 v * i-1/2,j-1 = 1 4 (v i,j-3/2 2 = 1 ∆x m * i+1/2,j-1/2 u * i+1,j-1/2 -u * ∆x -m * i-1/2,j-1/2 u * i,j-1/2 -u * i-1,j-1/2 , ∆x where we can approximate the remaining m term as v m * m * 1 i+1/2,j-1/2 = 4	,

* i-1/2,j -v * i-1/2,j-1 ∆y , where (•) * x (m ∂ x v)| i,j-1/2 = 1 ∆x m * i+1/2,j-1/2 v i+1,j-1/2 -v i,j-1/2 ∆x -m * i-1/2,j-1/2 v i,j-1/2 -v i-1,j-1/2 ∆x , ∂ x (m ∂ x u)| i,j-1/

Acknowledgments. MC would like to acknowledge the support of Antoine Lemoine for the development of this project within the open-source massively parallel software notus CFD (https://notus-cfd.org/).

This work was funded by the ANR CAPSEULERIAN-FSI.

model with area variation) and, anisotropic elastic bodies [START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF]. About the possible applications, the proposed method is very promising for the simulation of biological capsules and red blood cells (RBCs), modelled as thin membranes characterized by shear effects, and their interaction with complex geometries [START_REF] Hu | Characterizing the membrane properties of capsules flowing in a square-section microfluidic channel: Effects of the membrane constitutive law[END_REF][START_REF] Häner | Sorting of capsules according to their stiffness: from principle to application[END_REF].

Appendix A. Z is the local area variation. In this appendix, we prove that Z introduced in (2.7) measures the local area variation as it was firstly introduced in [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF]. In particular, in [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF] it was proven that for incompressible media the gradient of the level-set |∇ϕ| transports information about the local area variation. Afterwards, in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF] a new theory was introduced with the goal of having a single formulation (for both compressible and incompressible media) to deal with both area and shear variations. For this reason, two parameters have been introduced for three-dimensional membranes as functions only of the backward characteristics Y .

Here, we prove that to capture the area variation for two-dimensional membranes a different parameter is needed with respect to [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF], and we also show how to relate it to the theory introduced in [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF].

Proposition A.1. Given that J = det(∇ ξ X) = det(∇Y ) -1 is the volume ratio, different than one for general compressible media,

Proof. The local area variation defined in (2.