Modelling of damage of a liquid-core microcapsule in simple shear flow until rupture

Abstract Capsules, composed of a liquid core protected by a thin deformable membrane, offer high-potential applications in many fields of industry such as bioengineering. One of their limitations comes from the absence of models of capsule damage and/or rupture when they are subjected to an external flow. To assess when rupture is initiated, we develop a fluid–structure interaction (FSI) numerical model of a capsule in Stokes flow that accounts for potential damage of the capsule membrane. We consider the framework of continuum damage mechanics and model the membrane with an isotropic brittle damage model, in which the membrane damage state depends on the history of loading. The FSI problem is solved by coupling the finite element method, to solve for the membrane deformation, with the boundary integral method, to solve for the inner and outer fluid flows. The model is applied to an initially spherical capsule subjected to a simple shear flow. Damage initiates at a critical value of the capillary number, ratio of the fluid viscous forces to the membrane elastic forces and rupture occurs at a higher capillary number, when it reaches a threshold value. The material parameters introduced in the damage model do not influence the mode of damage but only the values of the critical and threshold capillary numbers. When the capillary number is larger than the critical value, damage develops in the two symmetric central regions containing the vorticity axis. It is indeed in these regions that the internal tensions are the highest on the membrane.


Introduction
Capsules consisting of a liquid droplet enclosed by a thin elastic membrane are commonly encountered in nature in the form of red blood cells, fish eggs and vesicles and in numerous industrial processes. The protection and controlled release of active agents is of great importance for diverse applications in the food, cosmetic, bioengineering and medical engineering industries, among others. In medicine, encapsulation has opened the way to new treatment techniques like targeted drug/gene therapy (Bhujbal, de Vosl & Niclou 2014). New-generation bioartificial organs/cells are being developed, for instance, by encapsulating islets of Langerhans to treat diabetic patients (Su et al. 2010) or haemoglobin to create artificial blood (Li, Nickels & Palmer 2005).
But when placed in suspension, capsules are subjected to intense stresses from the surrounding flow, which may cause the mechanical degradation of the membrane. In vivo tests have shown that artificial blood cells could be easily damaged in circulation depending on the particle shape and deformability (Li et al. 2005): this example illustrates the importance of controlling rupture. Depending on the application, capsule damage is to be prevented to preserve the inner substance, or, on the contrary, fostered and directed to allow a targeted release of the encapsulated substance. This necessitates a good understanding of the capsule damage mechanisms under low-inertia flow conditions and of the parameters that control the initiation of rupture.
Very few studies have been conducted on the rupture of capsules subjected to hydrodynamic stresses. The only results that currently exist are experimental. Early experimental studies showed the possibility of wrinkling formation at low shear rates (Walter, Rehage & Leonhard 2001), which could lead to fatigue mechanisms, and of capsule burst at high shear (Chang & Olbricht 1993). The results by Chang and Olbricht (Chang & Olbricht 1993) were obtained on macroscopic spherical capsules, produced through interfacial polymerization. Flow-induced rupture initiated from one of the major axis tips of the deformed ellipsoidal shape of the capsule, which corresponds to the point of minimum thickness. The crack then propagated in the shear plane. Rupture of microcapsules under simple shear flow was observed by Koleva & Rehage (2012) on thin polysiloxane capsules having a high degree of crosslinking. It was reached at small deformations, indicating a brittle behaviour of the capsule membrane. Increasing the shear rate, rupture typically occurred in the central region, close to the tips of the flow vorticity axis, which correspond to the zones of maximum tension. This study corroborated the results by Husmann et al. (2005), who showed that spherical and non-spherical polysiloxane microcapsules burst at the points of maximum elastic tensions, when placed in a spinning-drop apparatus. A similar breakup mechanism has been obtained in confined environments by Abkarian et al. (2008) for red blood cells flowing through a 5 μm wide channel, and by Le Goff et al. (2017) for artificial millimetric capsules and fish eggs trapped at a constriction within a cylindrical channel under a set pressure difference. In both studies, rupture initiated at the front of the capsule, where the tensile tension is the highest. Note that, in Le Goff et al. (2017), rupture could also occur at the point of contact between the capsule and the constriction, but this mode of rupture is different, as it is induced by contact and not by deformation under flow.
What is currently lacking is a model of capsule deformation under flow, capable of assessing when and where the initiation of rupture occurs. The objective of the present study is to develop the first fluid-structure interaction model accounting for membrane damage induced on a liquid-core microcapsule subjected to a simple shear flow. We will use continuum damage mechanics (CDM) to model the initiation and growth of microdiscontinuities (microcavities and microcracks), which lead to the local initiation of macrocracking as they accumulate and coalesce. Contrary to fracture mechanics, which accounts explicitly for the inherent geometrical discontinuity and the associated boundary conditions, the microdiscontinuities are not geometrically modelled in CDM. The local average damage state due to the microdiscontinuities is instead represented by a continuum variable: the damage variable. The CDM has benefited from numerous contributions to its theoretical development (e.g. Kachanov 1986;Lemaitre & Desmorat 2005) since the pioneering work of Kachanov (1958). It is based on the thermodynamics of irreversible processes with internal variables (Coleman & Gurtin 1967), and has been applied to model the damage mechanisms of a large spectrum of materials, from engineering materials (an overview of applications is presented in Lemaitre & Desmorat 2005) to biological tissues (Hokanson & Yazdani 1997;Natali et al. 2005;Holzapfel & Fereidoonnezhad 2017). We propose to incorporate a CDM model into a fluid-structure interaction framework, in order to investigate the time evolution of damage as the capsule deforms under flow.
In this study, we focus on the damage process of a capsule under intense hydrodynamic stresses induced by an external flow over a relatively short time. Due to the short time scale of the phenomena, fatigue or creep damage models are thus not presently relevant. Previous studies have shown that microcapsules may experience ductile (Ghaemi et al. 2016) or brittle (Koleva & Rehage 2012;Le Goff et al. 2017) damage depending on the material and history of loading (external thermo-mechanical stresses). We derive the damage model assuming a quasi-brittle behaviour of the capsule membrane, for which dissipation prior to cracking occurs with negligible irreversible strains (i.e. negligible plasticity). However, CDM provides a general framework: the present model will thus be straightforwardly extended to the other damage behaviours (ductile material, creep or fatigue).
After having detailed the formulation of the damage model of a capsule in infinite shear flow in § 2, we present the model discretization and numerical solver in § 3. We first investigate damage of a spherical capsule under isotropic inflation in § 4, as it provides insight into capsule damage and allows us to validate the numerical method by comparison of the results with the corresponding analytical solution. We then study damage in simple shear flow in § 5, and assess the effect that the dimensionless parameters of the model have on damage evolution and rupture initiation. We finally discuss the model and results in § 6 and analyse the potential of the model to identify the capsule membrane limit of elasticity by comparison with experiments.

