Axisymmetric Boundary Element Method for vesicles in a capillary

Journal of Computational Physics 289 (2015) 62–82

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Axisymmetric Boundary Element Method for vesicles in a capillary R. Trozzo a , G. Boedec b , M. Leonetti b , M. Jaeger a,∗ a b

Aix-Marseille Université, CNRS, Centrale Marseille, M2P2 UMR 7340, 13451, Marseille, France Aix-Marseille Université, CNRS, Centrale Marseille, IRPHE UMR 7342, 13384, Marseille, France

a r t i c l e

i n f o

Article history: Received 7 May 2014 Received in revised form 15 December 2014 Accepted 15 February 2015 Available online 23 February 2015 Keywords: Boundary Element Method Fluid-structure interaction Vesicle Poiseuille flow Computational model

a b s t r a c t The problem of a vesicle transported by a fluid flow can present a large range of length scales. One example is the case of a vesicle producing a tether, and eventually pearls, in an elongational flow. Another case occurs when a lubrication film is formed, such as during the short range interaction between two vesicles. Such problems are still challenging for 3D simulations. On the other hand, a good understanding could be obtained by first considering the axisymmetric regime when such a regime exists. An axisymmetric model could then be used, without the criticisms that can be made of a 2D approach. We propose such a model, primarily interested in flows through narrow cylindrical capillaries. Two options are compared, with and without explicit representation of the capillary boundaries by a mesh. The numerical effort is characterized as a function of the vesicle’s initial shape, the flow magnitude and the confinement. The model is able to treat typical configurations of red blood cells flowing through very narrow pores with extremely thin lubrication films. © 2015 Elsevier Inc. All rights reserved.

1. Introduction In studies focusing on giant vesicles in fluid flows, numerical models built with the boundary integral method (BIM) have proven to be the most accurate, as long as the Stokes flow regime and Newtonian fluids are concerned. 3D BIM-models have been developed and applied to many problems concerning the dynamics of red blood cells and their biomimetic counterparts in fluid flows; see [1–3] and references therein. However, some situations involving more than one characteristic length scale can become very challenging for computations in 3D, whereas they could be treated much more efficiently with an axisymmetric model for all axisymmetric regimes. The paradigm of this situation is the emergence of a tether in elongational flows [4,5] or sedimentation [6,7], a tube of submicroscopic size and huge length, almost a factor 100 with respect to the characteristic size. For instance, if the mother vesicle has a unit size, the tube has a radius of the order of 1/10, while its length may exceed 10 [7]. Another noticeable situation concerns a vesicle interacting with a wall. Experimental results on adhesion to a substrate show that the thickness of the lubrication film between the vesicle and the wall is of the order of 1/200 with respect to the vesicle typical size [8]. Although problems of this kind are often fully 3D, like the dynamics of a vesicle in a shear flow near a flat wall [9], some can present axial symmetry in a large number of flow regimes. A first example concerns deposition and spreading of a vesicle on a substrate as the final stage of sedimentation [10] or the unbinding of a vesicle from a substrate by pipette

*

Corresponding author. E-mail addresses: [email protected] (R. Trozzo), [email protected] (G. Boedec), [email protected] (M. Leonetti), [email protected] (M. Jaeger). http://dx.doi.org/10.1016/j.jcp.2015.02.022 0021-9991/© 2015 Elsevier Inc. All rights reserved.

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

63

suction [8]. A last example concerns vesicle motion in a capillary. In this case, for really deflated vesicles and large flow velocities, the lubrication film is not optically measurable and therefore strictly lower than 1/10 of the vesicle size. Depending on the value of the operating parameters, in all these examples the computation power needed can become significant. Indeed, a crude rule says that the lubrication film is properly described when the distance between two adjacent computational nodes is smaller than the thickness of the film. So, if one wishes to accurately determine the dynamics in the lubrication film, it is necessary to greatly increase the number of nodes and thus the computational cost. Of course, some 3D computations can be done, but an exhaustive study of the problem generally involves a large number of computations, to build at least a vesicle shape diagram for the considered flow. Moreover, these axisymmetric flow configurations can be of interest to develop mechanical properties characterizing methods [11]. In this case, the computations must be implemented with high efficiency. Here we will focus on the motion of a vesicle in a capillary. The understanding of vesicle behavior in confined capillary flows would lead to a better knowledge of red blood cells (RBCs) in human capillaries. Although RBCs behave differently from vesicles due to the shear elasticity of their plasma membrane, they also share some similar mechanical behaviors with vesicles, especially when capillary size is comparable to cell diameter [12–18]. For what concerns boundary integral approaches, to our knowledge most studies devoted to this problem are either 2D [19], either restricted to unbounded Poiseuille flow [20,21] and sometimes both [22–25]. Experimental studies have shown that the most commonly stationary vesicle’s shapes are axisymmetric bullet and parachute like shapes, until these shapes destabilize to 3D shapes like slippers [26,27]. These shapes have also been found in the few existing 3D simulations [21,27,28]. However for an exhaustive parametric study, some approximations to the real flow must be introduced to render the problem tractable: either stay in 3D and consider the related unbounded Poiseuille flow or drift to 2D [29,30]. Both have proven to be very useful, but none of these approaches can fully take into account the action of a capillary confinement on vesicle’s shapes, especially when confinement becomes significant. So, in other words, an efficient BIM axisymmetric model to study the transport of vesicles in a capillary would fill a miss in the class of BIM models. In fact axisymmetric BIM models have not attracted so much attention at all. The reason is that their development is more involved than for a 2D formulation and limited compared to a full 3D model. Considering only the case of soft deformable particles, BIM axisymmetric models have been developed first to study axisymmetric configurations of drops [31]. Later, the axisymmetric Green function for free-space flows has been used to develop an axisymmetric formulation for capsules in elongational flows in [32] and used by the same group to study extensively the deformation of capsules flowing through pores [11,33,34]. Its however only recently that axisymmetric models have been proposed for vesicles in free space flows [5,35]. An axisymmetric Green function taking into account the no-slip condition at the wall have been proposed in [36] for the case of a flat wall perpendicular to the flow, and in [37] for a circular pipe. Only the periodic version of this last has been used in practice to study the motion of a train of drops or red blood cells in a capillary [37,38]. To our knowledge no BIM-model for an isolated capsule or vesicle in an infinite capillary has been proposed so far. Based on the 3D Boundary Element Method (BEM) of [39], we propose and discuss its axisymmetric version here. In the frame of BEM, there are two ways to compute the flow confined in a capillary, either use a Green’s function that takes into account the confinement action of the channel wall, or use the free space Green’s function and mesh the wall. Each of these two approaches presents advantages and drawbacks. The first needs more development efforts but simplifies the job at the simulation phase by eliminating all questions concerning inflow and outflow boundary conditions as well as the mesh size for the wall. These questions must be answered for each new set of operating parameters for the second approach, but development efforts are limited. Moreover the second approach is more versatile since it is not limited to capillaries of constant section. Both approaches will be considered here and compared. This comparison will help us to define guides for the axisymmetric computation of vesicles in capillary flows. The paper is organized as follows. In Section 2 we first briefly recall the principle of the 3D model of [39]. We present then the axisymmetric extension for free space cases before taking into account the confined geometry, with particular attention to the numerical integration difficulties. In Section 3 the axisymmetric model is validated first for free space configurations and then for capillary flows. Finally the model is used to address some open questions concerning the influence of confinement on the stationary shape taken by a vesicle in a capillary flow. 2. Model 2.1. 3D BEM formulation in free space 2.1.1. Physical problem The description of vesicles in fluid flows is a fluid-structure interaction problem. Because of the typical microscopic scale, the assumption of Stokes flow regime is usually justified. Thus the flow of ambient fluids, internal (i) and external (e), is governed by Stokes equations:

∇ · σ i ,e + ρ i ,e g = 0 ∇·u

i ,e

=0

with the velocity continuity condition at the interface  (vesicle’s membrane):

(1) (2)

64

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

u e (x) = u i (x) = u 

(3)

i ,e

In Eq. (1) σ represents the stress tensor in the inside and outside bulk fluid, with inner and outer viscosities being supposed to be the same (η = ηi = ηe ). In case of non-zero density contrast between inside and outside, the contribution of buoyancy can be included into a modified pressure term p¯ i /e = p i /e + ρ i /e gz, such that the Stokes equations become:

∇ · σ¯ i ,e = 0 , ∇ · u i ,e = 0 with

(4)



σ¯ i,e = − p¯ i/e I + η ∇ u i/e + ∇ T u i/e

 (5)

The mechanical action of the phospholipid membrane on the surrounding fluids is accounted for by:

• A surface incompressibility condition ∇ s · u  = 0, leading to a surface tension γ and thus to a tension surface force density, which is purely normal only if γ takes a constant value on  : f γ = −2γ Hn + ∇ s γ

(6)

• A resistance to bending specified by a bending modulus κ

∼ 10−19 J (a typical value for vesicles and RBC [40]), leading

to a bending surface force density, which is purely normal:





f b = κ 2s H + 4H ( H 2 − K ) n

(7)

The total action of the membrane is given by the membrane surface force density f = f b + f γ in which n is the normal vector pointing out of the vesicle, H and K the mean and Gaussian curvatures, ∇ s and s the gradient and Laplace–Beltrami operators. Since the membrane is assumed to be in quasi-static mechanical equilibrium, the membrane force density f is balanced by the hydrodynamic stress action:

[[σ ]] · n + f = [[σ¯ ]] · n + f + f g = 0

(8)

[[σ¯ ]] · n = (σ¯ e − σ¯ i ) · n , f g = ρ gzn

(9)

with

where f g is due to the jump in pressure resulting from the density difference ρ = ρ i − ρ e . The problem can be expressed in terms of dimensionless variables (the tilde indicates a dimensionless variable):

x = R ref x˜ , t = t ref t˜ , u = R ref /t ref u˜ , f = f ref ˜f ,

γ = f ref R ref γ˜

(10)

using the following reference quantities:

• lref = R 0 , radius of the sphere of same enclosed volume; • t ref = η R 30 /κ viscous damping time of a bending disturbance; • f ref = κ / R 30 typical bending density forces; In the following, only dimensionless variables are used, and the tilde is omitted for clarity. The first characteristic dimensionless number of the problem is the excess area , measuring the additional area compared to the sphere of same volume:

A = (4π + ) R 20

(11)

Because of the surface area conservation constraint, a spherical vesicle ( = 0) with fixed volume and surface will act as a rigid sphere. In order to deform the vesicle, it must be “deflated”, and this √ deflation is measured by the parameter . Equivalently, instead of the excess area, the reduced volume defined as v = 3 4π V A −3/2 can be used. It defines the vesicle deflation as the ratio of the actual enclosed volume V over the volume of a sphere of same surface area A. 2.1.2. 3D BEM formulation We first briefly recall the principles of the 3D model of [39], which will serve as a base for the development of the axisymmetric model. For details see [39]. Owing to the linearity of Stokes equation, the velocity can be decomposed into an external flow field u ∞ and a perturbation flow field due to the force density f generated by the membrane. For an interface embedded in a viscous fluid of dynamic viscosity η , the ith component of the velocity field at any point x0 of the membrane is entirely determined by the force distribution on the interface [37]:

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

u i ( x0 ) = u ∞ i ( x0 ) +

1

65

 G i j (x0 , x) f j (x)dS (x)

8πη

(12)



where G i j are the components of the Stokeslet:

G i j ( x, x0 ) =

δi j r

+

xˆi xˆj

(13)

r3

with xˆ = x − x0 and r = ˆx, and δi j is the Kronecker symbol. For convenience, the problem is often solved in the vesicle’s referential frame. In that case, the imposed background velocity field u ∞ is modified by subtracting the vesicle’s mean velocity (velocity of the center of gravity). This last is a function of the vesicle’s shape and must thus be adjusted at each time step until the stationary regime is reached, if such a regime exists. Discretizing this relation over a surface mesh consisting of a set of planar triangular elements, the problem can be expressed as a set of linear algebraic relations between nodal velocities and nodal surface force densities:

⎛ . ⎞ ⎛ ⎞ ⎛ .. .. ⎜ ⎟ ⎜ . ⎟ ⎜ ⎜ un ⎟ ⎜ ∞ n ⎟ ⎜ i ⎟ = ⎝ (u )i ⎠ + G ⎜ ⎝ ⎝ ⎠ .. .. . .

⎞ .. . ⎟ f in ⎟ ⎠ .. .

(14)

where n is the node number in the triangulation, and G is a matrix resulting from the numerical integration of the integral in (12). We can further split the one column matrix containing the nodal values of the membrane forces in bending force and surface tension contributions, leading to the matrix operators F b acting on the vector containing the nodal positions and F γ acting on the γ nodal values [39]:

⎞ ⎛ ⎞ .. .. . . ⎟ ⎜ n⎟ ⎜ ⎜ f ⎟ = F b ⎜ xn ⎟ + F γ ⎝ i ⎠ ⎝ i ⎠ .. .. . . ⎛

where the tension form it reads:

⎞ .. . ⎜ n⎟ ⎜γ ⎟ ⎠ ⎝ .. . ⎛

(15)

γ is such that the total velocity field (14) satisfies the surface incompressibility constraint. In matricial

⎞ ⎞ ⎡⎛ ⎛ ⎞⎤ .. .. .. . . . ⎟⎥ ⎟ ⎟ ⎜ ⎢ ⎜ ⎜ n ⎟ ⎢⎜ (u ∞ )n ⎟ + G F b ⎜ xn ⎟⎥ = 0 γ DG F γ ⎜ D + i ⎠ ⎠ ⎝ ⎣⎝ ⎝ i ⎠⎦ .. .. .. . . . ⎛

(16)

where the matrix operator D defines the surface divergence operator applied to the velocity nodal values. Inverting this expression in order to eliminate γ , the velocity nodal values on the interface are finally given by:

⎞ ⎞ ⎛ ⎞⎤ ⎡⎛ .. .. .. . . ⎟⎥ .   ⎟ ⎜ n⎟ ⎜ ⎢ ⎜   − 1 ⎜ u ⎟ = Id − G F γ D G F γ ⎢⎜ (u ∞ )n ⎟ + G F b ⎜ xn ⎟⎥ D i ⎠ ⎝ i ⎠ ⎝ i ⎠⎦ ⎣⎝ .. .. .. . . . ⎛

(17)

In order to discretize in time the equation of the interface evolution, an explicit time stepping algorithm would result in severe constraints on the time step, growing with increasing resolution of the mesh. It is then necessary to use a semi-implicit time scheme to keep computations at a reasonable cost, computing stiff bending forces at the updated position x(t + dt ). Using a semi-implicit time stepping algorithm also improves long term stability of the algorithm: spurious modes which could lead to shape degeneracy are correctly damped by the implicit treatment of the bending forces. The semi-implicit time scheme that we have developed in [39] for 3D computations can be used for the axisymmetric model as well. Geometric properties such as the normal vector and mean curvature involve first and second order derivatives of the position, and direct computation is thus prohibited on a C 0 representation of the interface. However, this difficulty has been overcome in [39] by using discrete estimates of the operator acting on a smooth surface. The computation of the curvature vector (mean curvature times the normal vector) follows the approach of [41] and the principle has been extended to the computation of the bending forces that involves in turn the second derivative of the mean curvature. The surface Laplacian for a scalar field f on each node n of the triangulated surface mesh is expressed in terms of the flux of the surface gradient of f on the contour ∂ S of a small patch S of surface surrounding the node:

66

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

Fig. 1. Integration contour ∂ S (bold line) used for 3D computations. The gray zone is the patch attributed to the node n. Vector the boundaries.



1

<s f >n =

∇ s f · ν dl =

S

 1 

∂S

S

ν is the normal vector to

∇ s f e · ν e dl

(18)

e ∈ E n∂ S e

E n is the subset of elements of the C 0 triangulation connected to the node n, ∂ S e is the restriction of ∂ S on the element e and ν e is the tangent vector to the surface, normal to ∂ S e . In [39] contours ∂ S made of segments passing through the midpoints of the element’s edges have been used (see Fig. 1). With triangulated surfaces, surface parametrization is available only piecewise in the form of element local coordinates. Therefore, the computation of a surface gradient is performed with global Cartesian components of vectors and with derivatives with respect to element local coordinates. The surface gradient of a scalar field f is then given on each element by:

∂ f e ∂ xe ∂ξ α ∂ξ β

∇ s f e = aα β

(19)

∂x ∂x where ξ α (α = 1, 2) are the element coordinates, and aα β = ∂ξ α · ∂ξ β are the components of the inverse of the corresponding local metric. Moreover, with a C 0 piecewise linear approximation, all derivatives with respect to element coordinates take constant values, thus leading to constant values of ∇ s f e over each element e. In this way the i-th component of the curvature vector and the surface Laplacian of the mean curvature are computed by:

< Hni >n =

 1  S

∇ s xei · ν e dle

(20)

e∈ E n

∂ Se

<s H >n =

 1  S

∇ s H e · ν e dle

(21)

e∈ E n

∂ Se

where xei and H e are linear interpolations on the element of the corresponding nodal values of xi and H . Concerning the Gaussian curvature K , it is computed using a discrete version of the Gauss–Bonnet theorem [41] in the 3D model, but it will be computed in a different way in the axisymmetric one, see the end of Section 2.2.2 for details. Finally the expression of nodal bending forces as a function of nodal coordinates can be summarized by the matrix form used in (15):

⎞ ⎛ ⎞ .. .. . . ⎟ ⎜ n⎟ ⎜ ⎜ f ⎟ = F b ⎜ xn ⎟ ⎝ i ⎠ ⎝ i ⎠ .. .. . . ⎛

(22)

Similarly, the surface divergence of the velocity vector field u on an element e is also constant, given by:

∇ s · u e = aα β

∂ xe ∂ u e ∂ξ α ∂ξ β

(23)

and the approximated nodal values are simply obtained by area weighted average:

<∇ s · u >n =

 Ae e∈ E n

A

(∇ s · u )e

(24)

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

where A e is the area of the element e and A = matrix form used in (16):



67

A e . Thus the expression of nodal velocity divergence can be cast in the

⎞ ⎛ ⎞ .. .. . ⎟ ⎜ ⎜ .n ⎟ n ⎜ (∇ s · u ) ⎟ = D ⎜ u ⎟ ⎠ ⎝ ⎝ i ⎠ .. .. . . ⎛

(25)

The F γ matrix giving the nodal surface forces due to the surface tension

 Ae

<∇ s γ >n = (Id − n ⊗ n)

A

e∈ E n

γ is built in the same way, with:

(∇ s γ )e

(26)

where the projection operator ensures projections of the membrane tension force on the surface tangent plane at the node. 2.2. Axisymmetric BEM formulation in free space 2.2.1. Axisymmetric Green function in free space In this section we develop simplified versions of the boundary integral equation for axisymmetric flow configurations. Thus none of the variables is a function of the azimuthal angle φ . Moreover we will only consider axisymmetric flow with no swirling motion, thus velocity and force vector fields do not have any azimuthal component. The integral over  in (12) can then be explicitly integrated in the azimuthal angle φ , reducing the 3D problem to a 2D one, namely the determination of the radial and axial components of the velocity u r (r , z) and u z (r , z) in an azimuthal plane. Working with the xz plane, we have the correspondence r with x, u r with u x , y = u y = 0 and dS = r dφ dl, where dl is the differential arc length of the trace of the boundary  on the azimuthal plane. We then apply Eq. (12) at a point x0 in the xz azimuthal plane and obtain: ∞

u α ( x0 ) = u α ( x0 ) +



1

M α β (x0 , x) f β (x)dl(x)

8πη

(27)



where x0 = r0 e r + z0 e z , x = re r + ze z . The Greek subscripts α and β are either r or z, indicating the radial and axial components respectively. The term u ∞ α is the α component of the imposed external velocity field.  is now the curve of the membrane interface in the plane (r , z), dl the differential arc length on  , f β the β component of viscous strain jump at the interface and M α β is the azimuthal angle integration of the 3D Green function:

 M ( x0 , x) =

M rr M rz M zr M zz



2π =r

G xx cos φ + G yx sin φ G xz cos φ + G yz sin φ

G zx G zz

 dφ

(28)

0

Following the approach proposed by [37], the explicit expression of the matrix M is:



M =r

I 11 + (r 2 + r02 ) I 31 − rr0 ( I 30 + I 32 ) zˆ (r I 30 − r0 I 31 )

zˆ (r I 31 − r0 I 30 ) I 10 + zˆ 2 I 30



(29)

with zˆ = z − z0 and:

2π I mn (ˆz, r , r0 ) ≡ 0

cosn θ

[ˆz2

+ r2

+ r02 − 2rr0 cos θ]m/2

dθ

(30)

The I mn can be expressed in terms of complete elliptic integrals of the first and second kinds [37] and these elliptic integrals can be computed using polynomial approximations [42]. 2.2.2. Axisymmetric vesicle model In the axisymmetric case the membrane is discretized by C 0 piecewise linear 2D elements, leading to the following approximation of any scalar field f on each element:

f e (ξ ) =



N i (ξ ) f ie

(31)

i

where i is the i-th element node, f ie the value of f at node i and

N 1 (ξ ) =

1−ξ 2

N 2 (ξ ) =

1+ξ 2

(32)

68

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

are the 1D linear Lagrange interpolation functions, with the local coordinate ξ ∈ [−1, +1]. The derivative of f e (ξ ) with respect to ξ then takes a constant value on each element e:

df e dξ

=

f 2e − f 1e

(33)

2

The differential surface operators are computed on each element with the element local parametrization ξ1 = ξ along  and the axisymmetric natural global parametrization ξ2 = φ for the azimuthal direction. The components of the local metric are thus given by:

2  e 2  e 2  e 2  e r − r1e z − ze1 ∂ xe ∂ xe dr dz + = 2 + 2 · = ∂ξ ∂ξ dξ dξ 2 2 ∂ xe ∂ xe  e 2 a22 = · = r ∂φ ∂φ a12 = a21 = 0 a11 =

where:

j=



√

r2e − r1e

a11 =

2

2  + ze2 − ze1 2

=

(34)

Le

(35)

2

is the Jacobian of the transformation from the reference element to the real one, which is the ratio of the real element’s length L e to the reference element’s length. The surface gradient of f on an element e is then given by:

  df e dr e dze er + ez = a11 dξ ∂ξ dξ dξ dξ  e  e  e e e e f − f r − r z −2 2 1 2 1 2 − z1 =j er + ez 2 2 2     e f 2 − f 1e r2e − r1e ze − ze = er + 2 e 1 e z e e L L L   e f 2 − f 1e = νe e

∇ s f e = a11

df e ∂ xe

(36)

L

where ν e is the unit vector in the direction of the element e (see Fig. 2). The surface divergence of the velocity is:

∂ xe ∂ u e ∂ xe ∂ u e · + a22 · ∂ξ ∂ξ ∂φ ∂φ   e e dr du r dze du ez ue + er + = j −2 dξ dξ dψ dψ r  e   e  e   e r2 − r1e u r 2 − u r e1 z2 − ze1 u z 2 − u z e1 u re + + = e e e e e

∇ s · u e = a11

L

L

L

L

(37)

r

Note that in the axisymmetric case, the surface divergence of the velocity is not constant on an element because of the term u re /r e . However this term doesn’t introduce any singularity since u r = 0 if r = 0. So this term is simply set to zero for pole nodes. In order to compute the bending forces with the same approach as in 3D, the small patch to consider for a node n is now the ring formed by the rotation of the two segments that link the node n to the middle of the neighboring elements, involving the nodes (n − 1, n) and (n, n + 1), as shown in Fig. 2. Thus the integration contour ∂ S is made of the two circles ∂ S − and ∂ S + forming the boundary of the patch. If the node is on a pole (r = 0), then the patch is the cap coming from the rotation of the segment formed by the node and the middle of the nearest element. Moreover the border is symmetrized by introducing an additional fictitious element consisting of the node at the pole and the symmetric (r → −r) of the next node. The surface Laplacian of f at node n is then given by:

⎛

< s f >n = =

1 S 1 S

⎝ 



 ∇ s f · ν dl +

∂ S+

2π



⎞

∇ s f · ν dl⎠

∂ S−

r n + r n +1 2



+

∇ s f e · ν e + − 2π



r n −1 + r n 2



−

∇s f e · ν e

−



R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

69

