三維橢圓介面問題的數值研究

國

立 交 通 大 學

應 用 數

學 系

碩 士

論 文

三維橢圓介面問題的數值研究

An Immersed Interface Method for 3D

Elliptic Interface Problems

研

究

生

:

郭柏均

指

導 教 授

:

賴明治 教授

三維橢圓介面問題的數值研究

An Immersed Interface Method for 3D Elliptic Interface

Problems

研 究 生

_:

郭柏均

_{Student: Po-Chun Kuo}

指導教授

_:

賴明治 教授

_{Advisor: Prof. Ming-Chih Lai}

國 立 交 通 大 學

應 用 數

學 系

碩

士

論 文

A Thesis

Submitted to Department of Applied Mathematics

College of Science

National Chiao Tung University

in partial Fulfillment of the Requirements

for the Degree of

Master

三維橢圓介面問題的數值研究

學生:郭柏均

指導教授:賴明治教授

國立交通大學應用數學系碩士班

摘

要

在本篇論文中,系使用內嵌介面法把一個二維橢圓介面問題推廣

到三維。此方法為調整過的有限差分法，其修正項透過介面上的不連

續條件而得。藉由由法向量方向所展開的泰勒展開式，在介面附近的

網格點，離散差分方程會被調整。因為在介面上不連續條件型態的關

係， 在解橢圓介面問題的過程中，必須要解一個線性系統。首先用

重新初始化水平集方法來找網格點在界面上的正交投影，接下來我們

運用內嵌介面法，廣義最小殘量方法，和最小平方法來解橢圓介面問

題。

An Immersed Interface Method for

3D Elliptic Interface Problems

Student: Po-Chun Kuo Advisor: Prof. Ming-Chih Lai

Department (Institute) of Applied Mathematics

National Chiao Tung University

Abstract

In this thesis, we extend the immersed interface method for a 2D elliptic interface problem to 3D. The numerical method is a finite difference method modified with some correction terms from the jump conditions. By using Taylor’s expansions along the normal direction, the discrete difference equation is modified at the gird point close to interface. Because of the types of the jump condition, we have to solve a linear system in the process of solving the elliptic interface problem. We first use the reinitialization of level set method to find the orthogonal projection of the grid point and perform it to check its accuracy. Then, we solve the elliptic interface problem by the immersed interface method, GMRES and least squares method, and make some numerical tests to check the rate of convergence.

誌 謝

本論文的完成，首先要感謝我的恩師賴明治教授，在進行研究及

撰寫論文的過程中，老師細心的指導我，提供我許多寶貴的意見並協

助我擬定正確的方向，讓我對偏微分方程和科學計算有更多的認識與

理解。在未來的生涯規劃與求學方向上，老師更提供了各方面具建設

性的建議及看法，能當老師的學生真的是我的福氣。接著我要感謝楊

肅煜教授及舒宇宸教授，在碩士論文口試時指點我口頭報告上不足之

處，還提供許多論文需要修正的地方，老師們的建議可以使本論文更

加週延與完善，在此致上最深的謝意。此外，我要感謝陳冠羽學長、

胡偉帆學長和成澤仕軒學長，在我研究遇到問題或困難時，給了我許

多許多的幫助，認識您們真好。

在大學及研究所求學期間，非常感謝系(所)上的老師們，不只傳

授給我數學相關的知識，也教導我許多探究問題的思考模式和解決問

題的方法，並指點我許多為人處事應有的態度與道理。我也要感謝學

校及捐贈獎學金的團體與人士，讓我在求學期間獲得許多的肯定與鼓

舞。當然，我要感謝我求學期間的好夥伴，也就是系上的同學們，你

們帶給我許多美好的回憶，在我遇到困難時提供我許多的協助。

最後，我要感謝我的家人他們是我精神上最大的支柱，在我背後

默默的支持我、鼓勵我及幫助我，並時時給我鼓勵與叮嚀，讓我在無

後顧之憂的環境下能夠順利的完成這個階段的學業。

Contents

method 2

2.1 Elliptic interface equation . . . 2 2.2 The interpolation of the jump condition . . . 3 3 Find the foot of an irregular point using reinitialization

of level set method 6 3.1 The reinitialization of level set method . . . 6 3.2 Numerical results . . . 7 4 Solve the elliptic interface problem 9 4.1 Solve the elliptic interface equation . . . 9 4.2 Numerical results . . . 13 5 Transfer the elliptic interface equation with zero jump

condition [u] 20 6 Conclusion 22

List of Figures

1 The interface Γ divides the domain into Ω+ _{and Ω}− _{. 2}

2 The seven-point Laplacian of the irregular point xi,j,k 3

3 The foot X∗ _{of the irregular point x . . . . 5}

List of Tables

1 Example 1: The L∞_{errors and the order of accuracy}

of distance, foot and φ on irregular points . . . 7 2 Example 2: The L∞_{errors and the order of accuracy}

of Φ on irregular points . . . 9 3 Example 3: The L∞_{, L}1 _{and L}2 _{errors and the order}

of accuracy with σ+= 2, σ−_{= 23 by method 1 . . . . 14}

4 Example 3: The L∞_{, L}1 _{and L}2 _{errors and the order}

of accuracy with σ+= 23, σ−_{= 2 by method 1 . . . . 14}

5 Example 3: The L∞_{, L}1 _{and L}2 _{errors and the order}

