A simple implementation of the immersed interface methods for Stokes flows with singular forces

(1)

A simple implementation of the immersed interface methods

for Stokes ﬂows with singular forces

Ming-Chih Lai

*

_{, Hsiao-Chieh Tseng}

Department of Applied Mathematics, National Chiao Tung University, 1001, Ta Hsueh Road, Hsinchu 300, Taiwan Received 11 November 2006; received in revised form 15 April 2007; accepted 27 April 2007

Available online 21 June 2007

Abstract

In this paper, we introduce a simple version of the immersed interface method (IIM) for Stokes flows with singular forces along an interface. The numerical method is based on applying the Taylor’s expansions along the normal direction and incorporating the solution and its normal derivative jumps into the finite difference approximations. The fluid variables are solved in a staggered grid, and a new accurate interpolating scheme for the non-smooth velocity has been developed. The numerical results show that the scheme is second-order accurate.

1. Introduction

The fluid problems with interfaces have many applica-tions in science and engineering. For example, Peskin’s for-mulation for the blood flow in the heart valve leaflet (interface) involves exerting singular forces along the leaflet into the blood (fluid) [19]. Another example is the Hele-Shaw flow which can be formed by pumping a more viscous fluid into a less viscous surrounding fluid in two parallel thin plates. The shape of the interface is known to exhibit a fingering phenomenon. The governing equa-tions for the above problems involve solving elliptic partial differential equations with possible discontinuous coeffi-cients or with possible jump conditions for the unknown and its derivative across the interface. Several classes of practical finite difference methods have been developed in the past three decades. Those methods provide different treatments near the interface while keep the standard differ-ence discretization away from the interface. We summarize those different techniques as follows.

The immersed boundary (IB) method of Peskin[20] is the ﬁrst Cartesian grid method to provide a simple solution

to the fluid problems interacting with complicated inter-faces. The method is to treat the interface as a singular force (source) generator along which the force can be trans-mitted into the fluid smoothly. More precisely, the method uses the Lagrangian markers to track the interface and affects the fluid in Eulerian grid via a smooth version of dis-crete delta function. Since the equations are solved on a Cartesian grid, many fast fluid solvers can be applied easily. For instance, the original IB method exploits Fast Fourier Transform (FFT) in solving the Stokes equations in rectangular domain. The method is easy to implement but it is only first-order accurate. Besides, the solution is smoothing out across the interface. Although there are dif-ferent variants of formally second-order schemes recently

[9,7], the overall accuracy can be increased only for the suf-ficiently smooth problems. Despite this first-order accuracy issue, the IB method nowadays is getting more attention on applying to the flow interacting with rigid boundaries which has important applications in computational fluid dynamics. The reason for that is quite simple since it does not need to handle grid generations for unsteady problems which can save the computational effort significantly.

Unlike the IB method to have numerical smearing near the interface, the immersed interface method was invented to capture the solution and its derivative jumps sharply.

*

Corresponding author.

E-mail address:mclai@math.nctu.edu.tw(M.-C. Lai).

www.elsevier.com/locate/compﬂuid Computers & Fluids 37 (2008) 99–106

(2)

The idea is to incorporate those jump conditions into the finite difference discretization by modifying the difference approximations near the interface. Mayo[18]used the idea to solve the Poisson equation in an irregular domain based on an integral equation. Later, LeVeque and Li applied the similar idea to an elliptic interface problem with discontin-uous coefficients and singular sources[13]. Since the prob-lem is more difficult to solve than the Poisson probprob-lem, they introduce the local coordinates and derive more inter-face relations so that the finite difference scheme can match up to all second-order partial derivatives. Their method appears to be the first finite difference method that is sec-ond-order accurate for the problem. Dumett and Keener has extended the method to the anisotropic case of elliptic interface problems[6]. There are other applications to the fluid problems [12,16,15,23], just to name a few. We refer the interested readers to Li’s recent review article[10] and the newly published book[11].

There are several methods sharing the similar spirit as IIM in literature. For example, the boundary condition capturing method proposed by Liu et al.[17]can be imple-mented to solve the elliptic interface problems dimension by dimension. The method captures the solution and its normal derivative jumps sharply while smoothing the tan-gential derivative. The method is easy to implement; how-ever, it is only ﬁrst-order accurate. There are other variants of IIM such as the explicit jump [24]and the decomposed immersed interface method [4]. Beale and Lai [2] applied the boundary integral method to compute the solution for the Laplace problem near the interface and used it to form the correction terms. This method has been further extended to the Stokes problem with singular forces [3].

