Alexis Dalmon

Kalyani Kentheswaran

Guillaume Mialhe

Benjamin Lalanne

Sébastien Tanguy
email: tanguy@imft.fr

Guillaume Mialhe

Benjamin Lalanne

Fluids-membrane interaction with a full Eulerian approach based on the level set method

Keywords:

come

Introduction

ABSTRACT

A fully Eulerian approach to predict fluids-membrane behaviours is presented in this paper. Based on the numerical mode( proposed by li et al. (2012), we present a sharp methodology to account for the jump conditions due to hyperelastic membranes. The membrane is considered infinitely thin and is represented by the level set method. lts deformations are obtained from the transport of the components of the left Cauchy Green tensor throughout time. Considering the linear or a hyperelastic material law, the surface stress tensor is computed and gives the force exerted by the membrane on the surrounding fluids. The membrane force is taken into account in the Navier-Stokes equations as jump conditions on the pressure and on the velocity derivatives by imposing suitable singular source terrns in cells crossed by the interface. To prevent stability issues, an extension algorithm has been developed to remove the normal derivatives of the scalar fields specific to the membrane. ln particular, a subcell resolution at the interface of the extrapolated variable is proposed for increasing the accuracy of the extension algorithm. These improvements are validated by comparing our numerical results with benchmarks from the literature. Moreover, a new benchmark is proposed for fluids with both different viscosities and different densities to target applications where a gas and a liquid phase are separated by a membrane.

The fluids-membrane interaction study of this paper is part of a global project on propellant sloshing in satellite tanks. The latter phenomenon happens during a satellite manoeuvre and can be a major disturbance of the stability. The tanks contain Jiquid propellant and gas to maintain a sufficient pressure within the tank. During a manoeuvre, inertial forces lead to a motion of the fluids and thus of the centre of mass. This generates disturbing forces and torques on the whole structure which may deteriorate the quality of satellite imaging. Considering simple tanks, numerical methods have been developed in our home-made code DIVA (Dynamics of Interface for Vaporisation and Atomisation) to model propellant sloshing in micro-gravity conditions [START_REF] Dalmon | Direct numerical simulation of a bubble motion in a spherical tank under external forces and microgravity conditions[END_REF][START_REF] Lepilliez | On two-phase flow solvers in irregular domains with contact line[END_REF]. A parametric study has been done on typical rotational manoeuvres exerted by satellites in space and has been validated by comparisons with data from the FLUIDICS (FLUId DynamICs in Space) experiment, performed in the International Space Station (ISS) [START_REF] Dalmon | Direct numerical simulation of a bubble motion in a spherical tank under external forces and microgravity conditions[END_REF][START_REF] Dalmon | Comparison between the fluidics experiment and direct numerical simulations of fluid sloshing in spherical tanks under microgravity conditions[END_REF][START_REF] Mignot | Fluid dynamic in space experiment[END_REF]. The work presented in this paper is a first step towards the extension of the sloshing study to diaphragm tanks for which a hyperelastic membrane separates the liquid propellant and the gas. The modelling of the interaction between the membrane and the fluids within the tank is crucial. It must be predicted accurately and thus numerical developments are required to do so.

The fluid-membrane interaction is a challenging problem to solve numerically. Peskin [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF] developed the Immersed Boundary Method (IBM) to predict incompressible flows with moving elastic boundaries. These first applications were the modelling of blood flows in the heart [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF] and has been extended to many biological problems cited in [START_REF] Peskin | The immersed boundary method[END_REF]. With this method, the Navier-Stokes equations are solved on a Cartesian mesh and the membrane is described by Lagrangian markers. The force exerted by the membrane is deduced from the position of the markers and is interpolated onto the Cartesian mesh using Dirac delta functions. The Navier-Stokes equations are then solved with the forcing term corresponding to the elastic contribution from the membrane. Finally, the predicted fluid velocity is used to update the location of the marker points defining the membrane. This methodology is repeated for each time step. The smoothed Dirac functions induce a numerical smearing of the membrane in the fluid grid which may affect the accuracy of the method.

The Immersed Interface Method (IIM) of Leveque and Li [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficient and singular sources[END_REF] replaces the interpolated forcing term due to the elastic membrane with sharp jump conditions. This method has been introduced for elliptic equations and extended to Stokes flows [START_REF] Leveque | Immersed interface methods for stokes flow with elastic boundaries or surface tension[END_REF] and incompressible viscous flows [START_REF] Li | The immersed interface method for the Navier-Stokes equations with singular forces[END_REF]. The artificial smearing of the elastic force in the fluid grid generates spurious velocities at the membrane and a smoothing of the pressure field where a discontinuous jump must appear. With the IIM, the normal component of the elastic force is enforced through a jump condition on the pressure when solving the Poisson equation. The tangential component of the force induces jumps in the derivatives of the velocity across the membrane.

The Material Point Method (MPM) describes both the solid and the fluid phases using Lagrangian markers in the whole computational domain [START_REF] Hu | Material point method applied to fluid-structure interaction (fsi)/aeroelasticity problems[END_REF]. This method has been used in the case of fluid-membrane interaction [START_REF] York | Fluid-membrane interaction based on the material point method[END_REF] and allows large deformations of the linear elastic membrane.

Full Eulerian approaches do not use Lagrangian particles to follow the membrane motion and to compute the elastic forces. Cottet and Maitre [START_REF] Cottet | A level set method for fluid-structure interaction with immersed surfaces[END_REF] use the level set method to follow the motion of an elastic membrane immersed in a fluid. Without the reinitialisation algorithm of the level set function [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF], the variation of the level set gradient can be directly related to the membrane stretching. The authors use this property to enforce the elastic material law on the zero level set and to integrate the resulting force in the Cartesian mesh using delta functions. However, this method only gives the membrane stretching and limits the complexity of the membrane model. Moreover, the membrane is still artificially thickened by the delta functions.

Another full Eulerian method consists in defining the components of the deformation tensor of the membrane as scalar variables in the whole computational domain [START_REF] Ii | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach[END_REF]. This approach is based on the work of Sugiyama [START_REF] Sugiyama | A full Eulerian finite difference approach for solving fluid-structure coupling problems[END_REF] which developed a full Eulerian fluid-structure interaction model working with hyperelastic solid bodies. The position of each phase, fluid and solid, is described by the volume-of-fluid function (VOF) [START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF] and only one set of governing equations is solved in the whole domain. In the fluid region, the stress tensor contains the pressure term and the viscous stress tensor, in the solid region, the stress tensor is deduced from the solid deformation and the hyperelastic material law. The left Cauchy-Green deformation tensor is updated throughout time thanks to a transport equation in the whole computational domain. This allows to follow the deformation of the hyperelastic solid in a Eulerian manner. Then, the solid stress tensor can be computed following the hyperelastic material law and enforced in the Cartesian mesh. Ii et al. [START_REF] Ii | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach[END_REF] extend this method to fluid-membrane interactions, which means that the solid region is reduced to a codimension-one subspace. Following the methodology of Barthes-Biesel and Rallison [START_REF] Barthes-Biesel | The time-dependent deformation of a capsule freely suspended in a linear shear flow[END_REF], the deformation tensor of the membrane is defined as a solid deformation tensor projected onto the tangent plane of the membrane. The normal projection of the membrane deformation tensor is not considered because the membrane material is supposed to be incompressible. The membrane stresses are computed thanks to the deformation tensor and the hyperelastic material law and the resulting force is integrated in the Navier-Stokes equations with delta functions. Nevertheless, this method presents losses of accuracy due to the numerical dissipation of the smoothed Dirac function. Moreover, instability issues may appear over long time periods because the values of the deformation tensor far from the membrane may evolve chaotically with the fluid velocity and may influence the computation at the membrane.

In this paper, we propose to improve the full Eulerian method of Ii et al. [START_REF] Ii | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach[END_REF] by adding some aspects of the sharp methodology of the IIM. Moreover, we extend this model to handle different fluids on each side of the membrane. In section 2, the modelling of two-phase flows and the membrane model are described. The numerical methods implemented in the code are explained in section 3. More particularly, the discretisation of the jump conditions from the Ghost Fluid Method [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF][START_REF] Liu | A boundary condition capturing method for Poisson's equation on irregular domain[END_REF] is detailed. Furthermore, an extension algorithm based on a subcell resolution has been developed to reduce the numerical instability by removing the spurious normal derivatives. The validation of the methodology is done in section 4 with several benchmarks of growing complexity from the literature: from the simple case of a bubble rising due to a surface tension gradient until the much more complex cases of a stretched membrane separating different fluids and a capsule immersed in a shear flow. Finally, the paper is concluded by some remarks and perspectives in section 5.

Eulerian fluid membrane model

Two-phase flow model

We consider a domain with a boundary ∂ which contains two different fluids defined by + and -such as = + ∪ -. The interface between the two fluids regions is denoted and its outward normal vector is n. Each fluid is incompressible and Newtonian and follows the Navier-Stokes equations.

∇ • u = 0, (1)
ρ ∂u ∂t + (u • ∇) u = -∇ p + ∇ • 2μ D , ( 2
)
with u = (u, v, w) the velocity field, t the time, ρ the fluid density, μ the fluid viscosity, p the pressure and D the rate of deformation tensor defined as

D = ∇u + ∇u T 2 . ( 3
)
Considering the entire domain , special care must be taken at the interface between the two fluids regions. The following jump conditions must be accounted for,

[ρ] = ρ + -ρ -, (4)
[μ] = μ +μ -, [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF] [n • σ ] = f , [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF] with σ the stress tensor and f the local force density at the interface. In the case of two-phase flows with surface tension, the local force density becomes f = γ κn with γ the surface tension and κ = ∇ • n the mean curvature of the interface.

When a membrane separates the two fluids, the local force is defined in section 2.2.

Considering the entire domain , we define the density and viscosity field as

ρ(x) = ρ -+ (ρ + -ρ -)H(x), (7)
μ(x) = μ -+ (μ +μ -)H(x), [START_REF] Leveque | The immersed interface method for elliptic equations with discontinuous coefficient and singular sources[END_REF] with H the Heaviside function defined as H(x) = 0 in -and H(x) = 1 in + .

The Navier-Stokes equations in the entire domain can therefore be written as

∇ • u = 0, (9)
ρ ∂u ∂t + (u • ∇) u = -∇ p + ∇ • 2μ D + δ f , ( 10
)
with δ the Dirac function located at the interface.

Membrane model

In this section, the large deformable membrane model of Barthes-Biesel and Rallison [START_REF] Barthes-Biesel | The time-dependent deformation of a capsule freely suspended in a linear shear flow[END_REF] is described. The local force density of the membrane on the fluids is computed with a full Eulerian formulation. This formulation is inspired by the work of Sugiyama et al. [START_REF] Sugiyama | A full Eulerian finite difference approach for solving fluid-structure coupling problems[END_REF] for 3D solid problems and adapted to the membranes by Ii et al. [START_REF] Ii | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach[END_REF].

In a stress-free state, each particle of a solid is defined by its coordinate vector X. We denote by x(X, t) the position, at time t, of a particle located in X in the stress-free state. The link between the current state and the stress-free state is the deformation gradient tensor F defined as

F = ∂x ∂X . ( 11
)
The material time derivative of the deformation gradient tensor F follows d

F dt = (∇u s ) F , ( 12
)
with u s the velocity of the membrane particles. This velocity vector is obtained by extending u from its values at the membrane toward the normal direction. The extension method is defined in section 3.6. Now, we consider the special case of a membrane which thickness is very small compared to its other dimensions. Therefore, we will neglect the thickness and represent the membrane as a surface in a 3D space. Each fibre dX of the membrane in the reference state belongs to the tangent plane of the membrane. Similarly, each deformed fibre dx belongs to the tangent plane of the membrane in the current state. This means that only the components perpendicular to the normal direction of the membrane must be considered by the deformation gradient tensor. Let n R be the outward normal of the membrane in the reference state and n the outward normal in the current state, the previous condition can be expressed respectively as F s .n R = 0 and n. F s = 0 with F s the surface deformation gradient tensor of the membrane. This tensor can then be written

F s = P F P R , ( 13
)
with P = Īn ⊗ n the surface projection tensor in the current state and P R = Īn R ⊗ n R in the reference state. The local deformation of the membrane can be obtained with the surface left Cauchy-Green tensor,

Bs = F s F T s . ( 14
)
The Cauchy-Green tensor Bs is symmetric and its eigenvalues λ 2 1 and λ 2 2 correspond to the square of the two principal strains of the membrane in its tangent plane and the third eigenvalue zero corresponds to the deformation in the normal direction which has no existence because of the projections performed previously. The scalar invariants of the surface left Cauchy-Green tensor are

I 1 = tr( Bs ) = λ 2 1 + λ 2 2 , ( 15
)
I 2 = 1 2 tr( Bs ) 2 -tr( B2 s ) = λ 2 1 λ 2 2 , ( 16
)
I 3 = det( Bs ) = 0. (17)
This strain tensor, considering ( 13) and ( 14), can be written thanks to an intermediate tensor Ḡ s , Bs = P Ḡs P with Ḡs = F P R F T .

(

Considering ( 12) and ( 18), the material time derivative of Ḡ s is given by d Ḡs dt = (∇u s ) Ḡs + Ḡs (∇u s ) T .

(

) 19
In the case of hyperelastic materials, the stress tensor σ can be written according to the surface strain energy function W (I 1 , I 2 , I 3 ) which depends on the left Cauchy-Green tensor and its scalar invariants. In the case of a membrane, Barthes-Biesel and Rallison [START_REF] Barthes-Biesel | The time-dependent deformation of a capsule freely suspended in a linear shear flow[END_REF] show that the surface stress tensor σ s is,

σ s = 2 √ I 2 ∂ W ∂ I 1 Bs + I 2 ∂ W ∂ I 2 P . ( 20
)
Several strain energy functions [START_REF] Gent | Forms for the stored (strain) energy function for vulcanized rubber[END_REF][START_REF] Mooney | A theory of large elastic deformation[END_REF][START_REF] Ogden | Large deformation isotropic elasticity -on the correlation of theory and experiment for incompressible rubberlike solids[END_REF][START_REF] Yeoh | Characterization of elastic properties of carbon-black-filled rubber vulcanizates[END_REF] exist to describe different hyperelastic behaviours. All functions depending on the scalar invariants of the Cauchy-Green tensor can be considered with this methodology. We will consider in this study the strain energy function of the neo-Hookean solid [START_REF] Treloar | The elasticity of a network of long chain molecules[END_REF], adapted to the membrane case by considering the incompressibility of the membrane as

W = E s 6 (I 1 + 1 I 2 -3), ( 21
)
with E s the surface elastic modulus defined as E s = Eh R with E the Young's modulus of the membrane and h R the thickness of the membrane in the reference state. The surface stress tensor becomes for this material model,

σ s = E s 3 √ I 2 Bs - P I 2 . ( 22
)
Finally, the local force density exerted by the membrane is given by the surface divergence of the membrane stress tensor, f = ∇ s • σ s , [START_REF] Ogden | Large deformation isotropic elasticity -on the correlation of theory and experiment for incompressible rubberlike solids[END_REF] with ∇ s the surface gradient operator defined by ∇ s = P .∇. Considering the distance property of the level set function, which implies (n • ∇)n = 0, and the geometric properties of the membrane, some simplifications can be made, P σ s = σ s , n • σ s = 0, [START_REF] Yeoh | Characterization of elastic properties of carbon-black-filled rubber vulcanizates[END_REF] and the local force becomes f = ∇ • σ s . [START_REF] Treloar | The elasticity of a network of long chain molecules[END_REF] The local force density can be divided into two components, respectively toward the normal and the tangent directions of the membrane such as [START_REF] Huber | Direct numerical simulation of nucleate pool boiling at large microscopic contact angle and moderate Jakob number[END_REF] because the stress tensor is symmetric, and

f n = f • n = -σ s : ∇n = -σ s ∇ • n,
f τ = P f = P (∇ • σ s ). ( 27
)

Numerical methods

This section describes specific numerical methods implemented to model the membrane behaviour and its interaction with the surrounding fluids. These developments have been integrated in the home-made code DIVA. This solver is based on several numerical methods dedicated to the computation of two-phase flows. The DIVA code can also consider liquid-vapour phase change [START_REF] Huber | Direct numerical simulation of nucleate pool boiling at large microscopic contact angle and moderate Jakob number[END_REF][START_REF] Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: application to the Leidenfrost effect[END_REF][START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF][START_REF] Tanguy | Benchmarks and numerical methods for the simulation of boiling flows[END_REF][START_REF] Urbano | Direct numerical simulation of nucleate boiling in micro-layer regime[END_REF] and compressible flows [START_REF] Huber | A time splitting projection scheme for compressible two-phase flows. Application to the interaction of bubbles with ultrasound waves[END_REF]. Complex geometry can be accounted for by using the irregular domain method proposed in [START_REF] Ng | An efficient fluid-solid coupling algorithm for single-phase flows[END_REF] for single phase flows and extended to two-phase flows in irregular domains in [START_REF] Lepilliez | On two-phase flow solvers in irregular domains with contact line[END_REF]. In this study, the interface between the two fluids corresponds to the membrane and is represented by the level set method [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]. The Ghost-Fluid method [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF][START_REF] Gibou | A level set based sharp interface method for the multiphase incompressible Navier-Stokes equations with phase change[END_REF][START_REF] Gibou | A second-order-accurate symmetric discretization of the Poisson equation on irregular domains[END_REF][START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF][START_REF] Liu | A boundary condition capturing method for Poisson's equation on irregular domain[END_REF][START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF] is used to consider the sharp jump conditions at the membrane. It should be pointed out here that Ghost Fluid Method is first order accurate for imposing jump conditions, but some recent works have proposed extensions to second order accuracy as in [START_REF] Cisternino | A parallel second order Cartesian method for elliptic interface problems[END_REF][START_REF] Guittet | Solving elliptic problems with discontinuities on irregular domains -the Voronoi interface method[END_REF][START_REF] Hou | A numerical method for solving variable coefficient elliptic equation with interfaces[END_REF].

Interface tracking method

The level set method [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] is used to track the infinitely thin membrane throughout time. We define the scalar field φ which corresponds to the signed distance from the membrane . Each fluid region corresponds to the sign of the level set function as + = {x : φ(x) > 0}, -= {x : φ(x) < 0} and the membrane corresponds to the zero level set = {x : φ(x) = 0}.

The motion of the membrane is updated by solving the following transport equation ∂φ ∂t

+ u • ∇φ = 0. ( 28
)
To keep the signed distance property of the level set function throughout time, we consider the reinitialisation algorithm [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF]. The following equation is solved iteratively to correct the distance between each level set

⎧ ⎨ ⎩ ∂d ∂τ + sign(φ)(1 -||∇d||) = 0, d(τ = 0) = φ, ( 29
)
with τ a fictitious time along which the reinitialised distance function d is corrected to maintain the distance property ||∇d|| = 1. The smoothed sign function sign(φ) is defined in [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF][START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF].

The level set method allows us to compute the geometric properties of the membrane such as its outward normal using simple differencing, n = ∇φ ||∇φ|| , [START_REF] Urbano | Direct numerical simulation of nucleate boiling in micro-layer regime[END_REF] and thus, the projection tensor P = Īn ⊗ n can be computed in the whole domain.

The spatial derivatives of the transport equation and the reinitialisation algorithm are computed using the WENO-Z scheme [START_REF] Borges | An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[END_REF]. The temporal scheme of the transport equation is the second order TVD Runge-Kutta scheme. The reinitialisation algorithm is solved at the end of each time step.

Two-phase flow solvers

In this section, we describe the methodology to solve the Navier-Stokes equations with two different ways to consider the membrane contribution. The first one is the "delta" formulation for which the elastic force is a source term in the right-hand side. The second method integrates the elastic force through sharp jump conditions on the pressure and on the velocity derivatives.

The "delta" formulation

The "delta" formulation [START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF][START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] solves the Navier-Stokes equations for two-phase flows by smoothing the interface on an artificial thickness. The method is based on a projection method inspired by Chorin [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] in the case of single phase flows. For two-phase flows, the method we use is derived from the Ghost-Fluid viscous Conservative Method with an Implicit scheme of Lepilliez et al. [START_REF] Lepilliez | On two-phase flow solvers in irregular domains with contact line[END_REF], inspired by the work of Sussman et al. [START_REF] Sussman | A sharp interface method for incompressible two-phase flows[END_REF]. First, given a velocity field u n at the time t n = n t, an intermediate velocity u * is computed without considering the pressure term,

ρ n+1 u * -t∇ • 2μ n+1 D * = ρ n+1 u n -t u n • ∇ u n + tδ ε f n . ( 31
)
The viscous term is considered implicitly to avoid its time step restriction. This leads to the resolution of a large linear system where the three components of the velocity are coupled. The advection term is computed with WENO-Z schemes [START_REF] Borges | An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[END_REF] and the viscous term with central differencing schemes [START_REF] Lepilliez | On two-phase flow solvers in irregular domains with contact line[END_REF]. The local force f n is multiplied by the smoothed Dirac distribution δ ε defined below and added to the right-hand side. As the resulting matrix is diagonally dominant, the system is solved using a few steps of the Gauss-Seidel algorithm.

Then, the pressure field p n+1 serves as the scalar potential function of the Hodge decomposition which satisfies the following Poisson equation

∇ • ∇ p n+1 ρ n+1 = ∇ • u * t , ( 32
)
with homogeneous Neumann boundary conditions on ∂ . The resolution of this equation is done with the Black Box

Multigrid method [START_REF] Dendy | Black box multigrid[END_REF][START_REF] Maclachlan | Fast and robust solvers for pressure-correction in bubbly flow problems[END_REF] to reduce the computation time.

Finally, the velocity field u n+1 is defined as the projection of the intermediate velocity u * onto the divergence-free space

u n+1 = u * - t ρ n+1 ∇ p n+1 . ( 33
)
In the framework of the Whole Domain Formulation, the density and viscosity fields are updated with the level set function using a smoothed Heaviside distribution [START_REF] Gibou | A second-order-accurate symmetric discretization of the Poisson equation on irregular domains[END_REF] with H ε (x) the smoothed Heaviside distribution defined as

ρ(x) = ρ -+ (ρ + -ρ -)H ε (x), (34) μ(x) = μ -+ (μ + -μ -)H ε (x),
H ε (x) = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 0 if φ(x) < -ε, 1 2 1 + φ(x) ε + sin(π φ(x)/ε) π if |φ(x)| < ε, 1 if φ(x) > ε, ( 36
)
with ε the fictitious thickness of the membrane, which is equal to two or three times the size of a cell [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF]. Similarly, the smoothed Dirac distribution δ ε (x) is defined as the derivative of the smoothed Heaviside function,

δ ε (x) = ⎧ ⎪ ⎨ ⎪ ⎩ 0 if |φ(x)| > ε, 1 2ε 1 + cos πφ(x) ε if |φ(x)| < ε. ( 37
)

The sharp formulation

Let us consider that n is the normal vector to the membrane and e τ 1 , e τ 2 are two orthogonal vectors belonging to the tangent plane of the membrane. The no-slip condition on both sides of the membrane enforces [u] = 0 and gives

∂u ∂e τ 1 = 0 and ∂u ∂e τ 2 = 0. ( 38
)
The incompressibility condition ∇ • u = 0 in both fluids ensures that there is no jump at the membrane [∇ • u] = 0 and gives the following relation from [START_REF] Tan | An immersed interface method for the incompressible Navier-Stokes equations with discontinuous viscosity across the interface[END_REF] ∂u ∂n

• n = 0. ( 39
)
The density and viscosity fields are piecewise constant and their jumps at the membrane correspond to equations ( 4) and [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF]. With the sharp formulation, the local density force corresponding to the membrane contribution is expressed as three primary jump conditions:

[p] = 2 μ ∂u ∂n • n + f • n = 2[μ] ∂u ∂n • n + f n , (40)
μ ∂u ∂n

• e τ 1 + μ ∂u ∂e τ 1 • n = -P f • e τ 1 = -f τ • e τ 1 , (41)
μ ∂u ∂n

• e τ 2 + μ ∂u ∂e τ 2 • n = -P f • e τ 2 = -f τ • e τ 2 .
(42)

These jump conditions have been proven by Xu and Wang [START_REF] Xu | Systematic derivation of jump conditions for the immersed interface method in three-dimensional flow simulation[END_REF] in single phase flows and extended to piecewise constant viscosity by Tan et al. [START_REF] Tan | An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane[END_REF][START_REF] Tan | An immersed interface method for the incompressible Navier-Stokes equations with discontinuous viscosity across the interface[END_REF]. In this study, we only consider these primary jump conditions which account for all the physical phenomena of interest for this study. It is noticeable that higher accuracy numerical methods have been achieved in [START_REF] Tan | An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane[END_REF][START_REF] Tan | An immersed interface method for the incompressible Navier-Stokes equations with discontinuous viscosity across the interface[END_REF][START_REF] Xu | Systematic derivation of jump conditions for the immersed interface method in three-dimensional flow simulation[END_REF] in the simpler case of linear elasticity, by imposing further secondary jump conditions on the second order velocity derivatives and first order pressure derivatives. However, these numerical methods have not been generalised to more complex configurations such as hyperelastic membranes.

Following the Ghost-Fluid Conservative viscous Method (GFCM), Lalanne et al. [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF] show that only the normal component of the force must be considered in the projection and the correction steps because the viscous component is already taken into account in the predictor step of the projection method. The remaining jump condition on the pressure is then

[p] = f • n = f n . ( 43
)
From equations [START_REF] Cisternino | A parallel second order Cartesian method for elliptic interface problems[END_REF], the jump conditions ( 41) and ( 42) can be written as

μ ∂u ∂n • e τ 1 + [μ] ∂u ∂e τ 1 • n = -f τ • e τ 1 , (44)
μ ∂u ∂n • e τ 2 + [μ] ∂u ∂e τ 2 • n = -f τ • e τ 2 . ( 45
)
Following a similar approach to that of Lalanne et al. [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF], the source terms due to the tangential components of the elastic force are considered as jumps that are only enforced on the first term:

μ ∂u ∂n • e τ k with k = [1, 2]
. The second terms of equations ( 44) and ( 45) will be zero if there is no viscosity jump, even in the case where the tangential component of the local force is not zero. This approach still stands with a viscosity jump and has been used in [START_REF] De Langavant | Level-set simulations of soluble surfactant driven flows[END_REF] where several simulations with surfactants were performed. Consequently, the splitting of the jump conditions allows a simple numerical implementation of the source terms in the cells cut by the membrane. Indeed, the jump conditions on the tangent components of the viscosity-scaled velocity gradient in the first step of the GFCM can be taken into account by imposing

μ ∂u ∂n • e τ k = -f τ • e τ k , (46)
with k = [START_REF] Dalmon | Direct numerical simulation of a bubble motion in a spherical tank under external forces and microgravity conditions[END_REF][START_REF] Lepilliez | On two-phase flow solvers in irregular domains with contact line[END_REF]. Because the jump of the normal component of the viscosity-scaled velocity gradient is only due to the jump on the viscosity, the following source term is considered in the first step of the algorithm

μ ∂u ∂n = -P f = -f τ . (47)
In section 4.1, a validation of this methodology to deal with jumps of both tangential and normal stresses is presented by comparing simulations and the theoretical solution of a droplet rising in a surface tension gradient in the Stokes regime. The second terms of equations ( 44) and ( 45) are zero in this case, with or without viscosity jump. As the simulations match to the theoretical solution, this comparison demonstrates the relevance of the presented approach. For a general case with

nonzero terms [μ] ∂u ∂e τ 1 • n and [μ] ∂u ∂e τ 2
• n and viscosity jump, we assume that all the terms of equations ( 44) and [START_REF] Maclachlan | Fast and robust solvers for pressure-correction in bubbly flow problems[END_REF] are well taken into account. To our knowledge, a demonstration for this general case has never been provided and further investigation would be appreciated. It appears from this discussion that the key-point is to introduce a consistent numerical approximation of the source term resulting from the jump condition. Finally, the jump condition will be automatically satisfied if a consistent approximation of the source term is added at the right place. Other general considerations on this specific point can be found in Appendix B. From this, some analogies can be found between jump condition formulations and δ formulations for two-phase flows Navier-Stokes equations.

To solve the Navier-Stokes equations for two-phase flows with the Ghost Fluid viscous Conservative Method, the same projection method inspired by [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF] is considered with jump conditions to describe the membrane force. First, the intermediate velocity is computed using the same implicit scheme for the viscous term,

⎧ ⎪ ⎨ ⎪ ⎩ ρ n+1 u * -t∇ • 2μ n+1 D * = ρ n+1 u n -t u n • ∇ u n , μ ∂u ∂n = -f τ , ( 48
)
with the source term due to the tangent component of the singular force imposed as a jump on the viscous stress. It is noteworthy that, unlike the delta formulation, the prediction step of the projection method only contains this tangent component, since the normal component is accounted for in the projection and correction steps. The analogy between sharp methods using jump conditions and the "delta" formulation is thoroughly explained in [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF]. The linear system resulting from ( 48) is a single linear system where all the velocity components are coupled. This allows a fully implicit discretization of the viscous term without additional time step constraints due to viscosity jump and it can be solved easily with a Gauss-Seidel algorithm as pointed out by Lepilliez et al. in [START_REF] Lepilliez | On two-phase flow solvers in irregular domains with contact line[END_REF]. The numerical implementation of the jump conditions on the viscous stress tensor is detailed in section 3.5.

Then the Poisson equation is solved with

∇ • ∇ p n+1 ρ n+1 = ∇ • u * t , ( 49
)
with the following pressure jump condition accounting for the normal component of the local force

[p] = f n , ( 50
)
and the following jump condition on the pressure normal gradient

n • ∇ p n+1 ρ n+1 = μn • u n ρ n+1 , ( 51
)
details of which can be found in Appendix B. The numerical discretization of this singular source term is identical to the one described by Liu et al. in [START_REF] Liu | A boundary condition capturing method for Poisson's equation on irregular domain[END_REF] and is detailed in section 3.4. These equations based on a jump condition formulation can also be rewritten in a similar way with singular source terms instead of jump condition. It will give a mathematically identical system that will express as [START_REF] Min | A second order accurate level set method on non-graded adaptive Cartesian grids[END_REF] to compute the velocity field in the predictor step. Next, the pressure field can be computed with the appropriate jump condition with the following equation

ρ n+1 u * -t∇ • 2μ n+1 D * = ρ n+1 u n -t u n • ∇ u n + tf τ δ ,
∇ • ∇ p n+1 ρ n+1 = ∇ • u * t + ∇ • f n nδ ρ n+1 , ( 53
)
which is identical to equations ( 49) and [START_REF] De Langavant | Level-set simulations of soluble surfactant driven flows[END_REF]. Indeed, the term

u n+1 = u * - t ρ n+1 (∇ p n+1 -f n nδ ). ( 54
)
By considering that modifications of the numerical schemes accounting for jump conditions act as sharp approximations of singular source terms, it can be understood that the two presented formulations are identical. See for instance [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF] for more details on this specific point.

Discretisation of the membrane model

Considering the Cartesian mesh, the standard MAC grid is used: the scalar variables are located in the centre of the mesh cells and the components of the velocity field are staggered at the cell faces in each direction. Unlike in [START_REF] Ii | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach[END_REF], all the components of the different solid tensors are located at the centre of the mesh cells in this study. The location of the different variables in a 2D mesh cell is illustrated in Fig. 1.

At each time step, the position of the interface is updated using [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF]. The six components of the intermediate strain tensor Ḡs in 3D are then computed using equation (19) Ḡn+1

s = Ḡn s -t (u s n • ∇) Ḡn s -(∇u s n ) Ḡn s -Ḡn s (∇u s n ) T . ( 55
)
The advection term is computed with the WENO-Z scheme [START_REF] Borges | An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws[END_REF] and the velocity derivatives are obtained with central differencing schemes. For example, in 2D, the advection term corresponding to the first component of the intermediate strain tensor is computed as

(u s n • ∇)G n 11 = u s i, j ∂ G 11 ∂x i, j + v s i, j ∂ G 11 ∂ y i, j (56)
with u s i, j = 1 2

u s i-1/2, j + u s i+1/2, j and v s i, j = 1 2 v s i, j-1/2 + v s i, j+1/2 . ( 57
)
Considering now the derivatives of the x-component of the velocity,

∂u s ∂x i, j = u s | i+1/2, j -u s | i-1/2, j
x and

∂u s ∂ y i, j = u s | i, j+1/2 -u s | i, j-1/2 y , (58)
with

u s i, j±1/2 = 1 4 u s i-1/2, j + u s i+1/2, j + u s i-1/2, j±1 + u s i+1/2, j±1 . ( 59
)
The computation of the components of the surface left Cauchy-Green tensor and the stress tensor at the cell centres follows equations [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] and [START_REF] Liu | A boundary condition capturing method for Poisson's equation on irregular domain[END_REF].

The components of the surface density force f are computed at the centre of the faces of the cells following [START_REF] Treloar | The elasticity of a network of long chain molecules[END_REF]. Considering the sharp methodology and the force decomposition [START_REF] Huber | Direct numerical simulation of nucleate pool boiling at large microscopic contact angle and moderate Jakob number[END_REF][START_REF] Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: application to the Leidenfrost effect[END_REF], the normal component of the force is computed at the centre of the cell and the tangent components are staggered at the cell faces. Considering a 2D example, the discretisation of the normal component of the force, following [START_REF] Huber | Direct numerical simulation of nucleate pool boiling at large microscopic contact angle and moderate Jakob number[END_REF], is

f n i, j = σ 11 i, j ∂n x ∂x i, j + σ 12 i, j ∂n x ∂ y i, j + σ 21 i, j ∂n y ∂x i, j + σ 22 i, j ∂n y ∂ y i, j = σ 11 i, j n x | i+1, j -n x | i-1, j 2 x + σ 12 i, j n x | i, j+1 -n x | i, j-1 2 y + σ 21 i, j n y | i+1, j -n y | i-1, j 2 x + σ 22 i, j n y | i, j+1 -n y | i, j-1 2 y . ( 60
)
The tangent component of the force in the x-direction from [START_REF] Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: application to the Leidenfrost effect[END_REF] gives

f τ • e x i+1/2, j = P 11 i+1/2, j ∂σ 11 ∂x i+1/2, j + ∂σ 12 ∂ y i+1/2, j + P 12 i+1/2, j ∂σ 12 ∂x i+1/2, j + ∂σ 22 ∂ y i+1/2, j = P 11 i+1/2, j σ 11 | i+1, j -σ 11 | i, j x + σ 12 | i+1/2, j+1/2 -σ 12 | i+1/2, j-1/2 y + P 12 i+1/2, j σ 12 | i+1, j -σ 12 | i, j x + σ 22 | i+1/2, j+1/2 -σ 22 | i+1/2, j-1/2 y , ( 61
)
with P kl i+1/2, j =

σ kl i, j + σ kl i+1, j + σ kl i, j+1 + σ kl i+1, j+1 . ( 62
)
With the "delta" formulation, the discretisation of the elastic force corresponds to equation [START_REF] Lee | An immersed interface method for incompressible Navier-Stokes equations[END_REF] with Ī instead of P to consider all the components of the force. This method is identical in the other directions and can easily be extrapolated in 3D.

Jump condition on the pressure

The pressure jump is taken into account following the Ghost Fluid Method [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF][START_REF] Liu | A boundary condition capturing method for Poisson's equation on irregular domain[END_REF]. The Poisson equation on the pressure (53) can be written as ∇ • (β∇ p) = R H S, [START_REF] Ramanujan | Deformation of liquid capsule enclosed by elastic membranes in simple shear flow: large deformation and the effect of fluid viscosities[END_REF] with β = 1/ρ the diffusion coefficient. The equation contains jump conditions on the pressure denoted a = f n and on the diffusion coefficient [β] .

The discretisation of the 2D Poisson equation at point (i, j) gives

β i+1/2, j p i+1, j -p i, j x -β i-1/2, j p i, j -p i-1, j x x + β i, j+1/2 p i, j+1 -p i, j y -β i, j-1/2 p i, j -p i, j-1 y y = R H S i, j + g i, j , ( 64
) with (β i+1/2, j , β i-1/2, j , β i, j+1/2 , β i, j-1/2
) the harmonic averages of the diffusion coefficient at the centre of the cell borders.

If the interface crosses the right mesh segment [x i, j , x i+1, j ],

β i+1/2, j = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ β + β - β -θ R + β + (1 -θ R )
if φ i, j < 0 and φ i+1, j > 0,

β + β - β + θ R + β -(1 -θ R ) if φ i, j > 0 and φ i+1, j < 0, ( 65
) with θ R = |φ i+1, j | |φ i+1, j | + |φ i, j | . ( 66
)
Similarly, if the interface crosses the left mesh segment [x i-1, j , x i, j ],

β i-1/2, j = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ β + β - β -θ L + β + (1 -θ L ) if φ i, j < 0 and φ i-1, j > 0, β + β - β + θ L + β -(1 -θ L ) if φ i, j > 0 and φ i-1, j < 0, (67)
with

θ L = |φ i-1, j | |φ i-1, j | + |φ i, j | . ( 68
)
The g i, j term in (64) corresponds to the jumps enforced when the interface crosses at least one of the four neighbouring mesh segments. The latter are denoted by the superscripts R, L, T and B which correspond respectively to the right, left, top and bottom borders. As a general rule,

g i, j = g R i, j + g L i, j + g T i, j + g B i, j . ( 69
)
Each one of these values exists only if the membrane crosses the mesh segments and then equals

g R i, j = ± β i+1/2, j a R x 2 , g L i, j = ± β i-1/2, j a L x 2 , g T i, j = ± β i, j+1/2 a T y 2 , g B i, j = ± β i, j-1/2 a B y 2 , ( 70
)
with ± corresponding to the opposite sign of φ i, j and

a R = a i, j θ R + a i+1, j (1 -θ R ), a L = a i, j θ L + a i-1, j (1 -θ L ), ( 71
)
a T = a i, j θ T + a i, j+1 (1 -θ T ), a B = a i, j θ B + a i, j-1 (1 -θ B ). ( 72
)
The extension of these schemes to 3D problems is straightforward.

Jump condition on the velocity derivatives

The discretisation of each component of the viscous-stress tensor is described in this section with a specific emphasis on how to enforce suitable jump conditions. Considering a 2D example, the divergence of the viscous-stress tensor gives

∇ • (2μ D) = ⎛ ⎝ ∂ ∂x 2μ ∂u ∂x + ∂ ∂ y μ ∂u ∂ y + ∂ v ∂x ∂ ∂x μ ∂u ∂ y + ∂ v ∂x + ∂ ∂ y 2μ ∂ v ∂ y ⎞ ⎠ . ( 73
)
With [START_REF] Tan | An immersed interface method for the incompressible Navier-Stokes equations with discontinuous viscosity across the interface[END_REF], the jump condition of each velocity derivative can be expressed as

μ ∂u ∂x = n x μ ∂u ∂n = -n x (f τ • e x ) , μ ∂u ∂ y = n y μ ∂u ∂n = -n y (f τ • e x ) , (74) μ ∂ v ∂x = n x μ ∂ v ∂n = -n x f τ • e y , μ ∂ v ∂ y = n y μ ∂ v ∂n = -n y f τ • e y . ( 75
)
Considering the first term of the viscous-stress tensor in (73), let β = 2μ be the diffusion coefficient and b = -n x (f τ • e x )

be the jump condition on the velocity derivative computed at the right and left cell borders. The discretisation of this term at point (i

+ 1/2, j) gives ∂ ∂x β ∂u ∂x i+1/2, j = β i+1, j u i+3/2, j -u i+1/2, j x -β i, j u i+1/2, j -u i-1/2, j x x -h i+1/2, j (76)
with β i, j and β i+1, j the harmonic averages of the diffusion coefficient at the centre of the cells which are computed with the same methodology as in equations ( 65) to (68).

The jump conditions to be enforced are h i+1/2, j = h R i+1/2, j + h L i+1/2, j if the membrane crosses respectively the interval

[x i+1/2, j , x i+3/2, j ] or [x i-1/2, j , x i+1/2, j ] defined as h R i+1/2, j = ± β i+1, j b R θ R β ± x and h L i+1/2, j = ∓ β i, j b L θ L β ± x (77) with ± the opposite sign of φ i+1/2, j , ∓ the sign of φ i+1/2, j , b R = b i+1/2, j θ R + b i+3/2, j (1 -θ R ) and b L = b i+1/2, j θ L + b i-1/2, j (1 -θ L ), ( 78
) with θ R = |φ i+3/2, j | |φ i+3/2, j | + |φ i+1/2, j | and θ L = |φ i-1/2, j | |φ i-1/2, j | + |φ i+1/2, j | . ( 79
)
The same method is applied for all the second derivatives of the same variable. Special care must be taken in the case of mixed derivatives. Considering the last term of the divergence of the viscous stress tensor in the x-direction, the discretisation at point (i

+ 1/2, j) gives ∂ ∂ y β ∂ v ∂x i+1/2, j = β i+1/2, j+1/2 v i+1, j+1/2 -v i, j+1/2 x -β i+1/2, j-1/2 v i+1, j-1/2 -v i, j-1/2 x y -h i+1/2, j . ( 80
)
The diffusion coefficient is now β = μ and the jump condition on the mixed derivative is b = -n x f τ • e y . The harmonic averages of the diffusion coefficient are similarly computed following equations ( 65) to (68). The jump conditions h i+1/2, j = h T i+1/2, j + h B i+1/2, j exist only if the membrane crosses respectively the interval [x i+1, j+1/2

, x i, j+1/2 ] or [x i+1, j-1/2 , x i, j-1/2 ]
and are computed similarly as in equation ( 77),

h T i+1/2, j = ± β i+1/2, j+1/2 b T θ T β ± y and h B i+1/2, j = ∓ β i+1/2, j-1/2 b B θ B β ± y , ( 81
)
with ± corresponding to the opposite sign of φ i+1/2, j and ∓ the same sign than φ i+1/2, j .

Considering that the membrane crosses the top mesh segment, the sign of the level set function at the calculation point

(i + 1/2, j) gives the expression of b T and θ T . If φ i+1/2, j φ i, j+1/2 < 0, the computation of the jump condition becomes b T = b i+1, j+1/2 θ T + b i, j+1/2 (1 -θ T ) with θ T = |φ i, j+1/2 | |φ i, j+1/2 | + |φ i+1, j+1/2 | . ( 82
) Conversely, if φ i+1/2, j φ i, j+1/2 > 0, the jump condition becomes b T = b i+1, j+1/2 (1 -θ T ) + b i, j+1/2 θ T with θ T = |φ i+1, j+1/2 | |φ i, j+1/2 | + |φ i+1, j+1/2 | . ( 83
)
Considering that the membrane crosses the bottom mesh segment, the same methodology is applied: if φ i+1/2, j φ i, j-1/2 < 0, the computation of the jump condition becomes

b B = b i+1, j-1/2 θ B + b i, j-1/2 (1 -θ B ) with θ B = |φ i, j-1/2 | |φ i, j-1/2 | + |φ i+1, j-1/2 | . ( 84
) Conversely, if φ i+1/2, j φ i, j-1/2 > 0, the jump condition becomes b T = b i+1, j-1/2 (1 -θ B ) + b i, j-1/2 θ B with θ B = |φ i+1, j-1/2 | |φ i, j-1/2 | + |φ i+1, j-1/2 | . ( 85
)
The jump conditions of each term of the viscous stress tensor are added to the right-hand side of the first step of the projection method without changing the coefficients of the coupled linear system. More details about the linear system to solve the velocity are presented in [START_REF] Lepilliez | On two-phase flow solvers in irregular domains with contact line[END_REF]. The extension of these schemes to 3D problems is straightforward.

Extension algorithm

The components of the solid tensors and the solid velocity u s are scalar fields defined in the whole computational domain. Nevertheless, they are only physically meaningful at the interface (where φ = 0). The values of these scalar fields in the fluids domain are not physical and must not influence the computation at the membrane position. Through equation [START_REF] Fedkiw | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF], shear flows or vortices, even far from the interface, may affect the solid tensors by generating parasitic strains. The tensors components may grow exponentially and lead to numerical instability [START_REF] Ii | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach[END_REF]. To avoid this problem, the authors replace the surface left Cauchy-Green tensor by the current projection tensor in the regions remote from the membrane. This method improves the computation near the membrane but stability problems still occur in the long run [START_REF] Ii | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach[END_REF].

In this paper, an extension algorithm has been developed to avoid these stability issues by extending the components of the velocity field and the solid tensor from the membrane toward its normal direction. Therefore, only the physical values are extended in the whole domain and the gradient of the scalar fields in the normal direction of the membrane equal zero. The methodology is inspired by the reinitialisation algorithm of level set functions [START_REF] Min | On reinitializing level set functions[END_REF][START_REF] Min | A second order accurate level set method on non-graded adaptive Cartesian grids[END_REF][START_REF] Russo | A remark on computing distance functions[END_REF] where a subcell resolution has been proposed. The extension of a scalar function h is described by the following equation

⎧ ⎨ ⎩ ∂h ∂τ + (sign(φ)n) • ∇h = 0, h(x, τ ) = h(x, 0) for τ > 0 and x ∈ , (86)
with τ a fictitious time. After the convergence, the normal derivative of the function tends toward zero in the whole domain and the scalar field h tends toward its value at the closer point of the membrane. The spatial derivatives are discretised using second order ENO finite differences in each direction. In the x-direction, two one-sided ENO finite differences exist,

∂h ∂x r = h i+1 -h i x - x 2 minmod h i+1 -2h i + h i-1 x 2 , h i+2 -2h i+1 + h i x 2 , ( 87
) ∂h ∂x l = h i -h i-1 x + x 2 minmod h i+1 -2h i + h i-1 x 2 , h i -2h i-1 + h i-2 x 2 , ( 88
)
with the minmod function defined as minmod(a, b) = min (|a|, |b|) if ab > 0, 0 else.

The choice of the one-sided finite difference depends on the sign of (sign(φ)n) • e x to propagate the function values from the membrane toward the fluid domain,

∂h ∂x = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ ∂h ∂x l if (sign(φ)n) • e x > 0, ∂h ∂x r else. ( 89
)
The method is identical in the other directions. The value of the function at the interface is of utmost importance because it will be extended to the whole domain. The interface does not conform with the mesh grid and special care is needed to compute the interface values within mesh cells. We propose in this paper a subcell method to improve the accuracy on the cells cut by the membrane.

First, considering that the membrane cuts the interval [x i , x i+1 ], we define θ q x the distance between x i and the membrane (cf. Fig. 2). This distance can be computed using the level set function and second order Taylor expansions

φ(x m ) = 0 = φ i + θ q x ∂φ ∂x i + (θ q x) 2 2 ∂ 2 φ ∂x 2 i + O ( x 2 ). ( 90
)
We construct the derivatives of φ using central finite differences and solve the quadratic polynomial to find θ q , .!_ ( (</,;+1 -<Pi-1) -sign(<f,;)J (</,;+1 -<Pi-1) 2 -8<f,;( </J ;+1 -2</,; + </J ;-1)

) if 1 8 2 t I; 1 > E, 0 _ 2 </,;+1 -2</,; +<Pi-1 ax q- 2<Pi else. <Pi+l -<Pi-1 (91) (92)
The value of the scalar field hm at the membrane location can then be computed. If the membrane cuts the interval [X;,Xi+tl, 1 1 hm= 2 (-0q + BJ)h;-1 + (1 -0J)h; + 2 (0q +0J)h;+1, and, if the membrane cuts the interval [x;_1 ,x;],

(94)

The choice of the data points allows us to not favour values of points from one region instead of the other. Even if there might have a ditference in the computation of h at the membrane, the computation is consistent because the error decreases with the mesh size. Finally, the gradient computation of the extension algorithm is modifi ed so that the point located at the membrane is considered in the numerical scheme. For example, if the membrane eut the interval [x;, X ; + 1 J, él h ct ( hi+1-2h;+hi-1 h;-2h;-1+hî-2

)

" t -0 + 2 m m m o 2 , 2 . oX q6.X 6.X 6.X (95) ( 96
)
This method is identical in the other directions. The example in section 4.2 justifies the use of the subcell method by comparing results from the extension algorithm with or without the subcell resolution.

To avoid the time step restriction due to the subcell method, the extension algorithm is solved implicitly with the Gauss-Seidel method. The extension algorithm can then be written as The fi ctitious time step 6. r is defined as 6.r =amin(6.x , 6.y , 6.z) , (97) (98) with a = 0.3 which minimises the computation time to reach the fully extended field. Few iterations of the Gauss-Seidel algorithm are needed to reach the extended state because the scalar field at the membrane faintly evolves at each time step.

In our fluid-membrane solver, the extension algorithm is applied at each time step on the solid velocity field u s which allows to compute the components of the Cauchy-Green tensor (see section 3.3). The solid velocity field corresponds to the fluid velocity at the membrane extended in the whole domain. This algorithm is also applied directly on the components of Ḡs every 10 temporal iteration. The execution of the extension algorithm on the tensor components at each time step would increase the computation time and decreasing the extension frequency would generate fictitious normal derivatives at the membrane location.

Time step constraints and stability condition

The effects of the convection and the membrane elasticity lead to time step restrictions to ensure the numerical stability:

t conv = 1 max (|u|) x + max (|v|) y + max (|w|) z and t elas = 1 2 max(ρ + , ρ -) E s min( x, y, z) 3/2 . ( 99
)
The restriction due to the elasticity is similar to the restriction due to the surface tension effects in two-phase flow cases [START_REF] Kang | A boundary condition capturing method for multiphase incompressible flow[END_REF]. There is no restriction due to viscosity because of the implicit computation of the viscous term. The global condition on the time step becomes

1 t > 1 t conv + 1 t elas . ( 100
)
All the numerical algorithms, from the transport of the level set function to the solving of the Navier-Stokes equations with the elastic contribution, are integrated in a second order TVD Runge-Kutta scheme.

To fully resolve the dynamics near the membrane, the mesh grid must be fine enough to compute accurately the elastic contribution on the flow field. This results in a stability condition depending on a mesh Reynolds number. This mesh Reynolds number is based upon the characteristic velocity due to the membrane elasticity defined as U elas = E s /(max(ρ + , ρ -)L) and on the mesh size x such as

Re x = U elas x min(ν + , ν -) = x min(ν + , ν -) E s max(ρ + , ρ -)L , ( 101
)
with ν = μ/ρ the kinematic viscosity. The stability condition requires a mesh Reynolds number around one or less.

Numerical results

In this section, we compute some well-known benchmarks and compare our numerical results with those from the literature. In a first time, we present a detailed comparison between the proposed numerical solver against a theoretical solution involving the droplet migration induced by a prescribed surface tension gradient. This preliminary test-case allows us demonstrating the correct behaviour of the proposed numerical discretization to impose the jump condition on the normal derivative of the tangential velocity. Next, the improved extension algorithm is implemented on the benchmark from [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF] to highlight the efficiency of the method and its convergence rate. Then, the stretched and pressurized membrane immersed in a fluid, first introduced in [START_REF] Tu | Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods[END_REF], is considered with these new numerical methods. After the grid sensitivity study, the results on the radii evolution, the pressure jump and the velocity fields are compared to the results from [START_REF] Tan | An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane[END_REF]. The test case is extended to different fluids in the computational domain. First, different viscosities, and then, different densities, are enforced on each side of the membrane. Finally, the capsule immersed in a shear flow from [START_REF] Pozrikidis | Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow[END_REF] is implemented and its deformation throughout time is compared to the results from the literature.

Test-case: rising of a droplet in a quiescent liquid due to Marangoni stresses

In this section, we provide a numerical example to test the accuracy of the jump condition computation on the normal derivatives of the tangential velocity, by performing a simulation of an axisymmetric droplet with surface tension gradient. This configuration is based on the Marangoni effects test-case presented in [START_REF] Muradoglu | A front-tracking method for computation of interfacial flows with soluble surfactants[END_REF]. Along a droplet interface, non uniform local temperature or surfactants concentration can trigger gradients of surface tension γ , which generates a tangential stress jump. Then, a flow, called Marangoni convection, occurs to balance the gradients of γ , leading to a rising motion of the droplet in the quiescent liquid. This viscous tangential stress jump is assimilable to the jump condition of equation ( 46) where We performed axisymmetric simulations in a domain of size [l r , l z ] in respectively radial and longitudinal directions, in a moving frame until a steady state. We consider a droplet of radius R, viscosity μ + and density ρ + , in a quiescent liquid (μ -, ρ -), with a gradient of surface tension imposed by

f τ = ∇ s γ . (102)
γ (z) = γ 0 1 -β z l z , ( 103
)
where γ 0 = 10, β = 0.27, l z = 16R, and l r = l z /2. To ensure that the motion is solely due to Marangoni stresses, we impose μ + = μ -= 0.1, and ρ + = ρ -= 0.02. From the theoretical studies of [START_REF] Hadamard | Motion of liquid drops (viscous)[END_REF] and [START_REF] Rybczynski | Uber die fortschreitende Bewegung einer flussigen Kugel in einem zahen Medium[END_REF], where the Stokes equations are solved around and within a spherical droplet, an analytical expression of the droplet rising velocity is proposed in [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF] accounting for the jump condition due to Marangoni stresses. In [START_REF] Muradoglu | A front-tracking method for computation of interfacial flows with soluble surfactants[END_REF], the authors have adapted this analytical velocity to propose a benchmark case without any gravity:

u ∞ = 2γ 0 β R l z (6μ -+ 9μ + ) . ( 104
)
This expression is valid for a spherical droplet, i.e. at Capillary number Ca = μ -u d /γ 0 1, and for the creeping flow limit at Reynolds number Re = ρ -u d D/μ -1, where u d is the droplet velocity and D its diameter. The parameters were chosen in a such a way that Re = 0.045 and Ca = 0.0023, based on the terminal velocity of equation ( 104). From the studies of [START_REF] Hadamard | Motion of liquid drops (viscous)[END_REF], [START_REF] Rybczynski | Uber die fortschreitende Bewegung einer flussigen Kugel in einem zahen Medium[END_REF] and [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF], and considering the tangential stress jump as a boundary condition, we have calculated the analytical expression of tangential velocity around and within the droplet (more details are provided in the Appendix C),

u - θ = -u ∞ 1 + 1 2 R 3 r -3 sin(θ ) (105) u + θ = - 3u ∞ 2R 2 2r 2 -R 2 sin(θ ) ( 106
)
with u - θ and u + θ the tangential velocities respectively outside and inside the droplet. From this solution, the tangential velocity at the interface is given by the following expression,

u - θ | r=R = u + θ | r=R = 3 2 u ∞ sin(θ ) . ( 107
)
In Fig. 3, are plotted the droplet shape and the pressure field inside the droplet which generates its migration. In Figs. 4 and5, the dimensionless rising velocity of the bubble is plotted until steady state for four meshes with respectively implicit 1. 2�--------------- 1. 2�--------------- and explicit resolutions for the temporal discretization of the viscous terms. The error rates were calculated by comparison with the analytical solution. The numerical results are in good agreement with the theoretical expression of the rising velocity provided by eq. ( 104), for both resolutions.

ln Figs. 6 and7, we plot the dimensionless tangential velocity profile along the droplet surface, compared to the analytical expression eq. ( 107) that we present in the appendix. The more the mesh is refined, the better the agreement is with the theoretical solution.

ln Figs. 8 and9, the tangential stress profile µ,ôu 8 /ôn is plotted along the axis z = 0 for the most refined mesh, at the droplet equator, and compared to the tangential stress calculated from the solutions eq. ( 105) and eq. ( 106). Then, the numerical accuracy of this jump condition is verified in Table 1 by comparing the tangential stress jump from simulations to the calculated v'sY value from eq. (103). At each mesh refinement, the error rate is divided by two, typical of a first order accuracy as it was expected for the numerical discretization of the singular source terms.

F.xtension algorithm

ln this section, we present numerical results of the extension algorithm coupled with the subcell resolution on an example inspired by [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF]. Let us consider a 2D computational domain [-rr, rr] x [-rr, rr]. Initially, a level set function is defined as <f>(x, y) = Jx 2 + y 2 -Ro, the zero level set of which corresponds to a circle of radius Ro = 2. A scalar field denoted by h is initialised as h(X, y) = COS(X) sin(y), (108) -------------------� � t and some contours are shown in the left plot of Fig. 10 with the zero level set line depicted in solid black.

The extension algorithm described in section 3.6 is enforced on the scalar field h and the result is plotted on the right plot of Fig. 10. We observe that the contours of the obtained scalar field respect well the condition of no normal derivative at the membrane. The exact solution of the extended field h is easily computable and equals Fig. 11 focuses on the scalar field h near the membrane, represented by a bold solid black line, in the specific region [O, 1.1) x [1.3, 2.4). The initial contours of h are plotted with solid black lines and the extended scalar field with dashed black lines. ln the left plot of Fig. 11, we obseive that the values of the scalar field h at the membrane do not change after the extension. The solid lines correspond exactly to the dashed lines at the membrane location. The same computation is done without the subcell method and the result is illustrated in the right plot of Fig. 11. We obseive that the contours do not correspond as well. The zero normal derivative at the membrane is well respected but the values of the scalar field at the membrane have been modified. The scalar field variation is more important where the initial scalar field is tangent to the membrane. lndeed, the contour h (x, y) = 0.9 after the extension without the subcell method exhibits an important variation at the membrane position.

h(x, y) = cos R 0 x x 2 + y 2 sin R 0 y x 2 + y 2 . ( 109
)
Table 2 summarises the errors of the extended scalar field with the subcell resolution compared to the exact solution of equation (109). Iwo errors are computed: in the whole domain, except from the centre of the circle where the normal directions to the membrane meet up, and in the region near the membrane. The extension algorithm with the subcell method is third order accurate near the zero level set and almost second order in the whole region.

Considering the extension algorithm without the subcell method, the accuracy decreases notably, as shown in Table 3. Near the membrane, the order of magnitude of the error is at least two times higher than with the subcell resolution. The convergence rate of the method is around one in the whole domain and near the zero level set. The study of the accuracy and the contour plots of the scalar field h explains that the subcell method must be considered to extend the right value of the scalar field at the membrane toward its normal direction. Now, the level set function is modified to obtain an ellipse of major axis a = 2.5 and minor axis b = 1.25 as the zero level set. The initialisation of the scalar field h(x, y) stays identical and the exact solution of the extended field is obtained

Table 4

Numerical accuracy of the extension method with the subcell resolution in the case of the ellipse.

Grid

Whole domain |y| > 0. [START_REF] Ii | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach[END_REF] Near the membrane |φ(x,

y)| < 1.2 x E 1 R 1 E ∞ R ∞ E 1 R 1 E ∞ R ∞ 64

2.95

Table 5

Numerical accuracy of the extension method without the subcell resolution in the case of the ellipse.

Grid

Whole domain |y| > 0. [START_REF] Ii | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach[END_REF] Near the interface |φ(x, by an iterative process. Fig. 12 shows the initial and final scalar field h(x, y) with the ellipse corresponding to the zero level set represented in bold solid black. The same study on the extension algorithm is done on the ellipse case with and without the subcell resolution and are detailed in Tables 4 and5 respectively.

y)| < 1.2 x E 1 R 1 E ∞ R ∞ E 1 R 1 E ∞ R ∞ 64
The errors with the exact solution are computed in the whole domain, except from the region |y| < 0.15 where singularities due to the computation of the normal direction to the membrane appear, and in the region near the membrane. With the subcell resolution, the errors in the whole domain and near the membrane are always at least one order of magnitude lower than without it. The extension algorithm is still third order accurate near the membrane with the subcell resolution and second order accurate in the whole domain. Without the subcell resolution, the computed order of convergence is important, around 2, for the coarser meshes but decreases significantly with finer meshes.

The study on the ellipse case shows that the third order accuracy obtained near the interface is due to the subcell resolution and is not an artefact due to the simple circular membrane. Without the subcell resolution, the computed order of convergence oscillates between 1 and 2 depending on the membrane shape and the mesh considered.

Stretched and pressurized membrane immersed in a fluid

The stretched and pressurized membrane immersed in a viscous fluid is implemented in this section. It has been used by several authors to test the immersed boundary method and the immersed interface method [START_REF] Lee | An immersed interface method for incompressible Navier-Stokes equations[END_REF][START_REF] Leveque | Immersed interface methods for stokes flow with elastic boundaries or surface tension[END_REF][START_REF] Tu | Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods[END_REF]. In this example, the elastic contribution will be implemented using the Ghost-Fluid viscous Conservative Method of section 3.2.2. It consists in a 2D membrane, the shape of which is an ellipse of major axis a = 0.75 and minor axis b = 0.5, relaxing in a fluid at rest. In its stress-free state, we consider that the membrane is a circle of radius r 0 = 0.5. The equilibrium shape of the initially stretched membrane is a circle with the same amount of liquid within. The equilibrium radius equals r e = √ ab ≈ 0.61237. The domain is a closed square of length L = 3 and wall boundary conditions are considered. The ellipse is initially centred in the domain. The three different states of the membrane are depicted in Fig. 13. Different densities and viscosities of the fluid are considered: ρ ∈ [1, 10, 100] and μ ∈ [0.01, 0.1, 1]. Depending on the fluid properties, the membrane will relax toward its circular shape differently.

In this benchmark, the membrane follows the linear Hooke law. However, the numerical methods presented previously in this paper deal with hyperelastic material laws. Fortunately, in 2D cases, the linear Hooke law can be considered using the following methodology because the membrane is only stretched along its tangent direction denoted e τ . The first invariant of the Cauchy-Green tensor gives directly the square of the principal strain λ τ . The linear Hooke law corresponds to the following surface stress tensor σ s = E s ( I 1 -1) P , (110) with I 1 = λ 2 τ , P = e τ ⊗ e τ and the value of the surface elastic modulus E s = 10.

Therefore, the surface stress tensor can directly be computed with the projection tensor and the first invariant of the surface Cauchy-Green tensor. The latter tensor is initialised by the projection tensor multiplied by the initial stretching of the membrane. Considering a uniform stretching, it corresponds to the square of the ratio of the membrane perimeter between its initial state P 0 and its stress-free state P R ,

Bs (t = 0) = λ 2 τ (t = 0) P (t = 0) = P 0 P R 2 P (t = 0) (111)
with P 0 /P R ≈ 1.262.

The evolution of the radii is investigated and compared to the results from [START_REF] Tan | An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane[END_REF]. First, to ensure the validity of our computation, a grid sensitivity study has been performed for ρ = 1 and two different viscosities μ = 0.1 and μ = 0.01. The computations have been performed on three different meshes for each viscosity. Fig. 14 shows that for μ = 0.1, the three meshes, 128 2 , 256 2 and 512 2 , give the same radii evolution. In this case, the 128 2 mesh is fine enough to do the comparison. With μ = 0.01, finer grids are needed to reach a steady radii evolution. The temporal evolutions are close to each other at early flow time but deflect with time. This deviation decreases when refining the mesh. For meshes containing more than 512 cells in each direction, the same solution is obtained and gives consistent results with [START_REF] Tan | An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane[END_REF]. In what follows, the 512 2 mesh will be considered in case of a viscosity of μ = 0.01. For lower values of the viscosity, the 128 2 mesh is considered.

We observe in Fig. 14 that in both cases, the membrane oscillates around its equilibrium value. The magnitude of the oscillations decreases with time until the membrane reaches the circular shape. The theoretical radius r e is obtained for the horizontal r x and vertical r y radii in the steady state. The damping of the oscillations increases with the viscosity but its frequency remains stable. Fig. 15 shows the pressure distribution in the computational domain at t = 0.5 and t = 2. The sharp jump on the pressure is clearly visible at the membrane location. At early flow time, the pressure field exhibits important variations in each fluid .._:.. region because of the fluid motion. At t = 2, the pressure field is almost constant inside and outside the membrane. In the steady state, only the pressure jump remains with constant values of the pressure in each phase.

Fig. 16 shows the velocity field at two different times with μ = 0.1 in the whole domain. At t = 0.5 and t = 2, four main vortices are located in the fluid outside the membrane leading the membrane to deform in the vertical direction. The pressure and velocity fields of Figs. 15 and 16 are in good agreement with the results from [START_REF] Tan | An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane[END_REF]. These simulations have also been implemented without the extension algorithm and the results are presented in Appendix A. This study shows that the extension algorithm is necessary to keep the computation stable throughout time by removing the normal derivatives of the solid scalar fields.

The membrane behaviour is investigated for higher values of the density. The viscosity is kept at μ = 0.1 and the density is increased to ρ = 10 and ρ = 100. A grid sensitivity study is done on several meshes. We observe in Fig. 17 that the radii evolution reaches the same solution from 128 cells onwards in each direction with ρ = 10 and from 512 cells onwards in each direction with ρ = 100. In what follows, the 512 2 mesh will be considered when the density reaches ρ = 100, the 128 2 mesh will be used otherwise.

We observe that the time to relax the membrane increases with the density. The frequency of the oscillations seems to decrease with the square root of the density. Fig. 18 exhibits two examples of the pressure field for the two higher densities. The sharp jump condition on the pressure is well predicted by the numerical simulations even in the right plot of Fig. 18 where the pressure variation is important in each region.

Stretched and pressurized membrane immersed in fluids with different viscosities

Now, we distinguish the fluids inside and outside the membrane. In this section, the two fluids have different viscosities μ + and μ -. As depicted in Fig. 13, the superscript + corresponds to the fluid outside the membrane and the superscript - to the fluid inside the membrane. The same density ρ + = ρ -= 1 is enforced and following the work of [START_REF] Tan | An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane[END_REF], the viscosity ratio takes the values of 0.1 and 10.

First, we consider two cases with large viscosities: μ + = 1, μ -= 0.1 and μ + = 0.1, μ -= 1. Fig. 19 shows the radii evolution throughout time in both cases. The membrane does not oscillate in the first case, it relaxes gradually to its final position. When the fluid with the higher viscosity is inside the membrane, the latter performs only one small oscillation. The vertical and horizontal radii tend toward the theoretical value r e in both cases.

Then, the same viscosity ratios are maintained but the fluid viscosities are divided by 10. Fig. 20 shows the radii evolution for these two cases. The vertical and horizontal radii oscillate around the equilibrium value with frequencies of the same order of magnitude. Similarly to Fig. 19, the damping is more important when the fluid with the higher viscosity is outside the membrane. The two plots of Fig. 19 and 20 are really close to the results from [START_REF] Tan | An immersed interface method for solving incompressible viscous flows with piecewise constant viscosity across a moving elastic membrane[END_REF].

An example of the jump on the velocity derivatives is shown in Fig. 21. The evolution of both components of the velocity at y = 0.3 is plotted at the same time t = 1 with μ + = 0.1 and μ -= 1. We observe that in both cases, the velocity is continuous but its derivatives are discontinuous at the membrane location (represented by the vertical dashed lines). This validates the sharp methodology of the jump condition on the viscosity-scaled velocity gradient.

Stretched and pressurized membrane immersed in fluids with different densities

In this section, we consider that the fluids inside and outside the membrane have the same viscosity μ + = μ -= 0.1 but different densities. First, we enforce densities of 1 and 10 on both sides of the membrane. The radii evolutions of these two cases are represented in Fig. 22. As observed in the left plot of Fig. 17, the relaxing time of the membrane increases with the density. The membrane oscillates around its circular shape, the radius of which corresponds to the theoretical value r e at the end of the simulation. It takes longer to reach the equilibrium state when the heaviest fluid is inside the membrane. The frequency of the oscillations has the same order of magnitude in both cases. Then, the fluids densities are multiplied by 10. Fig. 23 exhibits the radii evolutions throughout time. We observe that the radii evolution with ρ + = 10 and ρ -= 100 is close to the radii evolution with ρ = 100 in the whole domain (cf. Fig. 17). Similarly, the damping rate of the oscillations is less efficient when the heaviest fluid is inside the membrane. The oscillation frequency has the same order of magnitude in these two cases.

Because the viscosity stays constant in the whole domain, the jumps on the velocity derivatives are directly linked to the tangent component of the elastic force. Fig. 24 shows the profiles of the two components of the velocity at y = 0.3 with ρ + = 10 and ρ -= 1. The jump on the velocity derivatives is clearly visible at the membrane location.

Immersed capsule in a shear flow

In this section, we consider a 3D spherical membrane immersed in a shear flow. The membrane deforms until the equilibrium between the viscous and the elastic forces is reached. This study case has been first introduced by Pozrikidis [START_REF] Pozrikidis | Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow[END_REF] and has been used to validate several immersed boundary methods [START_REF] Eggleton | Large deformation of red blood cell ghost in a simple shear flow[END_REF][START_REF] Ramanujan | Deformation of liquid capsule enclosed by elastic membranes in simple shear flow: large deformation and the effect of fluid viscosities[END_REF]. This example aims at validating the use of hyperelastic laws with our model.

Initially, the membrane is a sphere of radius a. It is located in the centre of the computational domain, a cube with a side of 4a. The shear flow enforced in the domain is characterised by the shear rate k such as u = (ky, 0, 0) at the top and bottom sides. On the other sides of the domain, periodic boundary conditions are imposed. We consider that the same fluid is inside and outside the membrane. Its physical properties are its density ρ and viscosity μ.

D xy = L -l L + l , ( 113
)
with L and l the maximum and minimum radii of the contour of the membrane with the (x, y)-plane computed with the non-linear least square method.

In this example, we use the "delta" formulation described in section 3.2.1. Sharp methods cannot be used in this example because of the compression zones on the membrane [START_REF] Barthes-Biesel | Capsule motion in flow: deformation and membrane buckling[END_REF][START_REF] Lac | Spherical capsules in three-dimensional unbounded stokes flows: effect of the membrane constitutive law and onset of buckling[END_REF][START_REF] Ramanujan | Deformation of liquid capsule enclosed by elastic membranes in simple shear flow: large deformation and the effect of fluid viscosities[END_REF]. Indeed, where at least one of the two principal membrane stresses becomes negative, a buckling instability appears [START_REF] Lac | Spherical capsules in three-dimensional unbounded stokes flows: effect of the membrane constitutive law and onset of buckling[END_REF]. Because no bending stiffness is considered, a zone of the membrane which is not under tension may buckle and break down the computation. In this case, folds appear at the membrane surface which may overlap and lead to a stability issue. Fig. 25 shows an example of folds which appear in the equatorial area in our computation with G = 0.0125 after kt = 0.8. Sharp methods, based on the local equilibrium at the nodes crossed by the membrane, become unstable when compression zones appear [START_REF] Barthes-Biesel | Capsule motion in flow: deformation and membrane buckling[END_REF]. The use of the smoothed method introduces a numerical dissipation which stabilises the computation and allows the solver to tolerate negative stress locally.

Nevertheless, in specific cases such as G = 0.0125, stability issues appear in the long run whatever the numerical methods used.

Fig. 26 shows the evolution of the deformation parameter D xy throughout time for 5 different values of G and three different meshes. We observe that the membrane deformation increases with G. At early flow time, the membrane loses its spherical shape and deforms with different magnitudes until it reaches a steady shape. The computations have been done on a maximum time of kt = 2 to be able to compare with the results from the literature. For the two lowest G, the buckling instability described below appears before kt = 2 and prevents us to compute the deformation parameter in the whole time period.

With G = 0.0125 and G = 0.025, the membrane deformation stays low and can be approached by the linear elasticity theory for which theoretical values of the membrane deformation have been computed in [START_REF] Barthes-Biesel | The time-dependent deformation of a capsule freely suspended in a linear shear flow[END_REF]. The evolution of the Taylor parameter with the linear theory is depicted by triangle markers in Fig. 26. The Taylor deformation parameter obtained Considering higher values of G, the linear theory cannot be used because the non-linearity of the material cannot be neglected. Our numerical results are compared to simulations done with the Immersed Boundary method and the Boundary Element method [START_REF] Eggleton | Large deformation of red blood cell ghost in a simple shear flow[END_REF][START_REF] Pozrikidis | Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow[END_REF][START_REF] Ramanujan | Deformation of liquid capsule enclosed by elastic membranes in simple shear flow: large deformation and the effect of fluid viscosities[END_REF]. The Taylor parameter obtained by Pozrikidis [START_REF] Pozrikidis | Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow[END_REF] for the highest values of G are plotted in Fig. 26 with filled triangles. Pictures of the steady state of the membrane for G ∈ [0.0125, 0.05, 0.2] are depicted in Fig. 27 and show different deformations computed with our numerical methods.

Fig. 1 .

 1 Fig. 1. Location of the variables on a 2D mesh cell.

Fig. 2 . 8 ,

 28 Fig. 2. The membrane r crossing the x-grid.

  the membrane eut the interval [x;_ 1 ,x;],

Fig. 3 .

 3 Fig. 3. Pressure field and streamlines of the rising droplet. (For interpretation of the colours in the figures, the reader is referred to the web version of this article.)

Fig. 4 .

 4 Fig. 4. Rising velocity of the droplet as a function of rime, for the implicit resolution, with the error rate when comparing to the theoretical solution.

  Fig. S. Rising velociry of the droplet as a function of rime, for the explicit resolution, with the error rate when comparing to the theoretical solution.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Tangenàal velocicy profile a long the drop let surface, as a funccion of the angle 0. for the implicit resolucion.

Fig. 8 .

 8 Fig. 8. Tangencial viscous stress profile along the radius, ac z = 0 for the implicit resolution.

Fig. 9 .

 9 Fig. 9. Tangential viscous stress profile along the radius, at z = 0 for the explicit resolution.

Fig. 10 .

 10 Fig. 10. Initial and final scalar fields h(x, y) with the circular zero level set represented by the bold solid black line.

6 XFig. 11 .

 611 Fig. 11. Zoom on the scalar field h(x, y) before and afrer the extension with the subcell resolution on the left and without on the right.

Fig. 12 .

 12 Fig. 12. Initial and final scalar fields h(x, y) with the zero level set represented by the bold solid black line.

Fig. 13 .

 13 Fig. 13. The membrane in different states.

Fig. 14 .Fig. 15 . 2 • 1 : 1 : 5 Fig. 16 .

 1415211516 Fig. 14. Radii evolution with µ, = 0.1 on the left and µ, = 0.01 on the right for three different meshes.

Fig. 17 .

 17 Fig. 17. Radii evolution with ρ = 10 on the left and ρ = 100 on the right for three different meshes.

Fig. 18 .

 18 Fig. 18. Pressure distribution with ρ = 10 at t = 1.5 on the left and with ρ = 100 at t = 5 on the right.

Fig. 19 .

 19 Fig. 19. Radii evolution with μ + = 1 and μ -= 0.1 on the left and μ + = 0.1 and μ -= 1 on the right.

Fig. 20 .

 20 Fig. 20. Radii evolution with μ + = 0.1 and μ -= 0.01 on the left and μ + = 0.01 and μ -= 0.1 on the right.

Fig. 21 .

 21 Fig. 21. Components of the velocity u = (u, v) at t = 1 and y = 0.3 with μ + = 0.1 and μ -= 1.

Fig. 22 .

 22 Fig. 22. Radii evolution with ρ + = 10 and ρ -= 1 on the left and ρ + = 1 and ρ -= 10 on the right.

Fig. 23 .

 23 Fig. 23. Radii evolution with ρ + = 100 and ρ -= 10 on the left and ρ + = 10 and ρ -= 100 on the right.

Fig. 24 .

 24 Fig. 24. Components of the velocity at y = 0.3 with ρ + = 10 and ρ -= 1 at t = 1.5.

  The membrane is made of a hyperelastic material following the neo-Hookean model. The surface elastic modulus of the membrane is denoted E s . The membrane is initially at rest, leading to Bs = Ḡs = P at t = 0. We define the dimensionless parameter G which compares the effect of the viscosity of the fluid over the elasticity of the membrane of G are considered: G ∈ [0.0125, 0.025, 0.05, 0.1, 0.2]. The deformation of the membrane is expressed through the Taylor deformation parameter

Fig. 25 .

 25 Fig. 25. Folds due to the buckling instability at the membrane surface with G = 0.0125 at kt > 0.8.

Fig. 26 .

 26 Fig. 26. Temporal evolution of the Taylor deformation parameter for the five values of G and three different meshes: 32 3 , 64 3 and 128 3 .

Fig. 27 .

 27 Fig. 27. Steady state or final state of the membrane for a) G = 0.0125, b) G = 0.05 and c) G = 0.2.