Formulation of the problem
We consider a spherical microcapsule of radius a enclosed in an elastic envelope of very small thickness with respect to its radius. The capsule is thus modelled as a two-dimensional incompressible membrane with surface shear elastic modulus G s . It is placed in an infinite shear flow of shear rateγ . The problem is studied in the reference frame of centre O and Cartesian basis (e x , e y , e z ) corresponding to the barycentric reference frame of the capsule oriented such that the unperturbed velocity field is given by v ∞ (x) =γ ze x (figure 1). The inner and outer fluids are the same incompressible Newtonian fluids of dynamic viscosity μ and density ρ. Gravitational and inertial effects being negligible due to the microscopic capsule size, the fluid-structure interaction problem is governed by only one non-dimensional parameter: the capillary number Ca = μγ a/G s , the ratio of the viscous to the elastic characteristic forces. 2.1. Internal and external flows Inertial effects being neglected, the fluid problem is governed by the Stokes equations where σ designates the Cauchy stress tensor, v is the velocity vector and div(·) is the divergence operator. At a given point x of the membrane S, the boundary integral formulation of the Stokes equations gives the relationship between the velocity vector v and the stress tensor σ (Pozrikidis 1992) where n is the unit vector normal to S pointing towards the external fluid and [[σ ]] · n = (σ ext − σ int ) · n is the stress jump across the membrane. We denote as J the second-order Oseen-Burgers tensor defined by where r = x − y, r = r and 1 is the identity tensor.

Wall mechanics
The capsule wall is modelled as a membrane of mid-surface S. The curvilinear coordinates (ξ 1 , ξ 2 ) describe the position x(ξ 1 , ξ 2 , t) on S in the configuration at time t. The position x(ξ 1 , ξ 2 , 0) on the initial configuration S 0 of S is noted X . It is convenient to write the membrane equations in local tangent bases. In what follows, if not specified, indices written with Latin letters take values in {1, 2, 3}, while indices written with Greek letters are in {1, 2}. The covariant basis (a i ) attached to S is defined by The contravariant basis (a i ) is defined by a i · a j = δ j i , where δ j i designates the Kronecker symbol. On S 0 , the covariant and contravariant bases are denoted (A i ) and (A i ), respectively. The metric tensor is g on S and G on S 0 . The contravariant and covariant components of g are a αβ = a α · a β and a αβ = a α · a β , respectively (similar definitions for the components A αβ and A αβ of G).
and Y C being the model constants Table 1. Summary of the key ingredients of the present associated damage model.
The wall inertia being negligible (Walter et al. 2010), the motion of the membrane is governed by the local mechanical equilibrium where q is the surface external load, T is the tension (resultant of the internal Cauchy stress over the thickness), ∇ s · is the surface divergence operator. The dynamic boundary condition imposes that (2.6) The weak form of the membrane equilibrium equation is obtained applying the principle of virtual work For any virtual displacementû ∈ H 1 (S), where H 1 (S) designates the Sobolev space associated with the Lebesgue space L 2 (S) and ε(û) is the symmetric part of g · ∇û, the tensor ∇û being the gradient ofû. In terms of kinematics, the no-slip boundary condition holds on S and gives the relationship between the fluid velocity and the position x of the corresponding point of the membrane v = dx dt . (2.8)

Material behaviour
The model of the capsule wall behaviour is developed in the standard framework of CDM (Lemaitre & Desmorat 2005) to account for the progressive degradation of the membrane while staying in the field of continuum mechanics. More specifically, CDM is a branch of the thermodynamics of irreversible processes with internal variables, the focus of which is to model irreversible transformations associated with damage. The development of a damage model is thus based on four key concepts inherited from the thermodynamics of irreversible processes: state variables, state potential, damage criterion and damage evolution law. A short review of these concepts together with the details of how we developed the model are given hereafter. We specify them in the case of quasi-brittle damage which corresponds to the membrane deformation until the initiation of rupture without irreversible strains (see table 1 for a summary).