In this paper, we extend the immersed interface method developed in[21]for Poisson problems to Stokes problems with singular forces along an interface. The numerical method is based on applying the Taylor’s expansions along the normal direction and incorporating the solution and its normal derivative jumps into the finite difference approxi-mations. Unlike the traditional IIM, the present fluid vari-ables are solved in a staggered grid [8] instead. Another contribution of this work is to derive a new accurate inter-polation scheme for the non-smooth velocity field. The numerical results show that the scheme is indeed second-order accurate.

The paper is organized as follows. In Section 2, we review the immersed interface methods for Poisson prob-lems with an interface and introduce a new interpolation scheme from a grid function to the interface. Then the numerical scheme for Stokes ﬂows with singular forces is discussed and the numerical results are shown in Section

3. Some conclusions are given in Section4.

2. Poisson equations with jump discontinuities across the interface

In this section, we consider the Poisson problem on a computational domain X¼ ½a; b ½c; d with an immersed

interface C¼ fXðsÞ ¼ ðX ðsÞ; Y ðsÞÞ; 0 6 s < Lg, where s is the arc-length parameter and Xð0Þ ¼ XðLÞ. The interface C divides the domain X into two regions,; namely inside (X) and outside (Xþ) the interface. Across C, the solution and its normal derivative exhibit jump discontinuities. The Poisson problem can be written as

Du¼ f in X; ð2:1Þ ½uðsÞ ¼ gðsÞ; ou on ðsÞ ¼ hðsÞ on C; ð2:2Þ u¼ ub on oX; ð2:3Þ

where the right-hand side function f could be discontinuous across the interface C. The jumps ½uðsÞ ¼ uþ_u _and

ou on ðsÞ ¼ouþ on ou

on are deﬁned as the diﬀerence of the

limit-ing values from two diﬀerent sides of the interface. Here, we also assume the interface parametrization XðsÞ and the jump discontinuity functions g(s) and h(s) are fairly smooth (say C2ðC)) on the interface.

2.1. An IIM to incorporate the jumps in the normal direction In this subsection, we review the scheme developed in

[21] to solve Eqs. (2.1)–(2.3). Before we proceed, let us lay out a uniform Cartesian grid in X with mesh width h¼ Dx ¼ Dy, and let the discretized solution at the grid point xi;j¼ ðxi; yjÞ ¼ ða þ iDx; c þ jDyÞ be denoted by ui;j.

As in [18], we classify the grid point as either a regular or irregular point. For a regular grid point, we mean that the interface does not cut through the standard five-point Laplacian of that point. On the other hand, if the five-point Laplacian of a grid five-point involves using the grid points inside and outside the interface simultaneously, then we call such point as an irregular point. Since the solution is smooth either inside or outside the interface, the five-point Laplacian of a regular point does not need to modify in order to have the second-order accuracy. However, it has to be modified at an irregular point since the solution is not smooth across the interface, and the modification depends on the jumps [u] and ou

on

. We explain how to modify the discrete Laplacian in the following.

Let xi;j be an irregular grid point in X as shown in

Fig. 1. The ﬁve-point Laplacian at xi;j can be written as

(3)

Dhuðxi;jÞ ¼

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

h2 þ

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

h2 ¼u

þ

i1;j 2ui;jþ uiþ1;j

h2 þ uþ_i;jþ1 2u i;jþ ui;j1 h2 ¼u

i1;j 2ui;jþ uiþ1;j

h2 þ u_i;jþ1 2u i;jþ ui;j1 h2 þu þ i1;j ui1;j h2 þ uþ i;jþ1 ui;jþ1 h2 ¼ uxxðxi;jÞ þ uyyðxi;jÞ þ Oðh2Þ þ

uc i1;j h2 þ uc i;jþ1 h2 ¼ fi;jþ 1 h2 u c i1;jþ uci;jþ1 þ Oðh2_Þ: _ð2:4Þ

To derive the correction term uc

i1;j, we need to ﬁnd the

orthogonal projection of xi1;j at the interface (say X¼

XðsÞ in Fig. 1), and then apply the Taylor’s expansions

