Henri Bruneau
email: charles-henri.bruneau@math.u-bordeaux.fr

M Jedouaa
email: mjedouaa@yahoo.fr

C-H Bruneau

E Maitre

An efficient interface capturing method for a large collection of interacting bodies immersed in a fluid

Keywords: fluid/structure interaction, level set method, multiple bodies, collision model

published or not. The documents may come

An efficient interface capturing method for a large collection of interacting bodies immersed in a fluid

Meriem Jedouaa, Charles-

Introduction

Numerical simulations of fluid/structure and structure/structure interaction (FSI) have attracted an increasing interest and several methods have been proposed during the last decades.

A popular and wide spread method is the Arbitrary Lagrangian Eulerian approach (ALE) introduced by Donea in 1982 (see [START_REF] Donea | An arbitrary lagrangian-eulerian finite element method for transient dynamic fluidstructure interactions[END_REF]) and extensively studied by several teams [START_REF] Grandmont | Fluid-structure interaction: a theoretical point of view[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] Hron | A monolithic FEM/multigrid solver for an ALE formulation of fluid-structure interaction with applications in biomechanics[END_REF][START_REF] Richter | Finite elements for fluid-structure interaction in ale and fully eulerian coordinates[END_REF][START_REF] Wick | Fluid-structure interactions using different mesh motion techniques[END_REF]. The ALE strategy is an hybrid method that combines the Lagrangian and Eulerian descriptions using a mobile non structured grid that follows the normal displacement of the fluid/structure interface. The fluid and solid equations are solved individually and continuity conditions for the velocity and stress tensor are explicitly discretized at the interface. The main weakness of the ALE method is its difficulty of implementation, especially when dealing with large displacements in dimension three. In addition the added mass effect [START_REF] Conca | Added mass and damping in fluid-structure interaction[END_REF] has been a long standing difficulty which has been worked around only recently [START_REF] Fernández | A projection algorithm for fluid-structure interaction problems with strong added-mass effect[END_REF][START_REF] Fernández | Coupling schemes for incompressible fluid-structure interaction: implicit, semi-implicit and explicit[END_REF][START_REF] Fernández | Incremental displacement-correction schemes for incompressible fluid-structure interaction[END_REF][START_REF] Fernández | Fully decoupled time-marching schemes for incompressible fluid/thin-walled structure interaction[END_REF][START_REF] Fernández | Generalized robin-neumann explicit coupling schemes for incompressible fluid-structure interaction: Stability analysis and numerics[END_REF]. Moreover, the computational grid has to be remeshed when the elements get too distorted, which could be a very costly procedure in dimension three.

Another method, introduced by Cottet and Maitre in [START_REF] Cottet | A level-set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF][START_REF] Cottet | A level set method for fluid-structure interactions with immersed surfaces[END_REF], consists in using a purely Eulerian formulation for describing the fluid/structure interaction. This approach was inspired by the immersed boundary method of Peskin [START_REF] Peskin | The immersed boundary method[END_REF] where the forces at the interface were described in a Lagrangian manner. In the model [START_REF] Cottet | A level-set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF], a level set method is used to capture the interface. The level set method was developed in [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF] to treat problems involving interfaces. It is used in many domains because of its numerous advantages: its ease of implementation, topological changes are directly handled and one single level set function can capture an arbitrary number of interfaces. This last property is largely used as it can be a very efficient tool to capture several interfaces, for instance when merging and splitting of interfaces are allowed. The present work aims at dealing with a dense suspension of biological cells immersed in a fluid. In this kind of application, using one level set is not always sufficient as will be explained thereafter. The other existing level set models for capturing a large number of objects are either computationally expensive or cannot be used to handle collisions. Indeed, in [START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF] one level set function is required for each body, leading to a huge computational cost although contacts are easily avoided. In [START_REF] Wang | color level sets: a multi-phase method for structural topology optimization with multiple materials[END_REF], a formulation using log 2 N level set functions to represent N different regions is designed. This model, based on the four color theorem, substantially reduces the number of level set functions and can handle very easily complex topologies. However, reconstruction of all distance between bodies can not be achieved, as a consequence this model is not able to deal with several pairwise interacting bodies immersed in a fluid. Indeed, in that case, the investigation of fluid/structures interaction raises the problem of collisions between the bodies. These collisions can be handled if the distance between bodies can be computed from the localisation procedure. As confirmed by theoretical works [START_REF] Hillairet | Lack of collision between solid bodies in a 2d incompressible viscous flow[END_REF], the hydrodynamical forces between two smooth bodies evolving in a Navier-Stokes flow do prevent contact in finite time. However, during numerical calculations, the flow between too close particles is not accurately resolved which leads to contacts or inter-penetration of particles. The numerical handling of these contacts and overlaps hence is crucial for the simulation of particles suspensions. A natural approach proposed in [START_REF] Hu | Direct simulation of flows of solid-liquid mixtures[END_REF] is to refine the mesh in the inter particle gap in order to resolve accurately the flow fields. However, this strategy is highly expensive as several refinements are necessary and in the case of multiple particles the cost of the simulation is highly increased. Other techniques consists in overcoming this low resolution problem near contact by taking into account lubrication forces, when the gap between the bodies is very small [START_REF] Maury | A many-body lubrication model[END_REF][START_REF] Dance | Incorporation of lubrication effects into the force-coupling method for particulate two-phase flow[END_REF]. Due to the singular behavior of the forces, and the time discretization errors, this approach appears to be insufficient and might still lead to contacts and overlaps at low spatial resolution. Other numerical strategies, less respectful of the underlying physics, consist in imposing a constraint on the particle motion by means of artificial short range repulsive forces [START_REF] Glowinski | A distributed lagrange multiplier/fictitious domain method for particulate flows[END_REF][START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF] or by directly enforce a minimal distance between the particles [START_REF] Maury | A many-body lubrication model[END_REF]. By contrast to the refinement strategy, these collision methods allow, in addition of handling overlaps and contacts between particles, to use a coarser discretization, reducing substantially the computational cost compared to the method proposed in [START_REF] Hu | Direct simulation of flows of solid-liquid mixtures[END_REF].

In our work we consider such an algorithm, with short range repulsive forces, built on the multi geometric deformable model (MGDM) of interface capturing introduced by J. Bogovic [START_REF] Bogovic | A multiple object geometric deformable model for image segmentation[END_REF] for image segmentation. The proposed algorithm can handle multiple deforming bodies and avoid collision using a short range repulsive force depending on the distance to the closest interface, following [START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF]. The main advantage of this method is that it requires, whatever the number of interacting objects, only five fields to (1) locate and evolve each structure in the domain, (2) specify a speed or a force independently for each structure, (3) handle numerical contacts between the structures. This substantially reduces the computational cost, as will be illustrated below. More precisely, the level set defined in a neighborhood of the interfaces is transported with the fluid velocity. A multi-label fast marching method is then performed in a narrow-band around the interfaces allowing to update the label and distance functions. This combines the advantage of the MGDM method which efficiently captures a large number of bodies and their relative neighbours and of the collision model introduced in [START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF] using a level set decomposition. In order to validate the ability of this method to avoid numerical contacts and its efficiency to deal with a large number of structures, two applications are explored: rigid bodies and vesicles suspensions. This paper is organized as follows: Section 2 presents the proposed model with a careful description of the three label maps and the two distance functions. In section 3, the repulsive forces used to avoid collisions are described. Sections 4 and 5 shows how to apply the method to the case of rigid bodies immersed in an incompressible fluid. In section 6 is proposed a benchmark of numerical simulations starting with a qualitative study of the grid convergence on the sedimentation of rigid disks in two dimensions. Then, other simulations in two-and three-dimensions are performed and compared to the results of the literature. Section 7 is dedicated to the numerical simulations of immersed vesicles. Two-dimensional numerical results of vesicles under different flow conditions are presented. Finally, some conclusions are derived.

Domain labeling

In this section, we first recall some basic principles of the level set method. Then, we provide a description of the method used to capture multiple interfaces. This method, inspired by the multi geometric deformable model of J. Bogovic [START_REF] Bogovic | A multiple object geometric deformable model for image segmentation[END_REF], is introduced in the context of several bodies immersed in a fluid. The main idea is to partition the entire fluid/structures domain into several objects. In order to locate the different objects in the domain we introduce a set of label maps and distance functions. Let Ω ⊂ R d (d = 2 or 3) be a bounded domain which contains N bodies immersed in an incompressible fluid, we denote Ω i and Γ i the interior and the boundary of the i th body. By the term body we refer in this article to a simply connected domain of the ambiant space. We consider the fluid Ω N f as an object. Thus the fluid/structures domain Ω is partitioned into N + 1 objects as:

     ∀i = j, Ω i Ω j = ∅ Ω N f = Ω\{ N i=1 Ω i } Γ N f = N i=1 Γ i .
(1)

Outline of the level set method

Pioneered by Osher and Sethian in [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF], the level set method is very popular to address problems involving interfaces. It is widely used for numerical analysis of surfaces and shapes and in the context of fluid/structure interaction. The general idea of the level set method is to define a scalar function in the all computational domain that vanishes on the location of the interface to capture.

Let Ω be a bounded domain in R d (d = 2 or d = 3) partitioned into two sub domains Ω 1 and Ω 2 and Γ be the interface between Ω 1 and Ω 2 . The aim is to follow the evolution of the interface Γ that is defined as the zero level set of a smooth function φ, with a gradient modulus bounded away from zero in a neighborhood of this level set. At each time t, the interface Γ is characterized by:

Γ(t) = {x ∈ Ω, φ(x, t) = 0}
It is for example defined as:

   φ(x) < 0 x ∈ Ω 1 φ(x) = 0 x ∈ Γ φ(x) > 0 x ∈ Ω 2
