A FRACTIONAL STEP IMMERSED BOUNDARY METHOD FOR STOKES FLOW WITH AN INEXTENSIBLE INTERFACE ENCLOSING A SOLID PARTICLE

(1)

SIAM J. SCI.COMPUT. c 2012 Society for Industrial and Applied Mathematics Vol. 34, No. 5, pp. B692–B710

A FRACTIONAL STEP IMMERSED BOUNDARY METHOD FOR STOKES FLOW WITH AN INEXTENSIBLE INTERFACE

ENCLOSING A SOLID PARTICLE∗

MING-CHIH LAI†, WEI-FAN HU‡, AND WEN-WEI LIN‡

Abstract. In this paper, we develop a fractional step method based on the immersed boundary

(IB) formulation for Stokes flow with an inextensible (incompressible) interface enclosing a solid particle. In addition to solving for the fluid variables such as the velocity and pressure, the present problem involves finding an extra unknown elastic tension such that the surface divergence of the velocity is zero along the interface, and an extra unknown particle surface force such that the velocity satisfies the no-slip boundary condition along the particle surface. While the interface moves with local fluid velocity, the enclosed particle hereby undergoes a rigid body motion, and the system is closed by the force-free and torque-free conditions along the particle surface. The equations are then discretized by standard centered difference schemes on a staggered grid, and the interactions between the interface and particle with the fluid are discretized using a discrete delta function as in the IB method. The resultant linear system of equations is symmetric and can be solved by fractional steps so that only fast Poisson solvers are involved. The present method can be extended to Navier– Stokes flow with moderate Reynolds number by treating the nonlinear advection terms explicitly for the time integration. The convergent tests for a Stokes solver with or without an inextensible interface are performed and confirm the desired accuracy. The tank-treading to tumbling motion for an inextensible interface enclosing a solid particle with different filling fractions under a simple shear flow has been studied extensively, and the results here are in good agreement with those obtained in literature.

Key words. immersed boundary method, inextensible interface, solid particle, fractional step

method, Stokes ﬂow

AMS subject classifications. 65M06, 76D07 DOI. 10.1137/100818777

1. Introduction. Most biological cells consist of a lipid bilayer membrane en-capsulating the cellular content. For instance, a human leukocyte contain a nucleus occupying 18%–44% of the volume [21], which affects the cell adhesion during inflam-matory responses. It is important to study the dynamics of a compound vesicle (a lipid bilayer membrane enclosing a fluid with a suspended particle) [26]. For simplic-ity, we regard the lipid bilayer membrane as an inextensible interface. Thus, in this paper, we develop a fractional step method based on the immersed boundary (IB) formulation for Stokes flow with an inextensible (incompressible) interface enclosing a solid particle.

Most previous studies in the literature have focused on an inextensible interface enclosing a homogeneous fluid. Even in that case, there are still some numerical issues to be addressed. In additional to solving for the fluid variables such as the velocity and pressure, the above problem involves finding extra unknown elastic tension such that ∗_{Submitted to the journal’s Computational Methods in Science and Engineering section December}

20, 2010; accepted for publication (in revised form) August 22, 2012; published electronically October 23, 2012.

http://www.siam.org/journals/sisc/34-5/81877.html

†_{Corresponding author. Center of Mathematical Modeling and Scientiﬁc Computing &}

Depart-ment of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan (mclai@math.nctu.edu.tw). This author’s work was supported in part by National Science Council of Taiwan under research grant NSC-98-2115-M-009-014-MY3, and by the NCTS in Taiwan.

‡_{Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road,}

Hsinchu 300, Taiwan (weifanhu.am95g@g2.nctu.edu.tw, wwlin@math.nctu.edu.tw). B692

(2)

AN INEXTENSIBLE INTERFACE WITH A SOLID PARTICLE B693

the surface divergence of the velocity is zero along the interface. Once the velocity is found, the interface moves according to the local ﬂuid velocity as usual. Since the interface is inextensible, the total length of the interface should be conserved. This mathematical model is motivated by the simulation of vesicle dynamics [10], and the deformation of erythrocytes [14, 17] and drug-carrying capsules [20], just to name a few. In particular, the dynamics of moving vesicles have been extensively explored both experimentally [6, 9] (see also and the references therein) and computationally [10, 28, 25, 22, 7, 11]. Notice that the dynamics of vesicles are determined by their boundary rigidity, inextensibility, and the hydrodynamical forces.

In previous literature, most related work is based on boundary integral methods; see, for example, [28, 25, 22] and the references therein. However, boundary integral methods generally assume infinite domains and cannot be generalized to full Navier– Stokes equations since there is no corresponding Green function. Until recently, Kim and Lai [7] applied a penalty IB method to simulate the dynamics of inextensible vesicles. By introducing two different kinds of Lagrangian markers, the authors are able to decouple the fluid and vesicle dynamics so that the computation can be per-formed more efficiently. One potential problem with this approach is that the time step depends on the penalty numbers and must be chosen smaller as the penalty number becomes larger. In [11], a new finite difference scheme based on the immersed interface method (IIM) has been developed for solving the present problem in Navier– Stokes flow. The authors treat the unknown elastic tension as an augmented variable so that the augmented IIM can be applied.

We summarize the contributions of the present work as follows:

• We develop a linearly semi-implicit scheme (section 3) for the model of Stokes

ﬂow with an inextensible interface enclosing a solid particle, which has poten-tial applications to the dynamics of a compound vesicle and is less investigated in the literature. The present scheme can be extended to Navier–Stokes ﬂow with moderate Reynolds number by treating the nonlinear advection terms explicitly.

• We show that the spreading operator of the tension and the surface

diver-gence operator of the velocity are skew-adjoint mathematically (section 2). This skew-adjoint property is also preserved in the discrete form (section 4.1), which makes the resultant linear system symmetric (section 4.1 and also Appendix B). A fractional step method that exploits fast Poisson solvers can be eﬃciently applied to solving the linear system (section 4.3).

• Unlike our previous work using a penalty approach [7], here we are able to

estimate the local error of inextensibility for two successive time steps (see (21) in section 3). In addition, there are no penalty parameters introduced as in [7], so the time step size can be significantly increased (see the time step chosen in section 5.2). From those points of view, the present scheme does offer significant improvements or the accuracy and stability achieved in the previous work.

• We present numerical results (section 5) of the tank-treading to tumbling

motions for an inextensible interface enclosing a solid particle with different filling fractions under shear flow. We have found that, by the increase of the filling fraction, the interface tends to transit from tank-treading to tumbling (section 5.3). We also compute the critical filling fractions for different re-duced areas (section 5.3), and the results are qualitatively consistent with those in [26].

The rest of the paper is organized as follows. In the next section, we describe

(3)

B694 MING-CHIH LAI, WEI-FAN HU, AND WEN-WEI LIN

