An implicit immersed boundary method for three-dimensional fluid–membrane interactions

Journal of Computational Physics 228 (2009) 8427–8445

Contents lists available at ScienceDirect

Journal of Computational Physics journal homepage: www.elsevier.com/locate/jcp

An implicit immersed boundary method for three-dimensional fluid–membrane interactions D.V. Le a,b,e,*, J. White a,b, J. Peraire a,c, K.M. Lim a,d, B.C. Khoo a,d a

Singapore-MIT Alliance, 4 Engineering Drive 3, National University of Singapore, Singapore 117576, Singapore Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, 77 Massachusetts Avenue, MA 02139, USA c Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, 77 Massachusetts Avenue, MA 02139, USA d Department of Mechanical Engineering, National University of Singapore, Kent Ridge Crescent, Singapore 119260, Singapore e Institute of High Performance Computing, 1 Fusionopolis Way, #16-16 Connexis, Singapore 138632, Singapore b

a r t i c l e

i n f o

Article history: Received 14 October 2008 Received in revised form 30 April 2009 Accepted 20 August 2009 Available online 2 September 2009 MSC: 65M06 76D05 74F10 92C10 Keywords: Immersed boundary method Newton–Krylov method Navier–Stokes equations Membrane capsules Red blood cells

a b s t r a c t We present an implicit immersed boundary method for the incompressible Navier–Stokes equations capable of handling three-dimensional membrane–fluid flow interactions. The goal of our approach is to greatly improve the time step by using the Jacobian-free Newton–Krylov method (JFNK) to advance the location of the elastic membrane implicitly. The most attractive feature of this Jacobian-free approach is Newton-like nonlinear convergence without the cost of forming and storing the true Jacobian. The Generalized Minimal Residual method (GMRES), which is a widely used Krylov-subspace iterative method, is used to update the search direction required for each Newton iteration. Each GMRES iteration only requires the action of the Jacobian in the form of matrix–vector products and therefore avoids the need of forming and storing the Jacobian matrix explicitly. Once the location of the boundary is obtained, the elastic forces acting at the discrete nodes of the membrane are computed using a finite element model. We then use the immersed boundary method to calculate the hydrodynamic effects and fluid–structure interaction effects such as membrane deformation. The present scheme has been validated by several examples including an oscillatory membrane initially placed in a still fluid, capsule membranes in shear flows and large deformation of red blood cells subjected to stretching force. Ó 2009 Elsevier Inc. All rights reserved.

1. Introduction This paper considers an implicit immersed boundary method for simulating viscous incompressible flows with immersed elastic membranes. The immersed boundary (IB) method was originally developed by Peskin [36] to study the fluid dynamics of blood flow in the human heart. Peskin’s immersed boundary method has proven to be a very useful method for modeling fluid–structure interaction involving large geometry variations. The original method has been developed further and applied to many biological problems including platelet aggregation [14,15,52], the deformation of red blood cells in a shear flow [11], the swimming of bacterial organisms and others [10,13]. More details on the immersed boundary method can be found in [37] and the references therein. Typically, in the framework of the IB method, the elastic boundary is treated as a collection of elastic fibers. Alternatively, the elastic boundary is also modeled as an elastic membrane [11] or a thin shell [17]. These models are used to compute the * Corresponding author. Address: Institute of High Performance Computing, 1 Fusionopolis Way, #16-16 Connexis, Singapore 138632, Singapore. Tel.: +65 64191219; fax: +65 64674350. E-mail addresses: [email protected] (D.V. Le), [email protected] (J. White), [email protected] (J. Peraire), [email protected] (K.M. Lim), [email protected] (B.C. Khoo). 0021-9991/$ - see front matter Ó 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2009.08.018

8428

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

forces acting at the discrete nodes representing the immersed boundary. For most biological tissues, their membranes are stiff and therefore a small perturbation of the boundary can lead to large elastic forces. This causes a severe restriction in time step required to maintain stability of the immersed boundary method. Much effort has been made to analyze the stiffness of the IB method and remove this restriction [18,22,33,47,48]. Several semi-implicit and implicit methods have been developed to alleviate this problem [12,21,31,35,51]. Comparisons of the explicit method and the implicit methods in the context of moving immersed boundaries have been presented in [35,51]. In the context of immersed interface methods (IIM) a quasi-Newton method has been proposed [27–29] to improve time stability. However, for very stiff problems, small time steps are still required. In order to alleviate this problem, an unconditionally stable discretization of the immersed boundary equations has been proposed [34]. In this scheme, an approximate Newton solver is employed to advance the location of the boundary. The Newton method, however, requires a Jacobian matrix which is extremely expensive to compute explicitly for every time step. To avoid forming the Jacobian matrix explicitly, an iterative matrix-free method has been introduced in the context of the immersed continuum method [53]. Another strategy for solving the implicit IB equations involves deriving Schur complement equations by eliminating one or more of the unknown variables [33,35]. Several methods have been employed to solve the Schur complement equations such as fixed point methods, projection methods and Krylov-subspace methods. Recently, an efficient semi-implicit IB method with arclength–tangent angle formulation has been proposed for twodimensional Navier–Stokes equations [21]. In this formulation, an unconditionally stable semi-implicit discretization is derived and the small scale decomposition technique [20] is applied to the discretization. This semi-implicit scheme has much better stability property than the explicit scheme and therefore offers a substantial computational cost saving. We note that the stability of this scheme is weaker than the unconditionally stable scheme proposed in [34] because this scheme treats only the leading order term implicitly [21]. In the present paper, an implicit immersed boundary method is presented with vastly improved time step by using the Jacobian-free Newton–Krylov method (JFNK) [24] to advance the location of the elastic membrane implicitly. This matrixfree approach has many advantages. The most attractive is Newton-like nonlinear convergence without the cost of forming and storing the true Jacobian. In this Jacobian-free method, the search direction required for each Newton iteration is updated using the Generalized Minimal Residual method (GMRES) [44], which is a widely used Krylov-subspace iterative method. Each GMRES iteration only requires the action of the Jacobian in the form of matrix–vector products and therefore avoids the need of forming and storing the Jacobian matrix explicitly. The JFNK method has become established in computational fluid dynamics (CFD) to deal with the nonlinear convection term [3,6,16,25]. In this paper, we employ the JFNK method to advance the membrane location implicitly while still approximating the convection term explicitly. Once the location of the boundary is obtained, the elastic forces acting at the discrete nodes of the membrane are computed. In the present work, an elastic membrane model with bending stiffness proposed in [38,42] is employed. In our numerical studies, the immersed boundary method provides the means of calculating the hydrodynamics and fluid–structure interaction effects such as membrane deformation, and the finite element method with membrane model is used to calculate the elastic forces on the membrane. The remainder of this paper is organized as follow. In Section 2, we describe the governing equations, the immersed boundary algorithm, the discrete model of capsule membranes and the method for advancing the membrane evolution in time. We then present a detailed implementation of the present method in Section 3. In Section 4, some numerical examples are presented and finally, some conclusions are given in Section 5. 2. Numerical methods 2.1. Governing equations In a three-dimensional bounded domain X that contains an enclosed membrane CðtÞ, we consider the incompressible Navier–Stokes equations formulated in primitive variables, written as