Geometrical characteristics of the curve such as normal vectors n and curvature κ are obtained explicitly using the level set function:

n = ∇φ |∇φ| κ = ∇ • n.
The displacement of the interface is obtained by the evolution of the level set function φ. Let u be the velocity in the domain Ω, the level set function is the solution of the scalar transport equation:

∂ t φ + u • ∇φ = 0. ( 2
)
The velocity field u may depend on the space, the time, the geometric properties of the curve and/or the physics of the problem. For instance, if we consider a problem of fluid/structures interaction this velocity field is the velocity of the fluid or of the structures (assuming continuity of velocity at interfaces). To avoid situation where gradient modulus could become close to zero, we usually define the level set function as a signed distance function that is regular in each corresponding domain:

φ(x) = -d(x, Γ) x ∈ Ω 1 d(x, Γ) x ∈ Ω 2 (3)
where d(x, Γ) = min y∈Γ xy .

We then perform redistancing during the algorithm moving the interfaces.

Using an implicit function to capture the interface tells us directly to which region belongs a point x with the help of a Heaviside function H and its corresponding Dirac function ζ applied to φ(x). In practice, a regularized version of these functions H ε and ζ ε is used on the interface in order to reduce grid effects:

H ε (φ(x)) =          0 φ(x) ≤ -ε 1 2 (1 + φ(x) ε + sin( πφ(x) ε ) π ) |φ(x)| ≤ ε 1 φ(x) ≥ ε and ζ ε (φ(x)) = H ε (φ(x))
, where ε represents half of the interface thickness. Note that a classically mentioned advantage of the level set method is to handle automatically changes of topology. This property is problematic in the case of a collections of objects which are not supposed to merge (such as red blood cells), as splitting or merging of interfaces are directly taken into account by the level set function. Our algorithm will be built to prevent this unwanted merging.

Level set functions for multiple interfaces

Let us consider N pairwise disjoints bodies Ω i , i ∈ {1, ..., N } with Ω i ⊂ Ω, and let us denote Γ i the frontier of the cell Ω i and Ω N f the fluid background subdomain with N f = N + 1. A major advantage of the level set method is that one level set function can capture an arbitrary number of interfaces between the bodies and the fluid. Let φ be the level set function which captures the union of the N bodies Ω i , i ∈ {1, ..., N }, defined by:

φ(x) = -d(x, ∪ N i=1 Γ i ) x ∈ ∪ N i=1 Ω i d(x, ∪ N i=1 Γ i ) elsewhere
Using one single level set function is very time efficient for computations, since it is independent on the number of cells. However, it is not possible to specify a different velocity model and/or a force for each cell. If two bodies are too close the curvature and the normal are not well computed and a strong drawback is that we do not have any information on the distance between bodies. So this method can lead to collision and merging of bodies. Another way to capture multiple bodies is to use one level set function for each cell. Then each cell interface Γ i is captured by one level set function φ i . Thus, it is possible to specify a different speed or force to a body Ω i by using the level set function φ i and we get the distance between all the bodies. Indeed, the distance between two bodies Ω j and Ω i is given by the level set functions φ i and φ j as follows:

∀x ∈ Ω i , φ j = d(x, Γ j ) or ∀x ∈ Ω j , φ i = d(x, Γ i )
This multiple level set decomposition has been widely used in image segmentation and in fluid/structures interaction. A major disadvantage is its computational cost when the number of bodies increases.

Label maps

At every point x of the fluid/structures domain Ω, we define the label functions L 0 , L 1 , L 2 as:

∀x ∈ Ω, ∀i ∈ {1, ..., N + 1},          L 0 (x) = i if x ∈ Ω i L 1 (x) = arg min j =L0(x) d(x, Γ j ) L 2 (x) = arg min j / ∈{L0(x),L1(x)} d(x, Γ j ).
The label map L 0 provides a partition labelling of the whole computational domain Ω into N + 1 different objects. The label map L 1 identifies the index in this labelling of the first closest object at all points in Ω.

The label map L 2 identifies the index of the second closest object at all points in Ω.

     L 0 (x) = i if x ∈ Ω i L 1 (x) = j if the first closest object to x is Ω j L 2 (x) = k if the second closest object to x is Ω k .
In particular, for each point of a solid the label map L 1 corresponds to the index of the fluid. The label map function L 2 gives the index of the first closest body for any x in the structure domain, and the second closest body for x in the fluid domain.

∀x ∈ Ω i ,      L 0 (x) = i L 1 (x) = N f L 2 (x) = k
where Ω k is the first closest structure to x.

In cases where two structures are at equal distance to a point, the label is one of these structures. We will see thereafter that this is not a problem in our algorithm, since repulsion forces acting on two equally distant objects will move them so that is equality case does not hold anymore after one time step. Figure 1 shows an illustration of the three label maps in the case of five structures immersed in a fluid. Taking advantage of this local configuration of the closest object, one can define two related distance functions.

Distance functions

The distance functions associated to the first and second closest object are given by:

∀x ∈ Ω, ϕ 1 (x) = d(x, Γ L1(x) ) ϕ 2 (x) = d(x, Γ L2(x) ). ( 4
)
The distance function ϕ 1 (x) is the distance from x to the first closest object's boundary Γ L1(x) and ϕ 2 (x) is the distance from x to the second closest object's boundary Γ L2(x) . At any point of the domain Ω, ϕ 1 captures the union of all bodies interfaces and ϕ 2 provides the distance to the first closest body. As a consequence, on each point of a structure, we have the distance to the next closest one allowing to define a collision model to the closest interface, since avoiding contacts is equivalent to require:

∀x ∈ Ω , ϕ 2 (x) > 0.
Figure 2 shows an illustration of the two distance functions related to the label maps of the Figure 1.

Collisions strategy

As discussed before, it is crucial to develop a collision model to avoid contacts between bodies. In this section, we present the proposed collision model that consists in a short range repulsive force taking into account the interactions between the closest bodies. This short range repulsive force is inspired by the collision model introduced in [START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF].

Collision model using several level set functions

The collision model [START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF] was developed within a level set framework in the context of fluid/rigid bodies interaction. The level set decomposition is used in such a way that each body interface is captured by one level set function. We consider N bodies immersed in a fluid and we denote by F j,i the force applied by the body Ω j on the body Ω i and φ i the level set function which captures the boundary Γ i of the body Ω i . The distance of a point x of Ω i to the body Ω j is provided by φ j (x) and the direction of the force F j,i is obtained directly by ∇φ j . Moreover, to locate the interface Γ i we use a cut off function regularized on a thickness ε on each part of the interface. The short range repulsive force is formulated as:

∀x ∈ Ω, F j,i (x) = k ε ρ(x)ζ φ i (x) ε ∇φ j (x) φ j (x) exp - φ j (x) ε b (5)
Consequently, we obtain the following collision model:

∀x ∈ Ω, F global (x) = N i,j=1 i =j ρ(x) k ε ζ( φ i (x) ε ) ∇φ j (x) φ j (x) exp - φ j (x) ε b ( 6
)
where k is a repulsive constant proportional to the square of the relative velocities of the corresponding bodies just before collision and ρ denotes the density. The coefficient ε represents the half thickness of the interface on which the repulsive forces are applied and ε b corresponds to the rebound coefficient. In [START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF], this value is fixed to ε. The interaction forces decrease exponentially fast for far away structures, reducing the number of influent interacting neighbours. Nethertheless, this collision model accounts for all possible interactions between the N bodies. Consequently for N bodies captured by N level set functions, N 2 computations of the repulsive forces are required, which represents a huge computational effort for large N .

The proposed collision model using two distance functions

To reduce the high computational cost of (6), we adapt F global so that it depends only on the two distance functions ϕ 1 and ϕ 2 . We set:

∀x ∈ Ω, F label (x) = k ε ρ(x)ζ ϕ 1 (x) ε ∇ϕ 2 (x) ϕ 2 (x) exp - ϕ 2 (x) ε b (7)
The term ζ ϕ1(x) ε provides the location of the union of all interfaces. Thus, this force has its support on a subset Γ ε = {x ∈ Ω, ϕ 1 (x) ≤ ε}. In numerical computations, the intensity of the force is controlled by the repulsive parameter k = K(L 0 , L 2 ), this parameter is proportional to the relative velocities between closest bodies just before collision. This coefficient is computed as:

K(L 0 , L 2 ) = max(K min , | (U L2 -U L0 ) • ∇ϕ 2 |, ((U L2 -U L0 ) • ∇ϕ 2 ) 2 ) ( 8
)
where U L0 denotes the velocity of the body on which forces are applied, whereas U L2 is the velocity of the first closest body. The term K min corresponds to the minimal intensity of the force which is used when some particles are almost at rest and hence, relative velocities are negligible.

This collision model takes into account the interaction between the closest bodies at all points. Indeed, as ϕ 2 is the distance to the second closest object at all points of the fluid/structures domain, if a body is surrounded by other bodies the interaction of the other structures are taken into account on different part of its interface.

For instance, in Figure 1, the repulsive force applied to the blue surrounded body comes from the four surrounding bodies as the label map L 2 has partitioned it into four parts. On the yellow part, the force applied is the one exerted by the yellow body and so forth. On the four surrounding bodies, only the forces exerted by the blue body is applied. The advantage of this formulation is that we get rid of the sum in (6) leading to a considerable saving of computational cost. In this section, we consider the case ρ = 1 to compare the two collision forces.

Comparison of the two collision models

