Florian Desmons

Thomas Milcent

Anne-Virginie Salsac

Mirco Ciallella
email: mirco.ciallella@ensam.eu

Fully Eulerian models for the numerical simulation of capsules with an elastic bulk nucleus

Keywords: Biological capsules, Fluid-structure interaction, Eulerian elasticity, full membrane elasticity, Incompressible Navier-Stokes

Introduction

Biological capsules have revolutionized the fields of medicine and biotechnology, offering remarkable advancements and opportunities for research, development, and practical applications. For example, in medicine, the use of biological capsules have greatly improved drug delivery systems [START_REF] Del Mercato | Biological applications of lbl multilayer capsules: From drug delivery to sensing[END_REF]. By encapsulating drugs within protective coatings, scientists can enhance their stability, solubility, and targeted release. Another field of application in biotechnology concerns the encapsulation for the efficient confinement of enzymes within protective matrices, enhancing their reusability, and enabling their application in diverse industrial processes [START_REF] Li | Fabricating covalent organic framework capsules with commodious microenvironment for enzymes[END_REF].

In the context of studying the interaction between such capsules and a surrounding fluid, numerical simulation of fluid-structure interaction (FSI) is a very promising tool to understand their behavior and deformation. FSI simulations involve the coupled analysis of fluid flow and structural deformation, allowing researchers to investigate the complex interactions between structures and their surrounding fluid environments. When applied to biological capsules, FSI simulations provide valuable insights into their mechanical behavior, fluid dynamics, and overall performance. Nowadays, several methods are used to tackle the coupled problem coming from the interaction between fluid and structure. In the seminal work [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF], Peskin introduced an innovative approach where Lagrangian markers are used to track the membrane deformation. The elastic force appears as a source term in the fluid equations and is spread on the fluid grid with discretized Dirac mass. This work was established as a reference for future developments of the so-called immersed boundary method [START_REF] Mittal | Immersed boundary methods[END_REF][START_REF] Lee | An immersed interface method for incompressible navier-stokes equations[END_REF][START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Griffith | Immersed methods for fluid-structure interaction[END_REF][START_REF] Sotiropoulos | Immersed boundary methods for simulating fluid-structure interaction[END_REF]. Other popular approaches to solve such problems are based on arbitrary Lagrangian-Eulerian (ALE) methods [START_REF] Liu | A numerical study of insect flight[END_REF][START_REF] Sahin | An arbitrary lagrangian-eulerian formulation for the numerical simulation of flow patterns generated by the hydromedusa aequorea victoria[END_REF][START_REF] Donea | An arbitrary lagrangian-eulerian finite element method for transient dynamic fluid-structure interactions[END_REF][START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF][START_REF] Fernández | A newton method using exact jacobians for solving fluid-structure coupling[END_REF], which consider body-fitted grids to follow the displacement of the interface. Although ALE method has been widely used in many FSI problems, they become extremely cumbersome to apply to problems with large deformations, due to the challenging re-meshing process.

More recently, Eulerian models for fluid-structure interaction have been derived for computational purposes as they enable complex multi-dimensional applications, with large deformation, to be simulated quite easily on fixed Cartesian meshes [START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF][START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF][START_REF] Cottet | Méthodes Level Set pour l'interaction fluidestructure[END_REF]. In this framework, both the fluid and the solid are described by an Eulerian approach. The interface and deformations of the media are captured by Eulerian fields, which are advected by the fluid velocity in a level-set fashion [START_REF] Osher | Level set methods and dynamic implicit surfaces[END_REF]. The fluid-structure problem is recast as a complex flow: the Navier-Stokes equations with an elastic source term is coupled with a transport equation on the Eulerian interface and deformation. In the last decade, this approach has been applied to both incompressible and compressible materials [START_REF] Maitre | Applications of level set methods in computational biophysics[END_REF][START_REF] Bergmann | An eulerian finite-volume approach of fluidstructure interaction problems on quadtree meshes[END_REF][START_REF] Gorsse | A simple cartesian scheme for compressible multimaterials[END_REF][START_REF] De Brauer | A cartesian scheme for compressible multimaterial hyperelastic models with plasticity[END_REF] proving its high flexibility, especially when dealing with complex geometries.

When dealing with microcapsules, Stokes flows are sometimes considered due to the low Reynolds numbers, which corresponds to slow and highly viscous flows. Furthermore, the Stokes equations are simpler than the Navier-Stokes equations and can be solved more efficiently using the boundary integral method [START_REF] Ramanujan | Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of fluid viscosities[END_REF][START_REF] Pozrikidis | Interfacial dynamics for stokes flow[END_REF][START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a stokes flow[END_REF]. The boundary integral method directly solves the fluid flow problem on the surface of the capsule, eliminating the need for a volumetric discretization of the entire fluid domain. This reduction in dimensionality significantly reduces the computational cost. Difficulties remain for the simulation of large deformations, interactions with complex geometries and applications where the viscosity of the inner fluid is different from that of the outer fluid [START_REF] Foessel | Influence of internal viscosity on the large deformation and buckling of a spherical capsule in a simple shear flow[END_REF].

The novelty of this paper is to present a framework based on Eulerian models for FSI for the realistic simulation of capsules undergoing large deformations. In particular, the objective is to study the influence of an internal nucleus, modeled as an elastic solid [START_REF] Deborde | Numerical simulations of the interaction of solitary waves and elastic structures with a fully eulerian method[END_REF], within the capsule for long-time simulations. The introduction of an internal nucleus is extremely promising for numerous applications since many biological cells, e.g. white blood cells, are often modeled as a thin membrane with a nucleus, also modeled as a thin elastic membrane [START_REF] Banaei | Numerical simulations of elastic capsules with nucleus in shear flow[END_REF]. Once the proposed model and method are validated against the literature of microcapsules, a new configuration of a relaxation phenomenon [START_REF] Gires | Transient behavior and relaxation of microcapsules with a cross-linked human serum albumin membrane[END_REF] is simulated. Due to the difficulties in simulating both experimentally and numerically this phenomenon, up to our knowledge this is the first simulation of this kind, where the nuclues is modeled as an elastic bulk, and it was possible thanks to the flexibility of the fully Eulerian framework.

The paper is organized as follows. First, in Section 2 we describe the physical problem, represented by a capsule immersed within an incompressible fluid, and its main parameters. Next, in Section 3, we briefly recall the formulation of the Eulerian model used to simulated this complex fluid-structure problem. In Section 4, we present the general numerical scheme used to discretize the FSI model along with a novel method to mitigate the degradation of the deformation vector in the case of large deformation and long-time simulations. The mathematical model, along with the proposed numerical method, is then validated with respect to other approaches present in the literature, in Section 5. In Section 6, the relaxation phenomenon is deeply analyzed considering the impact of an internal solid, and its stiffness, on the capsule elastic behavior. Finally, in Appendix A, we present a test case that is a modification of the sheared elastic sphere relaxation proposed in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF].

Problem statement

We consider an elastic capsule (see Figure 1) modeled as a thin elastic membrane, with an internal nucleus of radius a b and stiffness χ, immersed within an incompressible Newtonian fluid with viscosity µ and density ρ. The same viscosity and density are considered inside and outside the membrane. The radius a m of the capsule membrane is taken as a characteristic length of the problem. The capsule immersed in the fluid undergoes large deformations and its elastic properties are characterized by a non linear constitutive law depending on the shear modulus G s and the dilatation modulus K s . The fluid and the elastic capsule are strongly coupled and thus this phenomenon is a fluid-structure interaction problem. Therefore, the phenomena arising are governed by two dimensionless numbers: the Reynolds number

Re = ρU a m µ , (1)
and the capillary number

Ca = µU G s ( 2
)
where U is the fluid reference velocity.

The additional dimensionless parameters a b /a m and χ/G s are considered to study the effect of the internal nucleus on the capsule profile and deformation. In this work, we propose the original idea of using a fully Eulerian model for the simulation of capsules. The elastic deformations of the capsule are described by Eulerian fields that are advected by the fluid velocity and the associated elastic force appears as an Eulerian source term in the fluid equations.

Fully

Eulerian fluid-structure model

Eulerian description of the fluid

Herein, the incompressible fluid modeled through the incompressible Navier-Stokes equations: where

ρ(∂ t u + (u • ∇)u) -div(2µD(u)) + ∇p = F div(u) = 0 (3) Y (x, t) = ξ x = X(t, ξ) n(x, t) Initial configuration Ω 0 Deformed configuration Ω t
D(u) = 1 2 ([∇u] + [∇u] T
) is the strain rate tensor, p the pressure and F is a source term. When no structure is considered, the source term F is set to zero retrieving the original Navier-Stokes equations. In later sections, more details are given on the formulation of the elastic force term for both thin membranes and elastic solids.

Eulerian description of solid deformations

Let Ω 0 ⊂ R 3 be the reference configuration of a continuous medium and Ω t ⊂ R 3 the deformed configuration at time t. In order to describe the evolution of this medium in the Lagrangian frame we define the forward characteristics X(t, ξ) as the image at time t in the deformed configuration Ω t = X(t, Ω 0 ) of a material point ξ belonging to the initial configuration, i.e., X : R + × R 3 -→ R 3 (see Fig. 2). The corresponding Eulerian velocity field is defined as u : R 3 × R + -→ R 3 where

∂ t X(t, ξ) = u(X(t, ξ), t), X(0, ξ) = ξ, ξ ∈ Ω 0 . (4)
To describe the continuous medium in the Eulerian frame, we introduce the backward characteristics Y (x, t) that for a time t and a point x in the deformed configuration, gives the corresponding initial point ξ in the initial configuration, i.e., Y : R 3 × R + -→ R 3 (see Fig. 2). Since Y (X(t, ξ), t) = ξ, differentiating with respect to time and space in turn we have:

∂ t Y + (u • ∇ x )Y = 0, Y (x, 0) = x, x ∈ Ω t , (5)
and

[∇ ξ X(t, ξ)] = [∇ x Y (x, t)] -1 , for x = X(t, ξ). (6)
The relation [START_REF] Cottet | Eulerian formulation and level set models for incompressible fluid-structure interaction[END_REF] is the Eulerian equivalent of the characteristic equation (4). In addition, equation [START_REF] Cottet | Méthodes Level Set pour l'interaction fluidestructure[END_REF] allows to compute the gradient of the deformation in the Eulerian frame via Y . The next sections are devoted to the description of bulk and membrane elastic deformations with the backward characteristics Y .

Hyperelastic bulk models

For hyperelastic bulk material, the internal energy E b is a function of the deformation tensor ∇ ξ X. We focus in this paper on materials that are Galilean invariant and isotropic. With these assumptions it can be proven (see [START_REF] Holzapfel | Nonlinear Solid Mechanics. A continuum approach for engineering[END_REF]) that E b is expressed as a function of the invariants of the left Cauchy-Green tensor B = [∇ ξ X][∇ ξ X] T . This tensor is written in the Eulerian description with ( 6)

B(x, t) = [∇ x Y ] -1 [∇ x Y ] -T . (7)
In this work, we consider the neo-Hookean constitutive law (E b depends only on Tr(B)). The associated Cauchy stress tensor is given by σ b = 2χB where χ is the elastic modulus coefficient.

In this work, we aim at developing a numerical tool able to simulate elastic capsules with a bulk kernel inside. The neo-Hookean model allows to describe, with one parameter, elastic media subjected to large but moderate deformations. To model various aspects of the resulting non linear relation between stress and deformation, other elastic constitutive laws can be considered: Mooney-Rivlin, Saint Venant Kirchhoff, Ogden. In the final model the bulk solid will be immersed in a incompressible fluid and the fluid-structure interface will be captured by a level set function ϕ b : R 3 × R + -→ R which is advected by the Eulerian velocity field u:

∂ t ϕ b + u • ∇ϕ b = 0. ( 8
)
The force associated to the neo-Hookean solid is then given by

F b = div H ϕ b ε σ b , σ b = 2χB, ( 9
)
where H is a smoothed Heaviside function and 2ε is the width of the interface. Therefore ϕ b > ε corresponds to the fluid domain and ϕ b < -ε to the solid domain.

Hyperelastic membrane models

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

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

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

A = B - (Bn) ⊗ (Bn) (Bn) • n . ( 12
)
where B is given by [START_REF] De Brauer | A cartesian scheme for compressible multimaterial hyperelastic models with plasticity[END_REF]. This tensor measures surface deformations by projecting the 3D deformations (measured by B) on the surface Γ t (represented locally by n). The vector n is an eigenvector of A associated to the eigenvalue 0 so det(A) = 0. We denote by (λ 1 ) 2 and (λ 2 ) 2 the two other eigenvalues. The other invariants are used to define the following quantities

Z 1 = Tr(Cof(A)) = |λ 1 λ 2 |, (13)
Z 2 = Tr(A) 2 Tr(Cof(A)) = 1 2
λ 1 λ 2 + λ 2 λ 1 . (14)
The eigenvalues λ 1 and λ 2 correspond to the local deformation of the surface, therefore it is intuitive that Z 1 measures the local area variation whereas Z 2 measures the local shear variation. These results are demonstrated in details in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF]. Note also that Z 1 ≥ 0 whereas Z 2 ≥ 1. Now we introduce a constitutive law E that depends on the invariants Z 1 and Z 2 and the associated energy which is localized in a neighborhood of the membrane,

