Gwenn Boedec
email: gwenn.boedec@univ-amu.fr

Marc Leonetti
email: leonetti@irphe.univ-mrs.fr

Marc Jaeger
email: marc.jaeger@centrale-marseille.fr

Isogeometric FEM-BEM simulations of drop, capsule and vesicle dynamics in Stokes flow

Keywords: Drops, Capsules, Vesicles, Stokes flow, Boundary element method, Finite element method

We develop an algorithm for the three dimensional simulation of the dynamics of soft objects (drops, capsules, vesicles) under creeping flow conditions. Loop elements are used to describe the shape of the soft objects. This surface representation is used both for membrane solver based on finite element method (FEM) and for the fluid solver based on the boundary element method (BEM). This isogeometric analysis of the low Reynolds fluid-structure interaction problem is then coupled to high-order explicit time stepping or second-order implicit time stepping algorithm. For vesicles simulation, a preconditioner is designed for the resolution of the surface velocity incompressibility constraint, which is treated by the use of a local Lagrange multiplier. A mesh quality preserving algorithm is introduced to improve the control mesh quality over long simulation times. We test the proposed algorithm on capsule and vesicle dynamics in various flows, and study its convergence properties, showing a second-order convergence O(N -2 ) with mesh number of elements.

Introduction

Dynamics of small soft objects in viscous flows have been the subject of numerous studies in the past decades. One fundamental motivation is to better understand the behavior of complex fluids, down to the microscopic level of their elementary constituents: polymers, flexible fibers, drops and bubbles are the most studied soft particles. Over the two last decades however, attention has started to shift towards more complex artificial particles such as capsules and vesicles, which can be thought of as drops enclosed by a complex interface. Vesicles are drops enclosed by an incompressible fluid lipid membrane which resists bending, while capsules are drops enclosed by a solid polymerized membrane resisting shear and area deformation (and bending).

The studies of these particles were first motivated by the desire to model blood flow down to individual cells that compose it: red blood cells (RBCs), platelets, white blood cells. The common model for RBCs from a fluid dynamics point of view is to consider it as a combination of a vesicle and a capsule: an elastic network of polymers (cytoskeleton) is anchored into a lipid membrane, giving the RBCs membrane the properties of surface incompressibility, shear elasticity, as well as bending resistance [START_REF] Li | Continuum-and particle-based modeling of shapes and dynamics of red blood cells in health and disease[END_REF][START_REF] Freund | Numerical simulation of flowing blood cells[END_REF]. Thus both capsules and vesicles makes excellent candidates as biomimetic models of RBCs, allowing to focus on the influence of specific mechanical properties.

The numerical simulation of the dynamics of these particles requires a fluid solver to describe internal and external flows, a membrane solver to describe interfacial forces, as well as a proper description and tracking of the interface position to couple those solvers. Several computational models have been developed for capsules, vesicles and RBCs, and we refer to recent reviews [START_REF] Petia | Vesicles and red blood cells in flow: From individual dynamics to rheology[END_REF][START_REF] Li | Continuum-and particle-based modeling of shapes and dynamics of red blood cells in health and disease[END_REF][START_REF] Freund | Numerical simulation of flowing blood cells[END_REF][START_REF] Barthès-Biesel | Motion and deformation of elastic capsules and vesicles in flow[END_REF] for a complete overview of the existing literature.

For low Reynolds number flows typical of these objects, the Boundary Element Method (BEM) simplifies greatly the coupling of fluid and membrane solvers since only the interface needs to be meshed, and thus the same interface representation/mesh can be used for both solvers. However, BEM usually uses C 0 meshes while a direct computation of membrane bending forces requires at least a C 4 representation of the position, since fourthorder derivatives of the position appears in the bending forces. To circumvent this difficulty, special methods can be designed to compute membrane bending forces on C 0 meshes, using either a local reconstruction of the interface [START_REF] Yazdani | Three-dimensional numerical simulation of vesicle dynamics using a front-tracking method[END_REF][START_REF] Farutin | 3d numerical simulations of vesicle and inextensible capsule dynamics[END_REF][START_REF] Mendez | An unstructured solver for simulations of deformable particles in flows at arbitrary reynolds numbers[END_REF], or designing computations of Laplace-Beltrami operator on C 0 based on discrete geometry [START_REF] Thierry Biben | Three-dimensional vesicles under shear flow: Numerical study of dynamics and phase diagram[END_REF][START_REF] Boedec | 3d vesicle dynamics simulations with a linearly triangulated surface[END_REF][START_REF] Trozzo | Axisymmetric boundary element method for vesicles in a capillary[END_REF]. On the other hand, methods based on a spectral representation have been developed [START_REF] Zhao | A spectral boundary integral method for flowing blood cells[END_REF][START_REF] Shravan | A fast algorithm for simulating vesicle flows in three dimensions[END_REF][START_REF] Zhao | The dynamics of a vesicle in simple shear flow[END_REF][START_REF] Rahimian | Boundary integral method for the flow of vesicles with viscosity contrast in three dimensions[END_REF], but their global basis function implies that the whole interface must be refined even if the deformation is highly localized. It would thus be desirable to design a method allowing for a local representation of the interface that could be used both in fluid and membrane solver. It is possible to lift the C 4 requirement for interface representation if the forces are computed in the weak sense, using a finite element method, which only requires the approximation of the position to be in H 2 . [START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF] have shown that subdivision surfaces meet this requirement and can be used to accurately model thin shells. Recently, schemes based on subdivision surfaces have been developed for capsules [START_REF] Le | Subdivision elements for large deformation of liquid capsules enclosed by thin shells[END_REF][START_REF] Le | A front-tracking method with catmullclark subdivision surfaces for studying liquid capsules enclosed by thin shells in shear flow[END_REF][START_REF] Huang | Three-dimensional simulation of elastic capsules in shear flow by the penalty immersed boundary method[END_REF] or vesicles [START_REF] Feng | Finite element modeling of lipid bilayer membranes[END_REF][START_REF] Ma | Viscous regularization and r-adaptive remeshing for finite element analysis of lipid membrane mechanics[END_REF][START_REF] Zhao | The dynamics of a vesicle in a wall-bound shear flow[END_REF][START_REF] Spann | Loop subdivision surface boundary integral method simulations of vesicles at low reduced volume ratio in shear and extensional flow[END_REF], but [START_REF] Le | Subdivision elements for large deformation of liquid capsules enclosed by thin shells[END_REF][START_REF] Le | A front-tracking method with catmullclark subdivision surfaces for studying liquid capsules enclosed by thin shells in shear flow[END_REF][START_REF] Huang | Three-dimensional simulation of elastic capsules in shear flow by the penalty immersed boundary method[END_REF] use immersed boundary method where the forces (and thus the interface representation) are spread to fluid solver nodes via approximate delta-function, while in [START_REF] Zhao | The dynamics of a vesicle in a wall-bound shear flow[END_REF][START_REF] Spann | Loop subdivision surface boundary integral method simulations of vesicles at low reduced volume ratio in shear and extensional flow[END_REF] bending forces are computed by first variation of the bending energy, computed from a surrogate surface based on Loop subdivision, but every other physical quantity is computed on linear basis, especially the boundary integral fluid solver. Finally, [START_REF] Feng | Finite element modeling of lipid bilayer membranes[END_REF][START_REF] Ma | Viscous regularization and r-adaptive remeshing for finite element analysis of lipid membrane mechanics[END_REF] computes equilibrium shape of vesicles but do not consider coupled dynamics under flow.

On the other hand, the rapid development of isogeometric analysis [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF] has shown that using a consistent geometry representation throughout the whole computation can greatly improve the quality of the results obtained. The idea stems from the fact that in engineering analysis CAD softwares used for designing pieces use a different surface representation like Bézier curves, splines, or Non-Uniform Rational Basis Spline (NURBS) than the software used for finite element analysis of these pieces (generally linear and quadratic Lagrange elements). Thus the standard finite element analysis starts with a modeling error in the sense that shape from CAD software cannot be represented exactly in the FEM model: before doing any computation, it is necessary to discretize the shape with a mesh which is typically based on Lagrange elements. In two dimensions, the simplest mesh would use linear triangle which will generally represent only imperfectly the real geometry of the computational domain (e.g. a disk-shape domain). Only in the limit in infinite refinement the mesh will converge to the true geometry. On the other hand, representing simple geometrical shape such as circle or ellipsoids is possible with NURBS, even with a finite number of elements. The idea behind isogeometric analysis is to use the same approximation for both the geometry discretization and the solution discretization. This leads to a significant improvement of solution quality [START_REF] Hughes | Isogeometric analysis: Cad, finite elements, nurbs, exact geometry and mesh refinement[END_REF]. While this approach has been used in a boundary element context for elastostatics, aeroacoustics, electrostatics, to our knowledge, its use for viscous flow simulations has been so far limited to [START_REF] Heltai | Nonsingular isogeometric boundary element method for stokes flows in 3d[END_REF] where no deformation of the interface is involved, and [START_REF] Joneidi | Isogeometric boundary integral analysis of drops and inextensible membranes in isoviscous flow[END_REF] for drops and inextensible membranes.

In this paper, we propose an algorithm based on the idea of using Loop subdivision surfaces to represent the shape of deformable objects under flow. Starting from this assumption, a natural question is the following: how to evolve dynamically the control values of Loop elements such that the resulting evolution of the (limit) physical surface they control is a discretized version of the fluid structure interaction problem. To do so, we derive a consistent framework to discretize the problem: the boundary element method extends naturally to Loop elements to compute viscous flows, while a finite element framework can be used for membrane forces computations with a variety of physical properties of the interface. Main originality of this work is that Loop elements are used to discretize every physical quantity: position, velocity, forces, tension field. This algorithm has been already used to simulate capsule deformation in extensional flow [START_REF] De Loubens | Stretching of capsules in an elongation flow, a route to constitutive law[END_REF] and drop dynamics in shear flow [START_REF] Gounley | Influence of surface viscosity on droplets in shear flow[END_REF] with surface viscosity effects considered. Here, we describe in greater details the algorithm which was briefly outlined in [START_REF] Gounley | Influence of surface viscosity on droplets in shear flow[END_REF], especially discretization of Loop boundary element as well as update of control values (including a mesh quality preserving algorithm). We also discuss extension of this algorithm to vesicle simulations which requires treatment of surface incompressibility, and introduce an implicit-time stepping algorithm. To accelerate the resolution of surface incompressibility constraint, a preconditioner based on a linear discretization is introduced. The paper is organized as follows : the system that we aim to simulate is described in section 2, the numerical description of this system is presented in section 3. The algorithm is then validated with comparison to literature results in section 4, together with a discussion of its properties.

Physical system description

Fluid flows

The internal flow (denoted with superscript int) and the external flow (ext) are described by the Stokes equation:

η int,ext ∆v int,ext -∇p int,ext = 0 ∇ • v int,ext = 0 ( 1
)
where η is the viscosity of the fluid, v the velocity and p the pressure.

Membrane mechanics

The membrane is a closed surface S, which is described by its position x(t) at time t, as well as its reference position x 0 (surface S 0 ) if any. The surface density of force exerted by the membrane f onto surrounding fluids is given by the first variation of its surface energy

f = f (x) = - 1 √ a δE δx , E = S w s dS (2)
where w s is the surface energy density which completely describes the mechanical properties of the interface, and a is the determinant of the local metric.

For example, a drop with surface tension γ can be modeled with a surface energy w s = γ, which upon application of equation ( 2) leads to surface density of forces

f = f γ (x) = ∇ s γ + 2γHn (3)
where ∇ s is the surface gradient operator, H is the mean curvature of the interface, and n is the outward pointing normal vector. For a vesicle membrane, the most used model for surface energy is Helfrich energy, together with a local Lagrange multiplier γ for local surface incompressibility constraint (∇ s • v = 0):

w s = γ + w H s , w H s = κ 2 (2H) 2 (4)
Note that the Helfrich model may contain a spontaneous curvature H 0 which is set to zero here. It also contains a Gaussian curvature term which can be ignored here since its integral is constant by Gauss-Bonnet theorem (no topological changes considered). For a capsule membrane, several modeling based on thin-shell equations [START_REF] Green | Theoretical Elasticity[END_REF][START_REF] Green | Large Elastic Deformations[END_REF][START_REF] Naghdi | The Theory of Shells and Plates[END_REF] have been developed [START_REF] Barths-Biesel | Motion of a spherical microcapsule freely suspended in a linear shear flow[END_REF][START_REF] Barths-Biesel | The time-dependent deformation of a capsule freely suspended in a linear shear flow[END_REF][START_REF] Pozrikidis | Effect of membrane bending stiffness on the deformation of capsules in simple shear flow[END_REF]. We refer to recent reviews [START_REF] Pozrikidis | Interfacial dynamics for stokes flow[END_REF][START_REF] Barthès-Biesel | Motion and deformation of elastic capsules and vesicles in flow[END_REF] for more details and simply note that the most commonly used model for numerical simulations of capsules is a membrane model, neglecting bending resistance, and modeling the interface as a two-dimensional elastic material. Within this framework, specific interfacial properties are described by the choice of an appropriate constitutive law [START_REF] Diaz | Effect of constitutive laws for two-dimensional membranes on flow-induced capsule deformation[END_REF] for the membrane, for which there is two common alternatives: either a strain-softening (Neo-Hookean law, noted NH) or strain-hardening (Skalak [START_REF] Skalak | Strain energy function of red blood cell membranes[END_REF] law, noted Sk). For these two laws the surface density of energy is defined on the reference configuration S 0 as:

w N H s 0 = G s 2 I 1 -1 + 1 I 2 + 1 , w Sk s 0 = G s 4 I 2 1 + 2I 1 -2I 2 + CI 2 2 (5)
with G s the surface shear modulus, C a parameter controlling the area dilatation resistance in Skalak law, and I 1 , I 2 strain invariants related to the principal extension ratios λ 1 , λ 2 by [START_REF] Skalak | Strain energy function of red blood cell membranes[END_REF] :

I 1 = λ 2 1 + λ 2 2 -2 I 2 = λ 2 1 λ 2 2 -1 (6)

Coupling conditions

Membrane and fluids surrounding it are coupled by the following conditions

• Quasi-static mechanical equilibrium: Membrane's inertia is negligible, thus, the surface density of forces due to membrane deformation is at equilibrium with fluid stresses, leading to the condition

[[σ]] • n + f = 0 (7)
where

[[σ]] = [σ ext -σ int ]
is the jump of fluid stresses at the interface and f is a density of force exerted by the membrane, see equation [START_REF] Barthès-Biesel | Motion and deformation of elastic capsules and vesicles in flow[END_REF].

• Continuity of the velocities at the interface:

v int (x) = v ext (x) = v S ∀x ∈ S (8)
where v S is the velocity of fluids at the interface.

• Advection of the membrane by the flow: The membrane is assumed to be impermeable, thus the velocity of the membrane is exactly the velocity of the fluid at the interface

dx dt = v S ∀x ∈ S (9)
• Surface incompressibility constraint (vesicles): The velocity field of a vesicle membrane is required to satisfy a surface divergence-free constraint:

∇ s • v S = 0 (10) 6 Control mesh 1 iteration of subdivision Limit surface (a) (b) (c)

Numerical discretization

Surface representation : Subdivision surfaces

The surface is discretized using Loop subdivision elements [START_REF] Teorell | Smooth subdivision surfaces based on triangles[END_REF], which are assembly of linear triangles refined using a subdivision process. Starting from an initial coarse mesh, called control mesh, new elements are created and vertices positions are updated according to Loop subdivision rules [START_REF] Teorell | Smooth subdivision surfaces based on triangles[END_REF]. In the limit of infinite subdivisions, the mesh converges toward a limit surface, which is C 2 continuous, except at some extraordinary vertices where it is C 1 . Because of the regularity of the limit surface, numerical computation of geometrical quantities such as normal vector or curvature tensors is more straightforward and accurate that when using classical linear or quadratic Lagrange elements which are only C 0 . The ordinary or regular vertices are linked to six elements, while any vertex linked to a different number of element is termed "extraordinary" since the subdivision rule must be adapted here. The number of extraordinary vertices depends on the initial mesh (level 0 of subdivision). Starting from an icosahedron, there are 12 vertices linked to only 5 elements. This number does not change with refinement of mesh, since any vertex introduced by Loop subdivision is regular. Note also that it is not possible to discretize a spherical surface without introducing irregular vertices : the reason behind this is that the Euler characteristic of a sphere is constant and can related to discretization parameters (number of faces, vertices, edges) by Euler's formula. An example of a simulation of a capsule under shear flow using such elements is presented in figure 1.

An interesting property of Loop subdivision surface is that it is possible to analytically compute the limit position of the original vertices composing the control mesh. Moreover, since the seminal work of [START_REF] Stam | Evaluation of loop subdivision surfaces[END_REF], it is also possible to compute the limit position of any point inside the initial triangles by using shape functions N i :

x e (s 1 , s 2 ) = p∈one-ring N p (s 1 , s 2 )X p (11)
where (s 1 , s 2 ) is a local parametrization of the element. The shape functions are polynomials of the local parametrization (see [START_REF] Stam | Evaluation of loop subdivision surfaces[END_REF]), and X p are the nodal values of the position field. Note that these values have no real physical meaning : they are neither the real position of the vertex, nor the limit position of this vertex. They are just numerical values controlling the approximation of the shape. Contrary to classical Lagrange triangular elements, the support of the shape function for a given node p extends outside the triangles containing p as a vertex, up to every triangle which is linked to p by at most an edge (see figure 2 (b)). Thus, the position inside a given triangle depends not only on nodal values of the vertices of this triangle, but also of nodal values X p of the nodes related to this triangle by at most one edge. These elements are called the one-ring of element e.

We refer to Stam [START_REF] Stam | Evaluation of loop subdivision surfaces[END_REF] and Cirak et al. [START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF] for a detailed discussion of the evaluation of the shape functions inside an element. In this work, this local basis (N 1 , N 2 , ...) is used to represent every function f which is defined on the interface, e.g. membrane force density or membrane velocity. Thus, f is written as:

f (x) = f (x(s 1 , s 2 )) = f (s 1 , s 2 ) = p∈one-ring F p N p (s 1 , s 2 ) ( 12
)
where F p is the p-th nodal value of the field representing the function f .The series on the right-hand side of equation ( 12) is an approximation of the function f using Loop subdivision functions. While equation [START_REF] Dahl | Experimental observation of the asymmetric instability of intermediate-reduced-volume vesicles in extensional flow[END_REF] gives a way to compute f at any point of the membrane if the nodal values F p are known, one also needs to convert a field defined onto the interface f (x) into its nodal values F p , that is, given f find the approximation under the form [START_REF] Dahl | Experimental observation of the asymmetric instability of intermediate-reduced-volume vesicles in extensional flow[END_REF] such that the approximation error is minimized. Two main alternatives could be considered : Galerkin formulation in which the approximation error is measured in the L 2 -norm over the surface, or a collocation formulation in which the approximation error is measured pointwise. In this work, we choose to collocate the known field at vertices. This choice will be discussed in section 3.2 when considering velocity evaluation. Note that the two formulations were compared in the context of drop isogeometric simulations by [START_REF] Joneidi | Isogeometric boundary integral analysis of drops and inextensible membranes in isoviscous flow[END_REF], which found that while the Galerkin formulation has better convergence properties, collocation is both more stable and more fast. Thus, assuming f is known at each point of the interface, we search the nodal values F p of a field f such that:

f p = f (x = x p ) = f (x p ) = q∈one-ring F q N q (s 1 (x p ), s 2 (x p )) ∀p ∈ {1, ..., N nodes } (13)
where x p is the position of the p-th control vertex on the limit surface, and N q (x p ) are the shape functions evaluated at parameters corresponding to x p . Assembling all the nodal values of the field f into a global vector {F p }, equation ( 13) is rewritten in matrix form as:

{f p } = f (x p ) = C {F p } ( 14
)
where C is a collocation matrix in the sense that it enforces pointwise the relationship between the known limit values f p and the unknown nodal values F p of the field f which interpolates f at positions x p . This matrix is sparse and thus the linear system is solved with the sparse solver UMFPACK [START_REF] Davis | Algorithm 832: Umfpack, an unsymmetric-pattern multifrontal method[END_REF][START_REF] Davis | Direct Methods for Sparse Linear Systems[END_REF]. This linear system is solved each time a conversion from limit value to nodal values is needed. However, for system sizes considered here, the computational cost is always a small fraction (typically, less than 10%) of the total cost of one time step. In the following, the notation {f } will be used to designate the N-dimensional vector containing the values of f evaluated at control vertices limit positions {x} while the notation {F } will designate the N-dimensional vector containing the nodal values defined on the control mesh.

Fluid solver : Boundary element method

For a soft liquid particle suspended in a viscous fluid, the velocity at a point y of the particle interface S is given by the integral relation [START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF]:

v(y) = 2 1 + λ v ∞ (y) + 1 4πη ext (1 + λ) S G ij (x, y)f i (x)dS(x)e j + (1 -λ) 4π(1 + λ) S v i (x)T ijk (x, y)n k (x)dS(x)e j (15)
where v ∞ is the imposed far-field velocity, G, T are the free space Stokeslet and Stresslet, f is the surface density of force due to the membrane, e i is the unit i-th Cartesian vector and λ = η int /η ext is the viscosity ratio. Roman indices (i, j, k, l) can vary between 1 and 3 and denote the Cartesian components of vectors and tensors. To lighten the expressions, Einstein summation convention is used e.g.

G ij f i = G 1j f 1 + G 2j f 2 + G 3j f 3 .
The integrals in [START_REF] De Loubens | Stretching of capsules in an elongation flow, a route to constitutive law[END_REF] needs special numerical treatment to deal with singularities of G, T occurring when x → y. A powerful method to treat such integrals consists in subtracting the singularities by using analytical integral identities [START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF]. This has been largely used for Stresslet as well as Stokeslet when the membrane force was purely normal (as is the case for a clean drop). However, more general expression for the forces (with tangential components) required different treatment like singularity cancellation by special quadratures. Only recently was derived a new integral identity allowing to treat general expression of membranes forces [START_REF] Farutin | 3d numerical simulations of vesicle and inextensible capsule dynamics[END_REF] without resorting to special integration rule for singular elements. Using both singularity subtraction for Stresslet [START_REF] Pozrikidis | Boundary integral and singularity methods for linearized viscous flow[END_REF] and Stokeslet [START_REF] Farutin | 3d numerical simulations of vesicle and inextensible capsule dynamics[END_REF], the boundary integral relation is rewritten as:

v(y) = v ∞ (y) + 1 8πη ext S G ij (y, x) fj (y, x)dS(x)e i + 1 4π S R(y, x)dS(x) ∧ [n(y) ∧ f (y)] + (1 -λ) 8π S [v i (x) -v i (y)]T ijk (x, y)n k (x)dS(x)e i ( 16
)
where the expression for the modified force density f and R are [START_REF] Farutin | 3d numerical simulations of vesicle and inextensible capsule dynamics[END_REF]:

f (y, x) = f (x) + [-n(x)(f (y) • n(y)) + n(x) ∧ (n(y) ∧ f (y))] R(y, x) = [(y -x) • n(x)] |y -x| 3 (y -x) (17)
In the expression ( 16) every integrands is now non singular and can thus be integrated by standard quadratures. In this work, Gauss points are used. Note that while there is an additional integral in equation ( 17) and thus an additional computational cost when compared to classical formulation, this doesn't change the overall O(N 2 ) computational complexity of the algorithm, only the prefactor. On the other hand, since all the integrals can be treated with standard quadratures without distinction between singular and regular elements, it eases the parallelization of the computation, when compared to using special quadratures on singular elements (e.g. Duffy quadrature or introducing polar coordinates).

If the internal and external viscosities are equal (λ = 1), the integral relation ( 16) is simplified and directly gives the velocity if the membrane forces are known. If λ = 1, the integral relation [START_REF] De Loubens | Tank-treading of microcapsules in shear flow[END_REF] needs to be solved for the velocity. This equation is enforced at each node of the mesh. Equation ( 16) is a linear integral relation between physical velocities v and forces f , which depends linearly on the nodal values of velocities V and forces F via interpolation [START_REF] Dahl | Experimental observation of the asymmetric instability of intermediate-reduced-volume vesicles in extensional flow[END_REF]. Equation ( 16) is thus rewritten as :

I + (λ -1)TC -1 {v} = {v ∞ } + G{F } ( 18
)
where I is the identity matrix, G, T are the matrix corresponding to the discretization of ( 16) and {v}, {v ∞ } are vectors containing the limit values of velocities, and {F } contains the nodal values of forces. This linear system is solved by GMRES [START_REF] Saad | Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] without explicitly computing the G, T matrices, but only the matrix-vector products. These matrix-vectors products corresponds to linear integral operators e.g.:

S G ij (y, x) fj (y, x)dS(x) = N el e=1 S e G ij (y, x) fj (y, x)dS e (x) (19)
where S e is the surface of element e. Integrals are then transformed in the parameter space (ξ, η) and then performed by numerical quadrature, using Gauss-Hammer quadrature points on triangular elements:

S e G ij (y, x) fj (y, x)dS e (x) = 1 s 1 =0 1-s 1 s 2 =0 G ij (y, x) fj (y, x) √ ads 1 ds 2 = N quad k=1 w k G ij (y, x k ) fj (y, x k ) a(s 1 k , s 2 k ) (20)
with

x k = x(s 1 k , s 2 k
) the position of k-th quadrature point. In this work, 12 Gauss quadrature points per element are used for the computation of boundary element integrals.

After obtaining of limit values of velocities, the resulting values of velocities on the limit surface are converted into nodal values of the velocity by resolution of the collocation equation [START_REF] Davis | Direct Methods for Sparse Linear Systems[END_REF]. Note that a Galerkin formulation of equation ( 16) would require an additional integration over the surface, meaning that the boundary integral equation would need to be evaluated at each quadrature point instead of mesh vertices. Since there is always more elements than vertices in the mesh, even using a one-point quadrature rule would require more integral evaluations than in the collocation case, and the gap only widens if one uses higher-order quadrature rules. Thus, in terms of computational complexity, it seems advantageous to use a collocation formulation.

Membrane solver : Finite element method

Starting from a reference shape x 0 (s 1 , s 2 ), the membrane is deformed into its current shape x(s 1 , s 2 ). To describe elastic deformations, it is useful to introduce a local basis on the deformed configuration composed of the tangent vectors a 1 , a 2 defined by the parametrization and the outward pointing normal vector n :

a 1 = x ,1 = ∂x ∂s 1 , a 2 = x ,2 = ∂x ∂s 2 , n = a 1 ∧ a 2 |a 1 ∧ a 2 | (21)
The local basis on reference configuration (a 0 1 , a 0 2 , n 0 ) is defined similarly, with x 0 replacing x. In the following, Greek indices (α, β, γ, δ) will be used to denote components of vectors and tensors in the local basis, and can take the value 1 or 2. With this notation, the deformed configuration metric a αβ and curvature b αβ tensors are defined by :

a αβ = a α • a β , b αβ = a α,β • n = ∂a α ∂s β • n (22)
with a similar definition for reference configuration metric a 0 αβ and curvature b 0 αβ . The determinant of the metric is noted a = det(a αβ ) and is related to the elementary area by dS = √ ads 1 ds 2 . Finally, it is useful to introduce the inverse metric a αβ defined by a αγ a γβ = δ α β with δ α β the Kronecker symbol. Following [START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF], we write the principle of virtual work for the membrane as :

- S σ αβ δ(E αβ ) + µ αβ δ(B αβ ) dS + S f ext • δxdS = 0 (23)
The first integral is the internal virtual work for a virtual displacement δx, where σ αβ , µ αβ are the effective membrane and bending stresses. The second integral is the external virtual work. E αβ is the Green-Lagrange strain tensor measuring membrane deformations

E αβ = 1 2 a α • a β -a 0 α • a 0 β = 1 2 a αβ -a 0 αβ ( 24
)
and B αβ is a measure of bending strains defined as

B αβ = a α,β • n -a 0 α,β • n 0 = b αβ -b 0 αβ (25) Noting that δ(E αβ ) = 1 2 δ(a αβ ), δ(B αβ ) = δ(b αβ ) and f ext = -f , equation (23) writes: S 1 2 σ αβ δ(a αβ ) + µ αβ δ(b αβ ) dS + S f • δxdS = 0 ( 26
)
The above equality holds whether it is a drop, vesicle or capsule (with Kirchoff-Love hypothesis), and thus constitutes an interesting unified framework to study both types of objects. Before describing the specific interfacial properties that comes into play via membrane (σ αβ ) and bending stresses (µ αβ ), we describe how this equality is discretized. From ( 21) and [START_REF] Feng | Finite element modeling of lipid bilayer membranes[END_REF] we have

δ(a αβ ) = [δa α • a β + a α • δa β ] = [δx ,α • x ,β + x ,α • δx ,β ] δ(b αβ ) = n • δ(a α,β ) + a α,β • δ(n) = n • δ(a α,β ) + b αβ n + Γ γ αβ a γ • δ(n) = n • δ(a α,β ) -Γ γ αβ n • δ(a γ ) (27)
where Γ γ αβ = a γδ a α,β • a δ are the Christoffel symbols. We used the identities

δ(||n|| 2 ) = 2n • δn = 0 and δ(a γ • n) = δ(a γ ) • n + a γ • δ(n) = 0 to express δ(b αβ ) as a function of δx.
Using the same basis functions (Loop elements) for the Cartesian components of membrane force, position and virtual displacement, we write:

f i (s 1 , s 2 ) = p∈one-ring F p i N p (s 1 , s 2 ) x i (s 1 , s 2 ) = p∈one-ring X p i N p (s 1 , s 2 ) δx i (s 1 , s 2 ) = p∈one-ring δX p i N p (s 1 , s 2 ) ( 28
)
The principle of virtual work [START_REF] Green | Large Elastic Deformations[END_REF] is thus discretized as:

N el e=1 S e [I σ + I M + I f ] √ ads 1 ds 2 = 0 (29)
with

I σ = p∈one-ring 1 2 σ αβ (N p,α x i,β + N p,β x i,α ) δX p i I M = p∈one-ring µ αβ n i N p,αβ -Γ γ αβ n i N p,γ δX p i I f = p∈one-ring q∈one-ring N p N q F q i δX p i ( 30
)
where N p,α = ∂Np ∂s α . Note that a double sum appears explicitly in I f because nodal values F are unknowns, while a double sum is "hidden" in I σ , I M because nodal values of positions are known and thus can be used to evaluate every terms inside the triangle (see discussion below).

Equation ( 29) can then be written in matrix vector form as

{RHS} + M{F} = 0 (31)
The right hand side vector {RHS} and the mass matrix M are formed by numerical integration of equation ( 29) using Gauss quadrature points, for instance, for membrane terms :

S e I σ √ ads 1 ds 2 ≈ N quad k=1 w k I σ (s 1 k , s 2 k ) a(s 1 k , s 2 k ) ( 32
)
where w k and (s 1 k , s 2 k ) are weights and parameters values of the k-th quadrature point, N quad is the total number of quadrature points, and the dependency of I σ and √ a on parameters has been made explicit for clarity. To evaluate I σ and √ a at parameters values (s 1 k , s 2 k ) one needs to evaluate shape functions and their derivatives, membrane stresses and moments, and also tangent and normal vectors. In this work, 12 Gauss quadrature points are used.

The evaluation of shape functions and their derivatives inside regular and irregular elements follows the strategy of [START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF]. The evaluation of the membrane stresses and moments and of virtual variations of membrane and bending strains is discussed in the following.

For drops and zero-thickness capsule model (for example, Neo-Hookean or Skalak law), the surface energy does not depend on curvatures, thus

µ αβ = 0 (33)
It then remains to compute σ, which for drops takes the form

σ αβ = γa αβ (34)
while for capsules it takes the form [START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a stokes flow[END_REF]:

σ αβ = 2 J s ∂w s ∂I 1 a 0,αβ + 2J s ∂w s ∂I 2 a αβ (35)
with J s = √ a √ a 0 the Jacobian of the transformation from reference to deformed position. Derivatives of the energy w s with respect to invariants is done analytically, and the invariants of the deformation are computed as:

I 1 = a αβ a 0,αβ -2 I 2 = a/a 0 -1 (36)
For a membrane obeying the Helfrich energy (4), the membrane and bending stresses write [START_REF] Boedec | Mechanics of a helfrich membrane as a two-dimensional cosserat medium[END_REF]:

σ αβ = 2 √ a ∂( √ aw H s ) ∂a αβ = κ 2 4H 2 a αβ -8Hb αβ + γa αβ µ αβ = ∂w H s ∂b αβ = κ 2 4Ha αβ (37)
In all cases, one needs to compute inverse metrics and curvature tensors and local coefficients (γ, H, J s , ∂ws ∂I 1 , ∂ws ∂I 2 ) at quadrature points. The metric tensors in the reference state a 0 αβ and the deformed state a αβ are readily computed from the interpolation of the position, e.g.:

a αβ (s 1 , s 2 ) = x ,α • x ,β = p∈one-ring q∈one-ring X p i X q i N p,α N q,β (38)
The inverses tensors (a 0,αβ , a αβ ) are computed explicitly at each quadrature point.

Since the Loop subdivision scheme guarantees C 2 continuity everywhere except at extraordinary vertices, it is possible to evaluate the curvature tensor at any quadrature point by direct differentiation of the shape functions, to yield:

b αβ (s 1 , s 2 ) = x ,αβ • n = p X p i N p,αβ e i • n (39)
where the normal vector is computed by

n(s 1 , s 2 ) = a 1 ∧ a 2 |a 1 ∧ a 2 | = p X p i ∂Np ∂s 1 e i ∧ q X q j ∂Nq ∂s 2 e j |a 1 ∧ a 2 | (40)
For drops and vesicles, the local value of the tension (or Lagrange multiplier) is computed from the interpolation of the tension field:

γ(s 1 , s 2 ) = p∈one-ring Γ p N p (s 1 , s 2 ) ( 41
)
For a drop, the surface tension is a physical property of the interface, which may depends on other variables such as temperature or surfactant concentration. For a vesicle, the tension is a Lagrange multiplier ensuring the constraint of surface divergence free velocity field, and is thus an additional unknown of the system. The determination of tension in this case is discussed in the following section.

Surface incompressibility constraint solver

For a vesicle, the surface incompressibility constraint is satisfied by the use of a local Lagrange multiplier γ, which is the solution of :

∇ s • v γ (x) + ∇ s • v ∞ (x) + ∇ s • v κ (x) = 0 ( 42
)
where v γ is the velocity due solely to tension, and v ∞ , v κ are the imposed external velocity and bending induced velocity. This equation is actually linear in the unknown tension field γ, and is thus rewritten as:

D γ {Γ} = -{∇ s • v ∞ (x) + ∇ s • v κ (x)} (43)
where D γ is a linear operator transforming the nodal values of tension {Γ} to the limit value of the surface divergence of the velocity field due to tension. Equation ( 43) is thus a collocation at nodes of surface incompressibility constraint [START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF]. Equation ( 43) is solved by GMRES [START_REF] Saad | Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. The operator D γ is thus never build explicitly, only its action on {Γ} is computed as follows. Given a particular vector {Γ}, the nodal values of the forces due to tension are solutions of the system [START_REF] Joneidi | Isogeometric boundary integral analysis of drops and inextensible membranes in isoviscous flow[END_REF], with the right-hand side computed with equation [START_REF] Le | A front-tracking method with catmullclark subdivision surfaces for studying liquid capsules enclosed by thin shells in shear flow[END_REF]. The limit values of velocity field {v} are then obtained by application of the BEM equation ( 18) without v ∞ . After resolution of the system [START_REF] Dodson | Dynamics of strain-hardening and strainsoftening capsules in strong planar extensional flows via an interfacial spectral boundary element algorithm for elastic membranes[END_REF], the nodal values {V } of the velocity field are obtained by application of the collocation equation [START_REF] Davis | Direct Methods for Sparse Linear Systems[END_REF]. Finally, the limit value of surface divergence is obtained by using the following formula

∇ s • v = a αβ ∂x i ∂s α ∂v i ∂s β ( 44
)
This equation is computed at each node. For regular vertices, the application is straightforward by injection of the interpolation of position and velocities in [START_REF] Pozrikidis | Effect of membrane bending stiffness on the deformation of capsules in simple shear flow[END_REF]. For irregular vertices, first derivatives are computed by using the tangent masks of Loop subdivision scheme [START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF]. The inverse metric is obtained by analytic inversion of metric.

Note that equation ( 43) is however badly conditioned and leads to an increasing number of GMRES iteration while decreasing the mesh size, as can be seen in table 1. To overcome this problem, a preconditioner has been designed. The idea is to use a coarse approximation of the operator D γ by approximating tension forces, velocity and surface divergence computation onto a C 0 mesh composed of linear triangles. The construction of these operators is based on our previous work using linear triangles for vesicle simulations [START_REF] Boedec | 3d vesicle dynamics simulations with a linearly triangulated surface[END_REF] and is briefly recalled here for completeness. The basic idea is to interpolate linearly fields such as tension and velocity on the triangles describing the shape

γ(s 1 , s 2 ) = γ 0 (1 -s 1 -s 2 ) + γ 1 s 1 + γ 2 s 2 (45)
where (s 1 , s 2 ) are local element coordinates and γ i is the nodal value at node i. Then, it is possible to compute gradients for tension or divergence for velocity with respect to local element coordinates:

∇ s γ = ∂γ ∂s α a α = (γ 1 -γ 0 )a 1 + (γ 2 -γ 0 )a 2 (46)
where a α = a αβ a β . Equation ( 46) gives a constant value by element, which is then converted into a nodal value by a weighted-area average, and then projected on tangent plane:

(∇ s γ) n = [I -n ⊗ n] • 1 A e∈E n A e (γ 1 -γ 0 )a 1 + (γ 2 -γ 0 )a 2 (47)
Similar construction is applied for the computation of surface divergence of the velocity field, and the normal component of the tension force is obtained by computing mean curvature at nodes using a discrete surface divergence theorem (see [START_REF] Boedec | 3d vesicle dynamics simulations with a linearly triangulated surface[END_REF] for details). It is then straightforward to construct matrix form of the relations between on the one hand limit values of tension ({γ}) and limit values of tension forces ({f γ }), and on the other hand, between limit values of velocities ({v}) and limit values of surface divergence ({∇ s • v}). To finish the approximation of D γ , the linear relation between limit values of forces ({f }) and limit values of velocities ({v}) expressed by [START_REF] De Loubens | Tank-treading of microcapsules in shear flow[END_REF] without viscosity contrast (λ = 1) is discretized over the linear mesh as in [START_REF] Boedec | 3d vesicle dynamics simulations with a linearly triangulated surface[END_REF] by using Gauss-Hammer quadrature for non-singular elements and Gauss-Legendre quadrature with polar coordinates for singular elements. Thus, with the matrices F γ , D, G defined by

{f γ } = F γ {γ} , {∇ s • v} = D{v} , {v} = G{f } (48)
No preconditioner the preconditioner D γ is defined as

preconditioner N el λ = 1 Ca = 1 λ = 1 Ca = 1 Ca λ Ca λ 0.
D γ = DC -1 GF γ (49)
and we solve:

D γ D γ -1 {y} = -{∇ s • v ∞ + ∇ s • v κ } , D γ {Γ} = {y} (50)
with {y} an intermediate variable of calculus. In theory, the C -1 is not necessary since for preconditioning, every variable is interpolated on the linear mesh. However, numerical tests have shown that including this matrix in the definition of D γ lowered the number of iteration. While we lack a rigorous explanation for this, we propose the following heuristic argument: Note that on a linear mesh, the value of the surface divergence of velocity depends only on the values of velocity at adjacent nodes. This is not the case with Loop subdivision, where the value of surface divergence on one node depends on the values of velocity at each node separated by at most two edges. There is thus an important loss of information in the C 0 approximation of F γ , D, G, which is balanced by the inclusion of the inverse of the collocation matrix (C -1 ). Note that the inverse of collocation matrix (C -1 ) depends only on mesh connectivity and not on the real shape of the mesh. Thus, it has to be computed only once at the beginning of the simulation. This preconditioner is computed once at the beginning of each time step with the current position x n of the interface, an then an LU-factorization is performed using LAPACK [START_REF] Anderson | LAPACK Users' Guide[END_REF]. This factorization is then used for all the intermediate steps if Runge-Kutta scheme is used or for all the Newton iterations if trapezoidal scheme is used. This renders simulations of vesicles much costlier than drop or capsule simulations, since precondition factorization can represent as much as 50% of the computational cost of a time step. Also, viscosity contrast was not included in the G matrix to avoid the need to solve an additional linear system. However, as shown in table 1, the preconditioning is still efficient even with large viscosity constrast. For the coupled system given by ( 58), the preconditioner is simply

P = I 0 0 DC -1 GF γ (51)

Time stepping algorithm

Two solvers are used for integrating the time evolution of the interface given by equation ( 9) : an explicit high-order solver and an implicit solver. Note that, for a given position of the interface x n , membrane forces can be computed directly from equation ( 31) and then the velocity can be computed by equation ( 18), possibly with a projection stage described in section 3.4 in order to satisfy the surface incompressibility constraint. We note by v(x n ) the velocity computed at position x n .

The explicit solver is a Runge-Kutta Fehlberg fourth-fifth order solver [START_REF] Fehlberg | Low-order classical runge-kutta formulas with step size control and their application to some heat transfer problems[END_REF], whose advantages are high-order and adaptive time stepping, allowing for very good preservation of invariants such as enclosed volume of fluid. Given a position and tension (x n , γ n ) at time t n , the position and tension (x n+1 , γ n+1 ) at time t n+1 are such that:

x n+1 = x n + ∆t 4 i=0 c i v (i) n (52)
where the intermediate velocities are computed by:

v (0) n = v(x n , γ n ) , ∇ s • v (0) n = 0 v (k) n = v x (k) n , γ (k) n , ∇ s • v (k) n = 0 with x (k) n = x n + ∆t k-1 i=0 β ik v (i) n ( 53
)
and c i , β ik the coefficients of the Runge-Kutta Fehlberg method [START_REF] Fehlberg | Low-order classical runge-kutta formulas with step size control and their application to some heat transfer problems[END_REF]. The intermediate velocities v n and using this intermediate position to compute all integrals needed to build systems given by ( 31) and [START_REF] Dodson | Dynamics of strain-hardening and strainsoftening capsules in strong planar extensional flows via an interfacial spectral boundary element algorithm for elastic membranes[END_REF]. Note that if this scheme is used for vesicles, the surface incompressibility constraint is enforced at each Runge-Kutta stage (see section 3.4 for details). A fifth order combination is also used to form another solution

xn+1 xn+1 = xn + ∆t 5 i=0 ĉi v (i) n (54)
which permits to estimate numerical error and thus select adaptively the time step.

For capsules without bending, the stability condition ∆t ≤ O(∆x/G s ) allows reasonable time steps to be selected, while for vesicles, the stability condition ∆t ≤ O(∆x 3 /κ) prevents the use of such an explicit scheme in practical conditions.

The implicit solver is the trapezoidal time scheme. Given a position and tension (x n , γ n ) at time t n , the position and tension (x n+1 , γ n+1 ) at time t n+1 are such that:

x n+1 = x n + ∆t 2 [v(x n , γ n ) + v(x n+1 , γ n+1 )] ∇ s • v(x n+1 , γ n+1 )) = 0 (55)
To solve these implicit equations, we proceed by iteration, and let

x n+1 (0) = x n + ∆t(v(x n , γ n )) γ (0) n+1 = γ n (56)
For the first time step, the tension γ 0 is found by solving ∇ s • v(x 0 , γ 0 ) = 0. Now, given a value of (x n+1 (r) , γ

n+1 ) at the r-th iteration, the residual writes:

res x = x (r) n+1 -x n + ∆t 2 (v(x n , γ n ) + v(x (r) n+1 , γ (r) n+1 )
)

res γ = ∇ s • v(x (r) n+1 , γ (r) n+1 ) (57)
If both residuals are less than a prescribed tolerance, the iteration finishes.

If not, we seek a correction (δx, δγ) such that (x

(r) n+1 +δx, γ (r)
n+1 +δγ) satisfies the linearization of (55): where J x is the Jacobian with respect to x, and v γ is the velocity due solely to tension, that is, the velocity given by equation ( 16) with f = f γ and v ∞ = 0. Note that tension-induced velocity is a linear operator in γ and thus v γ (x, δγ) is exactly the Jacobian with respect to γ applied to δγ , that is J γ δγ = v γ (x, δγ). Equation ( 58) is solved with the Jacobian-free Newton-Krylov method [START_REF] Knoll | Jacobian-free newtonkrylov methods: a survey of approaches and applications[END_REF]. Thus, in the GMRES solver, the Jacobian is not constructed explicitly, but only evaluated numerically through its action on vectors :

I - ∆t 2 J x δx - ∆t 2 v γ (x (r) n+1 , δγ) = -res x ∇ s • (J x δx) + ∇ s • v γ (x (r) n+1 , δγ) = -res γ (58)
J x δx ≈ v(x (r) n+1 + δx, γ n+1 ) -v(x (r) n+1 , γ n+1 ) ( 59
)
with a small parameter [START_REF] Knoll | Jacobian-free newtonkrylov methods: a survey of approaches and applications[END_REF]. Again, intermediate velocities are computed by moving the interface position to e.g. x (r) n+1 + δx and using this intermediate position to compute all integrals for [START_REF] Joneidi | Isogeometric boundary integral analysis of drops and inextensible membranes in isoviscous flow[END_REF] and [START_REF] Dodson | Dynamics of strain-hardening and strainsoftening capsules in strong planar extensional flows via an interfacial spectral boundary element algorithm for elastic membranes[END_REF].

The average number of GMRES iteration required to solve [START_REF] Zhao | A spectral boundary integral method for flowing blood cells[END_REF] and the average number of Newton iteration required to solve [START_REF] Yazdani | Three-dimensional numerical simulation of vesicle dynamics using a front-tracking method[END_REF] are reported in table 2. It can be noticed that, due to the ill-conditioning of J x , keeping the time step fixed while increasing the mesh resolution increases the number of iterations required. A possible remedy would be to design a preconditioner for the J x operator, as has been done in [START_REF] Shravan | A fast algorithm for simulating vesicle flows in three dimensions[END_REF] for vesicle simulations with spectral methods. Note that both the explicit and implicit solver are described for vesicle, but can be applied for capsules or drops, in which case the equations for tension γ are not solved.

Mesh quality preserving algorithm

For capsules, the shear elasticity prevents severe mesh distortion while for drops or vesicles, the membrane fluidity allows velocity fields that may strongly distort the mesh, as can be seen in figure 3 (c). Interestingly, even with badly shaped elements, the simulation is still stable for some time, and gives acceptable results. The most noticeable effect of poor mesh quality is an order of magnitude increase in the enclosed volume variation as shown in figure 3 (b), and increased variations in the Taylor parameter and inclination angle as shown in figure 3 (a), for a clean drop submitted to simple shear flow (v ∞ = γye x ). Note that the oscillation amplitude in the Taylor parameter is less than 0.5% and thus that the global picture is still correct, even without any treatment on the control mesh. One may argue that Taylor deformation parameter is a global quantity which thus hides the local details. To check the impact of the control mesh quality on local quantities we take four cross-section of the shape at z = 0; 0.5; 0.75; 0.85 as shown in figure 4 (a) and compare results for a control mesh without any treatment (left part of the shape) and a control mesh with a mesh quality preservation algorithm (described in the following) applied (right part of the shape). One can first notice that the cross-sections almost overlap, showing that the local limit shapes are very similar. Then, we measure the membrane force norm (thus essentially the curvature) along these cross sections and plot the results as a function of the polar angle in figure 4 (b). Again, results are almost identical, with the most noticeable difference being that the shape without mesh treatment overestimates the curvature at the tip of the drop. This example illustrates that control mesh quality seems not too critical to obtain acceptable results, and that Loop elements are certainly more robust than linear or quadratic Lagrange elements. A more systematic study would be needed to investigate this point further. Moreover, it is still worthwhile to preserve mesh quality, as it improves volume conservation, reduces parasitic oscillations of the Taylor parameter, and more importantly, allows to keep a reasonable time step : indeed, if mesh elements become to small, the stable time step will decrease dramatically, leading to increased computational cost.

A mesh quality preserving algorithm has thus been designed to redistribute mesh nodes onto the interface.

This algorithm is similar in principle to existing algorithms such as [START_REF] Loewenberg | Numerical simulation of a concentrated emulsion in shear flow[END_REF][START_REF] Zinchenko | A novel boundary-integral algorithm for viscous interaction of deformable drops[END_REF][START_REF] Cristini | An adaptive mesh algorithm for evolving surfaces: Simulations of drop breakup and coalescence[END_REF][START_REF] Wang | Elastic mesh technique for 3d {BIM} simulation with an application to underwater explosion bubble dynamics[END_REF][START_REF] Zinchenko | Emulsion flow through a packed bed with multiple drop breakup[END_REF] : using the fact that, for vesicles or drops, there is no tangential reference configuration means that the tangent velocities can be described in an Eulerian manner, that is, decoupled from mesh velocities. In other words, it is possible to move grid points along the surface in order to preserve/improve mesh quality, without affecting the physics. In their pioneering work, [START_REF] Loewenberg | Numerical simulation of a concentrated emulsion in shear flow[END_REF] update the shape with only normal component of the physical velocity, and introduces a tangential "velocity" field aiming to preserve mesh quality. Subsequent developments [START_REF] Zinchenko | A novel boundary-integral algorithm for viscous interaction of deformable drops[END_REF][START_REF] Cristini | An adaptive mesh algorithm for evolving surfaces: Simulations of drop breakup and coalescence[END_REF][START_REF] Wang | Elastic mesh technique for 3d {BIM} simulation with an application to underwater explosion bubble dynamics[END_REF][START_REF] Zinchenko | Emulsion flow through a packed bed with multiple drop breakup[END_REF] are all based on introducing somehow a "mesh energy", either "kinetic energy" [START_REF] Zinchenko | A novel boundary-integral algorithm for viscous interaction of deformable drops[END_REF][START_REF] Zinchenko | Emulsion flow through a packed bed with multiple drop breakup[END_REF], or "elastic energy" [START_REF] Cristini | An adaptive mesh algorithm for evolving surfaces: Simulations of drop breakup and coalescence[END_REF][START_REF] Wang | Elastic mesh technique for 3d {BIM} simulation with an application to underwater explosion bubble dynamics[END_REF]. In [START_REF] Zinchenko | A novel boundary-integral algorithm for viscous interaction of deformable drops[END_REF], passive stabilization algorithm seeks to preserve the edges current length, while introducing an elastic energy [START_REF] Wang | Elastic mesh technique for 3d {BIM} simulation with an application to underwater explosion bubble dynamics[END_REF] may permit to improve mesh quality over time, for instance when starting with a not-optimal mesh. Adaptivity to local features [START_REF] Cristini | An adaptive mesh algorithm for evolving surfaces: Simulations of drop breakup and coalescence[END_REF][START_REF] Zinchenko | Emulsion flow through a packed bed with multiple drop breakup[END_REF] is done by introducing local reference length scales based on e.g. curvature or element area. The idea of the present algorithm is to treat the mesh as a network of springs, similar to elastic mesh approach, but with varying stiffness aiming to take into account local properties. Doing so allows to keep the network under tension and avoids numerical instabilities that we found using a local reference length approach when this length scales varies too rapidly in the mesh or when the current configuration is too far from optimal. Elastic energy is then minimized by tangential displacements of the nodes along the shape that relax forces in the network. This is done by overdamped dynamics, that is computing on each node the velocity u MQ such that f MQ -µu MQ = 0 where µ is a fictitious viscosity and f MQ is the sum of forces due to network of springs. Since we do not want the redistribution of nodes to change the shape, only the tangent part of the mesh quality velocity is kept, thus u MQ = (I -n ⊗ n) • f MQ /µ. On a given node of position x i , the mesh quality force is computed as:

f MQ i = j∈neighbours (k i + k j ) l ij -l 0 ij t ij l ij (60)
where t ij = x j -x i , l ij = ||t ij || and l 0 ij is an equilibrium length, and the sum is taken over all "neighbors" nodes, that is, nodes connected by one edge to node i. This equilibrium length is defined as α × L ij , where L ij was the distance between x i and x j the last time that the mesh quality algorithm was applied. In practice, mesh quality algorithm is generally applied at the end of each time step, in which case L ij is thus simply the distance between x i and x j at the beginning of the time step. α ≤ 1 is a parameter of the mesh quality algorithm algorithm to allow nodes displacement even if nodes have not moved between two time steps (e.g. a stationary state has been reached). The constants k j are defined such that the mesh quality algorithm tries to maintain an uniform distribution of area of elements (variable A e in 61), while also having a more refined mesh in highly curved regions (variable (2H) 2 e in 61). Thus, k j is computed as:

k j = c A < A > e∈E j A e 3 + c H < (2H) 2 > e∈E j (2H) 2 e 3 (61)
where c A , c H are constants selecting the relative importance of the terms.

A e is the area of the element, (2H) 2 e is the value of the curvature energy integrated over the element, and

< A >= 1 N S dS , < (2H) 2 >= 1 N S (2H) 2 dS (62)
are the average value of area and curvature energy by node, with N the number of nodes in the mesh. The nodes positions are then evolved following:

dx dt * = u remesh (63)
Equation ( 63) is integrated in fictitious time t * by an adaptive third order Runge-Kutta time integrator, until the residual of forces i ||f i || is below a prescribed parameter . The effect of this mesh quality algorithm algorithm is illustrated with parameters (α = 0.9, /c A = 0.05, c H = c A ) which yielded good results for all cases simulated so far, as exemplified for a drop in: shear flow (see 3 (c)), planar extensional flow (see figure 3.6 (b)) and Poiseuille flow (see figure 3.6 (b)). Note that the value of the fictitious viscosity µ only changes the time scale necessary to reach the equilibrium.

Numerical validations and applications of the coupled algorithm

In this section, the algorithm is used to simulate several test cases, in order to validate its implementation and demonstrate possible applications. Since an extensive validation of the numerical algorithm applied to drops has been performed in [START_REF] Gounley | Influence of surface viscosity on droplets in shear flow[END_REF], only an example of drop simulation is presented in section 4.1 and then the numerical algorithm and its implementation are tested with respect to reference results of the literature for capsules and vesicles dynamics. The slice of the shape in the x -y plane is shown, as well as the control mesh vertices projected onto their limit position. For comparison, the slice of the simulation with 5120 elements is superposed (red curve).

To illustrate the application of Loop subdivision elements to drop dynamics, we consider the evolution of a drop in simple shear flow (v ∞ = γye x ), with a capillary number Ca = 0.47 above the critical capillary number for breakup: numerous works have considered similar simulations (in particular, in the context of BEM methods, [START_REF] Cristini | An adaptive mesh algorithm for evolving surfaces: Simulations of drop breakup and coalescence[END_REF][START_REF] Zinchenko | Emulsion flow through a packed bed with multiple drop breakup[END_REF]), but it is still a good illustration of the potential of the method. Starting from an initially spherical shape, the drop elongates and develops a thin neck connecting two bulges, which will then leads to breakup by neck pinching. Since the mesh used here is of fixed topology, the breakup process is not described. For such elongated shapes, a sphero-cylindrical mesh [START_REF] Farutin | 3d numerical simulations of vesicle and inextensible capsule dynamics[END_REF][START_REF] Spann | Loop subdivision surface boundary integral method simulations of vesicles at low reduced volume ratio in shear and extensional flow[END_REF] is more adapted than the classical spherical mesh obtained by regular subdivision of an icosahedron, and is used here. Simulations are performed with three different number of elements : 320, 1280 and 5120 elements. The evolution of the shape at different times is presented in figure 7. The difference between coarse mesh and fine mesh starts to be visible only in the latter stages of the elongation, while shapes obtained with medium and fine resolution cannot be distinguished for the whole duration of the simulation. Note also that at t = 100, the volume drift is around 0.7% for 320 elements and 0.03% for 1280 elements. Thus, even with low resolution and elongated shapes, it is possible to obtain acceptable results, which could be interesting for simulations of suspensions. Capsules dynamics and steady state deformation in simple flows have been studied by several independent groups and thus constitute an excellent validation case. We thus consider the steady state deformation of a Neo-Hookean (NH) capsule and of a Skalak (Sk) capsule, placed either in planar 29 In both cases, an initially spherical capsule of radius R 0 will deform into a roughly ellipsoidal shape. The deformation of the capsule is measured by the Taylor parameter D = L-l L+l , where L and l are semi-axis of the ellipsoid having the same tensor of inertia, and is represented as a function of the capillary number Ca = γR 0 ηext Gs in figures 8 (a) for NH capsule in extensional flow and 8 (b) for Sk capsule in shear flow. As can be seen in these figures, the agreement with results of [START_REF] Dodson | Dynamics of strain-hardening and strainsoftening capsules in strong planar extensional flows via an interfacial spectral boundary element algorithm for elastic membranes[END_REF][START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a stokes flow[END_REF][START_REF] Foessel | Influence of internal viscosity on the large deformation and buckling of a spherical capsule in a simple shear flow[END_REF] is excellent. The meshes used to obtain these results contains between 320 and 20480 elements.

Capsules

A convergence study is carried out for a Sk capsule in shear flow (Ca = 1.5, λ = 1) to assess the dependency on the number of elements, which is varied from 80 to 81920. As shown in figure 9 (a), when the Taylor parameter D or the inclination angle θ are plotted as a function of dimensionless time γt, the results from the different simulations overlap almost perfectly, the only noticeable difference being seen with the coarsest mesh (80 elements). A more important effect of mesh refinement is seen on the preservation of enclosed volume: the relative error v = V 0 -V (t) V 0 is represented as a function of dimensionless time in figure 9 (b), where V 0 = 4 3 πR 3 0 is the volume of the initially spherical capsule. This figure shows that one mesh refinement (multiplying the number of elements by four) leads to roughly gaining one order of magnitude in the volume preservation. However, one can also no- For spatial convergence (a), the RK45 scheme is used and the reference solution is computed with 81920 elements. The spatial convergence of a capsule with bending resistance modeled by Helfrich membrane is also included for two values of the dimensionless bending modulus κ = 10 -4 , 10 -2 . For temporal convergence of trapezoidal scheme, a mesh with 5120 elements is used, and the reference solution is computed with RK45 scheme.

tice that even with the coarsest mesh (80 elements) volume preservation is still acceptable : the volume drift after a dimensionless time of γt = 30 is only 1%.

To estimate the convergence order of the numerical method, the relative error = 10 (a), with the reference solution computed with 81920 elements. As can be seen in the figure 10 (a), the relative error is second order in number of elements O(N -2 ). Note that in [START_REF] Farutin | 3d numerical simulations of vesicle and inextensible capsule dynamics[END_REF], the authors mention that using the singularity-subtraction with standard quadrature will leads asymptotically to a O(N -1 ) convergence because the gradient of the Green function is still singular. For the number of elements used here, this seems not to affect the convergence rate of our method, but could possibly degrade it for finer resolutions, in which case adaptive integration scheme as described in [START_REF] Farutin | Exact singularity subtraction from boundary integral equations in modeling vesicles and red blood cells[END_REF] could be used to restore the O(N -2 ) convergence. To check whether inclusion of bending rigidity affects the convergence order of the method, we also run the same convergence test with two dimensionless bending modulus κ = κ GR 2 0 = 10 -4 , 10 -2 . Note that in this test case, the bending rigidity is added by simply taking into account Helfrich energy for the membrane. Results displayed in figure 10 (a) shows that the convergence order of the method is not affected by inclusion of bending rigidity, still displaying a O(N -2 ) convergence rate. Finally, we also use this test case to check the time convergence order of the trapezoidal time scheme as described in section 3.5 : a mesh containing 5120 elements is used to compute the solution, both with a Runge-Kutta scheme (52) and with trapezoidal scheme (55) for various time step. As expected, the relative error on Taylor parameter converges as O(∆t 2 ) when using the trapezoidal rule. A possible application of our algorithm is to simulate capsules with bending resistance, which is necessary to describe capsule dynamics in the parameter region where compressive stresses are found. The simplest inclusion of bending resistance is to describe the interface energy as the simple sum of a two-dimensional membrane model (e.g., NH, Sk) and a Helfrich membrane. Note that more elaborate model of thin-shells can easily be implemented within the framework presented here, but we choose this ap-proach for its simplicity. We thus simulate the stretching by an uniaxial extensional flow (v ∞ = γ -x 2 e x -y 2 e y + ze z ) of a NH membrane supplemented with Helfrich energy. After an initial transient, the capsule reach a steady stated showed in figure 4 We then test our numerical algorithm applied to vesicles on the prototypical configuration of a vesicle immersed in shear flow : v ∞ = γye x . In this configuration, several dynamics may be observed, the simplest being a tank-treading configuration where the vesicle reach a steady deformation and inclination angle similarly to drops or capsules. This steady angle depends on the viscosity contrast λ and on the excess area ∆ which measures the deformability of the vesicle. Indeed, since for a vesicle both enclosed volume and total surface must be preserved, an initially spherical vesicle will not deform, whatever the flow strength is. To be deformed, the vesicle must thus be deflated, and a measure of this deflation is provided either by the reduced volume v = . The evolution of the inclination angle as a function of excess area for a vesicle without viscosity contrast is reported in figure 12 (a) where it is compared to spectral simulations of [START_REF] Zhao | The dynamics of a vesicle in simple shear flow[END_REF] and simulations of [START_REF] Yazdani | Three-dimensional numerical simulation of vesicle dynamics using a front-tracking method[END_REF]. The agreement with [START_REF] Zhao | The dynamics of a vesicle in simple shear flow[END_REF] is excellent, even with 320 elements.

If one includes viscosity contrast, different dynamics can be observed, such as tumbling, where the vesicle inclination axis undergoes a full rotation. In between tank-treading and tumbling, for sufficiently high capillary number, a more complex dynamics termed vacillating-breathing [START_REF] Misbah | Vacillating breathing and tumbling of vesicles under shear flow[END_REF] or trembling [START_REF] Deschamps | Phase diagram of single vesicle dynamical states in shear flow[END_REF][START_REF] Zhao | The dynamics of a vesicle in simple shear flow[END_REF] can be observed, where the inclination axis undergoes oscillation without full rotations, while the shape also oscillates. The evolution of the VB-TR / TU period is plotted as a function of viscosity contrast in figure 12 (b) and compared to spectral simulations of [START_REF] Zhao | The dynamics of a vesicle in simple shear flow[END_REF] and simulations of [START_REF] Farutin | 3d numerical simulations of vesicle and inextensible capsule dynamics[END_REF], showing an excellent agreement with both, except for λ = 9 where results deviates slightly from [START_REF] Zhao | The dynamics of a vesicle in simple shear flow[END_REF], but seems consistent with [START_REF] Farutin | 3d numerical simulations of vesicle and inextensible capsule dynamics[END_REF] in doing so. For vesicles, both the enclosed volume and the total surface should be preserved. The relative error for a tank treading vesicle is shown in figure 13 (a), and for a tumbling vesicle in figure 13 (b). In figure 13 (a), influence of 34 mesh size and time step on reduced volume preservation is also studied. The reduced volume drifts for long simulation times (stationary state is reached around γt ≈ 5), but even with the coarsest mesh the relative error is below 0.5% after several full tank-treading period. As in figure 9 (a), refining the mesh once leads to an order of magnitude decrease in the variation of reduced volume, because the slope of the drift depends on the mesh used.

There is also an important effect of the time step used, because for 1280 elements, the major part of the variation of the enclosed volume occurs in the transient part where the vesicles deforms from its initial ellipsoidal shape towards the stationary tank treading shape. Significant gains in reduced volume preservation could thus be obtained by using an adaptive implicit time scheme, as has been done recently by [START_REF] Quaife | Adaptive time stepping for vesicle suspensions[END_REF]. We also simulate a deflated v = 0.65 vesicle in uniaxial extensional flow (v ∞ = γ -x 2 e x -y 2 e y + ze z ), which is known [START_REF] Zhao | The shape stability of a lipid vesicle in a uniaxial extensional flow[END_REF][START_REF] Spann | Loop subdivision surface boundary integral method simulations of vesicles at low reduced volume ratio in shear and extensional flow[END_REF][START_REF] Vivek Narsimhan | The mechanism of shape instability for a vesicle in extensional flow[END_REF][START_REF] Dahl | Experimental observation of the asymmetric instability of intermediate-reduced-volume vesicles in extensional flow[END_REF] to exhibit instability. For comparison, we follow the procedure described in [START_REF] Zhao | The shape stability of a lipid vesicle in a uniaxial extensional flow[END_REF][START_REF] Spann | Loop subdivision surface boundary integral method simulations of vesicles at low reduced volume ratio in shear and extensional flow[END_REF]: starting from the equilibrium shape at v = 0.65, we first simulate the stationary state under flow, with Ca A = 4. Here, Ca A = γη κ A 4π 3 2 is the capillary number with a reference length based on the area. After reaching the stationary state, the shape is then perturbed with a sinusoidal perturbation, and the L 2 -norm of the radial velocity is monitored as a function of time, as in [START_REF] Zhao | The shape stability of a lipid vesicle in a uniaxial extensional flow[END_REF][START_REF] Spann | Loop subdivision surface boundary integral method simulations of vesicles at low reduced volume ratio in shear and extensional flow[END_REF]. Snapshots of the shape at different dimensionless times γt as well as evolution of L 2 norm of radial velocity are shown in figure 14, and show an excellent agreement with reference results of [START_REF] Zhao | The shape stability of a lipid vesicle in a uniaxial extensional flow[END_REF][START_REF] Spann | Loop subdivision surface boundary integral method simulations of vesicles at low reduced volume ratio in shear and extensional flow[END_REF]. Finally, we simulate vesicles settling due to a density contrast between the inner and outer fluids, and compare with results obtained with an axisymmetric code using linear elements [START_REF] Boedec | Sedimentation-induced tether on a settling vesicle[END_REF][START_REF] Trozzo | Axisymmetric boundary element method for vesicles in a capillary[END_REF]. The dimensionless parameter controlling the shape is the Bond number, defined as Bo = ∆ρgR 4 0 κ . Starting with a prolate ellipsoid of excess area ∆ = 0.93, shape evolves toward a rounded vesicle at the front with a protrusion at the rear, whose length reaches a plateau at stationary state. Depending on the strength of the forcing, shapes varies continuously from pear-like shapes to tethered shapes. Comparison of profiles showed in figure 15 show an excellent agreement with [START_REF] Boedec | Sedimentation-induced tether on a settling vesicle[END_REF][START_REF] Trozzo | Axisymmetric boundary element method for vesicles in a capillary[END_REF].

Conclusion

In this paper, we have presented an algorithm to simulate the dynamics of soft objects such as drops, capsules or vesicles under Stokes flow. The algorithm is based on the use of Loop subdivision surfaces to describe the interface of these objects. This surface representation is used both for the finite element solver for membrane forces and for the boundary element method for the fluid solver. The surface incompressibility constraint for vesicles is treated by the introduction of a local Lagrange multiplier which is solved by an iterative method. To avoid the growth of the number of iterations, a preconditioner based on a linear approximation of the surface [START_REF] Boedec | 3d vesicle dynamics simulations with a linearly triangulated surface[END_REF] is introduced. The whole algorithm spatial discretization error decreases as O(N -2 ), where N is the typical mesh number of elements (or vertices), and is coupled with either a second order implicit time scheme for simulation with bending rigidity, or with an adaptive RK45 time scheme for capsule or drop simulations.

The combination of an high order spatial scheme with second order time algorithm allows a very good preservation of invariants such as enclosed volume or area for vesicles, without any correction on the velocity field, introduction of a reference/target volume or area, or rescaling of the shape after each time step. Note that it is of course possible to include these corrections into the proposed algorithm without difficulties to eliminate possible drift of the solution at long simulation times, but these corrections cannot be used for simulations of multiple objects where the rescaling of the shape will lead to a modification of inter-objects distances, which might crucially affect the dynamics. We plan to improve the conservation of invariants by using adaptive quadratures to limit volume drift and adaptive time stepping with implicit algorithm to limit both volume and area variations, such as the recent algorithm of [START_REF] Quaife | Adaptive time stepping for vesicle suspensions[END_REF] designed for 2D vesicles suspensions.

Finally, the algorithm presented here can be extended to account for more complex interfacial rheology (e.g. surface viscosity effects, Marangoni effects). Its extension to drops with viscous interfaces has been discussed in [START_REF] Gounley | Influence of surface viscosity on droplets in shear flow[END_REF]. The inclusion of surface viscosity for other soft objects (capsules, vesicles) is likely to strongly influence the dynamics, and might be necessary to have a more accurate representation of real capsules, as recent experiments suggest [START_REF] De Loubens | Tank-treading of microcapsules in shear flow[END_REF].

Figure 1 :

 1 Figure 1: Simulation of a Skalak capsule in shear flow v = γye x , showing (a) the control mesh, composed of 80 elements, (b) the mesh obtained after one iteration of Loop subdivision and (c) the limit surface corresponding to an infinite number of iterations of Loop subdivision. The original control mesh is also represented in black wireframe, and the limit position of control vertices on the limit surface is represented by black dots. First row is the view in the xy plane, and second row the view in the xz plane.

Figure 2 :

 2 Figure 2: (a) Scheme of an element (gray triangle) with its local parametrization. The position inside the element depends on the nodal values of position of every node connected to the central element by at most one edge, denoted the one-ring. (b) Support of the shape function associated with one node extends up to every triangle which is at most one edge afar from the node: the control mesh is composed of a flat mesh, except one node whose height is set to 1.

Figure 3 :

 3 Figure 3: (color online) Effect of mesh quality preserving algorithm on a tank-treading drop Ca = 0.3 : (a) evolution of Taylor parameter as a function of dimensionless time. (b) relative error on volume conservation (c) Snapshots of the shape and mesh at different dimensionless times (from top to bottom) : γt = 0; 3; 6; 15; 30. Left (respectively right) part of the shape is a simulation without (respectively with) mesh quality preservation algorithm.

Figure 4 :

 4 Figure 4: Influence of control mesh quality on the solution at γt = 37.5 (parameters are the same as in figure 3. (a) control mesh nodes projected onto the limit shape, and cross section of the shape in the z = 0; 0.5; 0.75; 0.85 planes. Red (resp. black) curve is with (respectively without) mesh quality preserving algorithm. (b) comparison of the membrane force norm on the different cross-section. The curves are plotted as a function of θ, the polar angle defined with respect to x-axis.

Figure 5 :

 5 Figure 5: (a) Evolution of Taylor parameter for a clean drop in extensional flow as a function of dimensionless time. (Inset) Evolution of relative error on enclosed volume as function of dimensionless time. (b) Snapshots of the shape of a clean drop in extensional flow at different dimensionless times.

Figure 6 :

 6 Figure 6: Effect of mesh quality preserving algorithm on a drop in Poiseuille flow (Ca = 0.5). Left part of the drop comes from a simulation without any node redistribution, right part comes from a simulation with application of algorithm described in section 3.6.

Figure 7 :

 7 Figure 7: (color online) Evolution of a drop in shear flow with Ca = 0.47, with (a) 320 and (b) 1280 elements. From top to bottom, times are : t = 0; 25; 50; 75; 95; 100.The slice of the shape in the x -y plane is shown, as well as the control mesh vertices projected onto their limit position. For comparison, the slice of the simulation with 5120 elements is superposed (red curve).

  λ=1 Present work, λ=5 Walter et al., 2010 (λ=1) Foessel et al., 2011 (λ=5)

Figure 8 :

 8 Figure 8: Comparison of present work with literature results for Taylor parameter D as a function of capillary number Ca for (a) neo-Hookean capsule in extensional flow [18, 55] (b) Skalak capsule in shear flow[START_REF] Walter | Coupling of finite element and boundary integral methods for a capsule in a stokes flow[END_REF][START_REF] Foessel | Influence of internal viscosity on the large deformation and buckling of a spherical capsule in a simple shear flow[END_REF]

Figure 9 :

 9 Figure 9: (a) Evolution of Taylor parameter and inclination angle of a Skalak capsule in shear flow (Ca = 1.5) as a function of time, for different number of elements. (b) Corresponding evolution of the relative error on enclosed volume as a function of time, for different number of elements.

Figure 10 :

 10 Figure10: Convergence study of the numerical algorithm : relative error on the steady state value (taken at t = 26.25) of the Taylor parameter of a Skalak capsule under shear flow (Ca = 1.5) with respect to (a) number of elements (b) time step size for trapezoidal scheme. For spatial convergence (a), the RK45 scheme is used and the reference solution is computed with 81920 elements. The spatial convergence of a capsule with bending resistance modeled by Helfrich membrane is also included for two values of the dimensionless bending modulus κ = 10 -4 , 10 -2 . For temporal convergence of trapezoidal scheme, a mesh with 5120 elements is used, and the reference solution is computed with RK45 scheme.

D

  ref -D D ref on Taylor parameter D at a dimensionless time γt = 26.25 (thus in stationary state) is represented as a function of the number of elements in figure

Figure 11 :

 11 Figure 11: Neo-Hookean capsule in uniaxial extensional flow (Ca = 0, 07; κ = 10 -4 ) : (a) Top -view in the x-y plane . Bottom -view in the y-z plane (b) Cross section in the z = 0 plane.

. 2 .

 2 Runs at three different resolutions (N = 5120, 20480, 81920 elements) have been performed, and the comparison of the cross-section in the z = 0 plane is shown in figure 4.2 (b), showing an excellent agreement for all resolutions.

Figure 12 :

 12 Figure 12: (a) Steady state inclination angle θ for a tank-treading vesicle in simple shear flow (λ = 1, Ca = 10) as a function of excess area ∆. (b) Dimensionless VB-TR/TU period (T ) as a function of viscosity contrast λ for a quasi-spherical vesicle (v = 0.95, Ca = 6.0)

  or by the excess area ∆ = A 4π 3V 2/3 -4π. The ratio between external flow time scale and bending time scale defines the capillary number : Ca = γηextR 3 0 κ

Figure 13 :

 13 Figure 13: (a) Relative error on reduced volume for a tank treading vesicle in simple shear flow (Ca = 10, λ = 1, v = 0.95), comparison of two meshes and two different time steps. (b) Relative error on volume and area for a tumbling vesicle in shear flow (Ca = 6, λ = 20), for 1280 elements and a time step of ∆t = 10 -2 .

Figure 14 :

 14 Figure14:(a) Snapshots of shape at different dimensionless times γt = 0; 18.1; 23.3; 25.9; 27.2; 28.5 (b) Evolution of L 2 norm of the radial velocity. The dashed line corresponds to an exponential growth with a growth rate of ω = 8.57 × 10 -2 as computed by[START_REF] Zhao | The shape stability of a lipid vesicle in a uniaxial extensional flow[END_REF][START_REF] Spann | Loop subdivision surface boundary integral method simulations of vesicles at low reduced volume ratio in shear and extensional flow[END_REF].

Figure 15 :

 15 Figure 15: (a) Stationary shape of a settling deflated vesicle (= 0.93, Bo = 216) (b) Stationary profiles for different Bond numbers (Bo = 59; 118; 196). Symbols are present results and solid lines are comparison with axisymmetric code [9, 52].

Table 1 :

 1 Tension solver preconditioning : number of GMRES iterations needed to ensure a = 10 -9 residual for the surface divergence free constraint (equation 43), for a vesicle under shear flow v = γxe y , under various flow conditions measured by capillary number Ca = γηextR 3

	1 1 10 0.1 10 0.1 1 10 0.1 10
	320 18 18 18 18 18 10 12 12 13 13
	1280 39 42 44 46 43 13 15 16 17 16
	5120 68 76 78 81 77 14 17 18 19 13

κ and viscosity contrast λ = ηint ηext .

Table 2 :

 2 Average GMRES iteration count for the resolution of equation[START_REF] Zhao | A spectral boundary integral method for flowing blood cells[END_REF], the average number of Newton iterations (r) necessary to solve[START_REF] Yazdani | Three-dimensional numerical simulation of vesicle dynamics using a front-tracking method[END_REF] is indicated in parenthesis. The average is performed over the 10 first time steps of a simulation of a vesicle in shear flow with Ca = 1.

Acknowledgements

The authors would like to thank Paul Gang Chen for careful reading of the manuscript. This work has benefited from financial support from the ANR Polytransflow (grant no.13-BS09-0015-01), from Labex MEC (grant no. ANR-11-LABX-0092), from A*MIDEX (grant no. ANR-11-IDEX-0001-02) and from CNES. This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01).