Unconditionally energy stable schemes for the inextensible interface problem with bending

(1)

UNCONDITIONALLY ENERGY STABLE SCHEMES FOR THE

INEXTENSIBLE INTERFACE PROBLEM WITH BENDING^∗

MING-CHIH LAI^† AND KIAN CHUAN ONG^‡

Abstract. In this paper, we develop unconditionally energy stable schemes to solve the inextensible interface problem with bending. The fundamental problem is formulated by the immersed boundary method where the nonstationary Stokes equations are considered, with the elastic tension and bending forces expressed in terms of Dirac delta function along the interface. The elastic tension is one of the solution variables which plays the role of Lagrange multiplier to enforce the inextensibility of the interface. Both the backward Euler and Crank–Nicolson methods are introduced and it can be proved that the total energy, i.e., kinetic energy and bending energy, is discretely bounded. The numerical results show that both schemes are unconditionally energy stable without any time-step restriction. The backward Euler scheme is also applied to study the dynamics of vesicles suspended in a shear flow.

Key words. immersed boundary method, unconditionally energy stable scheme, inextensible interface, bending

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

1. Introduction. The immersed boundary (IB) method proposed by Peskin [1] has been successfully applied to various fluid-structure interaction problems. A comprehensive review of the IB method is presented in [2]. The IB formulation utilizes an Eulerian frame for the fluid velocity and a Lagrangian frame for the configuration of the immersed elastic structure (immersed boundary or interface). The immersed structure exerts certain forces into the fluid that drive the fluid flow, and at the same time, the fluid flow carries the immersed structure to a new configuration. The fluid- structure interaction between the fluid and the immersed structure is coupled through a force spreading and a velocity interpolating operator by the usage of the smoothed version of the Dirac delta function [2]. The IB formulation is efficient and easy to implement because the immersed structure (regardless of complexity) is regarded as a force generator to the fluid so the fluid variables can be solved in a static Eulerian domain without generating any structure-conforming grid. Therefore, a variety of efficient fluid solvers can be exploited.

Although substantial success has been achieved in the practical applications using the IB method, it still possesses several deficiencies from the numerical point of view.

First, the IB method is restricted to first-order accurate even though second-order accurate fluid solvers are used. Since the immersed elastic structure is usually one- dimensional lower than the fluid space, the exerted force is singular (delta function like) and smoothing the delta function in the regular finite difference scheme limits the

∗Submitted to the journal’s Computational Methods in Science and Engineering section August 28, 2018; accepted for publication (in revised form) April 29, 2019; published electronically July 2, 2019.

https://doi.org/10.1137/18M1210277

Funding: The work of the first author was partially supported by the Ministry of Science of Technology of Taiwan under research grant MOST-107-2115-M-009-016-MY3 and by the National Center for Theoretical Sciences.

†Corresponding author. Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan ([email protected]).

‡National Center for Theoretical Sciences, No. 1, Sec. 4, Roosevelt Road, Astronomy-Mathematics Building, National Taiwan University, Taipei 10617, Taiwan ([email protected]).

Downloaded 07/04/19 to 140.113.22.164. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

(2)

method to being first-order accurate only. Although several attempts have been made to improve the accuracy, including adaptive local mesh refinements near the immersed boundary, those methods remain formally second-order accurate [3, 4, 5, 6, 7].

Another numerical issue arises from the stability. As revealed in the literature [8, 9, 10, 2], the IB method suffers a stringent time-step restriction to maintain its numerical stability. This restriction becomes severe when the elastic force is stiff, and the force spreading is evaluated at the beginning of each time-step (explicit IB scheme).

Such a time-step restriction cannot be elevated effectively even if the fluid solver is discretized in a semi-implicit manner, i.e., explicit differencing of the advection term and implicit differencing of the diffusion term. Rather than performing the force spreading at the beginning of each time-step, one of the remedies is to perform the procedure at each intermediate (semi-implicit IB scheme) stage or at the end of each time-step (implicit IB scheme). However, there is always a trade-off between the stability and the efficiency of those algorithms. The fully implicit IB schemes lead to a nonlinear system of equations which limit the range of applicable solvers, while the semi-implicit IB schemes result in a linear system of equations which provide the possibilities of efficient solvers. In the past decade, there have been many attempts to reduce the stiffness or to overcome the severe time-step restriction in the development of the fully implicit IB scheme and the semi-implicit IB scheme [11, 12, 14, 15, 16, 17, 18]. It was ultimately identified that numerical instability is attributed to a lack of energy conservation [11, 12]. The semi-implicit IB schemes proposed by [11] are proven to be unconditionally stable such that the discrete energy is bounded, considering the discrete force spreading and interpolation operators are discretized at the same location in space and time. Subsequently, within the context of unconditionally stable semi- implicit IB schemes, several efficient algorithms and linear solvers are explored. For instance, approximate projection methods [11], the double Schur complement [12], the geometric multigrid method [18], and the projection method with preconditioners [19].