Table 1

 1 Numerical accuracy of tangential stress jump compared to the calculated value

	Grid	Implicit	Explicit
	64 × 128	15.98%	17.21%
	128 × 256	8.68%	9.22%
	256 × 512	4.49%	4.68%
	512 × 1024	2.28%	2.37%

of ∇ s γ from eq. (

103

).

1.82 1.76 X 10-2 1.25 6.56 X 10-6 3.03 4.37 X 10-S 2.87 256 2 9.01 X 10-S 1.91 5.71 X 10-J 1.62 8.24 X 10-7 2.99 6.20x 10-6 2.82 512 2 2.30 X 10-S 1.97 1.65 X 10-J 1.79 1.01 X 10-7 3.03 7.27x 10-7 3.09 Table3 Numerical accuracy of the extension method without the subcell resolution in the circular case. Grid Whole domain <f

	f1	R1	foo	Roo	f1	R1	foo	Roo

,(x, y)> -1.6 Near the interface l,f,(x, y)I < 1.2.1.x

642 1.11 X 10-2 6.12 X 10-2 1.06 X 10-2 2.82 X 10-2 128 2 3.47 X 10-J 1.68 2.38 X 10-2 1.36 3.50 X 10-J 1.59 8.72 X 10-J 1.69 256 2 1.66 X 10-J 1.06 9.30 X 10-J 1.35 1.68 X 10-J 1.06 5.77 X 10-J 0.60 512 2 6.93 X 10 -4 1.26 3.68 X 10-J 1.34 6.91 X 10 -4 1.28 2.74 X 10-J 1.07

