Simulating the dynamics of inextensible vesicles by the penalty immersed boundary method

(1)

Simulating the dynamics of inextensible vesicles by the penalty

immersed boundary method

Yongsam Kim

a,1

, Ming-Chih Lai

b,*

a

Department of Mathematics, Chung-Ang University, Dongjak-Gu, Heukseok-Dong 221, Seoul 156-756, South Korea b

Center of Mathematical Modeling and Scientiﬁc Computing, Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan

a r t i c l e

i n f o

Article history:

Received 21 September 2009

Received in revised form 24 February 2010 Accepted 12 March 2010

Available online 17 March 2010 Keywords:

Penalty immersed boundary method Inextensible vesicles

Binary-component vesicles

a b s t r a c t

In this paper, we develop an immersed boundary (IB) method to simulate the dynamics of inextensible vesicles interacting with an incompressible fluid. In order to take into account the inextensibility constraint of the vesicle, the penalty immersed boundary (pIB) method is used to virtually decouple the fluid and vesicle dynamics. As numerical tests of our cur-rent pIB method, the dynamics of single and multiple inextensible vesicles under shear flows have been extensively explored, and compared with the previous literature. The method is also validated by a series of convergence study, which confirms its consistent first-order accuracy on the velocity field, the vesicle configuration, the vesicle area and the perimeter errors. In addition, the method is also applied to study a binary-component vesicle problem.

1. Introduction

We introduce an immersed boundary (IB) method[25]to simulate the dynamics of inextensible vesicles in an incom-pressible viscous fluid. The vesicle dynamics is important for numerous biophysical systems that involve deforming fibers interacting with a surrounding fluid. The inextensible vesicles have attracted a great deal of attention recently as they are present in many biological systems[18]and can be considered as models of red blood cells [24,29] and drug-carrying capsules[34]. The dynamics of moving vesicles, which are determined by their boundary elasticity, inextensibility, and hydrodynamical force, have been extensively explored both experimentally[7,12](and the references therein) and compu-tationally[18,38,36].

Vesicle simulations have been conducted mostly by molecular dynamics models[23]or by continuum mechanics models of the fluid and the vesicle membrane using the boundary integral method[36]. Here, we use a Cartesian-grid based IB meth-od[25]. The IB method is generally a useful computational tool for the problems in which a thin elastic boundary is im-mersed in and interacting with a surrounding fluid[4,3,26,11,13,17]. In[14], the authors have introduced a new version of the IB method, the penalty immersed boundary (pIB) method, which can handle efficiently the case of massive elastic boundary, and have shown that the new method can be applied to many problems in which the mass of the immersed boundary plays a crucial dynamical role, see also[2,9,15,16]for more examples of the application of the pIB method. This method has the virtue of simplicity; one can easily implement it in the context of an existing IB method code for the massless boundary case.

* Corresponding author. Tel.: +886 3 5731361; fax: +886 3 5724679.

E-mail addresses:kimy@cau.ac.kr(Y. Kim),mclai@math.nctu.edu.tw(M.-C. Lai). 1

Tel.: +82 2 820 5219; fax: +82 2 823 5214.

Contents lists available atScienceDirect

Journal of Computational Physics

(2)

The basic idea of the pIB method[14]is to use two Lagrangian immersed boundaries to represent the real immersed boundary for different purposes. The first immersed boundary is massless and interacts with the fluid directly as in the tra-ditional IB computation, i.e., it moves at the local fluid velocity and applies the force locally to the fluid. The other immersed boundary is massive which carries the fiber mass and is connected to the previous massless boundary by a system of stiff springs. However, this massive boundary does not interact with the fluid directly and it moves according to Newton’s law of motion with the forces generated by the system of stiff springs that connects the two boundary representations. Roughly speaking, as the spring stiffness tends to infinity, the massive boundary will track the motion of the massless boundary and endow it with the mass. This behavior has been verified quantitatively in[14].

In order to apply the penalty IB idea to the inextensible vesicle problems, as in[14], we represent the vesicle boundary by two Lagrangian immersed boundaries. One vesicle boundary interacts with the ﬂuid dynamics directly and the other vesicle boundary follows the equations of the vesicle dynamics, including the inextensibility constraint, without a direct interaction with the ﬂuid dynamics. The two boundaries are connected by the penalty forces which act on both boundaries by Newton’s third law. This idea of separation of dynamical systems using a penalty idea has also been used in[10].

To show that the present developed pIB method is a robust and efficient numerical tool to handle the inextensible bound-ary in a fluid, we simulate several benchmark two-dimensional problems with single and multiple vesicles under a shear flow. From these examples, we show that the simulation results obtained by the pIB method are comparable to those ob-tained from previous literature[7,12,38,36,23]. In order to provide further verification that the pIB method correctly solves the inextensible vesicle dynamics, we perform a convergence study, which confirms the consistent first-order accuracy on the velocity field, the vesicle configuration, the vesicle area error, and the vesicle perimeter error. Finally, we extend the method to solve a binary-component vesicle model which is obtained by simplifying the model in[32].

2. Equations of motion

We begin by stating the mathematical formulation of the equations of motion for a system comprised of a two-dimen-sional viscous incompressible ﬂuid in a domainXcontaining an immersed, inextensible, massless vesicle boundaryC. We assume that the ﬂuid is the same inside and outside of the vesicle boundary, and the equations of motion in immersed boundary formulation are

q

@u @t þ u

r

u ¼

r

p þ

l

r

2_{u þ f;} ð1Þ

r

u ¼ 0; ð2Þ fðx; tÞ ¼ Z C Fðs; tÞdðx Xðs; tÞÞds; ð3Þ @X @tðs; tÞ ¼ Uðs; tÞ ¼ Z X uðx; tÞdðx Xðs; tÞÞdx; ð4Þ Fðs; tÞ ¼ @ @sð

r

r

ðs; tÞ is an unknown function which plays the role of a Lagrange multiplier to en-force the inextensibility constraint of the vesicle, Eq.(7). This inextensibility constraint is equivalent to the fact that the local stretching factor of the vessel satisﬁes d

dtj @X

@sj ¼ 0. To see how the boundary force can be expressed as Eq.(5), we explain as

(3)