In past decades, the study of the morphological dynamics of biological cell mem- branes and vesicles has been an active source of mathematical modeling in biology, biophysics, and bioengineering. The biological system poses challenges to develop efficient and accurate numerical schemes, and the IB method provides a versatile approach to tackle the problem [20, 21, 22, 23, 24, 25, 26]. In the present study, we shall propose new semi-implicit IB schemes for solving the inextensible interface problem with bending that were proven to be unconditionally stable, i.e., bounded discrete energy. The present method and stability analysis differ from the first author’s pre- vious work [13] in which the nearly inextensible approach is adopted and no bending effect is considered there. The rest of the paper is organized as follows. In section 2, we introduce the formulation of incompressible Stokes equations with an immersed inextensible interface. Then, we provide the continuous energy estimate in section 3.

In section 4, we develop semi-implicit IB schemes based on the backward Euler (BE) and Crank–Nicolson (CN) schemes and show that those developed schemes are unconditionally energy stable. The details of numerical methods are described in section 5, followed by numerical results consisting of energy stability check, grid convergence studies, and an application of vesicle dynamics in section 6. Finally, some conclusions are presented in section 7.

2. Equations of motion. We begin by stating the governing equations for an enclosed inextensible interface with bending immersed in a two-dimensional viscous incompressible fluid domain Ω. This model is motivated by the vesicle problem that has attracted significant attention in past years. Throughout this paper, the inextensible

(3)

interface is also called the vesicle boundary. Here, we assume that the vesicle and the fluid have the properties of matched density and viscosity. The equations of motion in IB formulation can be written as

ρ∂u

∂t + ∇p = µ∆u + S(Fσ) + S(Fb) in Ω, (2.1)

∇ · u = 0 in Ω, (2.2)

S(Fσ) = Z

Γ

∂

∂s(στ )δ(x − X(s, t)) ds, (2.3)

S(Fb) = − Z

Γ

cb

∂⁴X

∂s⁴ δ(x − X(s, t)) ds, (2.4)

∂X

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

Ω

u(x, t)δ(x − X(s, t))dx on Γ, (2.5)

∇s· U = ∂U

∂s · τ = 0 on Γ, (2.6)

where ρ is the density, u the velocity, p the pressure, and µ the viscosity. Here, the vesicle boundary configuration Γ is an enclosed interface represented by Lagrangian markers X(s, t) with s the arc-length parameter. Since the vesicle is inextensible (2.6), the total arc-length of Γ is conserved. The notation τ (s, t) = ^∂X_∂s is the unit tangent vector and cb is the constant bending rigidity. The forces arising from the vesicle boundary Γ have two terms, namely, the elastic tension force Fσ in (2.3) and the bending force Fb in (2.4), which can be derived by taking the variational derivative of Helfrich type energy as shown in [27]. Note that the elastic tension σ(s, t) is a self-adjusting unknown function (Lagrange multiplier) to enforce the inextensibility constraint (2.6) which is exactly similar to the role of pressure to enforce the fluid incompressibility (2.2). Certainly, the above governing equations should be ac- companied with initial conditions for the velocity and vesicle boundary configuration, and suitable velocity boundary condition on ∂Ω. For simplicity, the periodic boundary condition is employed throughout the study, particularly in the following energy estimates. As far as we know, the solvability of the whole continuous system of (2.1)–

(2.6) has not yet been rigorously investigated. Nevertheless, the well-posedness of the steady problem, associated to the system in the absence of bending force (2.4), has been studied recently by Liu, Song, and Xu [28] using weak formulation, and the inf- sup condition has been established successfully. Certainly, the proof of well-posedness of (2.1)–(2.6) is a very important theoretical problem and we shall leave it for future work.

3. Continuous energy estimate. We define the inner product of functions on Ω and Γ in the following:

(3.1) hu, viΩ= Z

Ω

u(x, t) · v(x, t) dx, hf, giΓ = Z

Γ

f (s, t) g(s, t) ds.

Thus, the L² norm of functions on Ω and Γ can be defined as kuk²_Ω = hu, ui_Ω and kf k²_Γ = hf, f i_Γ, respectively. Note that the inner product on Γ can be extended to the vector-valued functions by hf , gi_Γ =R

Γf (s, t) · g(s, t) ds so that kf k²_Γ = hf , f i_Γ. Taking the inner product of (2.1) by u and integrating in Ω yields

(3.2) ρ ∂u

∂t, u

Ω

+ h∇p, ui_Ω= hµ∆u, ui_Ω+ hS(F_σ) + S(F_b), ui_Ω.

(4)

Invoking the incompressibility condition (2.2) and periodic boundary condition for the velocity, the above equation can be simplified as

(3.3) d

dt ρ

2kuk²_Ω= −µk∇uk²_Ω+ hS(Fσ) + S(Fb), uiΩ.

In the IB formulation, the spreading operator acting on the elastic tension force Fσand the surface divergence operator of the velocity are skew-adjoint to each other;

see the references in [21, 30]. For completeness, we provide the derivation herein.

The spreading operator acting on the elastic tension force reads hS(Fσ), ui_Ω=

Z

Ω

Z

Γ

∂

∂s(στ )δ(x − X(s, t)) ds

· u(x, t) dx

= Z

Γ

∂

∂s(στ ) ·

Z

Ω

u(x, t)δ(x − X(s, t)) dx

ds

= − Z

Γ

σ

τ · ∂U

∂s