In this subsection, we compare the two collision models introduced above in two-dimensions and for ρ = 1. We consider the case of N circular rigid bodies having the same radius R and we set ρ = 1 to compare the two forces. The following result provides an estimate of the difference introduced by using the label force [START_REF] Conca | Added mass and damping in fluid-structure interaction[END_REF] instead of the collision model [START_REF] Wick | Fluid-structure interactions using different mesh motion techniques[END_REF]. As expected, the difference is smaller for large bodies and/or small ε.

Proposition 1. Assuming that the N bodies are disks, such as:

∀x ∈ Ω, ϕ 2 (x) ≥ ε then: F global -F label L 1 ≤ N 2πk α + O ε R exp - αR ε where α = √ 13-3 2 ≈ 0.3.
Proof. Let L i (x) denote the i th closest object to x and φ Li(x) the distance function associated to the object Ω Li(x) , it comes:

∀x ∈ Ω, ∀i ∈ {1, ..., N }, φ Li(x) (x) = d(x, Γ Li(x) ).
Using the assumption on the distance between disks, since the supports of the cut-off functions do not intersect, it holds:

∀x ∈ Ω, ζ ε (ϕ 1 (x)) = N j=1 ζ ε (φ j (x)),
Then, the model ( 6) can be written using these functions:

∀x ∈ Ω, F global (x) = N i=2 k ε ζ ε (ϕ 1 (x)) ∇φ Li(x) (x) φ Li(x) (x) exp - φ Li(x) (x) ε .
As defined in subsection 2.4, we have ϕ 2 = φ L2 and using (7) leads to:

∀x ∈ Ω, F global (x) -F label (x) = N i=3 k ε ζ ε (ϕ 1 (x)) ∇φ Li(x) (x) φ Li(x) (x) exp(- φ Li(x) (x) ε ). (9)
Note that these short-range repulsive forces tends to zero out of a cut-off radius, thus only the interactions with the N n = O(1) closest bodies of a given body will be numerically significant. However if we consider all interactions we get:

∀x ∈ Ω, F global (x) -F label (x) N i=3 k ε ζ ε (ϕ 1 (x)) ∇φ Li(x) (x) φ Li(x) (x) exp(- φ Li(x) (x) ε ). ( 10
)
Considering that ∀i ∈ {2, .., N n }, φ Li are distance functions we get:

F global -F label L 1 ≤ k ε ζ ε (ϕ 1 ) L 1 (Ω) Nn i=3 1 φ Li exp(- φ Li ε ) L ∞ (Ω) .
Moreover by definition we have:

∀x ∈ Ω, ∀i ∈ {4, .., N }, φ Li(x) (x) ≥ φ L3(x) (x).
leading to :

F global -F label L 1 ≤ (N -2) k ε ζ ε (ϕ 1 ) L 1 (Ω) 1 φ L3 exp(- φ L3 ε ) L ∞ (Ω) .
It is a classical result from level-set theory that

1 ε ζ ε (ϕ 1 (x)
) is an approximation of the length of the zero level-set of ϕ 1 . Namely, in the case of N disks, we can prove that

1 ε ζ ε (ϕ 1 ) L 1 (Ω) = N (2πR + O(ε)),
Moreover, therefore, We are looking for a lower bound of φ L3 depending on the radius R. To that aim, we consider the worst case of three disks which are enlarged by ε and are in contact like shown in Figure 3. This configuration gives the minimum distance required between two disks. Taking a point x on the boundary of one disk (without loss of generality, Γ 1 ) as shown in Figure 3, we compute the minimal distance of x to the second closest disk. Once again, the worst case is obtained in the situation where x is as in Figure 3. Through a simple calculation, we obtain as distance:

F global -F label L 1 (Ω) ≤ N (2πRk + O(ε)) 1 φ L3 exp(- φ L3 ε ) L ∞ (11)
d = √ 13 -3 2 R
that is a lower bound for φ L3 . As r → 1 r exp(-r ε ) is decreasing, we obtained the announced estimation.

Remark 1.

1. Let us point out to the reader that this difference tends to zero when R tends to ∞ or when ε tends to zero. Consequently, we can adjust ε depending upon R such as this difference becomes negligible. Numerically, the coefficient ε depends on the spatial discretization step and is taken equal or larger than R 10 . Moreover, the influence of the repulsive force imposed by the first closest cell at all points F label , is the most influential on the dynamics of the cells as it is the largest.

2. Moreover, as previously stated, due to the exponential decay the repulsion forces magnitudes, only the closest bodies will have significant contribution to the global repulsion at some point of a given body. Therefore conceptually the leading coefficient N in ( 11) could be replaced by N n , where N n is the maximal number of bodies that could be in contact with a given body. This is dependent of the dimension and could depend also on the body if those have different sizes.

For a more in-depth study of the collision model proposed, we refer to Meriem Jedouaa PhD Thesis ( [START_REF] Jedouaa | Interface capturing methods for interacting immersed objects[END_REF], in English) where the full fluid-structure models are compared in the case of Stokes flows.

Evolution of the label maps and distance functions in the general case

This section is devoted to the evolution of the three label maps and the two associated distance functions. We present the evolution algorithm in two-dimensions, the extension to the 3D case is straightforward.

The evolution is based on the transport of one level set function that captures the union of all interfaces, then a multi label fast marching method is performed enabling a re-initialization of the label maps and distance functions. As the proposed collision model is a short-range repulsive force, we only evolve L 1 and ϕ 1 in a thin band close to the interfaces, we denote by d L1 the size of this narrow-band from each side of the interface. The computation of L 2 and ϕ 2 is activated only if these thin bands meet (just before collision) and is achieved in a thin band of size d L2 . The final step consists in updating the level set function φ as a signed distance function using the new updated distance function ϕ 1 close to the interfaces. This procedure is summarized as:

1. Transport a level set function φ which captures the union of all interfaces with the fluid velocity U, 2. Evolve the label maps L 0 , by redefining its values near the interfaces, 3. Perform a multi label fast marching to update L 1 , ϕ 1 in a thin band around the interfaces, 4. If bodies are nearby, perform a multi label fast marching to update L 2 , ϕ 2 .

Global level set field φ

We define the level set function φ which captures the union of all interfaces as:

φ(x) =          -C x ∈ ∪ N i=1 Ω i , d(x, ∪ N i=1 Γ i ) > d L1 -d(x, ∪ N i=1 Γ i ) x ∈ ∪ N i=1 Ω i , d(x, ∪ N i=1 Γ i ) ≤ d L1 d(x, ∪ N i=1 Γ i ) x / ∈ ∪ N i=1 Ω i , d(x, ∪ N i=1 Γ i ) ≤ d L1 C x / ∈ ∪ N i=1 Ω i , d(x, ∪ N i=1 Γ i ) > d L1 (12)
where C >> 1. This level set field is a signed distance function close to the interfaces, it is negative inside each structure and positive outside. Denoting by U the fluid velocity in the entire computational domain Ω, the level set function is solution of the following scalar transport equation:

∂ t φ + U • ∇φ = 0 φ(x, 0) = φ 0 (x) (13)
in a neighboorhood of the interfaces.

Redefinition of the label L 0

The evolution of the label map L 0 depends on the sign of the level set function. By definition the region where the level set function is non negative corresponds to the fluid:

∀x ∈ Ω, ∀i ∈ {1, ..., N }, L n 0 (x) = i = N f if φ n (x) < 0 N f if φ n (x) ≥ 0 ( 14
)
At each time step, we change the label value L 0 , near the interface at the points where the condition ( 14) is not verified. Namely, if the level set function is non negative, we set the new value of L 0 to the label of the fluid (L 0 = N f ) and if the level set function is negative and the label function L 0 was the label of the fluid, we assign the new value of L 0 to the value of the neighbour of x whose label is different from N f . We consider here that the level set function is transported with an explicit scheme under a CFL condition. For all point x = (i, j) we denote by Neighbours(x) the set which contains the four closest points to x:

Neighbours(x) = {(i + 1, j), (i -1, j), (i, j), (i, j -1), (i, j + 1)}.
The whole process is sketched as:

Algorithm 1 Redefinition of L 0 for (x ∈ Ω) do if φ n+1 (x) ≥ 0 then L n+1 0 (x) = N f else if (L n 0 (x) = N f ) then for (y ∈ Neighbours(x)) do if (L n 0 (y) = N f ) then L n+1 0 (x) = L n 0 (y) end if end for end if end if end for
This simple strategy is allowed only if the distance between two bodies is strictly greater than two spatial discretization steps. It will have to be taken into account by our collision model.

Multiple label fast marching method

We describe here the local fast marching algorithm that allows to update ϕ 1 , ϕ 2 and the label maps L 1 and L 2 in a vicinity of the interfaces. This local fast marching is an extension of the fast marching method [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF] that was introduced in [START_REF] Sifakis | Moving object localisation using a multi-label fast marching algorithm[END_REF] and [START_REF] Bogovic | A multiple object geometric deformable model for image segmentation[END_REF]. To this end, we solve two eikonal equations of the form:

∀x ∈ Ω, |∇d(x)| = 1
This equation is solved using a first order numerical scheme [START_REF] Rouy | A viscosity solutions approach to shape-from-shading[END_REF] and has to be supplemented by a boundary condition which depends on ϕ 1 or ϕ 2 .

Computation of L 1 and ϕ 1

The first step of the proposed procedure is the computation of the distance function ϕ 1 and the values L 1 by solving the following equation:

   ∀x ∈ Ω, |∇ϕ 1 (x)| = 1 ∀x ∈ Γ, ϕ 1 (x) = ϕ 1 (x) | ∇φ(x) |
This process is executed in a narrow-band of size d L1 .

At initialization, we set ϕ 1 (x) to

ϕ 1 (x) | ∇φ(x) |
close to the interfaces, this allows to propagate all the interfaces simultaneously.

