Unconditionally energy stable immersed boundary method with application to vesicle dynamics

(1)

Unconditionally Energy Stable Immersed Boundary Method with Application to Vesicle Dynamics

Wei-Fan Hu and Ming-Chih Lai^∗

Center of Mathematical Modeling and Scientific Computing & Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan.

Received 26 July 2013; Accepted (in revised version) 12 August 2013 Available online XXX 2013

Abstract. We develop an unconditionally energy stable immersed boundary method, and apply it to simulate 2D vesicle dynamics. We adopt a semi-implicit boundary forcing approach, where the stretching factor used in the forcing term can be computed from the derived evolutional equation. By using the projection method to solve the fluid equations, the pressure is decoupled and we have a symmetric positive definite system that can be solved efficiently. The method can be shown to be unconditionally stable, in the sense that the total energy is decreasing. A resulting modification benefits from this improved numerical stability, as the time step size can be significantly increased (the severe time step restriction in an explicit boundary forcing scheme is avoided). As an application, we use our scheme to simulate vesicle dynamics in Navier-Stokes flow.

AMS subject classifications: 65M06, 76D07

Key words: Immersed boundary method, unconditionally energy stable, inextensible vesicle, Navier- Stokes flow.

1. Introduction

The immersed boundary (IB) method proposed by Peskin [26] has been successfully applied to many fluid-structure interaction problems — cf. the review [27]. The IB method employs an Eulerian description for the fluid velocity and a Lagrangian description for the configuration of the immersed elastic structure (immersed boundary or interface). The immersed structure exerts some force into the fluid that drives the fluid flow, and at the same time the fluid flow carries the immersed structure to a new configuration. This interaction between the fluid and the immersed structure is linked through a force spreading and velocity interpolating operator, on using a smoothed version of the Dirac delta function [27].

The IB method is easy to implement and efficient, simply because the immersed structure

∗Corresponding author. Email addresses: weifanhu.am95gg2.n tu.edu.tw (W.-F. Hu), ^m
lai

math.n tu.edu.tw(M.-C. Lai)

http://www.global-sci.org/eajam 1 201x Global-Science Pressc

(2)

(no matter how complex) is regarded as a force generator to the fluid, so that the fluid variables can be solved in a fixed Eulerian domain without generating any structure-fitting grid. Many fast efficient fluid solvers can therefore be applied.

Despite substantial success with practical applications using the IB method, it still has some drawbacks from the numerical point of view. Firstly, the method is only first-order accurate, whereas second-order accurate fluid solvers are used. The immersed elastic structure is usually represented one-dimensionally lower than the fluid space so that the exerted force is singular (delta function like), and smoothing the delta function in a regular finite difference scheme causes the method to be first-order accurate only. Although there have been several attempts to improve accuracy, even some including adaptive local mesh re- finements near the immersed boundary, formally those methods still remain to be made second-order accurate [6, 7, 17, 22, 28].

Another issue is numerical stability. As is well known, the IB method suffers a time step restriction to maintain numerical stability [21, 27, 29, 30]. This restriction becomes more stringent when the elastic force is stiff and the force spreading occurs at the beginning of each time step (an explicit scheme). It is notable that such a time step restriction cannot be alleviated even when the fluid solver is discretised in a semi-implicit manner — i.e. with explicit differencing of the advection term and implicit differencing of the diffusion term.

Rather than performing the force spreading at the beginning of the time step, one might consider doing so at an intermediate stage (a semi-implicit scheme) or even at the end of the time step (an implicit scheme). In the past decade, there have been many attempts to reduce the stiffness or to overcome this time step restriction [3, 4, 8, 10, 11, 23, 24].

However, there is always a trade-off between the stability and efficiency of the algorithms involved. In this article, we propose a new semi-implicit scheme that can be solved quite efficiently, where the resultant linear system is symmetric positive definite and the time step size can be significantly increased.

In Section 2, we introduce the formulation for the incompressible Navier-Stokes equations with an immersed elastic interface. We then develop semi-implicit immersed boundary schemes based on the projection method for the fluid solver in Section 3, and show that these developed schemes are unconditionally energy stable. Then we modify these semi-implicit schemes for efficient implementation, with the resultant linear system symmetric positive definite. Numerical results from simulations of vesicle dynamics are given in Section 4, followed by our conclusions and discussion of future work in Section 5.

2. Governing Equations

We begin by stating the mathematical formulation of the Navier-Stokes flow with an immersed boundary (or interface). We consider a moving, immersed, elastic boundary Γ(t), which exerts forces into an incompressible fluid in a fixed fluid domain Ω. We assume that the fluids inside and outside of the boundary are the same, so the governing equations can be written as follows:

ρ

∂ u

∂ t + (u· ∇)u

+∇p = µ∆u + Z

Γ

F(s, t)δ(x− X(s, t)) ds in Ω , (2.1)

(3)

∇ · u = 0 in Ω , (2.2)

∂ X

∂ t(s, t) = U(s, t) = Z

Ω

u(x, t)δ(x− X(s, t)) dx on Γ . (2.3) Here ρ and µ are constant fluid density and viscosity, u = (u, v) and p are the velocity field and the pressure, both described in Cartesian coordinates. The immersed boundary Γ is represented by a parametric form X(s, t) = (X (s, t), Y (s, t)), where s is the Lagrangian parameter of the initial configuration. The unit tangent vector τ is defined as τ =X_s/|X_s|.

