A note on pressure accuracy in immersed boundary method for Stokes flow

Short note

A note on pressure accuracy in immersed boundary method for Stokes flow

Kuan-Yu Chen

a

, Ko-An Feng

b

, Yongsam Kim

c

, Ming-Chih Lai

d,⇑

a

Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 30010, Taiwan

b_{National Center for Theoretical Sciences, National Tsing Hua University, No. 101, Sec. 2, Kuang Fu Road, Hsinchu 30043, Taiwan} c_{Department of Mathematics, Chung-Ang University, Dongjak-Gu, Heukseok-Dong 221, Seoul 156756, South Korea}

d

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

a r t i c l e

i n f o

Article history:

Received 16 November 2010

Received in revised form 15 February 2011 Accepted 9 March 2011

Available online 16 March 2011 Keywords:

Immersed boundary method Pressure accuracy Discrete delta function Stokes flow

a b s t r a c t

In this short note, we provide a simplified one-dimensional analysis and two-dimensional numerical experiments to predict that the overall accuracy for the pressure or indicator function in immersed boundary calculations is first-order accurate in L1norm, half-order

accurate in L2norm, but has O(1) error in L1norm. Despite the pressure has O(1) error near

the interface, the velocity field still has the first-order convergence in immersed boundary calculations.

Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction

The immersed boundary (IB) method was first proposed by Peskin in[8]as a computational tool to study the blood flow in the heart. Over the past few decades, the IB method has become a useful computational framework for solving fluid–struc-ture interaction problems, see a recent review[9]. Whereas the IB method has an impressive ease of the implementations, it is well-known that the method is only first-order accurate since it is a smoothing method rather than the sharp capturing method. Recently, a rigorous convergence proof of the velocity for Stokes flow in the immersed boundary formulation has been provided by Mori[7]. The author has proved that the velocity is roughly first-order accurate in the L1norm; however, it has no conclusion about the accuracy of the pressure.

In[2], Beyer and LeVeque have analyzed the accuracy of one-dimensional model for the immersed boundary method. They gave some insight about the accuracy of solving one-dimensional heat equation with a delta source term by choosing appropriate discrete delta functions. Tornberg and Engquist[11]used the regularization technique to analyze the numerical accuracy of some elliptic PDEs. They have verified numerically that, for two-dimensional Poisson equation with singular del-ta source term, the sdel-tandard centered difference approximation with smoothing discrete deldel-ta function is first-order accu-rate in L1norm and second-order accurate in L1norm. They also showed that, away from the interface, the scheme has a better accuracy which is expected in the case of smooth problems without singular source term. Another more accurate ap-proach for solving PDEs with singular sources is to incorporate the solution or its derivative jumps across the interface into the finite difference scheme such as the immersed interface method[6]. Readers who are interested in IIM can refer to a recent survey book by Li and Ito[5].

0021-9991/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2011.03.019

⇑Corresponding author. Tel.: +886 3 5731361; fax: +886 3 5724679. E-mail address:mclai@math.nctu.edu.tw(M.-C. Lai).

Contents lists available atScienceDirect

Journal of Computational Physics

In this paper, we shall provide a simplified one-dimensional analysis and two-dimensional numerical experiments to pre-dict that the pressure in immersed boundary calculations is first-order accurate in L1norm, half-order accurate in L2norm, but has O(1) error in L1norm. Notice that, it has no surprise that the pressure has O(1) error near the interface, since it is discontinuous across the interface, and a smoothing method such as the immersed boundary method will not be able to cap-ture the right discontinuity. A rigorous convergence proof for the two-dimensional problems will need a further investiga-tion in the future.

Since the pressure appears in its gradient form in Stokes equations, one might wonder how this O(1) point-wise error near the interface affects the overall accuracy of the velocity field. The pressure gradient and the immersed boundary force terms have the same singular behavior (a delta function singularity), finding the velocity field involves solving Poisson equation with a singular delta function source. It is known that the numerical approximation of Poisson equation with a singular delta function source only has first-order accuracy in L1norm[11]. Even though the calculated pressure has O(1) error near the interface, its gradient has the same discretization error as the singular delta function force near the interface. Thus, the over-all accuracy of velocity field in IB method is first-order accurate in L1norm.

2. Model problems

The problem which we are interested in arises from solving the stationary Stokes flow in a two-dimensional domainX

with one-dimensional boundary (or interface)Cimmersed inXas



r

p þ

l

D

v

þ Z C FðsÞ d2 ðx  XðsÞÞ ds ¼ 0; ð1Þ

r



v

¼ 0: ð2Þ