our governing equations for the Stokes flow with an inextensible interface enclosing a suspended solid particle based on an IB formulation. We also show the skew-adjoint property between the spreading operator acting on the tension and the surface divergence operator acting on the velocity. In section 3, we discretize the governing equations by centered difference schemes on a staggered grid. Some implementation details including the proof of the symmetric property of the resultant matrix and the existence of a solution in the linear system, and the fractional step method for solving the resultant linear system are all given in section 4. Numerical results including the convergence and efficiency tests and the simulations for tank-treading to tumbling motion with different filling fractions for the present model are shown in section 5. Some conclusions are given in section 6.

2. Governing equations. Consider a moving, immersed, inextensible interface Γ enclosing a suspended solid particle P in a two-dimensional fluid domain Ω, as shown in Figure 2.1. We assume that the fluids inside and outside of the interface are the same and governed by the incompressible Stokes equations, and the particle’s gravity is neglected. Using the IB formulation of the model [16, 7], the inextensible interface and the solid particle surface are regarded as singular force generators in the fluid equations. Let us describe the inextensible interface Γ and the particle surface ∂P by parametric forms X(s, t) = (X₁(s, t), X₂(s, t)) and Y(α, t) = (Y₁(α, t), Y₂(α, t)), respectively, where s and α are the corresponding Lagrangian parameters. The gov-erning equations in dimensionless form can be written as follows:

(2.1) −∇p + Δu + Γ ∂ ∂s(στ )δ(x − X(s, t)) ds + ∂P F(α, t)δ(x− Y(α, t)) dα = 0 in Ω, (2.2) ∇ · u = 0 in Ω, (2.3) ∇_s· U = ∂U ∂s · τ |∂X/∂s| = 0 on Γ, (2.4) U(s, t) = Ω u(x, t)δ(x− X(s, t)) dx, ∂X ∂t(s, t) = U(s, t) on Γ, (2.5) V(α, t) = Ω u(x, t)δ(x−Y(α, t)) dx = V_c(t)+ω(t) −(Y2(α, t)− Y2c(t)) Y₁(α, t)− Y1c(t) on ∂P, Γ P Ω

Fig. 2.1_{. A diagram of an inextensible interface enclosing a solid particle in a shear ﬂow.}

(4)

AN INEXTENSIBLE INTERFACE WITH A SOLID PARTICLE B695 (2.6) ∂P F(α, t) dα = 0, where F(α, t) = (F₁(α, t), F₂(α, t)), (2.7) ∂P F₁(α, t) (Y₂(α, t)− Y_2c(t))− F₂(α, t) (Y₁(α, t)− Y_1c(t)) dα = 0. Equations (2.1) and (2.2) are the familiar incompressible Stokes equations with singular force terms arising from the inextensible interface and the particle surface. Here, u = (u, v) and p are the velocity ﬁeld and the pressure, both described in Cartesian coordinates. Equation (2.3) represents the inextensibility constraint of the interface which is equivalent to the zero surface divergence of the velocity along the interface due to the fact that the local stretching factor [13],

∂ ∂t ∂X ∂s  = (∇s· U) ∂X ∂s  on Γ,

must be zero along the interface. Equation (2.4) simply states that the interface moves along with the fluid and that the velocity U is the interpolation of the fluid velocity at the interface. The interaction between the fluid and the interface or the particle is linked by the two-dimensional Dirac delta function δ(x) = δ(x)δ(y).

Since the particle behaves like a rigid body, its motion is governed by the trans-lational and rotational velocities around the center of the particle. Moreover, the velocity of the particle surface must equal the ﬂuid velocity (no-slip boundary con-dition). Equation (2.5) describes the particle surface velocity in which V_c(t) is the translational velocity of the center of particle and ω(t) is the angular velocity compo-nent of the particle. The particle centerY_c(t) = (Y_1c(t), Y_2c(t)) thus moves according to dYc

dt =Vc(t). Since the particle center velocity Vc(t) and the angular velocity

component ω(t) are also unknown, the system of equations must be closed by the force-free (2.6) and free (2.7) conditions for a rigid body motion. The torque-free condition is derived in detail in Appendix A.

Like the pressure in incompressible flow, the elastic tension σ and the particle surface force F in the present IB formulation are not known a priori, and they must be determined as parts of the solution. In fact, the tension and the particle surface force play the role of Lagrange multipliers in enforcing the local inextensible constraint (2.3) along the interface and the no-slip boundary condition (2.5) along the particle surface, respectively. So the major difficulty of solving the above interfacial problem enclosing a solid particle is that one needs to find those solution variables simultaneously. In this paper, we discretize (2.1)–(2.7) directly and use a fractional step method to solve the resultant linear system of equations. The detailed numerical algorithm will be given in the next section.

In the IB formulation, the force spreading operator and the velocity interpolating operator are self-adjoint in both continuous and discrete senses [16]. Here, we would like to show that the spreading operator acting on the function σ and the surface divergence operator of the velocity are skew-adjoint with each other. To proceed, let us deﬁne the spreading operator S of σ and the surface divergence operator∇_sof U as follows: (2.8) S(σ) = Γ ∂ ∂s(στ )δ(x − X(s)) ds, (2.9) ∇_s· U = ∂U ∂s · τ |∂X/∂s| = ∂ ∂s Ωu(x)δ(x− X(s)) dx · τ |∂X/∂s|.

(5)

We also deﬁne the inner product of functions on Ω and Γ in the following: (2.10) u, v_Ω= Ω u(x)· v(x) dx, f, g_Γ= Γ f (l) g(l) dl,

where l in (2.10) is the arc-length parameter. Then we have

S(σ), uΩ= Ω Γ ∂ ∂s(στ )δ(x − X(s)) ds · u(x) dx = Γ ∂ ∂s(στ ) · Ω u(x)δ(x− X(s)) dx ds =− Γ σ τ ·∂U ∂s ds

(integration by parts and the closed interface) = Γ σ −∂U ∂s · τ |∂X/∂s| ∂X ∂s  ds =σ, −∇_s· UΓ=σ, S∗(U)Γ. (2.11)

It follows that the spreading operator and the surface divergence operator are skew-adjoint.

The reason for showing the skew-adjointness of those two operators is twofold. First, since the surface divergence of the velocity is zero, the above derivation leads toS(σ), uΩ= 0. That is, the present elastic tension does not do work to the fluid, which is not surprising since it is merely the Lagrange multiplier for the inextensible constraint. However, if we add bending force along the interface, such as the one in vesicle problems [25, 7], then the bending force does do work to the fluid. Second, the skew-adjointness is also satisfied in the discrete sense (see the next section) so that the resultant matrix is symmetric.

3. Numerical discretization. We now are ready to discretize (2.1)–(2.7) by the IB method. For simplicity, we assume that the computational domain Ω = [a, b]×[c, d] is a rectangular domain and that the fluid variables are defined on the staggered marker-and-cell (MAC) grid [4]. That is, as shown in Figure 3.1, the pressure is defined on the grid points labelled as x = (x_i, y_j) = (a + (i−1/2)Δx, c+(j −1/2)Δy),

