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

Nicolas Grandmaison, Delphine Brancherie, Anne-Virginie Salsac

► To cite this version:

Nicolas Grandmaison, Delphine Brancherie, Anne-Virginie Salsac. Modelling of damage of a liquid-core microcapsule in simple shear flow until rupture. Journal of Fluid Mechanics, 2021, 914, 10.1017/jfm.2020.652. http://dx.doi.org/10.1017/jfm.2020.652

HAL Id: hal-03274554

https://hal.science/hal-03274554

Submitted on 30 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modelling of damage of a liquid-core microcapsule in simple shear flow until rupture

Nicolas Grandmaison<fnref n="1" />, Delphine Brancherie<fnref n="2" /> and Anne-Virginie Salsac<fnref n="1" />†

(Received xx; revised xx; accepted xx)

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 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.

## 1. 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 et al. 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 hemoglobin to create artificial blood (Li et al. 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

<footnote n="1">¹Biomechanics and Bioengineering Laboratory (UMR 7338), Université de technologie de Compiègne – CNRS, CS 60319, 60203 Compiègne cedex, France.</footnote>

<footnote n="2">²Roberval Laboratory (FRE 2012), Université de technologie de Compiègne – CNRS, CS 60319, 60203 Compiègne cedex, France.</footnote>

<footnote>† Email address for correspondence: a.salsac@utc.fr</footnote>
depending on the particle shape and deformability (Li et al. 2005): this example illustrates the importance to control rupture. Depending on the applications, 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 to gain 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 et al. 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 bursted 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. 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 (Holzapfel & Fereidoonnezhad 2017; Hokanson & Yazdani 1997; Natali et al. 2005). 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.
<figure><img src="image_2.png" /><figcaption>Figure 1: Capsule suspended in the unbounded simple shear flow.</figcaption></figure>

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 Section 2, we present the model discretization and numerical solver in Section 3. We first investigate damage of a spherical capsule under isotropic inflation in Section 4, as it provides insight on 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 section 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 Section 6 and analyze the potential of the model to identify the capsule membrane limit of elasticity by comparison with experiments.

## 2. 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 $\dot{\gamma}$. The problem is studied in the reference frame of center $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^\infty(\underline{x}) = \dot{\gamma}z e_x$ (Figure 1). The inner and outer fluids are the same incompressible Newtonian fluids of dynamic viscosity $\mu$ and density $\rho$. 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 = \mu\dot{\gamma}a/G_s$, ratio of the viscous to the elastic characteristic forces.
## Internal and external flows

Inertial effects being neglected, the fluid problem is governed by the Stokes equations:

$$ \mathrm{div}(\underline{\sigma}) = 0, \quad \mathrm{div}(\underline{v}) = 0, \qquad (2.1) $$

where $\underline{\sigma}$ designates the Cauchy stress tensor, $\underline{v}$ is the velocity vector and $\mathrm{div}(\cdot)$ is the divergence operator. At a given point $\underline{x}$ of the membrane $S$, the boundary integral formulation of the Stokes equations gives the relationship between the velocity vector $\underline{v}$ and the stress tensor $\underline{\sigma}$ (Pozrikidis 1992):

$$ \forall \underline{x} \in S, \quad \underline{v}(\underline{x}) = \underline{v}^{\infty}(\underline{x}) - \frac{1}{8\pi\mu} \int_S \underline{J}(\underline{x}, \underline{y}) \cdot [\underline{\sigma}] \cdot \underline{n}(\underline{y}) \, dS_{\underline{y}}, \qquad (2.2) $$

where $\underline{n}$ is the unit vector normal to $S$ pointing towards the external fluid and $[\underline{\sigma}] \cdot \underline{n} = (\underline{\sigma}_{\mathrm{ext}} - \underline{\sigma}_{\mathrm{int}}) \cdot \underline{n}$ is the stress jump across the membrane. We denote as $\underline{J}$ the second order Oseen-Burgers tensor defined by:

$$ \underline{J}(\underline{x}, \underline{y}) = \frac{1}{r} \underline{1} + \frac{1}{r^3} \underline{r} \otimes \underline{r}, \qquad (2.3) $$

where $\underline{r} = \underline{x} - \underline{y}$, $r = \|\underline{r}\|$ and $\underline{1}$ is the identity tensor.

## Wall mechanics