Eqs. (2.1) and (2.2) are the familiar incompressible Navier-Stokes equations, involving a force from the immersed boundary Γ represented by the line integral. Eq. (2.3) defines the position of the immersed boundary moving with the local fluid velocity (the interfacial velocity). Thus the interfacial velocityU is simply an interpolation of the fluid velocity at the immersed boundary, and the interaction between the fluid and the immersed boundary is linked by the two-dimensional Dirac delta functionδ(x) = δ(x)δ( y). We envisage that the force arises from elastic tension according to Hooke’s law — i.e.

F(s, t) = ∂

∂ s(στ) where σ = σ₀(|X_s| − r₀) , (2.4) withσ0is the elastic coefficient and r₀the rest length. Note that the rate of local stretching can then be considered to satisfy the following equation [18]:

∂

∂ t|X_s| = ∇_s· U

|X_s| = ∂ U

∂ s · τ . (2.5)

3. Numerical Discretization

We now proceed to discretise Eqs. (2.1)–(2.5) by the IB method, in a computational rectangular domain Ω = [a, b ]× [c, d ] for simplicity. The fluid variables are defined on the staggered marker-and-cell (MAC) grid introduced by Harlow & Welsh [9] — i.e. the pressure is defined at the cell centre labelled

x = (x_i, y_j) = (a + (i− 1/2)∆x, c + ( j − 1/2)∆ y) , i = 1, 2, · · · , m , j = 1, 2, · · · , n , while the velocity components u and v are defined at the cell edges

(x_i−1/2, y_j) = (a+(i−1)∆x, c+( j−1/2)∆ y) , (x_i, y_j−1/2) = (a+(i−1/2)∆x, c+( j−1)∆ y) . We assume a uniform mesh width h = ∆x = ∆ y in the fluid, although that is not necessary.

For the immersed boundary, we use a collection of discrete points s_k= k∆s, k = 0, 1,· · · M where the mesh width ∆s is chosen to be roughly half of the fluid mesh width h. The Lagrangian markers are denoted byX_k= X(s_k). We also assume the immersed boundary is closed, so that X₀ = X_M. The elastic tension is written σ_k−1/2, since it is defined at the "half-integer" points s_k−1/2 = (k− 1/2)∆s. Without loss of generality, for any function

(4)

r r r r r

→

↑

u_0,1 u_1,1 u_m−1,1 u_m,1

u_0,n u_1,n u_m−1,n u_m,n

v_1,0 v_1,1 v_1,n−1 v_1,n

v_m,0 v_m,1 v_m,n−1 v_m,n

p_1,1 p_m,1

p_1,n p_m,n

Figure1: The omputationaldomainΩ_h^with^a^staggered^grid.

φ(s) defined on the immersed boundary, we approximate the partial derivative ∂ φ/∂ s by central differences — i.e.

D_sφ =φ(s + ∆s/2)− φ(s − ∆s/2)

∆s .

The unit tangent approximated by τ = D_sX/

D_sX

is consequently defined at the half- integer points.

Let ∆t be the time step size, and the superscript index denote the time step level.

At the beginning of each time level n, the values of the boundary configuration Xⁿ_k and the unit tangent τⁿ_k−1/2 are all given. To simplify our notation, we define the discrete spreading operator acting on the tension and discrete surface divergence operator acting on the velocity as follows:

S_h^m[σ](x) =S_h^m[σ₀(|X_s| − r₀)](x) = XM k=1

D_s σ₀(|X_s| − r₀)τ^m

kδ_h(x− X^m_k)∆s , (3.1) T_h^m[u](s_k−1/2) = U_k− U_k−1

∆s · τ^m_k−1/2 U_k=X

x

u(x)δ_h(x− X^m_k)h², (3.2)

(5)

where the superscript m denotes values for the immersed boundary position at the time level m∆t and the discrete operators are skew-adjoint as〈S_h^m[σ], u〉_Ω

h=〈σ, −T_h^m[u]〉_Γ

h

— cf. the derivation in Ref. [19]. The skew-adjointness property plays an important role in our energy stability analysis later.

Our numerical discretisation and energy stability analysis are based on unsteady Stokes equations in this article. The nonlinear advection term in the Navier-Stokes equation (2.1) can be treated explicitly during the time evolution with a moderate CFL condition. How- ever, the Navier-Stokes equations can also be split into an advection part and unsteady Stokes part, where the advection part can be solved by an alternating direction implicit (ADI) method to maintain the unconditionally numerical stability [11].

We introduce two time integrations, the backward Euler (BE) and Crank-Nicholson (CN) scheme, in conjunction with the projection method as follows [5]:

• Backward Euler (BE) scheme

ρu^∗− uⁿ

∆t = µ∆_hu^∗+S_hⁿ[σ₀(|X_s|ⁿ⁺¹− r₀)] , in Ω_h; (3.3)

∆_hφ = 1

∆t∇_h· u^∗, in Ω_h; ∂ φ

∂ n

∂ Ω_h= 0 ; (3.4)

uⁿ⁺¹= u^∗− ∆t∇_hφ , in Ω_h; (3.5)

∇_hpⁿ⁺¹= ρ∇_hφ− µ∇_h(∇_h· u^∗) , in Ω_h; (3.6)