Fig. 2. Integration contour used for axisymmetric computations. The gray zone is the patch attributed to the node n.

where e − is the element joining nodes n − 1 and n, whereas e + is the one joining nodes n and n + 1. Using (36), this gives:

<s f >n =

π



(rn+1 + rn )( f n+1 − f n ) Le

S

with:

S = S+ + S− =

π

+

+

−

(rn + rn−1 )( f n − f n−1 ) Le



−

(38)



−

L e (r n+1 + 3r n ) + L e (3r n + r n−1 )

4

As in the 3D case, the approximated value of the divergence of the velocity at node n is obtained by area weighted average:

<∇ s · u >n =

1 A

+

+

−

A e (∇ s · u )e + A e (∇ s · u )e

−

 (39)

with: +

−

+

−

A = A e + A e = π L e (r n+1 + r n ) + π L e (r n + r n−1 )

(40)

A simpler length weighted average can also be used without noticeable difference on the results, thus with: +

−

+

A = Ae + Ae = Le + Le

−

(41)

Using (37), it gives:

<∇ s · u >n =

Ae

+



 (rn+1 − rn )(ur n+1 − ur n ) + ( zn+1 − zn )(u z n+1 − u z n ) +

A ( L e )2

+ +

Ae

−

 n  (r − rn−1 )(ur n − ur n−1 ) + ( zn − zn−1 )(u z n − u z n−1 ) −

u  r

r

A ( L e )2 (42)

n

The same area or length weighted average is used for the surface gradient of tension γ , followed by a projection on the surface tangent plane as in 3D. Finally, the Gaussian curvature K is computed very easily from the radial component of the normal vector since in axisymmetric configuration they are related by:

K =−

nr r

(43)

70

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

Fig. 3. Illustration of the highly oscillating nature of the integrands in Eq. (45) in the case ( zˆ = 1.5, rˆ = 0.2, R c = 2).

2.3. Axisymmetric BEM formulation for vesicles in capillary flows 2.3.1. Axisymmetric Green function for capillary flows In order to study the motion of a vesicle confined in a circular pipe of radius R c , we are interested in the axisymmetric Green’s function M T that takes into account the confinement action of the channel wall. In other words, we require that M T and hence the associated velocity field vanishes on the wall at r = R c . This kind of Green function is simply the sum of two contributions: a free space (FS) contribution M αF βS developed in the c previous section and a correction M α β for the wall action:

M T = M F S + Mc

(44)

Following [37,43], the latter is given by:

∞

c

M =r 0

− F rr cos(ˆzζ )

F rz sin(ˆzζ )

F zr sin(ˆzζ )

F zz cos(ˆzζ )

 dζ

(45)

where F α β are functions of ζ but depends also on the pipe radius R c , on the point force location (r , z) and on r0 , in form of the combinations R c ζ , r0 ζ , and r ζ . Their computation involves the computation of the modified Bessel functions of zeroth and first order I 0 , I 1 , K 0 , K 1 . They vanish as ζ tends to zero, except for the F zz component, which evolves like F zz ∼ 8 ln ζ . However, integrability at the origin can be obtained using a classical regularization, subtracting the logarithmic singularity. Moreover, the exponential decay for large ζ guarantees the integrability of all the F α β at any point within the pipe (see [43,37] for details). 2.3.2. Numerical integration The big difference with the free-space case is that, while in (29) all the terms in the matrix are given coefficients, in (45) a numerical integration on a semi-infinite domain is needed to obtain them. This integration is made difficult by the highly oscillatory nature of the integrands. Furthermore the order of magnitude of the M czz component is really high compared to the other ones, as shown in Fig. 3. So we need a very accurate and adaptive method. We have found that a good way to overcome the difficulty is to carefully divide the domain of integration into several zones with a choice of quadrature adapted to the oscillation frequency of integrands. The Filon formula [42], which is a modified Simpson’s rule, is such a quadrature (see Appendix A), developed especially for oscillatory integrals of the form:

b

b f (ζ ) sin(kζ )dζ or

a

f (ζ ) cos(kζ )dζ a

The Filon formula is advantageous over usual numerical integration formulas for smooth f (ζ ), especially for large k, since the number of points which need to be tabulated depends on the behavior of f (ζ ) rather than on f (ζ ) sin(kζ ). However here we have to deal with the extra problem of a semi-infinite integration domain. To overcome this problem, we split the integration domain into four zones and in each one we apply a different kind of quadratures. For the first interval from ζ ∈ [0, 10−5 ] we use Gaussian quadratures to avoid the singularity at zero. The choice of this interval is made in order to minimize the error on the integration of the fist part of the integral. We have observed numerically that taking a smaller interval does not improve further accuracy. The second and the third ones are computed using the Filon method with an adequate number of integration points n. The last semi-infinite interval is treated using the Gauss–Laguerre quadrature.

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

71

Fig. 4. Error on the absolute value of the four M c components integrated along the wall for a point of force at varying distances from the channel wall.

To test this integration method we compute the flow produced by a ring of point forces whose axis is collinear with the center line of the channel. The associated velocity field should vanish over the surface of the pipe (r = R c ). Fig. 4 shows the absolute value of the four M c components integrated on the wall along the interval z ∈ [−5 : 5] due to a unitary ring of point forces with a varying distance from the channel wall, which is a direct estimation of the error on the velocity at the wall. The error rapidly decreases as the distance of the point forces from the wall increases and it seems to saturate for very small distances, never exceeding 10−3 . 3. Results and discussion 3.1. Validation of the free space version of the axisymmetric model 3.1.1. Comparison to 3D computations In order to validate the free space version of the axisymmetric formulation, we consider first the sedimentation of a vesicle in open space, comparing the final stationary vesicle’s shape to that already obtained with the 3D version of the model [39] and validated by a theoretical study [44]. The vesicle settles under the action of gravity in a surrounding fluid at rest, u ∞ = 0. To fully describe the problem, a new dimensionless parameter must be introduced to characterize the magnitude of the gravitational action: the Bond number (ratio of gravitational to bending energies)

B0 =

ρ g R 40

(46)

κ

where ρ is the difference of density between inner and outer fluids, g the gravity, R 0 the typical size of the vesicle and κ the bending modulus. This validation is illustrated in Fig. 5 for ( = 0.6, B 0 = 47), using 260 linear elements, uniformly distributed on a meridian, for the axisymmetric simulation. The comparison is made to the 3D model with 642 triangular elements uniformly distributed on the whole surface of the vesicle. That corresponds to a mesh accuracy about 5 times lower for a meridian plane compared to the axisymmetric simulation. Computations are performed on an Intel Core i7 x980 and the ratio of computational cost between the axisymmetric computation and the three-dimensional one is 0.1. It demonstrates as expected that the axisymmetric code offers the possibility to refine the mesh (by a factor of about 5 here) with a much smaller computation time. Fig. 6 illustrates numerical convergence as a function of the number of elements in the mesh. We notice that a good description of the final vesicle shape is obtained even with a relative small number of mesh elements. Furthermore the solution quickly converges with mesh refinement. For what concerning computational time we notice that, using as a reference the axisymmetric mesh composed by 260 elements, the ratio of computational cost for less refined meshes is 0.2 for N = 130, 0.06 for N = 64 and 0.02 for N = 32. Comparing the N = 64 axisymmetric case with the previous 3D mesh, both giving almost the same precision of a meridian discretization, we thus obtain a speed up of about 500 times on the computational time. A second test, more related to the capillary flows we are interested in, concerns the motion of a vesicle in an unbounded Poiseuille flow, thus with an imposed background flow given by:



u ∞ = u P 1 − (r / R p )2



(47)

The background flow speed is thus adjusted with the parameter u P that corresponds to the maximal velocity at center line (r = 0). Another relevant parameter is the curvature of the parabolic velocity profile

72

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

Fig. 5. Comparison of stationary settling shape obtained by the 3D and the axisymmetric versions of the vesicle code ( = 0.6, B 0 = 47).

Fig. 6. Comparison of stationary settling shapes obtained by the axisymmetric code for different mesh refinements (N is the number of mesh elements).

C = 2u P / R 2p

(48)

