內嵌介面問題的數值方法與應用

國

立

交

通

大

學

應用數學系

碩 士 論 文

內嵌介面問題的數值方法與應用

Numerical Methods and Applications

for Immersed Interface Problems

研 究 生：曾孝捷

指導教授：賴明治 教授

內嵌介面問題的數值方法與應用

Numerical Methods and Applications

for Immersed Interface Problems

研 究 生：曾孝捷 Student：Hsiao-Chieh Tseng

指導教授：賴明治 Advisor：Ming-Chih Lai

國 立 交 通 大 學

應 用 數 學 系

碩 士 論 文

Submitted to Department of Applied Mathematics College of Science

National Chiao Tung University in partial Fulfillment of the Requirements

for the Degree of Master

in

Applied Mathematics July 2006

Hsinchu, Taiwan, Republic of China

內嵌介面問題的數值方法與應用

學生：曾孝捷

指導教授：賴明治 教授

國立交通大學應用數學系﹙研究所﹚碩士班

摘

要

本論文介紹計算有不連續條件的二維 Poisson 方程的一個簡單快速

的數值方法，並應用本方法來解 Stokes 方程和 Navier-Stokes 方程。在

本文中，提供在計算區域上，處理介面的幾何形狀以及在介面上方程的

不連續條件。我們使用三次雲線內插來近似不連續條件所在之介面，以

及使用二階精度的外插公式來處理不連續條件。在效率的考量上，使用

快速的 Poisson 算法。我們給出一些計算數據，來觀察此方法的精度以

及應用。

Numerical Methods and Applications for

Immersed Interface Problems

Student: Hsiao-Chieh Tseng

Advisor: Dr. Ming-Chih Lai

Department (Institute) of Applied Mathematics

National Chiao Tung University

Abstract

In this paper, we introduce a simple and efficient second order

numerical method to solve Poisson equation with jump conditions in

two-dimension (2-D). We also apply the method to simulate Stokes

and unsteady Navier-Stokes flow dynamics. One of the essential

components of the method is to handle the geometric of the

immersed interface and jump discontinuities along it in the

computational domain. Cubic spline interpolation is implemented to

deal with the interface and a second order accurate extrapolation is

applied to incorporate these jump conditions. We also apply the

marker-and-cell (MAC) method cooperating with the interface. In

the matter of efficiency, we can apply some fast Poisson solvers such

as Fishpack. We verify the accuracy, efficiency of our method by

showing some numerical experiments.

誌

謝

本篇論文的完成，首先要感謝我的指導老師 賴明治教授。記得剛

上研究所，老師帶著懵懂的我，一步一步地踏入科學計算這個領域，指

引我學習與研究的方向，並不時給予我鼓勵。從中，不僅學到許多相關

的知識與研究的方法，更讓我在這片領域中，找到了興趣與成就感。除

了指導老師之外，實驗室的夥伴們也給予我不少助力。曾昱豪學長是我

常請教的對象，他豐富的經驗與資源，讓我在研究的過程中，找到明確

的方向不至於迷惘；吳恭儉學長以他深厚的數學分析能力，提供我許多

在理論方面的驗證與推導。還有李雅羚同學與陳冠羽學弟和我在課業與

研究上的切磋，以及許多以畢業學長姐留下來的資料。學生在此謹致最

誠摯的謝意。

在論文口試期間，承蒙牛仰堯教授、楊肅煜教授，張書銘教授費心

審閱並提供許多寶貴意見，使本論文之完稿得以更加齊備，學生永銘在

心。

在研究生活外，我也要感謝我的研究室好夥伴：陳偉國、涂芳婷、

楊川和。除了課業上的討論之外，在閒暇之餘和他們話家常，交換各種

心得與想法，使我增廣見聞及開拓生活圈，在彼此的相互鼓勵下，一同

度過這段歲月。

最後，要感謝我的家人給我的關心與祝福，以及女友蔡琬琦的呵護

及支持。這是我心靈的寄託與避風港。他們對我的付出，讓我可以順利

地，專心完成我的學業。願與他們，以及所有關心我的人，分享此篇論

文完成之喜悅與榮耀。

目

錄

中文提要

………

i

英文提要

………

ii

誌謝

………

iii

目錄

………

iv

一、

Introduction ………

1

二、

Poisson Equation with Jump Discontinuities ……

2

2.1

Introduction ………

2

2.2

Discontinuous Poisson Problem ………

2

2.3

Cubic Spline Representation ………

3

2.3.1 A Brief Introduction to Cubic Spline

Representation ………

4

2.3.2 Solving the Linear System Ax = b with Specific

Matrix A ………

5

2.3.3 Flag Construction ………

6

2.4

Interpolation from Grid to Interface ………

7

2.5

Numerical Example ………

7

三、

Stokes Flow with Singular Force ………

11

3.1

Stokes Equation ………

11

3.2

Marker and Cell Method ………

11

3.3

Numerical Scheme ………

12

3.4

Handling the Derivatives with a Crossing

Interface ………

13

3.5

Numerical Example

14

四、

Navier-Stokes Equation ………

18

4.1

Introduction to Navier-Stokes Equation …………

18

4.2

Numerical Scheme ………

21

4.3

Numerical Example ………

23

4.3.1 Accuracy Check without Interfaces ………

23

4.3.2 Cavity Flow ………

24

4.3.3 Circular Flow with a Line Force ………

26

五、

Concluding Remarks ………

28

1

Introduction

In problems with interfaces, the unknowns or their derivatives may have jump discontinu-ities. For example, fluid flow or wave propagation an interface between different regions exerts a force on the material, or an interface separates regions of different material prop-erties. The interface may be a fixed obstacle, or be considered as an elastic membrane. The interface is often embedded in a regular computational domain, and varieties of fast methods can apply with some modifications. On the other hand, it also need to be considered that the accuracy, especially near the interface.

There are lots of treatments for the immersed interface problems. Mayo in this paper [9] used the fundamental solutions of Poisson equations to give some integral formulas on the interface. These seems to be some relations between Mayo’s treatment and the boundary element method.

In the 1970-80’s, C. S. Peskin developed the Immersed Boundary Method (IBM) to simulate the blood flow through heart valves. The valve leaflet (the interface) exerts some force into the blood (the fluid) and t the same time moves along with the local fluid velocity simultaneously. The interaction between the fluid and immersed boundary can be modeled by a well chosen discretized approximation to the Dirac delta function, which is called a discrete delta function. However, it is also known that generally the IBM is only first-order accurate for nonsmooth but continuous quantities. Different variations of IBM have been proposed to improve the accuracy. For example, Lai and Peskin [10] proposed a new formally second-order IBM with reduced numerical viscosity, and Griffith and Peskin [11] provided another higher order method for sufficiently smooth problems.

In contrast to the discrete delta function approach, Leveque and Li proposed Immersed Interface Method (IIM). It translates the jump discontinuities into the local coordinate system of the interface, and then find the correction terms to modify the system in Carte-sian grid. They apply this method to solve some fluid problems and obtain second-order accurate [8]. On the other hand, Jane Wang and her group computed the correction terms by extrapolation the jump discontinuities from the interface to the Cartesian grid [12]. They also made some applications in fluid mechanics [13].