of accuracy with σ+_{= 2, σ}−_{= 23 by method 2 . . . . 14}

6 Example 3: The L∞_{, L}1 _{and L}2 _{errors and the order}

of accuracy with σ+= 23, σ−_{= 2 by method 2 . . . . 14}

7 Example 4: The L∞_{, L}1 _{and L}2 _{errors and the order}

of accuracy with σ+= 2, σ−_{= 23 by method 1 . . . . 16}

8 Example 4: The L∞_{, L}1 _{and L}2 _{errors and the order}

of accuracy with σ+_{= 23, σ}−_{= 2 by method 1 . . . . 16}

9 Example 4: The L∞_{, L}1 _{and L}2 _{errors and the order}

of accuracy with σ+= 2, σ−_{= 23 by method 2 . . . . 16}

10 Example 4: The L∞_{, L}1 _{and L}2 _{errors and the order}

of accuracy with σ+= 23, σ−_{= 2 by method 2 . . . . 16}

11 Example 5: The L∞_{, L}1 _{and L}2 _{errors and the order}

of accuracy with σ+_{= 2, σ}−_{= 23 . . . . 18}

12 Example 5: The L∞_{, L}1 _{and L}2 _{errors and the order}

of accuracy with σ+= 23, σ−_{= 2 . . . . 18}

13 Example 6: The L∞_{, L}1 _{and L}2 _{errors and the order}

of accuracy with σ+= 2, σ−_{= 23 . . . . 19}

14 Example 6: The L∞_{, L}1 _{and L}2 _{errors and the order}

1

Introduction

Elliptic interface problems have many applications in science and en-gineering, such as multi-material problems, two-phase problems, fluid problems. We focus on elliptic interface problems with variable coeffi-cients. Along the interface, there exist the jump conditions of solution and flux across the interface. For example, this can applied in electro-hydrodynamics [2, 15].

Peskin used the immersed boundary method [1] to simply describe the fluid motion interacting with complicated interface at the Carte-sian grid. Many methods can solve elliptic interface problems, such as finite element method [5, 10], boundary integral method [13], ghost fluid method [3, 12], sharp interface method [7].

In this thesis, we use the immersed interface method [8, 15] to solve the problem from 2D to 3D. We have to find the orthogonal projec-tion to the interface. In 2D, we can use polynomials to interpolate the interface by picking some point and use these approximation to find the orthogonal projection, such as cubic spline representation [4]. However, it is difficult to use this method in 3D. Therefore, we use reinitialization of level set method [9, 11, 14].

In Section 2, we introduce the elliptic interface equation to study how to use the immersed interface method to solve this equation. In Section 3, we can see how to find the orthogonal projection to the interface by reinitialization of level set method. In Section 4, since we do not have [un], we have to get the solution through Poisson

solver and solve a linear system. In Section 5, a method to modify the equation s.t. [u] = 0 is introduced. Finally, some conclusions are displayed in Section 6.

2

Elliptic interface equation and the

immersed interface method

2.1

Elliptic interface equation

We consider a rectangular domain Ω = [xℓ, xr] × [yℓ, yr] × [zℓ, zr] with

a closed interface Γ. Since the interface Γ divides the domain Ω into two regions, we denote the inside region by Ω−_{and the outside region}

by Ω+, and with piecewise constant coefficient, σ− _{on Ω}− _{and σ}+ _on

Ω+. Along Γ, there exist the jump conditions of solution and flux across Γ.

We solve the following elliptic interface equation        ▽ · (σ ▽ u) = f in Ω − Γ, [u](X) = ω(X) on Γ, [σun](X) = ν(X) on Γ, u = ub on ∂Ω (1) where un= ∂u_∂n.

Since σ is a piecewise constant, we can rewrite the elliptic interface equation as ∆u = f in Ω − Γ where f = f /σ, i.e.        ∆u = f _{in Ω − Γ,} [u](X) = ω(X) on Γ, [σun](X) = ν(X) on Γ, u = ub on ∂Ω. (2)

2.2

The interpolation of the jump condition

To solve this problem, we use the immersed interface method devel-oped in [8, 15]. We use a uniform Cartesian grid in Ω with the mesh width h = ∆x = ∆y = ∆z, and use xi,j,k = (xl+ i∆x, yl+ j∆y, zl+

k∆z) to denote the mesh grid point and denote discretized solution at the mesh grid point by ui,j,k. Since we have the jump condition

on the interface Γ, we determine whether the grid point is a regular point or an irregular point. If the seven-point Laplacian of a grid point uses some points in Ω− _{and some points in Ω}+ _{simultaneously,}

the grid point is called the irregular point. Otherwise, when the grid point uses either points in Ω− _{or points in Ω}+_{, it is the regular point.}

At the irregular point, because of the jump conditions, we have to change the approximations according to the jump condition on the interface. For example, we let xi,j,k, xi+1,j,k, xi,j+1,k, xi,j,k+1 be in the

inside region, xi−1,j,k, xi,j−1,k, xi,j,k−1 in the outside region, then xi,j,k

is an irregular point, and the modified seven-point Laplacian ∆hu is

written as follows.

∆hu(xi,j,k) :=

ui−1,j,k− 2ui,j,k+ ui+1,j,k

h2

+ui,j−1,k− 2ui,j,k+ ui,j+1,k

h2

+ui,j,k−1− 2ui,j,k+ ui,j,k+1