ds (applying the integration by parts)

= Z

Γ

σ

−∂U

∂s · τ

ds = hσ, −∇s· UiΓ= 0 (by (2.6)).

(3.4)

As mentioned, analogous to pressure, the unknown elastic tension σ solely acts as a Lagrange multiplier to enforce the inextensibility constraint (2.6) on the vesicle boundary Γ.

The spreading operator acting on the bending force reads hS(Fb), uiΩ=

Z

Ω

− Z

Γ

cb

∂⁴X

∂s⁴δ(x − X(s, t)) ds

· u(x, t) dx

= − Z

Γ

cb

∂⁴X

∂s⁴ ·

Z

Ω

u(x, t)δ(x − X(s, t)) dx

ds

= − Z

Γ

cb

∂⁴X

∂s⁴ · U(s, t) ds

= − Z

Γ

cb

∂²X

∂s² · ∂²U

∂s² ds (applying integration by parts twice)

= − Z

Γ

cb

∂²X

∂s² · ∂²Xt

∂s² ds (by (2.5))

= −d dt

Z

Γ

cb

2

∂²X

∂s² ·∂²X

∂s² = −d dt

Z

Γ

cb

2

∂²X

∂s²

2

ds.

(3.5)

By substituting (3.4)–(3.5) into (3.3), we obtain the following energy estimate:

(3.6) d

dt ρ

2kuk²_Ω+c_b 2

∂²X

∂s²

2

Γ

!

= −µk∇uk²_Ω.

where the total energy of the system consists of the fluid kinetic energy (the first term) and the vesicle bending energy (the second term).

4. Numerical discretization. Now, we are ready to discretize (2.1)–(2.6) by the IB method. For simplicity, we assume a computational rectangular domain Ω = [0, Lx] × [0, Ly]. The fluid variables are defined on the staggered marker-and-cell

(5)

(MAC) grid introduced by Harlow and Welsh [29]; precisely, the pressure is defined at the cell center labeled as x = (xi, yj) = ((i − 1/2)∆x, (j − 1/2)∆y), i = 1, 2, . . . , m and j = 1, 2, . . . , n, while the velocity components u and v are defined at cell edges (x_i−1/2, yj) = ((i − 1)∆x, (j − 1/2)∆y) and (xi, y_j−1/2) = ((i − 1/2)∆x, (j − 1)∆y), respectively. Here, though it is not necessary, we assume a uniform mesh width h = ∆x = ∆y. For the immersed vesicle boundary, we use a collection of M discrete points sk= k∆s, k = 0, 1, . . . M −1, with mesh width ∆s comparable to the grid mesh h. The Lagrangian markers are denoted by Xk = X(sk). Since the vesicle boundary is closed, we define XM = X0. The elastic tension is defined at the “half-integer”

point given by sk−1/2 = (k − 1/2)∆s so we denote it as σk−1/2. For any function φ(s) defined on the immersed interface, we introduce three different approximations to the partial derivative ^∂φ_∂s, namely, the forward, backward, and central difference schemes, as

D_s⁺φ(s) = φ(s + ∆s) − φ(s)

∆s , D⁻_sφ(s) = φ(s) − φ(s − ∆s)

∆s ,

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

∆s .

Thus, the unit tangent can be approximated by τ = D_sX, which in turn is defined at the “half-integer” points.

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

At the beginning of each time level n, the boundary configuration Xⁿ_k, and the unit tangent τⁿ_k−1/2 are all given. In the present study, our numerical discretization and energy stability analysis are based on the unsteady Stokes equations instead of the Navier–Stokes equations. For the latter case, the nonlinear advection term can be treated explicitly during the time evolution with a moderate CFL condition. Al- ternatively, the Navier–Stokes equations can be split into an advection part and an unsteady Stokes part in which the advection equation is solved by the alternating direction implicit method to maintain the unconditionally numerical stability [15].

Here we introduce two different time integrations, namely, the BE scheme,

ρuⁿ⁺¹− uⁿ

∆t + ∇_hpⁿ⁺¹= µ∆_huⁿ⁺¹+ S_hⁿ(Fⁿ⁺¹_σ ) + S_hⁿ(Fⁿ⁺¹_b ) in Ω_h, (4.2)

∇h· uⁿ⁺¹= 0 in Ωh, (4.3)

S_hⁿ(Fⁿ⁺¹_σ ) =

M

X

k=1

Ds σⁿ⁺¹τⁿ

kδh(x − Xⁿ_k)∆s, (4.4)

S_hⁿ(Fⁿ⁺¹_b ) = −c_b

M

X

k=1

D_s⁺D⁻_sD⁺_sD_s⁻Xⁿ⁺¹_k δ_h(x − Xⁿ_k)∆s, (4.5)

Uⁿ⁺¹_k =X

x

u(x)ⁿ⁺¹δh(x − Xⁿ_k) h² ∀k, (4.6)

∇sh· Uⁿ⁺¹_k = D⁻_sUⁿ⁺¹_k · τⁿ_k−1/2= 0 ∀k, (4.7)

Xⁿ⁺¹_k − Xⁿ_k

∆t = Uⁿ⁺¹_k ∀k, (4.8)

and the CN scheme

