arXiv:0801.2398v2 [] 18 Jan 2008

Removing the Stiffness of Elastic Force from the Immersed
Boundary Method for the 2D Stokes Equations
Thomas Y. Hou∗

Zuoqiang Shi

†

June 16, 2007

Abstract
The Immersed Boundary method has evolved into one of the most useful computational methods in studying fluid structure interaction. On the other hand, the Immersed
Boundary method is also known to suffer from a severe timestep stability restriction when
using an explicit time discretization. In this paper, we propose several efficient semiimplicit schemes to remove this stiffness from the Immersed Boundary method for the
two-dimensional Stokes flow. First, we obtain a novel unconditionally stable semi-implicit
discretization for the immersed boundary problem. Using this unconditionally stable discretization as a building block, we derive several efficient semi-implicit schemes for the
immersed boundary problem by applying the Small Scale Decomposition to this unconditionally stable discretization. Our stability analysis and extensive numerical experiments
show that our semi-implicit schemes offer much better stability property than the explicit
scheme. Unlike other implicit or semi-implicit schemes proposed in the literature, our semiimplicit schemes can be solved explicitly in the spectral space. Thus the computational
cost of our semi-implicit schemes is comparable to that of an explicit scheme, but with a
much better stability property.

1

Introduction

The Immersed Boundary method was originally introduced by Peskin in the 1970’s to model
the flow around heart valves. Now it has evolved into a general useful method in studying the
motion of one or more massless, elastic surface immersed in an incompressible, viscous fluid,
particularly in biofluid dynamics problems where complex geometries and immersed elastic
membranes are present. The method has been successfully applied to a variety of problems
including blood flow in the heart [25, 16, 17, 18, 26, 19, 20], vibrations of the cochlear basilar
membrane [2, 8], platelet aggregation during clotting [7, 34], aquatic locomotion [5, 6, 11, 35, 3],
flow with suspended particles [5, 31], and inset flight [21, 22], We refer to [27] for an extensive
list of applications.
The Immersed Boundary method employs a uniform Eulerian grid over the entire domain
to describe the velocity field of the fluid and a Lagrangian description for the immersed elastic
∗

Applied and Comput. Math, Caltech, Pasadena, CA 91125. Email: hou@acm.caltech.edu.
Applied and Comput. Math, Caltech, Pasadena, CA 91125 and Zhou Pei-Yuan Center for Applied Mathematics, Tsinghua University, Beijing 100084, China. Email: shi@acm.caltech.edu.
†

1

structure. The force generated by the elastic structure drives the fluid and the fluid moves
the elastic structure. This interaction is expressed in terms of the spreading and interpolation
operations by use of smoothed Delta functions.
One of the main difficulties that the Immersed Boundary method encounters is that it
suffers from a severe timestep restriction in order to keep the stability [27, 32, 30]. This has
been the major limitation of the Immersed Boundary method. This restriction is typically
much more severe than the one that would be imposed from using an explicit discretization for
the convection term in the Navier-Stokes equation. The instability is known to arise from large
boundary force and small viscosity [32]. Much effort has been made to remove this restriction.
Some implicit and semi-implicit methods have been proposed in the literature [33, 23, 15].
Despite of these efforts, the timestep restriction has not been resolved satisfactorily. The
computational cost of using an implicit or semi-implicit scheme is still too high to be effective
in a practical computation. To date, almost all practical computations using the immersed
boundary method have been performed using an explicit discretization.
In this paper, we develop several efficient semi-implicit schemes to compute the motion
of an elastic interface immersed in a two-dimensional, incompressible Stokes flow. There are
several important ingredients in deriving our semi-implicit schemes. The first one is to use the
arclength and tangent angle formulation to describe the dynamics of the immersed interface
[9]. We remark that Ceniceros and Roma have also used the arclength and tangent angle
formulation to alleviate the stiffness of the viscous vortex sheet with surface tension in [4]. The
second one is to obtain an unconditionally stable semi-implicit discretization of the immersed
boundary problem. Throughout this paper, we use the term “stability” to mean that the energy
norm of the solution can be bounded in terms of the energy norm of the initial data, which is
a weaker result than proving that the difference between two solutions in the energy norm can
be bounded in terms of the energy norm of their difference at time zero. The third ingredient
is to perform Small Scale Decomposition to the unconditionally stable discretization to obtain
our efficient semi-implicit schemes. An important feature of our small scale decomposition
is that the leading order term, which is to be discretized implicitly, can be expressed as a
convolution operator. This property enables us to solve for the implicit solution explicitly
using the Fourier transformation. Thus, the computational cost of our semi-implicit schemes
is comparable to that of an explicit method. This offers a significant computational saving in
using the Immersed Boundary method.
The Small Scale Decomposition was first developed by Hou, Lowengrub and Shelley [9, 10].
They applied this method to remove the stiffness from interfacial flow with surface tension,
which has proved to be very successful. Due to the coupling between the elastic boundary
with the fluid, it is more difficult to remove the stiffness induced by the elastic force in the
Immersed Boundary method. To remove the stiffness in the Immersed Boundary method,
we need to decouple the stiffness induced by the elastic force from the fluid flow in such a
way that the resulting semi-implicit discretization is still unconditionally stable. This is accomplished by using a semi-implicit discretization which preserves certain important solution
structures which exist at the continuous level. Without obtaining this unconditionally stable
semi-implicit discretization, a straightforward application of the Small Scale Decomposition
to the Immersed Boundary method would not provide an efficient semi-implicit scheme with
the desirable stability property. Very recently, Newren et al. have obtained an unconditionally stable discretization for linear force in [24]. However, they did not perform Small Scale
2

Decomposition to their unconditionally stable discretization. As we will demonstrate in this
paper, the unconditionally stable semi-implicit discretization without using the Small Scale
Decomposition is still very expensive and the gain over the explicit discretization is quite
limited.
We develop several efficient semi-implicit schemes for both the steady Stokes flow and the
unsteady Stokes flow respectively. In both cases, our semi-implicit schemes work very well.
In the steady Stokes flow, we also develop a fourth order semi-implicit scheme by using the
integral factor method. For the unsteady Stokes flow, we develop a second order semi-implicit
method by combining our Small Scale Decomposition with a well known second order temporal
discretization [13, 27]. To illustrate the stability properties of our semi-implicit schemes, we
apply our methods to several prototype problems and test our schemes for a range of elastic
coefficients and viscosity coefficients. Our numerical results confirm that the semi-implicit
schemes remove the high order stability constraint induced by the elastic force. In the case
of unsteady Stokes equation, we also confirm the second order accuracy of our semi-implicit
scheme.
This paper is organized as follows. First, we review the classical formulation of the Immersed Boundary method in Section 2. Then, we introduce the arclength and tangent angle
formulation in Section 3. In Section 4, we describe the spatial discretization of the Immersed
Boundary method. In Section 5-6, we develop the numerical schemes for steady Stokes flow
and unsteady Stokes flow respectively. The numerical results are presented in Section 7. Our
numerical studies will focus on the stability restriction and computational cost of our methods.
Some concluding remarks are given in Section 8.

2

Review of the Immersed Boundary method

For simplicity, we just consider a viscous incompressible fluid in a two dimensional domain
Ω, containing an immersed massless elastic boundary in the form of a closed simple curve
Γ. The configuration of the boundary is given in a parametric form: X(α, t), 0 ≤ α ≤ Lb ,
X(0, t) = X(Lb , t), α tracks a material point of the boundary. We consider only the Stokes
equations in this paper and would neglect the convection term. Then the governing equations
are given as follows:
∂u
= −∇p + µ△u + f (x, t) ,
∂t
∇·u = 0,
ρ

(1)
(2)

∂X
(α, t) = u(X(α, t), t) ,
(3)
∂t
where u is the fluid velocity, p is the pressure, ρ and µ are constant fluid density and viscosity,
f (x, t) is the force density, which is not zero only on the boundary and which is infinite there.
The force density can be expressed as below
f (x, t) =

Z

Lb
0

F(α, t)δ(x − X(α, t))dα,

δ denotes the two-dimensional Dirac delta function and
∂
F(α, t) =
(T τ ),
∂α
3

(4)

(5)



T =T

∂X
∂α



.

(6)

The choice of function T in this paper is computed by Hook’s law




T = Sb

∂X
−1 ,
∂α

(7)

where Sb is the elastic coefficient of the boundary, and τ is the unit tangent vector along the
boundary, which is defined as
∂X
τ =
∂s

,

∂X
.
∂s

(8)

This choice of force density has been used widely in the literature in both computational and
theoretical studies [12],[29],[33].
We can rewrite (3) in the following way:
∂X
(α, t) =
∂t

Z

Ω

u(x, t)δ(x − X(α, t))dx.

(9)

Next, we introduce the spreading and interpolation operations. The spreading and interpolation operators are defined as follows:
L(X)(g(α))(x) =
L∗ (X)(u(x))(α) =

Z

ZΓ

Ω

g(α)δ(x − X(α, t))dα,

(10)

u(x)δ(x − X(α, t))dx .

(11)

It is easy to show that L and L∗ are adjoint operators:
< u(x), L(X)(g(α)) >Ω
=
=
=
=

Z

u(x)

ZΩ Z

ZΩ Z Γ

ZΓ
Γ

Ω

Z

Γ



g(α)δ(x − X(α, t))dα dx

u(x)g(α)δ(x − X(α, t))dαdx
u(x)g(α)δ(x − X(α, t))dxdα

g(α)

Z

Ω



u(x)δ(x − X(α, t))dx dα

= < L∗ (X)(u(x)), g(α) >Γ ,

(12)

where the inner product are defined as follows:
< u, v >Ω =
< f, g >Γ =

Z

ZΩ

u(x)v(x)dx,

(13)

f (α)g(α)dα.

(14)

Γ

Equations (1),(2) are the familiar Stokes equations of viscous incompressible fluid. Equations (3),(4) represent the interaction of the fluid and the elastic boundary. The elastic boundary applies the force to the fluid, the fluid carries the immersed boundary, and the force density
is determined by the configuration of the boundary.
4

3

The arclength-tangent angle formulation

In studying the evolution of a curve, it is useful to represent the curve by its tangent angle θ
and local arclength derivative sα . Previously, Hou, Lowengrub and Shelley [9] exploited this
formulation and combined it with a so-called ”Small Scale Decomposition” reformulation to
remove the stiffness induced by surface tension.
Consider the evolution of a simply closed curve Γ with known normal and tangent velocity
fields, U, V . Assume the curve is represented by X(α, t), α ∈ [0, Lb ]. We define the arclength
derivative, sα , and the tangent vector, θ, as follows
sα (α, t) = |Xα (α, t)|,

(15)

(xα (α, t), yα (α, t)) = sα (α, t)(cos θ(α, t), sin θ(α, t)).

(16)

The closed curve Γ evolves according to
∂X
= u(X, t) = U n + V τ ,
∂t

(17)

where τ and n are the unit tangent and normal vectors of the curve respectively. According
τ = kn, ∂n = −kτ , here s is the arclength variable. It is easy
to the Frenet formula, we have ∂∂s
∂s
to see that sα and θ satisfy the following evolution equations [9]:
(sα )t = Vα − θα U,
Uα V θα
+
.
θt =
sα
sα

(18)
(19)

Given sα and θ, the curve Γ can be reconstructed up to a translation by integrating (16).
However, we also need a point on the boundary to provide the constant of integration.
Using the sα − θ formulation, we can reformulate the immersed boundary problem as
follows:
∂u
= −∇p + µ△u + L(X) (F(sα , θ)) ,
∂t
∇·u = 0
ρ

U

V

(20)
(21)

∗

(22)

∗

(23)

= L (X)(u(x)) · n,

= L (X)(u(x)) · τ ,

(sα )t = Vα − θα U,
Uα V θα
+
,
θt =
sα
sα

(24)
(25)

where
F(sα , θ) =

4

∂
(T τ ) = Sb (sα,α τ + (sα − 1)θα n) .
∂α

(26)

Spatial Discretization

We use the spectral method to discretize both the Stokes equations and the immersed boundary
equations in space since we work on periodic domains. We first discuss the discretization of the
5

Stokes equations in a regular N × N Cartesian grid with a uniform meshsize h. Let xj = jh
and yj = jh. The discrete Fourier transform and inverse Fourier transform are defined as
follows:
Fh,x (φ)(k, y) =
Fh,y (φ)(x, k) =

