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

J. Walter, A.-V. Salsac, D. Barthès-Biesel, Patrick Le Tallec

► To cite this version:

J. Walter, A.-V. Salsac, D. Barthès-Biesel, Patrick Le Tallec. Coupling of finite element and boundary integral methods for a capsule in a Stokes flow. International Journal for Numerical Methods in Engineering, 2010, 83 (7), pp.829-850. 10.1002/nme.2859 . hal-00542306

HAL Id: hal-00542306

https://hal.science/hal-00542306v1

Submitted on 24 Mar 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 CC BY 4.0 - Attribution - International License
Coupling of finite element and boundary integral methods
for a capsule in a Stokes flow

J. Walter¹, A.-V. Salsac¹, D. Barthès-Biesel¹,*,† and P. Le Tallec²

¹UMR CNRS 6600, Biomécanique et Bioingénierie, Université de Technologie de Compiègne, BP 20529, 60205 Compiègne, France

²UMR CNRS 7649, Laboratoire de Mécanique des Solides, École Polytechnique, 91128 Palaiseau Cedex, France

SUMMARY

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

KEY WORDS: fluid-structure interaction; finite element method; boundary integral method; membrane model

# 1. INTRODUCTION

Synthetic capsules are liquid drops or gels protected by a thin elastic membrane and have a very wide range of applications in cosmetic, food, and pharmaceutical industries. They are used to protect fragile, volatile, or active substances and to control the release rate as desired. Other instances are found in biomedical engineering, where microcapsules are used for cell encapsulation or artificial blood manufacturing [1, 2]. In many applications, what is at stake is the characterization

<footnote>*Correspondence to: D. Barthès-Biesel, Université de Technologie de Compiègne, BP 20529, 60205 Compiègne, France.</footnote>

<footnote>†E-mail: dbb@utc.fr</footnote>

<footnote>Contract/grant sponsor: Conseil Régional de Picardie through project µFlec</footnote>

<footnote>Contract/grant sponsor: French Ministère de la Recherche through project Pilcam2</footnote>
of the complex behavior of a deformable capsule in an external liquid flow, to control or prevent the breakup of the capsule and the liberation of its contents. Numerical models of the fluid-structure interactions are thus necessary to predict the capsule deformation in flow conditions and compute the stress level in the membrane.

In many respects, the fluid-structure interaction of microcapsules suspended in an external flow is unconventional. (1) The capsules are closed surfaces freely suspended in a fluid. They are, therefore, not subjected to any boundary condition on the displacement apart from spatial periodicity. (2) At the small scale of the capsule, the inertia forces of the internal and external flows are negligible compared with the viscous forces. The fluid flows can, therefore, be modeled using the Stokes equations (balance between viscous and pressure forces). The load on the membrane is then due to pressure forces and viscous tractions with normal and in-plane components. (3) The inertia force of the capsule wall can be equally neglected. Consequently, the problem does not contain any dynamics and can be seen as a succession of equilibrium states. The drawback is that the stability of the equilibrium is not necessarily satisfied. (4) The capsule wall is subjected to not only large displacements but also large deformations. For instance, average stretch ratios of up to 1.2 are reported in experimental studies of capsules flowing in a pore [3, 4]. The non-linear effects due to these large deformations must be taken into account in the membrane mechanics model.

An analytical solution can only be found in the approximation of small deformations, using perturbation methods [5]. In the case of large deformations, a numerical simulation is needed to solve the fully coupled equations governing the fluid and membrane mechanics. This problem has been studied for the past three decades and different techniques have been considered to derive numerical solutions (as summarized in Table I). Many studies have used a coupling strategy based on the boundary integral method to solve the Stokes flow and the membrane elasticity equations on the same mesh [6–9]. The velocity field at any position within the fluid domain is given by surface integrals calculated on the geometric boundaries. This method, therefore, has the advantage of reducing the geometric dimension of the problem by one, which largely decreases the total number of nodes. It also does not involve discretizing spatial derivatives (using for example finite differences) and is very accurate. In particular, it has been shown to be efficient, precise, and stable when modeling the deformation of capsules subjected to shear flows. Another coupling method used for capsule simulation is the immersed boundary method [10–12]. Two grids are used in the method: a stationary 3D grid for the fluid flow and a moving 2D boundary grid for the interface. The forces exerted by the membrane on the fluid and the flow velocity convecting the membrane are applied locally from one grid to the other using approximate Dirac functions [13, 14]. The Navier–Stokes equations are generally solved using finite difference schemes. The immersed boundary method has the advantage over the boundary integral method to be applicable to non-zero Reynolds numbers and non-Newtonian fluids. However, the method does not treat the membrane as a physical boundary of the fluid domain, which can result in a lack of precision in the Lagrangian tracking of the interface.

The most frequently used model for the capsule wall is that of a 2D hyperelastic surface: the wall is considered to be infinitely thin and to have a negligible bending stiffness. Two approaches may be considered to model the capsule membrane mechanics: the equations of the force equilibrium on the capsule wall may either be written locally at each point (strong form) or multiplied by a test function and integrated over the capsule surface to define a variational problem (weak form). Most capsule studies have used the strong form of the equations. Numerical simulations of such models were first performed in axisymmetric configurations with a collocation technique to determine the unknown forces and velocities at discrete locations in elongational flows [6, 15] and in pore
<table><caption>Table I. Summary of the previous work, classified according to the solid (rows) and fluid (columns) methods used to model the fluid-structure interaction of a capsule in a flow.</caption><thead><tr><th></th><th>Boundary element method</th><th>Immersed boundary method</th></tr></thead><tbody><tr><td>Local equilibrium</td><td></td><td></td></tr><tr><td>3D axisymmetric</td><td>[6, 15–17]</td><td></td></tr><tr><td>3D with piecewise constant loads</td><td>[7, 18]</td><td>[12]</td></tr><tr><td>3D with B-spline discretization</td><td>[8, 19]</td><td></td></tr><tr><td>3D with spectral discretization</td><td>[9]</td><td></td></tr><tr><td>Variational method</td><td></td><td></td></tr><tr><td>3D with finite elements</td><td></td><td>[10, 11, 20]</td></tr></tbody></table>