E m = Q E(Z 1 , Z 2 ) 1 ε ζ ϕ m ε dx. ( 15
)
Here Q is a box containing the membrane, ε the width of the interface and ζ is a cut-off function used to spread the interface near {ϕ m = 0}. Elastic energy is conservative so a variation of the membrane shape induces an elastic force (principle of virtual power) and is given by

F m = div (E ,1 (Z 1 , Z 2 )Z 1 C 1 + E ,2 (Z 1 , Z 2 )Z 2 C 2 ) 1 ε ζ ϕ m ε , (16)
where

C 1 = I -n ⊗ n, C 2 = 2A Tr(A) -(I -n ⊗ n), (17)
and E ,i represents the derivative of E with respect to Z i . Three constitutive laws to model membranes are used in this article:

• Evan-Skalak model:

E ,1 (r 1 , r 2 ) = K s (r 1 -1), E ,2 (r 1 , r 2 ) = G s , (18)
• Skalak model:

E ,1 (r 1 , r 2 ) = G s 2 (-r 3 1 +4r 1 r 2 2 -r 1 -2r 2 )+ K s 2 (r 3 1 -r 1 ), E ,2 (r 1 , r 2 ) = G s r 1 (2r 1 r 2 -1), (19)
• Membrane Neo-Hookean model:

E ,1 (r 1 , r 2 ) = G s r 2 - 1 r 3 1 , E ,2 (r 1 , r 2 ) = G s r 1 . (20)

Overall model

The elastic membrane and bulk media are immersed in a incompressible fluid modeled by the Navier-Stokes equations. The overall fully Eulerian model is given by

         ρ(∂ t u + (u • ∇)u) -div(2µD(u)) + ∇p = F b (ϕ b , Y b ) + F m (ϕ m , Y m ) div(u) = 0 ∂ t Y ℓ + (u • ∇)Y ℓ = 0 ℓ = b, m ∂ t ϕ ℓ + u • ∇ϕ ℓ = 0 ℓ = b, m (21)
The membrane force F m is given by ( 16) and the bulk force F b is given by [START_REF] Del Mercato | Biological applications of lbl multilayer capsules: From drug delivery to sensing[END_REF]. Note that a level set function ϕ ℓ and a backward characteristic field Y ℓ is needed for each media (bulk (b) and membrane (m)). These equations are completed with appropriate initial and boundary conditions that will be detailed in the numerical results sections.

p, Y, ϕ > > u 1 ∧ ∧ u 2 ∆x 2 ∆x 1

General scheme

The equations ( 21) are discretized with finite volume schemes on a staggered grid (see Figure 3 for a 2D configuration). In the following, the notation ϕ (respectively Y ) denote either ϕ b or ϕ m (respectively Y b or Y m ) because the same discretization is used for each variable. Let ∆t be the time step and u n , p n , ϕ n , Y n , ρ n , µ n the time discretization of the variables at t n = n∆t. The semi-discretization in time is given by the projection method (see [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] for an overview on projection methods) for the Navier Stokes equations (Steps 1-2-3) and an Euler explicit scheme for the advection equations (Step 4):

Step 1 :

ρ n u ⋆ -u n ∆t + div(u n ⊗ u ⋆ ) -div(2µ n D(u ⋆ )) + ∇p n = F (ϕ n , Y n )
Step 2 : div ∆t ρ n ∇ψ n+1 = div(u ⋆ )

Step 3 :

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