along the normal direction at X. (Note that, there

are two diﬀerent techniques of ﬁnding the orthogonal projection of a grid point on the interface based on the interface representation; namely, the arc-length parametri-zation and the level set representation. We refer the inter-ested readers to the review article[10].) So we have

uc

i1;j¼ uþi1;j ui1;j

¼ uþ þ dou þ on þ d2 2 o2uþ on2 þ Oðh 3 Þ u þ d ou on þ d2 2 o2u on2 þ Oðh 3_Þ ¼ ½u_Xþ d ou on X þd 2 2 o2u on2 X þ Oðh3_Þ: _ð2:5Þ

The value d is the signed distance between the grid point xi1;j and the orthogonal projection X. Notice that, if

the grid point xi1;j is inside the interface, then d must be

negative.

In Eq.(2.5), the jumps½uX and

ou on

X are simply given

from Eq. (2.2); however, the second normal derivative jump o2u

on2

h i

X

is not known. To ﬁnd the second normal derivative jump, one can ﬁrst rewrite the Laplace operator in Eq.(2.1)on the interface as in[22]

o2u on2ðXðsÞÞ þ jðXðsÞÞ ou onðXðsÞÞ þ o2u os2ðXðsÞÞ ¼ f ðXðsÞÞ; ð2:6Þ where the value jðXðsÞÞ is the local curvature of the inter-face. Then, we can easily compute the second normal deriv-ative jump at X by o2u on2 X ¼ ½f _X jX ou on X o 2 os2½uX: ð2:7Þ

Finally, the correction term can be approximated by

uc_i1;j¼ ½uXþ d ou on X þd 2 2 o2u on2 X ¼ ½u_Xþ d ou on X þd 2 2 ½f X jX ou on X o 2 os2½uX ! : ð2:8Þ The correction term uc

i;jþ1can be evaluated in a similar way.

It is worth mentioning that the correction term uc i1;j

must be computed up to third-order accurate since it appears in the discrete Laplacian term, that is, a term of Oð1=h2Þ. Thus, the local truncation error at an irregular point is OðhÞ. It has been proved by diﬀerent authors

[3,11] that despite the fact of the first-order truncation errors at those irregular points (but second-order accurate at regular points), the overall accuracy is still second-order. Since the correction terms are only involved some modifica-tions in the right-hand side of the finite difference equamodifica-tions, the resultant linear system(2.4)can be efficiently solved by a fast Poisson solver such as Fishpack[1]. The computational cost spent on the modifications at irregular grid points is just a small portion of that for a fast Poisson solver.

2.2. An interpolation scheme from a grid function to the interface

In the Stokes problems, once we have obtained the veloc-ity field on the grid, we need to find the velocveloc-ity on the inter-face. Thus, a comparably accurate interpolation formula from a grid function to the interface markers needs to be developed. When the velocity field is smooth, the standard bi-linear interpolation can be used to achieve the second-order accuracy. However, if the velocity is not smooth across the interface (½u ¼ 0; ou

on

6¼ 0, as in the next section), the bi-linear formula needs to be modified in order to have second-order accurate interpolations. In the following, we derive such modification for the case that even the jump of u is not zero. Notice that, this new scheme differs from the one used in[14]where the authors used a linear combi-nation of the three nearby grid values with a correction term that consists of the function jump and the partial derivative jumps in both x- and y-directions at the interpolated mar-ker. Here, however, we simply modify the standard bi-linear interpolation formula (involving four neighboring grid val-ues) with the incorporation of the function and its normal derivative jumps which are given directly from the problem. The detail can be explained as follows.

Suppose we want to interpolate the value of u

I at a

mar-ker XI¼ ðXI; YIÞ ¼ ðxi;jþ ah; yi;jþ bhÞ from the inside of

the interface using surrounding values ui;j; uiþ1;j; ui;jþ1 and

uiþ1;jþ1 as illustrated inFig. 2. Let CI be the interpolating

correction term that should be added to the standard bi-lin-ear formula so that the interpolation accuracy is second-order. That is

u_I ¼ uBIþ CI þ Oðh

(4)

where uBI is the bi-linear interpolation of u at XI. To derive C_I, we have uBI¼ buTþ ð1 bÞuB ¼ bðð1 aÞuþ i;jþ1þ auiþ1;jþ1Þ þ ð1 bÞðð1 aÞu i;jþ au iþ1;jÞ ¼ bðð1 aÞu