There are three sets of points: alive(A), narrow-band (NB) and far away (F). At initialisation, the Alive set contains all the points of the N interfaces. The Narrow-Band set consists of the closest points to the interfaces. All others points are considered as far away and have to be computed and added to the narrow band.

Moreover, each contour is associated to the number of objects. To this end, we introduce a function lab which corresponds to the number of the interface that spreads.

During the propagation of the interfaces, the point x m for which ϕ 1 (x m ) is the lowest value of ϕ 1 in the narrow-band, is definitely deleted from the narrow-band and this point is considered as alive. This minimum value corresponds to the distance ϕ 1 (x m ) and we assign the label lab(x m ) to L 1 (x m ). Then, the values of ϕ 1 at the neighbors of the point x m are computed using the alive points which have the same label value, and these points are added to the Narrow-Band. The algorithm stops when the distance ϕ 1 has reached the narrow-band value d L1 . The steps performed during the initialization phase are provided by the algorithm2.

• d(y) = ϕ 1 (y) | ∇φ(y)) | • lab(y) = L 0 (x) • alive(y) = true The values L 1 (y) = lab(y) and ϕ 1 (y) are fixed. end if else • lab(x) = -1 • ϕ 1 (x) = -1 • alive(x) = false end if end for end for Algorithm 3 Iterative phase: computation of ϕ 1 and L 1 dist=0 while ((N B = ∅) and (dist ≤ N B L1 )) do Find x m such as d(x m ) = min (x∈N B) d(x) dist = max(dist, d(x m )) Delete x m from the NB if (L 1 (x m ) = lab(x m )) then alive(x m )=true end if if (alive(x m )=false) then ϕ 1 (x m ) = d(x m ) L 1 (x m ) = lab(x m ) alive(x m )=true for (y ∈ Neighbours(x m )) do ajout=true if (L 0 (y) = lab(x m )) then ajout=false end if if (L 1 (y) = lab(x m )) then ajout=false end if if (ajout=true) then
Compute d(y) using Neighbors(y) such as L 1 (Neighbors(y))=lab(x m ) Add y in NB end if end for end if end while

Computation of the label map L 2 and ϕ 2

As the procedure to update L 1 and ϕ 1 is only performed in a narrow-band of size d L1 , the computation of L 2 and ϕ 2 will be activated between two particles only if two or more particles are at a distance less or equal to 2d L1 .

If that is the case, we solve the following equation:

∀x ∈ Ω, |∇d(x)| = 1 ∀x ∈ Γ L1 , d(x) = ϕ 2 (x) (15)
in a narrow-band of size d L2 . Γ L1 corresponds to the intersection of the closest objects in the fluid domain, each point x ∈ Γ L1 is defined by:

∀y ∈ Neighbours(x), L 1 (y) = L 1 (x), L 0 (x) = L 0 (y)
Only the points close to the interface Γ L1 are visited as the algorithm stops propagation when the distance function ϕ 2 reachs the distance value d L2 .

Figure 4 provides an illustration of the results obtained for the label map L 1 and L 2 for a configuration of 20 rectangular particles and d L1 = 10h and d L2 = 5h. We observe that the computation of L 2 is achieved only for the closest structures allowing to reduce the CPU time. Moreover, the short-range repulsive force is not exerted on some particles as the distance is greater than 20h thus if ε = h the exponential term of the force is lower than 10 -9 .

Computational Complexity

We provide here the computational complexity of the multi-label fast marching (MLFMM) algorithm presented above. To this end, we add particles in the computational domain Ω next to each other so that the particles are close enough leading to the computation of L 2 as soon as there are two bodies in the domain.

Application to rigid bodies

In this section, we propose an application of the model to rigid bodies immersed in an incompressible fluid. First, we describe the model adopted for the numerical simulations as well as its numerical resolution. This

Volume penalization models

We consider N bodies evolving in an incompressible fluid denoted Ω N f and the entire domain Ω is partitioned in the same way as [START_REF] Donea | An arbitrary lagrangian-eulerian finite element method for transient dynamic fluidstructure interactions[END_REF]. In this case each body Ω i , 1 ≤ i ≤ N represents a rigid body. A volume penalization method is adopted to take into account the immersed rigid bodies.

Flow configuration

We denote by ρ f and µ the constant density and the viscosity of the fluid and by U and p the flow velocity and the pressure. The flow is governed by the viscous incompressible Navier-Stokes equations:

ρ f (∂ t U + (U • ∇)U ) -∇ • (µ∇U ) + ∇p = 0 in Ω N f × (0, T ) ∇ • U = 0 in Ω N f × (0, T )
with the condition at the interfaces of the solid bodies:

∀i ∈ {1, ..., N }, U = u i on Γ i ( 16
)
where u i is the rigid motion of the body Ω i . In order to get rid of the boundary conditions ( 16) on the solid boundaries, we use a penalization method ( [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF][START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF][START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF]). This method consists in adding a penalization term in the Navier-Stokes equations to impose the rigid motion inside the solid and to solve the boundary value problem inside the whole domain Ω including the bodies. So, there is no need for an adapted mesh to the geometry of the bodies and a Cartesian mesh on a box domain can be used. If we denote by χ i the characteristic function of the body Ω i and by λ the penalization parameter we obtain the following model:

     ρ(∂ t U + (U • ∇)U ) -∇ • (µ∇U ) + ∇p = ρg + λ( N i=1 χ i (u i -U )) + F global + F wall in Ω T = Ω × (0, T ) ∇ • U = 0 in Ω T (17)
complemented with homogeneous Dirichlet boundary and initial conditions, where g is the gravity force, F global is the collision force [START_REF] Conca | Added mass and damping in fluid-structure interaction[END_REF] and F wall is a repulsive force carry on by the solid boundaries of the domain Ω on the rigid bodies. The rigid velocity u i is obtained by averaging translation and angular velocities over the solid body Ω i ( [START_REF] Patankar | A formulation for fast computations of rigid particulate flows[END_REF]):

u i = u t i + u θ i = 1 | Ω i | Ωi ρχ i udx + J -1 i Ωi ρχ i u × (z -x g i )dz × (x -x g i ) (18)
where J i is the inertial matrix of the body Ω i and x g i its center . As the rigid bodies move with the fluid, the displacement of χ i is:

∀x ∈ Ω, ∀i ∈ {1, ..., N }, ∂ t χ i + U • ∇χ i = 0. or ∂ t χ i + u i • ∇χ i = 0.
Finally, denoting by ρ i the density of the immersed bodies Ω i , we obtain the following density function:

ρ = ρ f + N i=1 (ρ f -ρ i )χ i .
Let us note that the coefficient λ >> 1 means that the regions occupied by the solid bodies are considered as porous media with a very small permeability. Numerically, the penalisation parameter λ is fixed to 10 10 .

The penalization model using several level set functions

The characteristic function χ i can also be deduced considering a level set function φ i that captures the interface Γ i . This level set function evolves as:

∂ t φ i + u i • ∇φ i = 0
and using the regularized Heaviside function χ i can be smoothed as:

χ ε i = 1 -H ε (φ i ).
So the model becomes:

         ρ(∂ t U + (U • ∇)U ) -∇ • (µ∇U ) + ∇p = ρg + λ( N i=1 χ ε i (u i -U )) + F global + F wall in Ω T ∇ • U = 0 in Ω T ∂ t φ i + u i • ∇φ i = 0 in Ω T ( 19
)
This system is dependent on the number of structures causing a high computational cost if one wishes to simulate a large number of bodies.

The proposed penalization model

In order to alleviate the high computational cost due to a large number of particles N we take advantage of the label maps and distance functions to reformulate the penalization model.

To this end, we define a label map L 0,1 which allows to partition the computational domain into N subdomains, each subdomain is associated to the first closest structure, this function is defined as:

∀x ∈ Ω, L 0,1 (x) = L 0 (x) if (L 0 (x) = N f ) L 1 (x) otherwise (20)
The regularized characteristic functions of the solid bodies are then expressed as:

∀x ∈ Ω, ∀y ∈ Ω, χ L0,1(x) (y) = 1 -H( φ L0,1(x) (y) ε )
Denoting by ρ L0,1(x) the density of the body Ω L0,1(x) we obtain the following density function:

ρ x (y) = ρ f (1 -χ L0,1(x) (y)) + χ L0,1(x) (y)ρ L0,1(x) (y)
To compute the penalization term, the values of the rigid velocities are only required inside the particles and in a vicinity of the particle at a distance ε. Taking this in consideration and the advantage of the formulation above, we define a global rigid velocity which includes the N rigid velocities U i , i ∈ {1, ..., N }. More precisely, for each x ∈ Ω, y → U L0,1(x) (y) is the rigid velocity of the solid body Ω L0,1(x) obtained by averaging the translation and angular velocities over the solid.

Setting:

∀x ∈ Ω, | Ω L0,1(x) |= Ω L 0,1 (x) ρ(z)dz = Ω ρ(z)χ L ε 0 (x) (z)dz (21)
we obtain the following formulation . ∀x ∈ Ω, ∀y ∈ Ω,

U L0,1(x) (y) = 1 | Ω L0,1 | Ω ρ L0,1(x) (z)χ L0,1(x) (z)U (z)dz + J -1 L0,1(x) Ω ρ L0,1(x) (z)χ L0,1(x) (z)U (z) × (z -x g L0,1(z) )dz × (y -x g L0,1(y) ). ( 22
)
where J L0,1(x) and x g L0,1(x) are the inertial matrix and center of gravity of solid Ω L0,1(x) . In a neighboorhood of a particle Ω i that is when L 0,1 = i, the defined global velocity corresponds exactly to the rigid velocity U i .