qðut þ ðu  rÞuÞ ¼ rp þ lDu þ F; ru¼0

ð1Þ ð2Þ

with boundary conditions

uj@ X ¼ ub ;

ð3Þ

where u is the fluid velocity, p is the pressure, q and l are constant density and viscosity of the fluid, respectively. The effect of the membrane CðtÞ immersed in the fluid results in a singular force F which has the form

Fðx; tÞ ¼

Z

f ðq; r; s; tÞdðx  Xðq; r; s; tÞÞ dq d rds;

ð4Þ

CðtÞ

where ðq; r; sÞ are curvilinear coordinates attached to the membrane at a material point, Xðq; r; s; tÞ is the position at time t in Cartesian coordinates of the material point whose label is ðq; r; sÞ; x ¼ ðx; y; zÞ is spatial position, and f ðq; r; s; tÞ is the force

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

8429

strength. Here, dðxÞ is the three-dimensional Dirac function. The motion of the boundary can be determined by integrating the equation

dXðq; r; s; tÞ ¼ dt

Z

uðx; tÞdðx  Xðq; r; s; tÞÞ dx:

ð5Þ

X

2.2. Description of the IB and projection methods The immersed boundary method uses a set of N control points X l ; l ¼ 1; . . . ; N to represent the immersed boundary. The force density is computed at these control points and is distributed to the Cartesian grid points using a discrete representation of the delta function,

Fðxði; j; kÞ; tÞ ¼

N X

f l ðq; r; s; tÞDh ðxði; j; kÞ  X l ðtÞÞDqDrDs;

ð6Þ

l¼1

where f l ðq; r; s; tÞ is the force density at the control point X l whose label is ðq; r; sÞ; xði; j; kÞ and Fðxði; j; kÞÞ are the coordinate of grid point ði; j; kÞ and the force at that point, respectively. Dh ðxÞ is a three-dimensional discrete delta function,

Dh ðxÞ ¼

1 h

3

u

 x   y  z  u u ; h h h

ð7Þ

where h is the grid size, x; y and z are the Cartesian components of x and u is a continuous function which is taken from [43] and given by

8  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > 2 > > 16 5  3jdj  1  3ð1  jdjÞ ; 0:5 6 jdj 6 1:5; > <  uðdÞ ¼ 1  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 > 1 þ 3d þ 1 ; jdj 6 0:5; > 3 > > : 0; otherwise:

ð8Þ

Once the force density is computed at the control points and distributed to the grid, the Navier–Stokes equations with the forcing terms are then solved for the pressure and velocity field at the Cartesian grid points using the projection method [4]. Our numerical algorithm is based on the pressure-increment projection algorithm for the discretization of the Navier–Stokes equations. The spatial discretization is carried out on a standard marker-and-cell (MAC) staggered grid similar to that described in Kim and Moin [23]. Given the velocity un , the pressure pn1=2 and the forcing term F nþ1=2 , we compute the intermediate velocity u as follows:

 1 u  un l 2  1 rh u þ r2h un þ F nþ1=2 ; ¼ qðu  ruÞnþ2  GMAC pn2 þ Dt 2 u j@ X ¼ unþ1 ; b

q

ð9Þ

where the advective term is extrapolated using the formula,

3 1 ðu  rh uÞn  ðu  rh uÞn1 : 2 2

1

ðu  ruÞnþ2 ¼

ð10Þ

This intermediate velocity field, in general, does not satisfy the divergence-free condition (2). Therefore, we compute a pressure-increment /nþ1 and update the pressure and velocity field as

DMAC u ; n  r/nþ1 j@ X ¼ 0; Dt 1 ¼ u  DtGMAC /nþ1 ;

r2h /nþ1 ¼ q unþ1

ð11Þ ð12Þ

q

pnþ1=2 ¼ pn1=2 þ /nþ1 

l 2q

ðDMAC u Þ:

ð13Þ

We note that the above projection method is analogous to the pressure-increment projection method presented in [4]. In the above expressions, rh and r2h are the standard central difference operators, GMAC and DMAC are the MAC gradient and divergence operators, respectively [26]. The velocity field is then interpolated to find the velocity at the control points as,

uðX l ; tÞ ¼

X

3

uðxði; j; kÞ; tÞDh ðxði; j; kÞ  X l ðtÞÞh ;

i;j;k

and this velocity is used to advance the position of the membrane.

ð14Þ

8430

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

2.3. Discrete model of the capsule membrane 2.3.1. Membrane tensions Following the standard procedure [2,42], in order to keep track of the position of a material point on the membrane, let X and XðX; tÞ be its positions in the unstressed state and after deformation at time t, respectively. Let C be the deformation gradient tensor defined as

C¼

@X : @X

ð15Þ

Following [2], the surface deformation gradient tensor is defined as

n  Þ; A ¼ ðI  nnÞ  C  ðI  n

ð16Þ

 and n are the unit normal vector to the undeformed and deformed membranes, respectively. The left Cauchy–Green where n strain tensor is then defined as

B ¼ K 2 ¼ A  AT :

ð17Þ

The membrane tension tensor s is related to the surface strain energy function WðK1 ; K2 Þ and B by

s ¼ eK1



 @W @W Pþ B ; @ K1 @ K2

ð18Þ

where P ¼ I  nn and K1 ; K2 are the strain invariants

  1 1 log ðtrðBÞ2  trðB2 ÞÞ ; 2 2  1 2 1 K2 ¼ k1 þ k22  1 ¼ trB  1: 2 2

K1 ¼ log k1 k2 ¼

ð19Þ ð20Þ

Here, the eigenvalues k1 ; k2 of K are the principle planar stretches. The strain energy function for a neo-Hookean membrane is given by

W¼

 Eh ð2K2 þ e2K1  1Þ; 6

ð21Þ

 is the membrane thickness. For a thin membrane, where E is the Young’s modulus and h

W¼

Es ð2K2 þ e2K1  1Þ; 6

ð22Þ

 is the surface elastic modulus, to be distinguished from the volume modulus of elasticity of a three-dimenwhere Es ¼ Eh sional material [40]. Alternative constitutive equations for hyperelastic materials such as biological membranes can be used. For example, the Yeoh form of strain energy function [54] given by

W ¼ C 10 ð2K2 þ e2K1  1Þ þ C 30 ð2K2 þ e2K1  1Þ3 ;

ð23Þ

where C 10 and C 30 are constants, can be used to model red blood cell (RBC) deformation. Another strain energy function for the red blood cell membrane which conforms to deformation measurements for human RBCs has been developed in [46]

W¼

C 11 C 21 2K1 ð2K2 ðK2 þ 1Þ þ 1  e2K1 Þ þ ðe  1Þ2 ; 4 8

ð24Þ

where C 11 and C 21 are model constants and typically C 11  C 21 . 2.3.2. Membrane bending moments The bending moments developing along the edges of a membrane patch in the deformed state depend on the instantaneous edge curvature, as well as on the edge curvature in the rest state. For sufficiently small bending deformations, but not necessarily small in-plane deformations, the bending moments may be approximated with the linear constitutive equation

  m ¼ jB j  jRm P ;