State variables
We assume that the transformations of the capsule wall correspond to isothermal elastic deformation and damage. The damage variable represents the irreversible growth of microdefects in a representative volume element (RVE) (figure 2).
To illustrate the definition of the damage variable, we consider a deformed RVE of the capsule wall containing microdefects in the form of microcavities and microcracks (figure 2). We define damage in direction k as the surface ratio δS D /δS, with δS D the maximum intersection of microdefects in a cross-section δS of normal k of the RVE. The stresses on this cross-section are thus transmitted on δS = δS − δS D . We assume that the microdefects have no preferential orientation: the δS D /δS ratio is thus independent of the direction k and corresponds to isotropic damage. The state variable is then the scalar damage variable d defined as (2.9) It ranges from 0, for the local sound (undamaged) state of the material, to 1, when a crack initiates having the size of the RVE. The other state variable is the standard elastic deformation, used in all the mechanical models. The capsule incompressible wall being modelled as a membrane, the in-plane deformation tensor on the mid-surface S is given by the Green-Lagrange strain tensor e: (2.10) The tensor F is the gradient of the transformation of S (2.11) in which, as in what follows, we adopt the convention of summation over repeated indices.
In conclusion, the state variables are d and e, which both depend only on x ∈ S.

State potential
Following the standard framework of CDM, the constitutive law of the membrane and the definition of the variable controlling d are derived from a unique state potential function of the state variables. We note φ(e, d) the specific membrane free energy per unit surface of S 0 . Knowing φ, one can derive the associated variables dual to e and d, using the state laws π = ∂φ ∂e , where π is the second Piola-Kirchhoff tension tensor and Y the specific elastic energy release rate. The Cauchy tension tensor T is related to π through where J is the Jacobian of the transformation of S. The undamaged wall is chosen to follow the neo-Hookean (NH) law, which was shown to model well the elastic behaviour of thin artificial proteic membranes (Chu et al. 2011;Gubspun et al. 2016). The corresponding specific free energy φ NH (Barthès-Biesel, Diaz & Dhenin 2002) is where the two invariants of the transformation I 1 and I 2 are defined by (2.15) What is classical in CDM is to obtain the free energy φ in the damage state using homogenization, which is based on the principle of strain equivalence. We propose to illustrate this concept on the three-dimensional (3-D) RVE shown in figure 2, in the case of a uniaxial traction of intensity δF trac which induces an elongation λ (figure 3).We look for the equivalent RVE (right) having the same cross-section δS, and being subjected to the same elongation λ and loading δF trac as the real RVE. The stress in the equivalent RVE is thus σ = δF trac /δS, which is related to the effective stressσ = δF trac /δS through whereσ is computed from the constitutive law of the undamaged material. The concept of (2.16) can be translated to our 2-D membrane and generalized to any in-plane stress state with (2.17) We thus choose to express the specific free energy φ as Note that the present homogenization process preserves the membrane properties observed in the undamaged case. The state laws defined by (2.12) and (2.13) then have the following expressions: where the Cauchy tension tensor is given through its contravariant components.

Damage criterion and damage evolution law
The last ingredients of the model are the damage criterion and the damage evolution law. We choose to adopt an associated model (Besson et al. 2010), which is numerically robust. It only requires the introduction of the damage threshold function f (Y; d) (d acts as a parameter) to derive the damage criterion and the evolution law through the admissibility condition (i) and the principle of maximum dissipation (ii).
(i) Admissibility condition To be admissible, the associated variable Y must satisfy the standard admissibility condition It defines a bounded domain for Y, illustrated in figure 4(a). Figure 4. Illustration in two dimensions of (a) the admissible domain of the associated variable Y, defined by f (Y) ≤ 0, (b) the case of damage evolution (η > 0) where the yield surface f = 0 moves due to hardening and where the rate of damageḋ is along the normal to the yield surface (normality rule) and (c) the case when damage ceases (η = 0). The thick black lines represent one example of loading cycle, which successively contains all the phases given in table 2: (a) elastic loading 1 , (b) loading with damage 4 , (c) neutral loading 3 followed by elastic unloading 2 + 1 .
(ii) Principle of maximum dissipation The damage evolution is accompanied by dissipation. The associated governing laws are based on the principle of maximum dissipation The solution of the maximization problem under constrain (2.21) is provided by the Kuhn-Tucker conditionsḋ the four of which constitute the evolution law of damage, whereη acts as a Lagrange multiplier. The three conditions within (2.22) 2 are known as the loading/unloading conditions. They provide the damage criterion (2.23) The interior of the admissible domain corresponding to f (Y) < 0 (figure 4a) is thus the elastic domain, in which damage remains constant (ḋ = 0). The domain boundary corresponds to f (Y) = 0 and thus to cases where damage evolves. The damage evolution follows (2.22) 1 which can be interpreted geometrically asḋ being along the normal to the yield surface f = 0 (figure 4b). It is thus referred to as the normality rule. Together, the admissibility condition (2.20) and the damage criterion (2.23) lead to the consistency conditionηḟ = 0.
(2.24) Different cases of loading may exist (see table 2 and figure 4). Whenη = 0, no damage occurs regardless the values of f andḟ . Damage only occurs whenη / = 0, the value of which is obtained by solvingḟ (Y; d) = 0 (imposed by (2.24)). Table 2. Loading case possibilities as a function of the values of f ,ḟ andη. An illustration is given in figure 4 for a loading/unloading cycle.
Note that from the inequality of Clausius-Duhem D ≥ 0, and given that Y ≥ 0, the damage variable d can only grow in timeḋ ≥ 0. (2.25) which restrains the choice of f . Since most artificial and natural microcapsules have been shown to be brittle, we choose to follow the model developed by Marigo (1981) ( 2.27) We presently define κ as a function of two parameters, the damage threshold Y D ≥ 0 and the hardening modulus Y C ≥ 0, such that (2.28) The size of the domain of admissible states f ≤ 0 increases with damage (figure 4b). It is due to the hardening of the material and is controlled by the parameter Y C . The damage evolution law (2.22) can be written equivalently in an explicit form where · + designates the Macaulay brackets defined by (2.31)