ϕ n+1 -ϕ n ∆t + u n+1 • ∇ϕ n = 0, Y n+1 -Y n ∆t + (u n+1 • ∇)Y n = 0.
In Step 1 a prediction of the velocity u ⋆ is computed with an Euler implicit scheme for the viscous term. The convection term is treated explicitly in time with NSSP3 Runge-Kutta [START_REF] Ketcheson | Highly efficient strong stability-preserving runge-kutta methods with low-storage implementations[END_REF] and in space with WENO5 [START_REF] Jiang | Efficient implementation of weighted eno schemes[END_REF], following an high-order momentum procedure [START_REF] Desmons | A generalized high-order momentum preserving (homp) method in the one-fluid model for incompressible two phase flows with high density ratio[END_REF], and the source term explicitly with a standard second-order discretization. In step 2, the Poisson equation for the pressure increment ψ n+1 is solved with appropriate boundary conditions [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. The resulting linear systems are solved with the GMRES algorithm of the HYPRE library [START_REF] Falgout | Hypre: a library of high performance preconditioners[END_REF][START_REF] Falgout | The design and implementation of hypre, a library of parallel high performance preconditioners. chapter in Numerical Solution of Partial Differential Equations on Parallel Computers[END_REF] with Jacobi preconditioning for Step 1 and with multigrid preconditioning for Step 2. In Step 3 the velocity is corrected to enforce the incompressibility condition and the pressure is updated. In Step 4 the transport equations are discretized with an explicit NSSP3 Runge-Kutta scheme in time and a WENO5 scheme in space for ϕ, Y . For the cut-off function, we considered the following expression ζ(r) = 1 2 (1 + cos(πr)) on [-1, 1] and ζ(r) = 0 elsewhere. The Heaviside function is defined as our simulations ε is fixed at 2∆x which is the standard value used in the literature to spread the interface.

H(r) = r -∞ ζ(x)dx.

Extrapolations, reinitilization and inner diffusion

It was already discussed in section 3.2 of [START_REF] Deborde | Numerical simulations of the interaction of solitary waves and elastic structures with a fully eulerian method[END_REF] that the numerical scheme associated to the fully Eulerian fluid-structure model with bulk (corresponding here to the terms indexed by b in ( 21)) can be unstable in severe situations. Indeed the backward characteristics Y and the level set ϕ are computed on the whole domain and these fields can become distorted and irregular in the vicinity of the fluid-structure interface. These perturbations can lead to numerical errors on the location and magnitude of the elastic force and therefore create instabilities. These limitations were overcome using a redistancing technique for the level set function [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF][START_REF] Russo | A remark on computing distance functions[END_REF] together with a linear Aslam extrapolation [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF] on the backward characteristics of the external fluid.

In numerical simulations of fluid-structure interaction with bulk materials, the backward characteristics can become distorted only in the outside fluid region. In the inside region, the characteristics are smooth because the elastic forces tend to bring them back to their initial configuration. However these characteristics can become irregular in the interior fluid in the case of a thin elastic membrane surrounded from either side by a fluid. To illustrate this property we consider a representative test case in capsule applications where an elastic membrane is immersed within a fluid where a linear shear configuration is considered. More information about the numerical simulation of this test case are given in Section 5.1. The membrane will elongate in one direction and then will rotate along itself in a stationary velocity field: this behavior is called the tank-treading motion. The isolines of Y are plotted in Figure 4a at the early stage of rotation. It is clear that for larger physical times, i.e when the membrane will perform several full rotation, the characteristics will become even more distorted. This eventually leads also to a bad representation of the membrane.

The first idea to address this problem was to use an Aslam extrapolation of Y inside the fluid with the information coming from the membrane. Since the membrane is closed, the characteristics will cross together and then give rise to a non smooth field. Instead, we propose in this paper to perform a method, we called inner diffusion: the backward characteristics are smoothed by using a diffusion equation in the internal region. In Figure 4b we present the same illustration where Aslam extrapolations, and the inner diffusion were used. We clearly see that the characteristics are smooth everywhere and this allows us to run more severe configurations for longer times.

For the reinitilization technique we solve the equation for a fictitious time τ ,

∂ τ φ + sgn(φ 0 )(|∇φ| -1) = 0
with a WENO5 scheme. For the linear Aslam extrapolation,

∂ τ Y n + H(φ)n • ∇Y n = 0, ∂ τ Y + H(φ)(n • ∇Y -Y n ) = 0
the normal derivative of Y is extrapolated firstly in a constant, and then in a linear, manner in the region φ > 0 (H is a Heaviside function) with the initial condition Y n (t = 0) = n • ∇Y . Both equations are discretized by an Euler explicit scheme in time and a WENO5 scheme in space.

For the inner diffusion approach, we solve the equation ∂ τ Y = ∆Y , in the region φ < 0, by discretizing it with an Euler explicit scheme in time and a centered second-order scheme in space. The fictitious time τ has been set for every test case depending on the configurations: more severe tests may need more diffusion.

Validation

In this section, two test cases have been setup with different physical parameters (size, elasticity) of the membrane in order to validate the proposed method with respect to other approaches present in the literature. One additional test case, relevant for the grid convergence of the proposed approach has been proposed in Appendix A starting from a similar test introduced in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF]. In particular, this academic test is shown to give an idea about the required mesh refinement needed to obtain a good enough resolution of the deformed capsule.

Simple shear flow