(6)

ρuⁿ⁺¹− uⁿ

∆t + ∇_hpⁿ⁺¹^/²= µ

2∆_h(uⁿ⁺¹+ uⁿ) + S_hⁿ⁻¹^/²(Fⁿ⁺_σ ¹^/²) + S_hⁿ⁻¹^/²(Fⁿ⁺_b ¹^/²), (4.9)

∇h· uⁿ⁺¹= 0, (4.10)

S_hⁿ⁻¹^/²(Fⁿ⁺_σ ¹^/²) =

M

X

k=1

D_s

σⁿ⁺¹^/²τⁿ⁻¹^/²

k

δ_h(x − Xⁿ⁻_k ¹^/²)∆s, (4.11)

S_hⁿ⁻¹^/²(Fⁿ⁺_b ¹^/²) = −cb M

X

k=1

D⁺_sD_s⁻D_s⁺D⁻_sXⁿ⁺_k ¹^/²δh(x − Xⁿ⁻_k ¹^/²)∆s, (4.12)

Uⁿ⁺_k ¹^/²=X

x

u(x)ⁿ⁺¹+ u(x)ⁿ

2 δh(x − Xⁿ⁻_k ¹^/²) h² ∀k, (4.13)

∇s_h· Uⁿ⁺_k ¹^/²= D_s⁻Uⁿ⁺_k ¹^/²· τⁿ⁻_k−₁¹_/^/₂²= 0 ∀k, (4.14)

Xⁿ⁺¹_k − Xⁿ_k

∆t = Uⁿ⁺_k ¹^/² ∀k.

(4.15)

The spatial operators ∇_hand ∆_hare the standard second-order central difference approximations to the gradient and Laplacian on the MAC grid. δ_h is the discrete delta function. Note that Xⁿ⁺_k ¹^/², Xⁿ⁻_k ¹^/², and τⁿ⁻_k−₁¹_/^/²₂ in the CN scheme are lin- early approximated by Xⁿ⁺_k ¹^/² = (Xⁿ_k + Xⁿ⁺¹_k )/2, Xⁿ⁻_k ¹^/² = (Xⁿ_k + Xⁿ⁻¹_k )/2, and τⁿ⁻_k−₁¹_/^/₂² = (τⁿ_k−₁_/₂+ τⁿ⁻¹_k−₁_/₂)/2, respectively. The spreading operator S_h^n−1/2 used in the CN scheme is for the linearization of S_h^n+1/2by lagging one discrete time-step ∆t for the interface position at time t = (n − 1/2)∆t instead of t = (n + 1/2)∆t. This is exactly like the BE scheme that uses the spreading operator S_hⁿ for the interface position at t = n∆t instead of t = (n + 1)∆t. To be consistent, the discrete inextensibility constraint in present schemes also uses the one time-step lagging interface tangents, which linearizes the constraint as well. Here, we want to emphasize that the lagged operator is popularly used in IB simulations. Lagging the spreading and interpolation operators results in a linear system of equations which provides a variety of linear solvers applicable to the problem. Not lagging those operators results in a set of nonlinear equations for the implicit system [7]. This time lagging technique might reduce accuracy; however, a predictor-corrector approach which uses the lagged interface position to construct the spreading operator S in the predictor step and the predicted interface to construct the same operator in the corrector step does not exactly recover to full second-order accuracy. (This kind of predictor-corrector approach is called formally second-order accurate, and it is a fully second-order method only when the sharp interface is replaced by an elastic shell with thickness [4, 7, 22].) In fact, it is well known that the IB method is only first-order accurate due to the discretization of the singular delta function force in the formulation [2]. The authors have implemented the predictor-corrector IB method as in [4] for the vesicle in shear flow and the results are not significantly better than the ones shown in Table 3, so we omit them here. Hence, the time lagging technique is crucial for the discretization of the IB problem because it not only yields a linear system for the resultant matrix but also provides flexibility in solving the discretized equations and facilitates the foundation of the present energy stability analysis, as seen below.

(7)

4.1. Energy stability analysis. In this subsection, we follow the similar energy analysis proposed in [11, 15] and perform the energy stability analysis for our present schemes. We shall show that both methods are unconditionally stable in the sense that total energy is decreasing.

To proceed, we define the kinetic energy K and bending energy B of the system as (4.16)

K =ρ

2kuk²_Ω_h= ρ

2hu, uiΩ_h, B =c_b

2kD⁺_sD_s⁻Xk²_Γ_h= c_b

2hD_s⁺D_s⁻X, D⁺_sD_s⁻XiΓ_h, where the associated discrete inner products are defined as

(4.17) hu, viΩh=X

x

u(x) · v(x) h², hφ, ψiΓh =

M

X

k=1

φ_kψ_k∆s, respectively. Thus, the total energy is E = K + B.

Theorem 4.1. The BE scheme of (4.2)–(4.8) is unconditionally energy stable;

that is, the scheme satisfies Eⁿ⁺¹≤ Eⁿ for each time-step n.

Proof. Taking the discrete inner product of (4.2) with uⁿ⁺¹+ uⁿ, we obtain Kⁿ⁺¹−Kⁿ =ρ