The immersed boundary force F(s) = (Fx_{(s), F}y_{(s)) is only exerted along the interfaceCso that the integral is along the} one-dimensional interfaceCwhile the Dirac delta function d2_{(x) = d(x)d(y) is two-dimensional. Thus, the above immersed} bound-ary formulation is a typical singular problem with a delta function source. By taking the divergence operator on Eq.(1)and using the incompressibility constraint(2), we obtain the pressure equation

D

pðxÞ ¼

r

 Z

C

FðsÞ d2

ðx  XðsÞÞ ds: ð3Þ

Notice that, this equation is a Poisson equation with a source term which involves derivatives of the Dirac delta function. Once pressure is found, we can find the velocity by solving Eq.(1)which results two Poisson equations for the velocity com-ponents. In the immersed boundary computations, the pressure equation often uses periodic or Neumann boundary condi-tion; however, throughout this paper, we simply use the Dirichlet boundary condition since we are more concerned with the accuracy caused by the derivatives of Dirac delta function near the interface.

Another way of solving Stokes equations(1) and (2)is to discretize them directly and solve the resultant linear system for the velocity and pressure simultaneously. This procedure is mostly adopted by the problem with periodic boundary condi-tions in which FFT can be conveniently applied or fast Stokes solver is available. It seems that the above two solution pro-cedures are different; however, they share the same accuracy behavior for the pressure. The reason is that the pressure in IB formulation is discontinuous across the interface and both solution procedures involve the derivatives of the pressure and singular terms. Therefore, the pressure accuracy predicted in the paper does not change despite different solution procedures.

Another example leading to the same type of equation as Eq.(3)appears when we use the idea of Unverdi and Tryggvason

[10]to find the indicator function which is needed to track the regions of two-phase flow. If the viscosity is discontinuous across the immersed boundary, it can be represented by the following:

l

ðxÞ ¼

l

outþ ð

l

in

l

outÞIðxÞ;

where

l

inand

l

outare the viscosity inside and outside the interface, respectively. Since the indicator function has the value one (I = 1) inside the immersed boundaryCand the value zero (I = 0) outside, it can be calculated as the following procedure

[10]. LetXinrepresent the interior region and n be the unit outward normal vector to the interface, then the indicator func-tion can be represented by

IðxÞ ¼ Z

Xin

d2ðx exÞdex:

By taking the gradient and then divergence operators on both sides, we have

r

IðxÞ ¼  Z C n d2ðx  XðsÞÞ ds;

D

IðxÞ ¼ 

r

 Z C n d2 ðx  XðsÞÞ ds: ð4Þ

Thus, the indicator function can be obtained by solving the same type of Poisson equation as Eq.(3)with the special singular forcing term F(s) = n(s). In this paper, our goal is to consider the standard finite difference scheme with smoothing discrete delta function to solve the Eqs.(3) and (4), and to investigate its numerical accuracy.

3. One-dimensional analysis

Let us consider the one-dimensional counterpart as

d2u dx2¼ c

d

dxdðx 

a

Þ; 0 6 x 6 1; ð5Þ

with boundary conditions u(0) = u(1) = 0, and the immersed boundary point is located at x =

a

2 (0, 1). As is known, the exact solution u(x) can be represented as

uðxÞ ¼

Z 1

0

Gðx; yÞ c d

dydðy 

a

Þdy; ð6Þ

where G(x; y) is the well-known Green’s function defined as Gxx(x; y) = d(x  y), and can be explicitly written as

Gðx; yÞ ¼ xðy  1Þ; 0 6 x 6 y;

yðx  1Þ; y < x 6 1: 

For convenience, we set c = 1. By formally applying the integration by parts to Eq.(6), we obtain

uðxÞ ¼ 

Z 1

0

d

dyGðx; yÞdðy 

a

Þdy: ð7Þ

Now we specify a uniform grid with grid points xj= jh, j = 0, 1, . . . , N with h = 1/N, and then discretize the Eq.(5)with c = 1 by the standard centered difference scheme

Uj1 2Ujþ Ujþ1

h2 ¼

dhðxjþ1

a

Þ  dhðxj1

a

Þ

2h ; ð8Þ

where dhis the discrete delta function[8]defined as