We observe in Fig. 26 the good convergence of the proposed numerical methods since the solution depends weakly on the grid size. In particular, the Taylor deformation parameter is the same for both finest grids for all values of G. Our numerical results are in close agreement with the reference results presented in Fig. 26. However, we report a slight overestimation on the Taylor deformation parameter whatever the considered case.

Acknowledgements

The authors thank Airbus Defence & Space and the CNES (the French national space agency) for the funding of the PhD thesis of Alexis Dalmon. The authors gratefully acknowledge the CNRS (Centre National pour la Recherche Scientifique) for funding compensation of teaching hours to S. Tanguy during the year 2017-2018. This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation P17035.

Conclusion

ln this paper, we propose a full Eulerian method to deal with fluids-membrane interaction inspired by the work of lto et al. in [START_REF] Ii | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach[END_REF]. The main advantage of the Eulerian approach is to replace the use of l.agrangian markers to follow the membrane motion and deformations by scalar fields transported in the whole domain. Significant improvements are proposed in this paper to extend the method firstly introduced in [START_REF] Ii | A full Eulerian fluid-membrane coupling method with a smoothed volume-of-fluid approach[END_REF]. A sharp methodology has been implemented to consider the elastic forces as jump conditions in the fluid equations. Specific developments have been done to be able to predict the behaviour of a membrane separating different fluids in the Eulerian framework. Because of important stability issues, an extension algorithm has been developed to remove the parasitic normal derivative of the scalar fields specific to the membrane. These achievements have been tested on benchmarks from the literature and ensure the suitable prediction of the proposed numerical methods. The jump conditions on the pressure and the velocity derivatives are well predicted, even in cases where the viscosity or the density are piecewise constant across the membrane.