The vesicle membrane energy can be modelled by the Helfrich type energy ðEbÞ[8]to resist the bending of the membrane.

To enforce the inextensibility constraint, this energy is replaced by a modiﬁed form using the method of Lagrangian multi-plier. More precisely, the total energy can be written as

E ¼ Eeþ Eb¼ Z C

r

ds þ Z C cb 2 @2X @s2 2 ds; ð8Þ

where

r

ðs; tÞ is an unknown surface tension and cbis the constant bending rigidity. By taking the variational derivative of the

above energy form, the boundary force Eq.(5)is obtained by

F ¼ Feþ Fb¼ dEe dXþ dEb dX : ð9Þ

In[38], Zhou and Pozrikidis consider the deformation of liquid capsules with incompressible (same as inextensible) inter-faces in simple shear flow. They assume only the membrane tension force (first term in the right-hand side of Eq.(5)) to enforce the inextensibility. By using the boundary element formulation for Stokes flow, the authors discretize the whole equations by solving nonlinear system of vesicle boundary positions and the associated tension. Recently, Veerapaneni et al.[36]used a similar boundary integral method to simulate the dynamics of 2D[36]and 3D axis-symmetric[37] inexten-sible vesicles suspended in Stokes fluid. Unlike[38], Veerapaneni et al.[36,37]take the bending force into account. In this paper, we use the immersed boundary method to study the dynamics of inextensible vesicles in Navier–Stokes flow. Further-more, in a later section, we consider the case of binary-component vesicle in which the bending rigidity cband spontaneous

curvature are all dependent on the local heterogeneous concentration of those two components corresponding to liquid-or-dered and liquid-disorliquid-or-dered surface phases[22]. We note that a thermodynamically consistent model has been proposed to simulate the dynamics of two-dimensional multi-component vesicles in Stokes ﬂow[32].

2.1. Penalty IB method

Since the system of Eqs.(1)–(7)is complete, one can directly solve the whole system using some numerical methods. To enforce the inextensibility constraint, however, the methods should handle this constraint realized by Eqs.(5)–(7)in an im-plicit way, which requires an iterative method to solve the whole system. To avoid this complication and make the scheme efficient, we modify the formulation of Eqs.(1)–(7)so that the penalty IB (pIB) method[14]can be applied. In the pIB meth-od, we use a dual representation of the immersed boundary. One of its representatives, denoted by Xðs; tÞ, interacts with the fluid and moves at the local fluid velocity, just as in the case of an immersed boundary in traditional IB computation. The other representative, denoted by Yðs; tÞ, is elastic and linked to Xðs; tÞ by a system of stiff springs. The boundary points of Yðs; tÞ are not coupled directly to the fluid, and they move as if in a vacuum, seeFig. 1.

In mathematical terms, we replace Eqs.(5)–(7)by the following:

F ¼ KðYðs; tÞ Xðs; tÞÞ þ RðVðs; tÞ Uðs; tÞÞ; ð10Þ KðYðs; tÞ Xðs; tÞÞ þ RðVðs; tÞ Uðs; tÞÞ ¼@ @sð

r

ðs; tÞ

s

ðs; tÞÞ cb @4Yðs; tÞ @s4 ; ð11Þ

s

ðs; tÞ ¼@Y @s _@Y @s ; ð12Þ @V @s @Y @s¼ 0; ð13Þ Spring boundary _Y boundary _X

(4)

where U ¼ @X=@t (defined in Eq.(4)) and V ¼ @Y=@t are the velocities of the two representations of the immersed boundary. Eq.(10)defines the force density F which is transmitted by the boundary Xðs; tÞ to the fluid. It includes only the force density generated by the stiff springs that link the two representations of the immersed boundary and K and R are the penalty parameters. The larger they are, the greater the energy penalty that must be paid to separate the positions and velocities of the two boundary representations by any given amount. Eqs. (11)–(13)are the equations of motion and velocity of the boundary Yðs; tÞ. Note that the penalty force is instantaneously balanced by the elastic force of the boundary Yðs; tÞ and this elastic force does not act directly on the fluid, i.e., the boundary Yðs; tÞ does not interact with the fluid directly.

Generally speaking, as the penalty parameters K and R go to infinity, the boundaries Xðs; tÞ and Yðs; tÞ and velocities Uðs; tÞ and Vðs; tÞ will tend to coincide, and F in Eq.(10)approaches the force F in Eq.(5). In practice, we choose K and R sufficiently large to keep the positions and velocities of the two representations of the boundary close to each other. Notice that, the present pIB idea shares the same spirit with the ‘virtual boundary method’[5,30]for simulating the fluid interacting with rigid boundaries. Both methods introduce the penalty parameters K and R to penalize the difference in the desired and numerical velocities and positions of the immersed boundary. The difference is that, the virtual boundary method is devel-oped to impose the no-slip boundary condition while the present pIB method is to impose the inextensibility condition on the immersed boundaries. Besides, an extra unknown surface tension must be solved in the present approach.

The pIB method is an extension of the original IB method to handle the case where the mass of the immersed boundary plays an important dynamical role; for example, flapping filament in a flowing soap film, falling parachute in a steady flow, and compliant vessel interacting with blood[14,17,16]. The present pIB method is different from the one in[14]in two as-pects; namely, (1) no mass is considered here, (2) an additional damping-like force RðVðs; tÞ Uðs; tÞÞ is added to the penalty force. Nevertheless, the massive boundary case can be handled simply by adding the D’Alembert force MðsÞ@2Yðs;tÞ

@t2 to the left

side of Eq.(11)where MðsÞ is the mass density of the boundary. Then Eq.(11)simply represents Newton’s law of motion. The addition of RðVðs; tÞ Uðs; tÞÞ to the penalty force makes it possible to apply the inextensibility constraint (Eq.(13)) to the force balance equation (Eq.(11)) by means of the velocity Vðs; tÞ. A similar idea was also used in[10]. The advantage of the pIB method is that we only need to solve Eqs.(10)–(13)to enforce the inextensibility constraint, which are decoupled from the ﬂuid equations.

3. Numerical implementation of the pIB method