In this section, we present the numerical simulations for the simple shear flow test case. Thanks to the results available in the literature [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] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF][START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a stokes flow[END_REF], this test case allows us to validate the proposed method with complex configurations, typical for capsule applications, arising when a membrane is immersed within a linear shear velocity field. The complexity comes from the fact that the steady state consists in a constant shear occurring on both the inside and the outside of the membrane. The numerical challenge of such test cases is related to evolution of the deformation vector Y close to the interface. For this reason, it is mandatory for the long-time simulation of this test case to use both the Aslam extrapolation and the inner diffusion in order to capture properly the tank-treading motion (mentioned in Section 4.2). The simulation is setup by considering a thin membrane of radius a m , modeled with several constitutive laws, immersed within a viscous flow, characterized by the physical parameters in Table 1

, in the domain [-4, 4] × [-2, 2] × [-2, 2]. The dimensionless parameters for this configuration read Re = ρa 2 m γ µ = 0.0625, Ca = µ γa m G s , ( 22
)
where Re is the Reynolds number, and Ca is the capillary number, which is modified according to the configuration chosen. The initial conditions are given by Equation (A.1) to initialize the level-set of a sphere, and

u 0 (x, y, z) =   γy 0 0   , Y 0 (x, y, z) =   x y z   , (23)
ρ µ a m shear rate γ 1 4.0 0.5 1.0 to impose the linear shear velocity field and no pre-deformation on the initial configuration. We impose Neumann boundary conditions along the x-axis, moving wall along the y-axis, and slip wall along the z-axis. The computational domain is discretized using a uniform grid of 512 × 256 × 256 = 33 554 432 cells. The time step chosen for these simulations is set to ∆t = 2 × 10 -3 s. Further details about the configurations simulated in this section are given in Table 2, where different constitutive law types and multiple values of the capillary number are tested. The constitutive laws employed for these simulations are the Neo-Hookean (20) and the Skalak ( 19) laws. We remind here that, in general, the stretching modulus can be set as K s = 3G s . In these contributions, the shape of the membrane is assumed to be an ellipsoid once the steady state is reached and, by taking the slice z = 0 of this geometry, it is possible to introduce a dimensionless parameter depending on the characteristics of the ellipse,

D 12 = L 1 -L 2 L 1 + L 2 . ( 24
)
This parameter, computed by taking the major and minor axes of the ellipse (see Figure 5), has been mainly introduced for validation purposes. Here, the proposed approach is compared to the literature [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] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF][START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a stokes flow[END_REF]. In Figure 6 the geometry of the capsule is no longer changing, the structure dynamics is still evolving. It is then clear from the visualization that the steady state configuration consists in a fixed but rotating membrane. Figure 7 presents the results of D 12 when varying the capillary number. The results presented in this section are in good agreement with those coming from the literature showing a reasonable increasing trend of D 12 , when a larger capillary number is considered. In agreement with the results in Figure 7, it is possible to observe also in Figure 6 that using smaller capillary numbers leads to a stiffer membrane that tends to elongate less.

Capsule in a square-section channel

The test cases presented in this section have the goal of simulating capsules flowing in a squaresection channel, even when the capsule size at rest is larger than the channel dimension. In these situations, an interaction between the membrane and the boundary layer occurs and the capsule deforms into a parachute shape. The simulation is setup by considering a thin membrane, modeled with several constitutive laws, immersed within a viscous flow in a channel of dimensions [-1, 1] × [-0.5, 0.5] × [-0. initial conditions for both the level-set field ϕ and deformation vector Y are computed so that the capsule has a pre-deformed shape. This choice is related to the fact that, for some test cases, the capsule is assumed to be larger than the channel cross section. In this case, the initial conditions read ϕ 0 (x, y, z) = (xe

-2t 0 ) 2 + (ye t 0 ) 2 + (ze t 0 ) 2 -a m , (25)
where a m is the radius of the relaxed sphere, and

u 0 (x, y, z) =   V p -V 0 0   , Y 0 (x, y, z) =   xe -2t 0 ye t 0 ze t 0   , ( 26
)
where t 0 is the pre-deformation time, V is the mean capsule velocity and V p is a two-dimensional Poiseuille on a square section [START_REF] Pozrikidis | Introduction to Theoretical and Computational Fluid Dynamics. EngineeringPro collection[END_REF] such that,

V p = πV 2β ∞ i=0   1 (2i + 1) 3 - cosh (2i+1)πz ℓ (2i + 1) 3 cosh (2i+1)π 2   sin nπ y ℓ + 1 2 ,
where

β = π 4 96 - ∞ i=0 tan (2i+1)π 2 (2i+1) 5 π 2 .
The equations on ϕ 0 and Y 0 represent those of a deformed sphere into an ellipsoid elongated, by a virtual divergence-free velocity field, along the x-axis and compressed along the y and z axes. We impose inlet and Neumann conditions along the x-axis, and moving wall on the other boundaries. Multiple configurations are studied by varying the constitutive laws used to model the elastic properties of the membrane and the ratio a/ℓ, where 2ℓ = 1 represents the channel's size. Further details about the numerical parameters chosen in this section are given in Table 3. All simulations are setup by considering the following dimensionless parameters,

Re = ρa m V µ = 0.05, Ca = µV G s = 0.1, ( 27
)
where ρ is set to one. The goal of this benchmark is analyzing the deformation of the parachute shape with respect to the influence of the constitutive law and the ratio a m /ℓ. The computational domain is discretized using a uniform grid of 128 × 64 × 64 = 524 288 cells and the time step is set to ∆t = 10 -3 s. In Figure 8 we present the numerical results obtained when modeling the elastic membrane with Neo-Hookean constitutive law and varying the size of the immersed capsule. Figure 9 shows similar numerical results obtained by using the Skalak law, also with a capsule larger than the channel (TC14) to assess our method for highly compressed capsules. Due to the small size of the channel, the capsule is highly influenced by the higher velocity profile occurring in the middle of the channel that brings the structure to take an elongated shape that resembles that of a phantom. Increasing the size of the capsule leads to an even more elongated parachute, whose shape even features corners close to the boundary layer (see Figure 8b). Little differences in the computation of the membrane profile with respect to the literature are probably related to the different numerical approaches used for the simulations. Moreover, in our approach the forces are spread near the interface, which is not the case in [START_REF] Hu | Characterizing the membrane properties of capsules flowing in a square-section microfluidic channel: Effects of the membrane constitutive law[END_REF], hence regularizing the corners behind the parachute. In general, for all configurations a very good agreement with the literature is observed.

Numerical simulations of a capsule with an internal nucleus with complex geometries

In this section, we present the numerical simulations of a highly compressed capsule travelling within a square-section channel until it opens to a much wider space bringing the capsule to feature a relaxation phenomenon. A novelty with respect to previous works is the introduction of a solid bulk structure within the thin membrane. In particular, we perform thorough numerical studies on the influence of the size and stiffness of the bulk on the capsule relaxation.

The domain is divided in two part: the first one being a square-section channel [-1, 3] × [-0.5, 0.5] × [-0.5, 0.5] where the compressed capsule is initialized and the second one being the wider part [START_REF] Bergmann | An eulerian finite-volume approach of fluidstructure interaction problems on quadtree meshes[END_REF][START_REF] De Brauer | A cartesian scheme for compressible multimaterial hyperelastic models with plasticity[END_REF]×[-1, 1]×[-0.5, 0.5] where the relaxation happens (see Figure 10). The simulation is setup by considering a thin membrane, of radius a m , with an internal bulk of radius a b , immersed within a viscous flow. The internal square-section channel is treated as an immersed boundary by imposing wall boundary conditions through a classical penalization approach. ρ µ V a m G s K s law type 1 11 1 0.55 110 330 SK The dimensionless parameters used for this simulation are computed considering,

Re = ρa m V µ = 0.05, Ca = µV G s = 0.1, (28)
and are shown in Table 4, where 2ℓ represents the channel dimension with ℓ = 0.5. The initial condition on the velocity is given by the two-dimensional Poiseuille on a square-section given in Section 5.2. The level-set ϕ m and deformation vector Y m are initialized as

ϕ m,0 (x, y, z) = (xe -2t 0 ) 2 + (ye t 0 ) 2 + (ze t 0 ) 2 -a m , Y m,0 (x, y, z) =   xe -2t 0 ye t 0 ze t 0   , (29)
where t 0 = 0.2. The level-set ϕ b and deformation vector Y b for the elastic bulk are initialized as a sphere with no pre-deformation,

ϕ b,0 (x, y, z) = x 2 + y 2 + z 2 -a b , Y b,0 (x, y, z) =   x y z   . (30)
As mentioned above, several simulations have been performed by varying the size a b and the elastic modulus χ of the internal bulk. Further details about the setup of both a b and χ are given in Table 5. The computational domain is discretized using a uniform grid of 512 × 128 × 64 = 4 194 304 cells, and the time step chosen is ∆t = 5 × 10 -4 s. In Figure 11, we plotted the evolution of the membrane profile, with and without an internal bulk. At the beginning, as presented in Section 5.2, when flowing within the square-section channel, the compressed capsule takes a phantom shape. Then, once the capsule is at the end of the channel, it starts relaxing and modifying its shape until it becomes again similar to a sphere. To give more insights on the simulations, in Figure 12 we define three characteristic lengths to study the evolution of the parachute during the relaxation phenomenon. Figure 13 shows the evolution in non-dimensional time tV /ℓ of L 1 /ℓ and L 3 /ℓ when varying the bulk's size, while keeping a constant elastic modulus χ = 2200 (TC14, TC15, and TC16 in Table 5). As expected, in all cases L 1 /ℓ suddenly decreases and features a small jump when the capsule is exiting the channel. After that, the relaxation process begins and L 1 /ℓ increases again. Figure 13 shows that the presence of an elastic bulk slightly influences the evolution of L 1 /ℓ. Similar trends are observed for L 3 /ℓ, however it should be noticed that in this case the presence of an elastic bulk introduces remarkable changes in the parachute shape, which appears much more elongated for larger bulks. In particular, it is observed that when increasing the bulk's size, the minimum of L 3 /ℓ increases too. Although the dynamics of the relaxation process is highly influenced by the internal bulk, it is interesting to observe that the minimum value of L 3 /ℓ always occurs around tV /ℓ = 6.2. In Figure 14 we plot the evolution of L 2 /ℓ where the strong correlation with L 1 /ℓ is observed: when L 1 /ℓ decreases, L 2 /ℓ increases and vice versa.

Figure 15 shows the time evolution of the L 3 /ℓ characteristic length when varying the elastic modulus of the internal bulk, while keeping a constant radius a b /a m = 0.5 (TC14, TC16, TC17, TC18, and TC19). Similar trends to Figure 13b are followed when performing analysis on the bulk stiffness. In particular, it should be noticed that the stiffer the internal bulk is the higher the minimum of L 3 /ℓ gets. As already observed in Figure 13b, even though the capsule dynamics is highly affected by the presence of the internal solid, the minimum value of L 3 /ℓ always occurs at the same time. More details about the influence of an internal bulk on the capsule's shape is given in Figure 16, where the capsule dynamics is depicted at different times for different configurations (TC15, TC16, TC20, and TC21). For this case, an internal bulk with χ = 2200 is considered for the simulations. This brings the solid not to deform that much during the relaxation process. However, the bulk's size remarkably influences the shape of the parachute that appears differently. In particular, the smaller the bulk is, the more the capsule resembles the phantom shape shown in Figure 9a. Instead, when a bigger bulk is considered the shape of the parachute takes a more deformed profile.

Finally, in Figure 17 a fixed ratio a b /a m = 0.5 is considered, and the stiffness of the bulk varies from 110 to 2200. In this case, the bulk now shows larger deformations because a lower elastic modulus is chosen. This brings the parachute to slightly modify its shape assuming a more, or less, elongated appearance. also seems to be in line with the time evolution of L 3 /ℓ presented in Figure 15.

Conclusion and discussion

In this paper, we have presented a numerical approach for the simulation of elastic capsules. The numerical model is based on a fully Eulerian formulation for both the fluid and the deformable structures, membrane and solid volume, that are here taken into account by the Navier-Stokes equation by means of a source term describing the elastic forces, which act on the fluid. The fluid-structure system of equations considered herein describes both the forces related to the area variation and shear by following the full membrane model introduced in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF]. As mentioned above, the evaluation of the backward characteristics used to compute the elastic deformations can degrade when complex simulations of capsules are addressed, especially when membrane are considered. In this case, distortions of backward characteristics may occur within the membrane, causing a less reliable prediction of the capsule dynamics. The novelty of our method is based on coupling the Aslam extrapolation to improve the isolines of the backward characteristics outside the capsule (as done in [START_REF] Deborde | Numerical simulations of the interaction of solitary waves and elastic structures with a fully eulerian method[END_REF]) with the inner diffusion approach to smooth the internal isolines. Unfortunately, the good performances provided by the Aslam extrapolation in the outer zone does not solve the same problem in the internal part. The new algorithm is validated on several complex 3D configurations where the capsule dynamics is simulated and compared to other numerical experiments present in the literature. In particular, the challenges of these test cases are such that the proposed algorithm is essential to simulate the correct evolution of the capsule. First, we focus on the validation of the proposed model and method for classical test cases related to capsule dynamics. We studied before the deformation of a capsule when immersed into a linear shear velocity field, which gives rise to the tank-treading motion. Then, we performed several simulations of a capsule in a square-section channel, which allows to study different parachute shapes that the capsule takes when highly compressed. The numerical ρ µ G s K s a m 1 0.01 0.1 1 0.5 The membrane is pre-deformed with a 3D circular shear until a fictitious time t 0 = π. corresponds to the following initial conditions on Y ,

Y 0 (x, y, z) = 1 1 + α   x cos(t 0 z) + y sin(t 0 z) y cos(t 0 z) -x sin(t 0 z) z   (A.2)
where α is the pressurization coefficient. Following the Equation ( 14), the initial area and shear variations reads,

Z 1,0 = (1 + α) 2 , Z 2,0 = 1 + t 2 0 (x 2 + y 2 ) 2 2(x 2 + y 2 + z 2 ) . (A.3)
This introduces a shear variation that is equal to zero at the poles and varies along the z-axis (see Figure A.18). It should be noticed that, α does not influence the shear variation which stays the same, but introduce an area variation on the initial conditions. The coefficient α has been introduced to obtain a sphere shape after the relaxation process. Moreover, when taking α = 0 we obtain the same configuration analyzed in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF]. The computational domain is discretized using a Cartesian mesh with N = {64, 128, 256} cells in each direction. The time step chosen for the simulations varies depending on the mesh refinement: ∆t = 4 × 10 -3 for N = 64, ∆t = 2 × 10 -3 for N = 128 and ∆t = 10 -3 s for N = 256. In this simulation the pressurization coefficient α is set to 0.05. The numerical results at different time steps are presented in Figure A.19, showing the deformation isocontours until the new equilibrium is reached. At the beginning of the simulation, only the force related to the shear coefficient is acting on the membrane. However, the membrane inertia introduces a shape modification that brings into play also the force term related to the area variation. As expected, the shear deformation gives rise to a notable deformation and inertia, along both the x and z axes, such that the membrane begins to swing to one side and then to the other. From t = 3s, the membrane starts stabilizing and slightly oscillating around the steady state where Z 2 is uniformly equal one, and Z 1 = Z 1,0 . Because of the incompressibility condition, the internal fluid is going to be characterized by the same volume, that translates into an area variation at the equilibrium equal to the initial one. Having introduced the pressure coefficient α, the equilibrium configuration has an area variation greater than one that implies an elastic force at the interface. The pressure jump resulting from the elastic force is the reason why we are able to return to an actual sphere configuration at the equilibrium (which was not achieved in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF] where α was virtually set to zero). In order to provide more insights about these oscillations, in Figure A.20, we plotted the evolution of the horizontal and vertical radii. These curves also give an idea about the mesh refinement needed to achieve good enough results. The dimensionless pressure along the x and the z axes is also studied, in Figures A.