2huⁿ⁺¹, uⁿ⁺¹iΩ_h−ρ

2huⁿ, uⁿiΩ_h

=ρ

2huⁿ⁺¹+ uⁿ, uⁿ⁺¹− uⁿiΩ_h

=ρ

2 −huⁿ⁺¹− uⁿ, uⁿ⁺¹− uⁿi_Ω_h+ 2huⁿ⁺¹, uⁿ⁺¹− uⁿi_Ω_h

= −ρ

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

h+ huⁿ⁺¹, ρ uⁿ⁺¹− uⁿiΩ_h

= −ρ

2kuⁿ⁺¹− uⁿk²_Ω_h+ ∆thuⁿ⁺¹, −∇hpⁿ⁺¹+ µ∆huⁿ⁺¹ + S_hⁿ(Fⁿ⁺¹_σ ) + S_hⁿ(Fⁿ⁺¹_b )i_Ω_h.

Due to the discrete incompressibility condition (4.3) and the discrete inextensibility condition (4.7), the terms of the pressure and the elastic tension are canceled out.

Precisely, we have huⁿ⁺¹, −∇hpⁿ⁺¹iΩ_h = 0 (the proof can be performed easily by using the condition (4.3) and the summation by parts) and huⁿ⁺¹, S_hⁿ(Fⁿ⁺¹_σ )i_Ω_h =

−hσⁿ⁺¹, ∇_s_h· Uⁿ⁺¹iΓh = 0 (the proof can be found in [21]). Thus, the above equation becomes

Kⁿ⁺¹−Kⁿ= −ρ

2kuⁿ⁺¹− uⁿk²_Ω_h+ µ∆thuⁿ⁺¹, ∆huⁿ⁺¹iΩ_h+ ∆thuⁿ⁺¹, S_hⁿ(Fⁿ⁺¹_b )iΩ_h

= −ρ

2kuⁿ⁺¹−uⁿk²_Ω_h−µ∆th∇huⁿ⁺¹, ∇_huⁿ⁺¹iΩh+∆thuⁿ⁺¹, S_hⁿ(Fⁿ⁺¹_b )i_Ω_h. (4.18)

Now we need to compute the last term huⁿ⁺¹, S_hⁿ(Fⁿ⁺¹_b )iΩ_h, huⁿ⁺¹, S_hⁿ(Fⁿ⁺¹_b )iΩ_h =X

x

−cb M

X

k=1

D_s⁺D⁻_sD⁺_sD_s⁻Xⁿ⁺¹_k δh(x − Xⁿ_k)∆s

!

· uⁿ⁺¹(x) h²

= −c_b

M

X

k=1

D_s⁺D⁻_sD⁺_sD⁻_sXⁿ⁺¹_k · X

x

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

!

∆s

= −c_b

M

X

k=1

D_s⁺D⁻_sD⁺_sD⁻_sXⁿ⁺¹_k · Uⁿ⁺¹_k ∆s

(8)

= −cbhD_s⁺D⁻_sD⁺_sD⁻_sXⁿ⁺¹, Uⁿ⁺¹iΓ_h

= −cbhD_s⁺D⁻_sXⁿ⁺¹, D⁺_sD⁻_sUⁿ⁺¹iΓ_h. (4.19)

Note that the last identity is obtained by hD_s⁺φ, ψiΓ_h = −hφ, D_s⁻ψiΓ_h using the fact of summation by parts and the periodicity of the Γh.

The discrete bending energy is Bⁿ⁺¹− Bⁿ= c_b

2hD⁺_sD⁻_sXⁿ⁺¹, D_s⁺D_s⁻Xⁿ⁺¹iΓh−c_b

2hD_s⁺D⁻_sXⁿ, D⁺_sD_s⁻XⁿiΓh

= cb

2hD⁺_sD⁻_s(Xⁿ⁺¹+ Xⁿ), D_s⁺D_s⁻(Xⁿ⁺¹− Xⁿ)iΓ_h

= −cb

2hD⁺_sD_s⁻(Xⁿ⁺¹− Xⁿ), D_s⁺D⁻_s(Xⁿ⁺¹− Xⁿ)iΓ_h

+ cbhD⁺_sD_s⁻Xⁿ⁺¹, D⁺_sD⁻_s(Xⁿ⁺¹− Xⁿ)iΓ_h

= −c_b

2kD_s⁺D⁻_s(Xⁿ⁺¹− Xⁿ)k²_Γ

h+c_bhD⁺_sD_s⁻Xⁿ⁺¹, D⁺_sD⁻_s(Xⁿ⁺¹− Xⁿ)i_Γ_h

= −cb

2kD_s⁺D⁻_s(Xⁿ⁺¹− Xⁿ)k²_Γ

h+ cb∆thD⁺_sD⁻_sXⁿ⁺¹, D⁺_sD_s⁻Uⁿ⁺¹iΓ_h

= −cb

2kD_s⁺D⁻_s(Xⁿ⁺¹− Xⁿ)k²_Γ_h−∆thuⁿ⁺¹, S_hⁿ(Fⁿ⁺¹_b )iΩ_h(by (4.19)).

Thus, the total energy between two successive time-steps can be written as (4.20)

Eⁿ⁺¹− Eⁿ= −ρ

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

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