Finally, the proposed penalized model is provided by:

         ρ(∂ t U + (U • ∇)U ) -µ∆U + ∇P = ρg + λ(χ L0,1 (U L0,1 -U )) + F label + F wall in Ω T ∇ • U = 0 in Ω T ∂ t φ + U • ∇φ = 0 in Ω T (23)
This model has to be completed with the initial and boundary conditions on the velocities and pressure.

On the implementation side, all rigid velocities are computed incrementally, involving only one iteration on the mesh grid. Thanks to the label maps, we have suppressed the dependence on the number of bodies in the repulsive force and in the penalization term. The only dependence on the number of bodies comes from the storage of the N level set functions. Therefore, this is a very desirable model to simulate a large number of interacting bodies. In the case of circular rigid bodies, the level set equation is replaced in the previous model ( 23) by the transport of the gravity center of the N bodies and a reconstruction of the associated level set function is performed:

∂ t x g i = u t i , 1 ≤ i ≤ N

Discretization and Numerical implementation

This part elucidates the numerical resolution of the system [START_REF] Glowinski | A distributed lagrange multiplier/fictitious domain method for particulate flows[END_REF], two algorithms are considered for the evolution of the rigid bodies and we provide a comparison of these two algorithms in terms of computational complexity.

Resolution of a multi-fluid problem in the fluid/structures domain

The Navier-Stokes equations are solved using an incremental projection method of Chorin type. Given a time step ∆t, we set t n = n∆t and U n ≈ U (., t n ). First, we compute an intermediate State U * from:

U * = U n -∆t(U n • ∇)U n + ∆tµ ρ n ∆U n + ∆tg - ∆t∇p n ρ n . ( 24
)
Then we solve the pressure from the equation:

∇ • ( ∇p n+1 ρ n ) = div(U * ) ∆t + ∇ • ( ∇p n ρ n ) ( 25
)
so that the velocity:

U = U * - ∆t ρ 0 (∇p n+1 -∇p n ) (26)
is divergence free. The pressure equation ( 25) can be solved directly by a conjugate gradient algorithm or approximated with a relaxation procedure as follows:

∆p n+1 = ρ 0 ∆t div(U * ) + ∆p n setting 1/ρ n = 1/ρ 0 -(1/ρ 0 -1/ρ n
) in order to get a Poisson equation much faster to solve. In that case, the step u satisfies Navier-Stokes equations with a modification of the pressure term (1/ρ n -1/ρ 0 )∇p n + 1/ρ 0 ∇p n+1 instead of 1/ρ n ∇p n+1 . It remains to choose ρ 0 in order to have the best approximation. We have tested three different values ρ 0 = ρ f , ρ 0 = ρ L0,1 and ρ 0 = (ρ L0,1 + ρ f )/2. The Poisson solver is based on a classical 5 or 7 points second order stencil according to the dimension, the viscous terms are discretized by a second order central scheme and the convection term is discretized by a 5 th order WENO scheme.

To find the best value for ρ 0 , we perform a numerical test with six rigid disks falling under gravity in two dimensions. The Figure 6 shows that the results obtained with ρ 0 = ρ f gives the worst approximation whereas the two other values are very close to the approximation of the exact solution computed with the conjugate gradient. In the following we set

ρ 0 = (ρ L0,1 + ρ f )/2.
Concerning the penalization term an implicit treatment is adopted in order to use larger penalization coefficient λ and therefore the interface boundary condition is satisfied with better accuracy.

Penalization algorithm 1

In the general case the algorithm performs the following steps:

1. Compute and add the repulsive force

U col = U + ∆tF label , 2. Compute the translation velocity U t L0,1 of the body Ω L0,1 , U t L0,1 = Ω ρ n U col • χ L n 0,1 dx Ω ρ n χ L n 0,1 dx , 3.
Compute the angular velocity W L0,1 of the body Ω L0,1

W L0,1 = J -1 L0,1 Ω ρ n r n × U col • χ L n 0,1 dx × r n , where r n = (x -x g L n 0,1 ) 4. Compute the rigid velocity U L0,1 of the body Ω L0,1 , U L0,1 = U t L0,1 + W t L0,1
5. Correct the velocity using an implicit treatment of the penalization term

U n+1 -U col ∆t = 1 λ χ L0,1 (U L0,1 -U n+1 ), 6.
Transport the solid bodies with the fluid velocity or the rigid velocity

φ n+1 = φ n -∆tU n+1 • ∇φ n , 7. Redefine the distance function ϕ n+1 1 = |φ n+1 |, Redefine L n+1 0 ,L n+1 1 using φ n+1
Perform the multi label fast marching method described in Chapter 3 8. Redefine the gravity centers as:

x g (L 0,1 ) = Ω L 0,1 ρX Ω L 0,1 ρ (27)
The implicit treatment of the penalization term allows to use a large penalization coefficient λ. Moreover with this algorithm, the incompressibility constraint is imposed before the rigidity constraint. The two constraints are better imposed as:

D(U ) = 0 ⇒ div U = 0
Penalization algorithm 2

The second algorithm consists in reconstructing the N level set functions using the fact that objects are rigid. The center of gravity is moved according to the velocity field, and the distance function the corresponding solid is recovered by a rigid motion of the initial one. Then, the steps 6 and 7 of the previous algorithm are modified as follows:

6. Transport the solid bodies by a rigid motion from the initial condition.

7. Redefine the label maps and distance functions:

∀x ∈ Ω, ∀i ∈ {1, ..., N + 1},          L n+1 0 (x) = i if φ n+1 i ≤ 0 L n+1 1 (x) = arg min j =L n+1 0 (x) φ n+1 j L n+1 2 (x) = arg min j / ∈{L n+1 0 (x),L n+1 1 (x)} φ n+1 j ∀x ∈ Ω,      ϕ n+1 1 (x) = d(x, Γ L n+1 1 (x) ) ϕ n+1 2 (x) = d(x, Γ L n+1 2 (x) )

Comparison of the two algorithms of evolutions

To compare the CPU time of both algorithms, we fill the computational domain with particles next to each other so that the computation of L 2 is active as soon as there are two particles in the domain. We fix the narrow-band sizes S L1 to 10h and S L2 to 5h. where h is the space discretization step. Table 1 shows the CPU time of both algorithms, we can see that the CPU time of algorithm 2 is faster. However, the difference is small and the Algorithm 2 is provided in the case of circular rigid disks so that no interpolations of the initial level set functions are required. For general shape of bodies, interpolations of the N initial level set functions have to be achieved leading to an increase of the CPU time. It is hence most efficient to use the algorithm 1 for general shape of bodies.

Numerical illustrations

In this section, we present the numerical results obtained with the proposed model and give some comparisons with existing methods of the literature. In a first part, we investigate the space convergence of the method. Then we compare our model [START_REF] Glowinski | A distributed lagrange multiplier/fictitious domain method for particulate flows[END_REF] to the model [START_REF] Hillairet | Lack of collision between solid bodies in a 2d incompressible viscous flow[END_REF]. In a third part, we provide a qualitative comparison of our model with those proposed in [START_REF] Lefebvre | Fluid-particle simulations with freefem++[END_REF]. Finally, some simulations of dense suspensions of rigid bodies in two and three dimensions are presented. For all of the simulations presented in this paper, the dynamic viscosity µ is set to 0.01, the density of the fluid is set to ρ f = 1, and the density of the rigid bodies is the same for all the bodies. Let k = (k x , k y , k z ) t be the repulsive coefficient between bodies and k wall = (k wall

x , k wall y , k wall z ) t be the repulsive coefficient exerted by the walls, the value of the components depends on the amplitude of the force.

Grid sensitivity

A grid convergence is carried out in dimension two and three. In dimension two, the grid convergence is performed on four grid levels (G 1 , G 2 , G 3 , G 4 ): (128 × 128), (256 × 256), (512 × 512) and (1024 × 1024) cells on a uniform mesh. The regularization parameter ε is fixed to ∆x G1 where ∆x G1 denotes the mesh size corresponding to the coarsest grid (128 × 128). We take as test case the sedimentation of 25 circular rigid bodies having the same radius R = 0.025 + ε. The gravity force g is set to -980 and the repulsive coefficients on the bodies and on the walls are:

k x = -g/10, k y = -g/10, k wall x = -g/40, k wall y = -g.
The repulsive forces are applied on a ring around the interface of thickness 1.5∆x G1 .

In Figure 7 the results obtained with the four different grids are presented, the white line corresponds to the numerical size of the rigid particles (φ = ε). We can see that the dynamic of the rigid bodies and the interaction between them is quite similar for the different resolutions. The collisions are avoided thanks to the repulsive force. We observe the same phenomenon of kissing and tumbling of the bodies. By t = 0.48, all simulations have reached static equilibrium which represents different local minimum of the sedimentation of the 25 disks. We observe that the two fine simulations are very close to each other until time t = 0.15 and keep a symmetrical distribution. At the end of the simulation, a packing of the bodies is formed at the bottom of the computational domain and is composed of three layers. Each layer contains the same number of structures for the four different grids: eleven bodies on the first one, ten bodies on the second one and four on the last one. The distribution of the four bodies of the third layer is different for the various grids but is much closer on the two finest simulations. For the three dimensional case, we study the grid sensitivity using three different grids which contain respectively : (64 3 ), ( 1283 ) and ( 2563 ) bodies on a uniform mesh. In the three dimensional case, = -g. The results obtained with the finest grids are almost similar, the particles fall down symmetrically no interaction occurs between the particles before hitting the bottom wall. On the contrary, with the coarser grid, the particles interacts and falls unsymmetrically, this is due to the fact that the short range repulsive forces are active from the beginning as there are less discretization points in the interstitial gaps between the particles. Thus a 3D qualitative grid convergence is achieved for a large number of particles.

Comparison of the method with level set decomposition

We give here a comparison of our penalization model ( 23) and the penalization model that uses a level set decomposition (19)

Numerical efficiency

We investigate here the numerical efficiency of the two proposed algorithms (5.1) and (5.2) in the case of N rigid disks. The computational time is compared to the N level set method.

The algorithm associated to our model is given in section 5.2. Instead of transporting one level set function and performing the multi label fast marching method, we transport the gravity centers of the N rigid structures and reconstruct their associated level set functions. Then, the label and distance functions are reinitialized by using their definition. It is obvious that this algorithm depends on the number of bodies because of the N level set functions and their associated center of gravity. However, this part of the algorithm is very fast. This algorithm will allow us to compare the saving computational time which is induced by changing the existing collision model ( 6) and the penalization term. We average the computational time on the ten first iterations.

The averaged CPU time of our algorithm (table 2) is compared to the method using N level set functions (table 1), according to the number of bodies. As noticed before the collision model (6) computes N 2 repulsive forces which induced a high computational cost as shown in the second column of table 1 whereas in the present algorithm, the CPU time of the collision model is constant as it does not depend on the number of bodies. The CPU time of the penalization model is larger in [START_REF] Hillairet | Lack of collision between solid bodies in a 2d incompressible viscous flow[END_REF] because it depends on the number of bodies. Indeed, N rigid velocities must be computed to get the right velocity of each body. In the present algorithm, we must add the label redefinition that depends almost linearly on the number of bodies and so is quite cheap. For a low number of bodies the CPU time of both methods is close but increasing this number from 2 to 400 bodies induces a total CPU time 8000 larger for the model [START_REF] Hillairet | Lack of collision between solid bodies in a 2d incompressible viscous flow[END_REF]

Numerical comparison of the two collision models

To highlight the differences between the two collision models we first focus on a a test case with three circular rigid bodies falling on each other. The simulations are performed on a grid of size (128 × 128) corresponding to a spatial discretization step ∆x = 7.8125.10 -2 in order to better see the difference between the two models. The bodies have the same radius R = 0.1 and the thickness of the interface is ε = 2∆x. The repulsive coefficients are:

k x = -g/10, k y = -g/10, k wall x = -g/40, k wall y = -g.
in Figure 11, black and white lines represents the interfaces of the three different objects. The black line stands for the collision model ( 6) and the white line for the collision model ( 7) corresponding also to the colors that represent the values of the level set function. We can see that the bodies have the same behaviour, as expected, because on the one hand the radius is large and on the other hand the forces applied on the bodies are very similar. So the difference is small. The Figure 12 shows the results obtained with 6 rigid bodies. Here we can see again that the behaviour is very similar for a large radius. The last test concerns the same configuration with six smaller disks with radius R = 0.03 leading to a higher difference in the dynamics even if the final State is the same. The difference between the two models is stronger when the number of body is larger or when the force coefficients are higher. The bottom plot in Figure 14 shows the vorticity inside the fluid domain, when the bodies reached the bottom they move to the right and induce a strong positive vortex that has a strong influence on the dynamics of the bodies in its turn. We observe that the vorticity increases when the bodies are close to each other as the repulsive force is higher. At time t = 0.5 the vorticity vanishes because the six rigid bodies have reached the static equilibrium.

Comparison with the model introduced in [31]

445

In this part, we compare qualitatively our model to an existing method proposed in [START_REF] Lefebvre | Fluid-particle simulations with freefem++[END_REF]. In that model, the solid bodies are taking into account by penalizing the strain tensor to enforce the rigid body motion (see [START_REF] Janela | A penalty method for the simulation of fluid-rigid body interaction[END_REF]). A scheme for inelastic collisions is implemented imposing a minimal distance between bodies and therefore avoiding contacts (see [START_REF] Maury | A time-stepping scheme for inelastic collisions[END_REF] for more details). Moreover, the contacts with the four walls are also handled in the same way. The test case is the sedimentation of 100 rigid particles of radius R = 0.01 subject to the gravity force g = -70. The corresponding repulsive coefficient are:

k x = -g/7, k y = -g/7, k wall x = -g/28, k wall y = -g.
The Figure 15 shows the results obtained with the FreeFem code implemented by A. Lefebvre (see [START_REF] Lefebvre | Fluid-particle simulations with freefem++[END_REF]) on a mesh with x elements. Using our model, the simulations are performed on a grid of size (512 × 512), the half interface thickness ε is set to the mesh size ∆x and the results are shown in Figure 16. The computation on the FreeFEM++ code was performed on a non-symmetrical mesh so that the results in Figure 15 are unsymmetric from the beginning whereas our results, performed on a cartseian grid stay symmetric until time t = 1.44. However, we can see that the dynamics is globally the same and that the static equilibrium is reached at the same time t = 4.8.

Comparison of the repulsive force with or without constant intensity

In this section is shown the dynamic of rigid bodies of various size, the numerical contacts are either handled by using the short-range repulsive force with a constant intensity or according to the relative velocity between closest bodies.

The simulations deal with the sedimentation of 30 rigid bodies for the case of disks the radii are R = 0.05 and R = 0.025 and we also considered rigid squares each side is equal to √ πR.

The simulations are performed on a grid of size (512×512) and the half thickness of the interface is ε = 1.5∆x.

The coefficient of gravity g is set to -980. The repulsive coefficients are fixed to 100 in the constant intensity force case.

Figures 18 and17 show the sedimentation of 30 rigid squares of different sizes. In the presence of rigid particles, the hydrodynamical instabilities are more prominent, as initially the particles are closed, the symmetry loss is almost immediate. Once this flow symmetry lost, we can not observe the same dynamics for these two simulations, as due to instabilities the flow exhibits a very complex behavior, and it seems obvious that the trajectories of rigid bodies can not be the same for the two simulations. This is even more visible in Figures 20,[START_REF] Hillairet | Lack of collision between solid bodies in a 2d incompressible viscous flow[END_REF] where the vorticity field is represented at different time. At first the vorticity is created by the falling of the bodies whereas, later, the vortices convect the bodies, the bodies are then driven by the fluid flow.

Nevertheless, we point out that globally the behavior of the rigid bodies is similar: we observe that until t = 0.17 for both collision model the obtained results are similar, the particles are carried by the fluid flow and turns in the domain convected by the created vortices. At t = 0.25 the minimal distance between the rigid bodies is smaller for the simulation represented in Figure 17 than the particles in Figure 18, there are more interactions between the particles. At t = 3.18, the steady state has been reached for both simulations, the configuration of bodies are different. With the collision model depending on the relative velocities, the particles have fulfilled the bottom of the computational domain whereas with the constant model we observe three layers of particles not completely filled.

The same simulation for rigid disks having the same area has been carried out. The results are presented in Figure 21. We only colored six particles to better track their motions, we can see that at time t = 0.03 the models provide exactly the same results, then at time t = 0.08, we observe a slight delay between the two simulations for the orange and light blue particles. This slight delay is sufficient to change all the particles trajectories. However, a collective motion of the big green particle and the smaller red disks is observed for both simulations. Moreover, globally the motion of the bodies is similar with the two models, the bodies are carried by the fluid, we observe the same process of kissing, tumbling of particles. At the end of simulations, the final configuration is different, here again with the model depending on the relative velocities, the particles have completely fulfilled the domain. The model using the relative velocities allow the particle to be closed, meaning that the repulsive constant coefficient may not be well calibrated.

Dense suspensions of rigid bodies in 2D and 3D