i = 1, 2, . . . , m and j = 1, 2, . . . , n, while the velocity components u and v are deﬁned

at (x_i−1/2, y_j) = (a+(i−1)Δx, c+(j−1/2)Δy) and (x_i, y_j−1/2) = (a+(i−1/2)Δx, c+ (j− 1)Δy), respectively. Here, we assume that a uniform mesh width h = Δx = Δy is used, although that is not necessary. For the immersed inextensible interface X, we use a collection of discrete points s_k = kΔs, k = 0, 1, . . . M , with the interface mesh width Δs so the Lagrangian markers of the interface are represented by X_k = X(s_k). Similarly, we use a collection of discrete points α_k = kΔα, k = 0, 1, . . . M_p, so that the particle surface points are represented by Y_k = Y(α_k). Both Δs and Δα are roughly chosen as a half of fluid mesh h. Both the interface and the particle surface are closed so that we have X₀= X_M and Y₀= Y_M_p. The elastic tension is defined at the “half-integer” points given by s_k−1/2 = (k− 1/2)Δs, so we denote it by σ_k−1/2. Without loss of generality, for any function defined on the interface ψ(s), we approximate the partial derivative ∂ψ

∂s by the centered diﬀerence scheme as

(3.1) D_sψ = ψ(s + Δs/2)− ψ(s − Δs/2)

Δs .

(6)

Fig. 3.1. The computational domain Ω with staggered grid.

Thus, the interface stretching factor and the unit tangent can be approximated by

|DsX| and τ = DsX/|DsX|, which in turn are also deﬁned at the half-integer points.

We denote them by|D_sX|_k−1/2 andτ_k−1/2, respectively.

Let Δt be the time step size, and the superscript of the variables denote the time step index. At the beginning of each time step l, the interface conﬁguration Xl_k, the particle surface Yl

k = (Y1kl, Y2kl ), and its center Ycl = (Y1cl, Y2cl) are all given. Despite

the fact that the problem is nonlinear, here we propose a linearly semi-implicit scheme for (2.1)–(2.7). The time step can be advanced as follows:

(3.2) −∇hpl+1+ Δhul+1+ M k=1 D_s σl+1τl kδh(x− X l k)Δs + Mp k=1 Fl+1_k δ_h(x− Y_kl)Δα = 0, (3.3) ∇_h· ul+1= 0, (3.4) ∇_s_h· Ul+1_k = U l+1 k − Ul+1k−1 Δs · τ l k−1/2DsXl_k−1/2= 0 for k = 1, 2, . . . , M, (3.5) Ul+1_k = x u(x)l+1δ_h(x− Xl_k) h2 for k = 1, 2, . . . , M,

(7)

B698 MING-CHIH LAI, WEI-FAN HU, AND WEN-WEI LIN Vl+1 k = x u(x)l+1δ_h(x− Yl_k) h2 =Vl+1_c + ωl+1 −(Yl 2k− Y2cl) Yl 1k− Y1cl for k = 1, 2, . . . , M_p, (3.6) (3.7) Mp k=1 Fl+1_k Δα = 0, (3.8) Mp k=1 F_1kl+1(Y_2kl − Y_2cl)− F_2kl+1(Y_1kl − Y_1cl)Δα = 0,

where the spatial operators ∇_h, Δ_h, and ∇_s_h· are the standard second-order cen-tered diﬀerence approximations to the gradient, Laplacian, and surface divergence. Here, δ_h is a smoother version of discrete delta function developed in [27] as δ_h(x) =

1 h2φ(x_h)φ(y_h), with φ(r) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 3 8+32π −r 2 4, |r| < 0.5, 1 4+1−|r|8 −2 + 8|r| − 4r2₋1 8sin−1 √ 2(|r| − 1), 0.5≤ |r| ≤ 1.5, 17 16 −64π −3|r|4 +r 2 8 +|r|−216 −14 + 16|r| − 4r2 +₁₆1 sin−1 √2(|r| − 2), 1.5≤ |r| ≤ 2.5, 0, 2.5≤ |r|.

One should notice that this discrete delta function has one grid support wider than the one often used in the community. The advantage of using the above smoothing delta function is that it reduces oscillations of the elastic tension caused by the present IB method. The computation for the elastic tension cause oscillations in other literature [28, 7, 11].

Once we obtain the new velocity ﬁeld ul+1 _{on the ﬂuid grid, we can interpolate}

the new velocity to the marker points by (3.5) and move the Lagrangian markers to new positions. That is,

(3.9) Xl+1_k = Xl_k+ Δt Ul+1_k . Therefore, we have Xl+1_k − Xl+1_k−1 Δs = Xl k− Xlk−1 Δs + Δt Ul+1_k − Ul+1_k−1 Δs .

By multiplying the above equation by itself and using the zero discrete surface diver-gence (3.4), we obtain the equality

Xl+1_k − Xl+1_k−1 Δs 2 = Xl k− Xlk−1 Δs 2 + (Δt)2 Ul+1_k − Ul+1_k−1 Δs 2 , which leads to (3.10) |D_sXl+1|2_k−1/2=|D_sXl|2_k−1/2+ (Δt)2|D_sUl+1|2_k−1/2.

(8)

Thus, we conclude that the pointwise error for the local stretching factor is ﬁrst-order accurate, which is comparable to the accuracy of the IB method.

For the solid particle motion, the particle center velocity V_c and the angular velocity component ω will be solved as parts of the solution. Instead of ﬁrst using the interpolation formula to ﬁnd the new velocity and then updating the particle surface points, we simply adopt the idea of rigid body motion to determine the particle surface position. That is, we compute the new particle center and the rotational angle by (3.11) Yl+1_c = Yl_c+ Δt V_cl+1, θl+1= θl+ Δt ωl+1,

and then use them to ﬁnd the new particle surface position as

(3.12) Y_kl+1= Y_cl+1+ cos θl+1 _{− sin θ}l+1 sin θl+1 _{cos θ}l+1 Y_k0. 4. Implementation details.

4.1. Discrete skew-adjoint operators. In this subsection, we show that the spreading operator S acting on the elastic tension and the surface divergence operator

∇sacting on the velocity are also skew-adjoint in the discrete sense. That is, we shall

prove the numerical identity for (2.11). To proceed, we first define the corresponding discrete inner product on the fluid grid Ω_h and the interfacial grid Γ_h as

(4.1) u, vΩh = x u(x)· v(x) h2, φ, ψΓ_h = M k=1 φ_k−1/2ψ_k−1/2|D_sX|_k−1/2Δs,