h2 = u + i−1,j,k− 2u − i,j,k+ u − i+1,j,k h2 +u + i,j−1,k− 2u−i,j,k+ u − i,j+1,k h2 +u + i,j,k−1− 2u − i,j,k+ u − i,j,k+1 h2 = u − i−1,j,k− 2u − i,j,k+ u − i+1,j,k h2 + u+_{i−1,j,k− u}−_i−1,j,k h2 +u − i,j−1,k− 2u−i,j,k+ u − i,j+1,k h2 + u+_{i,j−1,k− u}−_i,j−1,k h2 +u − i,j,k−1− 2u−i,j,k+ u − i,j,k+1 h2 + u+_{i,j,k−1− u}−_i,j,k−1 h2

= uxx(xi,j,k) + uyy(xi,j,k) + uzz(xi,j,k) + O(h2)

+u c i−1,j,k h2 + uc_i,j−1,k h2 + uc_i,j,k−1 h2 = f_i,j,k+ 1 h2 u c

i−1,j,k+ uci,j−1,k+ uci,j,k−1
+ O(h2)

(3) where uc

i,j,k = u+i,j,k− u

− i,j,k.

To compute the correction term uc_i,j,k, we define the foot X∗

i,j,k of

xi,j,k. It means that X∗_i,j,kis the orthogonal projection to the interface

of xi,j,k. By Taylor’s expression, we can write

uc_i,j,k = u+_{i,j,k− u}−_i,j,k

=  u++ α∂u + ∂n + α2 2 ∂2_u+ ∂n2     X∗ i,j,k + O(h3) −  u−_{+ α}∂u− ∂n + α2 2 ∂2_u− ∂n2     X∗ i,j,k + O(h3) = [u]_X∗ i,j,k + α  ∂u ∂n  X∗ i,j,k +α 2 2  ∂2u ∂n2  X∗ i,j,k + O(h3). (4)

Figure 3: The foot X∗ _{of the irregular point x}

The |α| is the distance between xi,j,k and X∗_i,j,k. If xi,j,k is in the

inside region, then α is positive. Otherwise, if xi,j,k is in the outside

region, then α is negative. (i.e., α is the signed distance between xi,j,k

and X∗ i,j,k) Instead of findingh∂_∂n2u2 i X∗ i,j,k

, we use this equation

∆u = ∂

2_u

∂n2 + κ

∂u

∂n + ∆ssu (5) where κ = 2H = ∇ · n, H is the mean curvature. ∆ssu is the surface

laplacian operator of u. Then, we can write  ∂2u ∂n2  X∗ i,j,k = f _X∗ i,j,k− κX ∗ i,j,k  ∂u ∂n  X∗ i,j,k − ∆ss[u]X∗ i,j,k. (6)

Hence, the correction term can be written as

uci,j,k = [u]X∗ i,j,k+ α  ∂u ∂n  X∗ i,j,k +α 2 2  ∂2_u ∂n2  X∗ i,j,k + O(h3) = [u]X∗ i,j,k+ α  ∂u ∂n  X∗ i,j,k

3

Find the foot of an irregular point

using reinitialization of level set method

3.1

The reinitialization of level set method

For finding irregular points, we use reinitialization of level set method developed in [9, 11, 14]. First, we have to find signed-distance function φ : Ω → R where φ(x) =    d(x, Γ) in Ω+, 0 on Γ, −d(x, Γ) on Ω−_. (8)

Now, we consider a level set φ0 : Ω → R where

φ0(x)    > 0 in Ω+_, = 0 on Γ, < 0 on Ω−_. (9)

using reinitialization equation

φt+ S(φ0)(| ▽ φ| − 1) = 0 (10)

where S(φ0) is signed function. (i.e. S(φ0(x)) = 1 when x ∈ Ω+,

S(φ0(x)) = −1 when x ∈ Ω− and S(φ0(x)) = 0 when x ∈ Γ.)

It will converge to signed-distance function. Since the equation con-verges steady state quickly very much at the points near the Γ, we only use a few time steps to get the results. For getting better results, we modify the S(φ0) as S(φ0) = φ0 pφ2 0+ (∆x)2 . (11) as a numerical approximation. In this thesis, we use TVD-RK and WENO to compute Eq. (10). Now, we find the foot X∗ _{of x. Since}

X∗ is the foot of x, we have the equation

X∗ _{= x − φ(x)∇φ(X}∗). (12) Since x is close to X∗_{, we use this following approximation for}

sim-plicity

3.2

Numerical results

In these numerical experiments, we consider the rectangular domain Ω = [−1, 1] × [−1, 1] × [−1, 1]. In Example 1 and 2, we perform the numerical tests for finding the foots by the reinitialization of level set method.

Example 1

In this example, we consider the interface Γ as a sphere. The surface equation is

Σ : x2+ y2+ z2 = r2₀ (14) where r0 =r 2

13. The initial level set function is chosen to be

φ0(x, y, z) = x2+ y2+ z2− r02. (15)

The following numerical results are the errors about finding the foot by reinitialization of level set method. The length of the domain is 2 so we choose h = _N2 with N = 32, 64, 128, 256. Since the interface is a sphere, the exact distance and the exact foot of the irregular points are easy to obtain. Besides, we can choose a level set function φ to describe the error of the foot to the interface. Since Γ is a sphere, the signed-distance function φ is known. φ is suited to be the level set function describing the error of the foot. Therefore, we use φ to be