For numerical implementation, we use a ﬁrst-order IB method, generalized to take into account the equations for the boundary Yðs; tÞ that is linked to the boundary Xðs; tÞ by stiff springs. We use a superscript to denote the time level; thus, XnðsÞ is a shorthand for Xðs; nDtÞ, whereDt is the size of the time step, and similarly for all other variables. Our goal is to compute updated unþ1_; _Vnþ1_; _Xnþ1_{, and Y}nþ1_{from given data u}n_; _Vn_; _Xn_{and Y}n_.

The step-by-step procedure of the numerical implementation proceeds as follows:

tVnþ1ðsÞ: ð23Þ

First, we calculate the force density Fn

ðsÞ deﬁned on Lagrangian grid points and spread it into the Eulerian grid force fnðxÞ in the Navier–Stokes equations. These are done by Eqs.(14) and (15)which are the discretizations of Eqs.(10) and (3), respec-tively. Here and throughout the paperPxdenotes the sum over the rectangular lattice in physical space on which the ﬂuid

variables are deﬁned. Similarly,P_swill denote the sum over the boundary nodes on which the two boundary positions X; Y, and the force density F are deﬁned.

(5)

dhðxÞ ¼ h 2 / x1 h / x2 h ; ð24Þ

where x ¼ ðx1;x2Þ and h ¼Dx ¼Dy is the mesh width of the uniform square lattice. The function / is determined by the

following conditions: 1. / is a continuous function. 2. /ðrÞ ¼ 0 for jrj P 2. 3. For all r; Pj¼even/ðr jÞ ¼

P

j¼odd/ðr jÞ ¼12.

4. For all r; P_jðr jÞ/ðr jÞ ¼ 0.

5. For all r; P_jð/ðr jÞÞ2¼ C, where C is a constant, independent of r.

The motivation of these postulates is discussed in[26,27]. It follows that C ¼3

8and / is uniquely determined as

/ðrÞ ¼ 32jrjþpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1þ4jrj4r2 8 if jrj 6 1 52jrjpffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi7þ12jrj4r2 8 if 1 6 jrj 6 2 0 if 2 6 jrj: 8 > > > < > > > : ð25Þ

This is an even, bell-shaped function, which is not only continuous but also has a continuous ﬁrst derivative.

With the force density fnðxÞ in hand, we can turn to solving the Navier–Stokes equations. This is done by Eqs.(16) and (17)in which D0

i is the standard central difference operator in the spatial direction denoted by i, where i ¼ 1; 2, and L is

the standard 5-point discrete Laplacian. Note that the skew-symmetric difference is used for the convection term[19,25]. The vector operator D ¼ ðD1;D2Þ that is used for the discrete gradient and divergence can be deﬁned as follows[28]:

ðxÞ ¼ Rh=2

h=2 dhðx þ XÞdX. This operator is designed for ‘‘improved volume conservation”

and is constructed according to a recipe introduced in[28]which ensures that the value of D u at any particular grid point is exactly equal to the average over an h h square centered on that grid point of the continuous divergence of the interpo-lated velocity ﬁeld in which the immersed boundary points move. This deﬁnition ‘‘tunes” the operator D to the function dh

that deﬁnes the interpolation scheme of the immersed boundary method.

Eqs.(16) and (17)form a linear system with constant coefﬁcients in the unknowns unþ1_;_pnþ1_{. We solve this linear system}

by taking its discrete Fourier transform (implemented by the FFT algorithm), which reduces it to a collection of 3 3 uncou-pled linear systems, one at each grid point (wave number) in the discrete Fourier space. These 3 3 linear systems are easily solved, and then the inverse discrete Fourier transform (again implemented by the FFT algorithm) is used to obtain the solu-tion in the physical space. The above procedure is most naturally applied in the context of periodic boundary condisolu-tions, and for that reason we choose our fluid domain as a periodic box. Note, however, that there is no fundamental requirement of using a periodic domain in conjunction with the pIB method. Any fluid solver based on finite element or finite difference methods can be used, with whatever boundary conditions that solver can accommodate. We have emphasized the periodic case because of the efficiency that follows from the usage of the FFT algorithm for solving the linear systems that arise in our numerical scheme.

Next, we update the boundary velocity Vnþ1_{by solving Eq.}₍₂₀₎_{which is the discretization of Eq.}₍₁₁₎_{. The boundary}

veloc-ity Unþ1ðsÞ and the preliminary boundary position eYnþ1_{ðsÞ can be obtained by Eqs.}_{(18) and (19)}_{, respectively. Since ~}

_s

nþ1_{ðsÞ is}

deﬁned as in Eqs.(19)–(21)must be solved simultaneously in order to obtain Vnþ1ðsÞ and ~

r

nþ1_{ðsÞ. To see how this is done, we}

ﬁrst deﬁne

eFn_{¼ KðY}n

XnÞ RUnþ1þ cb

@4Yn

nþ1_: _ð31Þ

Using the relations j~

s

nþ1_{j ¼ 1; ~}

_s

nþ1

s ~

s

nþ1¼ 0, and ~

nþ1ss j~

s

nþ1s j

2

_r

_~nþ1_{. Thus Eq.}₍₃₁₎_{is a modiﬁed Helmholtz equation, which can be solved by a standard}

ﬁnite difference method. The condition number of the resultant matrix equation is approximately OðM2Þ, where M is the number of the markers. Once we have obtained ~

r

nþ1_{, the velocity V}nþ1_{can be obtained using Eq.}₍₂₉₎_.

Finally we update the conﬁguration of the two representations of the boundary XðsÞ and YðsÞ by Eqs.(22) and (23), which completes the evolution over one time step.

Note that, the present numerical scheme shown in(14)–(23)is in fact a combination of the forward and backward Euler methods so the accuracy is only first-order. (This is verified in Section4.2below.) Even though there are second-order IB methods available in the literature[6,19], they are only second-order accurate for problems with sufficiently smooth solu-tions, say an immersed structure with no zero thickness. Since this is not the case considered here, we could have no gain by using the second-order methods.

4. Numerical results

4.1. Single vesicle in a shear ﬂow

In this section, we introduce a two-dimensional computational model of a single vesicle suspended in a simple shear ﬂow. This problem was explored by several researchers through numerical simulations[18,38,36]in Stokes ﬂow.

We present the initial setting of our model and display the physical and computational parameters which are used in our numerical experiments.