As emphasized in the introduction section, this configuration is often used to tackle the motion of a vesicle in a capillary to mimic a capillary wall at r = R p . Indeed, from Eq. (47) we note that the background flow velocity is required to vanish for r = R p . This is valid only for large enough capillaries, when the confinement effect is negligible. The curvature of the background Poiseuille flow is then the only parameter having a physical meaning. For more confined configurations, the presence of the wall is not a trivial constraint on the shape of the vesicle, which is then physically limited at the capillary radius R c . In this case it is the total velocity field that should vanish at r = R c . Thus, in the bounded case, both the confinement and the Poiseuille flow curvature assume a physical meaning and they should be considered separately. For this example, no density contrast is considered (ρ = Bo = 0) and, in order to characterize the strength of the flow stress exerted on the vesicle, we introduce:

• The bending capillary number Ca = η u P R 20 /κ which replaces the usual capillary number Ca = ηu P /γ for drops, where

(49)

γ is then the interfacial tension [45].

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

73

Fig. 7. Comparison of stationary shape obtained by the 3D code and the axisymmetric one for a vesicle immersed in an unbounded Poiseuille flow ( = 0.6, Caunb = 0.5).

Fig. 8. Time evolution of the vesicle’s height h for a vesicle sedimentation simulation at large Bond numbers ( = 0.93, B 0 = 392).

• The capillary confinement parameter

λ=

R0

(50)

Rc

Actually, for the unbounded Poiseuille flow, only one parameter is needed to characterize the curvature of the background flow velocity. This parameter is expressed in [27] as a combination of Ca and λ (with R c replaced by R p ) in the form of

Caunb =

C η R04

=

Ca λ2

(51) 2 However, the parameter λ will be used in the following to compare vesicles shapes in bounded Poiseuille flow with

κ

those obtained in the unbounded configuration in order to highlight the effect of confinement on the shape of the vesicle. As for the sedimentation example, for this case too, the comparison shows good agreement to the 3D simulation, as can be seen on Fig. 7. 3.1.2. Vesicle sedimentation The axisymmetric formulation presented here has already been used to study the formation of a membrane tube or tether at the rear of an initially prolate settling vesicle [7]. After a first growth, the elongation of the tether reaches a constant value, leading to the stationarity of the system. The stability of the tether is verified over long times as can be seen in Fig. 8. Moreover another common observed phenomena is the appearance of pearls along the tether [6]. And for high sedimentation velocities, vesicles exhibit several tens of microns long pearling tubes of sub-micron radius. Preliminary results indicate that pearls formation is only a transient phenomenon, disappearing when the system reaches the stationary state. The axisymmetric version of the code is able to tackle this kind of transient dynamics, out of reach of the 3D version, as illustrated in Fig. 9. Numerical results are in good agreement with theoretical linear stability analysis [46].

74

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

Fig. 9. Pearls formation for a vesicle sedimentation simulation at large Bond numbers ( = 0.73, B 0 = 9320). Selected dimensionless time are t = 0.39; 5.03; 12.80; 17.45; 19.60; 22.67(×10−4 ).

3.2. Validation of the axisymmetric model for capillary flows 3.2.1. Validation of the capillary axisymmetric formulation In order to validate the extension of the axisymmetric formulation to confined capillary flows, the simple case of a falling rigid sphere in a circular pipe is considered first with increasing confinement parameter λ = R 0 / R c . Since a spherical vesicle with a fixed enclosed volume and local surface incompressibility constraint behaves as a rigid body, this example corresponds actually to the sedimentation of a non-deflated thus spherical vesicle ( = 0). From the Stokes law, the drag force exerted on a spherical object of radius R 0 in an unbounded viscous fluid at rest is given by:

F S = 6πη u S R 0

(52)

where η is the fluid viscosity and u S is the settling velocity. Defining a correction factor K due to the presence of the cylinder wall, the drag force exerted on a sphere in a circular cylinder can then be expressed as:

F = K FS

(53)

Taking u ∞ = 0 and imposing an arbitrary settling velocity u S to all nodes of the membrane mesh, we can easily obtain the corresponding surface force density exerted by the induced flow on the membrane by inversion of the matrix M T of Eq. (44). Integrating then on the membrane surface, we obtain the resulting drag force and thus the correction factor K that can be compared to existing theoretical values [47,48]. Relative error on K estimation versus analytical solution is reported in Fig. 10 as a function of the number N of elements in the particle’s surface mesh. As shown in Fig. 10, third order convergence is obtained. Details on the influence of confinement on the axial component of the surface density force profile are reported in Fig. 11. 3.2.2. Wall meshing alternative The alternative of using the axisymmetric Green function that includes the capillary wall action (with the numerical integration difficulty mentioned in Section 2.3.2) is to consider a capillary of bounded length, using the free space axisymmetric Green function and an explicit representation of the wall by a mesh. The two questions that arises are:

• How long must the capillary be in the upwind and downwind directions? Or in other words, how long must the wall mesh be?

• In addition to the question of the mesh refinement needed for the vesicle’s membrane, what is the mesh refinement needed for the wall mesh? It is well known that in the Stokes regime the velocity perturbation induced by the presence of a particle in a fluid flow decays in r −1 under an external force. This would indicate that rather extended domains have to be used to ensure that the solution will not be influenced, or reasonably influenced, by boundary conditions. However, in [49] it has been shown that the perturbation flow generated by a point-force distribution decays exponentially when confined in a circular pipe. Fig. 12 illustrates the magnitude of the section mean perturbation velocity

 ( z) =

1 Rc

R c

|u − u ∞ | u∞

dr

(54)

0

as a function of axial distance from the sedimenting particle for increasing values of λ. In Fig. 12 the same phenomenon is also visible for a vesicle with λ = 0.5, although the decay profile is not fully upwind and downwind symmetric because of the flow adapted non-symmetric shape that can be taken by the vesicle, contrary to a solid particle (the respective shapes are reported in the middle bottom of Fig. 12).

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

75

Fig. 10. Correction factor K for a sedimenting spherical particle. Left: comparison to theoretical predictions of [47,48] in function of mesh refinement. Right: Numerical convergence as a function of the number N of mesh elements (the black line indicates convergence scaling in N −3 ).

Fig. 11. Influence of confinement (λ = R 0 / R c ) on the axial component f z of the surface density force profile along a meridian (−π /2  θ  π /2 from upwind to downwind) of a sedimenting spherical particle.

Fig. 12. Magnitude of the section mean perturbation velocity as a function of axial distance from the sedimenting particle or vesicle for increasing values of λ = R 0 / R c .

76

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

Fig. 13. Convergence of the particle’s surface density force profile increasing the number N w of wall mesh elements for the case 1/λ = R c / R 0 = 1.1.

From this observation it seems that a shorter capillary length could be used and thus that the meshing of the wall can be a good alternative. However the decay of the perturbative velocity field cannot be taken as the unique criteria to ensure an accurate solution. For studies focusing on vesicle shapes, it is the velocity field at the vesicle’s membrane which is determinant. But very often the preferential indicator used to characterize the problem is related to force quantities [45]. So, in order to objectively compare the two approaches, we proceed as follows:

• Using the first approach (no wall mesh), we find the minimum number N of uniformly distributed mesh elements on the particle’s surface to ensure convergence of the particle’s surface density force profile at a specified tolerated relative error. Note that the force profile like the perturbation flow profile has to be upwind/downwind symmetric for a non-deformable spherical particle. • From these simulations (first approach), we determine the influence length of the particle on the background flow (no background flow for a sedimenting sphere) at a specified tolerated relative error. • For the second approach, we use a capillary length as determined previously, and without changing the number N of elements for the particle’s surface mesh, we find the minimum number N w of mesh elements for the wall to retain the same accuracy on the particle’s surface density force profile. This number can be optimized by using a quadratic refinement rule from tube extremities toward its center. So it is refined where needed, that is in the particle region, since the particle’s referential is used with the particle located at mid length wall. The computations are then made imposing a zero perturbation velocity condition at inflow and outflow sections. In this way the two approaches can be considered as equivalent at the prescribed accuracy, and their computational costs compared. Considering that the influence length is reached when the perturbation flow profile decays to a relative value of less than 10−5 , Fig. 12 shows that the minimal capillary’s length in both upwind and downwind direction is respectively L / R 0 = 9, 7 and 5 for 1/λ = R c / R 0 = 3, 2 and 1.1. Convergence of the particle’s surface density force profile up to a relative error of less than 0.1% is reached for a number of mesh elements on the particle’s surface of 100, 100 and 200 in each case. The same accuracy for the surface density force is obtained with the second approach using at least 690, 1820 and 2921 elements for the wall mesh respectively. For example, convergence is illustrated in Fig. 13 for the case 1/λ = R c / R 0 = 1.1. Computations are performed on an Intel Core i7 x980 and the ratio of computational cost for the two approaches are then respectively 2, 0.4 and 0.3, showing that for large enough capillary radius the wall meshing method is faster, while for strong confinement this method is slower and less accurate. 3.2.3. Convergence study in function of vesicle’s mesh refinement The same procedure can be used to characterize and compare the two approaches for a vesicle in a confined Poiseuille flow. However, the vesicle shape becomes an additional convergence criteria. Fig. 14 compares the stationary shape taken by the vesicle for three uniform mesh refinements (in the case  = 1.59, 1/λ = R c / R 0 = 0.9, Ca = 100). As already observed for the free space configuration, a good description of the shape is obtained even with a relative small number of mesh elements and the solution quickly converges when increasing their number. Indeed a small difference can be observed between the shape obtained with 34 elements and that obtained with 65 elements, and this last shape superposes with that obtained with 120 elements. So a mesh of 65 elements seems sufficient for a good shape description at the vesicle’s length scale for this set of operating parameters. Fig. 15 top shows the corresponding f r and f z vesicle’s surface density force profiles. For acceptable accuracy on these quantities, the 65 elements mesh seems sufficient, in accordance with the shape convergence criteria. Fig. 15 bottom shows the shape and the surface density force profiles for a set of parameters representing a red blood cell flowing through a very narrow pore ( = 4.8, 1/λ = R c / R 0 = 0.55, Ca = 100). In this case, we have to use a minimum number of 130 elements for the vesicle’s mesh.

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

77

Fig. 14. Convergence of the vesicle’s stationary shape increasing the number N of mesh elements on the vesicle’s surface ( = 1.59, 1/λ = R c / R 0 = 0.9, Ca = 100).

Fig. 15. f r (circle labels) and f z (triangular labels) vesicle’s surface density force profiles as a function of vesicle’s surface mesh refinement, with corresponding vesicle shape in the background (dashed line represents the wall position). Top: ( = 1.59, 1/λ = R c / R 0 = 0.9, Ca = 100). Bottom: ( = 4.8, 1/λ = R c / R 0 = 0.55, Ca = 100).

The lubrication film has a mean thickness of about 0.2R 0 for ( = 1.59, 1/λ = R c / R 0 = 0.9, Ca = 100) and 0.03R 0 for ( = 4.8, 1/λ = R c / R 0 = 0.55, Ca = 100). This gives a minimum element size of about one third of the film thickness in the first case and about the double of it in the second case. Actually, in less confined configurations, as illustrated by the first case, it is the vesicle’s shape itself that imposes the minimum element size, especially in order to be able to describe the high curvature region in the corner of the vesicle’s back. In stronger confinements, as illustrated by the second case, the film thickness criteria becomes more constraining. Concerning the second approach (computation with a meshed wall), the perturbation flow decay criteria leads us for both cases to limit the capillary’s length to 5R 0 in both upwind and downwind directions. It must be underlined that, contrary to the spherical particle case, the quadratic refinement strategy to minimize the number of elements for the wall isn’t adequate for an elongated vesicle. So, for the wall mesh, we use instead a uniform distribution of elements, as for the vesicle’s surface. The minimum elements number determined by the accuracy criteria on the vesicle’s surface density force is then 103 for the first case and 1451 for the second. Comparing computational time, it appears that the first approach becomes more advantageous as the confinement increases. For example, for the second case we obtain a computational time ratio of 0.8. 3.3. Applications 3.3.1. Influence of confinement for Poiseuille flows A first non-trivial issue is the vesicle stationary shape dependence on initial conditions. In [24] it was already predicted the existence of different shapes, a bulletlike and a parachutelike shape, for a small curvature of the flow field. But, for large enough curvature of the flow, the two solutions becomes indistinguishable. In [50] the authors show that, for weak confinement and very low capillary number (Ca = 1), stationary shapes obtained starting from oblate and prolate configurations are not the same. They obtain a parachute-like shape with a concave rear part in the first case and a bullet-like shape with a convex rear part in the second one. In order to confirm the predictions of [24] for larger values of capillary numbers, we

78

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

Fig. 16. Time evolution of the vesicle shape for  = 0.2 and λ = 0.5 under the action of a Poiseuille flow with Ca = 100 in the case of prolate (top) and oblate (bottom) initial shape. Selected times are t = 0; 0.01; 0.03; 0.13.

Fig. 17. Comparison of stationary shapes for bounded (red line) and unbounded (blue line) configurations for several values of confinement parameter λ in the case of a vesicle with  = 1.59 and Ca = 100. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

repeat the experience of [50] for Ca = 100. Fig. 16 shows the shape evolution over time starting from a prolate (top) or an oblate (bottom) for  = 0.2 and 1/λ = R c / R 0 = 2). It clearly shows that the two dynamics tend toward the same stationary shape, and this is verified as well for the set ( = 0.9, λ = 0.3) used in [50]. We can therefore assume that for large enough capillary number the initial shape, even if it strongly modifies the dynamics, plays no role on the stationary shape. The second and main issue is the effect of confinement on the vesicle’s stationary shape, and this question can be answered by comparing results coming from simulations with bounded and unbounded Poiseuille flows. First we want to validate the key assumption of [23] that for λ  0.1 the presence of lateral walls can be disregarded. Effectively, as shown in Fig. 17(a), we obtain that the stationary shape at 1/λ = R c / R 0 = 10 is the same for both bounded and unbounded configuration. But decreasing the Poiseuille radius (and thus the wall distance in the bounded case) the stationary shape starts soon to differ for the two configurations (Fig. 17(b) and Fig. 17(c)). This result shows that for higher values of confinement (R c  5R 0 ) the influence of the tube wall on the vesicle’s shape cannot be investigated in the frame of the unbounded Poiseuille flow. Experiments with giant lipid vesicles flowing through cylindrical capillaries with diameters close to the vesicle size and a constant imposed flow rate [26] show that a significant deformation of the membrane occurs and increases when the velocity, confinement or deflation of the vesicle is increased. As they flow, isolated vesicles experience a lift force that pushes them away from the wall [9]. But for a very small density difference between the interior and exterior fluids, vesicles tend to align along the axis of the micropipette. The steady-state flow and vesicle shapes are then axisymmetric. Fig. 18

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

79

Fig. 18. Comparison of simulated shapes to the observed ones in the experimental study of [26]. Vesicles are flowing from left to right. The corresponding sets of operating parameters are from left to right ( = 0.45, 1/λ = 1.25, Ca = 10), ( = 1.7, 1/λ = 1.5, Ca = 30) and ( = 2, 1/λ = 2.5, Ca = 60).

Fig. 19. Stationary shapes of a vesicle sedimenting in a capillary for different values of the confinement parameter λ ( = 0.93, B 0 = 216).