φ(x, y, z) =px2_{+ y}2_{+ z}2_{− r}

0. (16)

The following table shows the L∞ _{errors at the irregular points and}

the steps of the iteration. The errors are the distance, the foot, and φ for different N . In this case, we can see that the method has the first order accuracy.

Example 2

In this example, we consider the interface Γ as an ellipsoid. The surface equation is Σ : x 2 a2 + y2 b2 + z2 c2 = 1 (17) where a =r 2 13, b = r 11 43, c = r 3 31.

A problem occurs when the initial level set function is chosen to be φ0(x, y, z) = x 2 a2 + y2 b2 + z2 c2 − 1. (18)

Using this φ0 to compute φ, a state happens in the iteration process.

For some iteration step n, there exists some xi,j,k, s.t.

φ0(xi,j,k)φn(xi,j,k) < 0, (19)

i.e. for some grid points in the outside region, they are translated to the inside region in the iteration process, or the opposite situation. To deal with this problem, we discuss |∇φ| − 1 near the interface. On Γ, |∇φ(x, y, z)| =   2x a2, 2 y b2, 2 z c2    = 2 r x2 a4 + y2 b4 + z2 c4 ≥ 2 s 1 b2  x2 a2 + y2 b2 + z2 c2  = 2 b r x2 a2 + y2 b2 + z2 c2 = 2r 43 11 ≃ 4. (20) Hence, the reason of this occasion may be that |∇φ| − 1 near the interface is too large.

Here, we scale the initial level set function as φ0(x, y, z) =  a2_{+ b}2_{+ c}2 3   x2 a2 + y2 b2 + z2 c2 − 1  . (21) We take h = _N2 with N = 32, 64, 128, 256. Since the exact distance and the exact foot of the irregular points can not be obtained directly,

we describe the error of the foot to the interface by level set function Φ. The level set function is chosen to be

Φ(x, y, z) = x 2 a2 + y2 b2 + z2 c2 − 1. (22)

Table 2 shows the L∞_{errors of Φ at the irregular points and the steps}

of the iteration for different N . In this table, we can see that the rate of convergence is first order.

Size Φ Order step 32 1.1443e-01 - 10 64 5.4037e-02 1.0824 10 128 2.5760e-02 1.0688 10 256 1.2700e-02 1.0203 10

Table 2: The L∞ _{errors and the order of accuracy of Φ on}

irregular points for different N

4

Solve the elliptic interface problem

4.1

Solve the elliptic interface equation

We follow these steps developed in [15] to solve the equation (2). Since we only know [σun] instead of [un], we use these equations

[un] = 1 σ− [σun] − [σ]u +
when σ−_{> σ}+_{, (23)} [un] = 1 σ+ [σun] − [σ]u −
when σ+> σ−. (24) We will explain why we choose the equation in Eq. (44).

We let Ci,j,k be the correction term, then we can write Eq. (3) as

∆hui,j,k = fi,j,k+ Ci,j,k. (25)

And we let Ci,j,k be the correction term without [un] term, i.e., Ci,j,k

u−_{from the grid points in Ω}−_{. The approximation method is written}

in the following.

Suppose we want to approximate u+(X∗

i0,j0,k0). First, we set a

square number SN . Then, we collect a set S written as S = {xi,j,k ∈ Ω+|i = i0− sn, . . . , i0+ (SN − sn − 1),

j = j0− sn, . . . , j0+ (SN − sn − 1),

k = k0− sn, . . . , k0+ (SN − sn − 1)} (28)

where sn is the largest integer not greater than SN −1₂ .

Figure 4: The set S of X∗ _{with SN = 7}

Given a polynomial P3(x, y, z) with degree 3, we determine the

coeffi-cients of P3(x, y, z) by least squares method approximating u+(x, y, z)

on S. So we approximate u+(X∗

i0,j0,k0) by

u+_n(X∗_{i0,j0,k0) ≈ ∇P}3(X∗i0,j0,k0) · n(X

∗

i0,j0,k0). (29)

Similarly, we can approximate u−_(X∗

i0,j0,k0) by the above method.

Since the approximation is linear for ui,j,k, we can obtain two linear

operators B+_{, B}− _s.t.

u+_{n ≈ B}+u, (30) u−n ≈ B−u. (31)

And with the approximation, we rewrite Eq. (23) and (24) as [un] = 1 σ− [σun] − [σ]B +_u
when σ−> σ+, (32) [un] = 1 σ+ [σun] − [σ]B −_u
when σ+> σ−_{. (33)}

B+u +σ − [σ][un] = [σun] [σ] , (34) B−u +σ + [σ][un] = [σun] [σ] . (35) From Eq. (27), (34) and (35), we get a linear system as

" ∆h −A B± σ∓ [σ] #  u [un]  = " f + C [σun] [σ] # . (36) To solve Eq. (36), we use Schur-complement method, i.e., we solve u and [un] separately. In this method, we have to use Poisson solver.

In this thesis, we use the public software package MUDPACK to solve Poisson equation. The details are written as follows.

First, we solve a Poisson equation u∗ _for

∆hu∗ = f + C, (37)

u∗ _{= u}

b on ∂Ω. (38)

Second, we compare u and u∗