The capsule wall is modelled as a membrane of mid-surface $S$. The curvilinear coordinates $(\xi_1, \xi_2)$ describe the position $\underline{x}(\xi_1, \xi_2, t)$ on $S$ in the configuration at time $t$. The position $\underline{x}(\xi_1, \xi_2, 0)$ on the initial configuration $S_0$ of $S$ is noted $\underline{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 $(\underline{a}_i)$ attached to $S$ is defined by:

$$ \underline{a}_{\alpha} = \frac{\partial \underline{x}}{\partial \xi_{\alpha}}, \quad \underline{a}_3 = \frac{\underline{a}_1 \times \underline{a}_2}{\|\underline{a}_1 \times \underline{a}_2\|}. \qquad (2.4) $$

The contravariant basis $(\underline{a}^i)$ is defined by $\underline{a}_i \cdot \underline{a}^j = \delta^j_i$, where $\delta^j_i$ designates the Kronecker symbol. On $S_0$, the covariant and contravariant bases are denoted $(\underline{A}_i)$ and $(\underline{A}^i)$, respectively. The metric tensor is $\underline{g}$ on $S$ and $\underline{G}$ on $S_0$. The contravariant and covariant components of $\underline{g}$ are $a^{\alpha\beta} = \underline{a}^{\alpha} \cdot \underline{a}^{\beta}$ and $a_{\alpha\beta} = \underline{a}_{\alpha} \cdot \underline{a}_{\beta}$, respectively (similar definitions for the components $A^{\alpha\beta}$ and $A_{\alpha\beta}$ of $\underline{G}$).

The wall inertia being negligible (Walter et al. 2010), the motion of the membrane is governed by the local mechanical equilibrium:

$$ \forall \underline{x} \in S, \quad \underline{\nabla}_s \cdot \underline{T} + \underline{q} = \underline{0}, \qquad (2.5) $$

where $\underline{q}$ is the surface external load, $\underline{T}$ is the tension (resultant of the internal Cauchy stress over the thickness), $\underline{\nabla}_s$ is the surface divergence operator. The dynamic boundary condition imposes that:

$$ \forall \underline{x} \in S, \quad \underline{q} = [\underline{\sigma}] \cdot \underline{n}. \qquad (2.6) $$

The weak form of the membrane equilibrium equation is obtained applying the principle of virtual work:

For any virtual displacement $\hat{\underline{u}} \in H^1(S)$,

$$ \int_S \hat{\underline{u}} \cdot \underline{q} \, dS = \int_S \underline{T} : \underline{\varepsilon}(\hat{\underline{u}}) \, dS, \qquad (2.7) $$
<table><caption>Table 1: Summary of the key ingredients of the present associated damage model.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>State / Associated variables</td><td>Green Lagrange strain <u>&epsilon;</u> / Second Piola-Kirchhoff tension <u>&pi;</u><br>Damage variable <i>d</i> / Energy release rate <i>Y</i></td></tr><tr><td>State potential</td><td>Free energy &phi;(<u>&epsilon;</u>, <i>d</i>) = (1 &minus; <i>d</i>)&phi;<sub>NH</sub>(&u;&epsilon;)<br>with &phi;<sub>NH</sub> the isochore neo-Hookean potential</td></tr><tr><td>Damage threshold function</td><td><i>f</i> = <i>Y</i> &minus; &kappa;(<i>d</i>)<br>with &kappa;(<i>d</i>) = <i>Y</i><sub>D</sub> + <i>Y</i><sub>C</sub><i>d</i>,<br><i>Y</i><sub>D</sub> and <i>Y</i><sub>C</sub> being the model constants</td></tr></tbody></table>

where $H^1(S)$ designates the Sobolev space associated to the Lebesgue space $L^2(S)$ and $\varepsilon(\hat{u})$ is the symmetric part of $g \cdot \nabla \hat{u}$, the tensor $\nabla \hat{u}$ being the gradient of $\hat{u}$.

In terms of kinematics, the no-slip boundary condition holds on $S$ and gives the relationship between the fluid velocity and the position $\mathbf{x}$ of the corresponding point of the membrane:

$$ \underline{v} = \frac{d\underline{x}}{dt}. \qquad (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 $\underline{k}$ as the surface ratio $\delta S_D/\delta S$, with $\delta S_D$ the maximum intersection of microdefects in a cross-section $\delta S$ of normal $\underline{k}$ of the RVE. The stresses on this cross-section are thus transmitted on $\delta \tilde{S} = \delta S - \delta S_D$. We assume that the microdefects have no preferential orientation: the $\delta S_D/\delta S$ ratio is thus independent of the direction $\underline{k}$ and corresponds to isotropic damage. The state variable is then the
<figure><img src="image_3.png" /><figcaption>Figure 2: Representation of a microcapsule of mid-surface S placed in an infinite shear flow (left). Zoom on a representative volume element (RVE) containing microcavities and microcracks (middle). Decomposition of the cross-section $\delta S$ of normal vector $\underline{k}$ into the effective load-bearing cross-section $\delta \tilde{S}$ and the total surface of the microdefects $\delta S_D$ (right).</figcaption></figure>

scalar damage variable $d$ defined as

$$d = \frac{\delta S_D}{\delta S} = 1 - \frac{\delta \tilde{S}}{\delta S} \qquad (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 $\underline{e}$

$$\underline{e} = \frac{1}{2}(\underline{F}^T \cdot \underline{F} - \underline{G}) \qquad (2.10)$$

The tensor $\underline{F}$ is the gradient of the transformation of $S$

$$\underline{F} = \frac{\partial \underline{x}}{\partial \underline{X}} = \underline{a}^{\alpha} \otimes \underline{A}_{\alpha} \qquad (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 $\underline{e}$, which both depend only on $\underline{x} \in 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 $\phi(\underline{e}, d)$ the specific membrane free energy per unit surface of $S_0$. Knowing $\phi$, one can derive the associated variables dual to $\underline{e}$ and $d$, using the state laws

$$\left\{ \begin{array}{l} \underline{\pi} = \frac{\partial \phi}{\partial \underline{e}}, \\ Y = -\frac{\partial \phi}{\partial d}, \end{array} \right. \qquad (2.12)$$
<figure><img src="image_4.png" /><figcaption>Figure 3: Illustration of the homogenization principle on a representative volume element (RVE) under an uniaxial traction of intensity $\delta F_{trac}$ and of the associated elongation $\lambda$. The heterogeneous damaged material (real RVE) is modelled as a homogeneous domain (equivalent RVE) with the same cross-section $\delta S$ and subjected to the same loading/elongation. The force equilibrium leads to $\sigma = \delta \tilde{S}/\delta S \tilde{\sigma} = (1 - d) \tilde{\sigma}$, where the effective stress $\tilde{\sigma}$ is the stress transmitted through the load-bearing cross-section $\delta \tilde{S}$ and determined with the constitutive law of the undamaged material.</figcaption></figure>

where $\underline{\underline{\pi}}$ is the second Piola-Kirchhoff tension tensor and $Y$ the specific elastic energy release rate. The Cauchy tension tensor $\underline{T}$ is related to $\underline{\underline{\pi}}$ through

$$ \underline{T} = \frac{1}{J} \underline{\underline{F}} \cdot \underline{\underline{\pi}} \cdot \underline{\underline{F}}^T. \qquad (2.13) $$

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 $\phi_{NH}$ (Barthès-Biesel et al. 2002) is

$$ \phi_{NH}(\underline{\underline{e}}) = \frac{G_s}{2} \left( I_1 - 1 + \frac{1}{I_2 + 1} \right), \qquad (2.14) $$

where the two invariants of the transformation $I_1$ and $I_2$ are defined by:

$$ \begin{aligned} I_1 &= \mathrm{tr}(\underline{\underline{F}}^T \cdot \underline{\underline{F}}) - 2 = A^{\alpha\beta} a_{\alpha\beta} - 2, \\ I_2 &= \mathrm{det}(\underline{\underline{F}}^T \cdot \underline{\underline{F}}) - 1 = \mathrm{det}(A^{\alpha\beta}) \mathrm{det}(a_{\alpha\beta}) - 1. \end{aligned} \qquad (2.15) $$

What is classical in CDM is to obtain the free energy $\phi$ in the damage state using homogenization, which is based on the principle of strain equivalence. We propose to illustrate this concept on the 3D RVE shown in Figure 2, in the case of a uniaxial traction of intensity $\delta F_{trac}$ which induces an elongation $\lambda$ (Figure 3). We look for the
equivalent RVE (right) having the same cross-section $\delta S$, and being subjected to the same elongation $\lambda$ and loading $\delta F_{trac}$ as the real RVE. The stress in the equivalent RVE is thus $\sigma = \delta F_{trac}/\delta S$, which is related to the effective stress $\tilde{\sigma} = \delta F_{trac}/\delta \tilde{S}$ through:

$$ \sigma(\lambda, d) = \frac{\delta \tilde{S}}{\delta S} \tilde{\sigma}(\lambda) = (1 - d) \tilde{\sigma}(\lambda). \qquad (2.16) $$

where $\tilde{\sigma}$ is computed from the constitutive law of the undamaged material.

The concept of equation (2.16) can be translated to our 2D membrane and generalized to any in-plane stress state with:

$$ \frac{\partial \phi}{\partial \underline{\epsilon}}(\underline{\epsilon}, d) = (1-d) \frac{\partial \phi_{NH}}{\partial \underline{\epsilon}}(\underline{\epsilon}). \qquad (2.17) $$

We thus choose to express the specific free energy $\phi$ as:

$$ \phi(I_1, I_2, d) = (1 - d)\phi_{NH}(I_1, I_2). \qquad (2.18) $$

Note that the present homogenization process preserves the membrane properties observed in the undamaged case. The state laws defined by equations (2.12) and (2.13) then have the following expressions:

$$ \left\{ \begin{array}{l} T^{\alpha\beta} = (1-d)G_s \left( \frac{1}{J} A^{\alpha\beta} - \frac{1}{J^3} a^{\alpha\beta} \right) \\ Y = \phi_{NH}, \end{array} \right. \qquad (2.19) $$

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

$$ f(Y; d) \le 0. \qquad (2.20) $$

It defines a bounded domain for $Y$, illustrated in Figure 4a.

### (ii) Principle of maximum dissipation

The damage evolution is accompanied by dissipation. The associated governing laws are based on the principle of maximum dissipation

$$ \mathcal{D}(Y, \dot{d}) = \max_{f(Y^*; d) \le 0} \{\mathcal{D}(Y^*, \dot{d})\}, \qquad (2.21) $$

where $\mathcal{D}(Y, \dot{d}) = Y \dot{d}$, and $\dot{d}$ designates the material temporal derivative of $d$.

The solution of the maximization problem under constrain (2.21) is provided by the Kuhn-Tucker conditions

$$ \left\{ \begin{array}{l} \dot{d} = \dot{\eta} \frac{\partial f}{\partial Y} \\ f \le 0, \dot{\eta} \ge 0, \dot{\eta} f = 0. \end{array} \right. \qquad (2.22) $$
<table><caption>Table 2: Loading case possibilities as a function of the values of $f$, $\dot{f}$ and $\dot{\eta}$. An illustration is given in Figure 4 for a loading/unloading cycle.</caption><tr><td>① elastic loading/unloading</td><td>$f < 0$</td><td>$\dot{\eta} = 0$</td></tr><tr><td>② elastic unloading</td><td>$f = 0$</td><td>$\dot{f} < 0$</td><td>$\dot{\eta} = 0$</td></tr><tr><td>③ neutral loading</td><td>$f = 0$</td><td>$\dot{f} = 0$</td><td>$\dot{\eta} = 0$</td></tr><tr><td>④ loading with damage</td><td>$f = 0$</td><td>$\dot{f} = 0$</td><td>$\dot{\eta} > 0$</td></tr></table>

<figure><img src="image_5.png" /><figcaption>Figure 4: Illustration in 2D of (a) the admissible domain of the associated variable $Y$, defined by $f(Y) \le 0$, (b) the case of damage evolution ($\dot{\eta} > 0$) where the yield surface $f = 0$ moves due to hardening and where the rate of damage $\dot{d}$ is along the normal to the yield surface (normality rule), and (c) the case when damage ceases ($\dot{\eta} = 0$). The thick black lines represent one example of loading cycle, which successively contains all the phases given in table 2: (a) elastic loading ①, (b) loading with damage ④, (c) neutral loading ③ followed by elastic unloading ② + ①.</figcaption></figure>

the four of which constitute the evolution law of damage, where $\dot{\eta}$ acts as a Lagrange multiplier.

The three conditions within equation (2.22)$_2$ are known as the loading/unloading conditions. They provide the damage criterion

$$ \begin{cases} f(Y) < 0 \Rightarrow \dot{\eta} = 0 \\ f(Y) = 0 \Rightarrow \dot{\eta} \ge 0. \end{cases} \qquad (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 ($\dot{d} = 0$). The domain boundary corresponds to $f(Y) = 0$ and thus to cases where damage evolves. The damage evolution follows equation (2.22)$_1$ which can be interpreted geometrically as $\dot{d}$ 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

$$ \dot{\eta} \dot{f} = 0. \qquad (2.24) $$

Different cases of loading may exist (see Table 2 and Figure 4). When $\dot{\eta} = 0$, no damage occurs regardless the values of $f$ and $\dot{f}$. Damage only occurs when $\dot{\eta} \neq 0$, the value of which is obtained by solving $\dot{f}(Y; d) = 0$ (imposed by equation (2.24)).

Note that from the inequality of Clausius-Duhem $D \ge 0$, and given that $Y \ge 0$, the damage variable $d$ can only grow in time

$$ \dot{d} \ge 0. \qquad (2.25) $$

Thus, during damage ($f = 0$)

$$ \frac{\partial f}{\partial Y} \ge 0, \qquad (2.26) $$

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) for quasi-brittle damage:

$$ f(Y; d) = Y - \kappa(d) \le 0. \qquad (2.27) $$

We presently define $\kappa$ as a function of two parameters, the damage threshold $Y_D \ge 0$ and the hardening modulus $Y_C \ge 0$, such that:

$$ \kappa(d) = Y_D + Y_C d. \qquad (2.28) $$

The size of the domain of admissible states $f \le 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 eq. (2.22) can be written equivalently in an explicit form:

$$ d = < \zeta(Y^{max}) >^+, \qquad (2.29) $$

where $< \cdot >^+$ designates the Macaulay brackets defined by

$$ \begin{cases} <x>^+ = x & \text{if } x \ge 0 \\ <x>^+ = 0 & \text{otherwise.} \end{cases} \qquad (2.30) $$

The function $\zeta(Y) = (Y - Y_D)/Y_C$ designates the reciprocal of the bijection $\kappa$, and $Y^{max}$ is defined by

$$ Y^{max}(t) = \max_{\tau \le t} \{Y(\tau)\}. \qquad (2.31) $$

### 3. 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 (eq. (2.2), (2.6), (2.7), (2.8), (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 forms of the equations, in which the lengths are non-dimensionalized by $a$, time by $1/\dot{\gamma}$ 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
<figure><img src="image_6.png" /><figcaption>Figure 5: Numerical method to solve the fluid-structure interaction problem over a time step.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Projection on the shear plane of the mesh of a damaged capsule captured (a) in the initial configuration and (b) at a steady deformed state. P2 elements, $N_E=1280$, $N_N=2562$.</figcaption></figure>

is computed explicitly at each node from eq. (2.2). Finally, eq. (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 conform 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 $\underline{q} \in L_h^2$ and the damage $d \in L_h^2$ that satisfy eq. (2.7) and eq. (2.29), where the subscript $h$ indicates the finite element space. The position $\underline{x}$ and the virtual displacement $\hat{\underline{x}}$ are searched in $H_h^1$. Using isoparametric elements, we restrict the solution for $\underline{q}$ in $H_h^1$. A field $\underline{v}(\underline{x}, t) \in H_h^1$ writes: $\underline{v}(\underline{x}, t) = N^{(p)}(\underline{x})\underline{v}^{(p)}(t)$, $p \in [1, N_N]$, where $N^{(p)}$ and $\underline{v}^{(p)}$ are the shape function and the nodal coordinates of $\underline{v}$ associated to the node $p$, respectively. Noting $v_{Xj}^{(p)}$ the coordinates of $\underline{v}^{(p)}$ in a Cartesian basis $(\underline{\epsilon}_{Xj})$, the left-hand side of the discretized form of eq. (2.7) writes:

$$ \sum_{el} \hat{u}_{Xj}^{(p)} \int_{el} N^{(p)} N^{(q)} dS \, q_{Xi}^{(q)} = \{\hat{u}_{Xj}^{(p)}\}^{T}[M]\{q_{Xj}^{(p)}\}, \quad (3.1) $$

and the right-hand side writes:

$$ \sum_{el} \int_{el} T^{\alpha\beta}(\underline{\epsilon}, d) \chi_{\alpha\beta}^{(p)Xj} dS \, \hat{u}_{Xj}^{(p)} = \{\hat{u}_{Xj}^{(p)}\}^{T}\{R\}(\underline{\epsilon}, d), \quad (3.2) $$

where $\{q_{Xj}^{(p)}\}$ and $\{\hat{u}_{Xj}^{(p)}\}$ are the vectors of size $3N_N$ of the nodal coordinates, and $\chi_{\alpha\beta}^{(p)Xj}$ is defined by:

$$ \chi_{\alpha\beta}^{(p)Xj} = \frac{1}{2} \left( \frac{\partial N^{(p)}}{\partial \xi_{\beta}} \underline{\epsilon}_{\alpha} + \frac{\partial N^{(p)}}{\partial \xi_{\alpha}} \underline{\epsilon}_{\beta} \right) \cdot \underline{\epsilon}_{Xj}, \quad p \in [1, N_N]. \quad (3.3) $$

Eq. (2.7) being satisfied for any virtual displacement, the discrete solid problem writes:

$$ \begin{cases} [M]\{q\} = \{R\}(\underline{\epsilon}, d) & (3.4a) \\ d = <\kappa^{-1}(Y_{max})>^+ & (3.4b) \end{cases} $$

The square and column matrices $[M]$ and $\{R\}$ are, respectively, computed at each time step by using 6 Hammer points on each element (Hammer et al. 1956). The new value of the damage variable is obtained from eq. (3.4b), solved locally at each integration point while computing $\{R\}$. Knowing the deformation, the variable $d$ is computed explicitly as $Y^{max}$ depends only on the deformation. The computation of $d$ ensures the admissibility condition (eq. 2.27) at each time step. Finally, $\underline{q}$ is computed by solving eq. (3.4a) with the Pardiso solver (Schenk & Gärtner 2004).

## Fluid solver

For a given deformed configuration of the capsule and knowing the stress exerted by the membrane on the fluid, the velocity field $\underline{v}$ is given explicitly by eq. (2.2). The velocity field $\underline{v}$ is computed at each node. The integral on the right-hand side of eq. (2.2) is computed with 12 Hammer points per element. To handle the singularity of the tensor $J$ at node $\underline{x}$, we use polar coordinates centred on $\underline{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 conform mesh with isoparametric elements, the loads $\llbracket \underline{\sigma} \rrbracket \cdot \underline{n}$ and $\underline{q}$ are in the same space $H_h^1$. Hence the dynamic coupling between the fluid and the solid (eq. 2.6) is verified in its strong form in this space. Considering the kinematic coupling, the no-slip
condition (eq. 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).

## 4. 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 + \alpha(t))$, where the inflation ratio $\alpha$ is such that $\alpha \ge 0$. We will study two cases: a monotonic increase of $\alpha$ and cyclic variations of $\alpha$ 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 (eq. 2.29) and the equilibrium of the membrane (eq. 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^{\underline{T}} \cdot F^{\underline{\underline{L}}}$, is simply $\lambda = 1 + \alpha$. The corresponding isotropic principal value $T$ of the tension is:

$$T = (1 - d)G_s(1 - \frac{1}{\lambda^6}), \qquad (4.1)$$

and the elastic energy release rate $Y$:

$$Y = \frac{G_s}{2}(2\lambda^2 + \frac{1}{\lambda^4} - 3). \qquad (4.2)$$

As $Y$ increases monotonically with $\alpha$, the evolution law for damage (2.29) writes:

$$d = < \kappa^{-1}(Y(\alpha^{\text{max}})) >^{+}, \qquad (4.3)$$

where $\alpha^{\text{max}}$ is defined similarly to $Y^{\text{max}}$ in eq. (2.31). Hence, the condition for $d$ to increase is that $\alpha$ is larger than any of its previous values.

The external load is $q = p\underline{n}$, where $p \ge 0$ is the difference between the internal and external pressures. Choosing test functions of the form $\underline{u} = \hat{u}\underline{x}$ in the equilibrium equation (2.7), we obtain the Laplace relation between $T$ and $p$:

$$T = \frac{a(1 + \alpha)p}{2}. \qquad (4.4)$$

We prescribe the radial displacements to the nodes and impose $\underline{x}^{(m)} = (1+\alpha)\underline{X}^{(m)}$, $\forall m \in [0, N_N]$. The pressure difference and damage variable $d$ are obtained analytically using eq. (4.1)-(4.4), and numerically using the solid solver presented in section 3. For the numerical solution, we compute $p$ and $d$ as surface averages, the pressure difference $p$ being given by $q \cdot \underline{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 $\alpha$ in the case of a monotonic inflation of the capsule (Figure 7). The numerical and analytical curves are perfectly superimposed (Figures 7 b-c) and comparison with the analytical solution of the undamaged capsule ($d = 0$) shows a clear effect of damage on the pressure difference (Figure 7 b). Damage is initiated at the critical value $\alpha = \alpha_c$, which corresponds to $Y(\alpha_c) = Y_D$. The loss in elastic properties of the damaged
<figure><img src="image_8.png" /><figcaption>Figure 7: Case of a monotonic inflation: for the stretch ratio $\alpha$ shown in (a), corresponding curves of the dimensionless pressure difference (b) and of the damage variable (c), computed for $Y_D = 0.2$, $Y_C = 2.0$.</figcaption></figure>

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 $\alpha = \alpha_\ell$.

We then compare the evolution of $ap/G_s$ and $d$ with $\alpha$ 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 $\alpha$ does not exceed the critical value $\alpha_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 (Figures 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 $\alpha$ 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 $\alpha$ reaches the limit value $\alpha_\ell$ (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 $\alpha$, 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(\alpha)$ decreases (Figure 8b), which means that damage reduces the stiffness of the structure. It is interesting to see how the values of $\alpha_c$ and $\alpha_\ell$ depend on the
<figure><img src="image_9.png" /><figcaption>Figure 8: Case of cycles of inflations and deflations with increasing maximum capsule size: for the stretch ratio $\alpha$ shown in (a), corresponding curves of the dimensionless pressure difference (b) and of the damage variable (c), computed for $Y_D = 0.2$, $Y_C = 2.0$.</figcaption></figure>

parameter values $Y_D$ and $Y_C$. Following equation (4.3), the values of $\alpha$ initiating damage and rupture are given respectively by the equations $Y(\alpha_c) = Y_D$ and $Y(\alpha_\ell) = Y_D + Y_C$, where the expression of $Y(\alpha)$ is obtained using equation (4.2). The critical inflation ratio $\alpha_c$ thus depends solely on $Y_D$, but the limit inflation ratio $\alpha_\ell$ depends on both $Y_D$ and $Y_C$. Furthermore, the higher the parameter values, the higher the two threshold inflation ratios.

## 5. 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_\ell$ is reached, rupture initiates putting a limit to the damage regime. In the last part of this section, we will finally study the influence of the material parameters $Y_D$ and $Y_C$ on the three identified regimes and on the values of $Ca_c$ and $Ca_\ell$.

### 5.1. *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_\ell$. Hence the capsule is damaged but the damage stabilizes and a steady state is reached.
<figure><img src="image_10.png" /><figcaption>Figure 9: Two principal ellipses of the ellipsoid of inertia of the capsule.</figcaption></figure>

$t = t_c$

$t_c < t < t_1$

$t = t_1$

$t = t^\infty$

<figure><img src="image_11.png" /><figcaption>Figure 10: Map of damage at the instant of initiation of damage $t_c$, at an intermediary instant between $t_c$ and the instant of maximum elongation $t_1$, at time $t_1$, and at steady state ($t^\infty$). The map of damage is represented on the current and reference configurations. The current configuration is observed in the shear plane ($O, e_x, e_z$) and in the plane ($O, e_1, e_y$) which is defined in Figure 9. The reference configuration of the capsule is observed in the shear plane ($O, e_x, e_z$). The points $P$ and $P'$ correspond to the intersection of the membrane with the vorticity axis $e_y$. The results are obtained for $Ca = 0.7$, $Y_D = 0.2$ and $Y_C = 2.0$. All the pictures are at the same scale. The capsule is delimited by a black line.</figcaption></figure>

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
<figure><img src="image_12.png" /><figcaption>Figure 11: Time evolution of different state quantities: elastic energy release rate $Y$ (first row), maximum principal tension $T_1$ (second row), normal load $\underline{q} \cdot \underline{n}$ to visualize wrinkling (third row) and negative part of principal tension $T_2$ (fourth row). The results are shown in the shear plane $(O, \underline{e}_x, \underline{e}_z)$ at the same instants as in Figure 10, for $Ca = 0.7$, $Y_D = 0.2$ and $Y_C = 2.0$.</figcaption></figure>

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_\infty$). 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 $\underline{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}^\infty$ are found at $P$ and $P'$ at steady state ($t = t_\infty$). 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 $\underline{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 localisation 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$ (first and second rows of Figure 11). Damage has no visible consequences on membrane wrinkling: the wrinkles visible on the normal load maps in the third row of Figure 11 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 (fourth row of Figure 11). 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
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure 12: Temporal evolution of (a) the lengths of the axes of the ellipsoid of inertia, (b) the Taylor parameter $D_{12}$, (c) the inclination angle $\beta$, and (d) the global surface expansion $\lambda_s$. Computed for $Ca = 0.7$, $Y_D = 0.2$ and $Y_C = 2.0$.</figcaption></figure>

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 is 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 $\beta$ between the flow direction $e_x$ and the direction of the major axis $e_1$. Figure 12c represents the temporal evolution of $\beta$ showing that the inclination angle decreases from the first measurable value near $\pi/4$. Figure 12d shows the evolution of the global surface expansion ratio $\lambda_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 $\infty$ hereafter. It is interesting to notice in Figure 12a 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
<figure><img src="image_17.png" /></figure>

Figure 13: Maximum damage value at steady state $d_{max}^{\infty}$ with respect to Ca for $Y_D = 0.2$ and $Y_C = 2$. The inserted images represent the map of damage in the shear plane at steady state for $Ca = 0.6$ (a), for $Ca = 0.7$ (b) and at the instant of initiation of rupture $t = t_{\ell}$ for $Ca = 0.8$ (c). The colormap for $d$ is saturated for values of $d$ larger than 0.2.

damage. But only small differences are observed on the principal lengths $L_i$ (Figure 12a), $D_{12}$ (Figure 12b) and $\beta$ (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 $\lambda_S$ (Figure 12d). Nevertheless, the difference at steady state is only of a few percent.

## 5.2. 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_\ell = 0.73$. The maximum value of damage at steady state $d_{max}^\infty$ 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 \sim 0.6$. Above, $d_{max}^\infty$ increases more rapidly with $Ca$ until $d_{max}^\infty = 0.4$. Around $Ca = Ca_\ell$, 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_\ell$.

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_\ell$ 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 \le Ca_\ell$. No results are shown above $Ca_\ell$, since no steady deformed shape exists any longer ($D_{12}$ diverges to infinity). Despite the large effect of
<figure><img src="image_18.png" /><figcaption>Figure 14: Evolution of the values (a) $D_{12}^\infty$ and (b) $\beta^\infty$, respectively the values of $D_{12}$ and $\beta$ at steady state, in relation to Ca. Computed for $Y_D = 0.2$ and $Y_C = 2.0$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 15: Influence of the parameters $Y_C$ and $Y_D$ on the evolution of the damage value $d_{max}^\infty$ at points P and $P'$ at steady state with Ca: (a) $Y_C = 2.0$ and $Y_D = 0.1, 0.2, 0.3$, (b) $Y_C = 1.0, 2.0, 3.0$ and $Y_D = 0.2$.</figcaption></figure>

Ca on $d_{max}^\infty$ for $Ca_c < Ca < Ca_\ell$, the $D_{12}^\infty$ and $\beta^\infty$ curves initially remain superimposed to the non-damaged case, and it is only close to $Ca_\ell$ that small differences become visible. No significant influence of damage is, thus, found on these global quantities. It is a consequence of the localisation 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_\ell$, which corresponding to when $D_{12}$ diverges to infinity (initiation of breakup).

### 5.3. 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}^\infty(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_\ell$ 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_\ell$). This relates to the facts that the criterion of initiation of damage is only a function of $Y_D$, whereas the
<figure><img src="image_20.png" /><figcaption>Figure 16: Evolution of the critical and limit capillary numbers $Ca_c$ and $Ca_\ell$ with $Y_D$. The solid lines represent the limit curves of $Ca_c$ and $Ca_\ell$ for $Y_C = 0.2$. They delimit three domains corresponding to three states of the capsule: undamaged, damaged and ruptured. We also show the limit curves of $Ca_\ell$ for $Y_C = 1.0, 3.0$ as dotted lines to show how the three domains evolve with the parameters.</figcaption></figure>

criterion of initiation of rupture is controlled by $Y_D + Y_C$, as already shown at the end of section 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_\ell(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_\ell$, ruptured for $Ca > Ca_\ell$. The only effect of $Y_C$ is to shift the $Ca_\ell$ 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$).

# 6. Discussion and conclusion

In response to the current need for a damage model of microcapsules in flow, we have developed the first fluid-structure interaction (FSI) numerical model accounting for the degradation of the capsule membrane till 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.
<figure><img src="image_21.png" /><figcaption>Figure 17: The RVE consists of a bundle of elastic initially unbroken fibres of specific elastic energy $\phi_{NH}$ and probability of rupture $P_f$ (eq. (6.3)). It is subjected to an elongation $\lambda$ up to the maximum elongation ratio $\lambda^{max}$. The zone where the microdefects have appeared upon the rupture of the fibres are represented in grey.</figcaption></figure>

### 6.1. Interpretation of the damage evolution law

We have explained the physical meaning of the damage model in section 2, but propose to further detail the interpretation of the damage evolution law eq. (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 $\lambda$.

Each fibre is associated to the specific elastic energy $\phi_{NH}$ and has a brittle behaviour given by the classical energetic criterion of rupture

$$ \begin{cases} \phi_{NH}(\lambda^{max}) < \phi_u \Rightarrow \text{ sound fibre} \\ \phi_{NH}(\lambda^{max}) = \phi_u \Rightarrow \text{ broken fibre,} \end{cases} \qquad (6.1) $$

where $\phi_u$ is a specific energy at rupture and $\lambda^{max}$ is defined similarly to $Y^{max}$ in eq. (2.31). The key ingredient of this model is to consider $\phi_u$ as a random variable with probability density $p(\phi_u)$ given by the following band-limited and uniform probability density

$$ p(\phi_u) = \begin{cases} \frac{1}{Y_C} & \forall \phi_u \in [Y_D, Y_D + Y_C] \\ 0 & \forall \phi_u \notin [Y_D, Y_D + Y_C]. \end{cases} \qquad (6.2) $$

where $Y_D$ and $Y_C$ are the parameters of the damage model introduced in eq. (2.28). Hence, from eqs. (6.1) and (6.2), the probability of rupture of a fibre is given by

$$ P_f(\lambda^{max}) = \int_0^{\phi_{NH}(\lambda^{max})} p(\phi_u)d\phi_u \qquad (6.3) $$

Consistent with eq. (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(\lambda^{max})N$, and thus $d = P_f(\lambda^{max})$. From eq. (6.3), we obtain

$$d = \begin{cases} 0 & \text{if } \phi_{NH}(\lambda^{max}) \le Y_D \\ \frac{\phi_{NH}(\lambda^{max}) - Y_D}{Y_C} & \text{if } Y_D \le \phi_{NH}(\lambda^{max}) \le Y_D + Y_C \\ 1 & \text{if } Y_D + Y_C \le \phi_{NH}(\lambda^{max}) \end{cases} \quad (6.4)$$

where $\phi_{NH}(\lambda^{max}) = Y^{max}$ (see eq. (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 (eq. (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). Eq. (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.

## 6.2. Capsule inflation test

We have first applied the model to a capsule under isotropic inflation, 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 $\alpha_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: it 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 $\alpha_\ell$ 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 et al. 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 et al. 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.
### 6.3. Capsule under simple shear flow

We have then considered a capsule under simple shear flow and similarly seen that there exist 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.

### 6.4. 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 et al. 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 the group of Prof. H. Rehage 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 section 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 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
<figure><img src="image_22.png" /><figcaption>Figure 18: Experimental results obtained by Koleva & Rehage (2012) on polysiloxane microcapsules: (a)Evolution of the Taylor parameter $D_{12}$ with the capillary number $Ca$, (b) Formation of wrinkles at $Ca = 0.0042$, (c) Formation of a crack at $Ca = 0.01$, (d) Divergence of the capsule shape for $Ca$ larger than $Ca = 0.01$. Numerical predictions given by the present damage FSI model for $Y_D/G_s = 5 \times 10^{-4}$ and $Y_C/G_s = 3.1 \times 10^{-3}$: (e) At $Ca = 0.005$, 3D rendering showing the presence of wrinkles. (f) At $Ca = 0.01$, map of damage at $t = t_\ell$ when rupture initiates ($d = 1$). (g) At $Ca = 0.012$, map of damage at an instant of time after $t = t_\ell$ while the capsule shape diverges due to infinite elongation, this is a case where damage initiates at the points on the vorticity axis but not rupture, which occurs in the nearby region. Pictures (a-d) are reproduced from (Koleva & Rehage 2012) with permission of The Royal Society of Chemistry.</figcaption></figure>

in Koleva & Rehage (2012) the capillary number is based on the surface Young 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_\ell(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 \in [0; 1 \times 10^{-3}]$ and $Y_C \in [0; 4.6 \times 10^{-3}]$. We have run simulations assuming $Y_D/G_s = 5 \times 10^{-4}$ and $Y_C/G_s = 3.1 \times 10^{-3}$, for which $Ca_\ell(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 \ge Ca_\ell$, 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 the group of Prof. Rehage. 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 (eq. (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.

## Acknowledgements

This work was supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (MultiphysMicroCaps, No. ERC-2017-COG 772191) and by the French Ministry of Higher Education and Research (fellowship of Nicolas Grandmaison).

## Declaration of Interests

The authors report no conflict of interest.

## REFERENCES

* ABKARIAN, M., FAIVRE, M., HORTON, R., SMISTRUP, K., BEST-POPESCU, C. A. & STONE, H. A. 2008 Cellular-scale hydrodynamics. *Biomed. Mater.* **3**, 034011.

* BARTHÈS-BIESEL, D., DIAZ, A. & DHENIN, E. 2002 Effect of constitutive laws for two-dimensional membranes on flow-induced capsule deformation. *J. Fluid Mech.* **460**, 211–222.

* BAŽANT, Z. P. & PIJAUDIER-CABOT, G. 1988 Nonlocal continuum damage, localization instability and convergence. *J. Appl. Mech.* **55**, 287–293.

* BENALLAL, A., BILLARDON, R. & GEYMONAT, G. 1993 Bifurcation and localization in rate-independent materials. some general considerations. In *Bifurcation and Stability of Dissipative Systems* (ed. Q. S. Nguyen), pp. 1–44. Springer.

* BESSON, J., CAILLETAUD, G., CHABOCHE, J.-L. & FOREST, S. 2010 *Non-Linear Mechanics of Materials*, pp. 127–194. Springer.

* BHUJBAL, S. V., DE VOSL, P. & NICLOU, S. P. 2014 Drug and cell encapsulation: Alternative delivery options for the treatment of malignant brain tumors. *Adv. Drug Deliv. Rev.* **67-68**, 142–153.
CHANG, K. S. & OLBRICHT, W. L. 1993 Experimental studies of the deformation and breakup of a synthetic capsule in steady and unsteady simple shear flow. *J. Fluid Mech.* **250**, 609–633.

CHU, T. X., SALSAC, A.-V., LECLERC, E., BARTHÈS-BIESEL, D., WURTZ, H. & EDWARDS-LÉVY, F. 2011 Comparison between measurements of elasticity and free amino group content of ovalbumin microcapsule membranes: Discrimination of the cross-linking degree. *J. Colloid Interf. Sci.* **355**, 81–88.

COLEMAN, B. D. & GURTIN, M. E. 1967 Thermodynamics with internal state variables. *J. Chem. Phys.* **47**, 597–613.

GHAEMI, A., PHILIPP, A., BAUER, A., LAST, K., FERY, A. & GEKLE, S. 2016 Mechanical behaviour of micro-capsules and their rupture under compression. *Chem. Eng. Sci.* **142**, 236–243.

GUBSPUN, J., GIRES, P.-Y., DE LOUBENS, C., BARTHÈS-BIESEL, D., DESCHAMPS, J., GEORGELIN, M., LEONETTI, M., LECLERC, E., EDWARDS-LÉVY, F. & SALSAC, A.-V. 2016 Characterization of the mechanical properties of cross-linked serum albumin microcapsules: effect of size and protein concentration. *Colloid Polym. Sci.* **294**, 1381–1389.

HAMMER, P. C., MARLOWE, O. J. & STROUD, A. H. 1956 Numerical Integration Over Simplexes and Cones. *Mathematical Tables and Other Aids to Computation* **10**, 130–137.

HOKANSON, J. & YAZDANI, S. 1997 A constitutive model of the artery with damage. *Mechanics Research Communications* **24**, 151–159.

HOLZAPFEL, G. A. & FEREIDOONNEZHAD, B. 2017 Modeling of damage in soft biological tissues. In *Biomechanics of Living Organs* (ed. Y. Payan & J. Ohayon), pp. 101–123. Elsevier.

HUSMANN, M., REHAGE, H., DHENIN, E. & BARTHÈS-BIESEL, D. 2005 Deformation and bursting of nonspherical polysiloxane microcapsules in a spinning-drop apparatus. *J. Colloid Interf. Sci.* **282**, 109–119.

KACHANOV, L. M. 1958 On time to rupture in creep conditions (in russian). *Isv Akad. Nauk. SSR Otd. Tekh. Nauk.* **8**, 26–31.

KACHANOV, L. M. 1986 *Introduction to continuum damage mechanics*. Dordrecht: Springer.

KOLEVA, I. & REHAGE, H. 2012 Deformation and orientation dynamics of polysiloxane microcapsules in linear shear flow. *Soft Matter* **8**, 3681–3693.

KRAJCINOVIC, D. 1989 Damage mechanics. *Mechanics of Materials* **8**, 117–197.

LAC, É., BARTHÈS-BIESEL, D., PELEKASIS, N. A. & TSAMOPOULOS, J. 2004 Spherical capsules in three-dimensional unbounded Stokes flows: effect of the membrane constitutive law and onset of buckling. *J. Fluid Mech.* **516**, 303–334.

LE GOFF, A., KAOUI, B., KURZAWA, G., HASZON, B. & SALSAC, A.-V. 2017 Squeezing bio-capsules into a constriction: deformation till break-up. *Soft Matter* **13**, 7644–7648.

LEMAITRE, J. & DESMORAT, R. 2005 *Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures*. Berlin Heidelberg: Springer-Verlag.

LI, S., NICKELS, J. & PALMER, A. F. 2005 Liposome-encapsulated actin-hemoglobin (LEAcHb) artificial blood substitutes. *Biomaterials* **26**, 3759–3769.

MARIGO, J.-J. 1981 Formulation d'une loi d'endommagement d'un matériau élastique. *C. R. Acad. Sci. Paris. Série II* **292**, 1309–1312.

MOËS, N. & BELYTSCHKO, T. 2002 Extended finite element method for cohesive crack growth. *Engng. Fract. Mech.* **69**, 813–833.

NATALI, A.N., PAVAN, P.G., CARNIEL, E.L., LUCISANO, M.E. & TAGLIALAVORO, G. 2005 Anisotropic elasto-damage constitutive model for the biomechanical analysis of tendons. *Medical Engineering & Physics* **27**, 209–214.

POZRIKIDIS, C. 1992 *Boundary Integral and Singularity Methods for Linearized Viscous Flow*. Cambridge: Cambridge University Press.

RICE, J. R. 1976 The localization of plastic deformation. In *Theoretical and Applied Mechanics* (ed. W. T. Koiter), pp. 207–220. Proceedings of the 14th ICTAM, Delft (Netherlands): North-Holland Publishing Company.

SCHENK, O. & GÄRTNER, K. 2004 Solving unsymmetric sparse systems of linear equations with PARDISO. *Future Generation Computer Systems* **20**, 475–487.

SIMO, J. C., OLIVER, J. & ARMERO, F. 1993 An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. *Comput. Mech.* **12**, 277–296.
SU, J., HU, B.-H., LOWE, W. L., KAUFMAN, D. B. & MESSERSMITH, P. B. 2010 Anti-inflammatory peptide-functionalized hydrogels for insulin-secreting cell encapsulation. *Biomaterials* **31**, 308–314.

SUKUMAR, N., MOËS, N., MORAN, B. & BELYTSCHKO, T. 2000 Extended finite element method for three-dimensional crack modelling. *Int. J. Numer. Meth. Eng.* **48**, 1549–1570.

UNVERFEHRT, A., KOLEVA, I. & REHAGE, H. 2015 Deformation, orientation and bursting of microcapsules in simple shear flow: Wrinkling processes, tumbling and swinging motions. *Journal of Physics Conference Series* **602**, 012002.

WALTER, A., REHAGE, H. & LEONHARD, H. 2001 Shear induced deformation of microcapsules: shape oscillations and membrane folding. *Colloids Surf.*, A **183-185**, 123–132.

WALTER, J., SALSAC, A.-V., BARTHÈS-BIESEL, D. & LE TALLEC, P. 2010 Coupling of finite element and boundary integral methods for a capsule in a Stokes flow. *Int. J. Numer. Meth. Eng.* **83**, 829–850.