shows the comparison between experimental and simulated shapes with our BEM axisymmetric model for capillary (first approach) for several values of deflation, confinement and flow velocity. Predicted stationary shapes match very well with the images resulting from experiments in a quite large range of parameters (, λ, Ca). 3.3.2. Influence of confinement for sedimentation Having presented sedimentation simulations to validate the free space version of the axisymmetric model, a natural question concerns the influence that confinement would have on sedimentation. To briefly answer this question, we present some results on the tether formation for vesicles sedimenting in a capillary. Starting from one of the cases of tethered vesicles settling in free space presented in [7], we study the evolution of the tether length when increasing the confinement from 1/λ = R c / R 0 = 10 up to 1/λ = R c / R 0 = 1.15. As expected, for 1/λ  5 no significant variations from the shape obtained in the free space configuration can be observed, but for lower confinements we observe a drastic reduction of the tether as shown in Fig. 19. In particular at 1/λ = 2 we tend again to the typical pear-like shape of smaller Bond numbers already studied in [39]. Finally at 1/λ = 1.15 the tube disappears and the vesicle takes a bullet like shape more typical of capillary flows. 4. Conclusions and perspectives When more than one characteristic length scale are involved, computations in 3D can be too expensive, and a 2D description is ever questionable. However, often there exists a large domain of existence of axisymmetric regimes for which an axisymmetric model would be fully appropriate. We propose such an alternative on the bases of our earlier 3D model [39]. With the axisymmetric assumption, this model presents the advantage to take a very simple form with all contributions to the membrane forces writing as explicit algebraic expressions in function of nodal quantities. Main efforts concern the implementation of an appropriate Green function. For free space flows, we have demonstrated the accuracy and efficiency of the axisymmetric formulation on the case of sedimentation and non-bounded Poiseuille flows, by comparing to 3D computations. We have shown the usefulness of the approach through results on tether and pearling phenomena that cannot be tackled with 3D computations. But the main motivation of the work concerns capillary flows, especially capillary Poiseuille flows. We have shown that for confinement parameters λ greater than 1/5, considering unbounded Poiseuille flows is no longer possible. Actually, the confinement action of the capillary wall must be considered either using an adapted Green function or by representing explicitly the wall by a mesh. For the first approach, we have proposed a way to overcome the difficulty of numerical integration of the highly oscillating additional terms to the free space Green function, and its accuracy has been demonstrated on the case of a sedimenting spherical particle in a capillary, for which theoretical results exist. For the second approach, we have used the computations with the first one to determine the minimal capillary length in order to impose non-perturbing inflow and outflow boundary conditions. The question of mesh refinement has been answered in the same way, determining the optimal mesh for

80

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

the particle with the first approach and that for the wall with the second. We have then extended this characterizing study to the case of vesicles in capillary bounded Poiseuille flows, determining optimal meshing rules in relation to the vesicle’s shape and/or to the lubrication film thickness. From the point of view of computational cost, the two approaches appear comparable with an advantage to the second for lower confinement configurations and to the first when strongly confined. This axisymmetric approaches seem promising to study the influence of confinement on vesicle’s shape in a very large range of confinement. The comparison to experimentally observed shapes is highly satisfactory, and we have already brought an answer on the open question concerning the importance of the vesicle’s initial shape. Moreover the results presented for a set of operating parameters typical of red blood cells flowing through very narrow pores show that we are able to consider highly confined situations with lubrication film thickness as low as five percent of the capillary’s radius. In the future it should be useful to have a better characterization of vesicles properties in strong confinement, such as the bullet-parachute shape transition for vesicles in little capillaries for large ranges of characteristic parameters λ,  and Ca. But other type of capillary flows are of interest too, as showed by our preliminary results on vesicle sedimenting in a capillary. Acknowledgements This work has benefited from financial support from the ANR CAPSHYDR (Grant No. 11-BS09-013-02), from Labex MEC (Grant No. ANR-11-LABX-0092), from A*MIDEX (Grant No. ANR-11-IDEX-0001-02), from the ANR Polytransflow (Grant No. 13-BS09-0015-01), and from CNES. This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01). Appendix A. Filon formula For a given function f (x) within a closed interval [a; b], the interval must be divided into 2p sections where each one has width h. Within each section, f (x) is replaced by a polynomial approximation, in particular a second-degree polynomial. The integration formula is applied to each panel. Finally the sum of the contributions from each section gives the desired quadrature formula. These formulas are [51]

S =h









C =h

α ( f 0 cos(kx0 ) − f 2p cos(kx2p )) + β S 2p + γ S 2p−1 + E s

(A.1)

α ( f 2p cos(kx2p ) − f 0 cos(kx0 )) + β C 2p + γ C 2p−1 + E c

(A.2)

where

S 2p =

p 

f 2i sin(kx2i ) −

i =0

S 2p −1 =

p 

1 2



f 0 sin(kx0 ) + f 2p sin(kx2p )

f 2i −1 sin(kx2i −1 )

(A.3)

(A.4)

i =1

C 2p =

p 

f 2i cos(kx2i ) −

i =0

C 2p −1 =

p 

1 2



f 0 cos(kx0 ) + f 2p cos(kx2p )

f 2i −1 cos(kx2i −1 )

(A.5)

(A.6)

i =1

α = 1/θ + sin(2θ)/2θ 2 − 2 sin2 (θ)/θ 3

  β = 2 (1 + cos2 (θ))/θ 2 − sin(2θ)/θ 3   γ = 4 sin(θ)/θ 3 − cos(θ)/θ 2 θ = kh

(A.7) (A.8) (A.9) (A.10)

where we call f i = f (xi ), xi +1 − xi = h, x0 = a, x2p = b and E s and E c are the errors associated with using the first term on the right of Eqs. (A.1) and (A.2) as an approximation for S and C . For small θ it is necessary to replace the expressions for α , β and γ by power series in θ to avoid the loss of significant figures due to cancellation in these expressions. So we have [42]:

α = 2θ 3 /45 − 2θ 5 /315 + 2θ 7 /4725 2

4

6

β = 2/3 + 2θ /15 − 4θ /105 + 2θ /567 2

4

6

γ = 4/3 − 2θ /15 + θ /210 − θ /11 340

(A.11) (A.12) (A.13)

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

81