∆hu − A [un] = f + C, (39) ∆hu∗ = f + C, (40) u = u∗ _{= u} b on ∂Ω. (41) ∆h(u − u∗) − A [un] = 0, (42) u − u∗ _{= 0 on ∂Ω. (43)} (u − u∗_{) − ∆}−1 h A [un] = 0, u − ∆−1h A [un] = u ∗_, B±_{u − B}±∆−1_h A [un] = B±u∗, −σ ∓ [un] + [σun] − B±∆−1_h A [un] = B±u∗,

this is, we do not have to construct the matrix −B±_∆−1

h A. We can

use the poission solver to replace ∆−1_h , and do not need to write the explicit form of A and B±_{. For the iteration, we let the stopping}

crite-ria be h2. We can do this setting for the following reasons. The error of [un] is O(h2) and d is O(h), so the error of [un] makes an error O(h)

in the correction term. The correction term originally has an error O(h). Hence, The error of [un] does not affect the overall accuracy.

Once we know [un], then we can solve u by

∆hu = f + C + A [un] , (45)

u = ub on ∂Ω. (46)

In summary, we write the following steps of the method solving Eq. (36).

Step 1.

Solving the Poisson equation u∗ _for

∆hu∗ = f + C, (47)

u∗ = ub on ∂Ω. (48)

Step 2.

Using GMRES with the stopping criteria h2 to solve [un] in the below

linear system  −B±∆−1_{h A −}σ ∓ [σ]I  [un] = B±u∗− [σun] [σ] (49) where ∆−1_h is poisson solver with zero Dirichlet boundary condition. Step 3.

Solving the Poisson equation u

∆hu = f + C + A [un] , (50)

4.2

Numerical results

Considering the rectangular domain Ω = [−1, 1] × [−1, 1] × [−1, 1], we display the numerical results of the immersed interface method for the elliptic interface equation (2).

Example 3

We consider an interface and an exact solution as

Σ : x2+ y2+ z2= r2₀, (52) u+= x2+ y2+ z2, (53) u−_{= sin(x}2_{+ y}2_{+ z}2_{) (54)} where r0 = r 2 13.

The force f , curvature κ and the boundary conditions can be calcu-lated as                u+_n = 2px2_{+ y}2_{+ z}2_, u− n = 2px2+ y2+ z2cos(x2+ y2+ z2), f+ = 6, f− = 6 cos(x2_{+ y}2_{+ z}2_{) − 4(x}2_{+ y}2_{+ z}2_{) sin(x}2_{+ y}2_{+ z}2_), κ =√ 2 x2_+y2_+z2. (55) Therefore, the equation (2) is written as

           ∆u = 6 in Ω+, ∆u = 6 cos(x2+ y2+ z2_{) − 4(x}2+ y2+ z2) sin(x2+ y2+ z2) in Ω−_,

[u] = r2 0− sin(r02) on Γ, [σun] = 2r0 σ+− σ−cos(r02)
on Γ, u = x2+ y2+ z2 on ∂Ω. (56) Here, [u] is constant on Γ, so ∆ss[u] = 0. We take h = _N2 with

errors, the order of accuracy and the iteration steps for different N , σ+_{, σ}− _{are displayed in the following tables.}

Method 1

Size L∞ _{Order L}1 _{Order L}2 _{Order step}

32 7.6125e-05 - 8.0125e-06 - 1.4937e-05 - 2 64 3.1055e-05 1.2936 3.4919e-06 1.1982 6.2259e-06 1.2625 3 128 1.0722e-05 1.5343 6.8645e-07 2.3468 1.4433e-06 2.1089 5

Table 3: The L∞_{, L}1 _{and L}2 _{errors and the order of accuracy for different}

N with σ+_{= 2, σ}− _{= 23 by method 1}

Size L∞ _{Order L}1 _{Order L}2 _{Order step}

32 1.0930e-03 - 1.7164e-04 - 2.9306e-04 - 2 64 2.6165e-04 2.0625 4.0319e-05 2.0898 6.7259e-05 2.1234 2 128 2.1928e-05 3.5768 2.9699e-06 3.7630 4.8823e-06 3.7841 2 256 1.1043e-06 4.3116 6.8618e-08 5.4357 1.2182e-07 5.3247 2

Table 4: The L∞_{, L}1 _{and L}2 _{errors and the order of accuracy for different}

N with σ+ _{= 23, σ}− _{= 2 by method 1}

Method 2

Size L∞ _{Order L}1 _{Order L}2 _{Order step}

32 1.1487e-04 - 1.4585e-05 - 2.6055e-05 - 2 64 3.7705e-05 1.6071 5.1938e-06 1.4896 8.9480e-06 1.5419 2 128 1.1889e-05 1.6652 1.5597e-06 1.7356 2.6672e-06 1.7463 3 256 3.3990e-06 1.8064 4.9160e-07 1.6657 8.2544e-07 1.6921 4

Table 5: The L∞_{, L}1 _{and L}2 _{errors and the order of accuracy for different}

N with σ+_{= 2, σ}− _{= 23 by method 2}

Size L∞ _{Order L}1 _{Order L}2 _{Order step}

32 7.4238e-04 - 1.1580e-04 - 1.9781e-04 - 1 64 1.8352e-04 2.0162 2.8322e-05 2.0316 4.7207e-05 2.0670 2 128 1.4103e-05 3.7018 1.8714e-06 3.9197 3.0544e-06 3.9500 2 256 1.1318e-06 3.6393 8.3163e-08 4.4920 1.6270e-07 4.2306 3