|X_s|ⁿ⁺¹− |X_s|ⁿ

∆t =T_hⁿ[uⁿ⁺¹] , on Γ_h ; (3.7)

Uⁿ⁺¹_k =X

x

u(x)ⁿ⁺¹δ_h(x− Xⁿ_k) h²; (3.8)

Xⁿ⁺¹_k = Xⁿ_k+ ∆t Uⁿ⁺¹_k . (3.9)

• Crank-Nicholson (CN) scheme

ρu^∗− uⁿ

∆t = µ

2∆_h(u^∗+ uⁿ) +S_hⁿ

σ₀

|X_s|ⁿ⁺¹+|X_s|ⁿ

2 − r₀

, in Ω_h; (3.10)

∆_hφ = 1

∆t∇_h· u^∗, in Ω_h; ∂ φ

∂ n

∂ Ω_h= 0 ; (3.11)

uⁿ⁺¹= u^∗− ∆t∇_hφ , in Ω_h; (3.12)

∇_hp^n+1/2= ρ∇_hφ− µ∇_h(∇_h· u^∗) , in Ω_h; (3.13)

|X_s|ⁿ⁺¹− |X_s|ⁿ

∆t =T_hⁿ[(uⁿ⁺¹+ubⁿ)/2] , on Γ_h; (3.14) Uⁿ⁺¹_k =X

x

u(x)ⁿ⁺¹δ_h(x− Xⁿ_k) h², Ubⁿ_k=X

x

u(x)ⁿδ_h(x− Xⁿ_k) h²; (3.15) Xⁿ⁺¹_k = Xⁿ_k+ ∆t(Uⁿ⁺¹_k + bUⁿ_k)/2 . (3.16)

(6)

The gradient and Laplacian spatial operators∇_hand ∆_habove are represented in terms of the standard second-order central difference approximations, andδ_his a smoother version of the discrete delta function developed in Ref. [32]. Unlike the explicit scheme, both of these backward Euler and Crank-Nicholson schemes treat the stretching factor|X_s|ⁿ⁺¹as an unknown solution with closure via the discretisation of Eq. (2.5), so in the immersed boundary force computation the two schemes are semi-implicit.

3.1. Energy stability analysis

We now perform the energy stability analysis for our present schemes in a similar fashion to Refs. [11, 23], and show that both methods are unconditionally stable in the sense that the total energy decreases. To proceed, we define the kinetic energy K and potential energy P of the system as follows:

K =ρ 2kuk²_Ω

h= ρ

2〈u, u〉_Ω_h, P =σ₀

2 k|X_s| − r₀k²_Γ

h= σ₀

2 〈|X_s| − r₀,|X_s| − r₀〉_Γ_h, (3.17) where the associated discrete inner products are defined as

〈u, v〉_Ω

h=X

x

u(x)· v(x) h², 〈φ, ψ〉_Γ

h = XM k=1

φ_k−1/2ψ_k−1/2∆s , (3.18)

respectively. The total energy is E = K + P.

Theorem 3.1. The backward Euler (BE) scheme of Eqs. (3.3)–(3.9) is unconditionally energy stable — i.e. the scheme satisfies Eⁿ⁺¹≤ Eⁿfor each time step n.

Proof. We first consider stability analysis for the BE scheme. For convenience, let us denote the singular force term in Eq. (3.3) by S_hⁿ[σⁿ⁺¹]. We first substitute u^∗ from Eq. (3.5) into Eq. (3.3) to get the equation for ρ(uⁿ⁺¹− uⁿ), and then take the discrete inner product withuⁿ⁺¹+ uⁿ to obtain the kinetic energy difference

Kⁿ⁺¹− Kⁿ=ρ

2〈uⁿ⁺¹,uⁿ⁺¹〉_Ω_h−ρ

2〈uⁿ,uⁿ〉_Ω_h

= ρ

2〈uⁿ⁺¹+ uⁿ,uⁿ⁺¹− uⁿ〉_Ω_h

= ρ 2

−〈uⁿ⁺¹− uⁿ,uⁿ⁺¹− uⁿ〉_Ω_h+ 2〈uⁿ⁺¹,uⁿ⁺¹− uⁿ〉_Ω_h

=−ρ

2kuⁿ⁺¹− uⁿk²_Ω

h+〈uⁿ⁺¹,ρ

uⁿ⁺¹− uⁿ

〉_Ω

h

=−ρ

2kuⁿ⁺¹− uⁿk²_Ω

h+ ∆t〈uⁿ⁺¹,−ρ∇_hφ + µ∆_huⁿ⁺¹+ µ∆t∇_h∆_hφ +S_hⁿ[σⁿ⁺¹]〉_Ω_h

=−ρ

2kuⁿ⁺¹− uⁿk²_Ω

h− ∆tρ〈uⁿ⁺¹,∇_hφ〉_Ω_h+ µ∆t〈uⁿ⁺¹, ∆_huⁿ⁺¹〉_Ω_h + µ∆t²〈uⁿ⁺¹,∇_h∆_hφ〉_Ω_h+ ∆t〈uⁿ⁺¹,S_hⁿ[σⁿ⁺¹]〉_Ω_h.

(7)

The potential energy difference is Pⁿ⁺¹− Pⁿ= σ₀

2