i;jþ1þ auiþ1;jþ1Þ þ ð1 bÞðð1 aÞui;j

þ au iþ1;jÞ þ bð1 aÞu c i;jþ1 ¼ uI þ bð1 aÞu c i;jþ1;

where the correction uc

i;jþ1 can be derived exactly the same

procedure as Eq. (2.5). However, since we only need to approximate the correction up to Oðh2Þ, one can neglect the second normal derivative jump, that is

uc i;jþ1¼ u þ i;jþ1 u i;jþ1¼ ½uXþ d ou on X þ Oðh2Þ; ð2:10Þ

where X is the orthogonal projection of the grid point

xi;jþ1 on the interface. Therefore, the interpolating

correc-tion term can be derived as

C_I ¼ bð1 aÞ ½u_Xþ d ou on X ! : ð2:11Þ 2.3. Numerical examples

We consider three diﬀerent exact solutions (shown in

Table 1)for Poisson equation in X¼ ½1; 1 ½1; 1 with jump discontinuities on an elliptical interface C¼ fx2

a2þ

y2

b2¼ 1g with a ¼ 0:8, and b ¼ 0:2. In [21], the

similar examples were used but with a circle interface so that the orthogonal projection point can be found straight-forwardly and the implementation is much simpler. In our present code, we have also implemented the uniformly

cho-sen markers to construct the interface and thus from which to ﬁnd the necessary orthogonal projections and the curva-tures. Our numerical scheme and interface construction can be applied to any smooth curve in principle.

Table 2shows the maximum errors of N N grid points in X. ForExamples 1 and 2, we can see that the present scheme is indeed second-order convergent. For Example 3, the discrete Laplacian approximates the Laplacian oper-ator exactly so that the error is roughly equal to the machine precision. Those numerical results behave simi-larly as the ones obtained in[21]for the circular interface problems. Throughout this paper, all of numerical results were obtained by running our Fortran 90 code in a PC of Intel Pentium 4 (3.00 GHz) processor with 4 GB RAM. We also list the CPU time (s) in Table 2 for our computations. One can see that even a grid size 320· 320, the overall CPU time takes less than a quarter of a second.

Once we obtain the solution on the grid, we can interpo-late the value on the interface from inside by using the tech-nique discussed before. Note that, the solution and its normal derivative are not continuous across the interface C. Table 3shows the maximum errors of the interpolating solution (from inside) ofExample 2on the interface mark-ersðXk; YkÞ ¼ ð0:8 cos hk;0:2 sin hkÞ; k ¼ 0; 1; . . . ; N 1 with

hk¼ kDh, Dh ¼ 2p=N . One can see the present

interpolat-ing accuracy is also close to second-order.

Fig. 2. A diagram of the interpolation from a grid function to an interface marker XI.

Table 1

Three test examples for Poisson problem (2.1)–(2.3) with an elliptic interface x2

ð0:8Þ2þ y2

ð0:2Þ2¼ 1 in X ¼ ½1; 1 ½1; 1

Example 1 Example 2 Example 3

u 1 expðxÞ cosðyÞ x2 y2 uþ 1þ lnð2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x2_{þ y}2 p Þ expðx2_{Þ cosðyÞ} ₀ f 0 0 0 fþ 0 ð1 þ 4x2Þ expðx2Þ cosðyÞ 0 Table 2

The maximum errors for the examples of Poisson problem on the N N grid points

N ku uek1 Ratio CPU time

Example 1 40 2.1577E03 – 0.015 80 6.3698E04 3.39 0.016 160 1.7153E04 3.71 0.046 320 4.0663E05 4.22 0.219 Example 2 40 1.1909E03 3.71 0.015 80 3.0901E04 3.85 0.016 160 7.8497E05 3.94 0.078 320 1.9776E05 3.97 0.235 Example 3 40 5.5511E16 – 0.015 80 7.9797E16 – 0.016 160 2.6749E15 – 0.047 320 1.5377E14 – 0.172 Table 3

The interpolation errors for Example 2 of Poisson problem on the interface markers ðXk; YkÞ ¼ ð0:8 cos hk;0:2 sin hkÞ; k ¼ 0; 1; . . . ; N 1

with hk¼ kDh, Dh ¼ 2p=N N ku_u ek1;C Ratio 40 1.6036E03 – 80 4.7394E04 3.38 160 1.2650E04 3.75 320 4.0435E05 3.13