where the second summation is nothing but the midpoint rule for the second integral of (2.10). We also deﬁne the discrete spreading operator S_h acting on the discrete elastic tension σ_h as (4.2) S_h(σ_h) = M k=1 D_s(στ )_kδ_h(x− X_k)Δs. Then we have Sh(σh), uΩh = x _M k=1 D_s(στ )_kδ_h(x− X_k)Δs · u(x) h2 = M k=1 D_s(στ )_k· x u(x)δ_h(x− X_k) h2 Δs = M k=1 D_s(στ )_k· U_kΔs = M k=1 σ_k+1/2τ_k+1/2− σ_k−1/2τ_k−1/2 Δs · UkΔs =− M k=1 σ_k−1/2 U_k− U_k−1 Δs · τk−1/2Δs (summation by parts) =σ_h,−∇_s_h· UΓ_h =σ_h, S_h∗(U)Γ_h.

Note that this discrete skew-adjoint property is crucial to our IB formulation for solving (3.2)–(3.8). Due to the fact that discrete surface divergence of the velocity is

(9)

zero in (3.4), we can rescale this constraint to make the resultant matrix obtained by (3.4) be the transpose of the matrix obtained by the discrete spreading operator of the tension. One can also verify this symmetric property by expressing those terms explicitly, which is given in Appendix B.

4.2. Existence of a solution. By using the staggered grid to discretize the ﬂuid variables, the matrix obtained by the discrete divergence of the ﬂuid velocity can be written as the transpose of the discrete gradient of the pressure. As derived before, the matrix obtained by the discrete surface divergence of the velocity can be written as the transpose of the matrix obtained by the discrete spreading operator of the tension. Similarly, the matrix obtained by the discrete spreading operator of the force arising from the particle solid boundary can be written as the transpose of the matrix obtained by the discrete interpolating operator of velocity. Thus, the linear system for (3.2)–(3.8) is symmetric and can be written as

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ L G S E 0 GT ₀ ₀ ₀ ₀ ST ₀ ₀ ₀ ₀ ET ₀ ₀ ₀ _R 0 0 0 RT 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ u p σ F Θ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ bc₁ bc₂ 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , (4.3)

where the submatrix L, G, S, E, and R represents the discrete Laplacian Δ_h, discrete gradient−∇_h, the discrete spreading operator S_h on tension, the discrete spreading operator E_h on the surface force, and the discrete rigid body motion equation. The submatrix sizes of L, G, S, E, and R are ((m−1)n+m(n−1))×((m−1)n+m(n−1)), ((m−1)n+m(n−1))×(mn), (m−1)n+m(n−1)×M, (m−1)n+m(n−1)×2M_p, and 2M_p× 3, respectively. The unknown Θ = [V_1c, V_2c, ω]T_{; the right-hand-side vector}

[bc₁, bc₂, 0, 0, 0]T _{of (4.3) consists only of the velocity boundary conditions, since the}

pressure does not need the boundary condition in the present grid arrangement. The detailed entities of the submatrices S and R can be found in Appendices A and B.

Let us discuss the existence of the solution for the linear system (4.3). For clarity, we denote the matrix in (4.3) by A. As known, without the eﬀects of the inextensible interface and the solid particle, the problem becomes pure Stokes ﬂow, and the linear system is (4.4) L G GT ₀ u p = bc₁ bc₂ .

Let us denote the matrix in (4.4) by A. It is wellknown that the nullity of A equals 1

since the pressure is unique up to a constant, and the existence of a solution can be veriﬁed by using the discrete incompressible constraint (3.3). To be precise, since the rank of deﬁciency of A is 1, based on the algebraic structure of the submatrix G, the

kernel of A is (4.5) ker( A) = span{[ 0· · · 0 (m−1)n+m(n−1) 1· · · 1 mn ]T}.

(10)

And for any vector z∈ ker( AT_{) = ker(}_{A), we have}

zT bc₁ bc₂ = mn k=1 (bc₂)_k= n j=1 u_0,j h − n j=1 u_m,j h + m i=1 v_i,0 h − m i=1 v_i,n h =− ⎛ ⎝−n j=1 u_0,jh + n j=1 u_m,jh− m i=1 v_i,0h + m i=1 v_i,nh ⎞ ⎠ h−2 = ⎛ ⎝m i=1 n j=1 u_i,j− u_i−1,j h + v_i,j− v_i,j−1 h h2 ⎞ ⎠ h−2

= 0 (by the discrete incompressibility (3.3)), (4.6)

which shows the compatibility condition for the existence of a solution.

If we add the eﬀects of the interface and the solid particle, then the matrix A is

augmented by additional submatrices to become A as in (4.3). Since the matrices S and E come from the discrete spreading operator, the entries of S and E depend on the positions of moving Lagrangian markers. It is unlikely to show rigorously that the nullity of A is exactly equal to 1; however, we have conﬁrmed the above statement to be true in our numerical experiments. So the apparent kernel can be

(4.7) ker(A) = span{[ 0· · · 0 (m−1)n+m(n−1) 1· · · 1 mn 0· · · 0 M 0· · · 0 2Mp 0, 0, 0 3 ]T},

and the existence of a solution for the linear system (4.3) follows the equality of (4.6) immediately.

4.3. Fractional step method. In this subsection, we apply the idea of the fractional step method developed by Taira and Colonius [24] to solve the resultant linear system of equations (4.3). In [24], the authors applied the IB method to simulate the incompressible ﬂow over solid bodies with prescribed body surface velocity. Unlike the previous feedback forcing approaches in [5, 19, 12], they introduced the boundary force as another Lagrange multiplier to enforce the no-slip condition for the velocity at the immersed boundary. From this point of view, the present approach shares the spirit of [24] by introducing the elastic tension as a new Lagrange multiplier to enforce the surface divergence-free constraint (3.4) along the interface. Thus, it is natural to group all Lagrange multipliers into a new column vector q = [p, σ, F, Θ]T and combine the submatrices G, S, and E as Q = [G, S, E, 0], so the linear system (4.3) becomes (4.8) L Q QT _D u q = bc₁ b , where D = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 0 0 0 0 0 0 0 0 R 0 0 RT ₀ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ and b = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ bc₂ 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .

As in [24], we perform a block LU decomposition to the above equation to obtain L 0 QT −QTL−1Q + D I L−1Q 0 I u q = bc₁ b . (4.9)

(11)

Note that the above decomposition is possible (L−1exists) since the matrix L arising from the discrete Laplacian operator is symmetric and negative deﬁnite. One can further split the above matrix equation into the following steps by introducing an intermediate velocity u∗: Lu∗= bc₁, (4.10) −QT_L−1_{Q + D}_{q =}_b_{− Q}T_u∗_, (4.11) u = u∗− L−1Qq. (4.12)

Now we are ready to describe the detailed numerical implementation for solving (4.10)–(4.12). It is common practice to avoid the direct computation of the inverse of the matrix L since it is too expensive. In [24], a second-order approximation for