Numerical method
Knowing the current position of the material points of the membrane, we perform a Lagrangian tracking of the nodes of the capsule to solve the fluid-structure interaction problem ((2.2), (2.6), (2.7), (2.8) and (2.29)). We use the strategy proposed by Walter et al. (2010) coupling the finite element method to solve for the solid and the boundary integral method to solve for the fluid (figure 5). The problem is solved using the dimensionless (1) Find d n solving (2.29) at the integration points.
(2) Find q n at the nodes solving (2.7) Find the position x n + 1 at the next time step solving (2.8) at the nodes using an explicit scheme x n ← x n + 1 Figure 5. Numerical method to solve the fluid-structure interaction problem over a time step.
forms of the equations, in which the lengths are non-dimensionalized by a, time by 1/γ and tensions by G s . The two parameters Y D and Y C are thus also non-dimensionalized by G s . The originality of our work consists in introducing a damage model in the solid problem. At the material level, the evolution of the damage variable d is determined for each integration point using the explicit equation (2.29). The external load q is then obtained by solving the global problem (2.7) and transferred to the fluid problem. The velocity is computed explicitly at each node from (2.2). Finally, (2.8) is integrated with a second-order explicit Runge-Kutta scheme to solve for the position of the membrane nodes at the next time step.

Mesh
A conformal mesh is used, the nodes on the capsule S being shared by the fluid and the solid problems. The mesh is composed of curved triangular elements containing six nodes with quadratic shape functions (P2 elements). The mesh is generated on the spherical shape corresponding to the initial configuration (figure 6). Following a previous study (Walter et al. 2010), the mesh contains N E = 1280 P2 elements corresponding to a total of N N = 2562 nodes.