〈|X_s|ⁿ⁺¹− r₀,|X_s|ⁿ⁺¹− r₀〉_Γ_h− 〈|X_s|ⁿ− r₀,|X_s|ⁿ− r₀〉_Γ_h

=σ₀

2 〈|X_s|ⁿ⁺¹+|X_s|ⁿ− 2r₀,|X_s|ⁿ⁺¹− |X_s|ⁿ〉_Γ_h

=σ₀

2 〈2|X_s|ⁿ⁺¹− 2r₀− ∆tT_hⁿ[uⁿ⁺¹], ∆tT_hⁿ[uⁿ⁺¹]〉_Γ

h

as|X_s|ⁿ=|X_s|ⁿ⁺¹−∆tT_hⁿ[uⁿ⁺¹]

=〈σ₀(|X_s|ⁿ⁺¹− r₀)−σ₀∆t

2 T_hⁿ[uⁿ⁺¹], ∆tT_hⁿ[uⁿ⁺¹]〉_Γ_h

= ∆t〈σⁿ⁺¹,T_hⁿ[uⁿ⁺¹]〉_Γ_h− σ₀∆t²

2 kT_hⁿ[uⁿ⁺¹]k²_Γ

h.

Thus the total energy difference between two successive time steps is Eⁿ⁺¹− Eⁿ

=−ρ

2kuⁿ⁺¹− uⁿk²_Ω

h− ∆tρ〈uⁿ⁺¹,∇_hφ〉_Ω_h+ µ∆t〈uⁿ⁺¹, ∆_huⁿ⁺¹〉_Ω_h+ µ∆t²〈uⁿ⁺¹,∇_h∆_hφ〉_Ω_h +∆t〈S_hⁿ[σⁿ⁺¹], uⁿ⁺¹〉_Ω_h+∆t〈σⁿ⁺¹,T_hⁿ[uⁿ⁺¹]〉_Γ_h−σ₀∆t²

h.

The second and fourth terms vanish due to the orthogonality of the discrete divergence-free velocity and the gradient, and the commutativity of∇_hand ∆_h. The fifth and sixth terms cancel out since both discrete operators are skew-adjoint as〈S_hⁿ[σ], u〉_Ω_h=〈σ, −T_hⁿ[u]〉_Γ_h. The third term is always negative since the Laplace operator is negative definite. Thus

Eⁿ⁺¹− Eⁿ=−ρ

2kuⁿ⁺¹− uⁿk²_Ω

h− µ∆tk∇_huⁿ⁺¹k²_Ω

h−σ₀∆t²

h, so the total energy is decreasing such that the present BE scheme is unconditionally energy stable.

Theorem 3.2. The Crank-Nicholson (CN) scheme of Eqs. (3.10)-(3.16) is unconditionally energy stable — i.e. the scheme satisfies Eⁿ⁺¹≤ Eⁿfor each time step n.

Proof. The proof is similar to that for Theorem 3.1. Let us again denote the singular force term in Eq. (3.10) byS_hⁿ[(σⁿ⁺¹+ σⁿ)/2], substitute u^∗in Eq. (3.12) into Eq. (3.10) to getρ(uⁿ⁺¹− uⁿ), and then take the discrete inner product with uⁿ⁺¹+ uⁿto obtain the kinetic energy difference

Kⁿ⁺¹− Kⁿ= ρ

2〈uⁿ⁺¹,uⁿ⁺¹〉_Ω_h− ρ

2〈uⁿ,uⁿ〉_Ω_h

= ρ

2〈uⁿ⁺¹+ uⁿ,uⁿ⁺¹− uⁿ〉_Ω_h

= ∆t

2 〈uⁿ⁺¹+ uⁿ,−ρ∇_hφ +µ

2∆_h(uⁿ⁺¹+ uⁿ) + µ∆t

2 ∇_h∆_hφ +S_hⁿ[(σⁿ⁺¹+ σⁿ)/2]〉_Ω

h.

(8)

The potential energy difference is Pⁿ⁺¹− Pⁿ=σ₀

2

〈|X_s|ⁿ⁺¹− r₀,|X_s|ⁿ⁺¹− r₀〉_Γ_h− 〈|X_s|ⁿ− r₀,|X_s|ⁿ− r₀〉_Γ_h

= σ0

2 〈|X_s|ⁿ⁺¹+|X_s|ⁿ− 2r₀,|X_s|ⁿ⁺¹− |X_s|ⁿ〉_Γ

h

= ∆t

4 〈σ₀(|X_s|ⁿ⁺¹− r₀+|X_s|ⁿ− r₀),T_hⁿ[uⁿ⁺¹+ubⁿ]〉_Γ_h using Eq. (3.14)

= ∆t

4 〈σⁿ⁺¹+ σⁿ,T_hⁿ[uⁿ⁺¹+buⁿ]〉_Γ_h.

Thus the total energy difference between two successive time steps is Eⁿ⁺¹− Eⁿ= −ρ∆t

2 〈uⁿ⁺¹+ uⁿ,∇_hφ〉_Ω_h+µ∆t

4 〈uⁿ⁺¹+ uⁿ, ∆_h(uⁿ⁺¹+ uⁿ)〉_Ω_h +µ∆t²

4 〈uⁿ⁺¹+ uⁿ,∇_h∆_hφ〉_Ω_h+∆t