dhðxÞ ¼ 1 4h 1 þ cosð px 2hÞ   ; if jxj 6 2h; 0; otherwise: ( ð9Þ

This discrete delta function satisfies

X

m¼N

m¼0

dhðxm

a

Þ h ¼ 1; ð10Þ

which is the corresponding basic requirement for the delta function. Other discrete delta functions can be found in[2,12]; however, the usage of other delta functions leads to the same conclusion as will be given in this paper. For clarity, we denote the first-order and second-order centered difference operators as Dhand D2h, respectively. Similar to the analytic solution in Eq.(6), the discrete solution Ujof Eq.(8)can also be written as

Uj¼ h

XN

m¼0

GjmDhdhðxm

a

Þ; ð11Þ

where Gjm= G(xj, xm) is the discrete version of Green’s function defined as

Gjm¼

xjðxm 1Þ; 0 6 j 6 m;

xmðxj 1Þ; m < j 6 N:



We can immediately check that G satisfies D2hGjm¼1hdjmwhere djmis the Kronecker delta function.

Now, by taking summation by parts and the property of the discrete delta function, the numerical solution Ujcan be rewritten as

Uj¼ h

XN

m¼0

DhGjmdhðxm

a

Þ: ð12Þ

Using the similar approach as in[11], for any discrete point xj, j = 0, 1, . . . , N, the point-wise error between Ujand u(xj) can be represented by

jUj uðxjÞj ¼ h XN m¼0 DhGjmdhðxm

a

Þ  Z 1 0 d dyGðxj;yÞdðy 

a

Þdy          6 h XN m¼0 DhGjmdhðxm

a

Þ      hX N m¼0 d dyGðxj;yÞ     y¼xm dhðxm

a

Þ     þ h XN m¼0 d dyGðxj;yÞ          _y¼x m dhðxm

a

Þ  d dyGðxj;yÞ     y¼a     ¼ E1þ E2:

Using the fact that the derivative of Green’s function is

d dyGðxj;yÞ     y¼xm ¼ xj; 0 6 xj6xm; xj 1; xm<xj61  ð13Þ

and that its discrete counterpart is

DhGjm¼ xj; j < m; xj1₂; j ¼ m; xj 1; j > m; 8 > < > : ð14Þ

we can conclude that the first part of the error E1becomes

E1¼ h 1 2dhðxj

a

Þ         ¼ Oð1Þ; when jxj

a

j 6 2h; 0; otherwise:  ð15Þ

The second part of the error E2is simply an interpolating error for the functiondydGðxj;yÞ   

y¼xm

. Using the formula in Eq.(13)

and the first moment condition in(10), since the discrete delta function has finite support 4h, we can obtain

E2¼ hP j1 m¼0 dhðxm

a

Þ        ; xj6

a

; hPN m¼j dhðxm

a

Þ          ; xj>

a

; 8 > > > > < > > > > :

which implies that

E2¼

Oð1Þ; as jxj

a

j 6 2h;

0; otherwise:



ð16Þ

From the above analysis, one can immediately see that the point-wise error appears only at some points around the sin-gular point

a

, which means that the maximum error kuh uk1is of order O(1). For the same reason, we can conclude that L1(kuh uk1) and L2(kuh uk2) errors are of order O(h) and O(h1/2), respectively. Our numerical results in next section will confirm this conclusion.

4. Numerical results 4.1. One-dimensional problem

In this subsection, we consider the following one-dimensional problem:

d2u dx2¼ c

d

dxdðx 

a

Þ þ g; 0 < x < 1; ð17Þ

with the interface at the point

a

=

p

/6. The exact solution is given as

uðxÞ ¼ x

3_{þ 2}

_a

_x2_{; if x 6}

_a

_;

7ðx3_{ 1Þ=3; if x >}

_a

_;

(

ð18Þ

where the jump size c of the solution u at the interface is set to be (2

a

3_{+ 7)/3. The regular source term g can be easily} com-puted by the analytic solution.

Throughout this section, the ratio between two consecutive errors is computed as kuuhk

kuuh=2k, where u is the exact solution and uhis the numerical solution with the mesh width h = 1/N. As the mesh is refined, the asymptotic ratios of 1, 1.414, 2 represent that the corresponding orders of accuracy are zeroth-order (log21), half-order log2

ffiffiffi 2 p

 

and first-order (log22).

Table 1shows the order of accuracy for our test which verifies exactly our one-dimensional analysis in the previous sec-tion. Notice that, when we implemented different versions of discrete delta function given in[2,12], we observed the same behavior of the errors.

4.2. Two-dimensional problems

For two-dimensional problem, we generally write the equation as

D

u ¼

r

 Z

C

FðsÞ d2_{ðx  XðsÞÞ ds þ g ð19Þ}

in a domainX= [a, b]  [c, d] with Dirichlet boundary conditions. We first divide the domainXinto M  N uniform grids with mesh widthDx =Dy = h, and ui,jdenotes the discrete solution at (xi, yj) where xi= ih, i = 0, 1, . . . , M and yj= jh, j = 0, 1, . . . , N. We also choose a collection of marker points X(sk) = (Xk, Yk) along the interfaceCwith the mesh points sk= kDs. Here, the marker mesh widthDs is about a half of h. Then we use the standard centered difference scheme to discrete Eq.(19)as

uiþ1;j 2ui;jþ ui1;j

h2 þ

ui;jþ1 2ui;jþ ui;j1

h2 ¼ fx iþ1=2;j f x i1=2;j h þ fy i;jþ1=2 f y i;j1=2 h þ gi;j; where fx iþ1=2;j¼ X k Fx ðskÞ dhðxiþ h=2  XkÞdhðyj YkÞ

D

s; f_i;jþ1=2y ¼X k FyðskÞ dhðxi XkÞdhðyjþ h=2  YkÞ

D

s and fx i1=2;j; f y

i;j1=2are defined in a similar fashion. The resultant matrix equation can be solved efficiently by the fast direct solver in Fishpack[1].

Example 1. For the first example, we test the accuracy of the indicator function which is obtained by solving Eq.(4). For completeness, we test three different interfaceCin the domain [1, 1]  [1, 1] as follows.

1. Cis a circle centered at (0, 0) with the radius 0.3.

2. Cis an ellipse centered at (0, 0) with the major radius 0.9 and minor radius 0.1. 3. Cis a simple closed curve written in polar coordinates: r = 0.5 + 0.25 cos(5h).

Tables 2–4show the errors and convergence ratios for these three different cases. One can immediately see that, whereas the indicator function does indeed have an O(1) error in maximum norm, it is first-order convergent in L1norm and half-order convergent in L2norm. The results are consistent to the one-dimensional analysis.

Example 2. In this example, we test an analytic solution which arises from the pressure Eq.(3)in Stokes flow developed in

[3]. This example is also used in [4]for a simple version of immersed interface method. The computational domain is

X= [2, 2]  [2, 2], and the interface is a unit circle centered at (0, 0), i.e., X(h) = (cos h, sin h). The exact solution is written in polar coordinates as uðr; hÞ ¼ r 3_{ðcosð3hÞ þ sinð3hÞÞ; if r 6 1;} r3_{ðcosð3hÞ  sinð3hÞÞ; if r > 1} ( Table 1

Order of accuracy for one-dimensional test.

N ku  uhk1 Ratio ku  uhk2 Ratio ku  uhk1 Ratio

32 8.8827E01 – 1.7057E01 – 4.2479E02 –

64 6.1709E01 1.4394 1.0736E01 1.5887 1.9004E02 2.2352

128 1.1847E00 0.5208 1.0708E01 1.0026 1.2980E02 1.4641

256 1.1579E00 1.0231 7.4120E02 1.4446 6.3791E03 2.0347

512 1.1030E00 1.0497 5.0171E02 1.4773 3.0741E03 2.0750

Table 2

Convergent test for indicator function in case 1: a circle.

M  N ku  uhk1 Ratio ku  uhk2 Ratio ku  uhk1 Ratio

32  32 3.6463E01 – 1.3162E01 – 6.3848E02 –

64  64 4.5555E01 0.8004 9.5529E02 1.3777 3.2182E02 1.9839

128  128 4.8736E01 0.9347 7.2764E02 1.3128 1.6837E02 1.9113

256  256 4.8610E01 1.0026 4.9738E02 1.4629 8.2361E03 2.0443

and the boundary force F(h) = 2 sin(3h)(X0_{(h) + X(h)).Table 5shows the convergence ratios which again verifies our} one-dimensional analysis.

Example 3. Finally we consider Eq.(19)for which the domain isX= [1, 1]  [1, 1] and the interfaceCis described as a simple closed curve r = 0.5 + 0.25 cos(5h) in polar coordinates. The analytic solution u is given by

u ¼ ðx 2_{ 1Þðy}2_{ 1Þ þ 1; if ðx; yÞ 2}

_X

out; ðx2_{ 1Þðy}2_{ 1Þ; if ðx; yÞ 2}

_X

in; (

whereXoutandXinrepresent the exterior and interior regions, respectively. Note that the boundary force F is simply the normal vector n along the interfaceC. The external source g can be easily computed from the exact solution u.

Table 6shows the convergence tests of the numerical solutions forExample 3. The results are consistent to what we expect, i.e. they are first-order accurate in L1norm, half-order accurate in L2norm, but have O(1) errors in L1norm.

Acknowledgment

Y. Kim is supported by the National Research Foundation of Korea, Grant 2010-0006165. M.-C. Lai is supported in part by National Science Council of Taiwan under research Grant NSC-97-2628-M-009-007-MY3, NSC-98-2115-M-009-014-MY3, and the support of NCTS in Taiwan.

Table 3

Convergent test for indicator function in case 2: an ellipse.

M  N ku  uhk1 Ratio ku  uhk2 Ratio ku  uhk1 Ratio

32  32 6.7302E01 – 1.9248E01 – 1.2276E01 –

64  64 5.0021E01 1.3454 1.4391E01 1.3374 6.5139E02 1.8845

128  128 4.9922E01 1.0019 9.7919E02 1.4697 3.1954E02 2.0385

256  256 4.9834E01 1.0017 6.8903E02 1.4211 1.5951E02 2.0033

512  512 4.9617E01 1.0043 4.8685E02 1.4152 7.9510E03 2.0061

Table 4

Convergent test for indicator function in case 3: a simple closed curve.

M  N ku  uhk1 Ratio ku  uhk2 Ratio ku  uhk1 Ratio

32  32 5.9986E01 – 2.5162E01 – 2.0860E01 –

64  64 5.5492E01 1.0810 1.8259E01 1.3780 1.0827E01 1.9266

128  128 5.3029E01 1.0464 1.2910E01 1.4143 5.4431E02 1.9891

256  256 5.1669E01 1.0263 9.0547E02 1.4257 2.7064E02 2.0111

512  512 5.1194E01 1.0092 6.4251E02 1.4092 1.3547E02 1.9977

Table 5

Convergent test forExample 2.

M  N ku  uhk1 Ratio ku  uhk2 Ratio ku  uhk1 Ratio

32  32 1.5643E00 – 9.5960E01 – 2.0421E00 –

64  64 1.7182E00 0.9104 6.9177E01 1.3871 1.1867E00 1.7209

128  128 1.8342E00 0.9367 4.9447E01 1.3990 6.4868E01 1.8292

256  256 1.9086E00 0.9610 3.4999E01 1.4128 3.4044E01 1.9054

512  512 1.9284E00 0.9897 2.4775E01 1.4126 1.7495E01 1.9459

Table 6

Convergent test forExample 3.

M  N ku  uhk1 Ratio ku  uhk2 Ratio ku  uhk1 Ratio

32  32 5.9986E01 – 2.5162E01 – 2.0860E01 –

64  64 5.5492E01 1.0809 1.8259E01 1.3780 1.0827E01 1.9266

128  128 5.3029E01 1.0464 1.2910E01 1.4143 5.4431E01 1.9891

256  256 5.1669E01 1.0263 9.0547E02 1.4257 2.7064E02 2.0111

References

[1] J. Adams, P. Swarztrauber, R. Sweet, Fishpack – A package of Fortran subprograms for the solution of separable elliptic partial differential equations. Available from:<http://www.netlib.org/fishpack>.

[2] R.P. Beyer, R.J. LeVeque, Analysis of a one-dimensional model for the immersed boundary method, SIAM J. Numer. Anal. 29 (1992) 332–364. [3] R. Cortez, The method of regularized Stokeslets, SIAM J. Sci. Comput. 23 (2001) 1204–1225.

[4] M.-C. Lai, H.-C. Tseng, A simple implementation of the immersed interface methods for Stokes flows with singular forces, Comput. Fluids 37 (2008) 99– 106.

[5] Z. Li, K. Ito, The Immersed Interface Method, SIAM, 2006.

[6] R.J. LeVeque, Z. Li, The immersed interface method for elliptic equations with discontinuous coefficients and singular sources, SIAM J. Numer. Anal. 31 (1994) 1019–1044.

[7] Y. Mori, Convergence proof of the velocity field for a Stokes flow immersed boundary method, Commun. Pure Appl. Math. LXI (2008) 1213–1263. [8] C.S. Peskin, Numerical analysis of blood flow in the heart, J. Comput. Phys. 25 (1977) 220–252.

[9] C.S. Peskin, The immersed boundary method, Acta Numerica 11 (2002) 479–517.

[10] S.O. Unverdi, G. Tryggvason, A front-tracking method for viscous incompressible multi-fluid flows, J. Comput. Phys. 100 (1992) 25–37. [11] A.-K. Tornberg, B. Engquist, Numerical approximations of singular source terms in differential equations, J. Comput. Phys. 200 (2004) 462–488. [12] X. Yang, X. Zhang, Z. Li, G.-W. He, A Smoothing technique for discrete delta functions with application to immersed boundary method in moving