Solid solver
For a given deformed configuration of the capsule, the discrete solid problem consists in finding the external load q ∈ L 2 h and the damage d ∈ L 2 h that satisfy (2.7) and (2.29), where the subscript h indicates the finite element space. The position x and the virtual displacementû are searched in H 1 h . Using isoparametric elements, we restrict the solution and v ( p) are the shape function and the nodal coordinates of v associated with the node p, respectively. Noting v ( p) Xj the coordinates of v ( p) in a Cartesian basis (e Xj ), the left-hand side of the discretized form of (2.7) writes where {q} and {û} are the vectors of size 3N N of the nodal coordinates, and χ ( p)Xj αβ is defined by Equation (2.7) being satisfied for any virtual displacement, the discrete solid problem writes Find q and d, such that, The square and column matrices [M] and {R} are, respectively, computed at each time step by using six Hammer points on each element (Hammer, Marlowe & Stroud 1956 3.3. Fluid solver For a given deformed configuration of the capsule and knowing the stress exerted by the membrane on the fluid, the velocity field v is given explicitly by (2.2). The velocity field v is computed at each node. The integral on the right-hand side of (2.2) is computed with 12 Hammer points per element. To handle the singularity of the tensor J at node x, we use polar coordinates centred on x when integrating on the elements sharing this node (for more details see e.g. Lac et al. 2004). We do not use penalty methods to impose the conservation of the volume of the fluids. Still, the maximum relative variation of the capsule volume is limited to 0.1 % of the initial volume.

Coupling
Using a conformal mesh with isoparametric elements, the loads [[σ ]] · n and q are in the same space H 1 h . Hence the dynamic coupling between the fluid and the solid (2.6) is verified in its strong form in this space. Considering the kinematic coupling, the no-slip condition (2.8) is solved at the nodes with an explicit second-order Runge-Kutta scheme to find the position of the nodes at the next temporal increment. Since the local problem of damage is solved in the solid problem with an implicit scheme, the condition of stability of the scheme of temporal integration of the fluid-structure interaction problem is the same as the one initially developed by Walter et al. (2010).

Capsule damage under isotropic inflation
We first analyse the damage of a spherical capsule under osmotic inflation. We impose radial displacements inflating the capsule from radius a to radius a(1 + α(t)), where the inflation ratio α is such that α ≥ 0. We will study two cases: a monotonic increase of α and cyclic variations of α with successive increase and decrease of the capsule diameter. We compare the solution given by the solid solver to the analytical solution.
The problem consists in finding the damage variable d and the external load q that satisfy the evolution law of damage (2.29) and the equilibrium of the membrane (2.7). An analytical solution of the problem exists. We look for it in the form of uniform fields that satisfy the spherical symmetry of the problem. The stretch ratio of the membrane, which is the square root of the isotropic principal value of the dilatation tensor F T · F , is simply λ = 1 + α. The corresponding isotropic principal value T of the tension is and the elastic energy release rate Y As Y increases monotonically with α, the evolution law for damage (2.29) writes where α max is defined similarly to Y max in (2.31). Hence, the condition for d to increase is that α is larger than any of its previous values. The external load is q = pn, where p ≥ 0 is the difference between the internal and external pressures. Choosing test functions of the formû =ûx in the equilibrium equation (4.1)-(4.4), and numerically using the solid solver presented in § 3. For the numerical solution, we compute p and d as surface averages, the pressure difference p being given by q · n. Between the numerical and analytical solutions, we always find relative errors lower than 10 −3 % for the pressure difference p and 10 −4 % for damage d.
We first compare how ap/G s (the dimensionless value of p) and d evolve with the inflation ratio α in the case of a monotonic inflation of the capsule (figure 7). The numerical and analytical curves are perfectly superimposed (figure 7b,c) and comparison with the analytical solution of the undamaged capsule (d = 0) shows a clear effect of damage on the pressure difference (figure 7b). Damage is initiated at the critical value α = α c , which corresponds to Y(α c ) = Y D . The loss in elastic properties of the damaged capsule leads to a reduction in pressure difference as compared to the undamaged case. The pressure difference returns to zero when d = 1, which occurs when α = α .
We then compare the evolution of ap/G s and d with α in the case of a capsule subjected to cyclic inflations and deflations with increasing maximum sizes (figure 8). During the first cycle corresponding to the inflation of the capsule until point A, the value of α does not exceed the critical value α c . Hence damage does not initiate and the curves of ap/G s for the damaged and undamaged capsules coincide during inflation and deflation. For the second cycle (inflation until point B), the curves of d and ap/G s coincide with the corresponding curves obtained for the monotonic size increase (figure 7b,c). During deflation from point B, damage remains constant and the curve of pressure difference ap/G s stays below the inflation curve when α decreases back to 0. For the third cycle, the inflation curves of ap/G s and d overlap the corresponding curves of the previous deflation until point B. But, between points B and C, damage increases during inflation, and the curve of ap/G s again coincides with the corresponding curve obtained for the monotonic size increase. The deflation from point C is then similar to that of the second cycle with constant damage and an ap/G s -curve below the inflation one. During the last inflation, capsule rupture occurs, when α reaches the limit value α (corresponding to d = 1). The case of the capsule under isotropic inflation illustrates the effects of damage on the behaviour of the capsule. For a given value of the inflation ratio α, the more the membrane is damaged, the lower the pressure difference (figure 8b), in other words, damage reduces the loading capacity of the membrane. For increasing d, the slope at the origin for the curve ap/G s (α) decreases (figure 8b), which means that damage reduces the stiffness of the structure. It is interesting to see how the values of α c and α depend on the parameter values Y D and Y C . Following (4.3), the values of α initiating damage and rupture are given respectively by the equations Y(α c ) = Y D and Y(α ) = Y D + Y C , where the expression of Y(α) is obtained using (4.2). The critical inflation ratio α c thus depends solely on Y D , but the limit inflation ratio α depends on both Y D and Y C . Furthermore, the higher the parameter values, the higher the two threshold inflation ratios.

Capsule damage under simple shear flow
We now study the damage of a capsule in simple shear flow. We first show the typical motion and evolution of damage of a capsule in a reference case and then study the influence of the capillary number on the capsule behaviour. We will see that, when the capillary number is increased, three different regimes are found. The capsule is first undamaged until a critical capillary number Ca c is reached, corresponding to the onset of damage. Above this value, the capsule reaches a steady-state deformed shape in which it is partly damaged. When the limit capillary number Ca is reached, rupture initiates putting a limit to the damage regime. In the last part of this section, we will finally study

Coupled kinetics of motion and damage on a reference case (Ca = 0.7)
As reference case, we choose Y D = 0.2, Y C = 2.0 and Ca = 0.7. The value of Ca is such that Ca c < Ca < Ca . Hence the capsule is damaged but the damage stabilizes and a steady state is reached.
Upon the start of the shear flow at t = 0, the initially spherical capsule rotates and takes an ellipsoidal deformed shape. It gets flattened while inclining towards the direction of the flow e x (figure 9). Figure 10 shows the evolution of the capsule state over time until steady state. Note that the membrane rotates around the vorticity axis e y and has a so-called tank-treading motion. We choose to show the capsule shape and damage distribution at different stages: at the onset of damage (t = t c ), at an intermediate instant while damage develops, at maximum elongation (t = t 1 ) and at steady state (t = t ∞ ). The capsule states are shown in the current configuration from two view points in the shear plane and in the transverse inclined plane containing the maximum principal direction e 1 (figure 9). Damage is initiated, at time t c , at the points P and P which are on the vorticity axis (O, e y ). As the capsule elongates, two symmetric disjoint damaged areas form around points P and P , which correspond to the locations of maximum damage d max at each instant. Due to the irreversibility of damage, the maximum values d ∞ max are found at P and P at steady state (t = t ∞ ). In order to see whether preferential direction of damage exists, we plotted the damage distributions on the initial capsule configuration (last row of figure 10). Damage initially develops preferentially along the direction of maximum elongation e 1 but the anisotropy decreases after time t 1 to reach a quasi-isotropic damage distribution at steady state. This may be induced by the tank treading of the capsule membrane around the vorticity axis. Figure 11 gives complementary information on the evolution of the state of the capsule over time until the steady damaged state. The localization of the energy release rate Y, and hence of damage, in the regions around the points P and P (see figure 10) is correlated with the maximum of the principal tension T 1 (figures 11a-d and 11e-h). Damage has no visible consequences on membrane wrinkling: the wrinkles visible on the normal load maps in figure 11(i-l) are the same as in Walter et al. (2010) in the case without damage. They are induced by the presence of negative T 2 tensions transverse to the direction of the wrinkles ( figure 11m-p). The capsule wall being presently modelled as a membrane devoid of bending stiffness, the wrinkle amplitude and wavelength are purely numerical. But the small amplitude of the negative part of T 2 tensions indicates that they hardly contribute to the energy release rate Y and thus to damage. Consequently, they do not lead to any numerical artefact and damage is well predicted by the present model.  We now investigate how the capsule shape and deformation are influenced by damage. In figure 12, we compare the time evolution of geometric parameters to the case of a capsule without damage. Since the shape of the capsule can be approximated by an ellipsoid of inertia, we define the principal lengths L 1 and L 2 of the major and minor axes (directions e 1 and e 2 ) in the shear plane (O, e x , e z ) and L 3 , the length along the vorticity axis e y (see figure 9). The capsule indeed elongates along the directions e 1 and e y (L 1 > L 3 > 2a) and shrinks along the direction of the minor axis (L 2 < 2a) (figure 12a). We quantify the deformation of the capsule with the Taylor parameter D 12 = (L 1 − L 2 )/(L 1 + L 2 ) which measures the distortion of the profile of the ellipsoid in the shear plane (figure 12b). The inclination of the capsule is measured by the angle β between the flow direction e x and the direction of the major axis e 1 . Figure 12(c) represents the temporal evolution of β showing that the inclination angle decreases from the first measurable value near π/4. Figure 12(d) shows the evolution of the global surface expansion ratio λ S = (S − S 0 )/S 0 . Figure 12 globally shows that a steady deformed shape is reached. All the quantities tend towards a plateau value which will be denoted with the symbol ∞ hereafter. It is interesting to notice in figure 12(a) that the onset of damage (t = t c ) is not visible on the L i curves. It is only close to t = t 1 that the curves slightly diverge from the case without damage. But only small differences are observed on the principal lengths L i (figure 12a), D 12 (figure 12b) and β (figure 12c) hereafter. In this reference case, we find that damage has no significant   effects on the motion and deformation of the capsule, suggesting that damage will be very difficult to detect experimentally. The geometrical parameter that is the most affected by damage ends up being the global surface expansion ratio λ S (figure 12d). Nevertheless, the difference at steady state is only of a few per cent.

Effect of Ca
We now study the effect of Ca for the same values of parameters (Y D = 0.2, Y C = 2.0) as in the reference case. The corresponding critical and limit capillary numbers are Ca c = 0.37 and Ca = 0.73. The maximum value of damage at steady state d ∞ max is shown as a function of the capillary number Ca in figure 13. For Ca > Ca c , it increases almost linearly with Ca until Ca ∼ 0.6. Above, d ∞ max increases more rapidly with Ca until d ∞ max = 0.4. Around Ca = Ca , it finally reaches the value of 1 at points P and P very sharply, with a slope close to infinity. It is for this reason that it is classical in damage mechanics to relax the criterion for rupture to d = 0.9 or even d = 0.8. Figure 13 indeed shows that they provide the same value for Ca = Ca .
The inserted images in figure 13 show that the damage maxima always lie at points P and P . They also provide an indication of the extent of the damaged zone for increasing values of Ca. Note that for Ca = 0.8 the damage distribution is given at the instant of initiation of rupture t = t and not at steady state, as it no longer exists. In these cases, the capillary number influences mainly the values of damage in the vicinity of points P and P and marginally the damaged surface.
The capsule deformation and inclination at steady state are compared in figure 14 with the non-damaged case for Ca ≤ Ca . No results are shown above Ca , since no steady deformed shape exists any longer (D 12 diverges to infinity). Despite the large effect of Ca on d ∞ max for Ca c < Ca < Ca , the D ∞ 12 and β ∞ curves initially remain superimposed to the non-damaged case, and it is only close to Ca that small differences become visible. No significant influence of damage is thus found on these global quantities. It is a consequence of the localization of damage around points P and P that occurs in the case of a shear flow. Although the evolution curve of D 12 with Ca does not provide information on when damage is initiated (i.e. on the value of Ca c ), it directly provides the value of Ca , which corresponds to when D 12 diverges to infinity (initiation of breakup).

Effect of Y D and Y C
We finally study the influence of the material parameters Y D and Y C on the damage of the capsule. The evolution of d ∞ max (Ca) is represented for different values of Y D and Y C in figure 15. We observe the same trend as in the reference case ( figure 13).
For a fixed value of Y C , Ca c and Ca increase with Y D (figure 15a). However, when Y D is fixed (figure 15b), increasing Y C does not impact when damage initiates (constant Ca c ) but delays when the capsule breaks up (increasing value of Ca ). This relates to the facts that the criterion of initiation of damage is only a function of Y D , whereas the criterion of initiation of rupture is controlled by Y D + Y C , as already shown at the end of § 4.
The results are synthesized in figure 16, which provides a phase diagram of the capsule states for a range of values of Y D and Y C . For a given Y C , the curves Ca c (Y D ) and Ca (Y D ; Y C ) delimit three domains in the parametric space (Ca, Y D ): undamaged for Ca < Ca c , damaged for Ca c < Ca < Ca , ruptured for Ca > Ca . The only effect of Y C is to shift the Ca delimiting curve to higher Ca values as the capsule is then more resistant. This is what is shown by the dotted lines in figure 16, which complete the base case (Y C = 2.0).

Discussion and conclusion
In response to the current need for a damage model of microcapsules in flow, we have developed the first FSI numerical model accounting for the degradation of the capsule membrane until the onset of rupture, when it is deformed by hydrodynamic forces. We have placed ourselves within the framework of continuum damage mechanics, and simulated microdefect development by degrading the elastic material parameters through the introduction of a damage variable d. We have used an isotropic brittle damage model, in which the damage evolution of the membrane depends on the history of loading. We have integrated it in a finite element method that solves for the membrane deformation, which we have coupled to a boundary integral method to solve for the Stokes flows inside and outside the capsule.
6.1. Interpretation of the damage evolution law We have explained the physical meaning of the damage model in § 2, but propose to further detail the interpretation of the damage evolution law (2.29). The capsule membrane being assumed to have a quasi-brittle behaviour, damage evolution is driven by Y max . As an illustration, we propose to introduce a toy model (figure 17), consisting of a bundle of parallel elastic fibres under uniaxial traction (Krajcinovic 1989;Lemaitre & Desmorat 2005). The RVE consists of N parallel elastic fibres initially unbroken and subjected to an elongation ratio λ. Each fibre is associated with the specific elastic energy φ NH and has a brittle behaviour given by the classical energetic criterion of rupture where φ u is a specific energy at rupture and λ max is defined similarly to Y max in (2.31).
The key ingredient of this model is to consider φ u as a random variable with probability density p(φ u ) given by the following band-limited and uniform probability density: where Y D and Y C are the parameters of the damage model introduced in (2.28). Hence, from (6.1) and (6.2), the probability of rupture of a fibre is given by Consistent with (2.9), the damage variable d corresponds to the ratio n b /N, where n b is the number of broken fibres. For a large number of fibres, we can postulate n b = P f (λ max )N, and thus d = P f (λ max ). From (6.3), we obtain where φ NH (λ max ) = Y max (see (2.19) 2 ). We thus retrieve the damage evolution law (2.29). This toy model thus shows that the damage evolution law (2.29) is dictated by local phenomena: each fibre has a binary state broken/unbroken (6.1), for which the transition is randomly triggered. By integrating the function of rupture probability over all the fibres, we obtain a deterministic damage model for the RVE, where d ranges from 0 (all the fibres are unbroken) to 1 (all the fibres are broken). Equation (6.2) shows that the model parameters Y D and Y C delimit the range of dispersion of the specific energy at rupture in the microstructure.

Capsule inflation test
We have first applied the model to a capsule under isotropic inflation, a case for which we derived an analytical solution. This has allowed us to validate the numerical simulations and to show the consequences of damage on the pressure difference p between the internal and external fluids. The main findings are that a given capsule remains sound up to a critical value of the inflation ratio α c , at which damage initiates. As the capsule further inflates above this critical value, the isotropic tension first increases with the isotropic strain, reaches a maximum and then decreases: this corresponds to what is generally defined as a softening behaviour. As damage builds up, the pressure difference decreases, as the global stiffness of the capsule is proportional to the local effective surface shear modulus (1 − d)G s . The pressure difference finally returns to p = 0, which occurs when the damage variable reaches d = 1: it corresponds to the moment when the membrane ruptures. A given capsule is thus characterized by a limit inflation ratio α at which it breaks up.
The inflation capsule test has shown how excellent the agreement is between the theoretical solution and the one obtained with the FSI damage model. If the problem had been solved in displacement (imposed pressure) as classically done in finite element numerical codes, the material softening behaviour resulting from damage would have induced a loss of stability of the uniform solution at the beginning of the regime of strain localization (Rice 1976;Benallal, Billardon & Geymonat 1993). In this regime, a small perturbation from the uniform solution would have localized damage and strain in a band of width of one element: the solution would have been strongly mesh dependent. To solve this issue, classical solid solvers require additional methods, called localization limiters, to obtain more objective solutions (Bažant & Pijaudier-Cabot 1988;Simo, Oliver & Armero 1993). However, it is interesting to notice that even for the case of a capsule in simple shear flow discussed below, where the solution is non-uniform, we did not observe the effect of strain localization by changing the mesh size (results not shown). This shows how advantageous it is to implement the damage model within an explicit FSI solver, where the node displacements are imposed by the fluid and the corresponding external loads exerted by the fluids on the membrane are solved for in the solid problem. Furthermore, the present FSI scheme is particularly robust and stable, thanks to the fact that the quantities are integrated over the surface in both the fluid and solid solvers.

Capsule under simple shear flow
We have then considered a capsule under simple shear flow and similarly seen that there exists a critical value of the capillary number Ca c , at which damage initiates, and a limit capillary number Ca , at which capsule rupture occurs. In the model, we have chosen to base the criterion for damage on the elastic energy release rate Y of the membrane and to use the evolution law developed by Marigo (1981) for quasi-brittle materials, in which damage evolves when Y = Y D + Y C d. The initiation of damage is then solely dictated by the threshold modulus Y D , to which Ca c is proportional. As for the hardening modulus Y C , it governs the rate at which damage occurs: the lower Y C , the faster rupture occurs. The initiation of rupture (d = 1) and the corresponding limit capillary number Ca are thus controlled by Y C + Y D .
For Ca c < Ca < Ca , irreversible damage appears on the flanks of the capsule at the points P and P located on the flow vorticity axis: it is at these locations that the internal membrane tension is the highest. As the capsule tank treads, the two damaged zones grow around these points, but remain confined in their vicinity, the maximum values remaining at P and P . The most striking results in this range of capillary numbers are that the capsule still reaches a steady deformed shape like in the case without damage, and that the effect of damage remains non-visible on the capsule deformed shape, inclination and dynamics. Indeed, damage concentrates around the capsule poles P and P in the case of a simple shear flow, without propagating to the entire capsule. Note that such would not be the case under other flows conditions with three-dimensional vorticity effects, as the capsule rotation would lead to an isotropic distribution of damage all over the capsule membrane. Still, at present, differences in shape and inclination with the no-damage case start to be visible, when Ca gets close to Ca . At Ca = Ca , rupture finally occurs at points P and P , and no steady deformed shape exists thereafter for the capsule.

Comparison with experiments of capsule damage
Damage models are phenomenological and require confrontation with experimental data to assess the relevance of choice of damage evolution law. It is interesting to observe that the present findings corroborate well the results of the few experimental studies present in the literature, which showed that rupture is initiated at the points of maximum elastic tension (Husmann et al. 2005;Abkarian et al. 2008;Koleva & Rehage 2012;Unverfehrt, Koleva & Rehage 2015;Le Goff et al. 2017). The damage model assumptions are thus relevant to study the dynamics of microcapsules in flow.
Comparing the results of the model with experiments also serves the purpose of identifying the values of the model parameters, namely Y D and Y C in the present model. We propose to look more closely at the results obtained by Professor H. Rehage's group on thin polysiloxane microcapsules subjected to a simple shear flow until breakup in a counter-rotating rheometer cell (Koleva & Rehage 2012;Unverfehrt et al. 2015). They followed a given capsule under increasing values of shear rate and found that wrinkles form on the capsule membrane (figure 18b) similarly to what was predicted by numerical models (Lac et al. 2004;Walter et al. 2010). Polysiloxane being very brittle and resistant to deformation, only a small increase in capsule deformation was observed as Ca increased (figure 18a), and rupture occurred at only 3 % of deformation. The crack formed in the region near the vorticity axis (figure 18c), in agreement with the prediction given by our model. Similarly to what we have shown in § 5.2, no influence of damage effects could be observed on the Taylor parameter curve (figure 18a). But, even though simple shear flow experiments do not allow us to identify the value of Ca c (and thus Y D ), Ca is easily identified from the point of divergence of the Taylor parameter curve. Note that in Koleva & Rehage (2012) the capillary number is based on the surface Young's modulus instead of the surface shear modulus as in the present study. However, for the capsules of figure 18, the authors estimated that the two moduli had practically the same values, indicating that the Poisson ratio of the membrane was negative. The polysiloxane capsules of figure 18 are thus found to have a limit capillary number Ca (Y D , Y C ) = 0.01, which provides an implicit relationship between Y D and Y C . Since we know from figure 16 that the damaged domain is delimited by the curves Ca = Ca c and Ca = 0.01, we deduce that Y D ∈ [0; 1 × 10 −3 ] and Y C ∈ [0; 4.6 × 10 −3 ]. We have run simulations assuming Y D /G s = 5 × 10 −4 and Y C /G s = 3.1 × 10 −3 , for which Ca (Y D , Y C ) = 0.01, and found a good fit between the numerical predictions (figure 18e-g) and the experimental results (figure 18b-d). For Ca ≥ Ca , we have continued the simulations after the critical instant t = t where rupture initiates, and found that the totally damaged state d = 1 of the membrane propagates in the plane perpendicular to the major axis (O, e 1 ) and that the capsule elongates indefinitely along its major axis (figure 18g). The divergence of the capsule shape in the simulations (figure 18g) is similar to what is observed experimentally (figure 18d). In retrospect, it is surprising that the experiments by Chang & Olbricht (1993) did not fit those by Professor Rehage's group. Chang & Olbricht (1993), who were the first to study the rupture of polyamide capsules using a counter-rotating rheometer cell, found that rupture initiated at the apex of the major axis, where the capsule is the thinnest. Although these results contradict what all the other studies of the literature have found, it could be interesting to use the damage FSI model to investigate for which damage threshold function (2.27) the model would predict an initiation of rupture at that location.
This study, based on an associated damage model with three ingredients (table 1), could be generalized to include other dissipative phenomena, such as irreversible strains. These have for instance been taken into account by Ghaemi et al. (2016), in the case of a capsule under compression. The model, however, does not include the gradual degradation of the membrane and information on rupture is obtained by post-processing the stress-strain results. The modularity of the framework that we are proposing represents a real advantage if one wants to generalize the use of the damage FSI model for crack nucleation prediction and damage property identification. Predicting crack propagation is, however, outside the scope of the model, as it would require the use of another approach. The extended finite element method (Sukumar et al. 2000;Moës & Belytschko 2002) could then be one option among others to provide answers on the subsequent events following crack nucleation.