L−1 based on Taylor expansion is implemented for solving equations similar to our (4.11)–(4.12), and the conjugate gradient method is applied to solve those equations iteratively. However, this leads to another time step constraint related to the viscosity and the eigenvalues of the discrete Laplacian. In this paper, since we are working on the Stokes ﬂow rather than the Navier–Stokes, we are unable to approximate L−1 using Taylor’s expansion. Although we do not approximate the L−1 directly, we still can solve (4.10)–(4.12) eﬃciently thanks to the fast Poisson solver developed in public software package FISHPACK [1]. (The present matrix L is nothing but the discrete Laplacian operator.) The detailed steps for solving (4.10)–(4.12) are as follows:

Step 1. Solve (4.10) by two fast Poisson solvers to obtain intermediate velocity ﬁeld

u∗.

Step 2. Solve (4.11) by the GMRES method. In each iteration, a matrix-vector

prod-uct (−QT_L−1_{Q) ϕ is needed; fortunately, this can be obtained by letting} z = L−1Qϕ and solving Lz = Qϕ. Once it is done, we multiply z by−QT _to

obtain the product needed. Again, the time-consuming cost in each iteration is one fast Poisson solver.

Step 3. Find the velocity ﬁeld u from (4.12). Since q is solved via Step 2, by solving Lw = Qq, we then obtain u = u∗− w. Again, this involves applying two fast

Poisson solvers.

Therefore, the overall cost in Steps 1–3 for our present numerical algorithm can be counted in terms of the number of fast Poisson solvers applied. In the next section, we shall show the numbers of fast Poisson solvers used in the Stokes ﬂow for diﬀerent grid resolutions.

5. Numerical results. In this section, we perform a series of numerical tests for the present scheme. Throughout this section, the computational domain is chosen as Ω = [−1, 1] × [−1, 1]. All numerical runs were carried out on a PC with 4G of RAM using double precision arithmetics.

5.1. Convergence and efficiency tests for the Stokes solver. In this sub-section, we perform the convergence test and evaluate the efficiency for the present Stokes solver. The numerical algorithm for solving this problem is exactly same as Steps 1–3 described in the previous section but in a simpler manner; that is, Q = G and q = p. Other efficient Stokes solvers can be found in [23, 2, 18] and the references therein.

Here, we use the following test example so that we can easily compute the errors between the exact and numerical solutions:

u_e(x, y) = sin x cos y, v_e(x, y) =− cos x sin y, and p_e(x, y) = exsin y.

(12)

Table 5.1

Numerical accuracy for the Stokes solver.

m = n 32 64 128 256 512

ue− uh∞ 1.578e-4 4.481e-5 1.206e-5 3.153e-6 8.120e-7

Rate - 1.82 1.89 1.94 1.96

ve− vh∞ 1.578e-4 4.481e-5 1.206e-5 3.153e-6 8.120e-7

Rate - 1.82 1.89 1.94 1.96

pe− ph∞ 9.615e-4 4.286e-4 2.052e-4 1.005e-4 4.970e-5

Rate - 1.17 1.06 1.03 1.02

Table 5.2

The cost in CPU time and iterations.

m = n 32 64 128 256 512

Iterations 12 14 15 16 18

CPU time(sec) 0.02 0.05 0.20 0.93 4.97

Note that the above solution does not satisfy the pure Stokes equations, so we need to add some external force ﬁeld (can be easily computed) into the equations. However, it does not change the method or algorithm since the extra force term appears on the right-hand side of equations. Along the boundary of computational domain, the Dirichlet boundary conditions for the velocity are provided, while no pressure boundary condition is needed in our setting.

It is worth mentioning that the pressure is unique up to a constant in Stokes equations. Rather than pinning a certain value to a particular discrete pressure as in [24], uniqueness can be guaranteed by setting up a constraint for the discrete pressure as (5.1) i,j p_i,jh2= Ω p_e(x) dx.

So our initial guess p0_ijin the GMRES iteration can be chosen as p0_ij=#_Ωp_e(x) dx/|Ω|. In those tests, the tolerance of the residual is chosen as 10−8.

Table 5.1 shows the maximum errors between the exact and numerical solutions for different grid resolutions. One can see that the velocity field has clear second-order accuracy, while the pressure has first-order accuracy. Table 5.2 shows the efficiency of the present Stokes solver. One can see that the number of iterations increases slightly even we double the grid size, and the CPU time for a 512× 512 grid is just a few seconds.

5.2. Convergence test for the Stokes flow with an inextensible interface enclosing a solid particle. In this subsection, we perform a convergence study for the present numerical algorithm to the Stokes flow with an inextensible interface enclosing a solid particle. Here, we put an inextensible interface Γ and particle P with initial configuration X(s) = (0.25 cos(s), 0.5 sin(s)) and Y(s) = (0.1 cos(s), 0.1 sin(s)) under a shear flow (u, v) = (χy, 0) in a fluid domain Ω. The dimensionless shear rate χ is chosen to be χ = 1. We also choose the different grid size as m = n = 32, 64, 128, 256, 512 in which the corresponding mesh width is h = 2/m. We also set the Lagrangian mesh widths to be Δs≈ h/2 and Δα ≈ h/2, and the time step size to Δt = h/4.

Since the analytical solution is not available in this test, we choose the result obtained from the ﬁnest mesh m = n = 512 as our reference solution and compute the

(13)

Table 5.3

The mesh reﬁnement results for the perimeter of the interfaceL_h, the interface conﬁguration

Xh, and the velocityu_handv_h.

m = 32 m = 64 rate m = 128 rate m = 256 rate

|Lh− L0|/L0 3.511e-2 2.089e-2 0.75 1.220e-2 0.78 6.749e-3 0.85

Xh− Xref∞ 3.413e-3 1.079e-3 1.66 4.438e-4 1.28 1.622e-4 1.45

uh− uref∞ 6.552e-2 2.592e-2 1.34 1.053e-2 1.30 4.047e-3 1.38

vh− vref∞ 1.133e-1 3.660e-2 1.63 1.554e-2 1.24 4.968e-3 1.64

maximum error between the reference and the numerical solutions. All the numerical solutions are computed up to time T = 0.0625. Since the interface is inextensible, the perimeter of the interface should remain a constant theoretically as time evolves. Let L₀ and L_h be the perimeters of the interface at the initial time and ﬁnal time

T = 0.0625, respectively. The relative error of the perimeter is defined as |L_h− L0|/L0. Table 5.3 shows the relative errors of the perimeter, the maximum errors of the interface configuration, and the maximum errors for the fluid velocity field. Note that the fluid variables are defined at the staggered grid, so when we refine the mesh, the numerical solutions will not coincide with the same grid locations. In these runs, we simply use a linear interpolation to compute the solutions at the desired locations. Due to the fact that the IB formulation has the singular forcing term in the equations, regularizing the singular term by smoothing the discrete delta function causes the method to be first-order accurate. As shown in [3], even though the calculated pressure has O(1) error near the interface, its gradient has the same discretization error as the singular delta force term near the interface. Thus, the overall accuracy of the velocity field in the IB method is first-order accurate in the