Figure 1 :

 1 Figure 1: Characteristic variables of the problem.

Figure 2 :

 2 Figure 2: Forward and backward characteristics.

Figure 3 :

 3 Figure 3: Staggered grid with position of unknowns

  In

Figure 4 :

 4 Figure 4: Deformation of the isolines of Y at early stage rotation in the linear shear illustration: (a) without Aslam extrapolations and inner diffusion; (b) with Aslam extrapolations and inner diffusion.

L 2 L 1 Figure 5 :

 215 Figure 5: Simple shear flow: steady state configuration of the deformed capsule, with L 1 and L 2 represent the principal directions.

Figure 6 :Figure 7 :

 67 Figure 6: Simple shear flow: three-dimensional visualization of the tank-treading motion at steady state.

  (a) TC9 : slice on the yz-plane. (b) TC10: slice on the yz-plane. (c) TC9 : slice on the xy-plane. (d) TC10: slice on the xy-plane.

Figure 8 :

 8 Figure 8: Capsule in a square-section channel: the membrane profiles obtained by the proposed method with Neo-Hookean constitutive law (dashed blue line) are compared to those presented in [22] (black dotted line).

  (a) TC11: slice on the yz-plane. (b) TC12: slice on the yz-plane. (c) TC13: slice on the yz-plane. (d) TC11: slice on the xy-plane. (e) TC12: slice on the xy-plane. (f) TC13: slice on the xy-plane.