−1
1 NX
b y),
φ(xj , y)e−ikxj = φ(k,
N j=0
−1
1 NX
b k),
φ(x, yj )e−ikyj = φ(x,
N j=0

−N/2 + 1 ≤ k ≤ N/2,

(27)

−N/2 + 1 ≤ k ≤ N/2,

(28)

N/2

−1 b
(φ)(xj , y)
Fh,x

=

X

k=−N/2+1
N/2

−1 b
(φ)(x, yj ) =
Fh,y

X

k=−N/2+1

b y)eikxj = φ(xj , y),
φ(k,

0 ≤ j ≤ N − 1,

(29)

b k)eikyj = φ(x, y ),
φ(x,
j

0 ≤ j ≤ N − 1.

(30)

Now we introduce the discrete differential operator using the discrete Fourier transform defined above. For a function φ(x, y) defined in the fluid domain Ω, we approximate its spatial
derivatives as follows:
−1
(ik (Fh,x φ) (k, y)) ,
(Dh,x φ) (x, y) = Fh,x

(Dh,y φ) (x, y) =

(31)

−1
(ik (Fh,y φ) (x, k)) .
Fh,y

(32)

Denote ∇h = (Dh,x , Dh,y ). The differential operators are discretized in terms of Dh :
∇p → ∇h p,

∇ · u → ∇h · u,
2

∇ u → ∇h · ∇h u ≡

(33)
(34)
∇2h u.

(35)

Next, we describe the discretization of the immersed boundary. We employ a Lagrangian
grid with grid space ∆α. The number of grid points along the boundary is Nb . For a function
ψ(α) defined on the interface Γ, we define the discrete Fourier transform and its inverse as
follows:
F∆α (ψ)(k) =

Nb −1
1 X
b
φ(αj )e−ikαj = ψ(k),
Nb j=0

αj = j∆α,

(36)

Nb

−1 b
(ψ)(αj ) =
F∆α

2
X

N
k=− 2b +1

ikαj
b
ψ(k)e
= ψ(αj ).

(37)

When the interface is a closed curve, we can approximate the derivative operator along the
interface by the spectral derivative:
−1
(ik (F∆α φ) (k)) .
(D∆α ψ) (α) = F∆α

(38)

When the solution is not periodic, we can also use a finite difference method to discretize the
derivative, we refer to [27] for more details.
6

Now we discuss the discretization of the spreading and interpolation operators. These two
operators both involve the use of a discrete delta function. The discrete delta function we use
is introduced by Peskin in [27]:
 

and
φ(r) =









1
8

1
8




3 − 2|r| +

5 − 2|r| −

 

1
x
y
φ
φ
,
h2
h
h

δh (x, y) =



p

p

1 + 4|r| − 4r 2 ,



(39)

|r| ≤ 1,

−7 + 12|r| − 4r 2 , 1 ≤ |r| ≤ 2,
0,
|r| > 2.

(40)

Using the above discrete delta function, we can discretize the spreading and interpolation
operator as follows
X

Lh (X)(g(α))(x) =

α∈GΓ

L∗h (X)(u(x))(α)

X

=

x∈GΩ

g(α)δh (x − X(α, t))∆α,

(41)

u(x)δh (x − X(α, t))h2 .

(42)

The summation above is over grid points in Γ in (41) and over grid points in Ω in (42). Operator
Lh and L∗h are still adjoint using the following discrete inner product:
X

< f, g >Γh =

f (α)g(α)∆α,

(43)

α∈GΓ

< u, v >Ωh =

X

u(x)v(x)h2 .

(44)

x∈GΩ

Using the inner product defined above, we have:
< u(x), L(X)(g(α)) >Ωh
=

X

u(x)L(X)(g(α))h2

x∈GΩ

=

X

x∈GΩ

=

X

x∈GΓ

u(x)h2

X

α∈GΓ

g(α)∆α

g(α)δh (x − X(α, t))∆α

X

α∈GΩ

u(x)δh (x − X(α, t))h2

= < L∗h (X)(u(x)), g(α) >Γh .

(45)

As we will see later, this discrete self-adjoint property is crucial in obtaining our unconditional
stable semi-discrete scheme for the immersed boundary problem.

5
5.1

Steady Stokes flow
Formulation

For simplicity, we study the steady Stokes flow first. The governing equations for the steady
Stokes flow are given as follows:
0 = −∇p + µ△u + L(X) (F(sα , θ)) ,
7

(46)

∇ · u = 0,
U

(47)

= u(X(α, t), t) · n,

V

(48)

= u(X(α, t), t) · τ ,

(49)

(sα )t = Vα − θα U,
Uα V θα
θt =
+
.
sα
sα

(50)
(51)

In this simple case, we can use a boundary integral method for the two dimension Stokes
flow (page 60 of [28]) to solve equations (46)-(47) to get the velocity on the boundary:
u(X(α, t)) =
v(X(α, t)) =

1
4πµ
1
4πµ

Z

Γ

Z

Γ

!

r2
r1 r2
−(ln r)F1 (α , t) + 12 F1 (α′ , t) + 2 F2 (α′ , t) dα′ ,
r
r

(52)

r1 r2
r2
−(ln r)F2 (α , t) + 22 F2 (α′ , t) + 2 F1 (α′ , t) dα′ ,
r
r

(53)

′

!

′

where r = |r| and
r = (r1 , r2 ) = X(α, t) − X(α′ , t),

5.2

F = (F1 , F2 ),

u = (u, v).

(54)

Small Scale Decomposition

As we can see from (52)-(53), the velocity field can be expressed as a singular integral with a
kernel ln(r). However, the singular velocity integral is nonlinear and nonlocal. It is difficult
to solve for the implicit solution if we treat the velocity integral fully implicitly. The main
idea of the Small Scale Decomposition technique introduced in [9] is to decompose the singular
velocity integral into the sum of a linear singular operator which is a convolution operator
independent of time t and the configuration of the curve, and a remainder operator which is
regular. Since the remaining operator, which is nonlinear and nonlocal, is regular, the simplified convolution integral operator captures accurately the high frequency spectral property
of the original velocity integral. Thus, if we treat only the leading order convolution operator
implicitly, but keep the regular remainder operator explicitly, we can effectively remove the
stiffness of the original velocity field which comes mainly from the high frequency modes of
the solution. In this subsection, we will show how to perform such Small Scale Decomposition
for the Immersed Boundary method applied to the Stokes flow.
Observe that in the integral representation of the velocity field, (52)-(53), the only singular
part of the kernel is ln(r). The other part of the kernel is smooth. Thus to the leading order
contribution of the velocity field can be expressed as follows:
1
u(X(α, t)) ∼
4πµ
V

Z

Γ

−(ln r)F(α′ , t)dα′ ,

= u(X(α, t), t) · τ (α)
Z
1
∼
−(ln r)F(α′ , t) · τ (α)dα′
4πµ Γ
Z

Sb
−(ln r) sα,α′ τ (α′ ) + (sα − 1)θα′ n(α′ ) · τ (α)dα′ .
=
4πµ Γ
8

(55)
(56)

(57)

Next, we perform a Taylor expansion for r, τ (α′ ) · τ (α) and n(α′ ) · τ (α) as a function of α′
around α. By keeping only the leading order term, we have
r ∼ sα (α)|α − α′ |,

(58)

′

(59)

′

(60)

τ (α ) · τ (α) ∼ 1,

n(α ) · τ (α) ∼ 0.

Substituting the above Taylor expansions to (56), we get
V

Sb
4πµ

∼

Z



ln sα (α)|α − α′ | sα,α′ dα′ .

Integrating by part, we obtain
V

Sb
4πµ

∼

where H is the Hilbert transform

Z

1
H[f ](α) =
π

(61)

Sb
sα′
dα′ = − H[sα ],
−α
4µ

(62)

α′
Z

+∞

−∞

f (α′ )
dα′ .
α − α′

(63)

Using the same method, we can get the leading order contributions of U and Uα as follows:
U

∼

Uα ∼

Z

Sb
− ln |α − α′ |(sα′ − 1)θα′ dα′
4πµ
Z
Sb
Sb
(sα′ − 1)θα′ ′
dα = − H[(sα − 1)θα ].
′
4πµ
α −α
4µ

(64)
(65)

Note that if f is a smooth function, then the commutator [H, f ]u ≡ H(f u) − f H(u) is a
smoothing operator for u. Thus we can factor a smooth function from the Hilbert transform
without changing its leading order spectral property. Suppose that sα is smooth, then we
obtain to the leading order that
Uα ∼ −

Sb
(sα − 1)H[θα ].
4µ

(66)

Applying the same analysis to the Eqs (115)-(116) gives




Sb
Sb
(sα )t = − H[D∆α sα ] + D∆α V − D∆α θU +
H[D∆α sα ] ,
(67)
4µ
4µ






V D∆α θ
Sb
D∆α U
1
1
Sb
+
+
H[D∆α θ] +
H[D∆α θ] . (68)
1−
1−
θt = −
4µ
sα
sα
sα
4µ
sα
Note that the leading order operator is linear. This suggests a natural semi-implicit discretization of the immersed boundary problem.
Since we are dealing with a closed immersed boundary, it is natural to work in the Fourier
space. Furthermore, the Hilbert operator has a very simple kernel under the Fourier transformation. Notice that θ is not a periodic function of α. Its value increases 2π every time α
increases Lb . Nevertheless, if we let
θ(α, t) =

2π
α + φ(α, t), α ∈ [0, Lb ],
Lb
9

(69)

then φ is periodic. It is more convenient to work with φ than θ. Substituting (69) into (68)
and taking the Fourier transform on both sides of (67),(68), we obtain
ŝα,t
φ̂t






Sb
Sb
= − |k|ŝα + F (D∆α V − D∆α θU ) +
|k|ŝα ,
4µ
4µ
 


Sb
D∆α U
V D∆α θ
Sb
= − γ|k|φ̂ + F
+
γ|k|φ̂ ,
+
4µ
sα
sα
4µ

(70)
(71)



1
b = −i sgn(k) with sgn(k) being
where γ = max 1 −
. We have also used the fact that H
k
α
sα
the signature function. The first term on the right hand side captures the leading order high
frequency contribution of the terms from the right hand side. An important property of this
leading order term is that it is linear in ŝα and θ̂ and has constant coefficient in space. This
provides a straightforward application of the implicit time discretization.
Since our small scale decomposition is exact near the equilibrium, we can use this result to
get the stability constraint of the explicit scheme by using a frozen coefficient analysis. The
stability constraint is given by
∆t < C

µ h
.
Sb γ

(72)

As we can see, the time step needs to be very small if Sb is large and µ is small. For example,
if Sb = 100, and µ = 10−2 , then the stability would require that ∆t ≤ C10−4 h.

5.3

Semi-implicit schemes

Based on the small scale decomposition presented in the previous subsection, we propose two
types of semi-implicit schemes in this section. The first implicit time discretization uses the
backward Euler method to discretize the leading order term while keeping the lower order term
explicit. This gives rise to the following semi-implicit scheme:




ŝn+1
− ŝnα
α
∆t

Sb
Sb
= − |k|ŝn+1
+ F (D∆α V n − D∆α θ n U n ) +
|k|ŝnα ,
α
4µ
4µ

φ̂n+1 − φ̂n
∆t

= −



Sb
γ|k|φ̂n+1 + F
4µ



D∆α U n V n D∆α θ n
+
sn+1
sn+1
α
α



+

(73)


Sb
γ|k|φ̂n .
4µ

(74)

We call the above discretization the semi-implicit method. Near equilibrium, the stability
constraint of this numerical method is ∆t < C(Sb , µ), independent of the meshsize h. Since
the small scale decomposition only captures the leading order contribution from the high
frequency components, this method can not eliminate the effect of Sb and µ completely. The
coefficients Sb and µ can still affect the time stability through the low frequency components
of the solution, which comes from the second term of the right hand side. In order to obtain a
semi-implicit discretization with better stability property, we can incorporate the low frequency
contribution from the second term in our implicit discretization. This scheme can be found in
the appendix A. We call it the semi-implicit scheme of second kind.
The accuracy of the semi-implicit schemes presented above is just first order. In order to
get a high order time discretization, we can use the integral factor method. The integral factor
method factors out the leading order linear term prior to time discretization. They usually
10

provide stable and high order time integration methods for stiff problems. To use the integral
factor method, we rewrite (70),(71) as




∂  ηt 
Sb
|k|t P (ŝα , φ̂),
= exp
e ŝα
∂t
4µ


Sb
∂  ξt 
γ|k|t Q(ŝα , φ̂),
e φ̂ = exp
∂t
4µ

where

Sb
|k|,
4µ

(75)
(76)

Sb
γ|k|,
(77)
4µ
Sb
|k|ŝα ,
(78)
P (ŝα , φ̂) = F (Vα − θα U ) +
4µ


Uα V θ α
Sb
Q(ŝα , φ̂) = F
+
γ|k|φ̂.
(79)
+
sα
sα
4µ
Now it is straightforward to discretize this system to high order. In particular, we can apply
the classical fourth order Runge-Kutta method to discretize the above system to obtain a
fourth order semi-implicit scheme.
We remark that although the fourth order semi-implicit scheme based on the integral factor
approach is much more accurate than the first order semi-implicit discretization, the stability of
the fourth order method is weaker than the first semi-implicit scheme based on the backward
Euler discretization. The fact that the higher order discretization gives a weaker stability
property is a phenomenon which has been observed for almost all time integration methods.
It is not a restriction of our semi-implicit schemes for the immersed boundary problem.
The semi-implicit schemes we describe above only update the θ and sα variables. We also
need to reconstruct the boundary at tn+1 from θ n+1 and sn+1
α . For this purpose, we need
to update a reference point of the boundary. This will be done by using an explicit time
integration method. The simplest one is the forward Euler method:
η =

ξ=

xn+1 (0) = xn (0) + ∆t (V n cos(θ n (0)) − U n sin(θ n (0))) ,
y

n+1

n

n

n

n

n

(0) = y (0) + ∆t (V sin(θ (0)) + U cos(θ (0))) ,

(80)
(81)

where U and V are evaluated at the reference point. A higher order integration method can
be also used. In the explicit update of the reference point, we can use the values of U and V
obtained using the semi-implicit discretization from the previous time steps to extrapolate the
values of U and V in the intermediate time steps in our explicit update of the reference point.
Once we have updated the reference point, we can obtain the configuration of the boundary
(x, y) from (sα , θ) by integrating (16)
x

n+1

(α) = x

n+1

(0) +

y n+1 (α) = y n+1 (0) +

Z

α

Z 0α
0

′
n+1 ′
sn+1
(α ))dα′ ,
α (α ) cos(θ

(82)

′
n+1 ′
sn+1
(α ))dα′ .
α (α ) sin(θ

(83)

We can use more than one reference point, then average them to get the last configuration.
This can improve the stability constraint significantly. Actually, in our computation, we use
two reference points X(0), X(Nb /2), then take the average to determine the position of the
interface at next time step. Since we update only two reference points, the extra cost in
updating the reference point is small compared with the overall computational cost.
11

6

Unsteady Stokes flow

6.1

Formulation

In this section, we will extend the semi-implicit discretization developed for the steady Stokes
flow to the unsteady Stokes flow. The governing equations of the immersed boundary method
for the unsteady Stokes flow are as follows:
∂u
= −∇p + µ△u + L(X) (F(sα , θ)) ,
∂t
∇·u = 0
ρ

U

V

(84)
(85)

= u(X(α, t), t) · n,

(86)

= u(X(α, t), t) · τ ,

(87)

sαt = Vα − θα U,
Uα V θα
+
.
θt =
sα
sα

(88)
(89)

It is much more difficult to solve the fluid velocity u analytically from (84)-(85). As for the
steady Stokes flow, we will first derive an unconditionally stable time discretization which will
be given in next section and then apply the Small Scale Decomposition to the unconditionally
stable time discretization to obtain our efficient semi-implicit schemes.

6.2

An unconditionally stable semi-implicit discretization

In this section, we will describe our unconditionally stable semi-implicit discretization of the
Immersed Boundary method for the incompressible unsteady Stokes equations and prove its
unconditional stability in the sense of total energy is non-increasing.
The unconditionally stable semi-implicit discretization is consisted of two steps. In the first
step, we update sα , u from tn to tn+1 , then we get θ n+1 in the second step.
Step 1: Update of un+1 and sn+1
α .
ρ



un+1 − un
n
n
n
= −∇h pn+1 + µ∇2h un+1 + Lh,n F(sn+1
,
θ
;
τ
,
n
)
,
α
∆t


n
n
n
∇2h pn+1 = ∇h · Lh,n F(sn+1
α ,θ ;τ ,n ) ,

V n+1 = L∗h,n (un+1 ) · τ n
U

n+1

− snα
sn+1
α
∆t

=

L∗h,n (un+1 ) ·

n

n ,

(90)
(91)
(92)
(93)

= D∆α V n+1 − D∆α θ n U n+1 ,

(94)

where τ n = (cos(θ n ), sin(θ n )), nn = (− sin(θ n ), cos(θ n )), Lh,n = Lh (Xn ), L∗h,n = L∗h (Xn ),
∇h and D∆α are discrete derivative operators for the Eulerian grid and the Lagrangian grid
respectively, and




n
n
n
n+1 n
n+1
F(sn+1
− 1)D∆α θ n nn .
α , θ ; τ , n ) = Sb D∆α sα τ + (sα

12

(95)

Step 2: Update of θ n+1 . After we have obtained un+1 , pn+1 and sn+1
α , we update θ at
n+1
t
using the following semi-implicit scheme:
ρ



u n+1 − un
n+1
= −∇h p n+1 + µ∇2h u n+1 + Lh,n F(sn+1
; τ n , nn ) ,
α ,θ
∆t


n+1
n
n
,
θ
;
τ
,
n
)
,
∇2h p n+1 = ∇h · Lh,n F(sn+1
α

V

U
θ n+1

n+1

= L∗h,n (u n+1 ) · τ n

n+1

−
∆t

θn

=

L∗h,n (u n+1 )

=



1
sn+1
α

n

·n ,

D∆α U

n+1

(96)
(97)
(98)
(99)

+ D∆α θ n V

n+1



.

(100)

where




(101)

< f, D∆α g >Γh = − < D∆α f, g >Γh ,

< u, ∇h g >Ωh = − < ∇h · u, g >Ωh .

(102)

n+1
n
n+1
− 1)D∆α θ n+1 nn
F(sn+1
; τ n , nn ) = Sb D∆α sn+1
α ,θ
α τ + (sα

It is important to note that the above discretization is not fully implicit. In fact, both the
spreading and interpolation operators are evaluated at the interface Xn from the previous time
and un+1 , in (90) - (94), we use θ n instead of θ n+1 to
step. Moreover, when solve the sn+1
α
evaluate the force density. This makes our semi-implicit discretization linear with respect to
the implicit solution variables, un+1 , θ n+1 , and sn+1
α . The above semi-implicit discretization
essentially decouples the stiffness induced by the elastic force from the fluid equations. This
enables us to remove the stiffness of the Immersed Boundary method effectively by applying
the Small Scale Decomposition and arclength/tangent angle formulation as was done in [9].
In the following, we will prove that this semi-implicit discretization is unconditionally stable
in the energy norm.
By using a discrete summation by parts, we can show that

First, we define the total energy of the physical system. The total energy includes the
kinetic energy K and the potential energy P , which are defined below:
K =

P

=

N
ρ X
1
uij · uij h2 ,
ρ < u, u >Ωh =
2
2 i,j=1

Nb
1
Sb X
Sb < sα − 1, sα − 1 >Γh =
(sα,j − 1)2 ∆α.
2
2 j=1

(103)

(104)

The total energy is then defined as
E = K + P.

(105)

Below we will prove the unconditional stability of our semi-implicit discretization. To
simplify the presentation, we still denote the discrete spectral derivative D∆α g of a function g
as gα .
13

Taking the discrete inner product defined by (44) of (90) with un+1 + un and using (102),
we obtain
2(K n+1 − K n ) = ρ < un+1 + un , un+1 − un >Ωh

= ρ < −un+1 + un , un+1 − un >Ωh +2ρ < un+1 , un+1 − un >Ωh

= −ρ < un+1 − un , un+1 − un >Ωh





n
+2∆t(< un+1 , −∇h p + µ∇2h un+1 + Lh,n F(sn+1
α , θ ) > Ωh )

= −ρ < un+1 − un , un+1 − un >Ωh −2∆t < un+1 , ∇h p >Ωh




n
+2∆t < un+1 , µ∇2h un+1 >Ωh +2∆t < un+1 , Lh,n F(sn+1
α , θ ) > Ωh

= −ρ < un+1 − un , un+1 − un >Ωh −2∆t < ∇h · un+1 , p >Ωh
−2µ∆t < ∇h un+1 , ∇h un+1 >Ωh




n
+2∆t < L∗h,n un+1 , F(sn+1
α , θ ) >Γh .

(106)

The second term on the right hand side of (106) is zero because the discrete velocity field is
divergence free, i.e. ∇h · un+1 = 0. The fourth term can be rewritten as




n
< L∗h,n un+1 , F(sn+1
α , θ ) >Γh





n
n+1
= < V n+1 τ n + U n+1 nn , Sb sn+1
− 1)θαn nn >Γh
α,α τ + (sα





n+1
= Sb < V n+1 , sn+1
, (sn+1
− 1)θαn >Γh .
α,α >Γh + < U
α

Combining (106) and (107), we can get

(107)

2(K n+1 − K n ) = −ρ < un+1 − un , un+1 − un >Ωh −2µ∆t < ∇h un+1 , ∇h un+1 >Ωh
n+1
+2Sb ∆t(< V n+1 , sn+1
, (sn+1
− 1)θαn >Γh ).
α,α >Γh + < U
α

(108)

+ snα − 2 and
Similarly, by taking the discrete inner product defined by (43) of (94) with sn+1
α
using (102), we get
2(P n+1 − P n ) = Sb < sn+1
+ snα − 2, sn+1
− snα >Γh
α
α

= Sb < −sn+1
+ snα , sn+1
− snα >Γh +2Sb < sn+1
− 1, sn+1
− snα >Γh
α
α
α
α

= −Sb < sn+1
− snα , sn+1
− snα >Γh
α
α

+2Sb ∆t < sn+1
− 1, Vαn+1 − θαn U n+1 >Γh
α

= −Sb < sn+1
− snα , sn+1
− snα >Γh
α
α

n+1
+2Sb ∆t(− < sn+1
>Γh − < sn+1
− 1, θαn U n+1 >Γh ).
α,α , V
α

(109)

Adding (108) to (109), we have
1
E n+1 − E n = − ρ < un+1 − un , un+1 − un >Ωh −µ∆t < ∇h un+1 , ∇h un+1 >Ωh
2
1
− snα , sn+1
− snα >Γh
− Sb < sn+1
α
α
2
≤ 0.
(110)
14

This proves that our semi-implicit discretization is unconditionally stable in the sense that the
total energy is non-increasing.
Remark 1. In our proof presented above, we have used two important properties of our
semi-implicit discretization. The first property is that the discrete spreading and interpolation
operators are adjoint. The second property is that the velocity field satisfies the discrete divergence free condition. It is clear from the above proof that as long as these two properties are
satisfied by our spatial discretization, the corresponding semi-implicit discretization introduced
in the previous subsection is unconditionally stable.
Remark 2. We remark that the proof above is similar in spirit to that of a semi-linear discretization obtained by Newrenn et al in [24]. There is some minor difference between the
unconditionally stable semi-implicit discretization obtained by Newren et al and our unconditionally stable semi-implicit discretization. In the problem considered by Newren et al., the
force is a linear function of the interface. On the other hand, in the problem we consider, the
force is a nonlinear function of the interface (the rest length of the boundary is not zero). By
using the sα − θ formulation, the force is a linear function of sα . By treating θ explicitly, we
obtain a semi-implicit discretization that is linear with respect to sα . Due to the decoupling
between sα and θ, we need to solve two Nb × Nb linear systems instead of one 2Nb × 2Nb linear
system in the semi-implicit discretization obtained by Newren et al.
Remark 3. For the steady Stokes flow, we can also prove the following semi-implicit discretization is unconditionally stable:
Step 1:
n
0 = −∇h pn+1 + µ∇2h un+1 + Lh,n (F(sn+1
α , θ )),

∇2h pn+1
n+1

=

n+1

=

V

U

sn+1
− snα
α
∆t
n+1
θ
− θn
∆t

=

= D∆α V n+1 − D∆α θ n U n+1 ,
=

(111)

n
n
n
∇h · Lh,n (F(sn+1
α , θ ; τ , n ))),
L∗h,n (un+1 (x)) · τ n ,
L∗h,n (un+1 (x)) · nn ,