L_∞ norm. The numerical results shown in Table 5.3 are consistent with what we expect from theory.

5.3. Tank-treading to tumbling motion under shear flow. Unlike the pre-vious subsection where we focus on the numerical convergence test for our present scheme, here we consider the physical transient deformation of an inextensible inter-face with or without an enclosing rigid particle in a simple shear flow. As mentioned before, the motivation of this test is to simulate the compound vesicle dynamics, which have a lot of applications in bio-fluid problems [26].

In the present numerical experiments, we choose the residual tolerance for GM-RES as 10−4, which is larger than the pure Stokes flow. This is because the elastic tension σ tends to oscillate along the interface, which makes the GMRES method solving (4.11) difficult to converge if the residual tolerance is too small. It will be more appealing if one can find a good preconditioner for solving linear system (4.11) so that GMRES method can converge faster. However, this issue has not yet been resolved.

We first consider the case of an inextensible interface without enclosing a solid particle in a simple shear flow. It is wellknown that the equilibrium dynamics of an inextensible interface or vesicle under a simple shear flow will undergo a tanking-treading motion if the viscosity contrast is under a certain threshold [8]. Here, by tank-treading motion we mean that the configuration of the interface remains stationary while there is a tangential motion along the interface. As studied in previous literature [28, 9, 10, 25, 7], the motion of this steady interface can be characterized by both the inclination angle θ between the long axis of interface and the flow direction, and the

(14)

AN INEXTENSIBLE INTERFACE WITH A SOLID PARTICLE B705 0.5 0.6 0.7 0.8 0.9 0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18 0.19 0.2 reduced area V θ /π χ=1 χ=5 χ=10 0.5 0.6 0.7 0.8 0.9 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 reduced area V f χ=1 χ=5 χ=10

Fig. 5.1_{. The inclination angles} θ/π (left) and the tank-treading frequency f = 2π/

Γu · τdl

(right) versus reduced areasV with diﬀerent dimensionless shear rates for the tank-treading motion

of an inextensible interface in a shear ﬂow.

tank-treading frequency f = 2π/#_Γ _udl

τ of the revolution, where uτ is the tangential velocity component. The inclination angle has been found to be strongly dependent on the reduced area V = 4πA0

L2₀ , where A0 is the enclosed area of the interface and L₀ is the total length of interface. (Notice that, by the above definition, a circle has the reduced area V = 1, while an ellipse with larger aspect ratio has a smaller reduced area.) However, the inclination angle is independent of the dimensionless shear rate χ. This behavior is verified in the left panel of Figure 5.1, which shows the steady inclination angle (θ/π) versus the reduced area (V ) with different shear rates

χ = 1, 5, 10, and it is in a good agreement with previous studies [28, 9, 10, 25, 7]. As

the reduced area increases, the inclination angle increases as well. The right panel of Figure 5.1 shows the tank-treading frequency f versus the reduced area V for diﬀerent shear rates. One can see that as the dimensionless shear rate increases, the tangential motion becomes stronger; thus, the frequency becomes larger. Moreover, by ﬁxing the shear rate, if the interface has larger reduced area, then it has larger frequency as well. Again, our numerical results are in good agreement with previous studies in the literature.

Now we consider the case of an inextensible interface enclosing a solid particle (a compound interface) under a shear flow. As shown in [26], the compound interfacial dynamics will have the transition from tank-treading (TT) to tumbling (TB) if the inclusion effect is strong enough. That is, if the filling fraction φ = a2π/A₀ of the particle (a is the inclusion particle radius, A₀is the interface enclosing area) is above some critical threshold, then the interface will start to tumble rather than being sta-tionary. This can be explained as follows. By including a solid particle, the energy dissipation enhances, so the compound interface behaves like an inclusion-free inter-face that encapsulates a more viscous fluid. The larger the inclusion is, the higher the viscosity will be. Figures 5.2 and 5.3 show the time evolutionary plots for interfacial configurations with different filling fractions, φ = 0.08 and φ = 0.42, respectively. One can indeed see that in the smaller filling fraction case in Figure 5.2, the compound interface has TT motion, while in the larger filling fraction case in Figure 5.3, it has the TB motion.

Meanwhile, in the TT regime, as the ﬁlling fraction increases, both the inclination angle and the TT frequency will decrease. Besides, as in the inclusion-free case, the

(15)

Fig. 5.2_. _{The motion of an compound interface in a shear ﬂow with initial conﬁguration}

X(s) = (0.25 cos(s), 0.5 sin(s)) and Y(s) = (0.1 cos(s), 0.1 sin(s)). φ = 0.08.

Fig. 5.3_. _{The motion of an compound interface in a shear ﬂow with initial conﬁguration}

X(s) = (0.25 cos(s), 0.5 sin(s)) and Y(s) = (0.23 cos(s), 0.23 sin(s)). φ = 0.42.

compound interface with larger reduced area has larger inclination angle and TT frequency when the ﬁlling fraction is small. Both behaviors can be conﬁrmed in Figure 5.4 and are consistent qualitatively with those in [26].

We also investigate the critical value of ﬁlling fraction versus the reduced area for the TT motion to TB transition, as in Figure 5.5. Above the critical value, the interface motion will transit from TT to TB. One can easily see that as the reduced area increases, the critical ﬁlling fraction increases too. This is exactly the similar behavior as the inclusion-free case in which the critical viscosity contrast increases with the reduced area.

(16)

AN INEXTENSIBLE INTERFACE WITH A SOLID PARTICLE B707 0 0.05 0.1 0.15 0.2 filling fraction φ V=0.84 V=0.95 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 filling fraction φ f V=0.84 V=0.95

Fig. 5.4. The inclination angles θ/π (left) and the tank-treading frequency f = 2π/

Γu · τdl

versus ﬁlling fractionφ the in tank-treading regime.

0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 reduced area filling fraction φ TB TT

Fig. 5.5_{. The critical ﬁlling fraction for the TT to TB transition versus the reduced area.}

6. Conclusions. Biological fluid mechanics often require one to study the dy-namics of a lipid membrane encapsulating some cellular contents. For instance, a compound vesicle usually consists of a lipid bilayer membrane enclosing a fluid with a suspended particle. Since the bilayer membrane resists area dilation and its thick-ness is on the order of a nanometer, one can regard the membrane as an inextensible (incompressible) interface. In this paper, we consider the problem of an inextensi-ble interface enclosing a solid particle in Stokes flow. We first write the governing equations in the immersed boundary (IB) framework, in which the interface and solid boundary are treated as force generators in the fluid. In additional to solving for the fluid variables, the present problem involves finding an extra unknown elastic ten-sion such that the surface divergence of the velocity is zero along the interface, and an extra unknown particle surface force such that the velocity satisfies the no-slip boundary condition along the particle surface. We show that the spreading operator of the tension and the surface divergence operator of the velocity are skew-adjoint mathematically. While the interface moves with local fluid velocity, the enclosed par-ticle hereby undergoes a rigid body motion, and the system is closed by the force-free and torque-free conditions along the particle surface.