4 〈uⁿ⁺¹+ uⁿ,S_hⁿ[σⁿ⁺¹+ σⁿ]〉_Ω_h +∆t

4 〈σⁿ⁺¹+ σⁿ,T_hⁿ[uⁿ⁺¹+buⁿ]〉_Γ_h=−µ∆t

4 k∇_h(uⁿ⁺¹+ uⁿ)k²_Ω

h,

The first and third terms vanish due to the orthogonality of discrete divergence-free velocity and the gradient, and the commutativity of∇_hand ∆_h. The fourth and last terms cancel out, since both discrete operators are skew-adjoint. The third term is always negative because the discrete Laplace operator is negative definite. Consequently, the total energy decreases such that the scheme is unconditionally energy stable.

3.2. Modified projection method

To avoid the necessity to solve a coupled linear system of equations in dealing with Eqs. (3.3)–(3.7) or Eqs. (3.10)–(3.14), we can modify the BE and CN schemes slightly. The key idea is to march the time integration of the stretching factor|X_s|ⁿ⁺¹by the intermediate velocityu^∗rather thanuⁿ⁺¹, and substitute it into the discrete spreading operatorS_hⁿsuxh that|X_s|ⁿ⁺¹is no longer among the unknowns in the following schemes..

• Modified BE scheme ρu^∗− uⁿ

∆t = µ∆_hu^∗+S_hⁿ[σⁿ] + σ₀∆tS_hⁿT_hⁿ[u^∗] , in Ω_h; (3.19)

∆_hφ = 1

∆t∇_h· u^∗, in Ω_h; ∂ φ

∂ n

∂ Ω_h= 0 ; (3.20)

uⁿ⁺¹= u^∗− ∆t∇_hφ , in Ω_h; (3.21)

∇_hpⁿ⁺¹= ρ∇_hφ− µ∇_h(∇_h· u^∗)− σ₀∆t²S_hⁿT_hⁿ[∇_hφ] , in Ω_h; (3.22) Uⁿ⁺¹_k =X

x

u(x)ⁿ⁺¹δ_h(x− Xⁿ_k) h², on Γ_h; (3.23)

Xⁿ⁺¹_k = Xⁿ_k+ ∆t Uⁿ⁺¹_k , on Γ_h. (3.24)

(9)

• Modified CN scheme

ρu^∗− uⁿ

∆t = µ

2∆_h(u^∗+ uⁿ) +S_hⁿ[σⁿ] +σ₀∆t

4 S_hⁿT_hⁿ[ubⁿ] +σ₀∆t

4 S_hⁿT_hⁿ[u^∗] , in Ω_h; (3.25)

∆_hφ = 1

∆t∇_h· u^∗, in Ω_h; ∂ φ

∂ n

∂ Ω_h= 0 ; (3.26)

uⁿ⁺¹= u^∗− ∆t∇_hφ , in Ω_h; (3.27)

∇_hp^n+1/2= ρ∇_hφ−µ

2∇_h(∇_h· u^∗)−σ₀∆t²

4 S_hⁿT_hⁿ[∇_hφ] , in Ω_h; (3.28) Uⁿ⁺¹_k =X

x

u(x)ⁿ⁺¹δ_h(x− Xⁿ_k) h², Ubⁿ_k=X

x

u(x)ⁿδ_h(x− Xⁿ_k) h², on Γ_h; (3.29) Xⁿ⁺¹_k = Xⁿ_k+ ∆t(Uⁿ⁺¹_k + bUⁿ_k)/2 , on Γ_h. (3.30)

The lower computational cost of the above numerical schemes is of interest. As we mentioned before, the discrete operatorsS_hⁿ andT_hⁿ are skew-adjoint with matrix forms related byT_hⁿ=−(h²/∆s)(S_hⁿ)^T. Thus the resultant matrix equations foru^∗in Eq. (3.19) and Eq. (3.25) are both symmetric and positive definite, and we can solve efficiently using the aggregation-based multigrid (AGMG) method recently developed by Notay [25]. In Eqs. (3.20) and (3.26), the Poisson equation for pressure increment φ can be solved by the public software package FISHPACK [1]. Despite the fact that these modified BE and CN schemes are not exactly unconditionally energy stable (unlike the original schemes), it is notable that they benefit from the stability such that the time step size ∆t can be tremendously alleviated, as our numerical results below show.

4. Application to Simulating Vesicle Dynamics

Kim & Lai [13] developed a penalty IB method to study the dynamics of inextensible vesicles by introducing a dual representationX(s, t) and Y(s, t) of the immersed boundary, where the first representation X(s, t) interacts with the fluid directly as in traditional IB computation and the second Y(s, t) follows the equations of vesicle dynamics (including the inextensibility constraint) without direct interaction with the fluid dynamics. Both IB representations are linked by stiff springs. The advantage of this penalty idea is that the solution of fluid and vesicle dynamics is decoupled at each time step, so the traditional IB implementation can be applied without much extra effort. Nevertheless, the tension must be solved in this approach, so the penalty numbers has to be chosen sufficiently large to keep the two IB representations close enough, and the time step size must also be small in their explicit forcing computation. Recently, we simplified the penalty approach by introducing a spring-like tension to keep the vesicle boundary nearly inextensible. As a result, the tension (Lagrange multiplier for inextensibility) no longer needs to be found as part of the solution, and we applied this simplified penalty approach to simulate the dynamics of