The capsule wall was each time modeled as a 2D membrane.

flows [16, 17]. The full 3D problem was then solved for infinite flow configurations. Capsules in simple shear flows have been considered by Pozrikidis [18], Ramanujan and Pozrikidis [7] and more recently by Li and Sarkar [12], who computed the membrane load as a piecewise constant function. Lac et al. [8] and Lac and Barthès-Biesel [19] used bi-cubic B-splines instead as interpolation functions to compute the loads with high accuracy, whereas Dodson and Dimitrakopoulos [9] opted for a spectral discretization of the problem. An alternative option is to write the equilibrium equations in their weak form and to use a finite element method. The local equilibrium equations are converted into a variational problem. Only two groups [10, 11, 20] have implemented a finite element method that was each time coupled with an immersed boundary method to compute the motion of a capsule in shear flows. However, their method [21] is based on the use of linear interpolation functions and lacks generality in the implementation.

Previous studies have shown that, for certain ranges of the governing parameters, the capsule wall undergoes in-plane compression and tends to buckle [8, 12]. Thus, bending effects must be taken into account to properly model the wall mechanics. However, no satisfying model for the bending stiffness of a capsule has been implemented yet. The finite element method has the advantage of offering a framework suitable for the modeling of thin shells. In any case, finite element models seem to be the right option if one wishes to have a general and versatile model involving realistic constitutive laws and post-buckling effects.

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

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

## 2. PROBLEM STATEMENT

We consider an initially spherical capsule (radius $a$) consisting of a liquid droplet enclosed by an infinitely thin membrane characterized by a surface shear modulus $G_s$. The capsule is suspended in a viscous liquid undergoing an unbounded shear flow. The frame of reference with Cartesian basis ($e_1, e_2, e_3$) is centered on the center of mass of the capsule (Figure 1). The velocity of the unperturbed flow is denoted $v^\infty$ with characteristic shear rate $\dot{\gamma}$. Both fluids are assumed to be incompressible, Newtonian and to have the same viscosity $\mu$ and density $\rho$. Dimensional analysis shows that an important non-dimensional parameter is the capillary number

$$Ca = \frac{\mu \dot{\gamma} a}{G_s}, \qquad (1)$$

which compares the viscous stresses exerted by the fluids to the elastic stiffness of the membrane.

### 2.1. Internal and external flows

The particle Reynolds number $Re = \rho a^2 \dot{\gamma} / \mu$ is assumed to be small compared to unity, so that both the internal and external flows follow the Stokes equations:

$$\nabla \cdot \mathbf{v} = 0, \quad \nabla \cdot \boldsymbol{\sigma} = 0, \qquad (2)$$

where $\mathbf{v}$ is the Eulerian velocity field and $\boldsymbol{\sigma}$ is the Cauchy stress tensor in the fluids.

The velocity of the points of the capsule can then be related to the tractions on the membrane through an integral equation over the deformed surface $S$ [22, 23]:

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

<figure><img src="image_2.png" /><figcaption>Figure 1. Capsule freely suspended in a Stokes flow.</figcaption></figure>
where $[\sigma] \cdot \mathbf{n} = (\sigma^{\text{ext}} - \sigma^{\text{int}}) \cdot \mathbf{n}$ is the traction jump across the interface ($\mathbf{n}$ is the unit normal vector pointing outwards and the superscripts indicate the internal and external flows), and

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

is Green's single layer kernel, with $\mathbf{r} = \mathbf{x} - \mathbf{y}$ and $r = \|\mathbf{r}\|$.

## 2.2. Membrane mechanics

Following Skalak et al. [24], we treat the membrane as a purely bidimensional sheet of hyperelastic material without bending moments. This is equivalent to considering a membrane with a finite thickness $\vartheta$ and supposing that the strains are constant across the thickness and that the external forces are exerted by the two flows on the midsurface $S$. Each material point of $S$ is characterized by two curvilinear coordinates $(\xi^1, \xi^2)$. The position of a given material point is thus given as $\mathbf{X}(\xi^1, \xi^2)$ in the reference state and $\mathbf{x}(\xi^1, \xi^2, t)$ in the deformed state. The displacement $\mathbf{U}(\mathbf{X}, t) = \mathbf{x}(\mathbf{X}, t) - \mathbf{X}$ is related to the velocity of the membrane through the kinematic condition:

$$ \frac{\partial}{\partial t} \mathbf{U}(\mathbf{X}, t) = \mathbf{v}(\mathbf{x}, t). \qquad (5) $$

Dimensional analysis shows that the ratio of the inertial effects in the membrane to the viscous stresses exerted by the flow is of order

$$ \frac{\text{membrane inertia}}{\text{viscosity}} \sim \frac{\rho_{\text{mb}}}{\rho} \frac{\vartheta}{a} Re, \qquad (6) $$

where $\rho_{\text{mb}}$ is the density of the membrane material. Thus, with the density ratio of order 1 for an organic membrane suspended in water, a membrane of negligible thickness and a vanishing Reynolds number, the inertia of the membrane can be neglected. The motion of the membrane, therefore, follows the local equilibrium equation:

$$ \nabla_s \cdot \mathbf{T} + \mathbf{q} = \mathbf{0}, \qquad (7) $$

where $\mathbf{T}$ are the bidimensional Cauchy tensions in the membrane and $\nabla_s \cdot$ is the surface divergence operator on the deformed configuration. The load $\mathbf{q}$ on the membrane is equal to the traction jump:

$$ \mathbf{q} = [\sigma] \cdot \mathbf{n}. \qquad (8) $$

Equation (7) can be rewritten as a variational problem, where the load $\mathbf{q}$ is unknown. Let $\mathcal{V}$ be the Sobolev space $H^1(S)$, and $\hat{\mathbf{u}} \in \mathcal{V}$ a virtual displacement field. The variational problem corresponding to Equation (7) is thus:

$$ \int_S \hat{\mathbf{u}} \cdot \mathbf{q} dS = \int_S \hat{\mathbf{e}}(\hat{\mathbf{u}}) : \mathbf{T}(\mathbf{U}) dS, \quad \forall \hat{\mathbf{u}} \in \mathcal{V}, \qquad (9) $$