References [1] X. Li, P.M. Vlahovska, G.E. Karniadakis, Continuum- and particle-based modeling of shapes and dynamics of red blood cells in health and disease, Soft Matter 9 (2013) 28–37, http://dx.doi.org/10.1039/C2SM26891D. [2] D. Abreu, M. Levant, V. Steinberg, U. Seifert, Fluid vesicles in flow, Adv. Colloid Interface Sci. 208 (2014) 129–141, http://dx.doi.org/10.1016/j.cis.2014.02, http://www.sciencedirect.com/science/article/pii/S0001868614000384. [3] P.M. Vlahovska, D. Barthes-Biesel, C. Misbah, Flow dynamics of red blood cells and their biomimetic counterparts, C. R. Phys. 14 (6) (2013) 451–458, http://dx.doi.org/10.1016/j.crhy.2013.05.001, http://www.sciencedirect.com/science/article/pii/S1631070513000765. [4] V. Kantsler, E. Segre, V. Steinberg, Critical dynamics of vesicle stretching transition in elongational flow, Phys. Rev. Lett. 101 (2008) 048101, http://dx.doi.org/10.1103/PhysRevLett.101.048101. [5] H. Zhao, E.S.G. Shaqfeh, The shape stability of a lipid vesicle in a uniaxial extensional flow, J. Fluid Mech. 719 (2013) 345–361. [6] Z.-H. Huang, M. Abkarian, A. Viallat, Sedimentation of vesicles: from pear-like shapes to microtether extrusion, New J. Phys. 13 (3) (2011) 035026, http://stacks.iop.org/1367-2630/13/i=3/a=035026. [7] G. Boedec, M. Jaeger, M. Leonetti, Sedimentation-induced tether on a settling vesicle, Phys. Rev. E 88 (2013) 010702, http://dx.doi.org/10.1103/PhysRevE. 88.010702. [8] S. Chatkaew, M. Georgelin, M. Jaeger, M. Leonetti, Dynamics of vesicle unbinding under axisymmetric flow, Phys. Rev. Lett. 103 (24) (2009) 248103, http://dx.doi.org/10.1103/PhysRevLett.103.248103. [9] M. Abkarian, C. Lartigue, A. Viallat, Tank treading and unbinding of deformable vesicles in shear flow: determination of the lift force, Phys. Rev. Lett. 88 (6) (2002) 068103, http://dx.doi.org/10.1103/PhysRevLett.88.068103. [10] I. Rey Suárez, C. Leidy, G. Téllez, G. Gay, A. Gonzalez-Mancera, Slow sedimentation and deformability of charged lipid vesicles, PLoS ONE 8 (7) (2013) e68309, http://dx.doi.org/10.1371/journal.pone.0068309. [11] C. Quéguiner, D. Barthès-Biesel, Axisymmetric motion of capsules through cylindrical channels, J. Fluid Mech. 348 (1997) 349–376. [12] T.W. Secomb, R. Skalak, A two-dimensional model for capillary-flow of an asymmetric cell, Microvasc. Res. 24 (1982) 194–203. [13] T.W. Secomb, Flow-dependent rheological properties of blood in capillaries, Microvasc. Res. 34 (1987) 46–58. [14] T.W. Secomb, R. Skalak, N. Ozkaya, J.F. Gross, Flow of axisymmetric red blood cells in narrow capillaries, J. Fluid Mech. 163 (1986) 405–423, http://dx.doi.org/10.1017/S0022112086002355. [15] H. Noguchi, G. Gompper, Shape transitions of fluid vesicles and red blood cells in capillary flows, Proc. Natl. Acad. Sci. USA 102 (40) (2005) 14159–14164, http://dx.doi.org/10.1073/pnas.0504243102, http://www.pnas.org/content/102/40/14159.full.pdf+html. [16] G.G.L. McWhirter, H. Noguchi, Flow-induced clustering and alignment of vesicles and red bloods cells in microcapillaries, Proc. Natl. Acad. Sci. USA 106 (2009) 6039. [17] J.B. Freund, M.M. Orescanin, Cellular flow in a small blood vessel, J. Fluid Mech. 671 (2011) 466–490, http://dx.doi.org/10.1017/S0022112010005835. [18] D.A. Fedosov, M. Peltomaki, G. Gompper, Deformation and dynamics of red blood cells in flow through cylindrical microchannels, Soft Matter 10 (2014) 4258–4267. [19] G. Coupier, B. Kaoui, T. Podgorski, C. Misbah, Noninertial lateral migration of vesicles in bounded Poiseuille flow, Phys. Fluids 20 (11) (2008) 111702, http://dx.doi.org/10.1063/1.3023159. [20] A. Farutin, C. Misbah, Analytical and numerical study of three main migration laws for vesicles under flow, Phys. Rev. Lett. 110 (2013) 108104, http://dx.doi.org/10.1103/PhysRevLett.110.108104, http://link.aps.org/doi/10.1103/PhysRevLett.110.108104. [21] A. Farutin, C. Misbah, Symmetry breaking and cross-streamline migration of three-dimensional vesicles in an axial Poiseuille flow, Phys. Rev. E 89 (2014) 042709, http://dx.doi.org/10.1103/PhysRevE.89.042709. [22] B. Kaoui, B. Ristow, I. Cantat, C. Misbah, W. Zimmermann, Lateral migration of a two-dimensional vesicle in unbounded Poiseuille flow, Phys. Rev. E 77 (2008) 021903. [23] B. Kaoui, G. Biros, C. Misbah, Why do red blood cells have asymmetric shapes even in a symmetric flow?, Phys. Rev. Lett. 103 (2009) 188101. [24] G. Danker, P.M. Vlahovska, C. Misbah, Vesicles in Poiseuille flow, Phys. Rev. Lett. 102 (2009) 148102, http://dx.doi.org/10.1103/PhysRevLett.102.148102. [25] P. Vlahovska, T. Podgorski, C. Misbah, Vesicles and red blood cells in flow: from individual dynamics to rheology, in: Complex and Biofluids, C. R. Phys. 10 (8) (2009) 775–789, http://dx.doi.org/10.1016/j.crhy.2009.10.001. [26] V. Vitkova, M. Mader, T. Podgorski, Deformation of vesicles flowing through capillaries, Europhys. Lett. 68 (3) (2004) 398, http://stacks.iop.org/ 0295-5075/68/i=3/a=398. [27] G. Coupier, A. Farutin, C. Minetti, T. Podgorski, C. Misbah, Shape diagram of vesicles in Poiseuille flow, Phys. Rev. Lett. 108 (2012) 178106, http://dx.doi.org/10.1103/PhysRevLett.108.178106. [28] H. Noguchi, G. Gompper, L. Schmidt, A. Wixforth, T. Franke, Dynamics of fluid vesicles in flow through structured microchannels, Europhys. Lett. 89 (2010) 28002. [29] B. Kaoui, N. Tahiri, T. Biben, H. Ez-Zahraouy, A. Benyoussef, G. Biros, C. Misbah, Complexity of vesicle microcirculation, Phys. Rev. E 84 (2011) 041906, http://dx.doi.org/10.1103/PhysRevE.84.041906. [30] A. Farutin, C. Misbah, Symmetry breaking of vesicle shapes in Poiseuille flow, Phys. Rev. E 84 (1) (2011) 011902, http://dx.doi.org/10.1103/ PhysRevE.84.011902. [31] C. Pozrikidis, The buoyancy-driven motion of a train of viscous drops within a cylindrical tube, J. Fluid Mech. 237 (1992) 627–648. [32] X.Z. Li, D. Barthes-Biesel, A. Helmy, Large deformations and burst of a capsule freely suspended in an elongational flow, J. Fluid Mech. 187 (1988) 179–196. [33] A. Diaz, D. Barthes-Biesel, Entrance of a bioartificial capsule in a pore, Comput. Model. Eng. Sci. 3 (2002) 321. [34] Y. Lefebvre, D. Barthes-Biesel, Motion of capsules in a cylindrical tube: effect of membrane pre-stress, J. Fluid Mech. 589 (2007) 157. [35] S. Veerapaneni, D. Gueyffier, G. Biros, D. Zorin, A numerical method for simulating the dynamics of 3d axisymmetric vesicles suspended in viscous flows, J. Comput. Phys. 228 (19) (2009) 7233–7249. [36] E.P. Ascoli, D.S. Dandy, L.G. Leal, Buoyancy-driven motion of a deformable drop toward a planar wall at low Reynolds number, J. Fluid Mech. 213 (1989) 287–311. [37] C. Pozrikidis, Boundary Integral and Singularity Methods for Linearized Viscous Flow, Cambridge University Press, 1992. [38] C. Pozrikidis, Axisymmetric motion of a file of red blood cells through capillaries, Phys. Fluids 17 (3) (2005) 031503. [39] G. Boedec, M. Leonetti, M. Jaeger, 3d vesicle dynamics simulations with a linearly triangulated surface, J. Comput. Phys. 230 (2011) 1020–1034. [40] L. Scheffer, A. Bitler, E. Ben-Jacob, R. Korenstein, Atomic force pulling: probing the local elasticity of the cell membrane, Eur. Biophys. J. 30 (2) (2001) 83–90. [41] M. Meyer, M. Desbrun, P. Schroeder, A. Barr, Discrete differential geometry operators for triangulated 2-manifolds, in: VisMath ’02 Proceedings, 2002, pp. 237–242. [42] M. Abramowitz, I.A. Stegun (Eds.), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover Publications, New York, 1972. [43] H. Tozeren, Boundary integral equation method for some stokes problems, Int. J. Numer. Methods Fluids 4 (2) (1984) 159–170.

82

R. Trozzo et al. / Journal of Computational Physics 289 (2015) 62–82

[44] G. Boedec, M. Jaeger, M. Leonetti, Settling of a vesicle in the limit of quasispherical shapes, J. Fluid Mech. 690 (2012) 227–261, http://dx.doi.org/ 10.1017/jfm.2011.427. [45] E. Lac, J.D. Sherwood, Motion of a drop along the centreline of a capillary in a pressure-driven flow, J. Fluid Mech. 640 (2009) 27–54. [46] G. Boedec, M. Jaeger, M. Leonetti, Pearling instability of a cylindrical vesicle, J. Fluid Mech. 743 (2014) 262–279, http://dx.doi.org/10.1017/jfm.2014.34. [47] J.J.L. Higdon, G.P. Muldowney, Resistance functions for spherical particles, droplets and bubbles in cylindrical tubes, J. Fluid Mech. 298 (1995) 193–210. [48] C. Pozrikidis, J.M. Davis, Resistance and pressure coefficients for a periodic array of spherical, spheroidal, and cylindrical particles inside a circular tube, IMA J. Appl. Math. (2015), in press. [49] N. Liron, R. Shahar, Stokes flow due to a stokeslet in a pipe, J. Fluid Mech. 86 (1978) 727–744, http://dx.doi.org/10.1017/S0022112078001366. [50] W.-F. Hu, Y. Kim, M.-C. Lai, An immersed boundary method for simulating the dynamics of three-dimensional axisymmetric vesicles in Navier–Stokes flows, J. Comput. Phys. A 257 (2014) 670–686, http://dx.doi.org/10.1016/j.jcp.2013.10.018, http://www.sciencedirect.com/science/article/pii/ S0021999113006943. [51] L.D. Fosdick, A special case of the Filon quadrature formula, Math. Comput. 22 (101) (1968) 77–81, http://www.jstor.org/stable/2004764.