h−cb

2kD⁺_sD⁻_s(Xⁿ⁺¹− Xⁿ)k²_Γ

h. Thus, the total energy is decreasing, which shows the present BE scheme is unconditionally energy stable.

Theorem 4.2. The CN scheme of (4.9)–(4.15) is unconditionally energy stable;

that is, the scheme satisfies Eⁿ⁺¹≤ Eⁿ for each time-step n.

Proof. Taking the discrete inner product of (4.9) with uⁿ⁺¹+ uⁿ, we obtain Kⁿ⁺¹− Kⁿ=ρ

2huⁿ⁺¹, uⁿ⁺¹iΩ_h−ρ

2huⁿ, uⁿiΩ_h

=ρ

2huⁿ⁺¹+ uⁿ, uⁿ⁺¹− uⁿiΩ_h

=∆t

2 huⁿ⁺¹+ uⁿ, −∇_hpⁿ⁺¹^/²i_Ω_h+µ∆t

4 huⁿ⁺¹+ uⁿ, ∆_h(uⁿ⁺¹+ uⁿ)i_Ω_h +∆t

2 huⁿ⁺¹+uⁿ, S_hⁿ⁻¹^/²(Fⁿ⁺_σ ¹^/²)iΩ_h+∆t

2 huⁿ⁺¹+uⁿ, S_hⁿ⁻¹^/²(Fⁿ⁺_b ¹^/²)iΩ_h

= −µ∆t

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

h+∆t

2 huⁿ⁺¹+ uⁿ, Sⁿ⁻_h ¹^/²(Fⁿ⁺_b ¹^/²)i_Ω_h, where the pressure term satisfies huⁿ⁺¹+ uⁿ, −∇hpⁿ⁺¹^/²iΩ_h = 0 by using (4.10) and

∇h· uⁿ= 0, and the elastic tension term satisfies

huⁿ⁺¹+ uⁿ, S_hⁿ⁻¹^/²(Fⁿ⁺_σ ¹^/²)iΩ_h= −

M

X

k=1

σ_k−1/2ⁿ⁺¹^/²τⁿ⁻_k−1/2¹^/² · D_s⁻Uⁿ⁺_k ¹^/²∆s

= −hσⁿ⁺¹^/², ∇s_h· Uⁿ⁺¹^/²iΓ_h = 0.

(4.21)

The proof of above equality can be easily extended from the derivation in [21].

(9)

Now we need to compute the last term huⁿ⁺¹+ uⁿ, S_hⁿ⁻¹^/²(Fⁿ⁺_b ¹^/²)iΩ_h, huⁿ⁺¹+ uⁿ, S_hⁿ⁻¹^/²(Fⁿ⁺_b ¹^/²)iΩ_h

=X

x

−c_b

M

X

k=1

D⁺_sD⁻_sD⁺_sD_s⁻Xⁿ⁺_k ¹^/²δ_h(x − Xⁿ⁻_k ¹^/²)∆s

!

· (uⁿ⁺¹(x) + uⁿ(x)) h²

= −cb M

X

k=1

D⁺_sD⁻_sD⁺_sD_s⁻Xⁿ⁺_k ¹^/² · X

x

(uⁿ⁺¹(x) + uⁿ(x))δh(x − Xⁿ⁻_k ¹^/²) h²

!

∆s

= −2cb M

X

k=1

D⁺_sD_s⁻D_s⁺D⁻_sXⁿ⁺_k ¹^/²· Uⁿ⁺_k ¹^/²∆s (by (4.13))

= −2cbhD⁺_sD_s⁻D_s⁺D⁻_sXⁿ⁺¹^/², Uⁿ⁺¹^/²iΓ_h

= −2c_bhD⁺_sD_s⁻Xⁿ⁺¹^/², D_s⁺D_s⁻Uⁿ⁺¹^/²iΓh. (4.22)

Note that the last identity is obtained by applying hD_s⁺φ, ψi_Γ_h= −hφ, D⁻_sψi_Γ_h twice again. The discrete bending energy is

Bⁿ⁺¹− Bⁿ =c_b

2hD⁺_sD⁻_sXⁿ⁺¹, D_s⁺D⁻_sXⁿ⁺¹iΓ_h−c_b

2hD_s⁺D⁻_sXⁿ, D_s⁺D_s⁻XⁿiΓ_h

=c_b

2hD⁺_sD⁻_s(Xⁿ⁺¹+ Xⁿ), D_s⁺D⁻_s(Xⁿ⁺¹− Xⁿ)i_Γ_h

= cb∆thD⁺_sD⁻_sXⁿ⁺¹^/², D⁺_sD⁻_sUⁿ⁺¹^/²iΓ_h

= −∆t

2 huⁿ⁺¹+ uⁿ, S_hⁿ⁻¹^/²(Fⁿ⁺_b ¹^/²)iΩ_h (by (4.22)).

Thus the successive difference of total energy reads (4.23) Eⁿ⁺¹− Eⁿ = −µ∆t

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

h,

which shows that the present CN scheme is unconditionally energy stable. Comparing (4.20) and (4.23), the present BE scheme is more energy dissipative than the CN scheme.

5. Numerical method.