where $\hat{\mathbf{e}}(\hat{\mathbf{u}}) = \frac{1}{2}(\nabla_s \hat{\mathbf{u}} + \nabla_s \hat{\mathbf{u}}^T)$.

The test function $\hat{\mathbf{u}}$ must be in $H^1(S)$ to allow for a direct calculation and integration of the strain tensor $\hat{\mathbf{e}}(\hat{\mathbf{u}})$. Note that neither the virtual work principle nor the boundary integral (Equation (3)) requires that the load $\mathbf{q}$ be in $H^1(S)$. In the virtual work principle, the load must only be in the dual space $\mathcal{V}'$. However, in the boundary integral equation, it is convenient to choose $\mathbf{q} \in L^2(S)$
allowing for an explicit calculation of the integral. A less regular load would require to integrate Equation (3) by part, which is very delicate as the gradient of the kernel **J** is highly singular for **x**=**y**.

2.3. Stress–strain relations

To close the problem, we need the (non-linear) relation between the displacement **U** and the
tensions **T**. Since we study the large deformation of a curved surface, it is convenient to use the
local covariant (**a**<sub>1</sub>, **a**<sub>2</sub>, **a**<sub>3</sub>) and contravariant (**a**<sup>1</sup>, **a**<sup>2</sup>, **a**<sup>3</sup>) bases constructed from the curvilinear
coordinates ξ<sup>1</sup> and ξ<sup>2</sup> [25, 26]. The vector **a**<sub>3</sub> = **a**<sup>3</sup> = **n** is the unit normal vector (pointing outward)
and the tangent vectors are defined by **a**<sub>α</sub> = **x**<sub>α</sub>, where ,<sub>α</sub> denotes derivation with respect to ξ<sup>α</sup>.
Throughout the paper, Greek indices are equal to 1 or 2, while Latin indices are equal to 1, 2, or
3, and summation is performed over repeated indices. A vector **v** can thus be written as:

$$
\mathbf{v} = v^i \mathbf{a}_i = v_i \mathbf{a}^i = v_{X_i} \mathbf{e}_i, \quad (10)
$$

where the subscript X indicates Cartesian components. The representations of the metric tensor
are defined as a<sub>αβ</sub> = **a**<sub>α</sub> ⋅ **a**<sub>β</sub> and a<sup>αβ</sup> = **a**<sup>α</sup> ⋅ **a**<sup>β</sup>. The same quantities are defined in the reference state,
using capital letters: the local bases are denoted (A<sub>1</sub>, A<sub>2</sub>, A<sub>3</sub>) and (A<sup>1</sup>, A<sup>2</sup>, A<sup>3</sup>), and the covariant
and contravariant representations of the metric tensor are A<sub>αβ</sub> and A<sup>αβ</sup>, respectively.

The gradient of the transformation **F**=**a**<sub>α</sub> ⊗ **A**<sup>α</sup> is introduced such that **dx**=**F**·**dX**. The Green–Lagrange strain tensor is **e**=1/2(**F**<sup>T</sup>·**F**–**1**). Invariants of the transformation can then be defined as:

$$
I_1 = \operatorname{tr}(\mathbf{F}^T \cdot \mathbf{F}) - 2 = A^{\alpha\beta} a_{\alpha\beta} - 2, \quad I_2 = \det(\mathbf{F}^T \cdot \mathbf{F}) - 1 = |A^{\alpha\beta}| |a_{\alpha\beta}| - 1. \quad (11)
$$

Assuming the membrane to be 2D isotropic, the Cauchy tensions are related to a strain energy
function w<sub>s</sub>(I<sub>1</sub>, I<sub>2</sub>) per unit undeformed surface area by:

$$
\mathbf{T} = \frac{1}{J_s} \mathbf{F} \cdot \frac{\partial w_s}{\partial \mathbf{e}} \cdot \mathbf{F}^T, \quad (12)
$$

where the Jacobian *J*<sub>s</sub>=det **F** represents the ratio between the deformed and undeformed surface areas. Using the chain rule, one obtains the following expression for the contravariant representation of **T** [25]:

$$
T^{\alpha\beta} = \frac{2}{J_s} \frac{\partial w_s}{\partial I_1} A^{\alpha\beta} + 2J_s \frac{\partial w_s}{\partial I_2} a^{\alpha\beta}. \quad (13)
$$

Several laws have been proposed to model the behavior of hyperelastic membranes [27].
A classical 3D law used to describe volume-incompressible, rubber-like materials is the
neo-Hookean (NH) law, with a 2D form given by:

$$
w_s^{\text{NH}} = \frac{G_s^{\text{NH}}}{2} \left( I_1 - 1 + \frac{1}{I_2 + 1} \right). \quad (14)
$$

Skalak *et al.* [24] proposed a law (Sk) with independent surface shear and area dilation moduli:

$$
w_s^{Sk} = \frac{G_s^{Sk}}{4} (I_1^2 + 2I_1 - 2I_2 + C I_2^2), \quad C > -1/2, \quad (15)
$$
where the parameter C is such that the ratio of the area dilation modulus to the surface shear modulus is 1+2C. For large values of C, the Sk law describes an almost area-incompressible material, which is appropriate for biological membranes. It has also been shown to model effectively the behavior of alginate artificial capsules when using C=0 [28]. Note that this law was developed directly from the perspective of 2D elasticity and only assumes in-plane isotropy. From a 3D perspective, it could describe a material either fully isotropic or transversely isotropic.

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

$$w_s^H = G_s \left( \text{tr}(\varepsilon^2) + \frac{v_s}{1-v_s} (\text{tr}\varepsilon)^2 \right), \quad (16)$$