ð25Þ R m

where jB is the bending modulus, j is the Cartesian curvature tensor and j are the position dependent mean curvatures of the resting configurations. More details on the derivation of Eq. (25) can be found in [40,41]. Finally, the local force density f exerted by the membrane is given by

f ¼ rs  ðs þ qnÞ;

ð26Þ

where rs is the surface gradient operator rs ¼ ðI  nnÞ  r and q is the transverse shear tension. The expression for the transverse shear tension vector in terms of the surface divergence of the tensor of the bending moments is

q ¼ ½ðP  rÞ  m  P ¼ tr½ðP  rÞm  P:

ð27Þ

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

8431

2.4. Advancing the membrane 2.4.1. Explicit method The simplest version of the immersed boundary method is the explicit scheme in which the forward Euler method is used to advance the location of the membrane. This explicit scheme includes the following steps:  Advance the control points as

X nþ1 ¼ X n þ Dtun ðX n Þ:

ð28Þ

nþ1

 Use X to compute the boundary force as described in Section 2.3. Distribute the boundary force to the nearby Cartesian grid points using discrete delta function.  Solve the Navier–Stokes equations using the projection method to calculate the velocity field and pressure field.  Interpolate the velocity field to determine the velocity at the control points unþ1 ðX nþ1 Þ. This method is simple, but for problems where the membrane is stiff, a very small time step is required to maintain stability. Details on the stability issues of the explicit immersed boundary method are discussed in [51]. 2.4.2. Quasi-Newton method In order to improve the time step, one could use an implicit method. The position of the control points is updated using the trapezoidal rule

1 X nþ1 ¼ X n þ Dtðun ðX n Þ þ unþ1 ðX nþ1 ÞÞ: 2

ð29Þ

Eq. (29) is implicit and couples the motion of the membrane with the solution at all grid points. Therefore at each time step, we need to solve a nonlinear system of equations for the position of the control points of the form

gðX nþ1 Þ ¼ 0;

ð30Þ

1 gðXÞ ¼ X  X n  Dtðun ðX n Þ þ unþ1 ðXÞÞ: 2

ð31Þ

where

Normally this nonlinear system of equations is solved by a Newton’s method in which the Jacobian matrix is required,

1 JðXÞ ¼ g 0 ðXÞ ¼ I  Dtu0 ðXÞ: 2

ð32Þ

Starting from X ð0Þ a Newton iteration involves a sequence of linear systems given by

JðX ðkÞ ÞdX ðkÞ ¼ gðX ðkÞ Þ;

ð33Þ

X ðkþ1Þ ¼ X ðkÞ þ dX ðkÞ ;

ð34Þ

k ¼ 0; 1; . . .

However, the evaluation of the Jacobian matrix and its inverse can be time consuming for many applications. A quasi-Newton method has been devised wherein the inverse Jacobian matrix J 1 is replaced by a suitable matrix H, which is easy to compute. The matrix H at iteration (k + 1)th, H kþ1 , is updated from the matrix H k at the previous iteration as follows

   qT H k q p pT ð1  hk Þ c hk  T H kþ1 ¼ ck H k þ 1 þ ck hk k T k kT k  ck T H k qk  qTk H k  Tk p q H k þ H k qk pTk ; qk H k qk pk qk pk qk p k qk k k

ð35Þ

where

pk :¼ X ðkþ1Þ  X ðkÞ ;

qk :¼ g ðkþ1Þ  g ðkÞ ;

and the parameters ck P 0; hk P 0. With ck ¼ 1; hk ¼ 1, we have the rank two method of Broyden, Fletcher, Goldfard and Shanno (BFGS method) [49]. In practice, the BFGS method requires only a few iterations to converge as the solution at the previous time step provides a very good initial guess for the iteration. This has been observed in [27–29] for two-dimensional problems in the context of the immersed interface method. Note that the BFGS method needs to form and store the inverse Jacobian matrix H, and consequently it requires a large storage in three-dimensional problems especially for those with multiple membranes in the fluid domain. The quasi-Newton method was not used in the present work for three-dimensional problems because of the limitations in computer memory. In order to avoid such large storage, we employ a Jacobianfree Newton–Krylov method as described in the next section. 2.4.3. Jacobian-free Newton–Krylov method An alternative approach to solve (30) is the Jacobian-free Newton–Krylov (JFNK) method [24]. In the JFNK approach, the GMRES method is used to solve the sequence of linear systems given by Eq. (33). The GMRES method requires the action of the Jacobian only in the form of a matrix–vector product, which may be approximated by [5,7]:

8432

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

J ðX Þy 

gðX þ yÞ  gðXÞ



ð36Þ

;

where  is a small perturbation. The error of this approximation is proportional to . In [7],  was set equal to something larger than the square root of machine error ðmach Þ. Another effective formula for the evaluation of  is

¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1 þ kXkÞmach : kyk

ð37Þ

More details on the Jacobian-free Newton–Krylov can be found in [24] and the references therein. 3. Implementation The present implicit immersed boundary method involves the following steps: 1. Trace the immersed boundary using a collection of control points that define an unstructured grid of triangular elements as shown in Fig. 1. 2. Compute the surface forces at the control points. 3. Solve the Navier–Stokes equations to compute the velocity and pressure fields. 4. Advance the position of the control points. In the numerical implementation, the immersed boundary is discretized into a mesh of six-node curved triangles. A function uðn; gÞ defined over a triangle is approximated as

uðn; gÞ ¼

6 X

ui Ni ðn; gÞ;

ð38Þ

i¼1

where ðn; gÞ are the local parametric coordinates, ui is the value of u at node i and N i ðn; gÞ are the basis functions for a quadratic six-node triangular finite element [39]. To evaluate the membrane tension tensor s, one needs to calculate the left Cauchy–Green strain tensor, which is determined from the surface deformation gradient tensor, A. For each element, the surface deformation gradient tensors at the nodes are obtained by solving the following system of equations,

A

@X @X ¼ ; @n @n

A

@X @X ¼ ; @g @g

 ¼ 0; An

ð39Þ

at each node of the element [42]. The components of A are then averaged over the elements sharing the node to obtain a smooth distribution. Once the deformation tensor has been determined, the strain and stress tensors follow from (17) and (18). The transverse shear tension is computed from the bending moment, which is determined from the Cartesian curvature tensor, j. To evaluate the curvature tensor j at a point, one needs to solve

@X @n j¼ ; @n @n

@X @n j¼ ; @g @g

n  j ¼ 0;

ð40Þ

at each element node and then average over the elements sharing that node [40]. Having computed the curvature tensor, one can evaluate the bending moments m using the constitutive equation (25). To evaluate the transverse shear tension, we need to compute the surface gradient of the bending moments denoted by K ðP  rÞm, using the relations

Fig. 1. (a) Discretization of a sphere. (b) Discretization of a biconcave shape.

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

@X @m K ¼ ; @n @n

@X @m K ¼ ; @g @g

n  K ¼ 0:

8433

ð41Þ