Consider an incompressible viscous ﬂuid in a square box ð0; DÞ ð0; DÞ where D ¼ 80

l

m. The fluid-filled domain contains an immersed elastic vesicle which is a closed curve. The interior and exterior of the vesicle consist of the same fluid. The left panel ofFig. 2shows the initial configuration of a vesicle of an arbitrary shape. In order to apply a shear flow, it is natural to impose the background velocity field u0¼

c

ððy D=2Þ; 0Þ to the whole domain, where y is from 0 to D and

c

is the shear rate.

In this paper, however, since we use the periodic boundary condition for the computational domain, we impose the periodic shear ﬂow ﬁeld, u0ðyÞ ¼

c

2Dpsinð2

p

y=DÞ; 0

instead, see the left panel ofFig. 2.

In our simulations, we use the following parameters: the ﬂuid density

q

is 1:0 g=cm3_{and the viscosity}

_l

_{is 0.01 g/(cm s),}

which is the same as those of water, the bending coefﬁcient cbis 1010dyne cm, the time stepDt ¼ 4 108s, and the spatial

mesh width h ¼Dx ¼Dy ¼ 80=256

l

m. Following the analysis and scaling in[18,36], we deﬁne the length scale R0¼ L=2

p

where L is the perimeter of the vesicle boundary and the time scale

s

¼

l

R3

0=cb. We generally use R0¼ 10

l

m, which induces

an intrinsic time scale

s

¼ 0:1 s. We vary the initial conﬁguration of the vesicle and the shear rate

c

to explore the difference of the motion of a vesicle. These variations introduce two important dimensionless parameters; namely, the reduced area V ¼ A

pR2 0

¼4Ap

L2 where A is the area of the vesicle, and the dimensionless shear rate

v

¼

c

l

R 3

0=cb¼

cs

¼ 0:1

c

.

One important computational issue is how to choose the penalty coefﬁcients K and R appropriately. Here, we offer a prac-tical guide as follows. First, we choose some allowable maximum deviations for two representations jYðs; tÞ Xðs; tÞj and

0 20 40 60 80 0 10 20 30 40 50 60 70 80 μm μ m 0 20 40 60 80 0 10 20 30 40 50 60 70 80 μm μ m

Fig. 2. The conﬁguration of a single vesicle suspended in the shear ﬂow u0¼c2Dpsinð2py=DÞ; 0

at t ¼ 0 (left) and t ¼ 0:04 (right) in dimensionless time. Here the reduced area of the vesicle is V ¼ 0:51, and the shear rate isv¼ 250. The right ﬁgure includes the velocity ﬁeld on the vesicle boundary which shows a tank-treading tangential motion of the steady state vesicle.

(7)

jVðs; tÞ Uðs; tÞj, respectively. Then, we adjust(increase) the values of K and R until they are large enough that these allowable deviations will not be exceeded. However, during each run, the time step sizeDt must be adjusted accordingly with the pen-alty parameters K and R to ensure the numerical stability. As discussed in[5,30], a useful time step stability constraint arising from solving Eqs.(14)–(17)can be chosen as

D

t < ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi R2þ 2KC p R K ; ð32Þ

where C is a problem dependent constant of order one. A more concrete stability analysis has been characterized in[21]. (Notice that, the roles K and R played in the present work are similar to the ones

a

and b in[5,30]. The integral feedback term in[5,30]involving the ﬂuid velocity and desired boundary velocity can lead to the present form of Eq.(10)easily.) One should notice that the present time step constraint is more restrictive than the above criterion since we solve the Eqs.(28)– (31)to obtain the surface tension mainly in an explicit way. Based on the above process, we choose K ¼ 2 109dyne=cm3

and R ¼ 4 103dyne s=cm3_.

It is known that the vesicle in a simple shear flow undergoes a tank-treading tangential motion at its equilibrium config-uration. The right panel ofFig. 2shows a simulation result in which the vesicle configuration and the velocity field on the vesicle boundary are drawn at time t ¼ 0:04. The top panel ofFig. 3shows a more detailed motion of the single vesicle including the streamlines of the velocity field at certain times. The surface tensions

r

ðs; tÞ at the corresponding times, which are computed counterclockwise along the vesicle boundary starting at a boundary point on the long axis, are plotted in the bottom panel. The surface tension also reaches the equilibrium state and remains almost unchanged as time evolves.

Notice that, inFig. 3, we have observed that there are some sawtooth oscillations for the surface tension along the vesicle boundary at t ¼ 0:036 when the vesicle conﬁguration is steady. This oscillatory behavior has also been found in previous study[38]where the point-wise collocation method is applied. In order to remove the oscillations, one might apply the numerical smoothing so that the high frequency modes can be eliminated. Here, however, these sawtooth oscillations are generally low amplitudes, and the surface tension remains with sufﬁcient accuracy as long as the vesicle maintains a com-pact shape. The convergence study shown below strongly suggests that the pIB method is convergent.

The motion of a steady state vesicle can be characterized by both the inclination angle h between the long axis of the ves-icle and the ﬂow direction, and the frequency

x

of the revolution. These two values have been found to be strongly depen-dent on the reduced area V but independepen-dent of the dimensionless shear rate

v

[18,36]. This is veriﬁed inFig. 4which shows the inclination angle (left) and the revolution angular velocity (right) in terms of reduced area when the vesicles are at their steady shapes. We observe that both the inclination angle and angular velocity increase with the reduced area V but are nearly independent of the shear rate

v

, which conﬁrms the observation from previous results in the literature[18,36].

Fig. 5shows the shape dynamics of a vesicle suspended freely in a quiescent flow. In the absence of the inextensibility constraint, the equilibrium shape of a vesicle is a circle. However, the equilibrium shape of an inextensible vesicle can be different from a circle. It is also known that the equilibrium shape is dependent only on the reduced area V of the vesicle but independent of the material properties of the vesicle and the surrounding fluid[31,36].Fig. 5shows the motion of ves-icles which have three different initial shapes (see the first column) but with the same reduced area V ¼ 0:33. The first row has

l

¼ 0:01 g=ðcm sÞ and cb¼ 1010dyne cm, the second row has

l

¼ 0:01g=ðcm sÞ and cb¼ 1011dyne cm, and the third