(5)

3. Stokes ﬂow with singular forces

In this section, we consider the following Stokes equations:

rp þ lDu þ f þ g ¼ 0 in X; ð3:1Þ r u ¼ 0 in X; ð3:2Þ

u¼ ub onoX; ð3:3Þ

where u¼ ðuðxÞ; vðxÞÞ is the ﬂuid velocity, p ¼ pðxÞ is the pressure, and l is the constant viscosity. There are two dif-ferent forces acting to the ﬂuid; namely, the singular force f exerted by the immersed interface (immersed boundary) and the external force g (could be discontinuous). More precisely, the singular force f has a support only along the interface C¼ fXðsÞ ¼ ðX ðsÞ; Y ðsÞÞ; 0 6 s 6 Lg and it has the form

fðxÞ ¼ Z

C

FðsÞdðx XðsÞÞds; ð3:4Þ

where d is the 2D Dirac delta function. The boundary force F deﬁned on the interface can be further decomposed to its normalðn ¼ ðn1; n2ÞÞ and tangential (s ¼ ðs1;s2Þ) directions

as

F¼ ðF1; F2Þ ¼ ðF nÞn þ ðF sÞs ¼ Fnnþ Fss; ð3:5Þ

where Fn and Fs represent the normal and tangential

forces, respectively.

Since the force f has a delta function singularity along the interface, one can expect that the velocity and the pres-sure will not be smooth across the interface. In fact, the boundary force F plays an important role to the jump con-ditions of u and p, and their derivatives. In IIM, we usually reformulate the Eqs. (3.1)–(3.5) into a Stokes problem without the delta function force term but with some known jump conditions across the interface. Thus, the governing equations are summarized as follows[14]:

rp þ lDu þ g ¼ 0 in X; ð3:6Þ r u ¼ 0 in X; ð3:7Þ ½u ¼ 0; l ou on ¼ Fss on C; ð3:8Þ ½p ¼ Fn; op on ¼oFs os þ ½g n on C; ð3:9Þ u¼ ub onoX: ð3:10Þ 3.1. Numerical scheme

Unlike the traditional IIM implementation, we use a staggered grid [8] with a uniform mesh width h¼ Dx ¼ Dy for the ﬂuid velocity and the pressure as in

Fig. 3. By applying the divergence operator to Eq. (3.6)

and using the incompressibility constraint of Eq. (3.7), the pressure satisﬁes the following equation:

Dp¼ r g in X; ð3:11Þ ½p ¼ Fn; op on ¼oFs os þ ½g n on C: ð3:12Þ This is exactly the type of Poisson equation that we have discussed in the previous section so the scheme discussed before can be applied directly. One remaining question is how to choose the pressure boundary condition on the computational boundaryoX. In this work, for testing pur-pose, the pressure boundary condition can be chosen either a Dirichlet (Examples 1 and 2in next subsection) or a zero Neumann boundary condition (Example 3 in next subsection).

Once we solve the pressure, the velocity ﬁeld can be obtained by solving other two Poisson-type equations with the interface as Du¼1 lðrp gÞ in X; ð3:13Þ ½u ¼ 0; ou on ¼ 1 lFss on C; ð3:14Þ u¼ ub on oX: ð3:15Þ

Once again, the numerical scheme in previous section can also be applied directly to ﬁnd the velocity. However, in order to compute the pressure gradient at the right-hand side of Eq.(3.13)accurately, one has to take its jump across the interface into account. To derive the jumps [px] and [py],

we ﬁrst recall that the unit tangent vector can be written as s¼ ðs1;s2Þ ¼ ðn2; n1Þ at each point on C. Since

op on ¼ ½pxn1þ ½pyn2¼ oFs os þ ½g n; ð3:16Þ o os½p ¼ ½pxs1þ ½pys2¼ ½pxn2þ ½pyn1¼ oFn os ; ð3:17Þ we can easily obtain the pressure gradient jump ½rp ¼ ð½px; ½pyÞ as

Fig. 3. A diagram of the interface cutting through a staggered grid with a uniform mesh width h in X. The pressure is located at the center of the cell, the velocity component u is at the left-right face, and v is at the top-bottom face of a cell.

(6)