(10)

three-dimensional axisymmetric vesicles in Navier-Stokes flows [12], but the time step size had to be sufficiently small because the IB force was computed at the beginning of each time step. Here we use our two semi-implicit modified schemes given above that permit much larger time steps, to simulate the vesicle dynamics in Navier-Stokes flows.

Vesicles are closed lipid membranes suspended in a viscous fluid. The membrane force consists of elastic and bending parts, according to

F(s, t) = ∂

∂ s(στ)− c_b∂⁴X

∂ s⁴ , (4.1)

where c_b is the bending rigidity. The tensionσ introduced here is an unknown function, acting as the Lagrange multiplier for enforcing the local inextensibility of the membrane given by

∇_s· U = ∂ U

∂ τ · τ = 0 , on Γ . (4.2)

As in Refs. [2, 12], we replace the unknown tensionσ by σ = σ₀

|X_s| − |X_s|⁰

, (4.3)

where the elastic coefficient σ₀ ≫ 1 and |X_s|⁰ is the local stretching factor of the initial vesicle configuration. This spring-like tension is intended to keep the stretching factor|X_s| close to its initial counterpart|X_s|⁰. (Note that exact inextensibility means|X_s| = |X_s|⁰ for all time.) The numerical schemes presented in Section 3 can be extended straightforwardly to handle the tension computation, by replacing the constant rest length r₀with the initial stretching factor |X_s|⁰. We leave the computational detail for the bending force to the Appendix.

Unless otherwise stated, throughout this section we set the fluid density ρ = 1, the viscosity µ = 1, and the computational domain Ω = [−1, 1] × [−1, 1]. The stopping tolerance of the iterative method AGMG is 10⁻⁶. All numerical runs were carried out on a PC with 16G RAM using double precision arithmetic.

4.1. Convergence study

The convergence of our modified BE and CN schemes is first investigated. We put a vesicle with an elliptical shape X(s, 0) = (0.2 cos(s), 0.5 sin(s)) in quiescent flow initially.

The elastic coefficient is chosen to be σ₀ = 10⁵, and the bending coefficient c_b = 0.01.

We use different mesh sizes where m = n = 64, 128, 256, 512 with the corresponding mesh h = 2/m. We also set the Lagrangian mesh width as ∆s≈ h/2, and the time step

∆t = h/4. Since the fluid is incompressible and the vesicle boundary is inextensible, the enclosed area and the total perimeter of the vesicle should be conserved as time evolves.

For the fluid variables, we choose the result obtained from the finest mesh m = 512 as our reference solution, and compute the maximal error between this reference solution and the numerical solution. All numerical solutions were computed up to time T = 0.125.

We list the results obtained by modified BE and CN schemes in Table 1, showing the relative errors of the vesicle area and the perimeter, the maximum errors of the vesicle

(11)

Table1: Themeshrenementresultsforthevesi leareaA_h^,^vesi
le^perimeterL_h^,^boundary onguration

X_h^,ând^velo
ityômponentsu_hândv_h^.

BE m = n = 64 m = n = 128 rate m = n = 256 rate

|A_h− A₀|/A₀ 1.470e-3 1.252e-3 0.23 7.162e-4 0.81

|L_h− L₀|/L₀ 7.444e-3 2.684e-3 1.47 1.154e-4 1.22 kX_h− Xrefk^∞ 4.209e-3 8.477e-4 2.31 4.762e-4 0.83 ku_h− urefk^∞ 5.453e-2 8.110e-3 2.75 1.578e-3 2.36 kv_h− vrefk^∞ 5.628e-2 1.098e-2 2.36 3.203e-3 1.78 CN m = n = 64 m = n = 128 rate m = n = 256 rate

|A_h− A₀|/A₀ 2.688e-3 2.082e-3 0.37 1.148e-4 0.86

|L_h− L₀|/L₀ 6.899e-3 2.447e-3 1.50 1.052e-4 1.22 kX_h− Xrefk^∞ 3.296e-3 9.643e-4 1.77 6.821e-4 0.50 ku_h− urefk^∞ 7.318e-2 4.919e-2 0.57 2.744e-2 0.84 kv_h− vrefk^∞ 7.561e-2 5.790e-2 0.39 1.935e-2 1.58

boundary configuration, and the fluid velocity field. Since the fluid variables are defined at the staggered grid, the numerical solution on a refined mesh does not coincide with the one obtained from a coarser mesh, and we simply used linear interpolation to compute the solutions at the same grid locations. From Table 1, we note that the numerical results show roughly first-order convergence for all solution variables.

4.2. Maximal time step comparison

Let us now investigate the numerical stability of the present schemes, by testing the maximal time steps for the earlier explicit scheme (EP) versus our modified BE and CN schemes. We used three different elastic coefficientsσ₀ = 10⁷, 10⁸, 10⁹ to study the numerical stability, by comparing the maximal time step that can be used in each scheme.

To determine the maximal time step, we ensured the vesicle boundary behaviours are rea- sonable, by requiring the relative errors of both the area and perimeter to be within 1%

and that there is no numerical instability. The numerical parameters were the same as those in the previous convergence study, except that the initial vesicle configuration was chosen to be X(s) = (0.1 sin(s), 0.5 cos(s)). Table 2 shows the maximal time steps for the three schemes, and it can be seen that the time step for both of the modified BE and CN schemes can be chosen to be 3 to 4 orders larger than the time step in the explicit scheme.