0 10 20 30 40 50 60 0 0.02 0.04 vesicle length (_μm) tension σ (s) (dyne) t=0.0004s 0.0012s 0.0036s

Fig. 3. The motion of a single vesicle and streamlines (top) and the surface tensionrðs; tÞ (bottom) along the vesicle at times 0.004, 0.012, and 0.036. The reduced area V ¼ 0:51, and the shear ratev¼ 250. We can observe that, once the vesicle reaches its steady conﬁguration, it goes through a tank-treading tangential motion.

(8)

row has

l

¼ 0:001 g=ðcm sÞ and cb¼ 1010dyne cm. We can observe that, independent of the bending rigidity of the vesicle

and the ﬂuid viscosity, the three vesicles have almost the same equilibrium shape. The only difference is the time elapsed to reach the equilibrium conﬁguration.

Even though a vesicle in a quiescent ﬂow has its equilibrium shape similar to those inFig. 5depending on the reduced area V, the equilibrium shape could be an ellipse, if the initial shape is an ellipse with a small aspect ratio and the bending rigidity is small[33]. We verify this fact fromFig. 6which shows the motion of vesicle with a shape of an ellipse with the aspect ratio 0.25. We use two different bending coefﬁcients: cb¼ 1010dyne cm (top) and 1014dyne cm (bottom). Whereas

the ellipse with a larger bending rigidity changes into a shape similar to those inFig. 5(top), the ellipse with a smaller bend-ing rigidity does not change its initial shape (bottom).

4.2. Convergence study

In order to provide further veriﬁcation that the pIB method correctly solves the inextensible vesicle dynamics, we perform a convergence study. We consider a single vesicle in a shear ﬂow with shear rates

v

¼ 10 and

v

¼ 250, seeFig. 2. Now we choose the mesh sizes of the domain N ¼ 128, 256, 512, and 1024 so that the corresponding mesh width is h ¼ D=N. We also

0.6 0.7 0.8 0.9 0.08 0.1 0.12 0.14 0.16 0.18 0.2 θ /π reduced area V 0.6 0.7 0.8 0.9 0.25 0.3 0.35 0.4 0.45 ωa reduced area V χ=100 χ=50 χ=10 χ=5 χ=1 χ=100 χ=50 χ=10 χ=5 χ=1

Fig. 4. The inclination angles (left) and the angular velocities (right) of vesicles with different reduced areas in shear ﬂows with different shear rates. We set the angular velocityxa¼

ﬃﬃ 2 p

x

cR0, where xis the mean angular velocity over the vesicle boundary at the equilibrium conﬁguration. We observe that both h and

xaare nearly independent of the dimensionless shear rate. This result is in agreement with[18,36].

t=0.0s t=0.032s t=0.064s t=0.16s

t=0.0s t=0.064s t=0.128s t=0.512s

t=0.0s t=0.016s t=0.032s t=0.096s

Fig. 5. The motion of three vesicles of the same reduced area with different bending rigidities and different ﬂuid viscosities at some chosen times. The equilibrium shapes of the three vesicles are almost the same, i.e., it is independent of ﬂuid viscosity and the material property.

(9)

chooseDs andDt proportional to h, so that each factor of two in refinement of the fluid mesh width is accompanied by the same factor of refinement for the boundary mesh and the time step duration. The penalty parameters are defined as K ¼ 104_{N dyne=cm}3_{and R ¼ 0:005 N dyne s=cm}3_{. i.e., both parameters increase as the mesh width and time step are refined.}

In order to make both penalty force terms in Eq.(10)equally important, the magnitude of K is obtained by scaling a time step factor of R. As we shall see inFig. 7, the distance between two represented markers X and Y gets smaller when the grid is reﬁned (the penalty parameters K and R thus get larger).

Since the vesicle boundary is inextensible and the ﬂuid is incompressible, the perimeter of the vesicle and the enclosed area should remain as constants as time goes on. Let L0ðXÞ and LTðXÞ be the perimeters of the vesicle boundary X at the initial

time and ﬁnal time T ¼ 0:04, respectively. The relative error of the perimeter is deﬁned as jLTðXÞ L0ðXÞj=L0ðXÞ, and the t=0.0s t=0.064s t=0.128s t=0.32s

t=0.0s t=0.064s t=0.128s t=0.32s

Fig. 6. The motion of vesicles with different bending rigidity. The initial ellipse has the aspect ratio 0.25. Whereas the ellipse with a larger bending rigidity cb¼ 1010dyne cm changes its shape (top), the ellipse with a smaller one cb¼ 1014dyne cm keeps its initial shape (bottom).

1 2 3 4 x 10−3 1.8 2 2.2 2.4 2.6 2.8 1 2 3 4 x 10−3 1.8 2 2.2 2.4 2.6 2.8 time (s) 1 2 3 4 x 10−3 1.8 2 2.2 2.4 2.6 2.8 time (s) ratio (velocity) 1 2 3 4 x 10−3 1.8 2 2.2 2.4 2.6 2.8 time (s) ratio (position)

Fig. 7. Convergence ratios of the computed velocity field uðx; tÞ of the fluid (top) and the configuration of the vesicle boundary Xðs; tÞ (bottom). The convergence ratios (defined in the text) are plotted as functions of time for two different cases:v¼ 10 (left) andv¼ 250 (right). In each panel, the dashed line is the convergence ratio obtained using the grids N ¼ 128; 256; 512, and the solid line is the ratio obtained with N ¼ 256; 512; 1024. The convergence ratio is nearly two (first-order accuracy) for all the different cases.

(10)

relative error of the area is jATðXÞ A0ðXÞj=A0ðXÞ where the functional A assigns to the vesicle X its enclosed area.Table 1

shows the relative errors of the perimeter and the area of the region bounded by the vesicle boundary in the two cases:

v

¼ 10 (left) and

v

¼ 250 (right). We observe that the two errors in both cases indicate ﬁrst-order convergence. Besides, those errors are reasonably small and comparable to the results in[36].

Fig. 7shows the convergence ratios of the computed ﬂuid velocity uðx; tÞ and the vesicle boundary Xðs; tÞ. Since we do not know the exact solution of the problem, the estimation of the convergence ratio requires three numerical solutions for three consecutive N’s. Let ðuN;

v