Declaration of competing interest

The authors declare that they have no known competing financial interests or persona! relationships that could have appeared to influence the work reported in this paper.

AppendixA. Stretched and pressurized membrane without the extension algorithm

The ob j ective of this appendix is to show through an example that the extension algorithm is necessary. We consider the benchmark of section 4.3, the stretched and pressurized membrane immersed in a fluid with p = 1 and µ, = 0.1 in the whole computational domain. The computation breaks down after a few hundreds of temporal iterations without the extension algorithm, whatever the mesh. The comparison is done at early flow time, at t = 0.1. We observe in the left plot that the first scalar invariant of the strain tensor exhibits important variations at the mem brane. Where the membrane crosses the axis y = 0, the first invariant increases sharply between the outside and inside fluid regions. Moreover, strong tangential variations of the scalar invariant are visible at the membrane. ln the right plot, the scalar field does not present variations in the normal direction and the tangential variation along the membrane is smooth and does not present singularities. Without the extension algorithm, nothing prevents the appearance of jumps in the solid variables. -u=0.15 1.5 ------------- This instability phenomenon appears for ail the benchmarks tested in section 4 without the extension algorithm. To conclude, the extension algorithm stabilises the computation by cancelling the normal derivatives of the solid variables which enables improving the computations stability.

Appendix B. About primary and secondary jump conditions