Moreover, as the grid becomes finer or the elastic coefficientσ₀ becomes larger, the time step must become smaller to maintain numerical stability in the explicit scheme. However, for both BE or CN schemes we can always set ∆t = O(h) to maintain the desired numerical stability.

Table 3 shows the average CPU time (in seconds) for each time step, and the total CPU time for the computation up to T = 1. The modified BE and CN schemes clearly outperform the explicit scheme in terms of the total CPU time.

(12)

Table2: Maximumtimestepsforthe expli it(EP)and themodiedBEandCN s hemes, forvarious

elasti oe ientsσ0^.

σ₀= 10⁷ σ₀= 10⁸ σ₀= 10⁹

m, n EP BE CN EP BE CN EP BE CN

128 2.22e-5 1.56e-2 1.56e-2 4.36e-6 1.56e-2 1.56e-2 1.02e-6 1.56e-2 1.56e-2 256 1.03e-5 7.81e-3 7.81e-3 1.74e-6 7.81e-3 7.81e-3 3.81e-7 7.81e-3 7.81e-3 512 4.74e-6 3.91e-3 3.91e-3 7.44e-7 3.91e-3 3.91e-3 1.53e-7 3.91e-3 3.91e-3

Table3: TheaverageCPU time(inse onds)ofea htimestepfordierents hemes. Thetotaltime

with"*"meanstheestimatedvalue.

EP BE CN

m, n ∆t Average Total ∆t Average Total ∆t Average Total

128 1.02e-6 0.07 68627 1.56e-2 0.13 8 1.56e-2 0.13 8

256 3.81e-7 0.16 419947^∗ 7.81e-3 0.39 50 7.81e-3 0.43 55 512 1.53e-7 0.45 2941176^∗ 3.91e-3 1.25 320 3.91e-3 1.45 371

4.3. Tank-treading motion under shear flow

In a more physically realistic application, we consider the transient motion of a vesicle with the initial configurationX(s) = (0.1 sin(s), 0.5 cos(s)) under shear flow u = (γ y, 0).

The elastic coefficient is chosen to beσ₀= 10⁵, and the bending coefficient c_b= 0.01. It is well-known that the equilibrium dynamics of a vesicle under a simple shear flow undergoes tanking-treading motion if the viscosity contrast is under a certain threshold [14], where by tank-treading motion we mean that the configuration of the vesicle remains stationary while there is a tangential motion along the vesicle boundary. Fig. 2 shows the evolutional motion of the vesicle at different times. One can see that after some time the vesicle shape does not change, although the Lagrangian point (marked by∗) along the interface moves with its tangential velocity. This tank-treading motion shown is in good agreement with previous studies [13, 20, 31, 33].

As discussed elsewhere [13, 15, 16, 31, 33], the motion of a steady vesicle under shear flow can be characterised by the inclination angle θ between the long axis of the vesicle and the flow direction and the tank-treading frequency f = 2π/R

Γd l/uτ, where u_τis the tangential velocity component. The inclination angle has been found to depend strongly on the reduced area V = 4πA₀/L₀², where A₀ is the vesicle area and L₀ is the total perimeter of vesicle. (Under this definition, a circle has the reduced area V = 1, while an ellipse with larger aspect ratio has a smaller reduced area.) As the reduced area increases the inclination angle also increases, but the angle is independent of the shear rate γ. This behaviour is verified in the left panel of Fig. 3, which shows the steady inclination angle (θ /π) versus the reduced area (V ) with different shear rates γ = 1, 5, 10. Our numerical results are again in good agreement with those obtained in previous studies [13, 31].

The right panel of Fig. 3 shows the tank-treading frequency f versus the reduced area V for different shear rates. It is evident that as the shear rate increases the tangential

(13)

T = 0.125 T = 0.375 T = 0.625 T = 0.875

T = 13.5 T = 16.5 T = 19.5 T = 22.5

T = 25.5 T = 28.5 T = 31.5 T = 34.5

Figure2: Thetank-treadingmotionofavesi leundershearow.

motion becomes stronger, so the frequency becomes larger. Moreover, on fixing the shear rate, if the vesicle has larger reduced area then it has larger frequency as well. Again, our numerical results are in good agreement with those obtained in previous studies [13, 31].

5. Conclusions and Future Work

We have introduced a semi-implicit boundary forcing approach, where the updated stretching factor is used in the forcing term so that the developed backward Euler and Crank-Nicholson schemes are unconditionally energy stable. By using the projection method to solve the fluid equations, our method decouples the pressure and produces a symmetric positive definite system that can be solved efficiently. The resulting modified method therefore benefits from improved numerical stability, so that the time step size can be significantly increased in contrast to the severe time step restriction for explicit boundary forcing schemes. As an application, we use the developed schemes to simulate vesicle dynamics in Navier-Stokes flow. The extension of our new method to axisymmetric 3D and other 3D immersed boundary problems is currently under investigation.

Acknowledgments

This work has been partially supported by National Science Council of Taiwan under research grant NSC-101-2115-M-009-014-MY3 and NCTS.

(14)

0.5 0.6 0.7 0.8 0.9 0.09

0.1 0.11 0.12 0.13 0.14 0.15 0.16 0.17 0.18

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