NÞ be the velocity ﬁeld, and let k k2be L2norm. The top panels ofFig. 7shows the convergence

versus time for each of the cases N ¼ 128 (dotted line) and 256 (solid line). The bottom panels ofFig. 7shows the corresponding convergence ratios of the configuration of the vesicle boundary Xðs; tÞ versus time. (Note that, the error estimation is calculated in a similar manner to that of the fluid velocity above.) One can see from this figure that the convergence ratio for the fluid velocity and the vesicle boundary are all around two, which indicates that the present method is first-order accurate. This is the typical behavior in accuracy for the IB method as applied to problems with thin elastic boundaries (for second-order IB method in the case of an immersed elastic structure of finite thickness, see[6,14]).

4.3. Multiple vesicles in a shear ﬂow

As in[36], we now consider the case of multiple vesicles suspended in a shear flow. Since the computational time of our method is mostly consumed by solving the fluid equations, and there is no direct interaction between two different vesicles (the interaction occurs only through the fluid equations), the total computational cost is almost independent of the number of vesicles. However, since the increase of the number of vesicles simulated requires to increase the size of the computational domain, in order to solve the fluid equations with the same resolution as in the case of one vesicle, the computing time will still increase proportionally to the size of the computational domain.

We first consider the dynamics of four vesicles under a shear flow in a fluid domain 2D D.Fig. 8shows the motion of the four vesicles including streamlines of the velocity field at some chosen times. The shear rate is

v

¼ 250 as in[36]. The group of vesicles turns together like one single vesicle initially (from 0.016 to 0.024), and then they start to separate (from 0.04 to 0.048). Once they separate, the shear ﬂow will drive them farther away from each other. The result is consistent with the one obtained in[36].

We now simulate the dynamics of a larger number (56) of interacting vesicles. As explained before, the computing time of our method is dependent mostly on the number of spatial grid. We use the domain ð0; 4DÞ ð0; 4DÞ which is sixteen times larger than the case shown inFig. 2to get the same resolution. This amount of the increasing domain requires about 64 times more computing time than the single vesicle case. The cost of our method is comparable to that of the boundary integral method incorporating fast multiple method which is Oðmn log mÞ, where n is the number of vesicles and m is the number of the marker points on each vesicle[36].

Fig. 9shows the motion of the vesicles at several times. While the group of vesicles behaves like a single vesicle and moves following the shear ﬂow, each vesicle changes its shape into a relaxed one. Note that, we allow the vesicles to move out of the domain, in which case we should handle the data outside the domain by duplicating them into the domain in a periodic way, see the vesicles leaving the left (or right) side of the domain and coming into the right (or left) side of the domain.

4.4. Binary-component vesicle

As a ﬁnal example, we consider a binary-component vesicle suspended freely in a viscous ﬂuid. The binary-component vesicle is also inextensible and the bending rigidity and the spontaneous curvature of the vesicle boundary are not constants, but depend on the concentration of the different lipid phases. We consider the two component system with mass density

q

_A

0ð/ÞÞ 2 ds; ð33Þ Table 1

The relative errors of the perimeter and the enclosed area measured at time 0.04. The left table is forv¼ 10 and the right table is forv¼ 250.

N v¼ 10 v¼ 250

jLT L0j=L0ðXÞ jAT A0j=A0ðXÞ jLT L0j=L0ðXÞ jAT A0j=A0ðXÞ

128 1.13e6 3.14e5 1.52e5 5.61e5

256 7.12e7 6.21e6 1.72e5 3.64e5

512 4.05e7 4.61e7 1.03e5 1.36e5

(11)

where

j

¼ @2_X @s2

is the curvature, and

j

0is the spontaneous curvature which represents the asymmetry nature of the

mem-brane. Note that here we use the arc-length parameter s. Taking the variational derivative of the above modiﬁed bending energy, we obtain the following bending force

t=0s t=0.0008s

t=0.0016s t=0.0024s

t=0.0032s t=0.004s

Fig. 8. The motion of four vesicles in a shear ﬂow at some chosen times with streamlines depicting the velocity ﬁeld.

t=0s 0.0008s 0.0016s

0.0024s 0.0032s 0.0040s

0.0048s 0.0056s 0.0064s

(12)

Fb¼ dEb dX¼ @2 @s2 cbð/Þ

j

0ð/Þ @2X @s2 1 0 B @ 1 C A@ 2_X @s2 0 B @ 1 C A ð34Þ

One can easily see that when the case of constant cband zero spontaneous curvature

j

0¼ 0, the above bending force

recov-ers to the original one in Eq.(8). So the total force F in Eq.(5)now becomes

Fðs; tÞ ¼ @ @sð

r

ðs; tÞ

s

ðs; tÞÞ þ @2 @s2 cbð/Þ

j

0ð/Þ @2X @s2 1 0 B @ 1 C A@ 2_X @s2 0 B @ 1 C A: ð35Þ

Since the phase concentration / affects the bending rigidity cbð/Þ and the spontaneous curvature

j

0ð/Þ of the membrane, in

this work, we simply adopt the same idea used in interfacial ﬂow in which the surfactant concentration affects the magni-tude of surface tension. Since the surfactant stays in the interface (insoluble), one has to solve the surfactant equation along the interface[20,35]. Here, the order parameter / plays a similar role as the surfactant in surface tension driven ﬂows. Thus, one needs to solve the phase concentration equation on the vesicle boundaryC. In order to separate those two components on the vesicle boundary, we adopt the Cahn–Hilliard type continuum approach[1], meaning that, the time evolution of /ðs; tÞ is governed by @/ @t ¼ M/

D

sw ð36Þ w¼1

e

G 0 ð/Þ

e

D

sw; ð37Þ

where the surface Laplacian is deﬁned asDs¼@s@ @s@

j@X @sj =@X @s

. Here, w represents the chemical potential. Notice that, a more complicated chemical potential form has been derived in[32]to preserve the thermodynamically consistent relationship be-tween the mass ﬂux term and the force term. The function Gð/Þ ¼1

4/ 2

ð1 /Þ2is a double well potential that describes the tendency of the phase concentration / into the two stable phase separation / ¼ 0 or / ¼ 1, and G0

D