½px ¼ n1 oFs os þ ½g n n2 oFn os ; ð3:18Þ ½py ¼ n2 oFs os þ ½g n þ n1 oFn os : ð3:19Þ We summarize our numerical algorithm as follows: Step 1: Solve the pressure Poisson equation Dhp¼ r g þ

pc

h2 in X ð3:20Þ

with modiﬁed right-hand side functions at those irregular points which the modiﬁcations are based on the jump con-ditions of the pressure. More precisely, the correction terms at those irregular points can be computed as

pc¼ ½p þ a op on þ1 2a 2 _{½Dp j} op on o 2 ½p os2 ¼ Fnþ a oFs os þ ½g n þ1 2a 2 _{½r g j} oFs os þ ½g n o 2_F n os2 : ð3:21Þ

Step 2: Solve the velocity equations

Dhu¼ 1 lðrhp gÞ þ uc h2 in X; ð3:22Þ u¼ ub on oX; ð3:23Þ

where the corrections are made only at irregular points as

uc¼ ½u þ a ou on þ1 2a 2 Du ½ j ou on o 2 ½u os2 ¼1 l aFssþ 1 2a 2_{ð½rp ½g þ jF} ssÞ : ð3:24Þ Note that, the above pressure gradient jump has been derived from(3.18) and (3.19).

The remaining question is how to approximate the pres-sure gradient in Eq. (3.22). At those regular points, the derivatives px and py at the locations of ui;j and vi;j can

be approximated by centered diﬀerences ðp_iþ1;j pi;jÞ=h,

and ðp_i;jþ1 pi;jÞ=h, respectively. However, a correction

term must be added when the p-grid points straddle across the interface. Suppose the locations of pi;jand ui;jfall in the

same side of the interface while the location of p_iþ1;jfalls in another side of the interface as illustrated inFig. 3, then we need to add a correction term pc

iþ1;j of piþ1;j so that the

approximation of px at ui;j grid can be computed by

ððpiþ1;jþ pciþ1;jÞ pi;jÞ=h. Note that, the correction term

pc

iþ1;jcan be computed by the formula(3.21). On the other

hand, if the locations of p_iþ1;jand ui;jfall in the same side of

the interface while the location of p_i;jfalls in another side of the interface, then the approximation of pxat ui;j grid can

be computed byðpiþ1;j ðpi;jþ pci;jÞÞ=h, where a correction

term pc

i;jof pi;jmust be added. The approximation of pyat

the vi;j grid point can be handled in a similar manner.

3.2. Numerical results

In this subsection, we perform three diﬀerent numerical tests to examine the accuracy of the present numerical schemes for Stokes problems. Throughout the following examples, the computational domain is chosen as X¼ ½2; 2 ½2; 2, and the interface is simply a circle with radius one so that it can be represented by C¼ fXðhÞ ¼ ðcosðhÞ; sinðhÞÞ; 0 6 h < 2pg where h is the arc-length parameter. For simplicity, we also choose the viscosity l = 1.

Example 1 (Normal force on an interface[5]). In the ﬁrst example, we consider the case in which the boundary force exerts only along the normal direction, and it has the form as