Figure 3: The in lination angles θ /π ^(left) ^and ^the tank-treading frequen y f = 2π/R

Γ

dl

u· τ ^(right)

versusredu edareasV ^with^dierent^shear^rates,^for^thetank-treadingmotionofaninextensiblevesi le inashearow.

Appendix

This appendix provides the numerical detail on how we compute the bending force in Eq. (4.1). We first define the discrete force spreading operator F_hⁿ and the velocity interpolating operatorI_hⁿby

F_hⁿ[F](x) = XM k=1

F(s_k)δ_h(x− Xⁿ_k)∆s , I_hⁿ[u](s_k) =X

x

u(x)δ_h(x− Xⁿ_k) h², (A.1) respectively. It is well-known [23, 27] that the above both operators are adjoint with each other — i.e.

〈F_hⁿ[F], u〉_Ω

h=X

x

XM k=1

F(s_k)δ_h(x− Xⁿ_k)∆s

!

· u(x) h²

= XM k=1

F(s_k)·

X

x

u(x)δ_h(x− Xⁿ_k) h²

∆s =〈F, I_hⁿ[u]〉_Γ_h. The singular immersed boundary force arising from the bending is written as

f_b(x) =−c_b Z

Γ

∂⁴X

∂ s⁴δ(x− X(s, t)) ds . (A.2) To simplify our notation, we define the discrete fourth-order centred difference opera- torA_has an approximation to the fourth-order derivative. The discretisation for Eq. (A.2) can thus be written as

fⁿ⁺¹_b (x) =−c_bF_hⁿ[A_hXⁿ⁺¹](x) .

(15)

By substitutingXⁿ⁺¹= Xⁿ+ ∆tI_hⁿ[uⁿ⁺¹] into the above equation, we have fⁿ⁺¹_b (x) =−c_bF_hⁿ[A_hXⁿ+ ∆tA_hI_hⁿ[uⁿ⁺¹]](x)

=−c_bF_hⁿ[A_hXⁿ](x)− c_b∆tF_hⁿA_hI_hⁿ[uⁿ⁺¹](x) . (A.3) Since the discrete operatorsF_hⁿ and I_hⁿ are self-adjoint and the discrete fourth differential operatorA_his symmetric positive definite, the above composite operatorF_hⁿA_hI_hⁿ is also symmetric positive definite. Thus in our modified BE scheme we need only add the additional singular force

f^∗_b(x) =−c_bF_hⁿ[A_hXⁿ](x)− c_b∆tF_hⁿA_hI_hⁿ[u^∗](x) . (A.4) It is notable that the symmetric positive definite matrix structure foru^∗in Eq. (3.19) does not change at all, even when we add this extra bending force term. The bending term for modified CN scheme can be derived similarly, so we omit that detail.

References

[1] J. Adams, P. Swarztrauber and R. Sweet, Fishpack – A Package of Fortran Subprograms for the Solution of Separable Elliptic Partial Differential Equations, (1980).

[2] J. T. Beale, Partially implicit motion of a sharp interface in Navier-Stokes flow, J. Comput. Phys.

231, 6159-6172 (2012).

[3] H. D. Ceniceros, J. E. Fisher and A. M. Roma, Efficient solutions to robust, semi-implicit dis- cretizations of the immersed boundary method, J. Comput. Phys.228, 7137-7158 (2009).

[4] H. D. Ceniceros and J. E. Fisher, A fast, robust and non-stiff immersed boundary method, J.

Comput. Phys.230, 5133-5153 (2011).

[5] J. L. Guermond, P. Minev and J. Shen, An overview of projection methods for incompressible flows, Comput. Methods Appl. Mech. Engrg.195, 6011-6045 (2006).

[6] B. E. Griffith and C. S. Peskin, On the order of accuracy of the immersed boundary method:

Higher order convergence rates for sufficiently smooth problems, J. Comput. Phys.208, 75-105 (2005).

[7] B. E. Griffith, R. D. Hornung, D. M. McQueen and C. S. Peskin, An adaptive, formally second order accurate version of the immersed boundary method, J. Comput. Phys.223, 10-49 (2007).

[8] R. D. Guy and B. Philip, A multigrid method for a model of the implicit immersed boundary equations, Commun. Comput. Phys.12, 378-400 (2012).

[9] F. H. Harlow and J. E. Welsh, Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface, Phys. Fluids8, 2181-2189 (1965).

[10] T.Y. Hou and Z. Shi, Removing the stiffness of elastic force from the immersed boundary method for the 2D Stokes equations, J. Comput. Phys.227, 9138-9169 (2008).

[11] T. Y. Hou and Z. Shi, An efficient semi-implicit immersed boundary method for the Navier-Stokes equations, J. Comput. Phys.227, 8968-8991 (2008).

[12] W.-F. Hu, Y. Kim and M.-C. Lai, An immersed boundary method for simulating the dynamics of three-dimensional axisymmetric vesicles in Navier-Stokes flows, submitted for publication.

[13] Y. Kim and M.-C. Lai, Simulating the dynamics of inextensible vesicles by the penalty immersed boundary method, J. Comput. Phys.229, 4840-4853 (2010).

[14] S. R. Keller and R. Skalak, Motion of a tank-treading ellipsoldal particle in a shear flow, J. Fluid Mech.120, 27-47 (1982).