We present in this section some developments to show how secondary jump conditions can be taken into account by imposing primary jump conditions. Indeed, it is well-known that primary jump conditions, especially the one on viscosity, will involve several secondary jump conditions on the pressure field or on the pressure gradient for instance. However, if one solves the Navier-Stokes equations for two-phase flows, these secondary jump conditions can be accounted for when imposing the primary jump conditions. For the sake of simplicity we will consider here a Stokes flow without elastic con straints and an explicit solver for viscous terms, but such a development can be generalized to more complex configurations. From [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF]

where the following term has been decomposed in a continuous part and a discontinuous part

Finally, following previous developments presented in [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF], the equation (B.6) enables showing that equation (B.1) for the pressure field contains both a numerical approximation for imposing the correct jump condition on the pressure,

and also for the jump condition on the pressure normal gradient such as

As a result of these developments, one can understand that secondary jump conditions due to the viscosity jump are involved by the resolution presented here. In particular, the jump condition on the normal derivative of the normal component velocity is imposed by keeping the viscosity inside the divergence operator in (B.1) as previously stated in [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF]. As a consequence, this term has not to be further imposed in [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF], since it is already contained in the divergence of u * . These elements are supported by theoretical and experimental benchmarks [START_REF] Lalanne | Non-linear shape oscillations of rising drops and bubbles: experiments and simulations[END_REF][START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF] for which space convergence toward the correct solution has been showed in several configurations involving a significant effect of the viscosity jump.