FðhÞ ¼ 2 sinð3hÞXðhÞ: ð3:25Þ For convenience, the pressure and the velocity are written in polar coordinates as pðr; hÞ ¼ r 3_{sinð3hÞ; r < 1;} r3_sinð3hÞ; _{r >}_1; ( uðr; hÞ ¼ 3 8r 2_{sinð2hÞ þ} 1 16r 4_sinð4hÞ 1 4r 4_sinð2hÞ; _{r <}_1; 1 8r 2_sinð2hÞ 3 16r 4_sinð4hÞ þ1 4r 2_sinð4hÞ; _r_P_1; 8 > > > > > < > > > > > : vðr; hÞ ¼ 3 8r 2_cosð2hÞ 1 16r 4_cosð4hÞ 1 4r 4_cosð2hÞ; _{r <}_1; 1 8r 2_{cosð2hÞ þ} 3 16r 4_cosð4hÞ 1 4r 2_cosð4hÞ; _r_P_1: 8 > > > > > < > > > > > :

Example 2 (Tangential force on an interface [5]). In this example, we consider the case in which the boundary force exerts only along the tangential direction, and it has the form as

FðhÞ ¼ 2 sinð3hÞX0_ðhÞ: _ð3:26Þ

The pressure and the velocity are given by pðr; hÞ ¼ r 3_cosð3hÞ; _{r <}_1; r3_{cosð3hÞ; r > 1;} uðr; hÞ ¼ 1 8r 2_{cosð2hÞ þ} 1 16r 4_cosð4hÞ 1 4r 4_cosð2hÞ; _{r <}_1; 1 8r 2_{cosð2hÞ þ} 5 16r 4_cosð4hÞ 1 4r 2_cosð4hÞ; _r_P_1; 8 > > > < > > > : vðr; hÞ ¼ 1 8r 2_{sinð2hÞ þ} 1 16r 4_sinð4hÞ þ1 4r 4_sinð2hÞ; _{r <}_1; 1 8r 2_{sinð2hÞ þ} 5 16r 4_sinð4hÞ 1 4r 2_sinð4hÞ; _r_P_1: 8 > > > < > > > :

(7)

After some careful calculations, one can obtain that the velocity and the pressure in above two examples satisfy the Stokes equations(3.6)–(3.10) without the external forcing term (g¼ 0). The pressure and velocity boundary condi-tions for those two examples are all Dirichlet.

Example 3 (Normal and tangential forces on an inter-face). In this example, we consider the boundary force that has both normal and tangential components as

FðhÞ ¼ 2 sinð3hÞX0ðhÞ cos3_ðhÞXðhÞ: _ð3:27Þ

The velocity components u and v are chosen exactly the same as inExample 2while the pressure written in Carte-sian coordinates (x¼ r cos h; y ¼ r sin h) is

pðx; yÞ ¼ x

3_{þ cosðpxÞ cosðpyÞ; r < 1;}

cosðpxÞ cosðpyÞ; rP1:

The velocity and the pressure satisfy the Stokes equations

(3.6)–(3.10) with a nonzero external force g6¼ 0 which can be calculated analytically. Here, we use the zero Neu-mann boundary condition (easily to be checked) for the pressure and the Dirichlet boundary conditions for the velocity.

Table 4 shows the maximum errors of the computed velocity u; v and the pressure p for those three examples. One can see from these numerical results that our scheme is indeed second-order accurate. In addition, using the computed values on the grid, we can interpolate the velocity from the grid to the interface. Table 5 shows

the interpolation errors for the velocity ﬁeld on the interface. One can also see that the interpolating accuracy is second-order on average.

4. Conclusion

In this paper, we introduce a simple version of the immersed interface methods (IIM) for Stokes flows with singular forces along an interface. The numerical method is based on applying the Taylor’s expansions along the normal direction and incorporating the solution and its normal derivative jumps into the finite difference approxi-mations. The fluid variables are solved in a staggered grid, and a new interpolation scheme for the non-smooth veloc-ity has been developed. The numerical results show that the scheme is second-order accurate. In the future work, we shall develop a fourth-order accurate scheme to the present problems by taking more jump relations into account and extend the same techniques to the three-dimensional problems.

Acknowledgments

The authors were supported in part by MOE-ATU program and National Science Council of Taiwan under research grant NSC-94-2115-M-009-004.

References

[1] Adams J, Swarztrauber P, Sweet R. Fishpack—A package of Fortran subprograms for the solution of separable elliptic partial diﬀerential equations. Available from:http://www.netlib.org/ﬁshpack.

[2] Beale JT, Lai M-C. A method for computing nearly singular integrals. SIAM J Numer Anal 2001;38:1902–25.

[3] Beale JT, Layton AT. On the accuracy of ﬁnite diﬀerence methods for elliptic problems with interfaces. Commun Appl Math Comput Sci 2006;1:91–119.

[4] Berthelsen PA. A decomposed immersed interface method for variable coeﬃcient elliptic equations with non-smooth and discon-tinuous solutions. J Comput Phys 2004;197:364–86.

Table 4

The maximum errors for the examples of Stokes problem on N N grid points

N ku uek1 Ratio kv vek1 Ratio kp pek1 Ratio

Example 1

32 2.9955E03 – 9.5555E03 – 1.4625E02 –

64 7.4576E04 4.02 2.1775E03 4.39 3.2027E03 4.57

128 2.1442E04 3.48 5.4344E04 4.01 8.2001E04 3.91

256 4.8445E05 4.43 1.3800E04 3.94 1.9358E04 4.24

Example 2

32 9.3164E03 – 5.5489E03 – 1.7579E02 –

64 2.2334E03 4.17 9.8214E04 5.65 3.5421E03 4.96

128 4.5329E04 4.93 2.6948E04 3.64 9.5814E04 3.70

256 1.2100E04 3.75 6.8943E05 3.91 2.1994E04 4.36

Example 3

32 9.9654E03 – 9.3837E03 – 2.5682E02 –

64 2.7483E03 3.63 1.8844E03 4.98 7.2394E03 3.55

128 5.2897E04 5.20 4.4803E04 4.21 1.8827E03 3.85

256 1.4410E04 3.67 1.2263E04 3.65 4.7359E04 3.98

Table 5

The interpolation errors forExample 3of Stokes problem on the inter-face markers ðXk; YkÞ ¼ ðcos hk;sin hkÞ; k ¼ 0; 1; . . . ; N 1, hk¼ kDh,

Dh¼ 2p=N N ku uek1;C Ratio kv vek1;C Ratio 32 1.0035E02 – 1.0923E02 – 64 2.3020E03 4.36 2.9853E03 3.66 128 4.5430E04 5.07 6.8889E04 4.33 256 1.2788E04 3.55 1.8553E04 3.71

(8)

[5] Cortez R. The method of regularized stokeslets. SIAM J Sci Comput 2001;23:1204–25.

[6] Dumett M, Keener JP. A numerical method for solving anisotropic elliptic boundary value problems in 3D. SIAM J Sci Comput 2003;25:348–67.

[7] Griﬃth BE, Peskin CS. On the order of accuracy of the immersed boundary method: higher order convergence rates for suﬃciently smooth problems. J Comput Phys 2005;208:75–105.

[8] Harlow FH, Welsh JE. Numerical calculation of time-dependent viscous incompressible ﬂow of ﬂuid with a free surface. Phys Fluids 1965;8:2181–9.

[9] Lai M-C, Peskin CS. An immersed boundary method with formal second-order accuracy and reduced numerical viscosity. J Comput Phys 2000;160:705–19.

[10] Li Z. An overview of the immersed interface method and its applications. Taiwan J Math 2003;7:1–49.

[11] Li Z, Ito K. The immersed interface method—numerical solutions of PDEs involving interfaces and irregular domains. SIAM Frontiers Appl Math 2006;33.

[12] Li Z, Lai M-C. The immersed interface method for the Navier–Stokes equations with singular forces. J Comput Phys 2001;171:822–42. [13] LeVeque R, Li Z. The immersed interface method for elliptic

equations with discontinuous coeﬃcients and singular sources. SIAM J Numer Anal 1994;31:1019–44.

[14] LeVeque R, Li Z. Immersed interface method for Stokes ﬂow with elastic boundaries or surface tension. SIAM J Sci Comput 1997;18:709–35.

[15] Le DV, Khoo BC, Peraire J. An immersed interface method for viscous incompressible ﬂows involving rigid and ﬂexible boundaries. J Comput Phys 2006;220:109–38.

[16] Linnick MN, Fasel HF. A high-order immersed interface method for simulating unsteady incompressible ﬂows on irregular domains. J Comput Phys 2005;204:157–92.

[17] Liu X-D, Fedkiw RP, Kang M. A boundary condition capturing method for Poisson’s equation on irregular domains. J Comput Phys 2000;160:151–78.

[18] Mayo A. The fast solution of Poisson’s and the biharmonic equations on irregular regions. SIAM J Numer Anal 1984;21:285–99. [19] Peskin CS. Numerical analysis of blood ﬂow in the heart. J Comput

Phys 1972;25:220–52.

[20] Peskin CS. The immersed boundary method. Acta Numer 2002:1–39. [21] Russell D, Wang ZJ. A Cartesian grid method for modeling multiple moving objects in 2D incompressible viscous ﬂow. J Comput Phys 2003;191:177–205.

[22] Xu J, Zhao H-K. An Eulerian formulation for solving partial diﬀerential equations along a moving interface. J Sci Comput 2003;19:573–94.

[23] Xu S, Wang ZJ. An immersed interface method for simulating the interaction of a ﬂuid with moving boundaries. J Comput Phys 2006;216:454–93.

[24] Wiegmann A, Bube KP. The explicit-jump immersed interface method: Finite diﬀerence method for PDEs with piecewise smooth solutions. SIAM J Numer Anal 2000;37:827–62.