1
sn+1
α



(112)
(113)
(114)
(115)



D∆α U n+1 + D∆α θ n+1 V n+1 .

(116)

Step 2:




n+1
0 = −∇h p n+1 + µ∇2h un+1 + Lh,n F(sn+1
; τ n , nn ) ,
α ,θ





n+1
; τ n , nn ) ,
∇2h p n+1 = ∇h · Lh,n F(sn+1
α ,θ

V
U
θ n+1

n+1

n+1

−
∆t

θn

= L∗h,n (u n+1 ) · τ n

= L∗h,n (u n+1 ) · nn ,
=

1

sn+1
α



D∆α U

n+1

(117)
(118)
(119)
(120)

+ D∆α θ n V

15

n+1



.

(121)

In this case the total energy is just the potential energy
E=

Nb
Sb
Sb X
< sα − 1, sα − 1 >Γh =
(sα,j − 1)2 ∆α.
2
2 j=1

(122)

As in the case of the unsteady Stokes flow, as long as the velocity field satisfies the discrete
divergence free condition and the discrete spreading and interpolation operators are adjoint,
we can prove that above semi-implicit discretization is unconditionally stable in the sense of
total energy is non-increasing.

6.3

Small Scale Decomposition

In order to apply the Small Scale Decomposition to our unconditionally stable time discretization, we would like to solve for the velocity field at time tn+1 from the space-continuous version
of (90) and (91) using an integral representation:








µ∆t 2 −1
∆t
n
∇
(1 − ∇(∇2 )−1 ∇·)Ln (F(sn+1
un +
α , θ ))
ρ
ρ

 

∆t
µ∆t 2 −1
n+1 n
n
∇
Ln (F(sα , θ ))
u +
=
1−
ρ
ρ



∆t
µ∆t 2 −1 2 −1 
n
(∇ )
∇∇ · Ln (F(sn+1
∇
−
1−
α , θ )) .
ρ
ρ

un+1 (x) =

1−

To solve for the velocity field at tn+1 , we need to use the following free space fundamental
solutions in two space dimensions which are defined as follows:




µ∆t 2
1−
∇ E1 = δ(x − x′ ),
ρ


µ∆t 2
2
∇ 1−
∇ E2 = δ(x − x′ ).
ρ

(123)
(124)

These two fundamental solutions can be expressed in terms of the modified Bessel function of
the second kind [1]:
E1 =
E2 =

λ2
K0 (λ|x − x′ |),
2π

1
K0 (λ|x − x′ | + ln(|x − x′ |) ,
2π

(125)
(126)

ρ
and K0 is a modified Bessel function of the second kind. By integrating by
µ∆t
part, we can further express the velocity un+1 as

where λ2 =

n+1

u

(x) =

1
2π

Z

Ω

λ2 K0 (λ|x − x′ |)un (x′ )dx′
Z

1 ∆t
n
′
λ2 K0 (λ|x − Xn (α′ )|)F(sn+1
+
α , θ )dα
2π ρ Γn
Z
1 ∆t
n
′
∇x ∇x (K0 (λ|x − x′ |) + ln(|x − x′ |)) · Ln (F(sn+1
−
α , θ ))dx
2π ρ Ω
16

=

1
2π

Z

Ω

λ2 K0 (λ|x − x′ |)un (x′ )dx′
Z

1 ∆t
n
′
+
λ2 K0 (λ|x − Xn (α′ )|)F(sn+1
α , θ )dα
2π ρ Γn
Z
1 ∆t
n
′
−
G(x − Xn (α′ )) · F(sn+1
α , θ )dα ,
2π ρ Ω
where G is defined as follows:
2ri rj
1
ri rj
δij
−
+ λ2 (K0 (λ|r|) + K2 (λ|r|)) 2
Gij (r) =
2
4
|r|
|r|
2
|r|


ri rj
δij
− 3 ,
−λK1 (λ|r|)
|r|
|r|

(127)

(128)

and K0 , K1 , K2 are all modified Bessel functions of the second kind [1].
In this subsection, we will perform a Small Scale Decomposition to the velocity field based
on the integral representation (127). Recall that in our semi-implicit discretization, the velocity
field at tn+1 is evaluated on the boundary Xn at tn . Thus we should perform our Small Scale
Decomposition for un+1 (Xn ). To this end, we first write down the integral expression of
un+1 (Xn ) as follows:
un+1 (Xn (α)) =

1
2π

Z

Ω

λ2 K0 (λ|Xn (α) − x′ |)un (x′ )dx′
Z

1 ∆t
n
′
λ2 K0 (λ|Xn (α) − Xn (α′ )|)F(sn+1
α , θ )dα
2π ρ Γn
Z
1 ∆t
n
′
G(Xn (α) − Xn (α′ )) · F(sn+1
−
α , θ )dα .
2π ρ Ω
+

(129)

To perform the Small Scale Decomposition to the above velocity integral, we would like to
decompose the singular velocity kernel as the sum of a linear singular operator of convolution
type and a remainder operator which is regular. Using the Taylor expansion for α′ around α,
we get the following decomposition:
V n+1 (α) = un+1 (Xn (α)) · τ n (α)
Z
Sb ∆t
′
∼
λ2 K0 (λsnα |α − α′ |)sn+1
α,α′ dα −
2πρ Γn
!
Z
Sb ∆t
1
1 2
n
′
n
′
λ (K0 (λsα |α − α |) + K2 (λsα |α − α |)) −
2πρ Γn 2
(snα )2 (α − α′ )2
′
sn+1
α,α′ dα ,

(130)

where snα | inside K0 (λsnα |α − α′ |) is evaluated at α. Notice that [1]
d2
dα′2

!


1
1
K0 (λsnα |α − α′ |) = λ2 K0 (λsnα |α − α′ |) + K2 (λsnα |α − α′ |) .
2
2
(snα (α))

(131)

Integrating the right hand side of (130) by parts twice, we get
V n+1 (α) ∼

Z

1
′
Sb ∆t
λ2 K0 (λsnα |α − α′ |)sn+1
α,α′ dα −
2πρ
Γn
Z

Sb ∆t
′
K0 (λsnα |α − α′ |) − ln(α − α′ ) sn+1
α,α′ α′ α′ dα .
2
2πρ (snα ) Γn
17

(132)

Similarly, we can obtain the leading order contribution of U

n+1

as follows:

n+1

(α) = u n+1 (Xn (α)) · nn (α)
Z


Sb ∆t
n+1
n
′
′ 
n+1
dα′ .
K
(λs
|α
−
α
|)
−
ln(α
−
α
)
(s
−
1)θ
∼
′
0
α
α
α
α′ α′
2πρ (snα )2 Γn

U

(133)

Using this decomposition, we obtain the following scheme:



− snα
sn+1
α
∗,n+1
n ∗,n+1
n
= T (sn+1
)
+
D
V
−
D
θ
U
−
T
(s
)
,
∆α
∆α
α
α
∆t
un+1 − un
n
ρ
= −∇h pn+1 + µ∇2h un+1 + Lh,n (F(sn+1
α , θ )),
∆t
n
∇2h pn+1 = ∇h · Lh,n (F(sn+1
α , θ )),

V

n+1

L∗h,n (un+1 ) · τ n ,
L∗h,n (un+1 ) · nn ,

1
S(θ n+1 )

=

U n+1 =
θ n+1

−
∆t

where
T (sn+1
α )



=

−
S(θ

n+1

θn

=

sn+1
α

1
Sb ∆t
2πρ

λ

Γn

Sb ∆t
2πρ (snα )2
Sb ∆t
2πρ (snα )2

) =

Z

Z

Z

2

Γn

Γn

+

(135)
(136)
(137)
(138)

sn+1
α

K0 (λsnα |α







D∆α U n+1 + D∆α θ n V n+1 − S(θ n ) ,

′

−α

′
|)sn+1
α,α′ dα



α



′
K0 (λsnα |α − α′ |) − ln(α − α′ ) sn+1
α,α′ α′ α′ dα

K0 (λsnα |α

(134)

′

′



− α |) − ln(α − α )

(sn+1
α

!

− 1)θα′

(139)

,
α



α′ α′

′

dα

!

,
α

and u∗,n+1 is the velocity at tn+1 which is calculated explicitly
ρ

u∗,n+1 − un
= −∇h p∗,n+1 + µ∇2h u∗,n+1 + Lh,n (F(snα , θ n ))
∆t
∇2h p∗,n+1 = ∇h · Lh,n (F(snα , θ n ))
V

∗,n+1

U

∗,n+1

=

=

L∗h,n (u∗,n+1 ) ·
L∗h,n (u∗,n+1 ) ·

τ

(140)
(141)

n

(142)

n

(143)

n .

The derivation of the above semi-implicit scheme is given in Appendix B.
However, the expressions of T and S are still too complicated and need to be further
simplified. The leading order linear operator, which contains K0 (λsnα (α)|α − α′ |), is not a
convolution operator. Thus, it does not have a simple kernel under the Fourier transform as
the Hilbert operator in the case of the steady Stokes flow. To further simplify the kernel,
we approximate snα (α) by min snα (α). With this approximation, the corresponding leading
α
order operator is a convolution operator and can be diagonalized under the Fourier transform.
Denote β = λ min snα (α). In Appendix C, we will show that
α

F



1
π

Z

+∞

−∞

K0 (β|α − α′ |)f (α′ )dα′
18



fb(k)
.
=p 2
β + k2

(144)

Using (144) and replacing snα (α) by min snα (α), we can simplify the leading order term T (sn+1
α )
α

and S(θ n+1 ) under the Fourier transform:

Tb(sn+1
α ) ∼ −
b n+1 ) ∼ −
S(θ

Sb ∆t



λ min snα
α


2 
 r
n

2ρ min sα
α

Sb ∆t max
α







2

λ min snα
α

k2

+

2

+ k2


−1 
|k|3 − r
2


n

sn+1
α

2ρ min sα
α

k4




n+1
− |k|3 
 ŝα ,

k4

λ min snα
α

2

+ k2

(145)


 n+1
 θ̂
.


(146)

1
≪ 1. By Taylor expanding (145) and (146) with respect to
µ∆t
λ and keeping only the first order term, we obtain the leading order term as follows:
When µ ≫ 1, we have λ = √

Sb
|k|ŝn+1
α ,
4µ


b n+1 ) ∼ − Sb max sn+1 − 1 |k|θ̂ n+1 ,
S(θ
α
4µ α

Tb (sn+1
α ) ∼ −

(147)
(148)

which is the same as the steady Stoke flow. This is also consistent with one’s physical intuition.
When the viscosity is very large, the flow changes very slowly. The inertial term can be
neglected.
1
≫ 1, the asymptotic expansion is
When µ ≪ 1, then λ = √
µ∆t
√
Sb ∆t
n+1
b

T (sα ) ∼ − 
k2 ŝn+1
(149)
α ,
n √
2 min sα
ρµ
α

b n+1 ) ∼ −
S(θ





Sb ∆t max sn+1
−1
α
α

2ρ



min snα
α

2

k3 θ̂ n+1 .

(150)

From the asymptotic expansion above, we can see that our small scale decomposition is also
consistent with the linearized stability analysis which Stockie and Wetton got in [30]. Using
the leading order term above, we can get the leading order term of the eigenvalue same with
the result in [30].
We can also obtain the corresponding stability constraint for the explicit scheme near the
equilibrium:
∆t < C(Sb , µ)hβ ,

(151)

where 1 ≤ β ≤ 3/2. The value of β depends on µ. If µ ≪ 1, then we have β ≈ 3/2. On the
other hand, if µ ≫ 1, we have β ≈ 1.

19

6.4

The numerical scheme

Based on the small scale decomposition we developed in the last subsection, we can now
describe our semi-implicit numerical scheme. Combining the time discretization (90) -(100)
with the decomposition (130)-(133) and using the approximation (145)-(146), we obtain the
following semi-implicit numerical scheme:
Step 1: Update of un+1 and sn+1
α .


− snα
sn+1
α
∗,n+1
n ∗,n+1
n
= T (sn+1
)
+
D
V
−
D
θ
U
−
T
(s
)
,
∆α
∆α
α
α
∆t
un+1 − un
n
ρ
= −∇h pn+1 + µ∇2h un+1 + Lh,n (F(sn+1
α , θ )),
∆t
n
∇2h pn+1 = ∇h · Lh,n (F(sn+1
α , θ )),

(152)
(153)
(154)

where
Tb(sn+1
α ) = −

Sb ∆t





snα
 λ min
α

2 
 r
n

2ρ min sα
α

2

λ min snα
α

k2 + k4
2

+

k2




n+1
− |k|3 
 ŝα ,

(155)
(156)

and u∗,n+1 is the intermediate velocity at tn+1 which is calculated by solving the unsteady
Stokes equations implicitly while evaluating the elastic force explicitly:
ρ

u∗,n+1 − un
= −∇h p∗,n+1 + µ∇2h u∗,n+1 + Lh,n (F(snα , θ n )),
∆t
∇2h p∗,n+1 = ∇h · Lh,n (F(snα , θ n )),
V

∗,n+1

=

U ∗,n+1 =

L∗h,n (u∗,n+1 )
L∗h,n (u∗,n+1 )

(157)
(158)

n

·τ ,

(159)

· nn .

(160)

Step 2: Update of θ n+1 . Once we have updated u, p, and sα at tn+1 , we update θ n+1
using the following semi-implicit scheme:
θ n+1 − θ n
∆t

=

S(θ n+1 )
+
minα sn+1
α



1
sn+1
α



D∆α U

n+1

n

+ D∆α θ V

n+1



n



− S(θ ) ,

(161)

where
V n+1 = L∗h,n (un+1 ) · τ n ,
U

n+1

b n+1 )
S(θ

=

L∗h,n (un+1 )

(162)

n

·n ,





Sb ∆t max (snα − 1) 
 n+1
k4
α
3
 θ̂

r
.
= −
|k|
−

2




n 2
2
2ρ min snα
λ min sα + k
α

(163)
(164)

α

This is our semi-implicit scheme for the unsteady Stokes flow. The spectral discretization
in space has the advantage of being high order accurate and the leading order operator has
20

a simple kernel under the Fourier transform. As it is, the time discretization is only first
order. Based on the first order semi-implicit scheme that we develop in this subsection, we
will develop a second order semi-implicit scheme in the next subsection.
A near equilibrium stability analysis shows that the stability constraint of this semi-implicit
scheme is of the form ∆t < C(Sb , µ), which is independent of the wave number, but still
dependent on Sb and µ. This is due to the fact that the Small Scale Decomposition does
not capture the low frequency components of the solution accurately. The low frequency
components of the solution can affect the stability of the time discretization in two ways.
The first one is through the small scale decomposition, which only captures the leading order
contribution of the solution at high wave numbers. The second one comes from the second
term of the right hand side of the dynamic equations for sα and θ. As in the case of the steady
Stokes flow, we can include the leading order contribution from the second term in our leading
order term and treat them implicitly. This treatment would significantly improve the stability
property especially when the elastic coefficient is large or the viscosity is small. This improved
stability is at the expense of solving a linear system for the implicit solution at each time step.
We call this semi-implicit discretization as the semi-implicit method of the second kind. More
discussions on the semi-implicit method of the second kind can be found in Appendix A.
Remark 4. The leading order term we derive above is calculated analytically using the spacecontinuous formulation with an unsmoothed Dirac delta function. As Stockie and Wetton
pointed out in [32], this analysis over-predicts the stiffness of the Immersed Boundary method
in a practical computation. If we use the leading order approximation directly, the semi-implicit
scheme with the leading order terms derived above tends to over-dissipate the solution. To
alleviate this effect in the practical implementation, we rescale the leading order term by a
coefficient which is calculated at the first time step in the following way:
maxα Vα1,∗
,
maxα T (s0α )
maxα U 1
CU =
,
maxα SU (θ 0 )
CV =

1

where SU (θ 0 ) is the leading order term of U , which can be computed from S(θ 0 ) via the
Fourier transform. The leading order term we use in a practical computation is actually
n+1 ).
CV T (sn+1
α ) and CU S(θ

6.5

A second order semi-implicit scheme

Based on the first order semi-implicit scheme we have developed in the previous subsection,
we will derive the corresponding second order semi-implicit scheme in this subsection.
First, we need to use a robust implicit second order temporal discretization. To simplify
the presentation, we will only describe the semi-discrete algorithm. The space discretization
is done in the same way as before. The second order temporal discretization we use consists
1
of two steps. In the first step, we take a fractional time step from tn to tn+ 2 . It is same with
the first order semi-implicit discretization (90)-(100) , except the timestep is ∆t
2 .
n
In the second step, we integrate the unsteady Stokes equations from t to tn+1 based on
the midpoint and the trapezoidal rules:
21

Step 1: Update of un+1 and sn+1
α .
ρ


 
1
1
1
un+1 − un
= −∇p̄ + µ∇2 ū + Ln+ 1 F s̄α , θ n+ 2 ; τ n+ 2 , nn+ 2 ,
2
∆t

 
1
1
n+ 21
2
; τ n+ 2 , nn+ 2 ,
∇ p̄ = ∇ · Ln+ 1 F s̄α , θ

(165)
(166)

2

1

= L∗n+ 1 (ū) · τ n+ 2 ,

V̄

(167)

2

1

= L∗n+ 1 (ū) · nn+ 2 ,

Ū

(168)

2

sn+1
− snα
α
∆t/2
1

n+ 21

= V̄α − θα
1

Ū ,

(169)

1

1

1

1

where τ n+ 2 = (cos(θ n+ 2 ), sin(θ n+ 2 )), nn+ 2 = (− sin(θ n+ 2 ), cos(θ n+ 2 )), Ln+ 1 = L(Xn +
1
∗
2 ), Ln+ 1

2

2

= L∗ (Xn + 21 ) and
1

1



1

1

1

1

F(s̄α , θ n+ 2 ; τ n+ 2 , nn+ 2 ) = Sb D∆α s̄α τ n+ 2 + (s̄α − 1)D∆α θ n+ 2 nn+ 2



(170)

n+1 using
Step 2: Update of θ n+1 . After we have obtained un+1 and sn+1
α , we update θ at t
the following semi-implicit scheme:

e n+1 − un
u
ρ
∆t

= −∇pe n+1 + µ∇2


e n+1 + u n+1
u
2
1

1

!

2



∇2 pe n+1 = ∇ · Ln+ 1 F(s̄α , θ̄; τ n+ 2 , nn+ 2 ) ,
2

Ve

e
U

θ n+1 − θ n
∆t

e n+1 + u n+1
u
2

= L∗n+ 1
2

=

L∗n+ 1
2
1

=

n+ 21

sα

where



e n+1 + u n+1
u
2
n+ 21

eα + θα
U
1



!
!

1

· τ n+ 2



1

(173)

1

· nn+ 2 ,

(174)
(175)



1

1

F(s̄α , θ̄; τ n+ 2 , nn+ 2 ) = Sb s̄α,α τ n+ 2 + (s̄α − 1)θ̄α nn+ 2
and
ū =

un+1 + un
,
2

s̄α =



(172)

Ve .

1

1

+ Ln+ 1 F(s̄α , θ̄; τ n+ 2 , nn+ 2 ) , (171)

sn+1
+ snα
α
,
2

θ̄ =

θ n+1 + θ n
.
2



(176)

(177)
1

Here, Ln+ 1 and L∗n+ 1 are the spreading and the interpolation operators evaluated at Xn+ 2 .
2

2

Using the same method of analysis, we can prove that the above second order semi-implicit
discretization is unconditionally stable in the sense that the total energy is non-increasing.
The first step in our second order method is identical to the first order method except
that the time step is ∆t
2 instead of ∆t. Thus we can use the first order semi-implicit scheme
introduced in last subsection to compute it directly.
22

In the second step, we can also apply the Small Scale Decomposition with some modifications. After applying the Small Scale Decomposition to the second step of the two-step
method, the second step of the semi-implicit scheme has the form
Step 1: Update un+1 and sn+1
α .
− snα
sn+1
α
∆t
ρ

+ snα
sn+1
α
2

= T

!



V̄α∗

+

−

n+ 1
θα 2 Ū ∗

−T



n+ 1
sα 2


 
1
un+1 − un
= −∇p̄ + µ∇2 ū + Ln+ 1 F s̄α , θ n+ 2 ,
2
∆t
 

n+ 21
2
∇ p̄ = ∇ · Ln+ 1 F s̄α , θ
,



,

(178)
(179)
(180)

2

The leading order terms, T is given by

Tb(s̄α ) =

−





n+ 1
λ̄ min sα 2
α



  s
1 2 
n+ 2


Sb ∆t

n+ 1
λ̄ min sα 2
α

4ρ min sα
α

2

k2

+

2

k4

+ k2

−





b̄
|k|3 
 sα ,


(181)

(182)

2ρ
, and u∗,n+1 is the intermediate velocity at tn+1 which is obtained by solving
µ∆t
the unsteady Stokes equations implicitly but with the forcing evaluated explicitly:

where λ̄2 =

u∗,n+1 − un
ρ
∆t

2 ∗

∗

 

1
n+ 1
sα 2 , θ n+ 2

= −∇p̄ + µ∇ ū + Ln+ 1 F
 

2

n+ 21

∇2 p̄∗ = ∇ · Ln+ 1 F sα
2

1

, θ n+ 2

1

V̄ ∗ = L∗n+ 1 (ū∗ ) · τ n+ 2 ,





,

(183)

,

(184)
(185)

2

1

Ū ∗ = L∗n+ 1 (ū∗ ) · nn+ 2 .

(186)

2

and
ū∗ =

u∗,n+1 + un
2

(187)

Step 2: After the un+1 and sn+1
are calculated in step 1, we update θ at tn+1 using the
α
following semi-implicit scheme:
θ n+1 − θ n
∆t

= S

θ n+1 + θ n
2

!

+

1
n+ 21

sα

where
V̄
Ū

=

L∗n+ 1
2

un+1 + un
2

=

L∗n+ 1
2

un+1 + un
2
23





1

Ūα + θ̄α V̄ − S θ n+ 2

!
!

1

· τ n+ 2 ,
1

· nn+ 2 .



,

(188)

(189)
(190)

and the leading order term is

b θ̄)
S(

Sb ∆t max
= −

α






−1 
 3
|k| − s



1 3
n+ 2


n+ 1
sα 2

α

k4

n+ 1
λ̄ min sα 2
α

4ρ min sα



2

+ k2


 b̄
 θ,



(191)

This completely defines our second order semi-implicit scheme.

7

Numerical results

In this section, we will perform a number of numerical experiments to test the stability of our
semi-implicit schemes for both the steady and unsteady Stokes equations. We also compare
the performance of our semi-implicit schemes with the explicit scheme and the fully implicit
scheme. Our numerical results indicate convincingly that our semi-implicit schemes has a much
better stability property that that of the explicit scheme. Moreover, the computational cost of
our semi-implicit schemes is comparable to that of an explicit scheme. Thus our semi-implicit
schemes offer significant computational saving over the explicit scheme, especially when the
number of grid points is large.

7.1

Model problem

The test problem we use is one typically seen in the literature, in which the immersed boundary
is a closed loop initially in the shape of an ellipse. We choose an ellipse initially aligned in the
coordinate directions with horizontal semi-axis a = 0.32 and vertical semi-axis b = 0.24. The
boundary can be parameterized as follows:
(

x(α, 0) = 0.5 + 0.32 cos α,
y(α, 0) = 0.5 + 0.24 sin α.

(192)

The fluid is initially at rest in a periodic domain Ω = [0, 1] × [0, 1]. We use a periodic boundary
condition for the fluid flow. For the initial condition defined above, the rest state of the
boundary is a circle with radius r = 0.2. For the unsteady Stokes flow, the immersed boundary
with the above initial condition evolves as damped oscillations around a circular equilibrium
state. The area is conserved during the time evolution since the flow is incompressible. For
the steady Stokes flow, the boundary converges to the circular state without oscillations.
We use a uniform N × N grid to discretize the fluid domain, Ω. We choose Nb = 2N
number of grid points to discretize the immersed boundary so that there are approximately
2 immersed boundary points per mesh width. We use the spectral method to discretize the
spatial derivatives both in the fluid domain and along the immersed boundary. The leading
order singular integral is also discretized by the spectral method.

24

7.2

Steady Stokes flow

First, in order to reduce the number of parameters in our test problem, we write the equations
in terms of the following dimensionless variables to get the nondimensional model [33],
t′ =

t
,
t0

x′ =

x
,
L

ut0
,
L

u′ =

p′ =

pt0
,
µ

f′ =

f Lt0
,
µ

where L is the size of computational domain, t0 is characteristic time.
Using these new variables, we have
0 = −∇p′ + △u′ + f ′ (x′ , t′ ),

(193)

′

0 = ∇·u.

(194)

For the equations of the elastic boundary, the dimensionless variables are
X′ =

X
,
L

s′α =

sα
,
L

θ ′ = θ,

α′ =

α
,
L

T′ =

T
,
Sb

F′ =

FL
,
Sb

τ′ = τ,

n′ = n.

Then the equations describe the interaction of the boundary and the fluid become


U ′ = u′ X′ (α′ , t′ ), t′ · n′ ,

V

′

′

′

′

′

−

θα′ ′

′

= u X (α , t ), t

s′α,t′

=

θt′ ′

=

Vα′ ′

U,

′

· τ ′,


1
Uα′ ′ + V ′ θα′ ′ ,
′
sα

(195)
(196)
(197)
(198)

where
f ′ (x′ , t′ ) =
u′ (X′ (α′ , t′ ), t′ ) =

S b t0
µL

Z

Ω

Z

Lb /L

0

F′ (α′ , t′ )δ(x′ − X′ (α′ , t′ ))dα′ ,

u′ (x′ , t′ )δ(x′ − X′ (α′ , t′ ))dx′ .

(199)
(200)

There are two nondimensional parameters in this problem:
S b t0
,
µL

Lb
.
L

µL
, then the only parameter in this dimensionless model is Lb /L which is fixed
Sb
in our test problem. So we can always fix Sb = µ = 1 in our numerical study.
The stability analysis in the steady Stokes flow suggests us to use the total energy as a
criterion to test the stability of different numerical methods. For the steady Stokes equations,
the total energy is equal to the potential energy. In Fig 1, we show that the energy for four
different numerical methods: the explicit scheme, the semi-implicit scheme of first kind, the
4th order semi-implicit scheme using the integral factor method and the unconditionally stable
semi-implicit scheme. In this figure and the subsequent figures, we use the the legend “semiimplicit” to denote the semi-implicit scheme of first kind, and the legend “integral factor” to

If we let t0 =

25

0.5

0.12

stable semi−implicit
semi−implicit
explicit
integral factor

0.115

stable semi−implicit
semi−implicit
explicit
integral factor

0.45
0.4
0.35

0.11

Energy

Energy

0.3
0.105

0.25
0.2

0.1

0.15
0.1

0.095

0.05
0.09
0

2

4

6

8

0
0

10

2

4

6

8

10

Time

Time

Figure 1: Energy of the system for four different schemes. N = 128, Sb = 1, µ = 1. Left one:
∆t = 0.1; Right: ∆t = 1. Here the legend “stable semi-implicit” stands for the unconditionally
stable semi-implicit scheme, “semi-implicit” for the semi-implicit scheme of first kind, and
“integral factor” for the semi-implicit scheme based on the integral factor method.
denote the semi-implicit scheme based on the integral factor method. We use two different
time steps, 0.1 and 1, respectively. When ∆t = 0.1, all the four methods are stable. They
give almost identical results. When ∆t = 1, the explicit scheme becomes unstable , but all
the semi-implicit schemes are stable. In fact, all the semi-implicit schemes remain stable with
much larger time steps. In Fig 2, we plot the energy of the system for semi-implicit schemes
of first kind and the semi-implicit scheme based on the integral factor method with ∆t = 10.
Fig 3 shows the configuration obtained by the two semi-implicit schemes at the final time with
∆t = 10. They both remain as a circle, but lose some area compared with the original state.
Next, we compare the performance of our semi-implicit schemes with the explicit and fully
implicit schemes. The fully implicit scheme we use here was originally proposed by Tu and
Peskin in [33]. In order to make a fair comparison, we run the implicit schemes (semi-implicit
and fully implicit) with a time step small enough to make sure that the computational results
have a reasonable accuracy. We take ∆t = 4 for the fully implicit and the semi-implicit
schemes. With this time step, the area loss is less than 5%. For the explicit scheme, we
take ∆t = 1/4, 1/8, 1/16, 1/32 which corresponds to N = 64, 128, 256, 512 respectively. These
time steps are the largest possible to keep the stability of the explicit scheme. The time we
compute is T = 20. The result is shown in Table 1. From this comparison, we can see that the
performance of our semi-implicit schemes is much better than the explicit, the fully implicit
scheme, and the unconditionally stable semi-implicit scheme in all cases. As we can see, the
larger the number of the spatial grid points is, the more computational saving we would get
using our semi-implicit schemes. Even for a modest grid size, our semi-implicit schemes still
give a significant computational saving compared with the explicit or the fully implicit scheme.
It is interesting to note that although the computational cost of the unconditionally stable semiimplicit scheme (labeled as s,s,i) is faster than the fully implicit method, the computational
cost of the unconditionally stable semi-implicit method is still more expensive than the explicit
scheme. This makes the unconditionally stable semi-implicit scheme not very practical.

26

0.15

integral factor
semi−implicit

Energy

0.12

0.09

0.06

0.03
0

20

40

60

80

100
Time

120

140

160

180

200

Figure 2: Energy for the semi-implicit scheme of first kind (labeled as “semi-implicit”) and
the semi-implicit scheme based on the integral factor method (labeled as “integral factor”).
∆t = 10, N = 128, Sb = 1, µ = 1.

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0
0

0.2

0.4

0.6

0.8

0
0

1

0.2

0.4

0.6

0.8

1

Figure 3: Dashed line: the initial boundary configuration; Solid line: the boundary configuration after 20 time steps with ∆t = 10, N = 128, Sb = 1, µ = 1. Left: semi-implicit ; Right:
integral factor method.

27

N
64
128
256
512

exp
1
5
30
344

s,i 1
0.4
0.7
1.3
4.2

s,i 2
2
3
7
19

s,s,i
7
30
139
611

f,i
9
39
206
1200

Table 1: Execution time for each computation in seconds. The legends are defined as follows:
“exp” stands for the explicit scheme, “s,i1” the semi-implicit method of first kind, “s,i2”
the semi-implicit scheme of the second kind , “s,s,i” the unconditionally stable semi-implicit
method, and “f,i” the fully implicit scheme . N is the number of grid points along each
dimension.

7.3

Unsteady Stokes flow

We can also get the nondimensional model for unsteady stokes flow. Similar as the steady
stokes case, we define the following dimensionless variables:
t′ =

t
,
t0

x′ =

x
,
L

ut0
,
L

u′ =

p′ =

pt0
,
µ

f′ =

f Lt0
.
µ

where L is the size of computational domain, t0 is characteristic time. Using these new variables, we have

µt0
∂u′
=
−∇p′ + △u′ + f ′ (x′ , t′ ) ,
′
2
∂t
ρL
0 = ∇ · u′ .

(201)
(202)

For the equations of the elastic boundary, the dimensionless variables are
X′ =

X
,
L

s′α =

sα
,
L

θ ′ = θ,

α′ =

α
,
L

T′ =

T
,
Sb

F′ =

FL
,
Sb

τ′ = τ,

n′ = n.

Then the equations describe the interaction of the boundary and the fluid become
U ′ = u′ (X′ (α′ , t′ ), t′ ) · n′ ,

V

′

s′α,t′

′

′

′

′

′

= u (X (α , t ), t ) ·
Vα′ ′

=

−

τ ′,

θα′ ′ U ′ ,


1
Uα′ ′ + V ′ θα′ ′ ,
′
sα

θt′ ′ =

(203)
(204)
(205)
(206)

where
′

′

′

f (x , t ) =
u′ (X′ (α′ , t′ ), t′ ) =

S b t0
µL

Z

Ω

Z

Lb /L

0

F′ (α′ , t′ )δ(x′ − X′ (α′ , t′ ))dα′ ,

u′ (x′ , t′ )δ(x′ − X′ (α′ , t′ ))dx′ .

28

(207)
(208)

From the nondimensional analysis, we can see that there are three nondimensional parameters
in this problem:
S b t0
,
µL

µt0
,
ρL2

Lb
.
L

µt0
µ2
µL
, then the parameters left in this dimensionless model is
and
=
Sb
ρL2
ρLSb
Lb
µ2
. Lb and L are fixed in our test problem only depends on the initial condition. So
is
L
ρLSb
the only parameter in our test model. For this reason, we always fix the elastic coefficient Sb
to 1, but vary µ.
In our computations, we use the following parameter values:
If we let t0 =

ρ = 1, Sb = 1, µ = 0.1, 0.01, 0.005.
We vary the number of the spatial grid points along each dimension in the following fashion:
N = 64, 128, 256, 512.
The criterion that we use to check whether one scheme is stable or not is that the total
energy of the system is non-increasing and the boundary configuration lies within the computational domain.
Next, we perform some numerical experiments to test the stability of our semi-implicit
schemes for the unsteady Stokes flow. Fig. 4 shows that the energy obtained by the explicit
scheme and the semi-implicit scheme of the first kind. We take two different timesteps, 0.005
and 0.05. With ∆t = 0.005, the explicit and semi-implicit schemes are all stable, and they give
nearly identical results. With ∆t = 0.05, the explicit scheme becomes unstable, but the semiimplicit schemes remain stable. Even if we increase the timestep to ∆t = 1, the semi-implicit
methods are still stable, as we can see from Fig. 5. For the semi-implicit scheme of the first
kind, we have used the Small Scale Decomposition and further simplification of the singular
integral kernel. Therefore, the total energy in our semi-implicit scheme is not guaranteed to
decrease monotonically in time. Nonetheless, we observe that the total energy still decreases
in time as is the case for the unconditionally stable semi-implicit scheme. In Fig. 5, we also
plot the boundary configuration at the final time step, which is an approximate circle.
We remark that the semi-implicit scheme is not unconditionally stable, although its stability
is much better than the explicit scheme. This is due to the fact that we have used the Small
Scale Decomposition and further approximation of the leading order singular integral operator
to simplify the computation of the implicit solution. As we mentioned before, the Small Scale
Decomposition captures only the high frequency contribution to the stiffness, but it does not
remove the stiffness of the system induced by the low frequency components of the solution.
Thus there is still some mild timestep stability constraint for the semi-implicit scheme. The
time step has a mild dependence on the elastic coefficient Sb and the viscous coefficient µ. On
the other hand, our numerical study shows that the time step for the semi-implicit scheme is
independent on the meshsize.
We also compare the performance of our semi-implicit schemes with the explicit scheme. We
do not compare the performance of our semi-implicit schemes with the fully implicit scheme
29

0.5

0.12

stable semi−implicit
semi−implicit
explicit

0.115

stable semi−implicit
semi−implicit
explicit

0.45
0.4
0.35
0.3

Energy

Energy

0.11

0.25
0.2

0.105

0.15
0.1

0.1

0.05
0.095
0

0.1

0.2

0.3

0.4

0
0

0.5

0.1

0.2

0.3

0.4

0.5

Time

Time

Figure 4: Total energy of the unsteady Stokes system for different schemes with two different
timesteps. N = 128, Sb = 1, µ = 0.01. Left: ∆t = 0.005; Right: ∆t = 0.05. The legend “semiimplicit” stands for the solution obtained by the semi-implicit scheme of first kind, “stable
semi-implicit” the unconditionally stable semi-implicit method.

0.16

stable semi−implicit
semi−implicit

0.14
0.12

Energy

0.1
0.08
0.06
0.04
0.02
0
0

5

10
Time

15

20

Figure 5: Total energy of the system for two semi-implicit schemes. Here “stable semiimplicit” stands for the unconditionally stable semi-implicit scheme, and “semi-implicit” the
semi-implicit scheme of first kind, ∆t = 1, N = 128, Sb = 1, µ = 0.01.

30

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0
0

0.2

0.4

0.6

0.8

0
0

1

0.2

0.4

0.6

0.8

1

Figure 6: Dashed line: the initial boundary configuration; Solid line: the boundary configuration after 20 time steps with ∆t = 1, N = 128, Sb = 1, µ = 0.01. Left: the unconditionally
stable semi-implicit scheme; Right: the semi-implicit scheme of the first kind.
here because the fully implicit scheme is quite expensive and is not competitive with the
explicit scheme. In order to keep the area loss is no more than 5%, we take ∆t = 14 for all of
the semi-implicit schemes . For the explicit scheme, we take ∆t = 1/64, 1/128, 1/256, 1/512
which correspond to the spatial mesh sizes N = 64, 128, 256, 512 respectively, when µ = 0.05.
When µ = 0.01 and 0.005, the time step is set to be ∆t = 1/128, 1/256, 1/512, 1/1024 and
t = 1/256, 1/512, 1/1024, 1/2048. These time steps are the largest ones we can take to keep the
stability. We compute the solution up to T = 2. The results are documented in Table 2. We
can clearly see that the semi-implicit scheme of the first kind gives a significant improvement
over the explicit scheme. The cost for the semi-implicit scheme of the second kind is higher
than that for the semi-implicit scheme of the first kind. This is because we need to solve
for a linear system at each time step for the semi-implicit scheme of the second kind. The
cost increases as the number of the spatial grid points increases. The semi-implicit scheme
of the second kind and the unconditionally stable semi-implicit scheme both need to solve a
linear system at each time step. Their complexity are same, both are O(Nb2 ). But for the
unconditionally stable semi-implicit scheme, the scaling constant in front of Nb2 is much larger
than the semi-implicit scheme of the second kind. The reason is that the cost of computing
the coefficient matrix of the linear system for the unconditionally stable semi-implicit scheme
is much higher. As we can see from Table 2, the unconditionally stable semi-implicit scheme
(labeled as s,s,i) is still quite expensive compared with our semi-implicit schemes that use the
Small Scale Decomposition. For µ = 0.05, the unconditionally stable semi-implicit scheme is
even more expensive compared with the explicit scheme. Although the unconditionally stable
semi-implicit scheme is slightly faster than the explicit scheme for smaller µ, the semi-implicit
scheme (labeled as s,i,1) which uses SSD to further simply the singular integral kernel gives a
much more efficient algorithm. It gives a factor of 242 times speed-up over the explicit scheme
in the case of µ = 0.005 and N = 512.

31

N
64
128
256
512

exp
1.8
9
58
738

µ = 0.05
s,i 1 s,i 2
0.5
4
1
10
2.4
25
12
99

s,s,i
11
48
229
980

exp
3.3
18
116
1461

µ = 0.01
s,i 1 s,i 2
0.5
4
0.9
10
2.4
25
12
98

s,s,i
12
47
228
982

exp
6.6
35
236
2910

µ = 0.005
s,i 1 s,i 2
0.5
4
0.9
10
2.4
25
12
98

s,s,i
12
48
226
977

Table 2: Execution times for each computation in seconds. The legends are defined as follows:
“exp” stands for the explicit scheme, “s,i1” the semi-implicit scheme of the first kind, “s,i2”
the semi-implicit scheme of the second kind, “s,s,i” the unconditionally stable semi-implicit
scheme.

7.4

Second order semi-implicit scheme for the unsteady Stokes flow

In this subsection, we perform numerical experiments to test the convergence rate and the
stability property of our second order semi-implicit scheme. To check the convergence rate in
1
1
to 128
. When µ = 0.005,
time, we set N = 256 and vary the time step in powers of 2 from 16
the solution becomes more singular. In order to fully resolve the spatial solution, we increase
the spatial resolution to N = 512. Following [23], we compute the time discretization error at
time t as follows:
eT (u; ∆t) = ku(T ; ∆t) − u(T ; ∆t/2)kl2 .

(209)

For a vector field u(x) = (u1 (x), u2 (x)) defined on the Cartesian grid with xi = ih, yj = jh,
the discrete l2 norm is defined as follows
kukl2



=

X

u21 (xi , yj )

+



u22 (xi , yj )

i,j

1

2

2

h

.

(210)

Similarly, the discrete l2 norm for a vector field w(α) = (w1 (α), w2 (α)) defined on the interface
Γ is defined below:
kwkl2 =

X



!1

w12 (αi ) + w22 (αi ) ∆α

i

2

.

(211)

We compute the solution up to T = 1 and evaluate the convergence rate based on the numerical
solution at T = 1 with different temporal resolutions. The results are shown in Fig. 7 and
Table 3. As we can see, the convergence rate is approximately second order.
Next we check the stability of the second order semi-implicit scheme. Fig 8 shows that
the total energy for the second order explicit and second order semi-implicit schemes with
different timesteps 0.002 and 0.02. We choose the same second order explicit scheme that was
used in [13]. With ∆t = 0.002, both the explicit and semi-implicit schemes are stable and
they give nearly identical results. With ∆t = 0.02, the explicit scheme becomes unstable, but
the semi-implicit scheme is still stable. The stability restriction of the semi-implicit scheme is
far less severe than the corresponding explicit scheme. Our numerical study also shows that
32

error of X

−2

error of velocity field u

−1

10

10

µ=0.01
µ=0.05
µ=0.005

−3

10

µ=0.01
µ=0.05
µ=0.005

−2

10

−3

10
−4

10

−4

10
−5

10

−5

10

−6

−6

10

1

10

2

10

1

10

2

10

10

1/∆ t

1/∆ t

Figure 7: Plot of l2 errors in time at time T = 1 for the second order semi-implicit scheme
of the first kind. We choose Sb = 1 and N = 256 in all computations except in the case of
µ = 0.005 where N is increased to 512. The line at the bottom of each graph is a reference
line which corresponds to the second order convergence rate.

µ
0.05
0.01
0.005∗

convergence rate of X
2.11
2.13
2.17

convergence rate of u
2.70
2.10
1.96

Table 3: Convergence rates for X and u fitted from the data shown in Fig 7. the case * is
computed using a refined mesh 512 × 512.

0.5

0.12

semi−implicit
explicit

semi−implicit
explicit

0.45
0.4

0.115

0.35
0.3

Energy

Energy

0.11

0.25
0.2

0.105

0.15
0.1

0.1

0.05
0.095
0

0.1

0.2

0.3

0.4

0
0

0.5

0.1

0.2

0.3

0.4

0.5

Time

Time

Figure 8: Total energy of the unsteady Stokes system for the second order semi-implicit scheme
and the second order explicit scheme. N = 128, Sb = 1, µ = 0.01. Left: ∆t = 0.002; Right:
∆t = 0.02.

33

N
64
128
256
512

µ = 0.05
explicit semi-implicit
7
0.8
37
1.6
249
4.4
3088
24

µ = 0.01
explicit semi-implicit
7
0.8
38
1.6
504
4.6
6182
25

µ = 0.005
explicit semi-implicit
7
0.8
38
1.6
506
4.5
6200
25

Table 4: Execution time for each computation in seconds. Here “explicit” stands for the second
order explicit scheme and “semi-implicit” the second order semi-implicit scheme.
the time step for the semi-implicit scheme is independent on the meshsize while the explicit
scheme requires finer time step as the spatial mesh is refined.
Finally, we compare the performance of our second order semi-implicit scheme with that
of the second order explicit scheme. As before, in order to keep the accuracy with 5%,
we take ∆t = 41 for our semi-implicit schemes. For the explicit scheme, we take ∆t =
1/128, 1/256, 1/512, 1/1024 which correspond to the spatial mesh sizes N = 64, 128, 256, 512
respectively, when µ = 0.05. When µ = 0.01 and 0.005, the time step is set to be ∆t =
1/128, 1/256, 1/1024, 1/2048. These time steps are the largest ones that we can take to keep
the stability of the explicit scheme. We compute the solution up to T = 2. The result is shown
in Table 4. Again, we observe the same qualitative behavior as the first order schemes we
reported earlier.

8

Concluding Remarks

In this paper, we have developed several efficient semi-implicit immersed boundary methods
for solving the immersed boundary problem for the steady and unsteady Stokes equations.
The immersed boundary method has emerged as one of the most useful numerical methods
in computing fluid structure interaction, and has found numerous applications. But it also
suffers from the severe time step stability limitation due to the stiffness of the elastic force.
Guided by our stability analysis, we have developed several efficient semi-implicit schemes
which remove the stiffness of the immersed boundary method. We have demonstrated both
analytically and computationally that our semi-implicit schemes have much better stability
property than the explicit scheme. More importantly, unlike most existing implicit or semiimplicit schemes, our semi-implicit schemes can be implemented very efficiently. In fact, our
semi-implicit scheme of the first kind has a computational cost that is essentially the same as
that of an explicit scheme in each time step, but with a much better stability property. The
saving in the computational cost is quite substantial. We have demonstrated this improved
stability for a range of parameters and numerical resolutions. Our computational results show
that the larger the spatial resolution is, the bigger the computational saving our semi-implicit
schemes can offer. Thus the semi-implicit schemes we develop in this paper provide an effective
alternative discretization to the explicit method.
One of the essential steps in developing our semi-implicit schemes is to obtain an unconditionally stable semi-implicit discretization of the immersed boundary problem. This provides
us with a building block to construct our efficient semi-implicit schemes. There are two im34

portant observations in constructing the unconditionally stable semi-implicit discretization.
The first one is that we need to preserve certain important solution structures at the discrete
level. Specifically, we need to ensure that the discrete velocity field is divergence free, and
the discrete spreading and interpolation operators are adjoint. Another essential step is to
decouple the stiffness of the elastic force from the fluid flow in some appropriate way. This
is difficult to achieve if we use the Cartesian coordinate. But it becomes easier if we use the
arclength and tangent angle formulation to describe the dynamics of the immersed interface
as was done in [9]. On the other hand, as we demonstrated in this paper, the unconditionally
stable semi-implicit scheme is still very expensive to implement, and the saving over the expicit
scheme is rather limited.
Based on this unconditionally stable semi-implicit discretization, we have developed several efficient schemes for both the steady and the unsteady Stokes flows. By applying the
Small Scale Decomposition to the unconditionally stable semi-implicit time discretization and
further simplifying the leading order singular kernel, we obtain our semi-implicit scheme. The
advantage of this semi-implicit scheme is that the leading order term can be expressed as
a convolution operator, which can be evaluated explicitly using the Fourier transformation.
This allows us to solve for the implicit solution explicitly in the spectral space, which offers
substantial computational saving over the explicit scheme.
It is a natural step to extend the the semi-implicit schemes developed for the unsteady
Stokes equations to the Navier-Stokes equations. The discretization of the Navier-Stokes equations shares many similar properties as the unsteady Stokes equations if we treat the convection
term explicitly. We have performed a number of numerical experiments to test the stability
and the robustness of our semi-implicit immersed boundary methods for the Navier-Stokes
equations. The results are qualitatively similar to those for the unsteady Stokes equations
which we have presented in this paper. These results will be reported in a subsequent paper.
Acknowledgments. We would like to thank Profs. Charles Peskin and Hector Ceniceros for
a number of stimulating discussions on the Immersed Boundary method. The research was in
part supported by DOE under the DOE grant DE-FG02-06ER25727 and by NSF under the
NSF FRG grant DMS-0353838, ITR grant ACI-0204932, and DMS-0713670.

A

The semi-implicit scheme of the second kind

In this appendix, we will derive the semi-implicit method of the second kind in more detail.
As we mentioned before, the small scale decomposition only captures the leading order contribution from the high frequency components, which can not remove the stiffness induced by Sb
and µ completely. The coefficients Sb and µ can still affect the time stability through the low
frequency components of the solution, which comes from the lower order term of the right hand
side. In order to obtain a semi-implicit discretization with better stability property, we can
incorporate the low frequency contribution from the second term in our implicit discretization.
For the steady Stokes flow, this gives rise to the following decomposition:
ŝn+1
− ŝnα
α
∆t



Sb
= − |k|ŝn+1
+ F θαn
α
4µ

Z

′

ln |α − α

35

|(sn+1
α

−

1)θαn dα′



,





+ F Vαn − θαn U n − θαn
φ̂n+1 − φ̂n
∆t

= −

Sb
γ|k|φ̂n+1 + F
4µ







Z





1



n n+1
+ F
n+1 V θα

sα





Uαn
sn+1
α

+



Sb
|k|ŝnα , (212)
4µ

ln |α − α′ |(snα − 1)θαn dα′ +



Sb
γ|k|φ̂n ,
4µ

(213)

1
. By replacing the continuous derivative by the discrete derivative,
α
sα
and discretizing the continuous integral by the trapezoidal rule, we obtain our second semiimplicit scheme. We call this semi-implicit scheme the semi-implicit scheme of the second
kind. Near equilibrium, we can show that the semi-implicit scheme of the second kind is
unconditionally stable. Moreover, the stability property is independent of the meshsize, elastic
coefficient Sb and viscosity coefficient µ. Our numerical study also confirms this. The trade-off
is that we need to solve a linear system at each time step to obtain the implicit solution at
tn+1 .
Similarly, in the case of unsteady Stokes flow, we can also include the second term of the
right hand side in the leading order term to derive a scheme with better stability property. In
this case, the leading order term becomes:
where γ = max 1 −



λ min snα
α

2

k2



k4

+
 n+1

Sb ∆t
r
− |k|3 
Tb (sn+1
2 
α )=− 

 ŝα

2
n
n
2
2 min sα
λ min sα + k
α
α

Sb ∆t
F
−
2π

θαn
(snα )2

b n+1 ) = −
S(θ

+F



1

sn+1
α

V

Z

K0 (λsnα |α

Sb ∆t max (snα
α


2 min snα
α

n+1 n+1
θα



′

′

− α |) − ln(α − α )


− 1) 

2



|k|3 − r



(sn+1
α

−



k4
λ min snα
α

2

1)θαn′

+ k2



α′ α′

′

dα

!

, (214)

 n+1
 θ̂


,

(215)

where the derivative will be discretized by the spectral method and the integration will be
discretized by the trapezoidal rule. We call the above scheme the semi-implicit scheme of the
second kind for the unsteady Stokes flow. Near equilibrium stability analysis suggests that the
semi-implicit scheme of the second kind is unconditionally stable.

B

Derivation of the semi-implicit scheme (140)-(143)

In this appendix, we will derive the semi-implicit scheme (140)-(143). We define the operator
G(sα ; un , θ n , Xn ) : sα → u by the following equations:
u − un
= −∇h p + µ∇2h u + Lh,n (F(sα , θ n )),
∆t
∇2h p = ∇h · Lh,n (F(sα , θ n )).

36

(216)
(217)
(218)

Given sα , we obtain u by solving above equations. From the definition of operator G, we have
n n
n
un+1 = G(sn+1
α ; u , θ , X ),

u∗,n+1 = G(snα ; un , θ n , Xn ),

where u∗,n+1 is calculated from equations (140)-(141).
Then the equation of sα can be rewritten as
sn+1
− snα
α
∆t

= D∆α V n+1 − D∆α θ n U n+1








= D∆α L∗h,n (un+1 ) · τ n − D∆α θ n L∗h,n (un+1 · nn )






n n
n
n
= D∆α L∗h,n G(sn+1
α ;u ,θ ,X ) · τ









n n
n
n
−D∆α θ n L∗h,n G(sn+1
α ;u ,θ ,X ) · n

n n
n
≡ R(sn+1
α ; u , θ , X ).



(219)

Using the small scale decomposition, we get
n n
n
n+1
R(sn+1
α ; u , θ , X ) ∼ T (sα ).

(220)

By treating the leading order term implicitly, we obtain our semi-implicit method as follows:
sn+1
− snα
α
∆t

n
n n
n
n
= T (sn+1
α ) + (R(sα ; u , θ , X ) − T (sα ))



∗
n
n n
n
n
= T (sn+1
α ) + D∆α Lh,n (G(sα ; u , θ , X )) · τ







−D∆α θ n L∗h,n (G(snα ; un , θ n , Xn )) · nn − T (snα )






∗
∗,n+1
) · τ n − D∆α θ n L∗h,n (u∗,n+1 ) · nn
= T (sn+1
α ) + D∆α Lh,n (u

−T (snα )





∗,n+1
= T (sn+1
− D∆α θ n U ∗,n+1 − T (snα ) .
α ) + D∆α V



(221)

This is exactly our semi-implicit scheme (134).
In the steady Stokes case, we can define the operator G(sα ; un , θ n , Xn ) : sα → u similarly:
0 = −∇h p + µ∇2h u + Lh,n (F(sα , θ n )),

∇2h p

n

= ∇h · Lh,n (F(sα , θ )).

(222)
(223)

In this case, we have
G(snα ; un , θ n , Xn ) = un .
Thus, the semi-implicit scheme becomes
− snα
sn+1
α
n
n n
n
= T (sn+1
α ) + (D∆α V − D∆α θ U − T (sα )) .
∆t
37

(224)

C

The derivation of the Fourier transform of K0

In this appendix, we derive the the Fourier transform of K0 , which is given in (144). By the
definition of the Fourier transform, we have
F

Z

1
π

Z

+∞

−∞

Z

K0 (β|α − α′ |)f (α′ )dα′





+∞
1 +∞
K0 (β|α − α′ |)f (α′ )dα′ e−ikα dα
π −∞
−∞
Z
Z
1 +∞ +∞
K0 (β|α − α′ |)f (α′ )e−ikα dα′ dα
π −∞ −∞
Z
Z
′
′
1 +∞ +∞
K0 (β|α − α′ |)f (α′ )e−ik(α−α ) e−ikα dα′ d(α − α′ )
π −∞ −∞
Z +∞

Z
1 +∞
′
′
K0 (β|α − α′ |)e−ik(α−α ) d(α − α′ ) f (α′ )e−ikα dα′
π −∞
−∞

=
=
=
=
=



Z

+∞

−∞

where F(f (α))(k) =

Z

fb(k)
f (α′ )
′
−ikα′
p
dα
=
e
,
β 2 + k2
β 2 + k2

p

+∞

(225)

f (α)eikα dα is the Fourier transform. In the calculation above, we

−∞

have used the expression of the Bessel function (p 376, [1])
K0 (x) =

Z

+∞
0

cos(tx)
√
dt,
1 + t2

(226)

and the identity below:
Z

+∞
−∞

K0 (β|x|)e−ikx dx =
=
=
=
=

This proves (144).

Z

+∞ Z +∞

cos(βtx) −ikx
√
e
dtdx
1 + t2
−∞
0
Z
Z
1 +∞ +∞ cos(βtx) −ikx
√
e
dtdx
2 −∞ −∞
1 + t2
Z
Z
1 +∞ +∞ eiβtx −ikx
√
e
dtdx
2 −∞ −∞
1 + t2
Z
Z
1 +∞ +∞ ei(βt−k)x
√
dtdx
2 −∞ −∞
1 + t2
Z +∞
δ(βt − k)
π
√
π
.
dt = p 2
β + k2
1 + t2
−∞

(227)

References
[1] M. Abramowitz and I. A. Stegun, Modified Bessel Functions I and K, New York:
Dover, 9th printing (1972).

38

[2] R. P. Beyer, A computational model of the cochlea using the immersed boundary method,
J. Comput. Phys., 98(1992) pp. 145-162
[3] R. Cortez, N. Cowen, R. Dillon, L. J. Fauci, Simulation of swimming organisms: Coupling
internal mechanics with external fluid dynamics, Comput. Sci. Eng. 6 (2004), pp. 38-45
[4] H. D. Ceniceros, A. M. Roma, Study of the long-time dynamics of viscous vortex sheet
with a fully adaptive nonstiff method, Physics of Fluid 16 (2004), pp 4285-4318
[5] L. Fauci and C. S. Peskin, A computational model of aquatic animal locomotion, J. Comput. Phys., 77 (1988), pp. 85-108
[6] L. J. Fauci, Interaction of oscillating filaments - A computational study, J. Comput. Phys,
86(1990) pp. 294-313
[7] A. L. Fogelson, A mathematical model and numerical methods for studying platelet adhesion and aggregation during blood clotting, J. Comput. Phys., 56 (1984), pp. 111-134
[8] E. Givelberg, Modeling elastic shells immersed in fluid, PhD thesis (1997), Courant Institute of Mathematical Sciences, New York University
[9] T. Y. Hou, J. Lowengrub and M. Shelley, Removing the stiffness from interfacial flows
with surface tension, J. Comput Phys., 114 (1994), pp. 312-338.
[10] T. Y. Hou, J. Lowengrub and M. Shelley, The long-time motion of vortex sheets with
surface tension, Phys. Fluid A, 9 (1997), pp. 1933-1954.
[11] M. Hopkins and L. J. Fauci, A computational model of the collective fluid dynamics of
motile microorganisms, J. Fluid Mech., 455(2002), pp. 149-174
[12] E. N. Jung and C. S. Peskin, Two-dimensional simulations of valveless pumping using the
immersed boundary method, SIAM J. Sci. Comput. 23 (2001), pp. 19-45.
[13] 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
[14] R. J. LeVeque and Z. Li, Immersed Interface methods for Stokes flow with elastic boundaries or surface tension, SIAM J. Sci. Comput. 18 (1997), pp. 709-735.
[15] A. A. Mayo and C. S. Peskin, An implicit numerical method for fluid dynamics problems
with immersed elastic boundaries, Contemporary Mathematics, 141 (1993), pp. 261-277.
[16] D. M. McQueen, C. S. Peskin, and E. L. Yellin, Fluid dynamics of the mitral valve:
Physiological aspects of a mathematical model, Amer. J. Physiol. 242 (1982), pp. H1095H1110.
[17] D. M. McQueen and C. S. Peskin, Computer assisted design of pivoting-disc prosthetic
mitral valves, J. Thorac. Cardiovasc. Surg., 86 (1983), pp. 126-135.
[18] D. M. McQueen and C. S. Peskin, Computer assisted design of butterfly bileaflet valves
for mitral position, Scand. J. Thorac. Cardiovasc. Surg., 19 (1985), pp. 139-148.
39

[19] D. M. McQueen and C. S. Peskin, A three-dimensional computational method for blood
flow in the heart: (II) Contractile fibers, J. Comput. Phys., 82 (1989), pp. 289-297.
[20] D. M. McQueen and C. S. Peskin, Heart simulation by an immersed boundary method
with formal second order accuracy and reduced viscosity. In Mechanics for a New Millennium, Proceedings of the International Conference on Theoretical and Applied mechanics
(ICTAM) 2000, Kluwer Academic Publishers.
[21] L. A. Miller and C. S. Peskin, When vortices stick: and aerodynamic transition in tiny
insects, J. Exp. Biol., 207(2004), pp. 3073-3088
[22] L. A. Miller and C. S. Peskin, A computational fluid dynamics of ’clap and fling’ in small
insects, J. Exp. Biol., 208(2005), pp. 195-212
[23] Y. Mori and C. S. Peskin, A second order implicit Immersed Boundary method with boundary mass. Preprint, 2007, submitted to Computational Mechanics.
[24] E. P. Newren, A. L. Fogelson, R. D. Guy and R. M. Kirby Unconditionally stable discretizations of the immersed boundary equations, J. Comput. Phys., 222 (2007), pp. 702-719.
[25] C. S. Peskin, Numerical study of blood flow in the heart, J. Comput. Phys., 25 (1977), pp.
220-252.
[26] C. S. Peskin and D. M. McQueen , A three-dimensional computational method for blood
flow in the heart: (I) Immersed elastic fibers in a viscous incompressible fluid, J. Comput.
Phys., 81 (1989), pp. 372-405.
[27] C. S. Peskin, The immersed boundary method, Acta Numerica, 11 (2002), pp. 479-517
[28] C. Pozrikidis, Boundary Integral and Singularity Methods for Linearized Viscous Flow, Cambridge University Press, 1992.
[29] M. E. Rosar and C. S. Peskin, Fluid flow in collapsible elastic tubes: a three-dimensional
numerical model, New York J. Math. 7 (2001), pp. 281-302.
[30] J. M. Stockie and B. R. Wetton, Stability analysis for the immersed fiber problem, SIAM
J. Appl. Math, 55 (1995), pp. 1577-1591.
[31] , J. M. Stockie and S. I. Green, Simulating the motion of flexible pulp fibers using the
immersed boundary method, J. Comput. Phys., 147(1998), pp. 147-165
[32] J. M. Stockie and B. R. Wetton, Analysis of stiffness in the Immersed Boundary method
and implications for time-stepping schemes, J. Comput Phys, 154 (1999), pp. 41-64.
[33] C. Tu and C. S Peskin, Stability and instability in the computation of flows with moving
immersed boundaries: A comparison of three methods, SIAM J. Sci. Stat. Comput, 13
(1992), pp. 1361-1376.
[34] N. T. Wang and A. L. Fogelson, Computational methods or continuum models of platelet
aggregation, J. Comput. Phys., 151 (1999), pp. 649-675.
40

[35] L. Zhu and C. S. Peskin, Simulation of a flexible flapping filament in a flowing soap film
by the immersed boundary method, J. Comput. Phys., 179(2002), pp. 452-468
[36] J. Zhu, X. Chen, T. Y. Hou, An efficient boundary integral method for the Mullins-Sekerka
problem, J Comput Phys., 127 (1996), pp. 246-267.

41