Appendix C. Capillary rising of a drop

The performed simulations are based on the theoretical study of [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF] and the test-case in [START_REF] Muradoglu | A front-tracking method for computation of interfacial flows with soluble surfactants[END_REF]. In the latter, we consider a drop of radius R in micro-gravity conditions (g = 0) immersed in another fluid where a linear profile of surface tension γ is imposed along z direction:

The constant surface tension gradient induces a movement of the drop at a steady-state velocity u ∞ . Both fluids are supposed to be incompressible, and their thermo-physics properties are chosen such that Capillary, Marangoni and Reynolds numbers are small enough to consider a Stokes flow around the droplet which remains spherical. Thus, the flow remains axisymmetrical, and the conservation equation can be written in the moving referential as:

with the asymptotic boundary conditions u → u ∞ e z and p → 0 as r → ∞. The solution of those equations has been derived by Hadamard in [START_REF] Hadamard | Motion of liquid drops (viscous)[END_REF], in the cylindrical coordinates:

where u - r and u + r are respectively the radial velocities outside and inside the droplet, u - θ and u + θ the tangential velocities, pand p + the pressures, μ -and μ + the dynamic viscosities, a -, a + and a + 0 three integration constants. These ones are calculated with the boundary conditions at the interface:

(C.12)

The no-slip condition gives the first equation. The second one takes into account the Marangoni effect through the discontinuity of the tangential viscous stress across the interface. The latter gives the pressure jump due to both viscosity jump and capillary pressure. By replacing the theoretical expressions in eq. (C.12) and eq. (C.11), and combining the equations, we deduce the integration constants: (C. [START_REF] Sugiyama | A full Eulerian finite difference approach for solving fluid-structure coupling problems[END_REF] This expression is actually the same as in [START_REF] Muradoglu | A front-tracking method for computation of interfacial flows with soluble surfactants[END_REF] with γ = -βγ 0 /l z . Finally, we obtain the following analytical expressions: