D Gloaguen
email: gloaguen@lamm.univ-nantes.fr

S Freour

M François

R Guillen

J Royer

Thermal stresses in hexagonal materials -heat treatment influence on their mechanical behaviour

Internal stresses due to anisotropic thermal and plastic properties were investigated in rolled zirconium and titanium. The thermal stresses induced by a cooling process were predicted using a selfconsistent model and compared with experimental results obtained by X-ray diffraction. The study of the elastoplastic response during uniaxial loading was performed along the rolling and the transverse direction of the sheet, considering the influence of the texture and the thermal stresses on the mechanical behaviour. An approach in order to determine the thermal behaviour of phases embedded in two-phase materials is also presented.

INTRODUCTION

Hexagonal close packed materials present anisotropic thermal, elastic and plastic properties at mesoscopic and macroscopic levels. These properties and the crystallographic texture explain the appearance and the development of internal stresses when a thermo-mechanical deformation is introduced during the manufacturing of the material. Two types of internal stresses can be distinguished: Thermal internal stresses, which may arise from thermal anisotropy. As single crystal thermal expansion coefficients are different in the direction of the c r axis and in the basal plane, intergranular thermal stresses are induced during the cooling process after a heat treatment because a full accommodation of the anisotropic thermal contraction of each grain is prevented by the neighboring grains. These stresses can modify the elastoplastic response of the polycrystalline material [1]. Mechanical internal stresses, induced by the elastoplastic anisotropy and the requirement of compatibility [2]. In this article, we present a systematic simulation of the thermal and mechanical response of rolled zirconium and titanium sheets, using the self-consistent elastoplastic formulation proposed by Hill [3]. We model the cooling process and the tensile elastoplastic response of samples, and calculate the evolution of internal stresses associated with these processes. We quantify precisely the level and the influence of the thermal stresses on the macroscopic behaviour. We combine experimental observations (X-ray diffraction, tensile tests) with the predictions of a self-consistent model in order to obtain more informations about the different factors responsible for the appearance of residual elastic mesostrains. This comparison allows a better understanding and interpretation of X-ray diffraction and mechanical results. X-ray diffraction technique is one of the most used to characterise the mechanical state at mesoscopic level of polycrystalline materials. It is performed from the elastic strains measurement and the subsequent calculation of the stress state by means of mechanics equations. The capacity to measure intergranular strains provides an experimental tool for understanding how intergranular strains are generated. Experimental data from X-ray diffraction offer a rigorous test at a microscopic level for the models [4][5].

In a second part, we present an approach developed in order to predict the internal stresses and strains due to thermal solicitations in multiphase metallic polycrystals. The procedure is based on the coupling of numerical simulations with scale transition models and experiments. The effects due to thermal interactions occurring in two-phases materials are analyse.

SELF-CONSISTENT SCHEME

The elasto-plastic self-consistent model predicts intergranular effects by accounting for the mechanical and thermal anisotropy and the texture of the polycristal. Thermal effects are accounted for here. We assume that thermal contractions do not introduce plastic deformation in the material during the cooling process [6]. In this case, a thermoelastic formulation is sufficient to simulate the evolution of thermal stresses associated with cooling. The isothermal uniaxial loading is an elastoplastic problem where thermal effects enter through the residual stresses left in the material by the previous cooling process.

Polycristal models [7] treat each grain as an inclusion embedded in a Homogeneous Effective Medium (HEM) exhibiting the same overall response of the aggregate. Each grain is treated as an inhomogeneity with elastic (c) and thermal (a) constants, slip and twin mechanisms characteristic of a single crystal. The texture is known, the polycristal is represented by a weighted discrete distribution of orientations The stress and strain rates, σ & and ε & , of each grain are linearly linked by an elastoplastic tensor l. This tensor accounts for anisotropy of the grains and the hardening effects of the deformation systems during plasticity. The number and type of active systems are defined through their critical resolved shear stress and some hardening behaviour. The volume averages of the local stress and strain tensors must coincide with the overall strain and stress, E and S. Using the Eshelby's inclusion formalism [8], the σ & and ε & of each grain are coupled to Σ & and E & through an interaction tensor L * . It is a function of the Eshelby tensor, the elastoplastic tensor of the polycristal L, and the shape and the orientation of the inclusion. L linked linearly Σ & and E & .

Since the properties of the medium derive from the average response of all grains, it is initially undetermined and must be solved by iteration. By specifying an overall stress or strain, the model can give the corresponding stress or strain tensors for each orientation. Thus, we can describe the mechanical response of the polycrystal at mesoscopic scale.

Data used in the simulations

For the cooling simulation, we supposed the sample is stress-free at 923 K [2]. The texture was introduced in the model by a set of 1000 grains characterized by Euler angles and weights which represent their volume fraction. The single crystal elastic constants and the thermal expansion coefficients used in the model are listed in [START_REF] Ortiz | Fourth European Conference on Residual Stresses[END_REF]. We are taking account the temperature dependence of these coefficients.

For the tensile loading, the temperature is kept constant, the initial internal stresses are the ones left by the cooling process. The main active slip systems are assumed to be : the 6 prismatics { } 0 twinning. We considered a linear hardening law [4], the coefficient H gr is equal to H gg for any deformation modes r :

å = r r γ g H g τ & &
, the CRSS for the prismatic systems is determined to reproduce the experimental yield stress in the rolling direction (RD). The CRSS and hardening coefficients values are adjusted to give a good description of the tensile test performed in the RD of the sheet. The same parameters are taken for the simulation of the tensile test in the transverse direction (TD) and X-ray diffraction results.

Stress measurements by X-ray diffraction technique

These experiments were carried out on a D500 SIEMENS goniometer with a Cr K a (Zr702) and Cu K a radiations (Ti40). An W goniometric assembly with a scintillation detector was used. Two plane families were studied : {

RESULTS AND DISCUSSION

Thermal stresses

The thermal strains due to a temperature variation from 923 to 298 K are calculated. Predicted results are compared with experimental data using X-ray diffraction method.

Using the self-consistent formulation, thermal residual strains projected on the measurement direction versus sin 2 y for the diffracting planes were simulated. These strains can be represented by a straight line A + B sin 2 y. Predicted results were compared with experimental data using X-ray diffraction technique. The experimental strain value are determined by : ( )

0 0 φψ φψ θ cotan 2θ 2θ 2 1 ε - - = (1)
where 0 q is the Bragg angle for a stress-free specimen of the studied a phase and fy q is the measured Bragg angle. Consequently, it is necessary to determine the stress-free lattice parameters of the sample to calculate the thermal strains. These parameters cannot be measured using polycrystalline sample, because they are distorted by intergranular strains induced by the anisotropic thermal behaviour of the grains. Then, one solution [6] is to find the stress-free parameters a 0 and c 0 chosen, such that the difference between predicted and experimental strains values were minimized. The self-consistent formulation was used in order to determine the stress free lattice parameters. It should be noted that this determination was performed using the model in the reverse way. In this case, single crystal properties are determined by the comparison with experimental measurements done on aggregates which are the only possible source of data.

Figure 1 shows the experimental (solid circles and triangles) and predicted (full line) thermal strain versus sin 2 y curves for { } (titanium) diffracting planes with f = 0°. We obtained a good agreement between the two approaches. To give an idea of the stress magnitude, the diagonal components

Mechanical stresses

The mechanical response in the elastoplastic region is calculated using the elastoplastic self-consistent formulation presented above and the initial residual stresses left by the previous cooling process. Uniaxial tensile test up to 2 % strain along the RD and the TD are simulated. The results indicate that the elastoplastic response is strongly influenced by the residual stresses and by the texture of the material. The experimental and predicted tests for zirconium and titanium are depicted in figures 2 and 3, respectively.

The experimental tests for zirconium show several differences between the two directions : the yield stress is higher, the elastoplastic transition is faster and the curve shows a plateau at 0.6%strain in the transverse direction. The measured yield stresses are 205 MPa and 290 MPa (stress differential of 85 MPa) in the RD an TD respectively. The corresponding Young's modulus are 93 GPa and 96 GPa. These values are similar because of the low elastic anisotropy. With the model developed here, the predicted Young's modulus are 93 GPa and 95 GPa in agreement with the experimental ones. The stress differential can be attributed to thermal residual stresses. This fact has been already analysed by Mac Ewen et al [START_REF] Mac Ewen | [END_REF] on Zircaloy-2. The predicted yield stresses are 205 MPa and 275 MPa in the RD an TD respectively. When no thermal stresses are assumed, the model predicts a low stress differential (<15 MPa) which could be attributed to texture. The model fails to reproduce the plateau in the TD because the stress relaxation (due to twinning for example) is not accounted for. Nevertheless, the simulated data exhibit the correct tendency.

For titanium, the predicted yield stress is about 251 MPa (0.23 % total strain) along the RD. This value is very close to the experimental one (253 MPa). The calculated Young's modulus (107 GPa) coincides with the experimental one. With a loading along the TD, the initiation of the elastoplastic transition begins at 271 MPa (0.27 % strain). The experimental yield stress is about 293 MPa. Although the predicted yield stress is lower than experimental one, the simulation matches the fact that the yield stress in TD is higher than in RD. To quantify the influence of thermal stresses on the macroscopic response, calculations were also performed without thermal stress in the initial conditions. The predicted stress-strain curves are identical to the previous case. The elastoplastic transition begins at 256 MPa (251 MPa with thermal stresses) in the RD anisotropy is very low and its influence on the macroscopic behaviour, contrary to zirconium, is negligible.

In the simulations, we have considered that prismatic glide is the main active mechanism, pyramidal slip is a secondary slip mode and twinning contributes very little to plastic deformation. Distribution was obtained in agreement with the published experimental observations.

THERMAL STRESSES IN TWO-PHASE MATERIALS

Let us consider a hypothetical two-phase system composed of one phase having the same properties than a-Zr, while the second one will be assume to behave as a-Ti. This material can be only designed in theory because titanium and zirconium are, in practice, miscible in any proportion. Nevertheless, this material has been supposed for the need of the simulation, due to the lack of data concerning the evolution with the temperature of the CTE and elastic properties of real hexagonal two-phase structures. The model used for our simulations was presented in a previous paper [11].

The behaviour has been simulated when the material is subjected to a cooling process from 998 K to 298 K. Five compositions (volume fractions) of the polycrystal have been considered : f Zr = {0 ; 0.25 ; 0.5 ; 0.75 ; 1}. In each case, the homogenisation procedure has been performed on 4000 crystallites per phase. A spherical shape has been chosen for the inclusions. The volume fraction of each crystallite was adjusted in order to satisfy the assumption on the volume fraction of the phase in the material. The texture is isotropic in each phase, so that the relative weight of the grains are the same in a given phase. Notice that the case when the volume fraction of a phase tends towards zero corresponds to the limit case of the inclusion embedded in an infinite matrix developed in Eshelby's model [24]. The measurement direction chosen in the point of view of the diffraction is defined by f=0° and y=0°.

The results obtained for the planes { } The amount of residual stresses generally increases with a decrease in the volume fraction of the considered phase. In a-Ti, the axial stress component simulated for the plane { } 2 2 30 vary from 0 MPa (Ti single phase) up to 286 MPa, (one inclusion of titanium in a zirconium matrix). Thus, in specific conditions, the induced residual thermal stresses may reach values very close to the yield limit of a phase. The residual strains computed owing to SC model increase in a-Ti as a function of the volume fraction of a-Zr. Conversely, in the same conditions, the strains decrease in a-Zr. These results show clearly the influence of the second phase on the residual stresses values in the studied phase.

CONCLUSION

In this study, X-Ray Diffraction and self-consistent modeling have been used to measure and simulate the thermal residual strains induced in either single-phase or two-phase polycrystalline materials. Comparisons of these two approaches performed on commercially pure a-Ti and a-Zr textured samples show a good agreement. The efficiency of SC estimates was also demonstrated. In particular, the model properly reproduces the behaviour of the diffracting volume of strongly thermal anisotropic materials as Zr.

For zirconium, the residual stresses due to thermal anisotropy are rather important (equivalent to 35% of the yield stress) and consequently they play an important role on the elastoplastic transition, contrary to titanium. We showed that the crystallographic texture and the thermal residual stresses can explain the stress differential between RD and TD of the sample the elastoplastic transition after an uniaxial loading. The study of two-phase material show the influence and the interaction of the second phase on the thermal behaviour in the studied phase.

  2q = 148.751° for titanium. Diffractograms were recorded for fifteen tilt angles y varying between -45 and 45° and for azimuth angles f = 0, 45 and 90°. The direction f = 0° corresponds to the tensile direction.

  The values obtained for the best fit are : Zr : a = 3.2365 Å and c = 5.1498 Å or 2

22 s and 33 sFigure 1 .

 331 Figure 1. Intergranular deformations as a function of sin 2 y at room temperature with f = 0° for a) -zirconium -• { 4 1 10 } and ▲ { 2 2 20 } diffracting plane -experimental result. ─ predicted results.

Figure 2 .Figure 3 .

 23 Figure 2. Predicted and measured stress-strain curves for uniaxial loading along the TD (a) and the RD (b) of the zirconium sheet.

  a-Ti are given on figure4. In each case, f denotes the volume fraction of the diffracting phase.

Figure 4 .

 4 Figure 4. predicted thermal induced residual strains on the diffracting volume of a-Ti embedded in a hypothetical Zr-Ti matrix, as a function of the temperature and phase volume fraction.

{30.2} plane -Ti in


J Walter

A.-V Salsac

D Barthès-Biesel

P Le Tallec

Coupling of finite element and boundary integral methods for a capsule in a Stokes flow

Keywords: fluid-structure interaction, finite element method, boundary integral method, membrane model

We introduce a new numerical method to model the fluid-structure interaction between a microcapsule and an external flow. An explicit finite element method is used to model the large deformation of the capsule wall, which is treated as a bidimensional hyperelastic membrane. It is coupled with a boundary integral method to solve for the internal and external Stokes flows. Our results are compared with previous studies in two classical test cases: a capsule in a simple shear flow and in a planar hyperbolic flow. The method is found to be numerically stable, even when the membrane undergoes in-plane compression, which had been shown to be a destabilizing factor for other methods. The results are in very good agreement with the literature. When the viscous forces are increased with respect to the membrane elastic forces, three regimes are found for both flow cases. Our method allows a precise characterization of the critical parameters governing the transitions.

INTRODUCTION

Synthetic capsules are liquid drops or gels protected by a thin elastic membrane and have a very wide range of applications in cosmetic, food, and pharmaceutical industries. They are used to protect fragile, volatile, or active substances and to control the release rate as desired. Other instances are found in biomedical engineering, where microcapsules are used for cell encapsulation or artificial blood manufacturing [START_REF] Kühtreiber | Cell Encapsulation Technology and Therapeutics[END_REF][START_REF] Chang | Future prospects for artificial blood[END_REF]. In many applications, what is at stake is the characterization of the complex behavior of a deformable capsule in an external liquid flow, to control or prevent the breakup of the capsule and the liberation of its contents. Numerical models of the fluid-structure interactions are thus necessary to predict the capsule deformation in flow conditions and compute the stress level in the membrane.

In many respects, the fluid-structure interaction of microcapsules suspended in an external flow is unconventional. [START_REF] Kühtreiber | Cell Encapsulation Technology and Therapeutics[END_REF] The capsules are closed surfaces freely suspended in a fluid. They are, therefore, not subjected to any boundary condition on the displacement apart from spatial periodicity. [START_REF] Chang | Future prospects for artificial blood[END_REF] At the small scale of the capsule, the inertia forces of the internal and external flows are negligible compared with the viscous forces. The fluid flows can, therefore, be modeled using the Stokes equations (balance between viscous and pressure forces). The load on the membrane is then due to pressure forces and viscous tractions with normal and in-plane components. [START_REF] Risso | Experimental investigation of a bioartificial capsule flowing in a narrow tube[END_REF] The inertia force of the capsule wall can be equally neglected. Consequently, the problem does not contain any dynamics and can be seen as a succession of equilibrium states. The drawback is that the stability of the equilibrium is not necessarily satisfied. [START_REF] Lefebvre | Flow of artificial microcapsules in microfluidic channels: a method for determining the elastic properties of the membrane[END_REF] The capsule wall is subjected to not only large displacements but also large deformations. For instance, average stretch ratios of up to 1.2 are reported in experimental studies of capsules flowing in a pore [START_REF] Risso | Experimental investigation of a bioartificial capsule flowing in a narrow tube[END_REF][START_REF] Lefebvre | Flow of artificial microcapsules in microfluidic channels: a method for determining the elastic properties of the membrane[END_REF]. The non-linear effects due to these large deformations must be taken into account in the membrane mechanics model.

An analytical solution can only be found in the approximation of small deformations, using perturbation methods [START_REF] Barthès-Biesel | The time-dependent deformation of a capsule freely suspended in a linear shear flow[END_REF]. In the case of large deformations, a numerical simulation is needed to solve the fully coupled equations governing the fluid and membrane mechanics. This problem has been studied for the past three decades and different techniques have been considered to derive numerical solutions (as summarized in Table I). Many studies have used a coupling strategy based on the boundary integral method to solve the Stokes flow and the membrane elasticity equations on the same mesh [START_REF] Li | Large deformations and burst of a capsule freely suspended in an elongational flow[END_REF][START_REF] Ramanujan | Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of capsule viscosity[END_REF][START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF][START_REF] Dodson | Spindles cusps and bifurcation for capsules in Stokes flow[END_REF]. The velocity field at any position within the fluid domain is given by surface integrals calculated on the geometric boundaries. This method, therefore, has the advantage of reducing the geometric dimension of the problem by one, which largely decreases the total number of nodes. It also does not involve discretizing spatial derivatives (using for example finite differences) and is very accurate. In particular, it has been shown to be efficient, precise, and stable when modeling the deformation of capsules subjected to shear flows. Another coupling method used for capsule simulation is the immersed boundary method [START_REF] Eggleton | Large deformation of red blood cell ghosts in a simple shear flow[END_REF][START_REF] Doddi | Effect of inertia on the hydrodynamic interaction between two liquid capsules in simple shear flow[END_REF][START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF]. Two grids are used in the method: a stationary 3D grid for the fluid flow and a moving 2D boundary grid for the interface. The forces exerted by the membrane on the fluid and the flow velocity convecting the membrane are applied locally from one grid to the other using approximate Dirac functions [START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Mittal | Immersed boundary methods[END_REF]. The Navier-Stokes equations are generally solved using finite difference schemes. The immersed boundary method has the advantage over the boundary integral method to be applicable to non-zero Reynolds numbers and non-Newtonian fluids. However, the method does not treat the membrane as a physical boundary of the fluid domain, which can result in a lack of precision in the Lagrangian tracking of the interface.

The most frequently used model for the capsule wall is that of a 2D hyperelastic surface: the wall is considered to be infinitely thin and to have a negligible bending stiffness. Two approaches may be considered to model the capsule membrane mechanics: the equations of the force equilibrium on the capsule wall may either be written locally at each point (strong form) or multiplied by a test function and integrated over the capsule surface to define a variational problem (weak form). Most capsule studies have used the strong form of the equations. Numerical simulations of such models were first performed in axisymmetric configurations with a collocation technique to determine the unknown forces and velocities at discrete locations in elongational flows [START_REF] Li | Large deformations and burst of a capsule freely suspended in an elongational flow[END_REF][START_REF] Pozrikidis | The axisymmetric deformation of a red blood cell in uniaxial straining Stokes flow[END_REF] and in pore 3D with finite elements [START_REF] Eggleton | Large deformation of red blood cell ghosts in a simple shear flow[END_REF][START_REF] Doddi | Effect of inertia on the hydrodynamic interaction between two liquid capsules in simple shear flow[END_REF][START_REF] Doddi | Lateral migration of a capsule in a plane Poiseuille flow in a channel[END_REF] The capsule wall was each time modeled as a 2D membrane.

flows [START_REF] Leyrat-Maurin | Motion of a deformable capsule through a hyperbolic constriction[END_REF][START_REF] Quéguiner | Axisymmetric motion of capsules through cylindrical channels[END_REF]. The full 3D problem was then solved for infinite flow configurations. Capsules in simple shear flows have been considered by Pozrikidis [START_REF] Pozrikidis | Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow[END_REF], Ramanujan and Pozrikidis [START_REF] Ramanujan | Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of capsule viscosity[END_REF] and more recently by Li and Sarkar [START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF], who computed the membrane load as a piecewise constant function. Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF] and Lac and Barthès-Biesel [START_REF] Lac | Deformation of a capsule in simple shear flow: effect of membrane prestress[END_REF] used bi-cubic B-splines instead as interpolation functions to compute the loads with high accuracy, whereas Dodson and Dimitrakopoulos [START_REF] Dodson | Spindles cusps and bifurcation for capsules in Stokes flow[END_REF] opted for a spectral discretization of the problem. An alternative option is to write the equilibrium equations in their weak form and to use a finite element method. The local equilibrium equations are converted into a variational problem. Only two groups [START_REF] Eggleton | Large deformation of red blood cell ghosts in a simple shear flow[END_REF][START_REF] Doddi | Effect of inertia on the hydrodynamic interaction between two liquid capsules in simple shear flow[END_REF][START_REF] Doddi | Lateral migration of a capsule in a plane Poiseuille flow in a channel[END_REF] have implemented a finite element method that was each time coupled with an immersed boundary method to compute the motion of a capsule in shear flows. However, their method [START_REF] Charrier | Free and constrained inflation of elastic membranes in relation to thermoforming-non-axisymmetric problems[END_REF] is based on the use of linear interpolation functions and lacks generality in the implementation.

Previous studies have shown that, for certain ranges of the governing parameters, the capsule wall undergoes in-plane compression and tends to buckle [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF][START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF]. Thus, bending effects must be taken into account to properly model the wall mechanics. However, no satisfying model for the bending stiffness of a capsule has been implemented yet. The finite element method has the advantage of offering a framework suitable for the modeling of thin shells. In any case, finite element models seem to be the right option if one wishes to have a general and versatile model involving realistic constitutive laws and post-buckling effects.

Our purpose in the present paper is to develop a strategy which combines an integral formulation for the fluid flow and a full finite element membrane model for the capsule wall. However, the integral form of the Stokes equations gives the flow velocity as a function of the force distribution on the domain boundary. Obtaining the force distribution as a function of velocity would require the inversion of a dense matrix. We thus consider that the fluids impose the membrane displacement and that the reaction forces of the membrane govern the fluid flow. This means that the finite element method is not used in its classical formulation and that it is impractical to use commercial or classical finite element programs. In this paper, we propose an explicit finite element formulation to model the large deformation of a capsule and couple it with a boundary integral formulation of the Stokes flow equations.

In Section 2, we detail the formulation of the coupling problem and, in Section 3, the numerical method used to solve it. In Section 4, we assess the non-linear behavior of the finite element method and show the convergence of the coupling procedure. We then study two complex test cases: a capsule in a simple shear flow and in a plane hyperbolic flow. First, we validate our results against previous studies. Second, we use the new numerical method to further the understanding of these two cases, as no consensus was yet reached between previous studies regarding the mechanical behavior of a capsule in these two flows. We determine whether some effects observed in previous studies are physical or byproducts of the numerical methods used. We conclude with a discussion of the advantages of the new coupling procedure and of possible ways to extend it.

PROBLEM STATEMENT

We consider an initially spherical capsule (radius a) consisting of a liquid droplet enclosed by an infinitely thin membrane characterized by a surface shear modulus G s . The capsule is suspended in a viscous liquid undergoing an unbounded shear flow. The frame of reference with Cartesian basis (e 1 , e 2 , e 3 ) is centered on the center of mass of the capsule (Figure 1). The velocity of the unperturbed flow is denoted v ∞ with characteristic shear rate ˙ . Both fluids are assumed to be incompressible, Newtonian and to have the same viscosity and density . Dimensional analysis shows that an important non-dimensional parameter is the capillary number

Ca = ˙ a G s , ( 1
)
which compares the viscous stresses exerted by the fluids to the elastic stiffness of the membrane.

Internal and external flows

The particle Reynolds number Re = a 2 ˙ / is assumed to be small compared to unity, so that both the internal and external flows follow the Stokes equations:

∇ •v = 0, ∇ •r = 0, ( 2
)
where v is the Eulerian velocity field and r is the Cauchy stress tensor in the fluids. The velocity of the points of the capsule can then be related to the tractions on the membrane through an integral equation over the deformed surface S [START_REF] Rallison | A numerical study of the deformation and burst of a viscous drop in extensional flow[END_REF][START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF]: where [r]•n = (r ext -r int )•n is the traction jump across the interface (n is the unit normal vector pointing outwards and the superscripts indicate the internal and external flows), and

∀x ∈ S, v(x) = v ∞ (x)- 1 8 S J(x, y)•[r]•n(y) dS y , (3)
J(x, y) = 1 r 1+ 1 r 3 r⊗r (4)
is Green's single layer kernel, with r = x-y and r = r .

Membrane mechanics

Following Skalak et al. [START_REF] Skalak | Strain energy function of red blood cell membranes[END_REF], we treat the membrane as a purely bidimensional sheet of hyperelastic material without bending moments. This is equivalent to considering a membrane with a finite thickness ϑ and supposing that the strains are constant across the thickness and that the external forces are exerted by the two flows on the midsurface S. Each material point of S is characterized by two curvilinear coordinates ( 1 , 2 ). The position of a given material point is thus given as X( 1 , 2 ) in the reference state and x( 1 , 2 , t) in the deformed state. The displacement U(X, t) = x(X, t)-X is related to the velocity of the membrane through the kinematic condition:

* *t U(X, t) = v(x, t). ( 5
)
Dimensional analysis shows that the ratio of the inertial effects in the membrane to the viscous stresses exerted by the flow is of order membrane inertia viscosity

∼ mb ϑ a Re, ( 6
)
where mb is the density of the membrane material. Thus, with the density ratio of order 1 for an organic membrane suspended in water, a membrane of negligible thickness and a vanishing Reynolds number, the inertia of the membrane can be neglected. The motion of the membrane, therefore, follows the local equilibrium equation:

∇ s •T+q = 0, (7)
where T are the bidimensional Cauchy tensions in the membrane and ∇ s • is the surface divergence operator on the deformed configuration. The load q on the membrane is equal to the traction jump:

q = [r]•n. (8)
Equation ( 7) can be rewritten as a variational problem, where the load q is unknown. Let V be the Sobolev space H 1 (S), and û ∈ V a virtual displacement field. The variational problem corresponding to Equation ( 7) is thus:

S û•qdS = S ê( û) : T(U) dS, ∀ û ∈ V, (9)
where ê( û) = 1 2 (∇ s û+∇ s ûT ). The test function û must be in H 1 (S) to allow for a direct calculation and integration of the strain tensor ê( û). Note that neither the virtual work principle nor the boundary integral (Equation ( 3)) requires that the load q be in H 1 (S). In the virtual work principle, the load must only be in the dual space V . However, in the boundary integral equation, it is convenient to choose q ∈ L 2 (S) allowing for an explicit calculation of the integral. A less regular load would require to integrate Equation (3) by part, which is very delicate as the gradient of the kernel J is highly singular for x = y.

Stress-strain relations

To close the problem, we need the (non-linear) relation between the displacement U and the tensions T. Since we study the large deformation of a curved surface, it is convenient to use the local covariant (a 1 , a 2 , a 3 ) and contravariant (a 1 , a 2 , a 3 ) bases constructed from the curvilinear coordinates 1 and 2 [START_REF] Green | Large Elastic Deformations (2nd edn)[END_REF][START_REF] Crisfield | Non-Linear Finite Element Analysis of Solids and Structures: Advanced Topics[END_REF]. The vector a 3 = a 3 = n is the unit normal vector (pointing outward) and the tangent vectors are defined by a = x , , where . , denotes derivation with respect to . Throughout the paper, Greek indices are equal to 1 or 2, while Latin indices are equal to 1, 2, or 3, and summation is performed over repeated indices. A vector v can thus be written as:

v = v i a i = v i a i = v Xi e i , ( 10
)
where the subscript X indicates Cartesian components. The representations of the metric tensor are defined as a = a •a and a = a •a . The same quantities are defined in the reference state, using capital letters: the local bases are denoted (A 1 , A 2 , A 3 ) and (A 1 , A 2 , A 3 ), and the covariant and contravariant representations of the metric tensor are A and A , respectively. The gradient of the transformation F = a ⊗A is introduced such that dx = F•dX. The Green-Lagrange strain tensor is e = 1 2 (F T •F-1). Invariants of the transformation can then be defined as:

I 1 = tr(F T •F)-2 = A a -2, I 2 = det(F T •F)-1 = |A ||a |-1. (11)
Assuming the membrane to be 2D isotropic, the Cauchy tensions are related to a strain energy function w s (I 1 , I 2 ) per unit undeformed surface area by:

T = 1 J s F• *w s *e •F T , ( 12
)
where the Jacobian J s = det F represents the ratio between the deformed and undeformed surface areas. Using the chain rule, one obtains the following expression for the contravariant representation of T [START_REF] Green | Large Elastic Deformations (2nd edn)[END_REF]:

T = 2 J s *w s *I 1 A +2J s *w s *I 2 a . ( 13
)
Several laws have been proposed to model the behavior of hyperelastic membranes [START_REF] Oden | Finite Elements of Non-Linear Continua[END_REF]. A classical 3D law used to describe volume-incompressible, rubber-like materials is the neo-Hookean (NH) law, with a 2D form given by:

w NH s = G NH s 2 I 1 -1+ 1 I 2 +1 . ( 14
)
Skalak et al. [START_REF] Skalak | Strain energy function of red blood cell membranes[END_REF] proposed a law (Sk) with independent surface shear and area dilation moduli:

w Sk s = G Sk s 4 (I 2 1 +2I 1 -2I 2 +C I 2 2 ), C > -1/2, (15)
where the parameter C is such that the ratio of the area dilation modulus to the surface shear modulus is 1+2C. For large values of C, the Sk law describes an almost area-incompressible material, which is appropriate for biological membranes. It has also been shown to model effectively the behavior of alginate artificial capsules when using C = 0 [START_REF] Carin | Compression of biocompatible liquid-filled HSA-alginate capsules: determination of the membrane mechanical properties[END_REF]. Note that this law was developed directly from the perspective of 2D elasticity and only assumes in-plane isotropy. From a 3D perspective, it could describe a material either fully isotropic or transversely isotropic.

In the domain of small deformation, both laws are equivalent to the 2D formulation of Hooke's law

w H s = G s tr(e 2 )+ s 1-s (tre) 2 , ( 16
)
where e is the 2D linearized Green-Lagrange strain tensor and s ∈]-1, 1[ is the 2D Poisson ratio.

The equivalence between the laws is obtained with

G s = G NH s , s = 1/2, ( 17
)
for the NH law and

G s = G Sk s , s = C 1+C , ( 18
)
for Sk law.

For large deformations, however, the two laws behave very differently, the NH law being strainsoftening whereas the Sk law is strain-hardening [START_REF] Skalak | Strain energy function of red blood cell membranes[END_REF]. Such large deformations can be expected to occur when dealing with thin capsules since the surface shear modulus G s decreases with the membrane thickness.

NUMERICAL METHOD

We present the numerical procedure we use to solve the coupled problem described in the previous section. The undeformed capsule is put in its initial position and the flow is started. At each time step, the load on the membrane is determined from the deformation by a finite element method. From the load, the velocity of the membrane points can then be computed explicitly using a boundary integral method. Finally, the velocity is integrated to obtain the new displacement and deformation. These steps will now be detailed.

Discretization of the interface

The surface of the capsule is discretized using triangular elements. Two types of elements are implemented: flat P 1 elements with three nodes (one at each vertex) and linear shape functions or curved P 2 elements with six nodes (one at each vertex and at the middle of each side) and quadratic shape functions. The total number of elements and nodes are denoted N E and N N , respectively. The characteristic mesh size is denoted x, or h = x/a in non-dimensional form.

The sphere at rest is meshed by first inscribing an icosahedron (regular polyhedron with 20 triangular faces) in the sphere. The mesh is constructed by placing a new node at the middle of each edge and dividing each element into four new elements; the new nodes are then projected onto the sphere. The procedure is repeated until the desired number of elements is attained [START_REF] Ramanujan | Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of capsule viscosity[END_REF]. The mesh is shown in Figure 2 for P 1 elements. The P 2 elements are obtained from the P 1 elements by cutting the edges in half in a last step and projecting the new nodes onto the sphere.

For the purpose of stability and convergence studies, we shall use two mesh size definitions

h 1 = 1 a 4 a 2 N E and h 2 = 1 2a 4 a 2 N E , (19)
for P 1 and P 2 elements, respectively. These definitions ensure that the value of h is the same for both element types when the number of nodes of the discretized interface is the same.

Coupling procedure

The whole problem can be viewed as a first-order differential equation in time (Equation ( 5)) with the displacement field U unknown. It describes the kinematic evolution of the membrane and hence of the finite element mesh. Equation ( 5) can be rewritten as:

*U *t = F(U), (20)
where F is a known function of U given by the fluid equation ( 3), in which the traction jump [r]•n = q is obtained by solving the membrane problem [START_REF] Dodson | Spindles cusps and bifurcation for capsules in Stokes flow[END_REF]. We use an explicit second-order Runge-Kutta method to solve this problem. Thus, time convergence can only be expected when the time step t is such that :

t < 3 2L , ( 21
)
where L is the Lipschitz constant of F [START_REF] Stuart | Dynamical Systems and Numerical Analysis[END_REF].

From Equations (3), [START_REF] Ramanujan | Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of capsule viscosity[END_REF], and (12), we observe that L varies approximately as:

L ∼ x × 1 x × G s x = a x ˙ Ca . ( 22
)
(Note that in Equation (3), J decreases as 1/r and the dominant terms in the integral correspond, therefore, to y close to x. The relevant length scale is then x rather than a). Using the order of magnitude of L found in Equation ( 22), the stability condition ( 21) can be rewritten in nondimensional form as

˙ t<O(hCa). ( 23
)
It follows that an explicit scheme can be used to study the motion of deformable particles in a Stokes flow, because viscous forces are comparable to elastic ones, i.e. Ca = O [START_REF] Kühtreiber | Cell Encapsulation Technology and Therapeutics[END_REF]. Thus the restriction on the time step is governed by the spatial discretization and allows to carry out computations in a reasonable time. Conversely, for very rigid structures, where Ca 1, condition [START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF] leads to very small time steps and unbearably long computation times. This is why explicit schemes are not used in such situations.

Such an explicit coupling strategy, which is only conditionally stable in time, is unusual in fluid-structure interaction. The classical reverse strategy, where the displacement is transferred from the structure to the fluid, involves two time differentiation steps: first to obtain the interface fluid velocity as a function of the structural displacement and second to compute the force as a function of the interface velocity. The present strategy, made simpler by the absence of inertia forces, can split the fluid and the structure problems, because the structure acts directly on the fluid displacement to yield a surface traction without time differentiation.

Finite element procedure

For a given deformed state of the capsule, we first solve the solid problem (Equation ( 9)), where V is discretized as a finite element space V h , using the mesh described above. This finite element discretization, used for the displacement fields U and û, is a convenient choice to approximate the unknown load q ∈ L 2 as well. This choice of a common finite element space also has the advantage of coupling neighboring nodes, which provides regularizing properties when solving linear systems.

For a given element, the shape function associated with node p ( p ∈{1 . . . 3} for P 1 elements, p ∈{1 . . . 6} for P 2 elements) is denoted N ( p) ( 1 , 2 ), and the value at node p of the jth Cartesian

component of v ∈ V h is noted as v ( p) X j .
The discretized problem leads to the following matrix system:

[M]{q} = {R}, (24)
where {q} corresponds to the Cartesian components of the nodal values of the discretized load.

The left-hand side of Equation ( 9) can be discretized element-wise as el û(p)

X j S el N ( p) N (q) dS q (q) X j , ( 25
)
where the integration is performed on the deformed state. The matrix [M], therefore, has the structure of a mass matrix.

The vector {R} corresponds to the right-hand side of Equation ( 9) and depends non-linearly on the displacement U. The tensions T are computed directly using Equation [START_REF] Peskin | The immersed boundary method[END_REF], after computing the metric tensor on the deformed geometry of the elements. The virtual strain tensor ε is related to the covariant representation ûi of û by

ε = 1 2 ( û , + û , -2 i ûi ), ( 26
)
where i = a , •a i are the Christoffel symbols (the curvature tensor b = a , •n is denoted 3 for the sake of brevity). The covariant and Cartesian representations of û are related through:

ûi = a X j i û X j , ( 27
)
where a X j i are the Cartesian components of a i . The virtual strain tensor can, therefore, be related to the nodal Cartesian components of the virtual displacement by

ε = ( 1 2 N ( p) , a X j + 1 2 N ( p) , a X j + N ( p) a X j , -i N ( p) a X j i ) û(p) X j . ( 28
)
To solve the solid problem, we first assemble {R} and [M] on the deformed state. In our numerical procedure, the deformed state is known when solving the structural problem and thus the integration on the deformed state is much easier than when using a more standard approach such as in [START_REF] Wriggers | A fully non-linear axisymmetrical membrane element for rubber-like materials[END_REF]. Surface integration is performed using 6 Hammer points [START_REF] Hammer | Numerical integration over simplexes and cones[END_REF] on the elements. Equation ( 24) is then solved using the sparse solver Pardiso [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF][START_REF] Schenk | On fast factorization pivoting methods for sparse symmetric indefinite systems[END_REF].

Fluid solver

Once the load q is known, the velocity field at the nodes is obtained explicitly from the boundary integral Equation (3), which is discretized on the same mesh as for the solid problem. We use [START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF] Hammer points for the integration. Note that the kernel J(x, y) varies as 1/r and becomes very large when y → x, which can lead to numerical errors. Following [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF], we switch to polar coordinates centered on x when y and x belong to the same element. This change of coordinates introduces a Jacobian which goes to 0 as fast as r , thereby eliminating the singularity in J. In this case, 6 Gauss points are used for integration along each of the polar coordinates.

VALIDATION

Large isotropic deformation

Before coupling the finite element method to the fluid solver, we first test it when the membrane of the capsule undergoes large static deformations. We consider a spherical capsule inflated from radius a to radius a p = (1+ )a by an internal pressure p. The membrane undergoes an isotropic traction characterized by a stretch ratio = 1+ . The isotropic principal tensions can be related to through the membrane strain energy function [START_REF] Barthès-Biesel | Effect of constitutive laws for two dimensional membranes on flow-induced capsule deformation[END_REF]. The analytical solution is:

T = G s 1- 1 6 , ( 29
)
for the NH law and for the Sk law. The analytical expression of the pressure can be deduced using Laplace's law:

T = G s ( 2 -1+C 2 ( 4 -1)), (30)
p = 2T a p = 2T (1+ )a . ( 31
)
In order to check the non-linear behavior of the finite element method, we impose a normal displacement field given by U(X) = X. [START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF] In Figure 3, we plot the computed non-dimensional pressure ap/G s as a function of the inflation factor . The results are shown for both laws, using P 1 and P 2 triangular elements. For between 0 and 0.5, we find a maximum relative error with respect to the analytical solution of 0.3% with P 1 elements and 0.09% with P 2 elements. This shows that both element types behave well when submitted to large isotropic deformation.

Convergence

The first question that arises is whether the coupling procedure is numerically stable and whether time and space convergence is achieved. To verify these points, we consider a classical test case: the motion of an initially spherical capsule in a simple shear flow. We focus on the case of a membrane following the NH law, at a capillary number of Ca = 0.6 and study the influence of the time step and of the mesh size. A full study of the motion of a capsule in a simple shear flow will be performed in Section 4.3.

The undisturbed velocity field is given by:

v ∞ (x) = ˙ x 2 e 1 . ( 33
)
In such a flow field, the capsule elongates along a direction contained in the shear plane. The deformed shape of the capsule is almost ellipsoidal and can be approximated by its ellipsoid of inertia. The deformation is then characterized by the Taylor parameter D ij :

D ij = |L i -L j | L i + L j (i, j = 1, 2, 3), ( 34
)
where L i are the lengths of the principal axes of the ellipsoid of inertia (due to flow symmetry, L 1 and L 2 are in the shear plane). At Ca = 0.6 using the NH law, a steady state can be reached. The coefficient D ∞ 12 represents the Taylor parameter D 12 in the shear plane at steady state. We first study the effect of the dimensionless time step ˙ t on the stability of the numerical procedure with a very fine mesh (N N = 10 242 nodes, mesh size: h = 2.5×10 -2 ). This mesh size corresponds to N E = 20 480 P 1 elements and to N E = 5120 P 2 elements. As predicted by Equation ( 23), the procedure is conditionally stable. The time step stability condition in this particular case is found to be ˙ t 6×10 -3 for P 1 elements and ˙ t 5×10 -3 for P 2 elements. Below this critical time step, a steady state is reached.

The theoretical value of the deformation is not known, as there is no analytical solution to this problem. We use as a reference value for each element type P i (i = 1, 2) the value of D ∞ 12 , denoted D ref i , obtained with this mesh size for a small time step ˙ t = 1×10 -4 . We first increase the time step and compute the error

i = |D ∞ 12 -D ref i |/D ref i
relative to the reference value, for each element type. As shown in Figure 4, the numerical procedure converges linearly in time for both P 1 and P 2 elements. It should be noted that the error remains small ( i <2×10 -4 ) for all values of the time step small enough to satisfy Equation [START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF].

We then study the influence of the mesh size by increasing h with a small time step ˙ t = 1×10 -4 . As shown in Figure 5, the numerical procedure converges with respect to h, approximately as h 4 , and the precision of the results remains good even for cruder meshes ( i <2×10 -2 ). Figure 5 also shows that the relative difference = |D 2 -D 1 |/D 2 between the values obtained with P 1 and P 2 elements decreases when h is reduced. The results given by the two types of elements thus converge toward a common solution.

The numerical procedure is thus found to be conditionally stable and to converge with respect to the time step and mesh size. In view of the convergence results, all following computations are based on P 2 elements with N N = 2562 nodes and N E = 1280 elements (h = 5×10 -2 ). For Ca 0.5, we use a time step of ˙ t = 5×10 -3 ; following Equation [START_REF] Pozrikidis | Boundary Integral and Singularity Methods for Linearized Viscous Flow[END_REF], the time step is reduced proportionally to Ca for lower values of the capillary number.

Simple shear flow

We study the behavior of an initially spherical capsule in a simple shear flow with our new coupling method. In order to validate the method and the non-classical explicit use of finite elements, the results are compared with three previous studies:

• the work of Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF], who used the same method for the fluid problem (boundary elements), but a local approach for the solid problem, based on a discretization of the interface by bi-cubic B-splines; • the work of Li and Sarkar [START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF], who used the immersed boundary method for the flow and a local approach with a piecewise-constant load for the solid problem; • the work of Doddi and Bagchi [START_REF] Doddi | Lateral migration of a capsule in a plane Poiseuille flow in a channel[END_REF], who also used the immersed boundary method for the fluid but finite elements to model the membrane, albeit with a method somewhat different from ours, as noted in the Introduction.

For the external fluid flow, the undisturbed velocity field is given by Equation [START_REF] Schenk | On fast factorization pivoting methods for sparse symmetric indefinite systems[END_REF]. The membrane behavior is simulated using either the NH law or Sk law (C = 1). As the flow is started, the capsule elongates while the membrane rotates along the deformed shape (tank treading motion). During the first stages of deformation, folding occurs in the vicinity of the plane orthogonal to the main extension direction (hereafter called the equatorial plane). This is due to the presence of compression tensions in the membrane. Depending on the value of the capillary number Ca, three steady regimes are found, as first described by Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF]. For Ca<Ca L , the membrane is still subjected to compression in the equatorial area and the folds persist at steady state (Figure 6(a)). For larger values of the capillary number, the capsule becomes more elongated by the shear flow; the value of the isotropic part of the tensions (related to the area dilation modulus) increases, which leads to positive tensions at the equator, and the disappearance of the folds (Figure 6(b)). For Ca>Ca H , the elongated capsule behaves like a slender body submitted to a torque (due to the flow vorticity); negative tensions are again present at steady state and buckling occurs at the tips (Figure 6(c)).

We compare the steady values of D ∞ 12 with the previous studies for the NH law (Figure 7) and the Sk law (Figure 8). We find a good agreement for all values of the capillary number. In particular, our results are always well within 1% of those of Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF] for both the NH and Sk laws, which validates the new method. A 5% difference is, however, found with Li and Sarkar [START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF] for both laws. It is probably related to the crude description of the load used by these authors. The results of Doddi and Bagchi [START_REF] Doddi | Lateral migration of a capsule in a plane Poiseuille flow in a channel[END_REF] are close to our results at low values of the capillary number but seem to diverge as Ca increases. This is somewhat surprising, as they also use finite elements to model the membrane.

The critical values Ca L and Ca H separating the three regimes is an unresolved question, about which no consensus was reached in previous studies. Those critical capillary numbers are defined as the values between which the capsule takes a steady shape, free of negative tensions. The values found with our coupling method are

Ca L = 0.45, Ca H = 0.63, (35)
for the NH law and

Ca L = 0.4, Ca H = 2.4, (36)
for the Sk law (C = 1). They are equal to those found by Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF]. Doddi and Bagchi [20] also found negative tensions at low Ca, but did not quantify the limiting capillary number Ca L . They did not consider large values of Ca and therefore, did not comment on the existence of Ca H . The three regimes were observed by Li and Sarkar [START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF]. They were, however, unable to characterize Li & Sarkar [START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF] Figure 8. Values of D ∞ 12 as a function of Ca for a capsule membrane following the Sk (C = 1) law in a simple shear flow. Results are compared with [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF][START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF]; this case was not studied in [START_REF] Doddi | Lateral migration of a capsule in a plane Poiseuille flow in a channel[END_REF]. Vertical lines indicate the critical capillary numbers Ca L and Ca H .

Ca L and Ca H , as they found negative tensions for all values of Ca, a phenomenon they ascribe to the lack of precision of their numerical method. Hence, our work is the only one beside [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF] that determines Ca L and Ca H ; despite using different numerical techniques, we find the same values as previously.

It should be noted that Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF] report that their numerical method becomes unstable and eventually fails when negative tensions appear: in the compression zones, the distance between nodes tends to become so small that it leads to oscillation of the interpolating B-spline polynomials and failure of the numerical computation. Using our coupling method, we also encounter folds for Ca<Ca L ; however, they do not alter the numerical stability of the computation. Similarly, Li and Sarkar [START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF] reported stability problems for Ca > Ca H . We encounter no such problems with our method, which gives numerically stable results for all the values of Ca that we tried, even far above Ca H (values of Ca up to 1.2 were tested for the NH law and 3.0 for the Sk law). Our method, therefore, presents the advantage of allowing a complete study of the capsule deformation, including at supercritical capillary numbers. We also find that the procedure remains stable even after long computation times. This feature makes it possible to study the tank-treading motion of a capsule, which is a very slow phenomenon: for instance, for a capsule following the NH law at Ca = 0.6, the steady deformed state is reached at ˙ t = 7 but the tank treading period is ˙ T = 22. We were able to run computations over two full tank-treading periods without encountering any stability problem (this corresponds to 8800 iterations).

Planar hyperbolic flow

We now consider the behavior of a capsule in a planar hyperbolic flow. The undisturbed velocity field is given by

v ∞ (x) = ˙ (x 1 e 1 -x 2 e 2 ). (37)
Two numerical studies of the deformation of a capsule in such a flow exist:

• the previously described work of Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF];

• the very recent work of Dodson and Dimitrakopoulos [START_REF] Dodson | Spindles cusps and bifurcation for capsules in Stokes flow[END_REF], who discretize the interface using spectral elements, on which they solve the local equilibrium equation for the membrane and who compute the velocity using the boundary integral technique.

Since Dodson and Dimitrakopoulos

have not yet published the full results of their study, the work of Lac et al. shall be our main reference in this section. We simulate the membrane behavior using the NH and Sk (C = 1) laws. Placed in a planar hyperbolic flow, the capsule elongates in the e 1 direction. As in a simple shear flow, compressive tensions and folds appear in the plane orthogonal to the direction of elongation during the transient phase. For values of the capillary number below a critical value Ca L , these folds remain at steady state. For higher values, a steady state free of compressive tensions can be reached. Similarly to the case of the simple shear flow, Ca L is defined as the value below which negative tensions exist at steady state. As Ca is further increased, we find that the capsule deformation depends on the law used to model the membrane. With the strain-softening NH law, a critical capillary number Ca ∞ exists above which the capsule keeps on extending without reaching a steady state. However, this critical capillary number Ca ∞ does not exist with the strain-hardening Sk law, for which we always find a stationary deformed shape.

We compare the values we obtain for D ∞ 12 to those of Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF] in Figure 9 (NH law) and Figure 10 (Sk law). In the range where they were able to obtain results, we find very good agreement with their values (again well within 1%). Note that the slight discrepancy at low Ca for the NH law is due to the presence of strong folds in the shape that makes it impossible to define D ∞ 12 accurately. As in the case of the simple shear flow, the critical capillary numbers between the three regimes remain an open problem. With our numerical method, the critical values that we find for the NH law are and for the Sk law (C = 1)

Ca L = 0.14, Ca ∞ = 0.35, (38)
Ca L = 0.15, Ca ∞ = +∞. (39)
We find the same values for Ca L as Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF]. However, the results for Ca ∞ : Lac et al. found a lower critical value for the NH law (Ca ∞ = 0.22) and a finite value (Ca ∞ = 0.7) for the Sk law. An explanation for this discrepancy can be found in the work of Dodson and Dimitrakopoulos [START_REF] Dodson | Spindles cusps and bifurcation for capsules in Stokes flow[END_REF], who have shown that capsules in a hyperbolic flow exhibit very high curvatures at the tips during the transient phase, before returning to smaller tip curvatures at steady state. We observe this phenomenon (as shown in Figure 11), although our discretization does not allow a precise quantification of the curvature at the tips. We believe that Lac et al. numerical method became unstable because of the high curvatures and that they misinterpreted the numerical failure as a sign of unbounded extension, which indeed occurs at higher Ca with certain strain-softening materials, but not with strain-hardening ones.

The results found in a simple shear flow and in a hyperbolic flow contribute to showing the accuracy and stability of the new numerical method. They also provide a better understanding of the three regimes that may exist for each flow and of the transition between those.

DISCUSSION AND CONCLUSION

We have developed a new numerical method to simulate the mechanical behavior of a capsule in an unbounded Stokes flow. This method couples a boundary integral method to model the internal and external flows with a non-classical explicit finite element model of the capsule membrane. Coupling those two techniques had never been attempted before and we have shown here that this coupling is feasible, accurate, and stable.

Limitations of the local approach

Writing the membrane equilibrium equation as a variational problem (Equation ( 9)) rather than a local equilibrium equation (Equation ( 7)) imposes different conditions on the discretization of the unknowns. When using the local equilibrium equation, the load q is equal to the surface divergence of the tensions T. Thus to be able to use q in the boundary integral (Equation (3)), T must be continuous over the entire surface of the capsule. The position (and displacement field) must, therefore, be of class C 1 over the surface. Unfortunately, such a C 1 map of the position is not possible on a sphere (this is a classical corollary of the Poincaré-Brouwer 'hairy ball' theorem [START_REF] Milnor | Analytic proofs of the 'Hairy ball theorem' and the Brouwer fixed point theorem[END_REF]). This topological problem is unavoidable when using the local equilibrium equation for a capsule.

All the studies that used the local equilibrium encountered the problem. Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF] meshed the sphere with cubic B-splines, along coordinates similar to standard spherical coordinates. While this method offers a C 2 description of the position on a large part of the capsule, it introduces two singular points-the poles-where the local covariant basis cannot be defined and where the tensions, therefore, cannot be computed. Furthermore, as the Christoffel symbols go to infinity when approaching the poles, numerical errors become important around those points. Lac et al. eventually alleviated this issue by introducing a second mesh, with the poles orthogonal to those of the first mesh [START_REF] Lac | Hydrodynamic interaction between two identical capsules in a simple shear flow[END_REF]. They solved the solid problem once on each mesh and computed a weighted average of the load to reduce the contribution of the areas close to the poles. This method, while giving very good results, is computationally costly as it requires to solve the solid problem twice and to regularly project the position and the load from one grid onto the other. Furthermore, it still exhibits some numerical instability in particular cases like the hyperbolic flow (Equation ( 37)) when one set of poles lies along the elongation axis.

Another technique to mesh the sphere while using a local approach is to use elements between which the first spatial derivatives are discontinuous. An ad hoc method then needs to be implemented to treat the discontinuities. An example of such a method is provided by the work of Wang and Dimitrakopoulos [START_REF] Wang | A three-dimensional spectral boundary element algorithm for interfacial dynamics in Stokes flow[END_REF]. They discretized the interface using a moderate number of spectral elements (between 6 and 14), which offer a highly regular representation of the position within each element, but have discontinuities between the elements. They added a first-order smoothing scheme to enforce the continuity of the position and of its first derivatives from one element to the next. In contrast with the low number of elements (and, therefore, limited amount of smoothing necessary) used in spectral methods, Ramanujan and Pozrikidis [START_REF] Ramanujan | Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of capsule viscosity[END_REF] used an unstructured mesh with 512 triangular P 2 elements, between which the tensions are discontinuous. They developed an intricate technique to average the tensor T at each node from the values on the connected elements and then re-interpolated it between the nodes to compute the load (which they treated as element-wise constant). The effects of so many averages are difficult to predict or control.

In contrast with the local equilibrium, methods relying on the variational problem are not subjected to these topological complications. Equation ( 9) allows u to be in H 1 , which means that the solid problem may be discretized using a finite element space where the first spatial derivatives of the position are discontinuous. The load can thus be obtained without resorting to any sort of uncontrolled smoothing or averaging technique. We believe that variational methods are a sounder way of handling capsule mechanics than methods based on the local equilibrium.

Comparison with other methods

We have compared the results provided by our new model in the case a simple shear flow with three previous studies and have found a good agreement. It may be noted, in particular, that our results are remarkably close to those of Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF] (within 1%). Both studies find the same critical capillary numbers Ca L and Ca H , although two different methods are used to model the behavior of the membrane. This seems to indicate that the values of the critical capillary numbers have a physical relevance and are not a consequence of the numerical method used.

We find results equally close to those of Lac et al. in a planar hyperbolic flow. The same limiting value of the capillary number Ca L is obtained in both studies, once again indicating that this value is physical rather than determined by the numerical method. However, Lac et al. encountered numerical problems at high capillary numbers, due to the presence of high curvatures during the transient phase. The higher numerical stiffness of our elements, compared with the cubic B-splines they used, allows us to remain stable during the transient phase and to reach a steady state when one exists. Our discretization, however, does not allow us to get a precise description of the local curvature at the tips when the capsule is very elongated. In this case, the spectral element method of Dodson and Dimitrakopoulos [START_REF] Dodson | Spindles cusps and bifurcation for capsules in Stokes flow[END_REF] seems better suited. However, with such high curvatures, bending effects should probably not be neglected. Further studies using different numerical methods would be needed to confirm our understanding of the different behaviors of the NH and Sk laws and the values of Ca ∞ that we find.

Bending stiffness and limits of the membrane model

A major difference between our results and those by Lac et al. [START_REF] Lac | Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling[END_REF] is the behavior of the numerical method when the membrane is undergoing in-plane compression. Whereas the technique used by Lac et al. eventually fails when negative tensions appear (Ca < Ca L or Ca > Ca H in a simple shear flow, Ca < Ca L in a planar hyperbolic flow), our numerical method remains stable and a steady equilibrium state is achieved. This is probably due to the larger numerical stiffness of the finite elements as compared to the bi-cubic B-spline functions used by Lac et al. Neither Doddi and Bagchi [START_REF] Doddi | Lateral migration of a capsule in a plane Poiseuille flow in a channel[END_REF] nor Li and Sarkar [START_REF] Li | Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane[END_REF], who also used lower order discretization than Lac et al., reported stability problems at low Ca, which would confirm the hypothesis. The stiffness introduced by the numerical method contributes to the stability of the problem because it enriches the membrane model with some bending stiffness. It allows the numerical procedure to remain stable during transient phases when in-plane compression and high curvatures may render other methods unstable. It also has the advantage that steady states with negative tensions can be computed and studied.

While it stabilizes the numerical procedure, the stiffness introduced by the elements is a byproduct of the numerical method and cannot be controlled or used to model the physical bending stiffness of a capsule. For instance the folds that are observed for Ca<Ca L (Figure 6(a)) depend on the mesh, with a wavelength equal to the element size. Experimentally, the wall of a capsule has a finite thickness and, therefore, a finite bending stiffness. However, it should be noted that, as long as the membrane is locally taut and the local curvature remains small, the bending effects are negligible. For a capsule that is initially spherical and must increase its surface area to deform, these hypotheses are frequently valid. This explains why our method, with its small non-physical stiffness, generally gives results identical to those obtained with Lac et al. method, which has a vanishing bending stiffness.

Wherever negative tensions appear, bending effects determine the behavior of the capsule wall. However, while not appropriate to compute the wrinkled or post-buckling state, a membrane model like ours is sufficient to determine the occurrence and location of negative tensions. A more thorough study of the wrinkles would require treating the capsule wall as a thin shell with a physical bending stiffness. The framework of finite elements seems appropriate to implement such a shell model, which we intend to do in the near future.

Figure 1 .

 1 Figure 1. Capsule freely suspended in a Stokes flow.

Figure 2 .

 2 Figure 2. Typical capsule shape meshed with flat triangles (P 1 elements, N E = 5120, N N = 2562): (a) at rest and (b) after deformation.

Figure 3 .

 3 Figure 3. Non-dimensional pressure as a function of the inflation factor for the NH and Sk (C = 1) laws.N N = 2562 for both types of elements.

Figure 4 .Figure 5 .

 45 Figure 4. Relative error on Taylor's parameter D ∞ 12 as a function of the time step for P 1 and P 2 elements at Ca = 0.6. In all cases, N N = 10 242, h = 2.5×10 -2 (i.e. 20 480 P 1 elements, 5120 P 2 elements). The reference value D ref i for each element type is taken for ˙ t = 1×10 -4 .

Figure 6 .

 6 Figure 6. Steady deformed shape of a capsule, following the NH law, for different values of Ca: Ca = 0.3 (a), 0.6 (b), and 1.2 (c). The color scale corresponds to the normal component of the load, q•n.

Figure 7 .

 7 Figure 7. Values of D ∞ 12 as a function of Ca for a capsule membrane following the NH law in a simple shear flow. Results are compared with [8, 12, 20]. Vertical lines indicate the critical capillary numbers Ca L and Ca H .

Figure 9 .Figure 10 .

 910 Figure 9. Values of D ∞ 12 as a function of Ca for a capsule membrane following the NH law in a planar hyperbolic flow. Results are compared with [8]. Vertical lines indicate the critical capillary numbers Ca L and Ca ∞ .

Figure 11 .

 11 Figure 11. Transient (˙ t = 1.5) and stationary (˙ t = 10) profiles in the (x 1 , x 2 ) plane of a capsule in a hyperbolic flow at Ca = 1.0. The capsule membrane follows the Sk law (C = 1). The curvature at the tips is higher during the transient phase than at steady state.

Table I .

 I Summary of the previous work, classified according to the solid (rows) and fluid (columns) methods used to model the fluid-structure interaction of a capsule in a flow.

		Boundary element Immersed boundary
		method	method
	Local equilibrium
	3D axisymmetric	[6, 15-17]
	3D with piecewise constant loads	[7, 18]	[ 12]
	3D with B-spline discretization	[8, 19]
	3D with spectral discretization	[9]
	Variational method

ACKNOWLEDGEMENTS

This project was supported in part by the Conseil Régional de Picardie through project Flec and by the French Ministère de la Recherche through project Pilcam2. We thank Jean-Louis Batoz (Université de Technologie de Compiègne) for many fruitful discussions at the beginning of the project.

Contract/grant sponsor: Conseil Régional de Picardie through project Flec Contract/grant sponsor: French Ministère de la Recherche through project Pilcam2