The evaluation of K at each element node requires solving nine systems of 3 linear equations [40]. The components of K are also averaged over the elements sharing the node, and the transverse shear tension can then be calculated using Eq. (27). To evaluate the force density, we adopt the method developed in [38] where the average value of the force is evaluated on an element. The average value of the force on the element En with area Sn enclosed by the contour C n is computed from the line integral

Df ¼

1 Sn

I

½b  s þ ðb  qÞn dl;

ð42Þ

Cn

where l is the arclength along the element contour, b t n is the cross-product of the unit tangential vector t along the contour C n and the surface unit normal vector n. Once the forces are computed over all elements and distributed to the grid, the Navier–Stokes equations with the forcing terms are then solved for the pressure and velocity at the Cartesian grid points. This velocity field is interpolated to find the velocity at the control points, and the position of the control points is updated implicitly using the trapezoidal rule (29). In summary, given the location of the control points, X n , the velocity field, un and the pressure field pn1=2 , the process of computing the new velocity unþ1 , pressure field pnþ1=2 and the location of the control points X nþ1 can be summarized as follows: Step 1: Set k :¼ 0 and make an initial guess for X nþ1 , i.e. X ð0Þ as

X ð0Þ ¼ 2X n  X n1 : Step 2: Evaluate the forces over all elements using Eq. (42). Step 3:  Distribute the forces to the nearby Cartesian grid points using (6).  Solve the Navier–Stokes equations using the projection method as described in Section 2.2.  Compute the velocity at the control points, unþ1 ðX ðkÞ Þ by interpolating from the velocity at the surrounding grid points using (14). Step 4:  Evaluate

1 gðX ðkÞ Þ ¼ X ðkÞ  X n  Dtðun ðX n Þ þ unþ1 ðX ðkÞ ÞÞ: 2  If kg ðkÞ k <  then X nþ1 ¼ X ðkÞ and stop the iteration. Otherwise, update X ðkþ1Þ by solving (33) using GMRES. Note that the matrix–vector product required by GMRES is approximated as show in (36). Set k :¼ k þ 1 and go to step 2. We note that solving the Navier–Stokes equations using the projection method involves solving three Helmholtz equations for the intermediate velocity and one Poisson equation for the pressure-increment. These equations are solved using a Fast Fourier Transform algorithm. For the numerical computations presented we employ the subroutines from the FISHPACK library [1,45] especially adapted to our grid arrangement. 4. Numerical results In this section, we present some numerical experiments involving membrane capsules immersed in a three-dimensional fluid domain. The biconcave-shape membranes are governed by the Yeoh form of the strain energy function (23) while the other capsules are modeled as neo-Hookean membranes. The relative tolerance for stopping the GMRES solver is chosen to be 103 and the Newton iteration converges when kg ðkÞ k=kXk < 108 . All the calculations in this section were performed on a Dual-core AMD Opteron 2 GHz Processor. Most simulations were completed within a couple of minutes to hours. 4.1. Oscillatory membrane Our first numerical experiment involves an elastic membrane initially placed in a still fluid. The initial configuration of the membrane is an ellipsoid,

X2 Y 2 Z2 þ þ ¼ 1; a2 b2 c2

ð43Þ

with semi-axes ða; b; cÞ ¼ ð0:25; 0:22; 0:2Þ. The domain of the experiment is the 1 1 1 cube. The domain is discretized with a 64 64 64 uniform grid. The membrane is located at the center of the domain and is discretized into 8192 quadratic triangular elements connecting 16,386 nodes distributed over the membrane. The initial configuration of the membrane is stretched from an undeformed reference configuration which is a spherical membrane of radius 0.2. Because of the

8434

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

incompressibility condition, wepexpect that the membrane should perform damped oscillations and relax to the spherical ffiffiffiffiffiffiffiffi 3 equilibrium state of radius r ¼ abc  0:2224. We start our simulation by setting the velocity and pressure fields to zero, and a homogeneous Dirichlet boundary condition for the velocity is applied at all boundaries. We assume that the surface is a neo-Hookean membrane with the surface elastic modulus Es ¼ 0:1. In this simulation, the fluid density is q ¼ 1, the fluid viscosity is l ¼ 0:01 and the time step is Dt ¼ Dx. With the Jacobian-free Newton–Krylov implicit method, we are able to increase the time step further without loss of stability. Fig. 2(a) and (b) show the number of iterations required at each time step to run the simulations to t ¼ 10 using Dt ¼ Dx and Dt ¼ 4Dx, respectively. It can be seen that we only need 1–4 Newton and GMRES iterations for both time steps. At the first 10 time steps, the Newton and GMRES iterations are zeros because we use a very small time step initially to start the simulation. Fig. 3 shows the configurations of the elastic membrane at different times. The evolution of the three axes of the ellipsoidal membrane is shown in Fig. 4. It can be seen that the membrane oscillates before settling down to the equilibrium state. Fig. 5 shows the time evolution of the three semi-axes of the membrane at l ¼ 0:1. In this case, the Reynolds number is smaller and the membrane relaxes gradually to the equilibrium state without oscillations. For sufficiently small time steps, the explicit and the JFNK implicit method are stable as expected. For larger time steps, we expect that the explicit method will become unstable. For example, when Dt ¼ 4Dx the JFNK method gives good results as 4

4

3

3

2

2

1

1

0

50

100

150

# GMRES iterations

# Newton iterations

# Newton iterations # GMRES iterations

0

timestep

5

5

4

4

3

3

2

2

1

1

0

50

100

150

# GMRES iterations

# Newton iterations

# Newton iterations # GMRES iterations

0

timestep

Fig. 2. Number of Newton iterations and average number of GMRES iterations performed at each time step to run the simulations to t ¼ 10. (a) Dt ¼ Dx and (b) Dt ¼ 4Dx.

8435

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

-0.1

0

0.1

0.1

0

0

-0.1

-0.1

-0.2

-0.2 -0.2

0.2

-0.1

X

-0.2

-0.1

0

Z

0.1

Z

0.2

0

0.1

0.2

X

0.1

0.2

0.1

0.1

0

0

Z

0.2

Z

-0.2

0.2

-0.1

-0.1

-0.2

-0.2 -0.2

0.2

X

-0.1

0

0.1

0.2

X

Fig. 3. Configurations of the membrane at different times.

that obtained at a smaller time step, but the explicit method is unstable. Fig. 6 shows the time evolution of the membrane for both methods before the explicit method goes unstable. The figure shows that the explicit method quickly diverges at step 12 and assumes an erroneous geometry. The JFNK method is stable and the ellipsoidal membrane gradually assumes a spherical configuration. In this example, we study the convergence rate of our scheme in time with l ¼ 0:01; 0:05. Following [21,33], we compute the time discretization error at time T as follows:

eT ðv ; DtÞ ¼ kv ðT; DtÞ  v ðT; Dt=2ÞkL2 : For a vector field wðxÞ ¼ ðw1 ðxÞ; w2 ðxÞ; w3 ðxÞÞ defined on the Cartesian grid, the discrete L2 norm is defined as follows

kwkL2 ¼

X

 3 w21 ðxði; j; kÞÞ þ w22 ðxði; j; kÞÞ þ w23 ðxði; j; kÞÞ h