Table 6: The L∞_{, L}1 _{and L}2 _{errors and the order of accuracy for different}

Example 4

In this example, an interface and an exact solution are chosen to be Σ : x2+ y2+ z2= r₀2, (57) u+= ex+y+z, (58) u−= sin(x + y + z) (59) where r0 = r 2 13.

The force f , curvature κ and the boundary conditions are calculated as                  u+ n = ex+y+z√_xx+y+z2_+y2_+z2, u− n = sin(x + y + z)√_xx+y+z2_+y2_+z2, f+ = 0, f− _{= −3 sin(x + y + z),} κ = √ 2 x2_+y2_+z2. (60)

The equation (2) is written as            ∆u = ex+y+z in Ω+, ∆u = −3 sin(x + y + z) in Ω−_,

[u] = ex+y+z _{− sin(x + y + z) on Γ,}

[σun] = x+y+z_r0 (σ+ex+y+z− σ−cos(x + y + z)) on Γ,

u = ex+y+z on ∂Ω.

(61)

Since the interface Γ is a sphere, ∆ss[u] is available.

∆ss[u] = 2 cos(x + y + z) x + y + z r2 0 + 2 sin(x + y + z)  1 −xy + yz + zx r2 0  + 2ex+y+z  1 −xy + yz + zx + x + y + z r2 0  . (62) We choose h = _N2 with different N = 32, 64, 128, 256. Since the

inter-tables show the L∞_{, L}1_{, L}2 _{errors, the order of accuracy and the}

it-eration steps. Method 1

Size L∞ _{Order L}1 _{Order L}2 _{Order step}

32 4.7880e-03 - 5.0982e-04 - 9.5825e-04 - 4 64 3.7523e-03 0.3516 4.8873e-04 0.0610 8.4660e-04 0.1787 6 128 1.9643e-03 0.9338 2.5085e-04 0.9622 4.2371e-04 0.9986 8 256 5.1709e-04 1.9255 6.4666e-05 1.9558 1.0849e-04 1.9655 10

Table 7: The L∞_{, L}1 _{and L}2 _{errors and the order of accuracy for different}

N with σ+_{= 2, σ}− _{= 23 by method 1}

Size L∞ _{Order L}1 _{Order L}2 _{Order step}

32 6.3817e-04 - 9.7196e-05 - 1.3414e-04 - 4 64 1.9205e-04 1.7325 2.4062e-05 2.0141 3.1806e-05 2.0764 4 128 3.2468e-05 2.5644 6.4305e-06 1.9037 9.0349e-06 1.8157 4 256 1.0633e-05 1.6104 1.6876e-06 1.9299 2.4299e-06 1.8946 5

Table 8: The L∞_{, L}1 _{and L}2 _{errors and the order of accuracy for different}

N with σ+_{= 23, σ}−_{= 2 by method 1}

Method 2

Size L∞ _{Order L}1 _{Order L}2 _{Order step}

32 2.5791e-03 - 2.0105e-04 - 4.2491e-04 - 3 64 1.3330e-03 0.9522 1.2981e-04 0.6311 2.4516e-04 0.7935 5 128 5.2778e-04 1.3367 4.1144e-05 1.6577 7.8304e-05 1.6465 7 256 1.6009e-04 1.7211 1.0226e-05 2.0085 2.0025e-05 1.9673 8

Table 9: The L∞_{, L}1 _{and L}2 _{errors and the order of accuracy for different}

N with σ+_{= 2, σ}− _{= 23 by method 2}

Size L∞ _{Order L}1 _{Order L}2 _{Order step}

32 7.0015e-04 - 1.0703e-04 - 1.4639e-04 - 3 64 1.8356e-04 1.9314 2.3551e-05 2.1842 3.2243e-05 2.1828 4 128 8.5021e-05 1.1103 5.0232e-06 2.2291 8.0645e-06 1.9993 5 256 7.4797e-05 0.1848 2.0646e-06 1.2827 4.6048e-06 0.8085 6

Table 10: The L∞_{, L}1_{and L}2 _{errors and the order of accuracy for different}

Example 5

Considering an interface and an exact solution as Σ : x 2 a2 + y2 b2 + z2 c2 = 1, (63) u+= x 2 a2 + y2 b2 + z2 c2, (64) u−= sin x 2 a2 + y2 b2 + z2 c2  , (65) where a =r 2 13, b = r 11 43, c = r 3 31.

We calculate the force f , curvature κ and the boundary conditions                                u+ n = 2 q x2 a4 + y2 b4 +z 2 c4, u− n = 2 q x2 a4 + y2 b4 +z 2 c4 cos  x2 a2 + y2 b2 +z 2 c2  , f+ = 2 _a12 +_b12 +_c12
, f− = 2 _a12 +_b12 +_c12
cos  x2 a2 + y2 b2 + z2 c2  −4xa24 + y2 b4 +z 2 c4  sinx_a22 + y2 b2 +z 2 c2  , κ =  1 a2+ 1 b2+ 1 c2  x2 a4+ y2 b4+ z2 c4  −  x2 a6+ y2 b6+ z2 c6   x2 a4+ y2 b4+ z2 c4 3₂ . (66)