Figure 9 :

 9 Figure 9: Capsule in a square-section channel: the membrane profiles obtained by the proposed method with Skalak constitutive law (red bold line) are compared to those presented in [22] (black dotted line).

Figure 10 :

 10 Figure 10: Capsule relaxation phenomenon: initial configuration of a membrane with an internal solid.

Figure 11 :

 11 Figure 11: Capsule relaxation phenomenon: evolution of the capsule profile at different times.

1 Figure 12 :

 112 Figure 12: Capsule relaxation phenomenon: L 1 , L 2 , L 3 are the characteristic length of deformed capsule.

  No bulk a b /a m = 0.6 a b /a = 0.5 No bulk a b /a = 0.6 a b /a m = 0.5 (b)

Figure 13 :

 13 Figure 13: Capsule relaxation phenomenon: time evolution of L 1 /ℓ and L 3 /ℓ when varying the size of the internal bulk, while keeping a constant elastic modulus χ = 2200. The numerical results of TC14 are depicted using a dashed line, TC15 in blue, and TC16 in red.

5 Figure 14 :

 514 Figure 14: Temporal evolution of the length L 2 /ℓ for three test case: elastic volume (T C15), a b /a = 0.6 (T C16) and a b /a = 0.5 (T C17). The elastic volume coefficient is χ 2200.