!1=2 :

i;j;k

Similarly, the discrete L2 norm for a vector field WðXÞ ¼ ðW 1 ðXÞ; W 2 ðXÞ; W 3 ðXÞÞ defined on the interface CðtÞ is defined as follows

kWkL2 ¼

 X 2 W 1 ðX l Þ þ W 22 ðX l Þ þ W 23 ðX l Þ DqDr Ds l

!1=2 :

8436

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

0.25

a b c

Semi-axes

0.24

0.23

0.22

0.21

0.2

0

2

4

6

8

10

Time Fig. 4. The evolution of the membrane axes with

l ¼ 0:01. The membrane oscillates as it converges to the equilibrium state.

0.25

a b c

Semi-axes

0.24

0.23

0.22

0.21

0.2

0

5

10

15

Time Fig. 5. The evolution of the membrane axes with

l ¼ 0:1. The membrane relaxes gradually to the equilibrium state without oscillations.

We calculate the convergence rate in time at T ¼ 1 by varying the time steps in powers of 2 from Dt ¼ 1=16 to Dt ¼ 1=128. The results are shown in Table 1. We can see the second-order convergence in time for both the velocity field and the interface location. 4.2. Membrane capsule in simple shear flow 4.2.1. Spherical membrane In this example, we study the deformation of a spherical neo-Hookean membrane in a shear flow given by the velocity u ¼ ðk z; 0; 0Þ, where k is the shear rate. We compare the results with the linear theory [2] and those obtained by the boundary element method (BEM) [42]. A spherical membrane of radius a is discretized into 5120 quadratic triangular elements connecting 10,242 nodes. This spherical shape is taken to be a strain free state and the shape factors for each element are obtained from this undeformed configuration. The radius of the spherical membrane is chosen to be sufficiently small so that the Reynolds number of the flow is effectively zero; in this work, it is set at Re ¼ 0:001. Hence, the inertia effect is negligible and the results can be compared with those obtained by the linear theory [2] or the BEM [40,42]. The center of the membrane is placed at the center of a cube whose side is 10a and the cube is discretized with 64 64 64 uniform grid. The computational domain is chosen to be large enough so that boundary effects are not important. Boundary conditions for the velocity

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

Fig. 6. Configurations of the membrane at different times with explicit method and implicit method.

8437

8438

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

are of the Dirichlet type at z ¼ 5a and periodic at other boundaries. The deformation of the membrane is described by the Taylor shape parameter Dxz ¼ ðL  BÞ=ðL þ BÞ, where L and B are the maximum and minimum radial distances from the origin in the plane of shear, respectively. We consider the deformation of the membrane for a wide range of dimensionless shear rates,

G¼

lka Es

ð44Þ

;

Table 1 Numerical errors of X and u with different time steps.

l X

u

Dt ¼ 1=16

Dt ¼ 1=32

0.01

1:3600 10

0.005

2:9613 104 4

0.01

4:2159 10

0.005

8:6409 104

Dt ¼ 1=64

Convergence rate

3:3555 10

6

8:3011 10

2.02

6:2315 105

1:4922 105

2.06

4

1:0016 10

5

2:4759 10

2.02

2:0554 104

5:0726 105

2.02

5

4

0.35

0.3 G = 0.05 0.25

Dxz

0.2 G = 0.025 0.15

0.1 G = 0.01 0.05

G = 0.005 G = 0.001

0

0

1

2

3

4

5

kt Fig. 7. The evolution of the Taylor shape parameter Dxz at a sequence of dimensionless shear rates G. The solid lines are the results obtained by the present algorithm, the dashed lines are obtained with the linear theory, and the circles are found using the BEM [42].

0.7 Sui et al. Present Ramanujan et al.

0.6 G = 0.2

0.5

Dxz

0.4 G = 0.1

0.3

0.2

0.1

0

0

5

10

15

20

kt Fig. 8. The evolution of the Taylor shape parameter Dxz for oblate spheroidal capsules with b=a ¼ 0:9 at different dimensionless shear rates G.

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

8439

which expresses the ratio of external viscous stresses to restoring elastic tensions. The time evolution of the Taylor shape parameter with several values of the dimensionless parameter G is shown in Fig. 7. There is excellent agreement with the linear theory at small G ðG < 0:025Þ since deformations are small. Simulation results for very small G ðG ¼ 0:001Þ show that the present implicit immersed boundary method can give good results for stiff problems at relatively large time steps. As the value of G increases the deformations increase and the linear theory does not apply. For larger values of G, we compare our results with those obtained in [42] and Fig. 7 shows good agreement between the two methods for G ¼ 0:05. 4.2.2. Oblate spherical Spherical capsules have been considered initially because of their simple geometry which allows analytical methods to yield deformations histories. Here, we simulate capsules whose initial shapes are oblate spheroids with different aspect ratios. To describe an oblate spheroid with aspect ratio of b=a, we use the mapping xobl ¼ Rx; yobl ¼ Ry; zobl ¼ ðb=aÞRz where ðx; y; zÞ is the coordinate of a point on the unit sphere and the radius R is adjusted to preserve the volume. For an oblate spheroid of volume V, the dimensionless shear rate G is defined in term of the radius of an isovolumic sphere, a ¼ ð3V=4pÞ1=3 . The strain energy function for the neo-Hookean membrane is employed. First, we consider an oblate spheroid of aspect ratio b=a ¼ 0:9, inclined at the angle h0 ¼ p=4 with respect to the streamlines of the unperturbed flow. These parameters were chosen to compare the evolution of the deformation parameter with 0.8 Present Sui et al. Ramanujan et al.

0.7 0.6

Dxz

0.5 0.4 0.3 0.2 0.1 0

0

5

10

15

20

kt Fig. 9. The evolution of the Taylor shape parameter Dxz for oblate spheroidal capsules with b=a ¼ 0:5 at G ¼ 0:2.

0.6

0.5

Dxz

0.4

64x64x64 grid; 4098 nodes 96x96x96 grid; 10242 nodes 128x128x128 grid; 16368 nodes

0.3

0.2

0.1

0

0

5

10

15

20

kt Fig. 10. Grid refinement study for oblate spheroidal capsules with b=a ¼ 0:9 at G ¼ 0:2.

8440

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

those calculate using the BEM [42] and the Lattice Boltzmann method (LBM) [50] under the same conditions. A 64 64 64 fluid grid and membrane discretization of 10,242 nodes and 5120 quadratic triangular elements was employed. The deformation parameter Dxz is calculated for two shear rates G ¼ 0:1 and 0.2 and is shown in Fig. 8 as solid lines. Unlike the spherical capsule, the oblate spheroid capsule undergoes oscillations in the deformation parameter. This has also been observed in [42,50]. Fig. 8 shows good agreement between our results and those obtained in [50]. The BEM [42] predicts smaller deformations than the present method and the LBM [50]. Next, we consider an oblate spheroid capsule with more-oblate unstressed shape of aspect ratio b=a ¼ 0:5, inclined at the angle h0 ¼ p=4 with respect to the streamlines of the unperturbed flow. The deformation parameter Dxz is calculated for G ¼ 0:2 and is shown in Fig. 9 as a solid line. The oblate spheroid capsule undergoes oscillation with larger amplitude in the deformation parameter. Good agreement between our result and that obtained in [50] can be seen in Fig. 9. And again, the BEM [42] predicts smaller deformations than the present method and the LBM [50]. To determine if numerical errors lead to the differences between the present method and the BEM, a grid refinement study was conducted for the oblate spheroid capsule with aspect ratio b=a ¼ 0:9 at G ¼ 0:2. The results are shown in Fig. 10. First a coarser fluid grid, 64 64 64 with a coarser membrane discretization (4098 nodes, 2048 elements) was used. Then, on the 96 96 96 fluid grid, a finer membrane discretization (10,242 nodes, 5120 elements) was employed.