Therefore, the equation (2) is written as                ∆u = 2 _a12 +_b12 +_c12
in Ω+, ∆u = 2 _a12 +_b12 +_c12
cos  x2 a2 + y2 b2 +z 2 c2  in Ω−_, [u] = 1 − sin(1) on Γ, [σun] = 2 q x2 a4 + y2 b4 + z2 c4 (σ+− σ−cos(1)) on Γ, u = x_a22 + y2 b2 + z2 c2 on ∂Ω. (67)

Here, [u] is constant on Γ, so ∆ss[u] = 0. We take h = _N2 with

Size L∞ _{Order L}1 _{Order L}2 _{Order step}

32 2.5842e-02 - 2.8891e-03 - 5.2446e-03 - 3 64 8.8586e-03 1.5445 9.8525e-04 1.5521 1.7558e-03 1.5787 6 128 3.2733e-03 1.4363 3.2026e-04 1.6213 5.8850e-04 1.5771 9 256 9.5019e-04 1.7844 1.0107e-04 1.6639 1.7713e-04 1.7323 11

Table 11: The L∞_{, L}1_{and L}2 _{errors and the order of accuracy for different}

N with σ+_{= 2, σ}− _{= 23}

Size L∞ _{Order L}1 _{Order L}2 _{Order step}

32 3.6022e-01 - 4.2028e-02 - 7.2799e-02 - 4 64 9.3078e-02 1.9524 9.8860e-03 2.0879 1.6791e-02 2.1162 6 128 7.7633e-03 3.5837 6.2775e-04 3.9771 1.0806e-03 3.9578 7 256 3.6076e-04 4.4275 3.2838e-05 4.2567 6.1427e-05 4.1368 8

Table 12: The L∞_{, L}1_{and L}2 _{errors and the order of accuracy for different}

N with σ+_{= 23, σ}−_{= 2}

Example 6

In this example, we consider an interface and an exact solution as Σ : x 2 a2 + y2 b2 + z2 c2 = 1, (68) u+= ex+y+z, (69) u−= sin(x + y + z) (70) where a =r 2 13, b = r 11 43, c = r 3 31.

We calculate the force f , curvature κ and the boundary conditions                            u+_n = x a2+ y b2+ z c2 q x2 a4+ y2 b4+ z2 c4 ex+y+z, u− n = x a2+ y b2+ z c2 q x2 a4+ y2 b4+ z2 c4 cos (x + y + z) , f+ = 3ex+y+z, f− _{= −3 sin(x + y + z),} κ =  1 a2+ 1 b2+ 1 c2  x2 a4+ y2 b4+ z2 c4  −  x2 a6+ y2 b6+ z2 c6   x2 a4+ y2 b4+ z2 c4 3₂ . (71)

Therefore, the equation (2) is written as                ∆u = 3ex+y+z in Ω+, ∆u = −3 sin(x + y + z) in Ω−_,

[u] = ex+y+z_{− sin(x + y + z)} on Γ, [σun] = x a2+ y b2+ z c2 q x2 a4+ y2 b4+ z2 c4 (σ+ex+y+z _{− σ}−_{cos(x + y + z)) on Γ,} u = ex+y+z on ∂Ω. (72) Since the interface Γ is an ellipsoid, ∆ss[u] is available.

∆ss[u] = ex+y+z 2 − 1 a2 +_b12 +_c12
α3+ 2α1 α2 +α3α4 α2 2 ! − 2  1 − α_α1 2  sin(x + y + z) −α_α32 2  1 a2 + 1 b2 + 1 c2  α2− α4  cos(x + y + z) (73) where α1= xy a2_b2 + yz b2_c2 + zx c2_a2, (74) α2= x2 a4 + y2 b4 + z2 c4, (75) α3= x a2 + y b2 + z c2, (76) α4= x2 a6 + y2 b6 + z2 c6. (77)

We choose σ+ = 2, σ− _{= 23 and σ}+ _{= 23, σ}− _{= 2 to be the cases}

in this example. Since we know the exact u , the error between the numerical uN and the exact u is available. The following tables show

the L∞_{, L}1_{, L}2 _{errors, the order of accuracy and the iteration steps}

for different N , σ+, σ−_.

Size L∞ _{Order L}1 _{Order L}2 _{Order step}

32 1.1158e-02 - 3.5443e-04 - 8.5945e-04 - 5 64 4.0653e-03 1.4566 1.7419e-04 1.0248 3.8375e-04 1.1632 5 128 1.6554e-03 1.2962 6.3940e-05 1.4459 1.3586e-04 1.4981 7 256 8.2041e-04 1.0127 2.9861e-05 1.0984 6.2598e-05 1.1179 8

Table 14: The L∞_{, L}1_{and L}2 _{errors and the order of accuracy for different}

N with σ+_{= 23, σ}−_{= 2}

5

Transfer the elliptic interface

equa-tion with zero jump condiequa-tion

[u]

Now, we consider an elliptic interface equation        ∆u = f _{in Ω − Γ,} [u](X) = ω(X) on Γ, [σun](X) = ν(X) on Γ, u = ub on ∂Ω (78) where ω ∈ C2(Γ).

In previous sections, if the equation is solved, we must calculate ∆ss[u].

When the interface is known, calculating ∆ss[u] is feasible. However,

for some problems, the interface will move with time. When the in-terface moves in the problem, the inin-terface is possible to become un-known. Then, calculating ∆ss[u] become a difficult problem.

If we do not want to calculate ∆ss[u], then we must transfer the