t M/ @ @s @wnþ1 @s , @Y0 @s !, @Y0 @s ; ð41Þ wnþ1¼1

e

G 0 ð/n_Þ

_e

@ @s @/nþ1 @s , @Y0 @s !, @Y0 @s : ð42Þ

The spatial discretization of these two equations induces a linear system for /nþ1_{and w}nþ1_{. Since} @Y0 @s

is initially deﬁned and does not change in time, we only need to invert the resulting matrix once and save it for the later usage. Once we obtain the new /nþ1_{, we can substitute it into the functions of c}

bð/nþ1Þ and

j

0ð/nþ1Þ so that Eq.(40)can be solved numerically using the

same form of discretization as before where Y in the right-hand side is treated explicitly.

In our tests, we choose M/¼ 106, and

e

¼ 104. The spontaneous curvature is

j

0ð/Þ ¼R10ð5ð1 /Þ þ 0:1/Þ, where the

length scale R0is 103, and the bending rigidity is cbð/Þ ¼ c0ðð1 /Þ þ 0:5/Þ, where the coefﬁcient c0is 1010. Notice that,

the above special form of cband

j

0are adopted from[32]except for some scaling factors. The initial shape of the

binary-component vesicle is the same as the one shown inFig. 2. As in[32], the initial concentration / is a mixture of two phases with the following perturbed form:

(13)

/ðs; 0Þ ¼ /þ 0:001ðcosð2sÞ þ cosð4sÞ þ cosð6sÞÞ; ð43Þ

where /is the average concentration.

Fig. 10shows the motion of the binary-component vesicle suspended freely with streamlines depicting the velocity field and the concentration of the lipid phase /ðs; tÞ at some chosen times. The cases of four different average concentrations have been simulated and ordered from top to bottom row as / ¼ 0 (first row), and /¼ 0:25; 0:5, and /¼ 0:75 (last two rows). Notice that, the case of / ¼ 0 (first row) shows the motion of a single component vesicle which has a constant concentration /¼ 0 all the time. The concentration of the lipid phase /ðs; tÞ is drawn in terms of s starting at the point marked by on each vesicle and in the counterclockwise way along the vesicle.

One can see that, as time evolves, the phase separation occurs and yields two larger regions: one region of / 1 at the vesicle top and bottom parts where the vesicle boundary has smaller curvature, and the other region / 0 phase at the two tips of the vesicle where the vesicle boundary has larger curvature. This is exactly what we expect due to the form of spon-taneous curvature. One can conclude from the ﬁgure that, as the average concentration gets lower, the curvature of the ves-icle gets larger due to the spontaneous curvature. That is, when the average concentration /becomes larger, the vesicle becomes ﬂatter in general.

5. Summary and conclusions

We have developed a penalty immersed boundary (pIB) method to simulate the dynamics of inextensible vesicle inter-acting with an incompressible ﬂuid. In order to take into account the inextensibility constraint of the vesicle, we virtually separate the dynamics of the inextensible vesicle from the ﬂuid equations using the penalty force, and then solve the system of equations governing the motion of the inextensible boundary. In this separate system for the vesicle motion, the inexten-sibility constraint simply produces a simple one-dimensional second-order elliptic equation for the unknown surface tension.

We have validated the method by simulating the motion of a single vesicle both in a shear flow and in a quiescent flow. We have shown that the vesicle in a simple shear flow undergoes a tank-treading tangential motion at its equilibrium state, and that its inclination angle h and frequency

x

of the revolution are strongly dependent on the reduced area V but inde-pendent of the shear rate

v

, as shown in other literature. The fact that the shape of a vesicle in a quiescent ﬂow depends only on the reduced area V of the vesicle but independent of the material properties of the interface and surrounding ﬂuid

0 2 4 6 0 0.5 1 0 2 4 6 0 0.5 1 0 2 4 6 0 0.5 1 0 2 4 6 0 0.5 1 0 2 4 6 0 0.5 1 0 2 4 6 0 0.5 1 0 2 4 6 0 0.5 1 0 2 4 6 0 0.5 1 0 2 4 6 0 0.5 1 0.00016s 0 2 4 6 0 0.5 1 0.00048s 0 2 4 6 0 0.5 1 0.0096s 0 2 4 6 0 0.5 1 0.016s

Fig. 10. The motion of binary-component vesicles suspended freely with streamlines depicting the velocity ﬁeld and the concentration of the lipid phase /ðs; tÞ at some chosen times. Four different concentration cases: constant concentration / ¼ 0 (ﬁrst row), /¼ 0:25 (2nd–3th rows), /¼ 0:5 (4th–5th rows) and /¼ 0:75 (bottom two rows).

(14)

has also been verified. Additional validation has been provided in the form of a convergence study, which confirms the ex-pected first-order accuracy of the scheme.

For multiple vesicle cases, the total computing time is dependent mostly on the number of the computational grids but is almost independent of the number of vesicles. The cost of our method is comparable to that of other numerical methods such as the boundary integral method. Moreover, we have also extended the method to solve a simpliﬁed binary-component ves-icle model.

Future work will include the extension of the methodology introduced here to the study of three-dimensional vesicle dynamics. The application of the present method to some biological systems such as red blood cells and drug-carrying cap-sules will also be the subject of future work. Unlike the boundary integral method, the present method can be applied to the high Reynolds number ﬂow in which the inertial effect cannot be neglected.

Acknowledgment

Y Kim is supported by the Chung-Ang University Research Grants in 2010. M.-C. Lai is supported in part by National Sci-ence Council of Taiwan under research Grant NSC-97-2628-M-009-007-MY3 and MoE-ATU project.

References

[1] J.W. Cahn, J.E. Hilliard, Free energy of a nonuniform system I interfacial free energy, J. Chem. Phys. 28 (1958) 258–267.

[2] C. Duncan, G. Zhai, R. Scherer, Modeling coupled aerodynamics and vocal fold dynamics using immersed boundary methods, J. Acoust. Soc. Am. 120 (2006) 2859–2871.

[3] L.J. Fauci, C.S. Peskin, A computational model of aquatic animal locomotion, J. Comput. Phys. 77 (1988) 85–108.

[4] A.L. Fogelson, C.S. Peskin, A fast numerical method for solving the three-dimensional Stoke’s equations in the presence of suspended particles, J. Comput. Phys. 79 (1988) 50–69.