The Ghost Fluid Method (GFM), which is also a sharp interface method that uses the jump conditions in the solution and the flux, has been developed by Fedkiw, Liu, Kang and their group for elliptic interface problems. The GFM is simpler than IIM and the resulting linear system of equations for the self-adjoint elliptic equations is symmetric positive definite. The trade-off is that the solution is first-order accurate in the maximum norm.

In this paper, we implement the extrapolation method in [12] and apply it to Poisson, Stokes, and full Navier-Stokes equation.

2

Poisson Equation with Jump Discontinuities

2.1

Introduction

Given a scalar function φ with real spatial variable x, we denote [φ(x∗)] = φ(x+∗) − φ(x

− ∗),

the right limit subtract the left limit at x∗. We call φ has a jump at x∗ if [φ(x∗)] 6= 0.

Our numerical method is based on finite difference method (FDM) with Cartesian Grid points. We denote the irregular points whose difference stencil crosses the jump, and the others are called regular points.

Taking an 1-D example. Consider the case that one wants to obtain some approxima-tion of derivatives of φ at xj with grid size h but there is a jump of φ at x∗ between the

irregular points xj and xj+1. Let xj+1 = x∗+ α, Taylor expansion tells us that

φ(xj+1) = φ(x+∗ + α) = φ(x+_∗) + φ0(x+_∗)α + φ00(x+_∗)α 2 2 + O(h 3₎ = φ(x−_∗) + φ0(x−_∗)α + φ00(x−_∗)α 2 2 + C (φ) j+1+ O(h 3₎ where C_j+1(φ) = [φ] + [φ0]α + [φ00]α 2 2

is the correction term of φ(xj+1) for xj. Hence, for example we can obtain

φ0(xj) = φ(xj+1) − C (φ) j+1
− φ(xj−1) 2h + O(h 2_), φ00(xj) = φ(xj+1) − C (φ) j+1
− 2φ(xj) + φ(xj−1) h2 + O(h).

2.2

Discontinuous Poisson Problem

We consider the interface problem

∆ψ± = ω± in Ω±, (2.1) [ψ] = q0 on Γ,  ∂ψ ∂n  = q1 on Γ, ψ = ψb on ∂Ω,

where ω±, ψ±are defined on Ω±and q0, q1are on Γ; with Γ ⊆ Ω, Ω− = Ω−∪Γ, Ω+= Ω+∪Γ.

We assume that ω±, q0, q1 and Γ are fairly smooth, with a possible jump in ω± at Γ.

We generate the 2-D Cartesian grid with grid size h = ∆x = ∆y, and overlay an irregular boundary Γ with known discontinuities in [ψ], ∂ψ_∂n. Here we assume that Γ can be represented by a periodic cubic spline. In our approach we find the foot of perpendicular on Γ for every irregular point. In some sense this is the most suitable location to be chosen

since we know the discontinuities in terms of the normal derivatives. Instead of finding h

∂2_ψ

∂n2

i

, we apply the identity of Laplace operator on a curve, i.e., ω = ∂ 2_ψ ∂n2 + κ ∂ψ ∂n + ∂2ψ ∂s2,

and the final correction formula is a simple Taylor series expansion C₄(ψ) = [ψ•] + α  ∂ψ• ∂n  +1 2α 2 ∂2ψ• ∂n2  = [ψ•] + α  ∂ψ• ∂n  +1 2α 2  [ω•] − κ  ∂ψ• ∂n  − ∂ 2_[ψ •] ∂s2 

where α is the signed distance to the foot of perpendicular and κ is the local curvature at that point. More specifically, α is positive if the grid point is outside Γ (in Ω+), and

vice versa. The final correction to ω is then

ω_♦ update= ω_♦− C

(ψ) 4

h2 .

After correcting the right hand side ω, we apply a fast Poisson solver called Fishpack to solve the Poisson equation ∆ψ = ω and then obtain ψ immediately.

Γ ♦ × × M M ZZ b b α s s Ω+ Ω−

Figure 1: A 5-point Laplacian at ♦. The two M are irregular points related to the irregular point ♦.

2.3

Cubic Spline Representation

We interpret the interface Γ by picking some points on Γ, called Lagrangian markers, and using these markers to generate a parameterized periodic cubic spline. We denote the cubic spline by X(s) = (X(s), Y (s)) where s, 0 < s < L, is the arc length parameter. Hence, |∂X/∂s| = 1 and then the local curvature can be expressed by κ = X0Y00− Y0_X00_.

The foot of perpendicular X(si,j) of the grid point xi,j = (xi, yj) can be found by bisection

method just like root finding since the inner product of tangential and normal vector on a curve is equal to zero, That is,

xi− X(si,j)
X0(si,j) + yj − Y (si,j)
Y0(si,j) = 0.

Thus, we have the signed distance α = ±

q

2.3.1 A Brief Introduction to Cubic Spline Representation

The cubic spline interpolation is significant since it yields C2 _{approximation of curves and}

it is sufficiently smooth in the presence of small curvatures. In the construction of the periodic parametric cubic spline X(s), we give the markers

Xk = (Xk, Yk), k = 1, 2, . . . , M

in order. We also define that X(sk) = Xk, and impose that XM +1= X1 for convenience.

For k = 1, 2, . . . , M , We define the length of each segment by hk+1 = kXk+1− Xkk2

and the accumulated length by

H1 = 0, Hk+1= Hk+ hk+1

at each marker, and then obtain the spline representation

X(s) = Xk+ a1(s − Hk) + a2(s − Hk)2+ a3(s − Hk)3, s ∈ [Hk, Hk+1] where a1 = (Xk+1− Xk) hk+1 − (2Mk+ Mk+1) hk+1 6 , a2 = Mk 2 , a3 = (Mk+1− Mk) 6hk+1 .

That is, the spline can be characterized by its moments Mi = X00(si). In fact, we

implement the parameterized cubic spline by using the cumulative chord length since the chord length is approximately equal to the arc length if the markers are close enough. For more detail, see [3] and [4].

By the continuity of X0(s) and periodicity of X(s), we have the moment equation µjMj−1+ 2Mj + λjMj+1= dj M0 = MM, MM +1= M1 for j = 1, 2, . . . , M , where λj = hj+1 hj+ hj+1 µj = hj hj+ hj+1 = 1 − λj dj = 6 hj+ hj+1  Xj+1− Xj hj+1 −Xj − Xj−1 hj 

The task of calculating these moments is to solve a nonsingular linear system with a tridiagonal matrix containing two nonzero elements at the right-top and left-bottom corner.

Hence, we need to solve the nonsingular linear system           2 λ1 µ1 µ2 2 λ2 µ3 . .. ... . .. ... ... . .. 2 λM −1 λM µM 2                     m1 m2 .. . .. . .. . mM           =           d1 d2 .. . .. . .. . dM           for each component of X(sj).

2.3.2 Solving the Linear System Ax = b with Specific Matrix A

We introduce Thomas algorithm, a simplified form of Gaussian elimination that can be used to solve tridiagonal system of equations

cjxj−1+ djxj+ ejxj+1 = bj, j = 1, 2, . . . , n,

c1 = 0, en = 0.

in O(n). In this case

A =           d1 e1 c2 d2 e2 c3 . .. ... . .. ... ... . .. dn−1 en−1 cn dn           .

Thomas algorithm is given here: !! Forward elimination phase do k = 2,n

m = c(k) / d(k-1) d(k) = d(k) - m*e(k-1) b(k) = b(k) - m*b(k-1) end do

!! Backward substitution phase x(n) = b(n)/d(n)

do k = n-1,1,-1

x(k) = (b(k) - e(k)*x(k+1)) / d(k) end do

For periodic case, i.e., system of equations gives that

cjxj−1+ djxj+ ejxj+1 = bj, j = 1, 2, . . . , n,

and A =           d1 e1 c1 c2 d2 e2 c3 . .. ... . .. ... ... . .. dn−1 en−1 en cn dn           .

We can solve the system by using Thomas algorithm twice. 1. Solve cjyj−1+ djyj + ejyj+1 = bj, j = 1, 2, . . . , n − 1 with y0 = yn= 0 2. Solve cjzj−1+ djzj+ ejzj+1 = 0, j = 1, 2, . . . , n − 1 with z0 = zn= 1 3. Obtaining xn= bn− eny1− cnyn−1 dn− cnzn−1− enz1 , and then we have

xj = yj+ xnzj, j = 1, 2, . . . , n − 1.

A subroutine DGLTSL in Linpack can solve this periodic type of linear system in O(N ), where N is the problem size.

2.3.3 Flag Construction

As we construct the spline, we also need to keep track of the irregular points. We allocate an integer array, called the flag, which denotes the Cartesian grid types. Since we assume that the interface is closed, the flag should record these types such as regular points inside/outside the interface and irregular points inside/outside/on the interface. By the way, the properties of the flag is somewhat like those of a level set function. We find the irregular point first by finding the neighbor point of the intersection of spline and grid line. Then, using some filling algorithm, the domain can be distinguished to be inside or outside. A basic pseudo code of filling algorithm is given here:

recursive subroutine FloodFill(ix, iy, oldColor, newColor)

implicit none; integer, intent(in) :: ix, iy, oldColor, newColor if(IsOutside(ix,iy)) return

if(iflag(ix,iy) == oldColor) then iflag(ix,iy) = newColor

call FloodFill(ix+1, iy, oldColor, newColor) call FloodFill(ix-1, iy, oldColor, newColor) call FloodFill(ix, iy+1, oldColor, newColor) call FloodFill(ix, iy-1, oldColor, newColor) endif

return

Corresponding to the spline representation, the flag also plays an important role in the further computation.

2.4

Interpolation from Grid to Interface

Let ψ∗ = ψ(x∗, y∗) be the value on the interface which we want to interpolate. First we

need to find the index (i, j) of grid point which is located on the left-bottom corner of the cell. Then, let ph = x∗− xi, qh = y∗− yj, with grid ratio 0 ≤ p, q ≤ 1. Hence

ψ+ = pψi+1,j+1+ (1 − p)ψi,j+1+ O(h2)

ψ× = pui+1,j + (1 − p)ψi,j+ O(h2)

ψ∗ = qψ++ (1 − q)ψ×+ O(h2)

= pqψi+1,j+1+ (1 − p)qψi,j+1+ p(1 − q)ψi+1,j+ (1 − p)(1 − q)ψi,j + O(h2).

The interpolation of an interface point from neighbor irregular points should be modified by correction terms. We define

χi,j = ( 1 xi,j ∈ Ω+ 0 xi,j ∈ Ω− and we have ψ∗ update = ψ∗+pqχi+1,j+1C (ψ) i+1,j+1+(1−p)qχi,j+1C (ψ) i,j+1+p(1−q)χi+1,jC (ψ) i+1,j+(1−p)(1−q)χi,jC (ψ) i,j . t t t t (i, j) (i + 1, j) (i, j + 1) (i + 1, j + 1) ∗ × + d d d d d ph (1 − p)h qh (1 − q)h

Figure 2: Interpolation on point ∗ using the four corner points.

2.5

Numerical Example

We consider three modified Poisson equations with exact solutions on Ω = [−1, 1]×[−1, 1] with jump discontinuities defined on a ellipse

Γ : x

2

a2 +

y2 b2 = 1.

Example. 1 Example. 2 Example. 3 ψ− 1 x2− y2 sin(x) cos(y)

ψ+ 1 + ln(2px2+ y2) 0 0

ω− 0 0 −2 sin(x) cos(y)

ω+ 0 0 0

Table 1: Three test cases for Eq. (2.1).

N k · k∞ ratio k · kL2 ratio

Example. 1 16 3.0347E-02 0.0000 9.1019E-03 0.0000 32 5.4561E-03 5.5621 6.4576E-04 14.0949 64 1.7111E-03 3.1887 2.2393E-04 2.8838 128 5.3475E-04 3.1997 1.6647E-04 1.3452 256 8.7802E-05 6.0904 1.6368E-05 10.1703 Example. 2 16 7.2378E-02 0.0000 1.0009E-02 0.0000 32 1.7418E-02 4.1553 1.3929E-03 7.1855 64 2.7228E-03 6.3972 4.1537E-04 3.3535 128 1.0008E-03 2.7206 1.3772E-04 3.0160 256 1.0402E-04 9.6208 1.6698E-05 8.2479 Example. 3 16 4.2294E-02 0.0000 6.2970E-03 0.0000 32 7.7546E-03 5.4540 7.7777E-04 8.0963 64 1.2982E-03 5.9732 1.0272E-04 7.5721 128 5.6115E-04 2.3135 1.1584E-04 0.8867 256 1.4999E-04 3.7413 2.3262E-05 4.9799 Table 2: Test results for a = 0.6, b = 0.4.

In the first case, the convergent rate fluctuates in a wide range. The accuracy may be affected by how the spline approximates the actual ellipse interface. General speaking, it has second order accuracy on average.

N k · k∞ ratio k · kL2 ratio

Example. 1 16 2.1539E+00 0.0000 5.7049E-01 0.0000 32 5.4730E-01 3.9355 9.5492E-02 5.9743 64 1.2800E-02 42.7584 6.9241E-04 137.9125 128 2.5374E-03 5.0446 1.9674E-04 3.5193 256 9.1706E-04 2.7668 4.2801E-05 4.5967 Example. 2 16 3.9062E-01 0.0000 4.0670E-02 0.0000 32 2.3950E-01 1.6310 4.0821E-02 0.9963 64 2.1336E-02 11.2252 1.0262E-03 39.7800 128 4.1392E-03 5.1545 3.9887E-04 2.5727 256 1.4158E-03 2.9237 5.4272E-05 7.3495 Example. 2 16 5.6932E-01 0.0000 7.1246E-02 0.0000 32 1.2497E-01 4.5556 1.3391E-02 5.3204 64 8.8370E-03 14.1420 5.5884E-04 23.9622 128 1.5545E-03 5.6848 9.7281E-05 5.7446 256 5.4722E-04 2.8408 3.2040E-05 3.0363 Table 3: Test results for a = 0.8, b = 0.2.

In the second case, the curvature varies severely. In our test, the cartesian grid cannot fetch the interface behavior as N < 64 in this table. As N exceeds this value, we also observe second order accuracy on average.

3

Stokes Flow with Singular Force

3.1

Stokes Equation

The steady Stokes equation is of this form

−∇p + µ∆u = f , ∇ · u = 0.

where u = u(x), v(x)
is the fluid velocity, p = p(x) is the fluid pressure, and the µ is the constant fluid viscosity. The force density f (x) = (f1(x), f2(x)) is applied by the

immersed interface Γ into the fluid, which the boundary force is F (s) = (F1(s), F2(s)).

More clearly, the singular force f has support only on Γ and has the form f (x) =

Z

Γ

F (s)δ(x − X(s)) ds

where δ is the Dirac delta function defined in 2-D. The boundary configuration of Γ is given by X(s) = (X(s), Y (s)), 0 ≤ s ≤ L, where the are-length parameter s tracks every material point of the initial configuration.

We write a 2-D vector form of this equation, which yields −px py  + µ∆u ∆v  =f1 f2  , ux+ vy = 0.

One can expect that the velocity and the pressure will not be smooth across the boundary. In fact the boundary force F plays an important role to the jump condition of u, p, and their derivatives. In the computational fluid dynamics, we usually impose the no-slip condition of fluid velocity, i.e., [u] = 0 along the interface. Thus, we have the following theorem.

Theorem. Let u, p be solution of Stokes’ equation. Then the jump across the interface satisfies [p] = F · n |∂X/∂s| , µ  ∂u ∂n  = −F · τ |∂X/∂s|τ ,  ∂p ∂n  =  1 ∂X_∂s   ∂ ∂s  F · τ |∂X/∂s|  , where n, τ are unit normal and tangent vector respectively on the interface.

We will give a proof in the next section.

3.2

Marker and Cell Method

We use marker and cell (MAC) method and take the computation domain Ω to be square (or rectangular) with staggered grid in Cartesian coordinates. For simplicity, we take the mesh size h = ∆x = ∆y to be equal and let Nx, Ny be the number of mesh in the

x-and y- directions respectively. Usually, we pick the same number of grid points in each direction, and denote that N = Nx = Ny.

r r r r r r r r r r r r r r r r r r r r r r r r r r r r r r 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × × u0,1 u1,1 uNx−1,1 uNx,1 u0,Ny uNx,Ny v1,0 v1,1 v1,Ny−1 v1,Ny vNx,0 vNx,1 vNx,Ny−1 vNx,Ny p1,1 pNx,1 p1,Ny pNx,Ny

Figure 4: The computational domain Ω using staggered Grid with mesh size h.

3.3

Numerical Scheme

In the following numerical scheme, we assume s be an approximation of arc-length pa-rameter such that |∂X/∂s| is nearly equal to 1. For simplicity, we denote Fn = F · n

and Fτ = F · τ . By translating the boundary force F to the jump condition, we have

−∇p + µ∆u = 0, ∇ · u = 0, [u] = 0, [p] = Fn, µ  ∂u ∂n  = −Fττ ,  ∂p ∂n  = ∂Fτ ∂s on Γ, u    ∂Ω= ub.

Consequently, by taking divergent to the above equations we have ∆p = 0, [p] = Fn,  ∂p ∂n  = ∂Fτ ∂s on Γ.

We need to impose a boundary condition for the above Laplace system. For the test example, we give the Dirchlet boundary condition with given solution p calculated on the boundary. In general, we consider the periodic boundary condition. Although any two solutions of this system may be up to a constant, it is fine since we need only ∇p.

Here we summarize the solution procedures. 1. Solve the Laplace equation

and obtain p, where the correction tern C(p) = 1 h2  [p] + α ∂p ∂n  +1 2α 2  [∆p] − κ ∂p ∂n  − ∂ 2_[p] ∂s2  = 1 h2  Fn+ α ∂Fτ ∂s + 1 2α 2  −κ∂Fτ ∂s − ∂2_F n ∂s2  .

2. Compute px, py on the u- and v-grid by center difference (pi+1,j − pi,j)/h, (pi,j+1−

pi,j)/h respectively. It needs some correction when the p-grid points straddle across

the interface. Let us explain the procedure more clearly by taking an example. Note that the location of ui,j is the middle point between those of pi,j and pi,j. Let

pi+1,j ∈ Ω− and pi,j, ui,j ∈ Ω+. The key point is to correct the point on the different

side with the middle point. We need to add the correction term if the point is in Ω−, and vice versa. In this case we need to add the correction term C_i+1,j(p) of pi+1,j,

which is defined as above. That is, the approximation of px can be expressed by

((pi+1,j + C (p)

i+1,j) − pi,j)/h. We will discuss the detail in the next section.

3. Compute [∇p] which appears in the correction term C(u) _{later. Recall that τ =}

(τ1, τ2), the unit tangent vector, is equal to (−n2, n1) at each point on Γ, where

n = (n1, n2) is the unit normal vector. Since

∂Fτ ∂s =  ∂p ∂n  = [px] n1+ [py] n2, ∂Fn ∂s = ∂ [p] ∂s = [px] τ1+ [py] τ2, we have [∇p] = ([px] , [py]) =  τ2 ∂Fτ ∂s + τ1 ∂Fn ∂s , −τ1 ∂Fτ ∂s + τ2 ∂Fn ∂s  . 4. Solve the Poisson equation

µ∆hu + C(u) = ∇hp where C(u)= 1 h2  µ [u] + αµ ∂u ∂n  +1 2α 2  µ [∆u] − κµ ∂u ∂n  − µ∂ 2_[u] ∂s2  = 1 h2  −αFττ + 1 2α 2_{([∇p] + κF} ττ )  .

3.4

Handling the Derivatives with a Crossing Interface

Consider the derivatives approximated by center difference scheme with a crossing inter-face. In general, there are four cases about locations of grid points and the interinter-face. We use a table to illustrate how to calculate the first order derivative at a certain point.

Case Figure Difference formula (I) r φL r φR b φM Ω+Ω− Γ dφM dx ' (φR+ CR) − φL h (II) r φL r φR b φM Ω+Ω− Γ dφM dx ' φR− (φL− CL) h (III) r φL r φR b φM Ω−Ω+ Γ dφM dx ' (φR− CR) − φL h (IV) r φL r φR b φM Ω−Ω+ Γ dφM dx ' φR− (φL+ CL) h

Table 4: The four cases about locations of grid points L, M, R and the interface. h denotes the length of the interval [L, R], and M is the middle point in [L, R] here. CL, CRare the

correction terms at grid points L, R respectively.

3.5

Numerical Example

Example (1). Normal forces on a circle. Consider a circular boundary of radius R = 1 parametrized by x(θ) = (cos(θ), sin(θ)) in Ω = [−2, 2] × [−2, 2] and define the force on the boundary by

f (θ) = 2 sin(3θ)x(θ)

so that the forces are normal to the boundary. The exact solution is given by p(r, θ) = ( −r3_{sin(3θ), r < 1} r−3sin(3θ), r > 1 u(r, θ) = ( 3 8r 2_{sin(2θ) +} 1 16r 4_{sin(4θ) −}1 4r 4_{sin(2θ), r < 1} 1 8r −2_{sin(2θ) −} 3 16r −4_{sin(4θ) +} 1 4r −2_{sin(4θ), r ≥ 1} v(r, θ) = (₃ 8r 2_{cos(2θ) −} 1 16r 4_{cos(4θ) −}1 4r 4_{cos(2θ), r < 1} 1 8r −2_{cos(2θ) +} 3 16r −4_{cos(4θ) −} 1 4r −2_{cos(4θ), r ≥ 1}

We note that on the interface R = 1, [p] = 2 sin(3θ),  ∂p ∂n  = 0, [∆p] = 0, ∂ 2_[p] ∂s2 = −18 sin(3θ) and [u] = 0,  ∂u ∂n 

= 0, [∆u] = −6 cos(3θ) sin(θ) 6 cos(3θ) cos(θ)

 , ∂

2_[u]

N kuerrk∞ Ratio kverrk∞ Ratio kperrk∞ Ratio

16 6.31625996E-02 - 5.71511520E-02 - 1.80064548E-01 -32 1.15154249E-02 5.4850 1.04312328E-02 5.4788 6.70713180E-02 2.6846 64 2.46761660E-03 4.6666 3.44416038E-03 3.0286 2.41626231E-02 2.7758 128 2.31427096E-04 10.6626 6.66326129E-04 5.1688 3.79474345E-03 6.3673

Table 5: A grid refinement analysis of u, v, p for Example(1).

Example (2). Tangential forces on a circle. We consider a similar case. The domain Ω and the circular boundary are defined as Example(1). Define the force on the boundary by

f (θ) = 2 sin(3θ)x0(θ)

so that the forces are tangential to the boundary. The exact solution is given by p(r, θ) = ( −r3_{cos(3θ), r < 1} −r−3_{cos(3θ), r > 1} u(r, θ) = ( 1 8r 2_{cos(2θ) +} 1 16r 4_{cos(4θ) −}1 4r 4_{cos(2θ), r < 1} −1 8r −2_{cos(2θ) +} 5 16r −4_{cos(4θ) −} 1 4r −2_{cos(4θ), r ≥ 1} v(r, θ) = ( −1 8r 2_{sin(2θ) +} 1 16r 4_{sin(4θ) +} 1 4r 4_{sin(2θ), r < 1} 1 8r −2_{sin(2θ) +} 5 16r −4_{sin(4θ) −} 1 4r −2_{sin(4θ), r ≥ 1}

We note that on the interface R = 1, [p] = 0,  ∂p ∂n  = 6 cos(3θ), [∆p] = 0, ∂ 2_[p] ∂s2 = 0 and [u] = 0,  ∂u ∂n  =  2 sin(3θ) sin(θ) −2 sin(3θ) cos(θ) 

, [∆u] = 3 cos(2θ) + 3 cos(4θ) −3 sin(2θ) + 3 sin(4θ)

 , ∂

2_[u]

∂s2 = 0.

N kuerrk∞ Ratio kverrk∞ Ratio kperrk∞ Ratio

16 5.56812172E-02 - 1.84446032E-02 - 1.92715009E-01 -32 9.06017465E-03 6.1457 5.85585435E-03 3.1497 2.08500578E-02 9.2429 64 3.34797871E-03 2.7061 1.43026464E-03 4.0942 5.88284456E-03 3.5442 128 6.04375427E-04 5.5395 3.19944420E-04 4.4703 1.32162086E-03 4.4512

Table 6: A grid refinement analysis of u, v, p for Example(2).

Example (3). Combine the normal and tangential forces on a circle. Since the Stokes equation is linear, we combine the above two example by taking summation.

N kuerrk∞ Ratio kverrk∞ Ratio kperrk∞ Ratio

16 1.17188492E-01 - 6.20299174E-02 - 2.02998893E-01 -32 1.22218170E-02 9.5884 1.25574945E-02 4.9396 8.55532629E-02 2.3727 64 3.73073530E-03 3.2759 3.67928831E-03 3.4130 2.61792711E-02 3.2679 128 5.23096557E-04 7.1320 7.33869339E-04 5.0135 4.55700360E-03 5.7448

Table 7: A grid refinement analysis of u, v, p for Example(3).

Example (4). A moving interface problem. This example is given in [8]. We consider a moving interface problem which has uniform density and viscosity. The initial interface is given by

r(θ) = r0(1 +  sin(kθ)), 0 ≤ θ ≤ 2π.

The initial velocity and pressure are all set to be zero. The only driven force of the fluid motion is the surface tension which is proportional to the curvature κ with a constant E, that is,

Fn = −Eκ, Fτ = 0.

The velocity is smooth but the pressure has a non-constant jump [p] = Fn,

 ∂p ∂n 

= 0. We simply take an interface moving formula which says

Xn+1= Xn+ Un+1∆t.

Figure 5: The numerical evolution of the moving interface using the IIM coupled with the interpolation formula introduced in previous section. Here we pick the number of grid: Nx = Ny = N = 64, number of interface marker: M = 42, and ∆t = 0.002.

Figure 6: The measure of the bubble area. The geometry of interface eventually converges to a circle with conserving measure of area.

4

Navier-Stokes Equation

4.1

Introduction to Navier-Stokes Equation

The Navier Stokes Equation is of this form ρ ∂u

∂t + (u · ∇)u 

+ ∇p = µ∆u + f + g, (4.2) ∇ · u = 0, (4.3) where u(x, t), p(x, t), µ, f (x, t) are similar as the above session, and g(x, t) is the external force defined on the computational domain. Note that the divergent-free condition ∇·u = 0 denotes the fluid is incompressible.

Here we have some important theorems for Navier-Stokes equation, which can also be applied to time-independent Stokes equation.

Lemma. If [u] = 0, then

[ut] + [∇u] · u = 0. Proof. 0 = d dt[u] =  ut+ ∇u · ∂X ∂t  = [ut+ ∇u · u] = [ut] + [∇u] · u.

Lemma. If [u] = 0 and ∇ · u = 0, then  ∂u

∂n 

· n = 0.

Proof. Denoting the unit normal and tangential vector be n = (n1, n2), τ = (−n2, n1),

and imposing ux+ vy = 0, we have

 ∂u ∂n  · n = [uxn1+ uyn2] n1+ [vxn1+ vyn2] n2 = [uxn1+ vxn2] n1+ [uyn1+ vyn2] n2 = [vy(−n1) + vxn2] n1+ [uyn1 + ux(−n2)] n2 = − [vyn1− vxn2] n1+ [uyn1− uxn2] n2 = − ∂u ∂τ  · τ = 0 since  ∂u ∂τ  = ∂ [u] ∂τ = 0.

Theorem. Let u, p be solution of Navier-Stokes equation without the external force g. Then the jump across the interface satisfies

[p] = F · n |∂X/∂s| , µ  ∂u ∂n  = −F · τ |∂X/∂s|τ ,  ∂p ∂n  =  1 ∂X_∂s   ∂ ∂s  F · τ |∂X/∂s| 

Proof. First we choose a banded domain Ω(t) enclosing the interface Γ(t) with the outer

and inner boundary Γ+

 (t) and Γ −

 (t). The  is a small parameter of the distance from Γ(t)

to Γ+

 (t) and to Γ −

(t), where Γ(t) = {X(s, t), 0 ≤ s ≤ L} be the parametric curve with

parameter s and dependence on time t. Let φ(x) be any smooth function with compact support defined on Ω(t). We multiply φ(x) on both sides of Navier-Stokes equation and

integrate over Ω(t), which yields

Z Ω(t) ρ ∂u ∂t + u · ∇u  φ dx + Z Ω(t) ∇pφ dx = Z Ω(t) µ∆uφ dx + Z Ω(t) Z L 0 F (s, t)δ(x − X(s, t))ds φ dx.

Using the relation of the transport theorem in the first term of the left hand side, and the fact of Ω(t) enclosing Γ(t) and the definition of Dirac delta function in the last term of

right hand side, we have d dt Z Ω(t) ρuφ dx + Z Ω(t) ∇pφ dx = Z Ω(t) µ∆uφ dx + Z L 0 F (s, t)φ(X(s, t)) ds. Let us consider the first component of the above equation, that is,

d dt Z Ω(t) ρuφ dx + Z Ω(t) pxφ dx = Z Ω(t) µ∆uφ dx + Z L 0 F1(s, t)φ(X(s, t)) ds.

Assume u is smooth enough and bounded, the first term of left hand side tends to zero as  → 0. Moreover Z Ω(t) pxφ dx = Z Γ+ pn1φ dσ − Z Γ− pn1φ dσ − Z Ω(t) pφx dx → Z Γ [p] n1φ dσ = Z L 0 [p] n1φ |∂X/∂s| ds and Z Ω(t) µ∆uφ dx = Z Γ+ µ∂u ∂nφ dσ − Z Γ− µ∂u ∂nφ dσ − Z Ω(t) µ∇u · ∇φ dx → Z Γ µ ∂u ∂n  φ dσ = Z L 0 µ ∂u ∂n  φ |∂X/∂s| ds as  → 0. Therefore, since φ is arbitrary, we have

[p] |∂X/∂s| n1 = µ

 ∂u ∂n 

|∂X/∂s| + F1.

Similarly, for the second component, we have [p] |∂X/∂s| n2 = µ

 ∂u ∂n 

|∂X/∂s| + F2.

Multiplying n1, n2 on both sides of the above equations respectively, summing them, and

by previous lemma, we obtain the first equality [p] = F · n

Now substituting this pressure jump into the above relations, we get µ ∂u ∂n  |∂X/∂s| = (F · n)n1− F1 = (−F1n2+ F2n1)n2 = (−F · τ )(−n2), µ ∂v ∂n  |∂X/∂s| = (F · n)n2− F2 = (F1n2− F2n1)n1 = (−F · τ )n1, that is, µ ∂u ∂n  = −F · τ |∂X/∂s|τ .

In order to derive the normal derivative jump of the pressure, let us apply the di-vergence operator to both side of Navier-Stokes equation and use the incompressibility constraint, we have ∇ ·  ρ ∂u ∂t + (u · ∇)u  + ∆p = ∇ · f .

As before, we multiply φ on both sides and integrate over Ω, which yields

Z Ω(t) ∇ ·  ρ ∂u ∂t + (u · ∇)u  φ dx + Z Ω(t) ∆pφ dx = Z Ω(t) ∇ · f φ dx. Observe that Z Ω(t) ∇ ·  ρ ∂u ∂t + (u · ∇)u  φ dx = Z Ω(t) ∇ ·  ρ ∂u ∂t + (u · ∇)u  φ  dx − Z Ω(t) ρ ∂u ∂t + (u · ∇)u  · ∇φ dx = Z Γ+∪Γ− ρ ∂u ∂t + (u · ∇)u  φ · n dσ − d dt Z Ω(t) ρu · ∇φ dx → 0 (as  → 0, by previous lemma, and u is bounded.) and, using twice the Green’s identity,

Z Ω(t) ∆pφ dx = Z Ω(t) ∆φp dx + Z Γ+∪Γ− φ∂p ∂n dσ − Z Γ+∪Γ− p∂φ ∂n dσ → Z Γ  ∂p ∂n  φ dσ − Z Γ [p]∂φ

∂n dσ (as  → 0, and p is bounded.) = Z L 0  ∂p ∂n  |∂X/∂s| φ ds − Z L 0 [p] |∂X/∂s| ∂φ ∂n ds.

delta function), we have Z Ω(t) ∇ · f φ dx = Z Ω(t) Z L 0  ∇ · (F δ(x − X(s, t))φ dsdx = Z L 0 F (s, t) Z Ω(t) (∇ · δ(x − X(s, t)))φ dx  ds = − Z L 0 F · ∇φ ds = − Z L 0 F · n∂φ ∂n ds − Z L 0 F · τ∂φ ∂τ ds = − Z L 0 F · n∂φ ∂n ds − Z L 0 F · τ |∂X/∂s| ∂φ ∂s ds = − Z L 0 [p] |∂X/∂s| ∂φ ∂n ds + Z L 0 ∂ ∂s  F · τ |∂X/∂s|  φ ds. Since φ is arbitrary, we obtain the normal derivative jump for the pressure

 ∂p ∂n  =  1 ∂X_∂s   ∂ ∂s  F · τ |∂X/∂s|  . Thus, the proof is complete.

Remark. The results can be generalized to the case with additional bounded piecewise continuous force g across the interface. The proof is straightforward without much effort needed, and the result of [u] ,∂u_∂n are still valid. However, the normal derivative jump condition for the pressure should be

 ∂p ∂n  =  1 ∂X_∂s   ∂ ∂s  F · τ |∂X/∂s|  + [g] · n.

We consider the fluid density ρ ≡ 1, and translate the interface force to the jump conditions. We have ∂u ∂t + (u · ∇)u + ∇p = µ∆u + g, ∇ · u = 0, [u] = 0 , µ ∂u ∂n  = −F · τ |∂X/∂s|τ , [p] = F · n |∂X/∂s| ,  ∂p ∂n  =  1 ∂X_∂s   ∂ ∂s  F · τ |∂X/∂s|  .

4.2

Numerical Scheme

There are 3 steps to approximate the solution numerically. Here we consider the scheme without interfaces first.

1 In the prediction step, we solve u∗ explicitly where u∗− un

∆t + (u · ∇)u

n+1_{2 + ∇p}n−1_{2 =} µ

2∆u

and ˜un+1 by a fast Poisson solver with jump discontinuity, ˜ un+1− u∗ ∆t = µ 2∆ ˜u n+1 , (4.5) ˜ un+1    ∂Ω= u n+1  ∂Ω,  ˜un+1 = 0, µ ∂ ˜u n+1 ∂n  = −F ∗_{· τ} |∂X/∂s|τ .

Note that the advection term (u · ∇)un+12 is approximated by the extrapolation

3

2(u · ∇)u

n₋1

2(u · ∇)u

n−1_.

We explain this step in detail. In the first component of the vector field we have these terms of the space discretization

 u∂u ∂x + v ∂u ∂y n+1₂ = 3 2  un_i,ju n i+1,j − uni−1,j 2∆x + v n stag un i,j+1− uni,j−1 2∆y  − 1 2 u n−1 i,j un−1_i+1,j− un−1 i−1,j 2∆x + v n−1 stag un−1_i,j+1− un−1 i,j−1 2∆y ! , ∂pn−1₂ ∂x = pn− 1 2 i+1,j − p n−1₂ i,j ∆x , ∆un = u n

i+1,j − 2uni,j+ uni−1,j

∆x2 +

un

i,j+1− 2uni,j+ uni,j−1

∆y2 ,

and the interpolation

vstag =

1

4(vi,j+ vi,j−1+ vi+1,j + vi+1,j−1) from v-grid to u-grid. The second component is similar.

2 We perform the projection step to correction the velocity to be divergent free and update the pressure. According to Helmholtz-Hodge decomposition, we have

˜

un+1 = un+1+ ∆t ∇φn+1, (4.6) ∇ · un+1 _{= 0.}

That is, we want to project the vector field ˜un+1 _{to a divergent free one, u}n+1_.

Taking divergent on both sides, we obtain ∇ · ˜un+1 ∆t = ∆φ n+1_, ∂φn+1 ∂n    ∂Ω = 0.

Note that we should calculate the divergent term ∇ · u∗ on the p-grid, solve the Laplace equation of φ, and then obtain un+1_.

3 Finally, the pressure should be corrected. See [6] which says that ∇pn+12 = ∇pn− 1 2 + ∇φn+1− µ∆t 2 ∇(∆φ n+1 ). Similarly, pn+12 = pn− 1 2 + φn+1−µ 2∇ · ˜u n+1_.

4.3

Numerical Example

4.3.1 Accuracy Check without Interfaces

Here we introduce some example with exact solutions. In the following examples, we put the time marching step ∆t = h/2 where h is again the mesh size of the staggered grid. Example (1). Consider the Navier-Stokes equation on a domain Ω = [0, 1] × [0, 1], and give the exact solution

u(x, y, t) = − cos(πx) sin(πy) exp(−2π2t/ Re), v(x, y, t) = sin(πx) cos(πy) exp(−2π2t/ Re),

p(x, y, t) = −0.25 (cos(2πx) + cos(2πy)) exp(−4π2t/ Re). where Re= 1/µ. We have

∂u

∂x = sin(πx)π sin(πy) exp(−2π

2_{t/ Re),} ∂u

∂y = − cos(πx) cos(πy)π exp(−2π

2_{t/ Re),} ∂2u ∂x2 = cos(πx)π 2_{sin(πy) exp(−2π}2_{t/ Re),} ∂ 2_u ∂y2 = cos(πx)π 2_{sin(πy) exp(−2π}2_{t/ Re),} ∂v

∂x = cos(πx) cos(πy)π exp(−2π

2_{t/ Re),} ∂v

∂y = − sin(πx)π sin(πy) exp(−2π

2_{t/ Re),} ∂2_v ∂x2 = − sin(πx)π 2_{cos(πy) exp(−2π}2_{t/ Re),} ∂2v ∂y2 = − sin(πx)π 2_{cos(πy) exp(−2π}2_{t/ Re),} ∂p ∂x = 1 2sin(2πx)π exp(−4π 2 t/ Re), ∂p ∂y = 1 2sin(2πy)π exp(−4π 2 t/ Re).

N kuerrork2 Ratio kverrork2 Ratio kperrork2 Ratio

16 2.5618E-05 - 1.7061E-05 - 6.1794E-03 -32 2.7016E-06 9.4825 1.5956E-06 10.6925 1.2137E-03 5.0913 64 2.8951E-07 9.3316 1.7090E-07 9.3364 2.6025E-04 4.6635 128 3.5376E-08 8.1837 2.0645E-08 8.2780 6.0763E-05 4.2830

Table 8: A grid refinement analysis in L2-Norm for Example(1).

N kuerrork∞ Ratio kverrork∞ Ratio kperrork∞ Ratio

16 8.2522E-05 - 6.6571E-05 - 1.0643E-02 -32 1.5195E-05 5.4308 9.4561E-06 7.0400 2.1667E-03 4.9120 64 2.4185E-06 6.2828 1.4674E-06 6.4441 4.8205E-04 4.4947 128 4.2419E-07 5.7014 2.4059E-07 6.0991 1.1171E-04 4.3151

Table 9: A grid refinement analysis in L∞-Norm for Example(1).

[0, 1], and give the exact solution u(x, y, t) = 1 4ye −t , v(x, y, t) = −1 4x(1 − x 2_)e−t , p(x, y, t) = (x2+ y2− 1)2e−t, g1(x, y, t) = − 1 4ye −t₋ 1 16x(1 − x 2_)e−2t + 4(x2+ y2− 1)e−tx, g2(x, y, t) = 1 4x(1 − x 2_)e−t + 1 4ye −t (−1 4(1 − x 2_)e−t +1 2x 2_e−t ) − 3 2xe −t + 4(x2+ y2− 1)e−ty. Here we have ∂u ∂x = 0, ∂u ∂y = 1 4e −t , ∂2_u ∂x2 = 0, ∂2_u ∂y2 = 0, ∂v ∂x = 1 4e −t_{(−1 + 3x}2_), ∂v ∂y = 0, ∂2_v ∂x2 = 3 2xe −t , ∂ 2_v ∂y2 = 0, ∂p ∂x = 4(x 2_{+ y}2_{− 1)e}−t x, ∂p ∂y = 4(x 2_{+ y}2_{− 1)e}−t y.

N kuerrork2 Ratio kverrork2 Ratio kperrork2 Ratio

16 5.3838E-04 - 8.1989E-04 - 1.3604E-02 -32 1.4644E-04 3.6764 2.2699E-04 3.6120 4.3058E-03 3.1594 64 3.8076E-05 3.8459 5.9219E-05 3.8330 1.4074E-03 3.0594 128 9.6009E-06 3.9658 1.4984E-05 3.9521 3.8226E-04 3.6817

Table 10: A grid refinement analysis in L2_{-Norm for Example(2).}

N kuerrork∞ Ratio kverrork∞ Ratio kperrork∞ Ratio

16 1.3922E-03 - 3.2717E-03 - 8.6817E-02 -32 4.0255E-04 3.4584 9.6117E-04 3.4038 4.8699E-02 1.7827 64 1.1010E-04 3.6562 2.5788E-04 3.7271 2.6233E-02 1.8564 128 2.8170E-05 3.9084 6.6552E-05 3.8748 1.3360E-02 1.9635

Table 11: A grid refinement analysis in L∞-Norm for Example(2).

4.3.2 Cavity Flow

Example (3). Let Ω = [0, 1]×[0, 1] be computational domain. Assume that flow is driven by the upper wall. That is, we have the boundary condition

and

u = (u, v) = (1, 0) on Ωlid,

where Ωlid = {(x, y) ∈ ∂Ω, y = 1}. Let N = 128 be the number of partition. We compare

our result with [14].

Figure 7: Figures of velocities u and v on x = 0.5 and y = 0.5 respectively. Re = 1000.

4.3.3 Circular Flow with a Line Force

Example (4). Let us assume the curve Γ to be a circle x2_+y2 ₌ 1

4 within a computational

domain Ω = [−1, 1] × [−1, 1]. That is, the interface can be represented by X(s), Y (s)
= 0.5 cos θ, 0.5 sin θ

where θ = 2s, 0 ≤ s ≤ π. Let the pressure and velocity components be u(x, y, t) = ( h(t) y_r − 2y
r > 1₂ 0 r ≤ 1₂ v(x, y, t) = ( h(t) −x_r + 2x
r > 1₂ 0 r ≤ 1₂ p(x, y, t) = ( 1 r > 1 2 0 r ≤ 1 2

where r = px2_{+ y}2_{. It is easy to verify the velocity satisfies the incompressibility}

constraint, and it is continuous but has finite jump across the interface. That is,  ∂u ∂n  = −2h(t) sin θ,  ∂v ∂n  = 2h(t) cos θ, [p] = 1 and then Fn = 1, Fτ = −2h(t).

The external force g can be computed to satisfy the Navier-Stokes equation, and yields g1(x, y, t) = ( h(t)2cos θ 4 − 4r − 1_r
+ h(t) sin θ_rµ2 + h 0_{(t) sin θ(1 − 2r) r >} 1 2 0 r ≤ 1₂ g2(x, y, t) = ( h(t)2sin θ 4 − 4r − 1_r
− h(t) cos θ_rµ2 − h 0_{(t) cos θ(1 − 2r) r >} 1 2 0 r ≤ 1₂ . One can verify that on Γ, (r = 1/2),

[g] · n = h(t)2  4 − 4r − 1 r  − 0 = 0.

Some results are given here. In our tests, we assume the boundary force is given. We observe that the accuracy in u, v is nearly second order.

N kuk∞ Ratio kvk∞ Ratio kpk∞ Ratio

16 6.2133E-03 0.0000 6.0783E-03 0.0000 1.3309E-01 0.0000 32 2.0049E-03 3.0990 1.9912E-03 3.0526 6.5489E-02 2.0322 64 5.6914E-04 3.5228 5.6695E-04 3.5121 4.4361E-02 1.4762 128 1.4834E-04 3.8367 1.4805E-04 3.8295 1.5033E-02 2.9509 N kukL2 Ratio kvk_L2 Ratio kpk_L2 Ratio

16 1.4411E-03 0.0000 1.4221E-03 0.0000 2.5676E-02 0.0000 32 3.8649E-04 3.7287 3.8360E-04 3.7073 1.0012E-02 2.5645 64 1.1245E-04 3.4370 1.1209E-04 3.4223 3.3265E-03 3.0099 128 3.2281E-05 3.4835 3.2235E-05 3.4772 1.2290E-03 2.7066

Table 12: The grid refinement analysis for u, v, and p. h(t) = 1.

N kuk∞ Ratio kvk∞ Ratio kpk∞ Ratio

16 1.8770E-03 0.0000 1.8741E-03 0.0000 3.7830E-02 0.0000 32 6.3666E-04 2.9483 6.4121E-04 2.9227 2.0364E-02 1.8577 64 1.8387E-04 3.4625 1.8356E-04 3.4932 1.4330E-02 1.4211 128 4.8301E-05 3.8068 4.8899E-05 3.7539 4.7389E-03 3.0238 N kukL2 ratio kvk_L2 ratio kpk_L2 ratio

16 4.3742E-04 0.0000 4.3582E-04 0.0000 7.8480E-03 0.0000 32 1.2249E-04 3.5712 1.2222E-04 3.5660 3.2169E-03 2.4396 64 3.5968E-05 3.4054 3.5932E-05 3.4013 1.0809E-03 2.9762 128 1.0295E-05 3.4937 1.0291E-05 3.4917 4.1351E-04 2.6139 Table 13: The grid refinement analysis for u, v, and p. h(t) = 1 − e−t.

5

Concluding Remarks

We have introduced a numerical technique for the computation of interface problems. The result shows the behavior near the immersed interface can be captured well, and in general second order accuracy can be achieved. Moreover, the method is efficient since fast solvers can still be applied.

One of the future works is to solve piecewise constant coefficient Poisson equation ∇ · (β±∇ψ±) = ω± on Ω±.

The immersed interface method can also handle this type of equation. Thus, we can extend this solver to solve piecewise constant coefficient Stokes equation or Navier-Stokes equation.

It also need to be considered that how the spline we constructed approximates the actual interface, and how accurate the curvature and the distance to the foot of perpen-dicular we obtain. The error generated by approximation of an interface by an spline may lead to the fluctuation of convergent ratio in grid refinement analysis.

It is a challenge that extending this method to three-dimension. Although there are some fast Poisson solvers in 3-D, constructing the interface of surface and finding foot of perpendicular seem to be a heavy task. It is important that maintaining the efficiency and accuracy when extending this method into 3-D.

References

[1] Alfio Quarteroni, Riccardo Sacco, Fausto Saleri, “Numerical Mathematics.” Texts in Applied Mathematics 37, Springer.

[2] Alexandre J. Chorin, Jerrold E. Marsden, “A Mathematical Introduction to Fluid Mechanics, 3rd. ed.” Texts in Applied Mathematics 4, Springer-Verlag.

[3] J. Stoer, R. Bulirsch: “Introduction to Numerical Analysis.” Springer-Verlag, 1983 [4] Bu-Qing Su, Ding-Yuan Liu: “Computational Geometry: Curve and Surfce

Model-ing.” Academic Press, INC.

[5] M. C. Lai, “Jump conditions and test examples for the 2D Navier-Stokes equations involving an immersed boundary.”

[6] R. Cortez, “The Method of Regularized Stokeslets.” SIAM Journal of Scientific Com-puting, vol. 23, No. 4, pp. 1204-1225.

[7] J. Thomas Beale and Anita T. Layton, “On The Accuracy of Finite Difference Meth-ods for Elliptic Problems with Interfaces.” Aug 19, 2005.

[8] Zhilin Li and M. C. Lai, “The Immersed Interface Method for the Navier-Stokes Equations with Singular Forces.” Journal of Computational Physics, Vol. 171, n. 2, pp. 822-842.

[9] A. McKenny and L. Greengard and A. Mayo, “A Fast Poisson Solver for Complex Geometries.” Journal of Computational Physics, Vol. 118, pp. 348-355.

[10] M. C. Lai and C. S. Peskin, “An immersed boundary method with formal second-order accuracy and reduced numerical viscosity. ” Journal of Computational Physics, Vol. 160, pp. 705-719. (2000)

[11] Boyce E. Griffith, Charles S. Peskin, “On the order of accuracy of the immersed boundary method: Higher order convergence rates for sufficiently smooth problems.” Journal of Computational Physics, Vol. 208, pp. 75-105.

[12] D. Russell. Z. Jane Wang, “A cartesian grid method for modeling multiple mov-ing objects in 2D incompressible viscous flow.” Journal of Computational Physics Vol. 191 (2003), pp.177-205.

[13] Sheng Xu, Z. Jane Wang, “An immersed interface method for simulating the in-teraction of a fluid with moving boundaries.” Journal of Computational Physics, in press.

[14] U. Ghia, K. N. Ghia, and C. T. Shin, “High-Re Solutions for Incompressible Flow Us-ing the Navier-Stokes Equations and a Multigrid Method.” Journal of Computational Physics Vol. 48, pp.387-411.