Figure 15 :

 15 Figure15: Capsule relaxation phenomenon: time evolution of L 3 when varying the elastic modulus of the internal bulk, while keeping a constant radius ratio a b /a = 0.5. The numerical results of TC14 are depicted using a dashed line, TC16 in blue, TC17 in red, TC18 in orange, and TC19 in green.

  (a) tV /ℓ = 5.75 (b) tV /ℓ = 6.25

Figure 16 :

 16 Figure 16: Capsule relaxation phenomenon: time evolution of the membrane and bulk profiles during the relaxation process when varying the size of the internal bulk, while keeping a constant elastic modulus χ = 2200. The numerical results of TC15 are depicted in blue, TC16 red, TC20 in orange, and TC21 in green. The walls of the square-section channel are represented in grey.

Figure A. 18 :

 18 Figure A.18: Relaxation of a sheared elastic sphere: initial condition of Z 2 .

  (a) t = 0.1s (b) t = 0.5s (c) t = 1.2s (d) t = 2s (e) t = 3s (f) t = 5s

Figure A. 19 :

 19 Figure A.19: Relaxation of a sheared elastic sphere: surface profile with the deformation isocontours Y 1 (red) and Y 3 (black) at different time steps.

21

 21

Figure A. 20 :Figure A. 21 :Figure A. 22 :

 202122 Figure A.20: Relaxation of a sheared elastic sphere: evolution of the x-radius and z-radius for N = {64, 128, 256}.

Table 1 :

 1 Simple shear flow: physical parameters.

  we show two three-dimensional visualizations of the tank-treading motion for different capillary numbers by plotting the velocity vectors of the fluid on the membrane surface. Here we show the flow velocity field over the membrane surface to point out that even if

	test case Ca	G s	K s law type
	TC1	0.15 13.33	-	NH
	TC2	0.4	5	-	NH
	TC3	0.6 3.33	-	NH
	TC4	0.9 2.22	-	NH
	TC5	0.15 13.33 40	SK
	TC6	0.4	5	15	SK
	TC7	0.9 2.22 6.66	SK
	TC8	2.0	1	3	SK

Table 2 :

 2

Simple shear flow: stretching and shear moduli, K s and G s , when varying the capillary number for the Neo-Hookean (NK) and Skalak (SK) constitutive laws.

Table 3 :

 3 Capsule in a square-section channel: stretching and shear moduli, K s and G s , when varying the ratio a m /ℓ for the Neo-Hookean (NK) and Skalak (SK) constitutive laws.

	5, 0.5]. The

Table 4 :

 4 Capsule relaxation phenomenon: simulation setup for the membrane.

Table 5 :

 5 Capsule relaxation phenomenon: simulation setup for different configurations of the internal bulk.

Table A . 6 :

 A6 Relaxation of a sheared elastic sphere: physical parameters.

results obtained for different capillary numbers and constitutive laws have been compared with other numerical techniques present in the literature, obtaining a good agreement in all cases. In Appendix A, we also present a grid refinement analysis for an academic test case to understand the needed level of refinement to achieve satisfying results. In particular, a modification of the circular shear benchmark introduced in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF] is proposed by introducing a pressurization coefficient that allows, with the right constitutive law, the correct recovery of the spherical shape after its relaxation. Finally, an original test case is presented where an elastic nucleus is considered within the thin membrane. We provide interesting insights about the influence of a nucleus on the deformation of a capsule during its relaxation after a compressed phase.

Appendix A. Relaxation of a sheared elastic sphere

This academic test case is a modification of the sheared elastic sphere relaxation proposed in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF]. Herein, we propose to use again this test case to validate the new method described in this paper by performing a mesh refinement analysis. As it was shown in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF], the relaxation of the pre-deformed sphere ends with a final configuration that does no longer resemble an actual sphere. This is mainly related to the fact that the sphere was initialized with no area variation (Z 1 = 1), which eventually brings the membrane to modify its aspect. The simulation is setup by considering a thin membrane, modeled with the Evans-Skalak constitutive law [START_REF] Gorsse | A simple cartesian scheme for compressible multimaterials[END_REF] immersed within a viscous flow in the domain [-1, 1] 3 . It should be noticed that also the constitutive law has changed with respect to the test case presented in [START_REF] Milcent | Eulerian model of immersed elastic surfaces with full membrane elasticity[END_REF]. The physical parameters used for this test case are given in Table A.6. We impose zero velocity as initial condition and Neumann boundary conditions on all boundaries of the domain. The initial conditions for the level-set field ϕ is given by a sphere of radius a m , ϕ 0 (x, y, z) = x 2 + y 2 + z 2 -a m .

(A.1)