0.3

0.25

Dxz

0.2

0.15

0.1

0.05

0

0

1

2

3

4

5

3

4

5

kt

0.6

0.5

Dxz

0.4

0.3

0.2

0.1

0

0

1

2

kt

^ B ¼ 0, 0.01, 0.025 Fig. 11. The evolution of the Taylor shape parameter Dxz for capsules with spherical initial shapes and reduced modulus of bending. (a) j ^ B ¼ 0, 0.04, 0.1 and 0.15 at G ¼ 0:2. The thickness of the lines increases as the bending modulus is raised. and 0.0375 at G ¼ 0:05. (b) j

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

8441

Finally, the simulation was performed on the 128 128 128 fluid grid with the finest membrane discretization (16,386 nodes, 8192 elements). The results show that the numerical error using the implicit immersed boundary method for the oblate spheroid capsule with aspect ratio b=a ¼ 0:9 at G ¼ 0:2 using these meshes, is small. 4.2.3. Effect of membrane bending stiffness Bending stiffness is believed to play a central role in determining the equilibrium configuration and the shape oscillations of biological membranes. Here, we consider the effect of membrane bending stiffness on the flow-induced deformation of spherical capsules in simple shear flows. We expect that the bending stiffness will restrict the overall capsule deformation. The deformation of the capsule is determined by the dimensionless shear rate G and the reduced ratio of the bending mod^ B jB =ða2 ES Þ. The reference state for computing the bending moment is the membrane with ulus to the elastic modulus, j flat resting shape, corresponding to zero reference curvature. The reference configuration for calculating the in-plane elastic tension is the membrane of an initially spherical capsule of radius a. We note that the reference state of the bending moments is not necessarily the same as that of the elastic tensions, reflecting differences in the physical mechanisms that are responsible for their respective development [40]. Simulations were performed for spherical capsules with different

20 Experiments µ0= 5.3 µN/m µ0= 7.3 µN/m µ0 = 11.3 µ N/m

18 16

Axial diameter

Diameter (µm)

14 12 10 8 6 4 Transverse diameter

2 0

0

50

100

150

200

Stretching Force (pN)

20 Experiments µ0= 3.7 µN/m µ0= 5.3 µN/m µ0 = 9.0 µN/m

18 16

Axial diameter

Diameter (µm)

14 12 10 8 6 4 Transverse diameter

2 0

0

50

100

150

200

Stretching Force (pN)

Fig. 12. Variation of the axial and transverse diameters of the red blood cell against the stretching force. (a) Computational results from [32] using the finite element program ABAQUS. (b) Computational results from [8] using the finite element code ADINA 8.3. The experimental results are taken from [32].

8442

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

shear rates and reduced bending modulus. Fig. 11(a) shows a family of Taylor shape parameters for reduced bending mod^ B ¼ 0, 0.01, 0.025 and 0.0375 at the low shear rate of G ¼ 0:05. Fig. 11(b) illustrates a family of Taylor shape parameters ulus j ^ B ¼ 0, 0.04, 0.1 and 0.15 at the higher shear rate of G ¼ 0:2. We note that the line thickness for reduced bending modulus j increases as the magnitude of the bending modulus increases. The Taylor shape parameters reveal that spherical capsules deform and reach stationary shapes and that raising the bending modulus reduces the capsule deformation.

4.3. Red blood cell deformed by optical tweezers Next, we consider the deformation of human red blood cells subjected to direct stretching by optical tweezers. Numerical simulations were performed to extract the large deformation elastic properties from the experimental results obtained in [9,30,32] during loading as well as upon relaxation of the load. Direct tensile stretching of the human red blood cell using optical tweezers to extract elastic properties was first reported by Henon et al. [19] who attached two silica beads non-specifically to diametrically opposite ends of the cell, trapped both beads with laser beams, and imposed tensile elastic deformation on the cell by moving the trapped beads in opposite directions. In [9,30,32], the optical tweezers system is designed to consist only of a single optical trap, one of the beads is adhered to the glass surface while the other is free to be trapped using the laser beam. By moving one of the beads with the laser beam, the cell is directly stretched. Fig. 12(a) shows the variation of axial and transverse diameters of the red blood cell against the stretching force obtained by both numerical simulations and experimental measurements [32]. In [32], the numerical simulations were performed using the commercially available general purpose finite element package ABAQUS for red blood cells with an initial diameter of 7.82 lm. The shape of the red blood cell is a biconcave shape given as

"



R ZðRÞ ¼ 0:5R0 1  R0

2 #12 "  2  4 # R R ; C0 þ C1 þ C2 R0 R0

ð45Þ

where R2 ¼ X 2 þ Y 2 6 R20 ; R0 ¼ 3:91 lm; C 0 ¼ 0:207161; C 1 ¼ 2:002558, and C 2 ¼ 1:122762. The contact dimension between the beads and the cell was taken to be 2 lm in the simulations. The strain energy function (23) was employed with 0 and C 30 ¼ l =30h  0 , where l is the initial in-plane shear modulus and h 0 is the initial membrane thickness. The C 10 ¼ l =2h 0

0

0

results for the in-plane shear modulus l0 ¼ 5:3—11:3 lN=m are shown in Fig. 12(a). These values are comparable to the range of 4.0–10 lN/m reported in the literature where the estimates have been principally based on micropipette aspiration experiments. Fig. 12(b) illustrates the variation of the axial and transverse diameters of the red blood cell against the stretching force obtained in [8] by numerical simulations. Computational results from [8] were obtained using the finite element code ADINA 8.3 with the same strain energy function (23). The model proposed in [8] consists of a deformable liquid capsule modeled as Newtonian fluid enclosed by a hyperelastic membrane with viscoelastic properties. The membrane shear modulus was estimated to range between 3.7 and 9.0 lN/m. This range also compares well with the reported results obtained from micropipette aspiration experiments.

20 Experiments µ0= 4.8 µN/m µ0= 7.3 µN/m µ0 = 10.0 µN/m

18 16

Axial diameter

Diameter (µm)

14 12 10 8 6 4 Transverse diameter

2 0

0

50

100

150

200

Stretching Force (pN) Fig. 13. Variation of the axial and transverse diameters of the red blood cell against the stretching force. The dotted, solid and dash-dotted lines represent computational results for the axial/transverse diameters with l0 ¼ 4:8, 7.3 and 10.0 lN/m, respectively. The experimental results are taken from [32].

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

8443