(17)

The governing equations are then discretized by standard centered difference schemes on a staggered grid, and the interactions between the interface and parti-cle and the fluid are discretized using discrete delta function as in the IB method. The discrete spreading operator of the tension and the discrete surface divergence op-erator of the velocity used in the present scheme preserve the skew-adjoint property numerically. The resultant linear system of equations is therefore symmetric and can be solved by fractional steps so that only fast Poisson solvers are involved. The present method can be extended to Navier–Stokes flow with moderate Reynolds number by treating the nonlinear advection terms explicitly for the time integration. It is also important to mention that, unlike our previous work using the penalty approach [7], we are able to estimate the local error of inextensibility for two successive time steps in the present scheme. In addition, since there are no penalty parameters introduced as in [7], the time step size can be significantly increased.

We study the tank-treading (TT) and tumbling (TB) dynamics for an inextensible interface enclosing a solid particle with different filling fractions under shear flow. We have found that by the increase of the filling fraction, the interface tends to transit from TT to TB. In the TT regime (when the filling fraction is small), the inclination angle and TT frequency will increase as the reduced area increases. Those angle and frequency will decrease as the filling fraction increases. We also have found that the critical filling fraction (from TT to TB) will increase as the reduced area increases, which is qualitatively consistent with the results in [26].

Appendix A. In this appendix, we ﬁrst give a simple derivation for the rigid particle motion equation (2.5) and the torque-free condition (2.7) in twodimensions. These can be derived directly from the following equations in threedimensions,

V = V_c+ω × (Y − Y_c),

∂P

F× (Y − Y_c) dP = 0,

by setting the third component of the vector-valued functions to be zero. For instance, we set V = (V₁, V₂, 0), so the vorticity vectorω = (0, 0, ω).

The symmetric property of the (3.6)–(3.8) can be seen as follows:

V_1k = V_1c− ω(Y_2k− Y_2c), V_2k = V_2c+ ω(Y_1k− Y_1c), Mp K=1 F_1kΔα = 0, Mp k=1 F_2kΔα = 0, Mp k=1 (F_1k(Y_2k− Y2c)− F2k(Y_1k− Y1c)) Δα = 0.

We can scale out the coeﬃcient Δα so that the above three equations can be written in the matrix form

0 R RT 0 F Θ , where R = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 −(Y₂₁− Y_2c) .. . ... ... 1 0 −(Y2M_p− Y2c) 0 1 (Y₁₁− Y1c) .. . ... ... 0 1 (Y_1M_p− Y1c) ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ and Θ = ⎡ ⎢ ⎢ ⎣ V_1c V_2c ω ⎤ ⎥ ⎥ ⎦ .

(18)

Appendix B. In this appendix, we derive that the matrix obtained from the discrete spreading operator S_h of σ_h and the matrix obtained from discrete surface divergence operator ∇_s_h of U are transpose with each other. To see this, we ﬁrst rewrite the operator S_h(σ_h) as

S_h(σ_h) = M−1 k=0 D_s(στ )_k δ_h(x− X_k) Δs = M−1 k=0 σ_k+1/2τ_k+1/2− σ_k−1/2τ_k−1/2 Δs δh(x− Xk)Δs = M k=1 σ_k−1/2τ_k−1/2δ_h(x− X_k−1)− M k=1 σ_k−1/2τ_k−1/2δ_h(x− X_k) − σ−1/2τ−1/2δh(x− X0) + σM−1/2τM−1/2δh(x− XM) = M k=1 (δ_h(x− X_k−1)− δ_h(x− X_k))τ_k−1/2σ_k−1/2.

Note that the last two terms are cancelled out due to the periodicity of the interface. Now we can write the discrete operator∇_s_h as

∇sh· Uk = U_k− U_k−1 Δs · τk−1/2|DsX|k−1/2 =− h 2 Δs |D_sX|_k−1/2 x (δ_h(x− X_k−1)− δ_h(x− X_k))τ_k−1/2· u_i,j.

Since the discrete surface divergence operator of the velocity is zero, we can scale out the coeﬃcient −_{Δs |D}_sX|h2

k−1/2 so that the resultant matrices obtained from Sh and

∇sh· are transpose to each other.

Acknowledgments. We thank Xiaofan Li for helping with the formulation of the motion of the solid particle.

REFERENCES

[1] J. Adams, P. Swarztrauber, and R. Sweet, Fishpack, a package of Fortran subpro-grams for the solution of separable elliptic partial diﬀerential equations, 1980; see http://www2.cisl.ucar.edu/resources/legacy/ﬁshpack.

[2] J. H. Bramble, J. E. Pasciak, and A. T. Vassilev, Analysis of the inexact Uzawa algorithm

for saddle point problems, SIAM J. Numer. Anal., 34 (1997), pp. 1072–1092.

[3] K.-Y. Chen, K.-A. Feng, Y. Kim, and M.-C. Lai, A note on pressure accuracy in immersed

boundary method for Stokes ﬂow, J. Comput. Phys., 230 (2011), pp. 4377–4383.

[4] F. H. Harlow and J. E. Welsh, Numerical calculation of time-dependent viscous

incompress-ible ﬂow of ﬂuid with a free surface, Phys. Fluids, 8 (1965), pp. 2181–2189.

[5] D. Goldstein, R. Handler, and L. Sirovich, Modeling a no-slip ﬂow boundary with an

external force ﬁeld, J. Comput. Phys., 105 (1993), pp. 354–366.

[6] K. H. de Haas, C. Blom, D. van den Ende, M. H. G. Duits, and J. Mellema, Deformation

of giant lipid bilayer vesicles in shear ﬂow, Phys. Rev. E, 56 (1997), pp. 7132–7137.

[7] Y. Kim and M.-C. Lai, Simulating the dynamics of inextensible vesicles by the penalty

im-mersed boundary method, J. Comput. Phys., 229 (2010), pp. 4840–4853.

[8] S. R. Keller and R. Skalak, Motion of a tank-treading ellipsoidal particle in a shear ﬂow, J. Fluid Mech., 120 (1982), pp. 27–47.

[9] V. Kantsler and V. Steinberg, Orientation and dynamics of a vesicle in tank-treading

motion in shear ﬂow, Phys. Rev. Lett., 95 (2005), 258101.

(19)

[10] M. Kraus, W. Wintz, U. Seifert, and R. Lipowsky, Fluid vesicles in shear ﬂow, Phys. Rev. Lett., 77 (1996), pp. 3685–3688.