where $\varepsilon$ is the 2D linearized Green-Lagrange strain tensor and $v_s \in ]-1, 1[$ is the 2D Poisson ratio. The equivalence between the laws is obtained with

$$G_s = G_s^{NH}, \quad v_s = 1/2, \quad (17)$$

for the NH law and

$$G_s = G_s^{Sk}, \quad v_s = \frac{C}{1+C}, \quad (18)$$

for Sk law.

For large deformations, however, the two laws behave very differently, the NH law being strain-softening whereas the Sk law is strain-hardening [24]. Such large deformations can be expected to occur when dealing with thin capsules since the surface shear modulus $G_s$ decreases with the membrane thickness.

### 3. NUMERICAL METHOD

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

#### 3.1. Discretization of the interface

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

The sphere at rest is meshed by first inscribing an icosahedron (regular polyhedron with 20 triangular faces) in the sphere. The mesh is constructed by placing a new node at the middle of each edge and dividing each element into four new elements; the new nodes are then projected onto the sphere. The procedure is repeated until the desired number of elements is attained [7].
<figure><img src="image_3.png" /><figcaption>Figure 2. Typical capsule shape meshed with flat triangles ($P_1$ elements, $N_E=5120$, $N_N=2562$): (a) at rest and (b) after deformation.</figcaption></figure>

The mesh is shown in Figure 2 for $P_1$ elements. The $P_2$ elements are obtained from the $P_1$ elements by cutting the edges in half in a last step and projecting the new nodes onto the sphere.

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

$$h_1 = \frac{1}{a} \sqrt{\frac{4\pi a^2}{N_E}} \quad \text{and} \quad h_2 = \frac{1}{2a} \sqrt{\frac{4\pi a^2}{N_E}}, \qquad (19)$$

for $P_1$ and $P_2$ elements, respectively. These definitions ensure that the value of $h$ is the same for both element types when the number of nodes of the discretized interface is the same.

## 3.2. Coupling procedure

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

$$\frac{\partial \mathbf{U}}{\partial t} = \mathcal{F}(\mathbf{U}), \qquad (20)$$

where $\mathcal{F}$ is a known function of $\mathbf{U}$ given by the fluid equation (3), in which the traction jump $[\sigma] \cdot \mathbf{n} = \mathbf{q}$ is obtained by solving the membrane problem (9).

We use an explicit second-order Runge-Kutta method to solve this problem. Thus, time convergence can only be expected when the time step $\Delta t$ is such that:

$$\Delta t < \frac{3}{2L}, \qquad (21)$$

where $L$ is the Lipschitz constant of $\mathcal{F}$ [29].

From Equations (3), (7), and (12), we observe that $L$ varies approximately as:

$$L \sim \frac{\Delta x}{\mu} \times \frac{1}{\Delta x} \times \frac{G_s}{\Delta x} = \frac{a}{\Delta x} \frac{\dot{\gamma}}{Ca}. \qquad (22)$$
(Note that in Equation (3), **J** decreases as $1/r$ and the dominant terms in the integral correspond, therefore, to **y** close to **x**. The relevant length scale is then $\Delta x$ rather than $a$). Using the order of magnitude of $L$ found in Equation (22), the stability condition (21) can be rewritten in non-dimensional form as

$$ \dot{\gamma}\Delta t < O(hCa). \qquad (23) $$

It follows that an explicit scheme can be used to study the motion of deformable particles in a Stokes flow, because viscous forces are comparable to elastic ones, i.e. $Ca = O(1)$. Thus the restriction on the time step is governed by the spatial discretization and allows to carry out computations in a reasonable time. Conversely, for very rigid structures, where $Ca \ll 1$, condition (23) leads to very small time steps and unbearably long computation times. This is why explicit schemes are not used in such situations.

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

### 3.3. Finite element procedure

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

For a given element, the shape function associated with node $p$ ($p \in \{1...3\}$ for $P_1$ elements, $p \in \{1...6\}$ for $P_2$ elements) is denoted $N^{(p)}(\xi^1, \xi^2)$, and the value at node $p$ of the $j$th Cartesian component of $\mathbf{v} \in \mathcal{V}_h$ is noted as $v_{Xj}^{(p)}$.

The discretized problem leads to the following matrix system:

$$ [M]\{q\} = \{R\}, \qquad (24) $$

where $\{q\}$ corresponds to the Cartesian components of the nodal values of the discretized load.

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

$$ \sum_{el} \hat{u}_{Xj}^{(p)} \left( \int_{S_{el}} N^{(p)} N^{(q)} dS \right) q_{Xj}^{(q)}, \qquad (25) $$

where the integration is performed on the deformed state. The matrix $[M]$, therefore, has the structure of a mass matrix.

The vector $\{R\}$ corresponds to the right-hand side of Equation (9) and depends non-linearly on the displacement $\mathbf{U}$. The tensions $T^{\alpha\beta}$ are computed directly using Equation (13), after computing
the metric tensor on the deformed geometry of the elements. The virtual strain tensor $\hat{\varepsilon}_{\alpha\beta}$ is related to the covariant representation $\hat{u}_i$ of $\hat{u}$ by

$$ \hat{\varepsilon}_{\alpha\beta} = \frac{1}{2}(\hat{u}_{\alpha,\beta} + \hat{u}_{\beta,\alpha} - 2\Gamma^i_{\alpha\beta}\hat{u}_i), \qquad (26) $$

where $\Gamma^i_{\alpha\beta} = \mathbf{a}_{\alpha,\beta} \cdot \mathbf{a}^i$ are the Christoffel symbols (the curvature tensor $b_{\alpha\beta} = \mathbf{a}_{\alpha,\beta} \cdot \mathbf{n}$ is denoted $\Gamma^3_{\alpha\beta}$ for the sake of brevity). The covariant and Cartesian representations of $\hat{u}$ are related through:

$$ \hat{u}_i = a_i^{Xj} \hat{u}_{Xj}, \qquad (27) $$

where $a_i^{Xj}$ are the Cartesian components of $\mathbf{a}_i$. The virtual strain tensor can, therefore, be related to the nodal Cartesian components of the virtual displacement by

$$ \hat{\varepsilon}_{\alpha\beta} = (\frac{1}{2}N_{,\beta}^{(p)} a_{\alpha}^{Xj} + \frac{1}{2}N_{,\alpha}^{(p)} a_{\beta}^{Xj} + N^{(p)} a_{\alpha,\beta}^{Xj} - \Gamma^i_{\alpha\beta} N^{(p)} a_i^{Xj}) \hat{u}_{Xj}^{(p)}. \qquad (28) $$

To solve the solid problem, we first assemble $\{R\}$ and $[M]$ on the deformed state. In our numerical procedure, the deformed state is known when solving the structural problem and thus the integration on the deformed state is much easier than when using a more standard approach such as in [30]. Surface integration is performed using 6 Hammer points [31] on the elements. Equation (24) is then solved using the sparse solver Pardiso [32, 33].

### 3.4. Fluid solver

Once the load $\mathbf{q}$ is known, the velocity field at the nodes is obtained explicitly from the boundary integral Equation (3), which is discretized on the same mesh as for the solid problem. We use 12 Hammer points for the integration. Note that the kernel $J(\mathbf{x}, \mathbf{y})$ varies as $1/r$ and becomes very large when $\mathbf{y} \to \mathbf{x}$, which can lead to numerical errors. Following [8], we switch to polar coordinates centered on $\mathbf{x}$ when $\mathbf{y}$ and $\mathbf{x}$ belong to the same element. This change of coordinates introduces a Jacobian which goes to 0 as fast as $r$, thereby eliminating the singularity in $J$. In this case, 6 Gauss points are used for integration along each of the polar coordinates.

## 4. VALIDATION

### 4.1. Large isotropic deformation

Before coupling the finite element method to the fluid solver, we first test it when the membrane of the capsule undergoes large static deformations. We consider a spherical capsule inflated from radius $a$ to radius $a_p = (1+\alpha)a$ by an internal pressure $p$. The membrane undergoes an isotropic traction characterized by a stretch ratio $\lambda=1+\alpha$. The isotropic principal tensions can be related to $\lambda$ through the membrane strain energy function [34]. The analytical solution is:

$$ T = G_s \left(1 - \frac{1}{\lambda^6}\right), \qquad (29) $$

for the NH law and

$$ T = G_s(\lambda^2 - 1 + C\lambda^2(\lambda^4 - 1)), \qquad (30) $$
<figure><img src="image_4.png" /><figcaption>Figure 3. Non-dimensional pressure as a function of the inflation factor for the NH and Sk (C = 1) laws. $N_N = 2562$ for both types of elements.</figcaption></figure>

for the Sk law. The analytical expression of the pressure can be deduced using Laplace's law:

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

In order to check the non-linear behavior of the finite element method, we impose a normal displacement field given by

$$U(\mathbf{X}) = \alpha \mathbf{X}. \qquad (32)$$

In Figure 3, we plot the computed non-dimensional pressure $ap/G_s$ as a function of the inflation factor $\alpha$. The results are shown for both laws, using $P_1$ and $P_2$ triangular elements. For $\alpha$ between 0 and 0.5, we find a maximum relative error with respect to the analytical solution of 0.3% with $P_1$ elements and 0.09% with $P_2$ elements. This shows that both element types behave well when submitted to large isotropic deformation.

## 4.2. Convergence

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

The undisturbed velocity field is given by:

$$\mathbf{v}^{\infty}(\mathbf{x}) = \dot{\gamma}x_2\mathbf{e}_1. \qquad (33)$$
In such a flow field, the capsule elongates along a direction contained in the shear plane. The deformed shape of the capsule is almost ellipsoidal and can be approximated by its ellipsoid of inertia. The deformation is then characterized by the Taylor parameter $D_{ij}$:

$$D_{ij} = \frac{|L_i - L_j|}{L_i + L_j} \quad (i, j = 1, 2, 3), \qquad (34)$$

where $L_i$ are the lengths of the principal axes of the ellipsoid of inertia (due to flow symmetry, $L_1$ and $L_2$ are in the shear plane). At $Ca=0.6$ using the NH law, a steady state can be reached. The coefficient $D_{12}^\infty$ represents the Taylor parameter $D_{12}$ in the shear plane at steady state.

We first study the effect of the dimensionless time step $\dot{\gamma}\Delta t$ on the stability of the numerical procedure with a very fine mesh ($N_N=10242$ nodes, mesh size: $h=2.5 \times 10^{-2}$). This mesh size corresponds to $N_E=20480$ $P_1$ elements and to $N_E=5120$ $P_2$ elements. As predicted by Equation (23), the procedure is conditionally stable. The time step stability condition in this particular case is found to be $\dot{\gamma}\Delta t \le 6 \times 10^{-3}$ for $P_1$ elements and $\dot{\gamma}\Delta t \le 5 \times 10^{-3}$ for $P_2$ elements. Below this critical time step, a steady state is reached.

The theoretical value of the deformation is not known, as there is no analytical solution to this problem. We use as a reference value for each element type $P_i$ ($i=1, 2$) the value of $D_{12}^\infty$, denoted $D_i^{\text{ref}}$, obtained with this mesh size for a small time step $\dot{\gamma}\Delta t = 1 \times 10^{-4}$.

We first increase the time step and compute the error $\epsilon_i = |D_{12}^\infty - D_i^{\text{ref}}|/D_i^{\text{ref}}$ relative to the reference value, for each element type. As shown in Figure 4, the numerical procedure converges linearly in time for both $P_1$ and $P_2$ elements. It should be noted that the error remains small ($\epsilon_i < 2 \times 10^{-4}$) for all values of the time step small enough to satisfy Equation (23).

We then study the influence of the mesh size by increasing $h$ with a small time step $\dot{\gamma}\Delta t = 1 \times 10^{-4}$. As shown in Figure 5, the numerical procedure converges with respect to $h$, approximately as $h^4$, and the precision of the results remains good even for cruder meshes ($\epsilon_i < 2 \times 10^{-2}$). Figure 5

<figure><img src="image_5.png" /><figcaption>Figure 4. Relative error $\epsilon_i$ on Taylor's parameter $D_{12}^\infty$ as a function of the time step for $P_1$ and $P_2$ elements at $Ca=0.6$. In all cases, $N_N=10242$, $h=2.5 \times 10^{-2}$ (i.e. 20480 $P_1$ elements, 5120 $P_2$ elements). The reference value $D_i^{\text{ref}}$ for each element type is taken for $\dot{\gamma}\Delta t = 1 \times 10^{-4}$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5. Relative error $\epsilon_i$ on Taylor's parameter $D_{12}^\infty$ as a function of the mesh size for $P_1$ and $P_2$ elements at $Ca=0.6$, and relative difference $\delta$ between the $P_1$ and $P_2$ values. In all cases, $\dot{\gamma}\Delta t = 1 \times 10^{-4}$. The reference value $D_i^{\text{ref}}$ for each element type is taken for $h=2.5 \times 10^{-2}$ (i.e. 20 480 $P_1$ elements, 5120 $P_2$ elements).</figcaption></figure>

also shows that the relative difference $\delta = |D_2 - D_1|/D_2$ between the values obtained with $P_1$ and $P_2$ elements decreases when $h$ is reduced. The results given by the two types of elements thus converge toward a common solution.

The numerical procedure is thus found to be conditionally stable and to converge with respect to the time step and mesh size. In view of the convergence results, all following computations are based on $P_2$ elements with $N_N=2562$ nodes and $N_E=1280$ elements ($h=5 \times 10^{-2}$). For $Ca \ge 0.5$, we use a time step of $\dot{\gamma}\Delta t = 5 \times 10^{-3}$; following Equation (23), the time step is reduced proportionally to $Ca$ for lower values of the capillary number.

### 4.3. Simple shear flow

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

* the work of Lac *et al.* [8], who used the same method for the fluid problem (boundary elements), but a local approach for the solid problem, based on a discretization of the interface by bi-cubic B-splines;

* the work of Li and Sarkar [12], who used the immersed boundary method for the flow and a local approach with a piecewise-constant load for the solid problem;

* the work of Doddi and Bagchi [20], who also used the immersed boundary method for the fluid but finite elements to model the membrane, albeit with a method somewhat different from ours, as noted in the Introduction.

For the external fluid flow, the undisturbed velocity field is given by Equation (33). The membrane behavior is simulated using either the NH law or Sk law ($C=1$).
<figure><img src="image_7.png" /><figcaption>Figure 6. Steady deformed shape of a capsule, following the NH law, for different values of Ca: Ca=0.3 (a), 0.6 (b), and 1.2 (c). The color scale corresponds to the normal component of the load, q·n.</figcaption></figure>

As the flow is started, the capsule elongates while the membrane rotates along the deformed shape (tank treading motion). During the first stages of deformation, folding occurs in the vicinity of the plane orthogonal to the main extension direction (hereafter called the equatorial plane). This is due to the presence of compression tensions in the membrane. Depending on the value of the capillary number Ca, three steady regimes are found, as first described by Lac et al. [8]. For Ca<Ca_L, the membrane is still subjected to compression in the equatorial area and the folds persist at steady state (Figure 6(a)). For larger values of the capillary number, the capsule becomes more elongated by the shear flow; the value of the isotropic part of the tensions (related to the area dilation modulus) increases, which leads to positive tensions at the equator, and the disappearance of the folds (Figure 6(b)). For Ca>Ca_H, the elongated capsule behaves like a slender body submitted to a torque (due to the flow vorticity); negative tensions are again present at steady state and buckling occurs at the tips (Figure 6(c)).

We compare the steady values of $D_{12}^{\infty}$ with the previous studies for the NH law (Figure 7) and the Sk law (Figure 8). We find a good agreement for all values of the capillary number. In particular, our results are always well within 1% of those of Lac et al. [8] for both the NH and Sk laws, which validates the new method. A 5% difference is, however, found with Li and Sarkar [12] for both laws. It is probably related to the crude description of the load used by these authors. The results of Doddi and Bagchi [20] are close to our results at low values of the capillary number but seem to diverge as Ca increases. This is somewhat surprising, as they also use finite elements to model the membrane.

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

$$Ca_L=0.45, \quad Ca_H=0.63, \qquad (35)$$

for the NH law and

$$Ca_L=0.4, \quad Ca_H=2.4, \qquad (36)$$

for the Sk law (C=1). They are equal to those found by Lac et al. [8]. Doddi and Bagchi [20] also found negative tensions at low Ca, but did not quantify the limiting capillary number Ca_L. They did not consider large values of Ca and therefore, did not comment on the existence of Ca_H. The three regimes were observed by Li and Sarkar [12]. They were, however, unable to characterize
<figure><img src="image_8.png" /><figcaption>Figure 7. Values of $D_{12}^{\infty}$ as a function of Ca for a capsule membrane following the NH law in a simple shear flow. Results are compared with [8, 12, 20]. Vertical lines indicate the critical capillary numbers $Ca_L$ and $Ca_H$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Values of $D_{12}^{\infty}$ as a function of Ca for a capsule membrane following the Sk (C=1) law in a simple shear flow. Results are compared with [8, 12]; this case was not studied in [20]. Vertical lines indicate the critical capillary numbers $Ca_L$ and $Ca_H$.</figcaption></figure>

$Ca_L$ and $Ca_H$, as they found negative tensions for all values of Ca, a phenomenon they ascribe to the lack of precision of their numerical method. Hence, our work is the only one beside [8] that determines $Ca_L$ and $Ca_H$; despite using different numerical techniques, we find the same values as previously.

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

### 4.4. Planar hyperbolic flow

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

$$v^{\infty}(\mathbf{x}) = \dot{\gamma}(x_1\mathbf{e}_1 - x_2\mathbf{e}_2). \qquad (37)$$

Two numerical studies of the deformation of a capsule in such a flow exist:

* the previously described work of Lac et al. [8];

* the very recent work of Dodson and Dimitrakopoulos [9], who discretize the interface using spectral elements, on which they solve the local equilibrium equation for the membrane and who compute the velocity using the boundary integral technique.

Since Dodson and Dimitrakopoulos have not yet published the full results of their study, the work of Lac et al. shall be our main reference in this section.

We simulate the membrane behavior using the NH and Sk ($C=1$) laws. Placed in a planar hyperbolic flow, the capsule elongates in the $\mathbf{e}_1$ direction. As in a simple shear flow, compressive tensions and folds appear in the plane orthogonal to the direction of elongation during the transient phase. For values of the capillary number below a critical value $Ca_L$, these folds remain at steady state. For higher values, a steady state free of compressive tensions can be reached. Similarly to the case of the simple shear flow, $Ca_L$ is defined as the value below which negative tensions exist at steady state. As $Ca$ is further increased, we find that the capsule deformation depends on the law used to model the membrane. With the strain-softening NH law, a critical capillary number $Ca_{\infty}$ exists above which the capsule keeps on extending without reaching a steady state. However, this critical capillary number $Ca_{\infty}$ does not exist with the strain-hardening Sk law, for which we always find a stationary deformed shape.

We compare the values we obtain for $D_{12}^{\infty}$ to those of Lac et al. [8] in Figure 9 (NH law) and Figure 10 (Sk law). In the range where they were able to obtain results, we find very good agreement with their values (again well within 1%). Note that the slight discrepancy at low $Ca$ for the NH law is due to the presence of strong folds in the shape that makes it impossible to define $D_{12}^{\infty}$ accurately.
<figure><img src="image_10.png" /><figcaption>Figure 9. Values of $D_{12}^{\infty}$ as a function of $Ca$ for a capsule membrane following the NH law in a planar hyperbolic flow. Results are compared with [8]. Vertical lines indicate the critical capillary numbers $Ca_L$ and $Ca_{\infty}$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10. Values of $D_{12}^{\infty}$ as a function of $Ca$ for a capsule membrane following the Sk (C=1) law in a planar hyperbolic flow. Results are compared with [8]. The vertical line indicates the critical capillary number $Ca_L$.</figcaption></figure>

As in the case of the simple shear flow, the critical capillary numbers between the three regimes remain an open problem. With our numerical method, the critical values that we find for the NH law are

$$Ca_L=0.14, \quad Ca_{\infty}=0.35, \qquad (38)$$
<figure><img src="image_12.png" /><figcaption>Figure 11. Transient (γt = 1.5) and stationary (γt = 10) profiles in the (x₁, x₂) plane of a capsule in a hyperbolic flow at Ca = 1.0. The capsule membrane follows the Sk law (C = 1). The curvature at the tips is higher during the transient phase than at steady state.</figcaption></figure>

and for the Sk law (C = 1)

$$Ca_L = 0.15, \quad Ca_\infty = +\infty. \qquad (39)$$

We find the same values for $Ca_L$ as Lac et al. [8]. However, the results differ for $Ca_\infty$: Lac et al. found a lower critical value for the NH law ($Ca_\infty$ = 0.22) and a finite value ($Ca_\infty$ = 0.7) for the Sk law. An explanation for this discrepancy can be found in the work of Dodson and Dimitrakopoulos [9], who have shown that capsules in a hyperbolic flow exhibit very high curvatures at the tips during the transient phase, before returning to smaller tip curvatures at steady state. We observe this phenomenon (as shown in Figure 11), although our discretization does not allow a precise quantification of the curvature at the tips. We believe that Lac et al. numerical method became unstable because of the high curvatures and that they misinterpreted the numerical failure as a sign of unbounded extension, which indeed occurs at higher $Ca$ with certain strain-softening materials, but not with strain-hardening ones.

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

# 5. DISCUSSION AND CONCLUSION

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

## 5.1. Limitations of the local approach

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

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

Another technique to mesh the sphere while using a local approach is to use elements between which the first spatial derivatives are discontinuous. An *ad hoc* method then needs to be implemented to treat the discontinuities. An example of such a method is provided by the work of Wang and Dimitrakopoulos [37]. They discretized the interface using a moderate number of spectral elements (between 6 and 14), which offer a highly regular representation of the position within each element, but have discontinuities between the elements. They added a first-order smoothing scheme to enforce the continuity of the position and of its first derivatives from one element to the next. In contrast with the low number of elements (and, therefore, limited amount of smoothing necessary) used in spectral methods, Ramanujan and Pozrikidis [7] used an unstructured mesh with 512 triangular $P_2$ elements, between which the tensions are discontinuous. They developed an intricate technique to average the tensor $T$ at each node from the values on the connected elements and then re-interpolated it between the nodes to compute the load (which they treated as element-wise constant). The effects of so many averages are difficult to predict or control.

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

## 5.2. Comparison with other methods

We have compared the results provided by our new model in the case a simple shear flow with three previous studies and have found a good agreement. It may be noted, in particular, that our results are remarkably close to those of Lac *et al.* [8] (within 1%). Both studies find the same critical capillary numbers $Ca_L$ and $Ca_H$, although two different methods are used to model the behavior of the membrane. This seems to indicate that the values of the critical capillary numbers have a physical relevance and are not a consequence of the numerical method used.

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

### 5.3. Bending stiffness and limits of the membrane model

A major difference between our results and those by Lac *et al.* [8] is the behavior of the numerical method when the membrane is undergoing in-plane compression. Whereas the technique used by Lac *et al.* eventually fails when negative tensions appear ($Ca < Ca_L$ or $Ca > Ca_H$ in a simple shear flow, $Ca < Ca_L$ in a planar hyperbolic flow), our numerical method remains stable and a steady equilibrium state is achieved. This is probably due to the larger numerical stiffness of the finite elements as compared to the bi-cubic B-spline functions used by Lac *et al.* Neither Doddi and Bagchi [20] nor Li and Sarkar [12], who also used lower order discretization than Lac *et al.*, reported stability problems at low $Ca$, which would confirm the hypothesis. The stiffness introduced by the numerical method contributes to the stability of the problem because it enriches the membrane model with some bending stiffness. It allows the numerical procedure to remain stable during transient phases when in-plane compression and high curvatures may render other methods unstable. It also has the advantage that steady states with negative tensions can be computed and studied.

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

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

## ACKNOWLEDGEMENTS

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

1. Kühtreiber WM, Lanza RP, Chick WL. *Cell Encapsulation Technology and Therapeutics*. Birkhäuser: Basel, 1998.

2. Chang TMS. Future prospects for artificial blood. *Trends in Biotechnology* 1999; **17**(2):61-67. DOI: 10.1016/S0167-7799(98)01242-6.

3. Risso F, Collé-Paillot F, Zagzoule M. Experimental investigation of a bioartificial capsule flowing in a narrow tube. *Journal of Fluid Mechanics* 2006; **547**:149-173.

4. Lefebvre Y, Leclerc E, Barthes-Biesel D, Walter J, Edwards-Levy F. Flow of artificial microcapsules in microfluidic channels: a method for determining the elastic properties of the membrane. *Physics of Fluids* 2008; **20**(12):123102-1-123102-10. DOI: 10.1063/1.3054128.

5. Barthès-Biesel D, Rallison JM. The time-dependent deformation of a capsule freely suspended in a linear shear flow. *Journal of Fluid Mechanics* 1981; **113**:251-267.

6. Li XZ, Barthès-Biesel D, Helmy A. Large deformations and burst of a capsule freely suspended in an elongational flow. *Journal of Fluid Mechanics* 1988; **187**:179-196.

7. Ramanujan S, Pozrikidis C. Deformation of liquid capsules enclosed by elastic membranes in simple shear flow: large deformations and the effect of capsule viscosity. *Journal of Fluid Mechanics* 1998; **361**:117-143.

8. Lac E, Barthès-Biesel D, Pelekasis NA, Tsamopoulos J. Spherical capsules in three-dimensional unbounded Stokes flow: effect of the membrane constitutive law and onset of buckling. *Journal of Fluid Mechanics* 2004; **516**:303-334.

9. Dodson III WR, Dimitrakopoulos P. Spindles cusps and bifurcation for capsules in Stokes flow. *Physical Review Letters* 2008; **101**(20):208102-1-208102-10. DOI: 10.1103/PhysRevLett.101.208102.

10. Eggleton CD, Popel AS. Large deformation of red blood cell ghosts in a simple shear flow. *Physics of Fluids* 1998; **10**:1834-1845.

11. Doddi S, Bagchi P. Effect of inertia on the hydrodynamic interaction between two liquid capsules in simple shear flow. *International Journal of Multiphase Flow* 2008; **34**:375-392.

12. Li X, Sarkar K. Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane. *Journal of Computational Physics* 2008; **227**(10):4998-5018. DOI: 10.1016/j.jcp.2008.01.034.

13. Peskin CS. The immersed boundary method. *Acta Numerica* 2002; **11**:1-39.

14. Mittal R, Iaccarino G. Immersed boundary methods. *Annual Review of Fluid Mechanics* 2005; **37**:239-261.

15. Pozrikidis C. The axisymmetric deformation of a red blood cell in uniaxial straining Stokes flow. *Journal of Fluid Mechanics* 1990; **216**:231-254. DOI: 10.1017/S0022112090000416.

16. Leyrat-Maurin A, Barthès-Biesel D. Motion of a deformable capsule through a hyperbolic constriction. *Journal of Fluid Mechanics* 1994; **279**:135-163.

17. Quéguiner C, Barthès-Biesel D. Axisymmetric motion of capsules through cylindrical channels. *Journal of Fluid Mechanics* 1997; **348**:349-376.

18. Pozrikidis C. Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow. *Journal of Fluid Mechanics* 1995; **297**:123-152.

19. Lac E, Barthès-Biesel D. Deformation of a capsule in simple shear flow: effect of membrane prestress. *Physics of Fluids* 2005; **17**:072105-1-072105-8.

20. Doddi S, Bagchi P. Lateral migration of a capsule in a plane Poiseuille flow in a channel. *International Journal of Multiphase Flow* 2008; **34**(10):966-986. DOI: 10.1016/j.ijmultiphaseflow.2008.03.002.

21. Charrier JM, Shrivastava S, Wu R. Free and constrained inflation of elastic membranes in relation to thermoforming—non-axisymmetric problems. *Journal of Strain Analysis for Engineering Design* 1989; **24**(2): 55-74.

22. Rallison JM, Acrivos A. A numerical study of the deformation and burst of a viscous drop in extensional flow. *Journal of Fluid Mechanics* 1978; **89**:191-200.

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

24. Skalak R, Tozeren A, Zarda RP, Chien S. Strain energy function of red blood cell membranes. *Biophysics Journal* 1973; **13**:245-264.

25. Green AE, Adkins JE. *Large Elastic Deformations* (2nd edn). Oxford University Press: Oxford, 1970.

26. Crisfield MA. *Non-Linear Finite Element Analysis of Solids and Structures: Advanced Topics*. Wiley: New York, NY, U.S.A., 1997.

27. Oden JT. *Finite Elements of Non-Linear Continua*. McGraw-Hill: New York, 1972.
28. Carin M, Barthès-Biesel D, Edwards-Lévy F, Postel C, Andrei D. Compression of biocompatible liquid-filled HSA-alginate capsules: determination of the membrane mechanical properties. *Biotechnology and Bioengineering* 2003; **82**:207-212.

29. Stuart AM, Humphries AR. *Dynamical Systems and Numerical Analysis*. Cambridge University Press: Cambridge, 1998.

30. Wriggers P, Taylor R. A fully non-linear axisymmetrical membrane element for rubber-like materials. *Engineering Computations* 1990; **7**:303–310.

31. Hammer P, Marlowe O, Stroud A. Numerical integration over simplexes and cones. *Mathematical Tables and Other Aids to Computation* 1956; **10**(55):130–137.

32. Schenk O, Gärtner K. Solving unsymmetric sparse systems of linear equations with PARDISO. *Future Generation Computer Systems* 2004; **20**(3):475–487. DOI: 10.1016/j.future.2003.07.011.

33. Schenk O, Gärtner K. On fast factorization pivoting methods for sparse symmetric indefinite systems. *Electronic Transactions on Numerical Analysis* 2006; **23**:158–179.

34. Barthès-Biesel D, Diaz A, Dhenin E. Effect of constitutive laws for two dimensional membranes on flow-induced capsule deformation. *Journal of Fluid Mechanics* 2002; **460**:211–222.

35. Milnor J. Analytic proofs of the 'Hairy ball theorem' and the Brouwer fixed point theorem. *American Mathematical Monthly* 1978; **85**(7):521–524.

36. Lac É, Morel A, Barthès-Biesel D. Hydrodynamic interaction between two identical capsules in a simple shear flow. *Journal of Fluid Mechanics* 2007; **573**:149–169.

37. Wang Y, Dimitrakopoulos P. A three-dimensional spectral boundary element algorithm for interfacial dynamics in Stokes flow. *Physics of Fluids* 2006; **18**(8):82–106. DOI: 10.1063/1.2337572.