In this part, we present some results of dense suspensions of rigid bodies evolving in a fluid which was performed using our numerical model. The first simulations deal with the sedimentation of 400 rigid bodies of radius R = 0.01 in the two dimensional case. The simulations are performed on a grid of size (512 × 512) and the half thickness of the interface is ε = 1.5∆x. The white line shows the real numerical size of the particles corresponding to the isoline φ = ε. The coefficient of gravity g is set to -980, the associated repulsive coefficients are -g/10. The 400 bodies fall down symmetrically to reach a dense repartition at the bottom as can be seen in Figure 22. The second simulation addresses the 3D case. Figure 23 show the simulation of 500 rigid spheres of radius R = 0.01 falling under gravity using a grid of size 128 3 . The half thickness of the interface is ε = ∆x. The coefficient of gravity g is set to -980. At initial step, there are five slices of 100 bodies at a distance d = 0.1 (distance of two closest bodies' centers). In that case the repulsive forces are negligible. The interactions between bodies occur at once after t = 1.5 on the fine mesh. As a consequence, the equilibrium state is reached fastly around t = 2.5 instead of t = 9.4.

The last simulation deals with the dynamics of 90 rigid disks suspended in a shear flow. The intensity of the short-range repulsive force is proportional to the relative velocities between the closest bodies. The computational domain Ω is a rectangular domain of size [0, 2] × [0, 1], the grid resolution is (512 × 256), the regularization parameter ε = 1.5∆ x 5.9.10 -3 . Figure shows the obtained results at different time, the disk are colored according to the values of the label map L 0 and the background color corresponds the magnitude of the velocity. Until time t = 0.25, the rigid bodies moves horizontally through the canal. Then, we observe a convective motion, the solids are carried by the fluid and rotates in the domain. Thanks to the collision model, the numerical contacts between particles are avoided. This test case confirms the ability of the proposed collision model to handle collisions only by using the relative velocities of the closest cells.

Application to vesicle simulation

The purpose of this section is to validate the proposed model which includes a short-range repulsive force in the case of deformable bodies. An application to immersed vesicles suspensions is proposed. We first describe the Eulerian fluid/elastic membrane coupling model. Then, we present the numerical procedure used for the simulations. Finally, a validation of the model in the case of multiple vesicles under various flow are presented.

Eulerian fluid/elastic membranes coupling model

We present here the Eulerian fluid/elastic membrane coupling model introduced in [START_REF] Cottet | A level-set formulation of immersed boundary methods for fluid-structure interaction problems[END_REF] in the case of multiple vesicles and using the label maps and distance functions.

Computation of the elastic and bending forces

In the context of multiple vesicles, the use of one level set function to capture the union of all interfaces requires to maintain a sufficiently large minimal distance between the vesicles in order to compute correctly the high order derivatives of the level set function which appear in the elastic and bending forces.

To reduce the number of grid points required in the narrow-gap between two interfaces, we use the distance function ϕ 2 to compute the forces of each vesicle Γ ε i in the regions where the level set function φ is not defined as the signed distance function to Γ i . By doing this, the stencil required to compute a geometrical property at a point x belonging to Γ ε i is prolongated to the regions where L 2 = i. Numerically, we introduce a function which indicate the points where this procedure has to be applied. This procedure is used in order to compute the curvature and other derivatives of order higher than two of the level set function φ.

Collision model

The short-range repulsive force is used here, to avoid numerical contacts between vesicles. Moreover, to compute accurately the curvature and elastic forces on the membranes, a minimal distance has to be maintained between the membranes which corresponds to the interfaces Γ ε i of thickness ε. We hence fix the activation distance of the repulsive force ε b to 2ε. The intensity of the short-range repulsive force is computed according to the relative velocities between closest vesicles just before collision. Each vesicle velocity is computed by averaging the fluid velocity on each membrane Γ ε i , so that setting:

∀x ∈ Ω, | Γ L0,1(x) |= Γ L 0,1 (x) dz = Ω 1 ε ζ ε (φ L0,1(x) (z))dz (28)
we obtain the following average velocity over the membrane boundary:

∀x ∈ Ω, U a (L 0,1 (x)) = 1 | Γ L0,1 | Γ L 0,1 (x) U (z)dz = Ω 1 ε ζ ε (φ L0,1(x) (z))U (z)dz (29)
then, at all points x ∈ Γ ε i , i = (1, .., N ) where the label map L 2 (x) is defined, the intensity of the force is computed as:

K L2(x),L0,1(x) =| (U a (L 2 (x)) -U a (L 0,1 (x))) • ∇ϕ 2 (x) |
where U a (L 2 (x)) denotes the averaged velocity associated to the vesicle Ω L2(x) , this term provides the averaged velocity of the second closest vesicle to x. Then, at each point x where Ω L0,1(x) and Ω L2(x) are at a distance ε b this value is fixed. Finally, the short range repulsive force is computed as:

∀x ∈ Ω, F label (x) = | (U a (L 2 (x)) -U a (L 0,1 (x))) • ∇ϕ 2 (x) | ε ρ(x)ζ ε (ϕ 1 (x)) ∇ϕ 2 (x) ϕ 2 (x) exp - ϕ 2 (x) ε b (30)

Complete model

The fluid-structure interaction system hence reduces to the incompressible Navier-Stokes equations with source terms (elastic, curvature forces) combined with a scalar transport equation of the level set function φ and the stretching e.

Find (U, P, φ) solution of :

         ρ(φ)(∂ t U + (U • ∇)U ) -div(µ(φ)D(U )) + ∇P = F e + F c + F label in Ω T = Ω × (0, T ) div u = 0 in Ω T = Ω × (0, T ) ∂ t φ + U • ∇φ = 0 in Ω T = Ω × (0, T ) ∂ t e + U • ∇e = -e(n(φ) ⊗ n(φ)) : D(U ) in Ω T (31)
Denoting by ρ L0,1(x) the density of the vesicle Ω L0,1(x) we obtain the following density function:

ρ x = ρ f (1 -χ L0,1(x) ) + χ L0,1(x) ρ L0,1(x) , µ x = ρ f (1 -χ L0,1(x) ) + χ L0,1(x) µ L0,1(x)
This model involved several parameters, in order to achieve numerical simulations, it is most convenient to introduce dimensionless parameters. For this purpose, we introduce reference values. Let L r , U r , ρ r and µ r denote the characteristic lenght, velocity, density and viscosity scales. Considering the following characteristic quantities :

             x = L r x , y = L r y , z = L r z U = U r U , P = ρ r U r L r P , t = L r U r t µ = µ r µ , ρ = ρ r ρ , ε = L r ε φ = L r φ , ϕ 1 = L r ϕ 1 , ϕ 2 = L r ϕ 2
Differentiating (and dropping '), we obtain the following dimensionless system:

           ρ(φ)(∂ t U + (U • ∇)U ) - 1 R e div(ν(φ)D(U )) + ∇P = 1 W e Fe (φ) + 1 W c Fc (φ) + F label in Ω T div U = 0 in Ω T ∂ t φ + U • ∇φ = 0 in Ω T ∂ t e + U • ∇e = -e(n(φ) ⊗ n(φ)) : D(U ) in Ω T (32)
where the physical parameters are given by: R e

= L r U r ρ r µ r , W e = µ r U r λ , W c = µ r U r L 2 r α 9.2.

Numerical procedure

The system is discretized by a finite difference method on a staggered grid (MAC type) with a finer resolution for the level set function. The fluid velocity and pressure are discretized on the coarse grid while all the other functions (label maps, distance functions and level set functions) are located on the finer grid. To compute the fluid velocity a Projection method is used, the diffusion term are treated implicitly and an explicit scheme is used for the forces and the convective term.

Stability conditions

The advection equation of the level set function is treated explicitly and to ensure the stability of the scheme, we impose the stability condition CFL:

∆t ≤ ∆x |U | ∞
The explicit treatment of the fluid/elastic coupling induces stability conditions. The Brackbill condition provided in [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] gives a stability condition in the case of perfect fluids separated by an interface with surface tension, expressed as:

∆t ≤ W e 2π (∆x) 3 2 (33)
where W e decreases while the surface tension coefficient increases. This condition is valid for high Reynolds number, where viscosity of the fluid could be neglected. In the case of small W e , which is the case of red blood cells, the time step has to be chosen very small. In the case of small Reynolds number, [START_REF] Galusinski | Level-set method and stability condition for curvature-driven flows[END_REF] exhibited the stability condition:

∆t ≤ W e R e ∆x (34)
An unifying condition has been proposed by [START_REF] Cottet | A semi-implicit level set method for multiphase flows and fluid-structure interaction problems[END_REF][START_REF] Bost | Méthodes level-set et pénalisation pour le calcul d'interactions fluide-structure[END_REF] where the stability of schemes for immersed elastic membranes have been studied.

For the curvature force, no condition exists, one can however take into account the typical velocity induced by the bending force and build a maximal time-step on the above conditions.

Algorithm

The level set function φ and the distance functions are updated at each time step using the multi-label fast marching on the finer grid. The narrow-band size set for the simulations are: d L1 = 10h where h corresponds to the discretization space step associated to the coarser grid (the fluid grid) and we fix a high narrow-band size for the computation of ϕ 2 and L 2 we set d L2 = 10h this allows to compute correctly the bending and elastic forces when the vesicles get closer by means of the reconstituted level set functions. At each time step, the algorithm performed the following steps:

1. Computation of the fluid velocity U n+1 using an implicit scheme and a projection method of Chorin type according to the forces F n e , F n c , F n

Numerical illustrations

This section is dedicated to the numerical results obtained with the proposed model. We first provide a grid convergence of the algorithm then numerical simulations of multiple vesicles in several type of external hydrodynamical flow are presented.

Grid sensitivity

A grid convergence is carried out in the case of two vesicles evolving in a Poiseuille type flow. The computational domain is a rectangle of size [0, 2] × [0, 1].

To study the grid convergence the simulations are performed on three grid levels (G 1 , G 2 ,G 3 ) which contain respectively : (256 × 128), (512 × 256), (1025 × 512) cells on a uniform mesh. The regularization parameter ε is fixed to ∆x G1 where ∆x G1 denotes the mesh size corresponding to the coarsest grid ∆x G1 0.0078. The three grids used to discretize the level set functions, label maps, distance functions and the stretching are twice finer than the three coarsest grid.

The velocity profile of the Poiseuille viscous fluid flow is parabolic and at initialization the shape of vesicles are Cassini Oval defined using the implicit function:

((x -a) 2 + y 2 )((x + a) 2 + y 2 ) = b 4
For both vesicles we set the parameters a = 0.18 and b = 0.1747, the associated reduced area is around 0.8. Using a fast marching procedure we initialize φ to a signed distance function to the interfaces defined as the minimum of the two implicit functions. The physical parameters fixed in this study are: R e = 1.10 -1 , W e = 0.05, W c = 200 At initialization, the vesicles are located at the entrance section of the domain. As illustrated in Figure 26, their shape progressively change during their motion in the canal. Figure 25 shows the results obtained with the three resolutions at different times. We can see that the deformation of the vesicles are similar for the three resolutions, we note however that at time t = 0.5 the deformation of the two finer grids are closer than the one obtained with the coarsest grid. At the end of the simulations, the shape of the right vesicle corresponds to the equilibrium shape of a vesicle having a reduced area around 0.8 whereas the second vesicle adopts a parachute shape. In light of these results, the grid convergence is achieved for the coarsest grid.

Vesicles passing through a bifurcation

We consider here the case of vesicles passing through a bifurcation, the velocity has a parabolic profile at the entrance and outlet sections.The computational domain Ω = [0, 2] × [0, 1]. In the case of a bifurcation, the proposed collision model failed at avoiding numerical contacts between vesicles and the wall. Instead, the intensity of the force is fixed, numerically we found that a coefficient k = 10 is sufficient to avoid contacts. We set the following parameters: R e = 0.01, W e = 0.005, W c = 30

The simulation represented in Figure 27 is performed on a grid of size (1024 × 512). The initial shape of the vesicle is an ellipse of size a = 0.1 and b = 0.4, the height of the bifurcation is around four times smaller than the vesicle height (b = 0.8). We can see that the vesicle progressively deforms while passing through the bifurcation. During the deformation, the vesicle adopts different shape in order to pass through the channel. Finally, at the end of the simulation the vesicle has a parachute shape which is a typical shape of vesicle deformed in a Poiseuille type flow. , as the size of the vesicles is small, the simulations are performed using a fine resolution, the grid contains (2048 × 1024) cells on an uniform mesh grid. As the level set grid is twice finer, the induced computational time of this simulation is very high, around 50 seconds per iteration. The obtained results are represented in Figure 28, colors indicate the value of the label map L 0 ,the white lines correspond to the isolines φ = 0 . At initialization, the region occupied by the vesicles represents around one-half of the computational domain, the configuration of the vesicle is seven layers of 15 vesicles. Each vesicle interface corresponds to a Cassini Oval with the parameters a = 0.076 and b = 0.08. We observe that depending on the layer, the vesicles adopt different shapes, these shapes are on one hand due to the pressure driven Poiseuille flow, to the elastic and bending forces and the interactions between the vesicles. At t = 0.14, we can see that the bottom blue left vesicle and the top red left vesicle advanced slowly than the other vesicles and hence interacts with their neighbors, the symmetry is lost and while the cells advances in the canal various shapes can be observed.

Conclusions

In this paper, we introduced a new model to simulate efficiently a large number of interacting bodies immersed in a fluid. This model involved three label maps and two distance functions which allow to locate the bodies and their closer neighbours in the domain. A collision model depending on the distance between the closest bodies is proposed. This model which is totally independent on the number of bodies, is compared both theoretically and numerically to the model introduced in [START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF]. We present an application to rigid structures with a penalisation model that only depends on five advected field functions. Numerical results are in good agreement with the results of the literature at least qualitatively. Compared to a model which is totally dependent on the number of bodies, our model substancially reduces the CPU time. A numerical test on bodies of various radii shows that the collision model is efficient even when the strength of the force is very different. An application of the method to immersed deformable vesicles have been proposed. Using the level set function φ, the stretching e and the label maps one elastic and bending forces are computed to impose the inextensibility of the whole set of vesicles membranes as well as their resistance to bending. To handle numerical contacts, the relative velocities are computed by using the averaged fluid velocity on each vesicle membrane. A convergence study in the case of two vesicles evolving in a Poiseuille is presented. The dynamic behavior of vesicles under different boundary conditions have been presented and confirms the capability of the proposed method to deal with numerical contacts between vesicles at low resolution.

Figure 1 :

 1 Figure 1: Illustration of the three label maps for a configuration of five bodies, from left to right: L 0 , L 1 and L 2 . Each object has a specific color and the red one corresponds to the fluid. The white contour represents the boundary of the bodies.

Figure 2 :

 2 Figure 2: Illustration of the two distance functions for a configuration of five structures, from left to right: ϕ 1 and ϕ 2 . The white contour represents the boundary of the circular bodies.

Figure 3 :

 3 Figure 3: Configuration of three bodies. The circles represented are of radius R + ε.

Algorithm 2

 2 Initialization phase for (x ∈ Ω) do for (y ∈ Neighbours(x)) do if (L 0 (x) = L 0 (y)) then points at the interfaces if (alive(y) = false) then • Add y to the Narrow-Band

Figure 4 :

 4 Figure 4: Illustration of the multi label fast marching algorithm for 20 objects. The black color corresponds to the undefined values of L 1 and L 2 . The white contours corresponds to the zero value of the level set function φ

  The CPU time is represented in Figure5 for different number of particles close to each other for d L1 = 10h 270 and d L2 = 5h where h denotes the spatial discretization step.

Figure 5 :

 5 Figure 5: CPU time of the proposed multi label fast marching algorithm.

Figure 6 :

 6 Figure 6: Flow field for six rigid disks falling under gravity (left) and a zoom of the location of the interface for different values of ρ 0 (right). The conjugate gradient method is in red, ρ 0 = ρ f in black, ρ 0 = ρ L 0,1 in white and ρ 0 = (ρ L 0,1 + ρ f )/2 in blue. The background color shows the level set amplitude.

Figures 8 ,

 8 9 and 10 show the simulation of 200 rigid spheres of radius R = 0.01 falling under gravity. The half thickness of the interface is ε = 2∆x. The coefficient of gravity g is set to -980. The repulsive coefficients are: k x = -g/10, k y = -g/10, k z = -g/10, k wall x = -g/40, k wall y = -g/40, k wall z

Figure 7 : 7 Figure 8 : 7 Figure 9 : 7 Figure 10 :

 77879710 Figure 7: Study of the grid convergence with a test case of 25 rigid disks falling under gravity. From left to right, grid 128 × 128, grid 256 × 256, grid 512 × 512 and grid 1024 × 1024. The background color shows the level set amplitude.

23 Figure 11 :Figure 12 : 5 Figure 13 : 5 Figure 14 :

 231112513514 Figure 11: Comparison of the two collision models for three disks of radius R = 0.1. The background colors show the level set amplitude.

8 Figure 15 : 8 Figure 16 :

 815816 Figure15: Simulation of 100 rigid particles submitted to gravity obtained with the FreeFem code[START_REF] Lefebvre | Fluid-particle simulations with freefem++[END_REF]

18 Figure 17 : 18 Figure 18 : 18 Figure 19 : 18 Figure 20 : 6 Figure 21 : 9 Figure 22 :Figure 23 :Figure 24 :

 18171818181918206219222324 Figure 17: Simulation of 30 square rigid bodies of different size falling under gravity. The intensity of the force depends on the relative velocities of the closest particles.

label 2 . 3 .• ∇φ n , 4 ., 5 . 1 =

 23451 Interpolation of the fluid velocity U n+1 on the finer grid providing U n+1 g Advection of the level set function φ n with the interpolated fluid velocity in a vicinity of the vesicles:φ n+1 = φ n -∆tU n+1g Transport of the stretching e with the interpolated fluid velocity:e n+1 = e n -∆tU n+1 g • ∇e ne n (n(φ n ) ⊗ n(φ n )) : D(U n+1 g )Update of the distance function ϕ n+1 |φ n+1 |, Update L n+1 0 ,L n+1 1 using φ n+1Perform the multi label fast marching method described in Chapter 3 6. Update the level set function as a signed distance function in thin narrow-bands around the vesicles 7. Update the N reconstituted level set functions

Figure 1 . 10 :

 110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :Figure 1 . 11 :

 26110111 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :

 26110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :

 26110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :Figure 1 . 11 :

 26110111 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :

 26110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :

 26110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :Figure 1 . 11 :

 26110111 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :

 26110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :

 26110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :Figure 1 . 11 :

 26110111 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :

 26110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :

 26110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :Figure 1 . 11 :

 26110111 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :

 26110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :

 26110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26 (Figure 1 . 10 :Figure 1 . 11 :

 26110111 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 10 :

 110 Figure 1.10: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Figure 1 . 11 :

 111 Figure 1.11: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

26

 26

Figure 25 :Figure 26 :

 2526 Figure 25: Simulation of two vesicles in a Poiseuille flow performed on three grid levels. From top to bottom, the associated discretization space steps are :h = 7.81 × 10 -3 , h 3.90 × 10 -3 , and h 1.95 × 10 -3 .

Figure 27 :

 27 Figure 27: Simulation of a vesicle passing through a bifurcation. The background color corresponds to the magnitude of the velocity. The white contour represents the isoline φ = ε.

Figure 28 :

 28 Figure 28: Simulation of 105 vesicles in a Poiseuille flow. The colors indicate the values of the label map L 0 from dark blue for the first body to dark orange for the 105 th body and red for the fluid that is the 106 th object.

Table 1 :

 1 Computational time of both algorithms.

	Number of Algorithm 1 Algorithm 2
	disks	CPU time	CPU time
	10	0.1	0.05
	50	0.5	0.6
	100	1.0	0.7
	200	2.1	1.4
	400	4.0	2.8

Table 2 :

 2 whereas with the present model the CPU time is only 25 times larger. So our method is around 340 times faster for 400 bodies. Averaged CPU time using the N level set decomposition

	Number of Collision model (6) Penalization model (19)	Total
	disks	CPU time		CPU time	CPU time
	2	0.02		0.06	0.2
	5	0.17		0.16	0.48
	10	0.72		0.35	1.24
	25	4.87		0.88	6
	50	19.25		1.75	21.5
	100	80.8		3.9	85.3
	400	1583.4		19.75	1605.3
	Number of Model (7) Model (23) Label redefinition	Total
	disks	CPU time CPU time	CPU time	CPU time
	2	0.015	0.05	0.008	0.2
	5	0.015	0.06	0.014	0.23
	10	0.015	0.09	0.02	0.25
	25	0.016	0.18	0.08	0.4
	50	0.016	0.3	0.16	0.6
	100	0.016	0.56	0.23	0.9
	400	0.016	2.52	2.06	4.7

Table 3 :

 3 Averaged CPU time using the algorithm of section 5.2.1

model will be compared numerically to a penalization model formulated in the context of a level set decomposition. The second part will be dedicated to numerical illustrations.

(g) t = 0.32 (h) t = 0.33 (i) t = 0.42

Acknowledgment

Most of the numerical simulations have been run on PLAFRIM platform supported by IMB University of Bordeaux and INRIA Bordeaux -Sud Ouest.