[11] Z. Li and M.-C. Lai, New ﬁnite diﬀerence methods based on IIM for inextensible interfaces

in incompressible ﬂows, East Asian J. Appl. Math., 1 (2011), pp. 155–171.

[12] M.-C. Lai and C. S. Peskin, An immersed boundary method with formal second order accuracy

and reduced numerical viscosity, J. Comput. Phys., 160 (2000), pp. 705–719.

[13] M.-C. Lai, Y.-H. Tseng, and H. Huang, An immersed boundary method for interfacial ﬂow

with insoluble surfactant, J. Comput. Phys., 227 (2008), pp. 7279–7293.

[14] H. Noguchi and G. Gompper, Shape transitions of ﬂuid vesicles and red blood cells in capillary

ﬂows, Proc. Natl. Acad. Sci. USA, 102 (2005), pp. 14159–14164.

[15] J. B. Perot, An analysis of the fractional method, J. Comput. Phys., 108 (1993), pp. 51–58. [16] C. S. Peskin, The immersed boundary method, Acta Numer., 11 (2002), pp. 479–517.

[17] C. Pozrikidis, The axisymmetric deformation of a red blood cell in uniaxial straining Stokes

ﬂow, J. Fluid Mech., 216 (1990), pp. 231–254.

[18] J. Peters, V. Reichelt, and A. Reusken, Fast iterative solvers for discrete Stokes equations, SIAM J. Sci. Comput., 27 (2005), pp. 646–666.

[19] E. M. Saiki and S. Biringen, Numerical simulation of a cylinder in uniform ﬂow: Application

of a virtual boundary method, J. Comput. Phys., 123 (1996), pp. 450–465.

[20] M. J. Stevens, Coarse-grained simulations of lipid bilayers, J. Chem. Phys., 121 (2004), pp. 11942–11948.

[21] G. W. Schmid-Schonbein, Y. Y. Shih, and S. Chien, Morphometry of human leukocytes, Blood, 56 (1980), pp. 866–875.

[22] J. S. Sohn, Y.-H. Tseng, S. Li, A. Voigt, and J. S. Lowengrub, Dynamics of

multicompo-nent vesicles in a viscous ﬂuid, J. Comput. Phys., 229 (2010), pp. 119–144.

[23] J. C. Strikwerda, An iterative method for solving ﬁnite diﬀerence approximations to the

Stokes equations, SIAM J. Numer. Anal., 21 (1984), pp. 447–458.

[24] K. Taira and T. Colonius, The immersed boundary method: A projection approach, J. Com-put. Phys., 225 (2007), pp. 2118–2137.

[25] S. K. Veerapaneni, D. Gueyffier, D. Zorin, and 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), pp. 2334–2353.

[26] S. K. Veerapaneni, Y.-N. Young, P. M. Vlahovska, and J. Blawzdziewicz, Dynamics of

a compound vesicle in shear ﬂow, Phys. Rev. Lett., 106 (2011), 158103.

[27] X. Yang, X. Zhang, Z. Li, and G.-W. He, A smoothing technique for discrete delta

func-tions with application to immersed boundary method in moving boundary simulafunc-tions, J.

Comput. Phys., 228 (2009), pp. 7821–7836.

[28] H. Zhou and C. Pozrikidis, Deformation of liquid capsules with incompressible interfaces in

simple shear ﬂow, J. Fluid Mech., 283 (1995), pp. 175–200.

(20)

SIAM J. SCI.COMPUT. c XXXX Society for Industrial and Applied Mathematics Vol. 34, No. 5, pp. B692–B710

A FRACTIONAL STEP IMMERSED BOUNDARY METHOD FOR STOKES FLOW WITH AN INEXTENSIBLE INTERFACE

ENCLOSING A SOLID PARTICLE∗

MING-CHIH LAI†, WEI-FAN HU‡, AND WEN-WEI LIN‡

Abstract. In this paper, we develop a fractional step method based on the immersed boundary

(IB) formulation for Stokes flow with an inextensible (incompressible) interface enclosing a solid particle. In addition to solving for the fluid variables such as the velocity and pressure, the present problem involves finding an extra unknown elastic tension such that the surface divergence of the velocity is zero along the interface, and an extra unknown particle surface force such that the velocity satisfies the no-slip boundary condition along the particle surface. While the interface moves with local fluid velocity, the enclosed particle hereby undergoes a rigid body motion, and the system is closed by the force-free and torque-free conditions along the particle surface. The equations are then discretized by standard centered difference schemes on a staggered grid, and the interactions between the interface and particle with the fluid are discretized using a discrete delta function as in the IB method. The resultant linear system of equations is symmetric and can be solved by fractional steps so that only fast Poisson solvers are involved. The present method can be extended to Navier– Stokes flow with moderate Reynolds number by treating the nonlinear advection terms explicitly for the time integration. The convergent tests for a Stokes solver with or without an inextensible interface are performed and confirm the desired accuracy. The tank-treading to tumbling motion for an inextensible interface enclosing a solid particle with different filling fractions under a simple shear flow has been studied extensively, and the results here are in good agreement with those obtained in literature.

Key words. immersed boundary method, inextensible interface, solid particle, fractional step

method, Stokes ﬂow

AMS subject classifications. 65M06, 76D07 DOI. 10.1137/100818777

1. Introduction. Most biological cells consist of a lipid bilayer membrane en-capsulating the cellular content. For instance, a human leukocyte contain a nucleus occupying 18%–44% of the volume [21], which affects the cell adhesion during inflam-matory responses. It is important to study the dynamics of a compound vesicle (a lipid bilayer membrane enclosing a fluid with a suspended particle) [26]. For simplic-ity, we regard the lipid bilayer membrane as an inextensible interface. Thus, in this paper, we develop a fractional step method based on the immersed boundary (IB) formulation for Stokes flow with an inextensible (incompressible) interface enclosing a solid particle.

Most previous studies in the literature have focused on an inextensible interface enclosing a homogeneous fluid. Even in that case, there are still some numerical issues to be addressed. In additional to solving for the fluid variables such as the velocity and pressure, the above problem involves finding extra unknown elastic tension such that ∗_{Submitted to the journal’s Computational Methods in Science and Engineering section December}

20, 2010; accepted for publication (in revised form) August 22, 2012; published electronically October 23, 2012.

http://www.siam.org/journals/sisc/34-5/81877.html

†_{Corresponding author. Center of Mathematical Modeling and Scientiﬁc Computing &}

Depart-ment of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan (mclai@math.nctu.edu.tw). This author’s work was supported in part by National Science Council of Taiwan under research grant NSC-98-2115-M-009-014-MY3, and by the NCTS in Taiwan.

‡_{Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road,}

Hsinchu 300, Taiwan (weifanhu.am95g@g2.nctu.edu.tw, wwlin@math.nctu.edu.tw). B692

(21)