[5] D. Goldstein, R. Handler, L. Sirovich, Modeling a no-slip ﬂow with an external force ﬁeld, J. Comput. Phys. 105 (1993) 354–366.

[6] B.E. Grifﬁth, C.S. Peskin, On the order of accuracy of the immersed boundary method: higher order convergence rates for sufﬁciently smooth problems, J. Comput. Phys. 208 (2005) 75–105.

[7] K.H. de Haas, C. Blom, D. van den Ende, M.H.G. Duits, J. Mellema, Deformation of giant lipid bilayer vesicles in shear ﬂow, Phys. Rev. E 56 (1997). [8] W. Helfrich, Elastic properties of lipid bilayers: theory and possible experiments, Z. Naturforsch C 28 (1973) 693–703.

[9] J.J. Heys, B. Knott, T. Gedeon, Y. Kim, Modeling arthropod ﬁliform hair motion using the penalty immersed boundary method, J. Biomech. Eng. 41 (2008) 977–984.

[10] W. Huang, S. Shin, H. Sung, Simulation of flexible filaments in a uniform flow by the immersed boundary method, J. Comput. Phys. 226 (2007) 2206– 2228.

[11] E. Jung, C.S. Peskin, Two-dimensional simulations of valveless pumping using the immersed boundary method, SIAM J. Sci. Comput. 23 (2001) 19–45. [12] V. Kantsler, V. Steinberg, Orientation and dynamics of a vesicle in tank-treading motion in shear ﬂow, Phys. Rev. Lett. 95 (2005).

[13] Y. Kim, C.S. Peskin, 2-D parachute simulation by the immersed boundary method, SIAM J. Sci. Comput. 28 (6) (2006). [14] Y. Kim, C.S. Peskin, Penalty immersed boundary method with an elastic boundary with mass, Phys. Fluids 19 (5) (2007).

[15] Y. Kim, C.S. Peskin, Numerical study of incompressible ﬂuid dynamics with nonuniform density by the immersed boundary method, Phys. Fluids 20 (2008) 062101.

[16] Y. Kim, S. Lim, A. Friedman, S. Raman, O. Simonetti, Blood ﬂow in a compliant vessel by the immersed boundary method, Ann. Biomed. Eng. 37 (5) (2009).

[17] Y. Kim, C.S. Peskin, 3-D parachute simulations by the immersed boundary method, Comput. Fluids 38 (2009) 1080–1090. [18] M. Kraus, W. Wintz, U. Seifert, R. Lipowsky, Fluid vesicles in shear ﬂow, Phys. Rev. Lett. 77 (1996).

[19] M.-C. Lai, C.S. Peskin, An immersed boundary method with formal second-order accuracy and reduced numerical viscosity, J. Comput. Phys. 160 (2000) 705–719.

[20] M.-C. Lai, Y.-H. Tseng, H. Huang, An immersed boundary method for interfacial ﬂows with insoluble surfactant, J. Comput. Phys. 227 (2008) 7279– 7293.

[21] C. Lee, Stability characteristics of the virtual boundary method in three-dimensional applications, J. Comput. Phys. 184 (2003) 559–591.

[22] J.S. Lowengrub, A. Ratz, A. Voigt, Phase-ﬁeld modeling of the dynamics of multi-component vesicles: spinodal decomposition, coarsening, budding and ﬁssion, Phys. Rev. E 79 (2009) 031926.

[23] A.J. Markvoort, R.A. Van Santen, P.A.J. Hilbers, Vesicle shapes from molecular dynamics simulations, J. Phys. Chem. B 110 (2006). [24] H. Noguchi, G. Gompper, Shape transitions of ﬂuid vesicles and red blood cells in capillary ﬂows, Proc. Natl. Acad. Sci. 102 (2005). [25] C.S. Peskin, The immersed boundary method, Acta Numer. 11 (2002) 479–517.

[26] C.S. Peskin, D.M. McQueen, Three dimensional computational method for flow in the heart: immersed elastic fibers in a viscous incompressible fluid, J. Comput. Phys. 81 (1989) 372–405.

[27] C.S. Peskin, D.M. McQueen, Fluid dynamics of the heart and its valves, in: Case Studies Math. Model.: Ecol. Physiol. Cell Biol., Prentice-Hall, Englewood Cliffs, NJ, 1996, pp. 309–337.

[28] C.S. Peskin, B.F. Printz, Improved volume conservation in the computation of ﬂows with immersed elastic boundaries, J. Comput. Phys. 105 (1993) 33– 46.

[29] C. Pozrikidis, The axis-symmetric deformation of a red blood cell in uniaxial straining Stokes ﬂow, J. Fluid Mech. 216 (1990) 231–254.

[30] E.M. Saiki, S. Biringen, Numerical simulation of a cylinder in uniform ﬂow: application of a virtual boundary method, J. Comput. Phys. 123 (1996) 450– 465.

[31] U. Seifert, Adhesion of vesicles in two-dimensions, Phys. Rev. A 43 (1991).

[32] J.S. Sohn, Y.-H. Tseng, S. Li, A. Voigt, J.S. Lowengrub, Dynamics of multi-component vesicles in a viscous ﬂuid, J. Comput. Phys. 229 (2010) 119–144. [33] P. Song, Numerical Simulations of Cylindrical Membranes, Ph.D. Thesis, Peking University, 2008.

[34] M.J. Stevens, Coarse-grained simulations of lipid bilayers, J. Chem. Phys. 121 (2004).

[35] H.A. Stone, A simple derivation of the time-dependent convective–diffusion equation for surfactant transport along a deforming interface, Phys. Fluids A 2 (1) (1990) 111–112.

[36] S.K. Veerapaneni, D. Gueyfﬁer, D. Zorin, G. Biros, A boundary integral method for simulating the dynamics of inextensible vesicles suspended in a viscous ﬂuid in 2D, J. Comput. Phys. 228 (2009) 2334–2353.

[37] S.K. Veerapaneni, D. Gueyfﬁer, G. Biros, D. Zorin, A numerical method for simulating the dynamics of 3D axis-symmetric vesicles suspended in viscous ﬂows, J. Comput. Phys. 228 (2009) 7233–7249.