In the present paper, we employ the implicit immersed boundary method to study the deformation of human red blood cells subjected to the direct stretching by optical tweezers in order to extract the membrane shear modulus from the experimental results. The tensile forces are applied uniformly on the contact areas of sizes 2 lm between the cell and the beads at opposite ends of the cell. The red blood cell of diameter 7.82 lm is discretized with 16,386 nodes and 8192 quadratic triangular elements. The red blood cell is placed initially at the center of a cube whose side is 8R0 . The computational domain is discretized with 128 128 128 uniform grid. The computational domain is filled with the fluid whose density and viscosity are given as 1050 kg m3 and 4:1 mPa s. These are typical values for the red blood cell cytoplasm density and viscosity. Comparisons of predicted and measured changes in the axial and transverse diameters of the red blood cell using the (Yeoh) strain energy function (23) are plotted in Fig. 13 for l0 ¼ 4:8, 7.3 and 10 lN/m. The simulation is able to capture the experimental trend over the entire range of experimental data well, including the small deformation range and the error bars. The in-plane shear modulus l0 ¼ 4:8—10 lN=m are well comparable to the range of 4.0–10 lN/m obtained from micropipette aspiration experiments. Comparing Figs. 12 and 13, it is apparent that our results are closer to those obtained in [32] than

38.9 pN

Initial undeformed shape

110 pN

194 pN

Fig. 14. The deformed shapes of the red blood cell under the stretching forces.

Fig. 15. Positions of the red blood cells at different times.

8444

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

those in [8]. Fig. 14 illustrates the initial undeformed shape and the deformed shapes of the red blood cell under the stretching force of 38.9, 110 and 194 pN. 4.4. Red blood cells in pressure-driven flow Using the shear modulus obtained from the previous example, we perform a simulation of the pressure-driven flow involving 32 red blood cells (RBCs) immersed in the fluid. The purpose of this simulation is to illustrate the performance of the present implicit immersed boundary method. Specifically, we compare the CPU time that is required for the present implicit scheme and the second-order Runge–Kutta scheme used in [37] to obtain a solution at a given time. In the simulation, the RBCs are initially placed in a box of size 24 lm 48 lm 12 lm as shown in Fig. 15(a). The computational domain is discretized with 64 128 32 uniform grid. The fluid viscosity and density are the same as those in the previous example. The pressure drop is maintained at 10 cm of water column. Each red blood cell is discretized with 10,242 nodes and 5120 quadratic triangular elements. We measure the computational times required for each scheme to obtain the solution at t = 0.05 s with the in-plane shear modulus of 10 lN/m. We used Dt ¼ 4:5 104 for our implicit scheme. We can actually choose larger time step while still maintain the stability, but this choice of Dt is due to the accuracy considerations. The maximum time step that is required for the second-order Runge–Kutta scheme to maintain the stability is Dt ¼ 1:875 106 . All the simulations are performed on a Dual-core AMD Opteron 2 GHz Processor. The positions of the RBCs at t ¼ 0:05 s are shown in Fig. 15(b), applicable for the two different time schemes. We consider periodic boundary conditions for the cells so that the cells which exit at one end of the channel will enter at the other end. We observe that the CPU time for our implicit scheme is about 22 min and this CPU time is about 43 times faster than that of the explicit scheme. The gain over the explicit scheme is even more substantial with larger shear modulus. 5. Conclusions In this paper, we have presented the implicit immersed boundary method for the incompressible Navier–Stokes equations capable of handling three-dimensional membrane–fluid flow interactions. The present study was undertaken with the objective of eliminating time step restrictions by using the Jacobian-free Newton–Krylov method (JFNK) to advance the location of the elastic membrane implicitly. We have shown the capability of the present technique by applying it to the viscous flow problems involving elastic membranes with bending stiffness in three-dimensional domains. Numerical simulations have been performed for the oscillatory membrane in viscous flow to illustrate that the method can maintain stability under relatively large time steps. Simulations have also been performed to reproduce some results for the membrane capsule in simple shear flow as a validation test for our method. It is found that our numerical results are in good agreement with those reported in literature. The present method was also used to simulate the large deformation of human red blood cells subjected to direct stretching by optical tweezers. We have compared the in-plane shear modulus of the normal human red blood cell with that determined by experimental measurements. The range of shear modulus obtained by the present method compares well to that reported previously in the literature. With the obtained shear modulus, the present implicit method has been used to perform simulations of red blood cells in the pressure-driven flow to further demonstrate the performance of the method. The present algorithm does not include any preconditioned techniques with the GMRES solver. We expect to improve the time step further without increasing the number of GMRES iterations by applying appropriate preconditioners. This will be the subject of future work. References [1] J. Adams, P. Swarztrauber, R. Sweet, FISHPACK: Efficient FORTRAN Subprograms for the Solution of Separable Elliptic Partial Differential Equations, 1999. . [2] D. Barthes-Biesel, J.M. Rallison, The time-dependent deformation of a capsule freely suspended in a linear shear flow, J. Fluid Mech. 113 (1981) 251– 267. [3] I. Borazjani, L. Ge, F. Sotiropoulos, Curvilinear immersed boundary method for simulating fluid–structure interaction with complex 3D rigid bodies, J. Comput. Phys. 227 (2008) 7587–7620. [4] D.L. Brown, R. Cortez, M.L. Minion, Accurate projection methods for the incompressible Navier–Stokes equations, J. Comput. Phys. 168 (2001) 464–499. [5] P.N. Brown, Y. Saad, Hybrid Krylov methods for nonlinear systems of equations, SIAM J. Sci. Stat. Comput. 11 (1990) 450–481. [6] X.-C. Cai, D.E. Keyes, V. Venkatakrishnan, Newton–Krylov–Schwarz: an implicit solver for CFD, in: Proceedings of the Eighth International Conference on Domain Decomposition Methods, Wiley, New York, 1997. [7] T.F. Chan, K.R. Jackson, Nonlinearly preconditioned Krylov subspace methods for discrete Newton algorithms, SIAM J. Sci. Stat. Comput. 5 (1984) 533– 542. [8] C.Y. Chee, H.P. Lee, C. Lu, Using 3D fluid–structure interaction model to analyse the biomechanical properties of the erythrocyte, Phys. Lett. A 372 (2008) 1357–1362. [9] M. Dao, C.T. Lim, S. Suresh, Mechanics of the human red blood cell deformed by optical tweezers, J. Mech. Phys. Solid 51 (2003) 2259–2280. [10] R. Dillon, L.J. Fauci, D. Graver, A microscale model of bacterial swimming, chemotaxis and substrate transport, J. Theor. Biol. 177 (1995) 325–340. [11] C.D. Eggleton, A.S. Popel, Large deformation of red blood cell ghosts in a simple shear flow, Phys. Fluid 10 (1998) 1834–1845. [12] L.J. Fauci, A.L. Fogelson, Truncated Newton methods and the modeling of complex immersed elastic structures, Commun. Pure Appl. Math. 66 (1993) 787–818. [13] L.J. Fauci, C.S. Peskin, A computational model of aquatic animal locomotion, J. Comput. Phys. 77 (1988) 85–108. [14] A.L. Fogelson, A mathematical model and numerical method for studying platelet adhesion and aggregation during blood clotting, J. Comput. Phys. 56 (1984) 111–134.