5.1. Linear solver. The resultant linear system for the BE scheme (equations (4.2)–(4.8)) is expressed as

(5.1)







ρ

∆tI − µ∆h

∇h −Sσ −Sb

∇^T_h 0 0 0

−S_σ^T 0 0 0

−Ih 0 0 _∆t¹ I











 uⁿ⁺¹ pⁿ⁺¹ σⁿ⁺¹ Xⁿ⁺¹







=







ρ

∆tuⁿ 0 0

1

∆tXⁿ





 .

Likewise, the resultant linear system for the CN scheme (equations (4.9)–(4.15)) is expressed as

(10)

(5.2)







ρ

∆tI −^µ₂∆h

∇h −Sσ −¹₂Sb

∇^T_h 0 0 0

−S^T_σ 0 0 0

−¹₂Ih 0 0 _∆t¹ I











 uⁿ⁺¹ pⁿ⁺¹^/² σⁿ⁺¹^/² Xⁿ⁺¹







=







ρ

∆tI +^µ₂∆h uⁿ+¹₂Sb(Xⁿ) 0

S_σ^T(uⁿ)

1

∆tXⁿ+¹₂Ih(uⁿ)





 .

In the staggered grid discretization framework, the discrete divergence operator of the fluid velocity can be written as the transpose of the discrete gradient operator of the pressure. Furthermore, the discrete surface divergence operator of the velocity can be written as the transpose of the discrete spreading operator of the elastic tension [21].

At each time-step, the operator Sσ, S_σ^T, Sb, and Ih are updated, as well as the corresponding source terms. Notice that in the case of no-slip boundary conditions being used, there will be an extra term in the right-hand side arising from the boundary conditions for uⁿ⁺¹.

There are several possible methods for solving the unsymmetric sparse linear systems (5.1) and (5.2). The straightforward strategy is to apply Krylov subspace methods directly, such as GMRES and BICGSTAB. However, it is very likely that these sparse iterative linear solvers may fail to converge without an efficient precon- ditioner. Another one is based on the double Schur complement [9, 7, 12] to which Krylov subspace methods can be applied subsequently. It is realized that this approach also suffers from a lack of efficient preconditioners on top of operator splitting errors. As a consequence, semi-implicit IB schemes with iterative linear solvers may be more computationally expensive than the explicit IB methods. For the same simulation time, they may need more iterations to solve the linear system at a given time-step than the number of time-steps that would be required by the explicit IB schemes, although they grant a larger time-step than explicit IB methods [7, 12].

In the present study, the unsymmetric sparse linear systems (5.1) and (5.2) are solved directly using the unsymmetric multifrontal method [31]. We assign the pressure value at the first grid point to remove the rank deficiency of the system due to the nonuniqueness of the pressure. As far as we are concerned, the option of semi-implicit IB schemes with sparse direct linear solvers has not been explored. The primary advantages of direct solvers are highly robust and accurate. They do not require any preconditioners, and they eliminate the need for any iterations within each time-step, and hence convergence analysis. On the other hand, the major deficiencies of sparse direct solvers are the programming complexity including explicit matrix entries at each time-step, and the numerical LU factorization would lead to a large amount of fill-in which occupies a large amount of computer memory. These shortcomings are outweighed by their advantages over sparse iterative linear solvers for the present unconditionally energy stable schemes to solve the inextensible interface problem with bending. Furthermore, semi-implicit IB schemes with sparse direct linear solvers could serve as an important benchmark for developing potential fast solvers and efficient preconditioners in the future.

5.2. Initial arc-length parametrization. In the present study, our initial interface is always chosen as an ellipse so the initial arc-length parametrization is needed.

Let us consider the ellipse be parameterized by Xk = (Xk, Yk) = (a cos θk, b sin θk), k = 0, 1, . . . M − 1, where a is the semiminor axis and b is the semimajor axis. Hence- forth, the total arc-length L of the ellipse is given by

(11)

L = 4b Z π/2

0

p1 − e²sin²θ dθ = 4bEm(e),

where E_m(e) is the complete elliptic integral of the second kind and e =p1 − (a/b)² is the eccentricity. In the current implementation, E_m(e) is computed numerically using the piecewise minimax rational function approximation [32]. Since an ellipse is fourfold symmetry, we first compute the position vector (Xk, Yk) of discrete points sk with uniform arc-length on the first quadrant. Then the discrete points sk on the other quadrants can be obtained subsequently. Let

qk = kEm(e)

M/4 , k = 0, 1, . . . M/4 − 1,

the corresponding θ_k is therefore given by θ_k = E_m⁻¹(q_k, e). The inversion of the incomplete elliptic integral of the second kind E_m⁻¹(qk, e) with respect to θk can be evaluated by solving the transcendental equation Em(θk, e)−qk= 0 with the Newton–

Raphson method [33], i.e.,

(5.3) θ^l+1_k = θ^l_k− E_m(θ^l_k, e) − q_k q

1 − e²sin²θ^l_k ,

where l is the iteration index and Em(θ, e) is the incomplete elliptic integral of the second kind which is evaluated numerically by the half and double argument trans- formations [34]. Practically, the solution of (5.3) converges to θ_k= E_m⁻¹(q_k, e) within an error tolerance of 10⁻¹² in five iterations. The overall procedure only requires it to be carried out at the initialization stage before the time-marching simulation, and it only imposes a negligible computation overhead.