equation s.t. [u] = 0. To do that, we can use this method in [16]. First, we get the signed-distance function φ(x) as follows

φ(x) =    d(x, Γ) in Ω+, 0 on Γ, −d(x, Γ) on Ω−_. (79)

Let us suppose φ ∈ C3_{(Ω) and the normal lines of the interface Γ do}

not intersect in the outside region Ω+, then every x in the outside region has only one foot X∗_{, i.e.}

We can define an extension ωe(x) of ω(x) in Ω+∪ Γ

ωe(x) = ωe(X∗+ φ(x)∇φ(X∗))

= ω(X∗) (81) and ωe(x) = 0 in Ω−.

Since φ ∈ C3_{(Ω) and ω ∈ C}2_{(Γ), ω}

e∈ C2(Ω \ Γ). Now, we define ˆu(x)

ˆ u(x) = H(φ(x))ωe(x) =    0 φ(x) < 0, 1 2ωe(x) φ(x) = 0, ωe(x) φ(x) > 0. (82) where H(·) is the Heaviside function.

Let q(x) = u(x) − ˆu(x), then we can get this equation        ∆q = f − H(φ)∆ˆu in Ω − Γ, [q] = 0 on Γ, [σqn](X) = ν(X) on Γ, q = ub− ˆu on ∂Ω. (83) Proof :

If x ∈ Ω−_{, then ˆu(x) = 0. Therefore, ∆ˆu(x) = 0, and H(φ(x)) = 0.}

Hence, ∆q(x) = ∆u(x) − ∆ˆu(x) = f (x) − 0 = f (x) − H(φ(x))∆ˆu(x). (84) If x ∈ Ω+_{, then H(φ(x)) = 1. Hence,} ∆q(x) = ∆u(x) − ∆ˆu(x) = f (x) − H(φ(x))∆ˆu(x). (85) Then, when X ∈ Γ, [q](X) is [q](X) = [u](X) − [ˆu](X)

Hence, [σ ˆun](X) = σ+ ∂ ˆu ∂n     X+− σ −∂ ˆu ∂n     X− = σ+∂ωe ∂n     X+− 0 = 0, (88) [σqn](X) = [σun](X) − [σˆun](X) = ν(X) − 0 = ν(X). (89)

6

Conclusion

In this thesis, we use the immersed interface method to solve the 3D elliptic interface problem, and reinitialization of level set method to find the foot of the irregular point. To deal with the problem, we only use a simple interpolation, Poisson solver, least squares method and GMRES. It is simple to solve the elliptic interface problem by the present method. There only needs a few GMRES iteration steps. About approximation of u±_{, using polynomial with degree 3 by least}

squares method makes the rates of convergence of the numerical re-sults in this thesis are better than first order. Even some rere-sults have high order accuracy.

References

[1] Charles S. Peskin. The immersed boundary method, Acta Nu-merica 11: 479-517 (2002)

[2] D. A. Saville. ELECTROHYDRODYNAMICS: The Taylor-Melcher Leaky Dielectric Model, Annual Review of Fluid Me-chanics Vol. 29: 27-64 (1997).

[3] Hiroshi Terashima , Gr´etar Tryggvason A front-tracking/ghost-fluid method for front-tracking/ghost-fluid interfaces incompressible flows, Journal of Computational Physics Volume 228, (2009).

[4] Hsiao-Chieh Tseng. Numerical methods and applications for im-mersed interface problem, Master thesis, National Chiao-Tung University, Taiwan, Republic of China, (2006).

[5] James H. Bramble, J. Thomas King. A finite element method for interface problems in domains with smooth boundaries and interfaces, Advances in Computational Mathematics 6: 109-138 (1996).

[6] Jian-Jun Xu, Hong-Kai Zhao. An eulerain formulation for solving partail differenctial equation along a moving interface, Journal of Scientific Computing 19(1-3)(2003) p.573-594.

[7] M. Sussman, K.M. Smith, M.Y. Hussaini, M. Ohta, R. Zhi-Wei. A sharp interface method for incompressible two-phase flows, Jour-nal of ComputatioJour-nal Physics Volume 221, (2007)

[8] Ming-Chih Lai, Hsiao-Chieh Tseng. A simple implementation of the immersed interface methods for Stokes flows with singular forces , Computers and Fluids , vol 37, 99-106, (2008).

[9] Narisawa Shinoki. An implicit closest point method for solving convection-diffusion equations on a moving surface, Master the-sis, National Chiao-Tung University, Taiwan, Republic of China, (2013).

[10] Songming Hou,Peng Song, Liqun Wang, Hongkai Zhao. A Weak Formulation for Solving Elliptic Interface Problems Without Body Fitted Grid Journal of Computational Physics Volume 249, (2013).

[11] Stanley Osher, Ronald Fedkiw. Level set methods and dynamic implicit surfaces, (2003)

[12] T.G. Liu, B.C. Khoo, K.S. Yeo. Ghost fluid method for strong shock impacting on material interface, Journal of Computational Physics Volume 190, (2003).

[13] T.Y.Hou, J.S.Lowengrub, and M.J.Shelley. Boundary Integral Methods for Multicomponent Fluids and Multiphase Materials, Journal of Computational Physics Volume 169, (2001)

[14] Yu-Lun Chang. A simple immersed interface method for 3D Pois-son equation with jump condition, Master thesis, National