D.V. Le et al. / Journal of Computational Physics 228 (2009) 8427–8445

8445

[15] A.L. Fogelson, Continuum models of platelet aggregation: formulation and mechanical properties, SIAM J. Appl. Math. 52 (1992) 1089–1110. [16] L. Ge, F. Sotiropoulos, A numerical method for solving the 3D unsteady incompressible Navier–Stokes equations in curvilinear domains with complex immersed boundaries, J. Comput. Phys. 225 (2007) 1782–1809. [17] E. Givelberg, Modeling elastic shells immersed in fluid, Commun. Pure Appl. Math. LVII (2004) 0283–0309. [18] Z. Gong, H. Huang, C. Lu, Stability analysis of the immersed boundary method for a two-dimensional membrane with bending rigidity, Commun. Comput. Phys. 3 (3) (2008) 704–723. [19] S. Henon, G. Lenormand, A. Richert, F. Gallet, A new determination of the shear modulus of the human erythrocyte membrane using optical tweezers, Biophys. J. 76 (1999) 1145–1151. [20] T.Y. Hou, J. Lowengrub, M. Shelley, Removing the stiffness from interfacial flows with surface tension, J. Comput. Phys. 114 (1994) 312–338. [21] T.Y. Hou, Z. Shi, An efficient semi-implicit immersed boundary method for the Navier–Stokes equations, J. Comput. Phys. 227 (2008) 8968–8991. [22] T.Y. Hou, Z. Shi, Removing the stiffness of elastic force from the immersed boundary method for the 2D Stokes equations, J. Comput. Phys. 227 (2008) 9138–9169. [23] J. Kim, P. Moin, Application of a fractional step method to incompressible Navier–Stokes equations, J. Comput. Phys. 59 (1985) 308–323. [24] D.A. Knoll, D.E. Keyes, Jacobian-free Newton–Krylov methods: a survey of approaches and applications, J. Comput. Phys. 193 (2004) 357–397. [25] D.A. Knoll, V.A. Mousseau, On Newton–Krylov multigrid methods for the incompressible Navier–Stokes equations, J. Comput. Phys. 163 (2000) 262– 267. [26] D.V. Le, B.C. Khoo, K.M. Lim, An implicit-forcing immersed boundary method for simulating viscous flows in irregular domains, Comput. Method Appl. Mech. Eng. 197 (2008) 2119–2130. [27] D.V. Le, B.C. Khoo, J. Peraire, An immersed interface method for viscous incompressible flows involving rigid and flexible boundaries, J. Comput. Phys. 220 (2006) 109–138. [28] D.V. Le, C. Rosales, B.C. Khoo, J. Peraire, Numerical design of electro-mechanical traps, Lap Chip 8 (2008) 755–763. [29] R.J. LeVeque, Z. Li, Immersed interface method for Stokes flow with elastic boundaries or surface tension, SIAM J. Sci. Comput. 18 (3) (1997) 709–735. [30] C.T. Lim, M. Dao, S. Suresh, C.H. Sow, K.T. Chew, Large deformation of living cells using laser traps, Acta Mater. 52 (2004) 1837–1845. [31] A. Mayo, C.S. Peskin, An implicit numerical method for fluid dynamics problems with immersed elastic boundaries, Contemp. Math. 141 (1993) 261– 277. [32] J.P. Mills, L. Qie, M. Dao, C.T. Lim, S. Suresh, Nonlinear elastic and viscoelastic deformation of the human red blood cell with optical tweezers, Mech. Chem. Biosyst. 1 (3) (2004) 169–180. [33] Y. Mori, C.S. Peskin, Implicit second-order immersed boundary methods with boundary mass, Comput. Method Appl. Mech. Eng. 197 (2008) 2049– 2067. [34] E.P. Newren, Unconditionally stable discretizations of the immersed boundary equations, J. Comput. Phys. 222 (2007) 702–719. [35] E.P. Newren, A.L. Fogelson, R.D. Guy, R.M. Kirby, A comparison of implicit solvers for the immersed boundary equations, Comput. Method Appl. Mech. Eng. 197 (2008) 2290–2304. [36] C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (1977) 220–252. [37] C.S. Peskin, The immersed boundary method, Acta Numer. 11 (2) (2002) 479–517. [38] C. Pozrikidis, Finite deformation of liquid capsules enclosed by elastic membranes in simple shear flow, J. Fluid Mech. 297 (1995) 123–152. [39] C. Pozrikidis, Numerical Computation in Science and Engineering, Oxford University Press, 1998. [40] C. Pozrikidis, Effect of membrane bending stiffness on the deformation of capsules in simple shear flow, J. Fluid Mech. 440 (2001) 269–291. [41] C. Pozrikidis, Numerical simulation of the flow-induced deformation of red blood cells, Ann. Biomed. Eng. 31 (2003) 1194–1205. [42] S. Ramanujan, C. Pozrikidis, Deformation of liquid capsules enclosed by elastic membrane in simple shear flow: large deformations and the effect of fluid viscosities, J. Fluid Mech. 361 (1998) 117–143. [43] A.M. Roma, C.S. Peskin, M.J. Berger, An adaptive version of the immersed boundary method, J. Comput. Phys. 153 (1999) 509–534. [44] Y. Sadd, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (1986) 856–869. [45] U. Schumann, R.A. Sweet, A direct method for the solution of Poisson’s equation with Neumann boundary conditions on a staggered grid of arbitrary size, J. Comput. Phys. 20 (1976) 171–182. [46] R. Skalak, A. Tozeren, R.P. Zarda, S. Chien, Strain energy function of red blood cell membranes, Biophys. J. 13 (1973) 245–264. [47] J.M. Stockie, B.R. Wetton, Stability analysis for the immersed fiber problem, SIAM J. Appl. Math. 55 (1995) 1577–1591. [48] J.M. Stockie, B.R. Wetton, Analysis of stiffness in the immersed boundary method and implications for time-stepping schemes, J. Comput. Phys. 154 (1999) 41–64. [49] J. Stoer, R. Bulirsch, Introduction to Numerical Analysis, third ed., Springer-Verlag, 2002. [50] Y. Sui, Y.T. Chew, P. Roy, H.T. Low, A hybrid method to study flow-induced deformation of three-dimensional capsules, J. Comput. Phys. 227 (2008) 6351–6371. [51] C. Tu, C.S. Peskin, Stability and instability in the computation of flows with moving immersed boundaries: a comparison of three methods, SIAM J. Sci. Statist. Comput. 13 (1992) 1361–1376. [52] N.T. Wang, A.L. Fogelson, Computational methods for continuum models of platelet aggregation, J. Comput. Phys. 151 (1999) 649–675. [53] X.S. Wang, An iterative matrix-free method in implicit immersed boundary/continuum methods, Comput. Struct. 85 (2007) 739–748. [54] O.H. Yeoh, Characterization of elastic properties of carbon-black-filled rubber vulcanizates, Rubber Chem. Technol. 63 (1990) 792–805.