5.3. Nonuniform fourth derivative operator. The exact inextensibility constraint (2.6) requires the initial uniform arc-length mesh to be always fixed as time proceeds. However, in practical simulations, the immersed interface is observed to elongate or shrink sluggishly, i.e., the distribution of discrete points sk with uniform arc-length will not be sustainable. Consequently, the uniform fourth derivative operator D⁺_sD⁻_sD⁺_sD_s⁻φk used in (4.5) for the BE scheme and in (4.12) for the CN scheme should be modified accordingly with respect to local discrete arc-length ∆sk.

Here, we define the nonuniform fourth derivative operator based on Taylor series expansion,

∂⁴φ_k

∂s⁴ ≈

2

X

p=−2

c_k+pφ_k+p

= ck−2φk−2+ ck−1φk−1+ ckφk+ ck+1φk+1+ ck+2φk+2, (5.4)

where the stencil coefficients are given by

ck−2= 24

∆sk−2(∆sk−2+∆sk−1)(∆sk−2+∆sk−1+∆sk)(∆sk−2+∆sk−1+∆sk+∆sk+1),

c_k−1= − 24

∆sk−2∆sk−1(∆sk−1+ ∆sk)(∆sk−1+ ∆sk+ ∆sk+1),

ck= 24

∆s_k−1∆s_k(∆s_k−2+ ∆s_k−1)(∆s_k+ ∆s_k+1),

ck+1= − 24

∆sk∆sk+1(∆sk−1+ ∆sk)(∆sk−2+ ∆sk−1+ ∆sk),

c_k+2= 24

.

(12)

At each time-step, the local arc-length ∆skis evaluated and the stencil coefficients are updated. Note that (5.4) is equivalent to D⁺_sD⁻_sD⁺_sD_s⁻φk when the local arc-length

∆sk = ∆s is uniform. More precisely, for the uniform arc-length mesh, (5.4) becomes (5.5) ∂⁴φk

∂s⁴ ≈ φk−2− 4φk−1+ 6φk− 4φk+1+ φk+2

∆s⁴ ≡ D⁺_sD_s⁻D_s⁺D⁻_sφk.

Since the interface is closed, the interfacial markers X_k, k = 0, 1, . . . M − 1, are periodic. Thus, the periodic boundary condition can be applied to the above difference formula as φk = φ_k(mod)M. For instance, at k = 0, the fourth-order derivative is calculated as

∂⁴φ0

∂s⁴ ≈ cM −2φM −2+ cM −1φM −1+ c0φ0+ c1φ1+ c2φ2. Similar treatments are also applied to the points at k = 1, M − 2, M − 1.

5.4. Discrete delta function. The two-dimensional discrete delta function is represented by a tensor product of a kernel ψ(r),

(5.6) δ_h(x) = 1

h²ψx h

ψy h

. We employ the following kernel ψ(r) [35]:

(5.7)

ψ(r) =











3

8+₃₂^π −^r₄², |r| < 0.5,

1

4+^1−|r|₈ p−2 + 8|r| − 4r²−¹₈sin⁻¹ √

2 (|r| − 1), 0.5 ≤ |r| < 1.5,

17

16−₆₄^π −^3|r|₄ +^r₈² +^|r|−2₁₆ p−14 + 16|r| − 4r² +₁₆¹ sin⁻¹ √

2 (|r| − 2), 1.5 ≤ |r| ≤ 2.5,

0, 2.5 ≤ |r|.

Equation (5.7) holds the advantage of suppressing numerical wiggles in the vicinity of an interface caused by the IB method. We remark that various discrete delta functions have also been utilized in the IB method, and the development of the discrete delta function is an area of active research [36].

6. Numerical results.

6.1. Energy stability check. We first conduct a simple energy stability check for the present BE and CN IB schemes developed in section 4. We set an initial vesicle shape as an ellipse X(θ) = (0.2 cos θ, 0.5 sin θ) which is immersed in a quiescent fluid domain Ω = [0, 2] × [0, 2]. For other parameters, we set ρ = 1, µ = 1, and cb = 0.01.

The simulations are time-marched from time t = 0 to t = 3.0 with a fixed mesh size of h = 1/128 but differ on time-step sizes ∆t = 2h, h, h/2, h² to check the numerical stability. We choose the number of Lagrangian markers M such that ∆s is proportional to h. One should be aware that the choice of the ratio ∆s/h is not exhaustive, and it is usually a localized problem-dependent parameter. Nevertheless, the current choice is realized to be one of the optimum ratios that are achieving the desired accuracy without compromising both the area conservation and total arc- length conservation, as well as alleviating the potential numerical wiggles, particularly on the solution of pressure and elastic tension in the vicinity of the interface.

Figure 1 shows the transient dissipation profiles of total energy E normalized by initial total energy E0at t = 0. Both BE and CN produce consistent and nearly identi- cal results. The total energy profiles are monotonically decreasing, which substantiates that present BE and CN schemes are unconditionally energy stable. Figure 2 depicts