模擬可溶性界面活性劑問題之沉浸邊界法

國

立

交

通

大

學

應用數學系

博

士

論

文

模擬可溶性介面活性劑問題之沉浸邊界法

An immersed boundary method for simulating the interfacial flows

with soluble surfactant

研 究 生： 陳冠羽

指導教授： 賴明治 教授

模擬可溶性界面活性劑問題之沉浸邊界法

An immersed boundary method for simulating the interfacial flows

with soluble surfactant

研 究 生：陳冠羽

Student：Kuan-Yu Chen

指導教授：賴明治 教授

Advisor：Dr. 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 Doctor of Philosophy

in

Applied Mathematics

November 2013

Hsinchu, Taiwan, Republic of China

模擬可溶性界面活性劑問題之沉浸邊界法

研究生：陳冠羽

指導教授

：賴明治 教授

國立交通大學

應用數學系 博士班

摘

要

本論文之第一部分在探討沉浸邊界法中求解壓力項或是指示函數之精度，我

們提供一維的理論證明與二維的數值結果來說明 L

1

範數為一階收斂，L

2

範數為

半階收斂，而 L

∞

範數則存在 O(1)的誤差。我們也將討論帶有其他種類界面奇異

項的帕松方程式(Poisson Equation)，求解時的精度預測。

第二部份我們探討兩項流中分子兩端極性不同的界面活性劑，這些分子通常

喜好駐於兩種液體的界面上，而且能透過吸收與釋放等過程跟可溶於液體中之界

面活性劑交流。這類問題牽涉到在可變形界面上或是複雜區域內求解偏微分方

程，因此在界面變動時如何精確計算界面與外在區域耦合之對流擴散方程實為本

問題重點所在。我們首先改寫可溶的複雜區域內界面活性劑濃度方程，透過前述

指示函數讓該方程鑲嵌於規則空間以方便計算，此外界面與外在區域之間的界面

活性劑交流，例如吸收與釋放等過程，則可視為界面的奇異項導入外在區域之濃

度方程中。在沉浸邊界法的模型之下，我們發展守衡數值格式求解界面與外在區

域耦合之濃度方程，在數值計算下依舊保持界面活性劑之總質量守衡。我們做了

一系列的數值測試來驗證我們提出的數值方法的正確性。我們也將過去針對不可

溶性界面活性劑的研究工作拓展到可溶性界面活性劑，並且將探討界面活性劑的

可溶性對於液體界面的形變造成的影響。

An immersed boundary method for simulating the interfacial flows

with soluble surfactant

Student：

Kuan-Yu Chen

Advisors：Dr.

Ming-Chih Lai

Department of Applied Mathematics

National Chiao Tung University

ABSTRACT

In the first part of this thesis, we provide a simplified one-dimensional analysis and

two-dimensional numerical experiments to predict that the overall accuracy for the

pressure problem or indicator function in immersed boundary calculations

is

first-order accurate in L

1

norm, half-order accurate in L

2

norm, but has O(1) error in

L

∞

norm. We also discuss the accuracy for another type of source terms for solving

Poisson problems with singular conditions on the interface.

In the second part, we consider the surfactant, which is an amphiphilic molecular,

under multi-phase fluids. These particles usually favor the presence in the fluid

interface, and they may couple with the surfactant soluble in one of bulk domains

through adsorption and desorption processes. This type of problem needs to solve

partial differential equations in deformable interfaces or complex domains. Thus, it is

important to accurately solve coupled surface-bulk convection-diffusion equations

especially when the interface is moving. We first rewrite the original bulk

concentration equation in an irregular domain (soluble region) into a regular

computational domain via the usage of the indicator function, which is described in

previous part, so that the concentration flux across the interface due to adsorption and

desorption processes can be termed as a singular source in the modified equation.

Based on the immersed boundary formulation, we then develop a new conservative

scheme for solving this coupled surface-bulk concentration equations which the total

誌

謝

本論文得以完成，首先要特別感謝本人的指導教授 賴明治老師在六年多的博士班生涯(以及之前兩年的碩 士班學業)之中，不管是在學業上、研究上或是生活上，都不斷給予在下各種幫助與指點。感謝老師多次安排到 國外（美國、香港等）參與研討會的機會，讓我能夠接觸到國外最新的研究發展，並且跟國外的學者溝通交流， 對研究工作與自身的語言能力都有不小的助益。另外也感謝老師經常帶領（或指派）我與其他師門學生參加國內 舉辦之各種不同的研討會，尤其是不少國際性的大型研討會，如東亞區工業與應用數學學會年會（EASIAM）等等， 其他像是國家理論科學中心的演講，以及物理、化學、機械、大氣等領域的活動，透過大量與不同學門的學者交 流，除了增廣見聞外，也能開拓研究工作在不同領域的應用。另外要感謝口試委員林文偉講座教授、楊肅煜教授、 吳金典教授與黃聰明教授的指點，使本論文更臻完善。 在博士班求學生涯中，感謝許多交大應數教授的指導與助理們的幫忙，如林琦焜教授對於偏微分方程理論的 指導，讓我對於流體力學的數學有更深的了解；感謝鄧君豪老師在數值方法分析上的指導，讓我在處理兩項流問 題得以找到盲點；其他還有葉立明教授、吳金典教授、白啟光教授、許元春教授、吳培元教授等等，在課業與其 他方面給我的提點。另外要感謝台大郭鴻基教授與中央(目前已轉往台大)蔡武廷教授對於流體力學在其他領域的 介紹與應用給我不同的體驗。由於在研究過程中需要大量運用電腦計算，感謝系上的資訊助理詹宗智先生，幫助 我處理電腦設備更新維護與其他軟硬體資源的採買，並且在各類設施操作與軟體使用上給予建議。感謝謝天長博 士在研究後期對於理論上的問題給予的協助。 沒有師門學生們互相幫忙，很難度過這漫長的博士生涯，感謝以下同門師兄與師弟妹：曾昱豪、曾孝捷、黃 仲尹、謝先皓、蔡修齊、郭乂維、胡偉帆、曾鈺傑、陳建明、張永潔、許哲維、陳昱丞、黃義閔、楊承翰、吳聲 華、林詩婷、成澤仕軒、關湘源、張毓倫、郭柏均等等。另外要感謝其他系上的學長姐、同學與學弟妹們，包括 吳恭儉、呂明杰、廖康伶、李信儀、郭君逸、黃韋強、方怡中、陳德軒、龔柏任、賴建綸、黃俊銘、蘇偉碩、邱 鈺傑等等，感謝他們在課業上與生活上的幫助。其他還有大學認識的同學與朋友們，包括黃致維、李仁喆、劉名 謙、胡仲軒、陳華軒、黃晟志、陳奕光、吳明道，還有其他透過網路認識的朋友們（族繁不及備載）。還要感謝 從小一起長大的死黨賴運辰，是到處遊山玩水的好伙伴。 最後要特別感謝我的家人，包括我的父親、母親、弟弟、祖母、姑姑與姑丈、大伯父、叔叔、姨媽、舅舅等 等，感謝他們始終不斷給我支持與鼓勵，以及無私的付出與關懷，讓我能夠專心在研究工作上持續邁進，不必為 了生活瑣事煩心。

目

錄

中文提要

………

i

英文提要

………

ii

誌謝

………

iii

目錄

………

iv

表目錄

………

vi

圖目錄

………

viii

符號說明

………

x

Chapter 1

Preliminaries………

1

Part

1

Numerical research on Poisson problems with interfacial

source terms………

3

Chapter 2

Introduction 1………

4

Chapter 3

One-dimensional analysis………

7

3.1

One-dimensional analysis for source terms as derivatives of

delta functions………

7

3.2

One-dimensional analysis for source terms as delta

functions………

11

Chapter 4

Numerical results 1………

15

4.1

Convergent test for source terms as derivatives of delta

functions………

15

4.1.1 One-dimensional problem 1………

16

4.1.2 Two-dimensional problems 1………

18

4.2

Convergent test for source terms as delta functions…………

27

4.2.1 One-dimensional problem 2………

28

4.2.2 Two-dimensional problems 2………

30

4.3

Applications to indicator functions and pressure………

37

4.3.1 One-dimensional problem 3………

38

4.3.2 Two-dimensional problems 3………

40

Part 2

Computational methods on interfacial flows with soluble

surfactant………

48

Chapter 5

Introduction 2………

49

Chapter 6

A coupled surface-bulk concentration model………

53

6.1

An embedding bulk concentration equation in a regular Cartesian

domain

………

57

Chapter 7

A conservative scheme for solving the coupled surface-bulk

concentration equations………

8.2

Algorithm for Navier-Stokes flow with soluble surfactant…

70

Chapter 9

Numerical results 2………

72

9.1

Bulk diffusion with a fixed interface………

72

9.2

Bulk convection-diffusion with a moving interface…………

76

9.3

Surface-bulk coupling with a moving interface………

79

9.4

A freely oscillating drop………

82

9.5

A drop under shear flow………

86

Chapter 10

Conclusion and future work………

92

表

目

錄

4.1

Order of accuracy for one-dimensional test withδ

hcos

………

16

4.2

Order of accuracy for one-dimensional test withδ

h √

………

18

4.3

Convergent test usingδ

h cos

for indicator function case 1 : a

circle………

19

4.4

Convergent test usingδ

h cos

for indicator function case 1 : an

ellipse………

20

4.5

Convergent test usingδ

hcos

for indicator function case 1 : a

simple closed curve………

20

4.6

Convergent test usingδ

h √

for indicator function case 1 : a

circle………

20

4.7

Convergent test usingδ

h √

for indicator function case 1 : an

ellipse………

21

4.8

Convergent test usingδ

h √

for indicator function case 1 : a

simple closed curve………

23

4.9

Convergent test usingδ

h cos

for Example 4.1.2.2………

23

4.10

Convergent test usingδ

h √

for Example 4.1.2.2………

25

4.11

Convergent test usingδ

h cos

for Example 4.1.2.3………

27

4.12

Convergent test usingδ

h √

for Example 4.1.2.3………

27

4.13

Order of accuracy for one-dimensional test withδ

h cos

_………

28

4.14

Order of accuracy for one-dimensional test withδ

h √

………

30

4.15

Convergent test usingδ

hcos

for Example 4.2.2.1………

31

4.16

Convergent test usingδ

h √

for Example 4.2.2.1………

31

4.17

Convergent test usingδ

h cos

for Example 4.2.2.2………

33

4.18

Convergent test usingδ

h √

for Example 4.2.2.2………

35

4.19

Convergent test usingδ

h cos

for Example 4.2.2.3………

37

4.20

Convergent test usingδ

h √

for Example 4.2.2.3………

37

4.21

Order of accuracy for one-dimensional test withδ

h √

………

39

4.22

Convergent test for Example 4.3.2.1………

42

4.23

Convergent test of u for Example 4.3.2.2………

47

4.24

Convergent test of v for Example 4.3.2.2………

47

4.25

Convergent test of ▽．u

h

for Example 4.3.2.2………

47

9.1

The L

2

and L

1

errors and their convergent rates for the bulk

diffusion with a fixed interface at T=0.5………

74

9.2

The convergent study of numerical leakage relative error

76

9.4

The L

2

errors and their convergent rates for the bulk and

surface concentrations at T = 0.5………

79

9.5

The L

2

errors and their convergent rates for the bulk and

surface surfactant concentrations, and the fluid velocity field at

T = 0.5………

83

9.6

The L

2

errors and their convergent rates for the bulk and

surface surfactant concentrations, and the fluid velocity field at

T = 0.5………

圖

目

錄

4.1

Comparison between numerical and analytic solutions for 1D

case………

17

4.2

Numerical solution for indicator function case 3 : a simple

closed curve………

21

4.3

Comparison between numerical and analytic solutions for

indicator function case 1 : a circle………

22

4.4

Comparison between numerical and analytic solutions for

Example 4.1.2.2………

24

4.5

Comparison between numerical and analytic solutions for

Example 4.1.2.3………

26

4.6

Comparison between numerical and analytic solutions for 1D

case………

29

4.7

Comparison between numerical and analytic solutions for 2D

case 4.2.2.1………

32

4.8

Comparison between numerical and analytic solutions for 2D

case 4.2.2.2………

34

4.9

Comparison between numerical and analytic solutions for 2D

case 4.2.2.3………

36

4.10

Comparison between numerical and analytic solutions for 1D

case………

40

4.11

Comparison between numerical and analytic solutions for 2D

case 4.3.2.1………

42

4.12

Comparison between numerical and analytic solutions for u in

case 4.3.2.2………

44

4.13

Comparison between numerical and analytic solutions for v in

case 4.3.2.2………

45

4.14

Divergence of velocity field of numerical solution in case

4.3.2.2………

46

5.1

Schematic diagram of surfactant particles with hydrophilic

heads and hydrophobic tails, cited from

”http://en.wikipedia.org/wiki/Surfactant”………

50

6.1

Illustration of domains………

54

7.1

The computational domain Ω using staggered grid with mesh

size h………

60

7.2

Illustration of Lagrangian markers………

60

9.1

The bulk concentration along the horizontal line y = 0 at

75

error.

Lower panel: Time evolutionary plot of leaking mass relative

error inside the interface………

9.3

Upper panel: Time evolutionary plot of total mass relative

error.

Lower panel: Time evolutionary plot of leaking mass relative

error inside the interface………

77

9.4

The bulk concentration at different times. Left column: The

bulk concentration contour plots and the interface positions.

Right column: The bulk concentration plots along the

horizontal line y = 0. The dashed line in the first right plot

denotes the initial bulk concentration along y = 0………

78

9.5

Upper panel: Time evolutionary plot of total mass relative

error.

Lower panel: Time evolutionary plot of leaking mass relative

error inside the interface………

80

9.6

Upper panel: the bulk concentration contour plots and the

interface positions.

Lower left: the bulk concentration plots along the horizontal

line y = 0. The dashed line in the first plot denotes the initial

bulk concentration.

Lower right: the surface concentration along the interface in

counter-clockwise way starting from the point marked

by ”o”………

81

9.7

Upper panel: Time evolutionary plot of total mass relative

error.

Lower panel: Time evolutionary plot of leaking mass relative

error inside the interface………

84

9.8

The comparison of insoluble (denoted by ”-.”) and soluble

(denoted by ”-”) interfacial flows for a freely oscillating drop

85

9.9

The comparison of clean (denoted by ”-.”) and soluble

(denoted by ”-”) interfacial flows for a drop under shear flow

87

9.10

Left: the bulk concentration plots along the horizontal line y =

0. The dashed line in the first plot denotes the initial bulk

concentration.

Right: the surface concentration along the interface in

counter-clockwise way starting from the point marked by ”o”

88

9.11

The comparison of aspect ratio between clean (denoted by ”-.”)

and soluble (denoted by ”-”) cases………

89

9.12

Upper panel: Time evolutionary plot of total mass relative

error.

Lower panel: Time evolutionary plot of leaking mass relative

error inside the interface………

符 號 說 明

Ω

：computational domain

Σ

：interface immersed in Ω

Ω

0

：inner region enclosed by Σ in Ω

Ω

1

：outer region of Ω

Δx ：discrete mesh width in x-direction

Δy ：discrete mesh width in y-direction

Δt ：discrete time step size

h

：unit mesh width

α

：surface parameterization

dα ：surface element under parameterization

Δα

：discrete surface mesh width

X(α) ：parameterized position of the interface



 

X

：stretching factor

τ ：unit tangent vector of interface

n

：unit outward normal vector of interface

C

：bulk surfactant concentration

Γ ：surface surfactant concentration

p

：pressure

u

：velocity field

u

：velocity component in x-direction

v

：velocity component in y-direction

U

：velocity field on the interface

f

：external force

F

：interfacial force

G(x;y)

：1D Green function

I(x)

：indicator function which indicates inner region

H(x)

：indicator function which indicates outer region

▽

2

：Laplacian ( orΔ)

▽

s

：surface gradient

▽

s

． ：surface divergence

▽

s2

：surface Laplacian ( orΔ

s

)

δ(x) ：1D delta function

δ

2

(x) ：2D delta function

δ

h

(x) ：1D discrete delta function

δ

h 2

(x) ：2D discrete delta function

δ

hcos

：discrete delta function with cosine type

δ

h

√

：discrete delta function with square root type

D

α

：surface first derivative difference operator

D

h

：1D first derivative centered difference operator

D

h2

：1D second derivative centered difference operator

▽

h

：2D gradient difference operator

▽

h

． ：2D divergence difference operator

▽

h 2

：2D Laplacian difference operator

|| ||

1

：(discrete) one norm

|| ||

2

：(discrete) two norm

|| ||

∞

：(discrete) infinity norm

μ

：viscosity

Re

：Reynolds’ number

Ca

：Capillary number

Cs

：bulk surfactant concentration adjacent to the interface

Pe

：Peclet number

Pe

s

：surface Peclet number

S

a

：absorption Stanton number

S

d

：desorption Stanton number

λ

：dimensionless absorption depth

σ

：surface tension

Chapter 1

Preliminaries

This thesis is a combination of our previous works. The first part is based on the article [12] ”K.-Y. Chen, K.-A. Feng, Y. Kim, M.-C. Lai, A note on pressure accuracy in immersed boundary method for Stokes flow, Journal of Computational Physics, 230 (2011), 4377– 4383.” We provide a simplified one-dimensional analysis and two-dimensional numerical experiments to predict that the overall accuracy for the pressure problem or indicator

function in immersed boundary calculations is first-order accurate in L1 norm, half-order

accurate in L2 norm, but has O(1) error in L∞ norm. We also discuss the accuracy for

another type of source terms for solving Poisson problems with singular conditions on the interface. In this case, we prove that the convergent rate is second-order accurate in

L1 norm, one and half-order accurate in L2 norm, and first-order accurate in L∞ norm.

Moreover, we will give some applications to solve second-order elliptic equations with piecewise-constant coefficients by indicator function, and compute the velocity in Stokes equations by using the solution of pressure equations we obtained.

The following part provides the details of the paper [13] ”K.-Y. Chen, M.-C. Lai, A con-servative scheme for solving coupled surface-bulk convection-diffusion equations with an application to interfacial flows with soluble surfactant, Journal of Computational Physics, 257(2014), 1–18.” We consider the surfactant, which is an amphiphilic molecular, under

and desorption processes. This type of problem needs to solve partial differential equa-tions in deformable interfaces or complex domains. Thus, it is important to accurately solve coupled surface-bulk convection-diffusion equations especially when the interface is moving. We first rewrite the original bulk concentration equation in an irregular do-main (soluble region) into a regular computational dodo-main via the usage of the indicator function, which is described in previous part, so that the concentration flux across the interface due to adsorption and desorption processes can be termed as a singular source in the modified equation. Based on the immersed boundary formulation, we then develop a new conservative scheme for solving this coupled surface-bulk concentration equations which the total surfactant mass is conserved in discrete sense. A series of numerical tests has been conducted to validate the present scheme. As an application, we extend our previous work [24] ”M.-C. Lai, Y.-H. Tseng and H. Huang, An immersed boundary method for interfacial flows with insoluble surfactant, Journal of Computational Physics, 227 (2008) 7279-7293” to the soluble case. The effects of solubility of surfactant on drop deformations in a quiescent and shear flow are investigated in detail.

Part I

Numerical research on Poisson

problems with interfacial source

Chapter 2

Introduction 1

In this part, we consider the problems of solving a Poisson equation in the computational domain Ω (either in one-dimension or in two-dimension) with a source term defined only on a boundary (one point for one-dimensional case, and one-dimensional interface for two-dimensional case, respectively) Σ immersed in Ω. This type of the problems arise from using Immersed Boundary Method to solve the stationary Stokes flow defined on irregular domain or containing interfacial singularities inside the regular domain. The Stokes problem in the immersed boundary formulation is defined as

−∇p + µ∆u + Z

Σ

F(s) δ2_{(x − X(s)) ds = 0,} (1)

∇ · u = 0. (2)

The two-dimensional Dirac delta function is defined as

δ2(x) = δ(x) δ(y), (3)

which is the combination of two one-dimensional Dirac delta functions. Since the im-mersed boundary force F is only exerted along the interface Σ, we use the integral with the Dirac delta function to keep the formulation is defined along the interface Σ. There-fore, the above immersed boundary formulation is a typical singular problem with a delta function source. To solve this problem, a simple ansatz is by taking the divergence oper-ator on Eq. (1), which means that

−∆p(x) + µ∇ · ∆u + ∇ · Z

Σ

Use the incompressibility constraint in Eq. (2) to take out u, we obtain the pressure equation ∆p(x) = ∇ · Z Σ F(s) δ2_{(x − X(s)) ds. (5)}

Notice that, the above equation is what we mention before, which involves solving Poisson equation with a source term which can be written as the derivatives of Dirac delta function. Furthermore, in the general immersed boundary computations, one often uses periodic or Neumann boundary conditions to solve pressure equation. For simplicity, throughout this work, we just use the Dirichlet boundary condition since we are more concerned about the accuracy caused by the derivatives of Dirac delta function near the interface.

After solving the pressure p, we have to solve the velocity field u by solving the equation (1) with substituting p. In this case, we need to solve another two Poisson equations with source terms on the interface, which can be described as Dirac delta function. Therefore, we deal with two kinds of source terms to the Poisson equations introduced by the Stokes flow problem.

Another example comes from two-phase flow problem. In the former work by Tryg-gvason et. al. [37], they introduced the indicator function in order to track the regions of two-phase flow. For any quantity q, which could be density, viscosity or others, is valued piecewise constant in the domain, i.e. it is discontinuous across the interface. It can be represented by the following:

q(x) = qout + (qin− qout)I(x), (6)

where qin and qout are the constant quantity inside and outside of the interface,

respec-tively. Note that, the indicator function has the value one (I = 1) inside the immersed boundary Γ and the value zero (I = 0) outside. Assume we use the immersed boundary

Σ to divide the domain Ω into two parts: Ω0, which is inside Σ; and Ω1, which is outside

Σ. The indicator function can be written as 

The indicator function can be calculated as the following procedure [37]. Let n be the unit outward normal to the interface, then the indicator function can be represented by

I(x) = Z

Ω0

δ2_{(x − ˜x) d˜x.}

By taking the gradient and then the divergence operators, we have ∇I(x) = − Z Σ nδ2_{(x − X(s)) ds,} ∆I(x) = −∇ · Z Σ nδ2_{(x − X(s)) ds.} (8)

Thus, the indicator function can be obtained by solving a similar equation as Eq. (5) with the special singular forcing term F(s) = −n(s).

In the following of this part, our goal is to use the standard finite difference scheme with smoothing discrete delta function to discretize the equations (5) and (8), and to investigate the numerical accuracy. Moreover, we will also discuss other types of singularity on the interface to the Poisson equations, and do further analysis and validations to these problems.

Chapter 3

One-dimensional analysis

In this chapter, we will discuss the accuracy of numerical solution by using simple finite-difference schemes to our problems.

3.1

One-dimensional analysis for source terms as

deriva-tives of delta functions

We consider the one-dimensional equation as

d2_u

dx2 = c

d

dxδ(x − α), 0 ≤ x ≤ 1 (9)

with boundary condition

u(0) = u(1) = 0, (10)

and the interface is located at x = α ∈ (0, 1). The exact solution of Eq. (9) can be expressed as u(x) = Z 1 0 G(x; y) c d dyδ(y − α) dy (11)

where G(x; y) is the well-known Green’s function, which solves

d2G

The Green’s function G(x; y) can be explicit written as G(x; y) =      x(y − 1), 0 ≤ x ≤ y, y(x − 1), y < x ≤ 1 . (13)

Without losing the generality, we set c = 1. By applying the integration by parts into Eq. (11) , we obtain u(x) = Z 1 0 G(x; y) c d dyδ(y − α) dy = G(x; y)δ(y − α)|y=10 − Z 1 0 d dyG(x; y) δ(y − α) dy = − Z 1 0 d dyG(x; y) δ(y − α) dy. (14)

The Boundary terms vanish by using the definition of Green’s function.

Based on uniform grid with grid points xj = jh, j = 0, 1, ..., N where h = 1/N, we use

the standard centered difference scheme to discretize Eq. (9) with c = 1 as following

U_{j−1− 2U}j + Uj+1

h2 =

δh(xj+1− α) − δh(xj−1− α)

2h . (15)

The discrete delta function in [31] defined as

δh(x)cos =      1 4h(1 + cos( πx 2h)), if |x| ≤ 2h, 0, otherwise . (16)

Although there are more different discrete delta functions can be found in [4, 44] and other papers, the usage of different delta functions cannot lead to different conclusions that will be given in this subsection.

The discrete delta functions above satisfies N

X

m=0

δh(xm− α)h = 1, (17)

which is the corresponding basic requirement for discrete delta function. For simplicity,

we denote the first-order and second-order centered difference operators as Dh and D2h,

respectively. Analog to analytic solution in (11), the discrete solution Uj of Eq. (15) can

also be written as Uj = h N X m=0 GjmDhδh(xm− α) (18)

where Gjm = G(xj; xm) is the discrete version of Green’s function defined as Gjm =      xj(xm− 1), 0 ≤ j ≤ m, xm(xj − 1), m < j ≤ N . (19)

We can immediately check that G satisfies D2

hGjm = _h1δjm where δjm is the Kronecker

delta function.

Now, by taking summation by parts and the property of discrete delta function, we

can rewrite the numerical solution Uj as

Uj = h N X m=0 GjmDhδh(xm− α) = h N −1 X m=1 Gjm δh(xm+1− α) − δh(xm−1− α) 2h = h N X m=2 G_j(m−1)δh(xm− α) 2h − h N −2 X m=0 Gj(m+1) δh(xm− α) 2h = −h N −1 X m=1 Gj(m+1)− Gj(m−1) 2h δh(xm− α) = −h N X m=0 DhGjmδh(xm− α) (20)

Hence the point-wise error between Uj and u(xj) can be expressed as

|Uj − u(xj)| =      h N X m=0 DhGjmδh(xm− α) − Z 1 0 d dyG(xj; y) δ(y − α) dy      ≤      h N X m=0 DhGjmδh(xm− α) − h N X m=0 d dyG(xj; y)     y=xm δh(xm− α)      +      h N X m=0 d dyG(xj; y)     y=xm δh(xm− α) − d dyG(xj; y)     y=α      = E1+ E2

First, E1 is the error from discretizing the differentiation of Green’s function. Using

the fact that the derivative of Greens function is

d 

 

and that its discrete counterpart is DhGjm =            xj, j < m xj −1₂, j = m xj − 1, j > m , (22) we obtain DhGjm− d dyG(xj; y)     y=xm =      1 2, if m = j 0, otherwise (23) By using (23), we compute E1 =      h N X m=0 DhGjmδh(xm− α) − h N X m=0 d dyG(xj; y)     y=xm δh(xm− α)      =      h N X m=0 DhGjm− d dyG(xj; y)     y=xm ! δh(xm− α)      =     h1 2δh(xj − α)     =      O(1), when |xj− α| ≤ 2h 0, otherwise (24)

The second part of the error E2 is simply an interpolating error for the function

d dyG(xj; y)    y=xm

. Using the formula in Eq. (21) and the first moment condition in (17), since the discrete delta function has finite support 4h, we can obtain

E2 =      h N X m=0 d dyG(xj; y)     y=xm δh(xm− α) − d dyG(xj; y)     y=α      =      h j−1 X m=0 d dyG(xj; y)     y=xm δh(xm− α) + h N X m=j d dyG(xj; y)     y=xm δh(xm− α) − d dyG(xj; y)     y=α      =      h j−1 X m=0 (xj− 1)δh(xm− α) + h N X m=j xjδh(xm− α) − d dyG(xj; y)     y=α     

For xj ≤ α, _dydG(xj; y)    y=α = xj, then E2 =      h j−1 X m=0 (xj − 1)δh(xm− α) + h N X m=j xjδh(xm− α) − xj      =      h j−1 X m=0 δh(xm− α)      =      O(1), as α − xj ≤ 2h 0, otherwise (25) Similarly, for xj > α, _dydG(xj; y)    y=α = xj − 1, then E2 =      h j−1 X m=0 (xj− 1)δh(xm− α) + h N X m=j xjδh(xm − α) − (xj − 1)      =      h N X m=j δh(xm− α)      =      O(1), as xj− α ≤ 2h 0, otherwise (26)

Therefore, we combine (25) and (26) to get

E2 =      O(1), as |xj − α| ≤ 2h 0, otherwise (27) From the above analysis, one can immediately see that the point-wise error appears only at some points around the singular point α, which means that the maximum

er-ror kuh− uk_∞ is of order O(1). For the same reason, we can conclude that L1(kuh− uk₁)

and L2(kuh− uk₂) errors are of order O(h) and O(h1/2), respectively. Our numerical

results in final section will confirm this conclusion.

3.2

One-dimensional analysis for source terms as delta

with the same boundary condition as Eq. (10), and the interface is located at x = α ∈ (0, 1). The exact solution of (28) can be expressed as

u(x) =

Z 1

0 G(x; y) cδ(y − α) dy

(29) we use the standard centered difference scheme to discretize Eq. (28) with c = 1 as following

U_{j−1− 2U}j + Uj+1

h2 = δh(xj − α), (30)

Analog to analytic solution in Eq. (29), the discrete solution Uj of Eq. (30) can also be

written as Uj = h N X m=0 Gjmδh(xm− α) (31)

Hence the point-wise error between Uj and u(xj) can be expressed as

|Uj − u(xj)| =      h N X m=0 Gjmδh(xm− α) − Z 1 0 G(x; y)δ(y − α) dy      (32) According to [38], [44] and [4], we consider two different discrete delta functions

δ_hcos=      1 4h(1 + cos( πx 2h)), if |x| ≤ 2h, 0, otherwise ; (33) δ_h√ =              1 8  3 − 2x h  + q 1 + 4x h  − 4  x h   2 if |x| ≤ h, 1 8  5 − 2x h  − q −7 + 12x h  − 4  x h   2 if h < |x| ≤ 2h, 0, otherwise . (34)

with the relative moment conditions h N X m=0 (xm− α)δhcos(xm− α) = O(h) (35) h N X m=0 (xm− α)δ √ h (xm− α) = 0 (36)

Suppose α lies in some interval (xi, xi+1). For xj < α − h, we have |Uj − u(xj)| =      h i+2 X m=i−1 Gjmδh(xm− α) − xj(α − 1)      =      h i+2 X m=i−1 xj(xm− 1)δh(xm− α) − xj(α − 1)      =      h i+2 X m=i−1 xj(xm− α + α − 1)δh(xm− α) − xj(α − 1)      =      h i+2 X m=i−1 xj(xm− α)δh(xm− α)      =     

O(h), for δcos

h 0, for δ_h√ (37) For α − h ≤ xj ≤ α, we compute |Uj − u(xj)| =      h i+2 X m=i−1 Gjmδh(xm− α) − xj(α − 1)      =      h x_i−1(xj − 1)δh(xi−1− α) + i+2 X m=i xj(xm− 1)δh(xm− α) ! − xj(α − 1)      = |h [(xj(xi−1− α + α − 1) + h) δh(xi−1− α) + i+2 X m=i xj(xm− α + α − 1)δh(xm− α) # − xj(α − 1)      =      h2δh(xi−1− α) + h i+2 X m=i−1 xj(xm− α)δh(xm− α)      = O(h). (38) For α = xj, we derive |Uj− u(xj)| =      h j+1 X m=j−1 Gjmδh(xm− α) − xj(xj − 1)      =     x_j−1(xj − 1) 4 + xj(xj − 1) 2 + xj(xj+1− 1) 4 − xj(xj − 1)     = h (39)

For α < xj ≤ α + h, we calculate |Uj − u(xj)| =      h i+2 X m=i−1 Gjmδh(xm− α) − α(xj− 1)      =      h xj(xi+2− 1)δh(xi+2− α) + i+1 X m=i−1 xm(xj − 1)δh(xm− α) ! − α(xj − 1)      = |h [((xi+2− α + α)(xj − 1) + h) δh(xi+2− α) + i+2 X m=i (xm− α + α)(xj− 1)δh(xm− α) # − α(xj− 1)      =      h2δh(xi+2− α) + h i+2 X m=i−1 (xj − 1)(xm− α)δh(xm− α)      = O(h). (40) For xj > α + h, we obtain |Uj− u(xj)| =      h i+2 X m=i−1 Gjmδh(xm− α) − α(xj − 1)      =      h i+2 X m=i−1 xm(xj− 1)δh(xm − α) − α(xj− 1)      =      h i+2 X m=i−1 (xm− α + α)(xj − 1)δh(xm− α) − α(xj − 1)      =      h i+2 X m=i−1 (xj− 1)(xm− α)δh(xm− α)      =     

O(h), for δcos

h

0, for δ√_h

(41)

Hence we combine Eq. (37) to Eq. (41) to get

for δ_hcos_{, |U}j − u(xj)| = O(h) ∀j (42)

for δ√_{h , |U}j − u(xj)| =      O(h), as |xj− α| ≤ h 0, otherwise (43)

From the above analysis, one can immediately see that, under proper moment condi-tions satisfied, the point-wise error appears only at some points around the singular point

α, which means that the maximum error kuh− uk_∞ is of order O(h). For the same

rea-son, we can conclude that L1(kuh− uk₁) and L2(kuh− uk₂) errors are of order O(h2) and

Chapter 4

Numerical results 1

Throughout this chapter, we validate our derivations of convergence analysis via a series of examples. For different kinds of source term, we check one-dimensional cases at first. We also test on two-dimensional examples to see if our analysis is hold numerically.

4.1

Convergent test for source terms as derivatives

of delta functions

In this section, we will show the convergent results of solving Poisson problems with source terms as derivatives of delta functions. Since the proof we derived is valid for any kinds

of discrete delta functions, we pick up δcos

h and δ

√

h from [4, 44] for comparison. The ratio

in all tables of this section means the convergent orders of accuracy, which is computed by

log( ku − uNk

ku − u2Nk

)/ log( 1/N

1/2N), (44)

4.1.1

One-dimensional problem 1

In this subsection, we consider the following one-dimensional problem to verify the proof in section 3.1. The equation is

d2_u

dx2 = c

d

dxδ(x − α) + g, 0 < x < 1. (45)

Here, the interface is set at the point α = π/6. The exact solution is given as

u(x) =      x3_{+ 2αx}2 _{if x ≤ α} 7 3(x 3_{− 1) if x > α} (46)

where the jump of solution u at the interface c is equal to −(2α3 + 7)/3. The regular

source term g can be computed by the analytic solution, which is written as

g(x) =      6x + 4α if x ≤ α 14x if x > α . (47)

Table 4.1 shows the order of accuracy for our test with using δcos

h which confirms our

one-dimensional analysis in previous section. In Table 4.2 we use δ√_h to verify the conclusion

we mention. We could observe the same behavior of the errors even if we change the discrete delta function. Figure 4.1 shows the maximum error always occurs across the interface. Even if we refine the mesh, there still exists an O(1) error, which matches our derivations of proof in previous section.

mesh _{ku − U}hk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32 8.8827E-01 – 1.7057E-01 – 4.2479E-02 –

64 6.1709E-01 0.5255 1.0736E-01 0.6678 1.9004E-02 1.1604

128 1.1847E-00 -.9409 1.0708E-01 0.0038 1.2980E-02 0.5500

256 1.1579E-00 0.0329 7.4120E-02 0.5307 6.3791E-03 1.0248

512 1.1030E-00 0.0700 5.0171E-02 0.5630 3.0741E-03 1.0531

Table 4.1: Order of accuracy for one-dimensional test with δcos

0 0.2 0.4 0.6 0.8 1 −2 −1.5 −1 −0.5 0 0.5 1D case comparison exact sol numer sol, N=32 numer sol, N=128 0.5 0.51 0.52 0.53 0.54 0.55 −2 −1.5 −1 −0.5 0 0.5 exact sol numer sol, N=512

mesh _{ku − U}hk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32 9.1311E-01 – 1.7356E-01 – 4.2955E-02 –

64 6.1970E-01 0.5592 1.0738E-01 0.6927 1.9005E-02 1.1764

128 1.1875E-00 -.9382 1.0731E-01 0.0008 1.3001E-02 0.5477

256 1.1634E-00 0.0295 7.4446E-02 0.5275 6.3993E-03 1.0226

512 1.1138E-00 0.0628 5.0603E-02 0.5569 3.0926E-03 1.0490

Table 4.2: Order of accuracy for one-dimensional test with δ_h√ .

4.1.2

Two-dimensional problems 1

For two-dimensional problem, we generally write the equation in the form : ∆u = ∇ ·

Z

Σ

F(s) δ2_{(x − X(s)) ds + g} (48)

in a rectangular domain Ω = [a, b] × [c, d] with given Dirichlet boundary conditions. We use an M × N uniform grid with mesh width ∆x = ∆y = h to divide the domain Ω. The

notation Ui,j represents the discrete solution at (xi, yj), where xi = ih, h = 0, ..., M and

yj = jh, j = 0, ..., N. We also use a group of marker points X(sk) = (Xk, Yk) along the

interface Σ with the mesh points sk = k∆s. In this case, the mesh width ∆s is about a

half of h. Then we use regular centered difference scheme to discrete Eq. (48) as Ui+1,j − 2Ui,j + Ui−1,j

h2 +

Ui,j+1− 2Ui,j + Ui,j−1

h2 = Fi+1/2,j − Fi−1/2,j h + Fi,j+1/2− Fi,j−1/2 h + gi,j (49) where Fi+1/2,j = X k F (sk) δh(xi+ h/2 − X(sk))δh(yj − Y (sk)) ∆s, (50) F_i−1/2,j =X k F (sk) δh(xi− h/2 − X(sk))δh(yj − Y (sk)) ∆s, (51) Fi,j+1/2= X k F (sk) δh(xi− X(sk))δh(yj + 1/2 − Y (sk)) ∆s, (52) and F_i,j−1/2=X k F (sk) δh(xi − X(sk))δh(yj − 1/2 − Y (sk)) ∆s. (53)

The resultant matrix equation can be solved efficiently by the fast direct solver in Fishpack [3].

Example 4.1.2.1: In the first example, we test the accuracy of the indicator function

which is described in Eq. (8). For completeness, we test three different interface Σ in the domain [−1, 1] × [−1, 1] as follows.

1. Σ is a circle with the radius 0.3, which is centered at (0, 0).

2. Σ is an ellipse with the major radius 0.9 and minor radius 0.1, which is centered at (0, 0).

3. Σ is a simple closed curve written in polar coordinates : r = 0.5 + 0.25 cos(5θ).

Table 4.3-4.5 show the convergence tests for those three different cases with using δcos

h ,

while Table 4.6-4.8 are the correspondent results with using δ_h√ . The convergent results

of the indicator function calculation are strongly supporting our conclusion, which shows

first-order convergence in L1 norm and half-order convergence in L2 norm, although there

is an O(1) error in maximum norm. The results are consistent with one-dimensional analysis. Figure 4.5 shows the cross section of the numerical and analytic solutions along the line y = 0. It implies that there exists an O(1) error at the transition area of the indicator function.

M × N kI − Ihk∞ ratio kI − Ihk2 ratio kI − Ihk1 ratio

32×32 3.6463E-01 – 1.3162E-01 – 6.3848E-02 –

64×64 4.5555E-01 -.3212 9.5529E-02 0.4622 3.2182E-02 0.9883

128×128 4.8736E-01 -.0974 7.2764E-02 0.3926 1.6837E-02 0.9345 256×256 4.8610E-01 0.0037 4.9738E-02 0.5488 8.2361E-03 1.0316 512×512 4.9805E-01 -.0620 3.4744E-02 0.5175 4.0955E-03 1.0078

Table 4.3: Convergent test using δcos

M × N kI − Ihk∞ ratio kI − Ihk2 ratio kI − Ihk1 ratio

32×32 6.7302E-01 – 1.9248E-01 – 1.2276E-01 –

64×64 5.0021E-01 0.4280 1.4391E-01 0.4194 6.5139E-02 0.9141

128×128 4.9922E-01 0.0027 9.7919E-02 0.5555 3.1954E-02 1.0275 256×256 4.9834E-01 0.0024 6.8903E-02 0.5070 1.5951E-02 1.0023 512×512 4.9617E-01 0.0061 4.8685E-02 0.5010 7.9510E-03 1.0043

Table 4.4: Convergent test using δcos

h for indicator function case 2 : an ellipse.

M × N kI − Ihk∞ ratio kI − Ihk2 ratio kI − Ihk1 ratio

32×32 5.9986E-01 – 2.5162E-01 – 2.0860E-01 –

64×64 5.5492E-01 0.1123 1.8259E-01 0.4625 1.0827E-01 0.9460

128×128 5.3029E-01 0.0654 1.2910E-01 0.5000 5.4431E-02 0.9211 256×256 5.1669E-01 0.0374 9.0547E-02 0.5116 2.7064E-02 1.0079 512×512 5.1194E-01 0.0132 6.4251E-02 0.4948 1.3547E-02 0.9983

Table 4.5: Convergent test using δcos

h for indicator function case 3 : a simple closed curve.

M × N kI − Ihk∞ ratio kI − Ihk2 ratio kI − Ihk1 ratio

32×32 3.6989E-01 – 1.3208E-01 – 6.3755E-02 –

64×64 4.5684E-01 -.3047 9.5640E-02 0.4657 3.2222E-02 0.9844

128×128 4.8714E-01 -.0926 7.2906E-02 0.3915 1.6858E-02 0.9346 256×256 4.8617E-01 0.0028 4.9854E-02 0.5483 8.2307E-03 1.0342 512×512 4.9809E-01 -.0350 3.4826E-02 0.5175 4.0799E-03 1.0124

Table 4.6: Convergent test using δ_h√ for indicator function case 1 : a circle.

Example 4.1.2.2: In this example, we solve a pressure equation from Stokes flow

problem in [14]. Since the analytic solutions to this problem are available, Lai et el. [23] had use a simplified version of this example to test immersed interface method. The pressure equation to be solved is described in Eq. (5). The computational domain is

Figure 4.2: Numerical solution for indicator function case 3 : a simple closed curve.

M × N kI − Ihk∞ ratio kI − Ihk2 ratio kI − Ihk1 ratio

32×32 6.7492E-01 – 1.9304E-01 – 1.2199E-01 –

64×64 5.0036E-01 0.4316 1.4485E-01 0.4142 6.4800E-02 0.9126

128×128 4.9957E-01 0.0021 9.8565E-02 0.5554 3.1684E-02 1.0322 256×256 4.9850E-01 0.0030 6.9379E-02 0.5065 1.5807E-02 1.0031 512×512 4.9660E-01 0.0054 4.8969E-02 0.5025 7.8913E-03 1.0022

Table 4.7: Convergent test using δ√_h for indicator function case 2 : an ellipse.

Ω = [−2, 2] × [−2, 2], and the interface is a unit circle centered at (0, 0), i.e. X(θ) = (cos θ, sin θ). The exact solution is written in polar coordinates as

p(r, θ) =       −r3 _{sin(3θ) if r ≤ 1} r−3 _{sin(3θ) if r > 1} (54)

−1 −0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 2D case comparison 0.25 0.3 0.35 0 0.2 0.4 0.6 0.8 1 exact sol numer sol, N=32 numer sol, N=128 exact sol numer sol, N=512

Figure 4.3: Comparison between numerical and analytic solutions for indicator function case 1 : a circle.

M × N kI − Ihk∞ ratio kI − Ihk2 ratio kI − Ihk1 ratio

32×32 5.9802E-01 – 2.5211E-01 – 2.0901E-01 –

64×64 5.5251E-01 0.1141 1.8301E-01 0.4620 1.0829E-01 0.9486

128×128 5.2782E-01 0.0658 1.2940E-01 0.5000 5.4356E-02 0.9942 256×256 5.1434E-01 0.0371 9.0733E-02 0.5120 2.7025E-02 1.0081 512×512 5.1359E-01 0.0020 6.4375E-02 0.4950 1.3529E-02 0.9982

Table 4.8: Convergent test using δ_h√ for indicator function case 3 : a simple closed curve.

Table 4.9 shows the convergence tests for this case with using δcos

h , and table 4.10

presents the computation results with using δ_h√ , respectively. Figure 4.6 shows the

com-parison of numerical and analytic solution along the cross section x = 0. One can see that the maximum error of pressure occurs at the interface.

M × N kp − Phk∞ ratio kp − Phk2 ratio kp − Phk1 ratio

32×32 1.5643E-00 – 9.5960E-01 – 2.0421E-00 –

64×64 1.7182E-00 -.1354 6.9177E-01 0.4720 1.1867E-00 0.7831

128×128 1.8342E-00 -.0943 4.9447E-01 0.4843 6.4868E-01 0.8712 256×256 1.9086E-00 -.0573 3.4999E-01 0.4985 3.4044E-01 0.9300 512×512 1.9284E-00 -.0149 2.4775E-01 0.4983 1.7495E-01 0.9604

Table 4.9: Convergent test using δcos

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 2D case comparison 0.95 1 1.05 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 exact sol numer sol, N=32 numer sol, N=128 exact sol numer sol, N=512

M × N kp − Phk∞ ratio kp − Phk2 ratio kp − Phk1 ratio

32×32 1.5054E-00 – 8.4134E-01 – 1.5983E-00 –

64×64 1.6069E-00 -.1676 6.4567E-01 0.3818 9.4537E-01 0.7576

128×128 1.8046E-00 -.0938 4.7711E-01 0.4364 5.2107E-01 0.8594 256×256 1.8913E-00 -.0676 3.4350E-01 0.4740 2.7398E-01 0.9274 512×512 1.9219E-00 -.0231 2.4531E-01 0.4857 1.4085E-01 0.9599

Table 4.10: Convergent test using δ_h√ for Example 4.1.2.2.

Example 4.1.2.3: As a last example, we consider Eq. (48) in the square domain

Ω = [−1, 1] × [−1, 1] with analytic solutions. The interface Σ dividing Ω into inner

part Ω0 and outer part Ω1, which is a simple closed curve written in polar coordinates

r = 0.5 + 0.25 cos(5θ). The solution u is given by

u =        (x2_{− 1)(y}2_{− 1) + 1 if (x, y) ∈ Ω} 0 (x2_{− 1)(y}2_{− 1) if (x, y) ∈ Ω} 1 , (55)

thus the boundary force F is simply the normal vector n along the interface Σ. The

external source g can be easily computed from the solution u as g = 2(x2_{+ y}2_{− 2).}

The convergence tests based on δcos

h are shown in table 4.11, while table 4.12 lists

the corresponding order of accuracy with using δ_h√ . Figure 4.7 shows the numerical and

analytic solution along the line y = 0. It indicates that refining mesh could not improve the maximum error.

One can observe that the convergent results from example 2 and 3 are consistent with

our conclusion, i.e. first-order accurate in L1 norm, half-order accurate in L2 norm, but

−1 −0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2D case comparison −0.3 −0.28 −0.26 −0.24 −0.22 −0.2 0.8 1 1.2 1.4 1.6 1.8 2 exact sol numer sol, N=32 numer sol, N=128 exact sol numer sol, N=512

M × N ku − Uhk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32×32 9.7976E-01 – 4.6428E-01 – 3.6083E-01 –

64×64 9.8231E-01 -.0037 3.4142E-01 0.4434 1.9150E-01 0.9140

128×128 9.8129E-01 0.0015 2.4376E-01 0.4860 9.7187E-02 0.9784 256×256 9.8105E-01 0.0003 1.7146E-01 0.5075 4.8266E-02 1.0097 512×512 9.8169E-01 -.0009 1.2165E-01 0.4951 2.4217E-02 0.9950

Table 4.11: Convergent test using δcos

h for Example 4.1.2.3.

M × N ku − Uhk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32×32 9.8024E-01 – 4.6421E-01 – 3.6094E-01 –

64×64 9.8190E-01 -.0024 3.4141E-01 0.4432 1.9136E-01 0.9155

128×128 9.8141E-01 0.0007 2.4374E-01 0.4861 9.7017E-02 0.9799 256×256 9.8090E-01 0.0007 1.7145E-01 0.5075 4.8187E-02 1.0095 512×512 9.8150E-01 -.0008 1.2164E-01 0.4952 2.4179E-02 0.9948

Table 4.12: Convergent test using δ_h√ for Example 4.1.2.3.

4.2

Convergent test for source terms as delta

func-tions

In this section, we will show the convergent results of solving Poisson problems with source terms as delta functions. Since the conclusion we obtained depends on the moment

condition of discrete delta functions, we pick up δcos

h and δ

√

h from [4, 44] for comparison.

The ratio in all tables of this section is the same as previous section, which means the orders of accuracy.

4.2.1

One-dimensional problem 2

In this subsection, we consider the following one-dimensional problem to verify the proof in section 3.2. The equation is

d2_u

dx2 = cδ(x − α) + g, 0 < x < 1. (56)

Here, the interface is fixed at the point α = π/4 − 0.15. The exact solution is given as u(x) =      (x − α)3 _{if x ≤ α} (x − α)3_{+ cx + d if x > α} (57)

where the jump of derivative of solution u at the interface c is set as c = −1 and the constant d for keeping the continuity of u is d = −cα The regular source term g can be derived from the analytic solution, which is represented by g(x) = 6(x − α).

By using δcos

h , table 4.13 shows the order of accuracy for our test. Since this type of

discrete delta function does not have the moment condition, the convergent rate is nearly

of first order. In Table 4.14 we use δ_h√ to verify the conclusion we mention. In this case,

the discrete delta function holds the moment condition, thus we can obtain the same result as we proved in previous subsection. Figure 4.8 shows that the maximum error is improving as the mesh is refining.

mesh _{ku − U}hk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32 3.5924E-03 – 6.8746E-04 – 2.7290E-04 –

64 1.8555E-03 0.9531 2.5911E-04 1.4077 1.0378E-04 1.3948

128 9.0120E-04 1.0419 9.3268E-05 1.4740 4.3089E-05 1.2681

256 4.5930E-04 0.9724 3.6476E-05 1.3544 1.8986E-05 1.1823

512 2.2936E-04 1.0018 1.4989E-05 1.2830 9.4484E-06 1.0068

Table 4.13: Order of accuracy for one-dimensional test with δcos

0 0.2 0.4 0.6 0.8 1 −0.3 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 1D case comparison exact sol numer sol, N=32 0.625 0.63 0.635 0.64 0.645 0.65 −15 −10 −5 0 5x 10 −3 exact sol numer sol, N=512 numer sol, N=128

mesh _{ku − U}hk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32 3.7517E-03 – 6.9307E-04 – 1.5281E-04 –

64 1.8672E-03 1.0066 2.4420E-04 1.5049 3.8151E-05 2.0019

128 9.4225E-04 0.9867 8.6932E-05 1.4900 9.5640E-06 1.9960

256 4.6251E-04 1.0266 3.0318E-05 1.5197 2.3780E-06 2.0078

512 2.3990E-04 0.9470 1.1019E-05 1.4602 6.0116E-07 1.9838

Table 4.14: Order of accuracy for one-dimensional test with δ_h√ .

4.2.2

Two-dimensional problems 2

In this subsection, the problems we need to solve can be written generally in the following formulation :

∆u = Z

Γ

f (s) δ2_{(x − X(s)) ds + g} (58)

in the computational domain Ω = [a, b] × [c, d] with an interface Σ inside, while the boundary conditions are still Dirichlet type. We use the same mesh structure as it was stated in previous section. Here, we use regular centered difference scheme to discrete Eq. (58) as

Ui+1,j− 2Ui,j+ Ui−1,j

h2 +

Ui,j+1− 2Ui,j+ Ui,j−1

h2 = fi,j+ gi,j (59) where fi,j = X k f (sk) δh(xi− X(sk))δh(yj− Y (sk)) ∆s. (60)

Example 4.2.2.1: The first example is coming from LeVeque and Li’s work [22]. The

computational domain is Ω = [−1, 1] × [−1, 1], and the interface Σ inside Ω is a circle centered at (0, 0) with radius 0.5. The analytic solution u is written in polar coordinates as u(r) =      1 _{if r ≤ 0.5} 1 + log(Cr) if r > 0.5 , (61)

where the jump on the interface f = C. Here, C is set to be 2, and the external source term g is zero, which can be verified easily.

The convergence tests based on δcos

h are shown in table 4.15, while table 4.16 lists the

corresponding order of accuracy with using δ_h√ . Figure 4.9 shows the difference between

numerical and analytic solutions near the interface. The maximum error is decreasing after the mesh refined.

M × N ku − Uhk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32×32 5.6933E-02 – 2.7227E-02 – 3.8063E-02 –

64×64 2.8740E-02 0.9862 1.3437E-02 1.0187 1.8874E-02 1.0120

128×128 1.4463E-02 0.9907 6.7304E-03 0.9974 9.4448E-03 0.9987 256×256 7.2452E-03 0.9972 3.3629E-03 1.0009 4.7169E-03 1.0016 512×512 3.6247E-03 0.9991 1.6816E-03 0.9998 2.3582E-03 1.0001

Table 4.15: Convergent test using δcos

h for Example 4.2.2.1.

M × N ku − Uhk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32×32 3.2541E-02 – 1.0381E-02 – 5.6095E-03 –

64×64 1.6039E-02 1.0206 3.4241E-03 1.6000 1.3934E-03 2.0092

128×128 8.0821E-03 0.9887 1.2542E-03 1.4489 3.6329E-04 1.9393 256×256 4.0789E-03 0.9865 4.3528E-04 1.5267 9.1119E-05 1.9952 512×512 2.0854E-03 0.9678 1.5457E-04 1.4936 2.3011E-05 1.9854

0.4 0.45 0.5 0.55 0.6 0.9 0.95 1 1.05 1.1 1.15 2D case comparison 0.48 0.485 0.49 0.495 0.5 0.505 0.51 0.515 0.52 0.95 0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04 1.05 exact sol numer sol, N=32 exact sol numer sol, N=512 numer sol, N=128

Example 4.2.2.2: The next example is a modification from Zhou and his collabo-rators’ work [47]. The domain Ω of the problem is [−1, 1] × [−1, 1], which contains an circular interface Σ centered at (0, 0) with radius 0.5. For convenience, the problem is described in the polar coordinates. The exact solution u is written as

u(r) =      r2 _{if r ≤ 0.5} 7/64 + r4_{/4 + r}2_{/2 if r > 0.5} . (62)

The right hand side f is given by

g(r) =      4 _{if r ≤ 0.5} 4r2_{+ 2 if r > 0.5} , (63)

with the source term on the interface f = −0.375.

The table 4.17 provides the convergence results based on δcos

h , while the corresponding

order of accuracy with using δ_h√ is listed in table 4.18. Figure 4.10 shows the maximum

error at the interface improves as the mesh number N increases.

M × N ku − Uhk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32×32 1.0692E-02 – 6.8852E-03 – 9.8518E-03 –

64×64 5.3116E-03 1.0092 2.9823E-03 1.2070 4.3218E-03 1.1887

128×128 2.6178E-03 1.0207 1.3325E-03 1.1622 1.9038E-03 1.1827 256×256 1.3075E-03 1.0015 6.5268E-04 1.0297 9.2571E-04 1.0402 512×512 6.7053E-04 0.9634 3.2934E-04 0.9867 4.6689E-04 0.9874

Table 4.17: Convergent test using δcos

0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.1 0.15 0.2 0.25 0.3 0.35 0.4 2D case comparison 0.48 0.485 0.49 0.495 0.5 0.505 0.51 0.515 0.52 0.23 0.235 0.24 0.245 0.25 0.255 0.26 0.265 0.27 exact sol numer sol, N=32 exact sol numer sol, N=512 numer sol, N=128

M × N ku − Uhk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32×32 5.2739E-03 – 5.2757E-03 – 8.4898E-03 –

64×64 1.8829E-03 1.4858 1.7842E-03 1.5641 2.8787E-03 1.5603

128×128 1.1681E-03 0.6887 5.0356E-04 1.8250 7.9509E-04 1.8562 256×256 6.2013E-04 0.9135 1.9113E-04 1.3976 2.9956E-04 1.4082 512×512 3.1040E-04 0.9984 1.0044E-04 0.9282 1.5913E-04 0.9126

Table 4.18: Convergent test using δ_h√ for Example 4.2.2.2.

Example 4.2.2.3: _{The computational domain Ω of the final example is [−1, 1] ×}

[−1, 1]. There is an interface Σ inside the domain Ω, which is a circle centered at (0, 0) with radius 0.5. The analytic solution of u is given in the polar coordinates as

u(r) =      exp(−r2_{) if r ≤ 0.5}

1/(er2_{) − 4/e + exp(−1/4) if r > 0.5}

. (64)

The source term on the interface f = −16/e + exp(−1/4), while the right hand side g could be derived as g(r) =      4(r2_{− 1) exp(−r}2_{) if r ≤ 0.5} 4/(er4_{) if r > 0.5} . (65)

The convergence tests based on δcos

h are shown in table 4.19, while table 4.20 lists the

corresponding order of accuracy with using δ√_h . Figure 4.11 shows the numerical solution

−1 −0.5 0 0.5 1 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 2D case comparison 0.45 0.5 0.55 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 exact sol numer sol, N=32 exact sol numer sol, N=512 numer sol, N=128

M × N ku − Uhk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32×32 3.7152E-01 – 1.9625E-01 – 2.8326E-01 –

64×64 1.7040E-01 1.1245 8.1905E-02 1.2606 1.1891E-01 1.2522

128×128 7.7697E-02 1.1330 3.6186E-02 1.1785 5.1663E-02 1.2026 256×256 3.7131E-02 1.0652 1.7712E-02 1.0306 2.5111E-02 1.0408 512×512 1.8568E-02 0.9997 8.9366E-03 0.9869 1.2662E-02 0.9878

Table 4.19: Convergent test using δcos

h for Example 4.2.2.3.

M × N ku − Uhk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32×32 2.1936E-01 – 1.3733E-01 – 2.0976E-01 –

64×64 9.0241E-02 1.2815 4.5147E-02 1.6049 6.9359E-02 1.5965

128×128 3.6605E-02 1.3017 1.2348E-02 1.8703 1.8486E-02 1.9076 256×256 1.6330E-02 1.1645 4.6871E-03 1.3974 7.1454E-03 1.3713 512×512 8.1021E-03 1.0111 2.5386E-03 0.8846 3.9730E-03 0.8467

Table 4.20: Convergent test using δ_h√ for Example 4.2.2.3.

4.3

Applications to indicator functions and pressure

The indicator function is not only providing a method to determine whether the position is in the inside region or the outside part of the closed interface, but also a useful tool to solve other kinds of partial differential equations. For example, this function can be applied when solving the following second-order variable-coefficient elliptic equation:

∇ · (b∇Ψ) = f + Z

Σ

g δ2_{(x − X(s)) ds.} (66)

Here b is a piecewise constant coefficient matrix, i.e. b = b0 inside Σ (in Ω₀), and b = b1

The source term on the surface g can be computed by

g = (b1∇Ψ₁− b0∇Ψ₀) · n, (68)

which is the jump of the normal derivative of the solution across the interface.

By using the indicator function, one can solve this type of equations without deal-ing with complex domain or complicated conditions on the boundary (interface). For

convenience, we only use δ_h√ in the following computation.

4.3.1

One-dimensional problem 3

The equation (66) in one dimension could be rewrite as d dx  bdΨ dx  = f + gδ(x − Ip) , (69)

where Ip is the point that the interface lies in the domain.

To calculate the solution of equation (69), we adapt the simple discretization bj+1/2Ψj+1_∆x−Ψj − bj−1/2Ψj−Ψ∆xj−1

∆x = fi,j + gi,j (70)

Here bj+1/2 = (bj+1+bj)/2, bj−1/2 = (bj+bj−1)/2. Let b be a piecewise constant coefficient

in the compuational domain, where b = bL if the position x < Ip, and b = bR for the case

x > Ip. By using the indicator function I(x) under the definition

I(x) =      1 if x > Ip 0 if x < Ip , (71)

which can be calculated by solving the following equation

d2_I

dx2 =

d

dxδ(x − Ip). (72)

We can write the coefficient bj as bj = bL+ I(xj)(bR− bL).

Example 4.3.1: In the first example, we consider the equation (69) in the domain

[0, 4] with the exact solution

Ψ(x) =      (x + 1)2_/b L if x < Ip ((x + 1)2_{+ cx + d)/b} R if x > Ip . (73)

Here the interface point Ip = 2π/3, the right hand side function f = 2, and the jump

condition g = c = 1. The coefficients bL= 1 and bR = 10000, while the constant d in the

exact solution can be derived as

d = (Ip+ 1)2

 bR

bL − 1



− cIp (74)

Table 4.21 shows the accuracy analysis of our scheme. One can find about first-order convergence obviously. Figure 4.12 provides the comparison between analytic and numer-ical solutions. The maximum error of numernumer-ical solution comes from the neighborhood of the interface point, but it will be improved after we refine the mesh.

mesh _{kΨ − Ψ}hk∞ ratio kΨ − Ψhk2 ratio kΨ − Ψhk1 ratio

32 1.3088E-00 – 1.1052E-00 – 1.3808E-00 –

64 5.7502E-01 1.1865 4.8300E-01 1.1942 6.0531E-01 1.1897

128 3.8643E-01 0.5733 3.2351E-01 0.5781 4.0480E-01 0.5804

256 1.9618E-01 0.9780 1.6407E-01 0.9795 2.0549E-01 0.9781

512 1.0032E-01 0.9676 8.3856E-02 0.9683 1.0508E-01 0.9676

0 0.5 1 1.5 2 2.5 3 3.5 4 1 2 3 4 5 6 7 8 9 10 1D case comparison exact sol numer sol, N=32 numer sol, N=128 numer sol, N=512

Figure 4.10: Comparison between numerical and analytic solutions for 1D case

4.3.2

Two-dimensional problems 3

In order to solve equation (66) in two-dimensional domain, we extend previous discretiza-tion in one-dimensional case as

bi+1/2,jΨi+1,j_∆x−Ψi,j − bi−1/2,jΨi,j−Ψ∆xi−1,j ∆x +bi,j+1/2 Ψi,j+1−Ψi,j ∆y − bi,j−1/2 Ψi,j−Ψi,j−1 ∆y

∆y = fi,j+ gi,j, (75)

where

gi,j =

X

k

g(sk) δh(xi− X(sk))δh(yj− Y (sk)) ∆s. (76)

Here we approximate the coefficients without integer index by taking averages of nearby values on the grid

bi+1/2,j = (bi+1,j+ bi,j)/2 b_i−1/2,j = (bi,j+ bi−1,j)/2 bi,j+1/2 = (bi,j+1+ bi,j)/2

Example 4.3.2.1: The following example can be found in Beale and Layton’s work [4]. The computational domain is [−1.3, 1.3] × [−1.31.3]. There is an interface Γ, which

is an ellipse, inside the domain with major radius ra = 0.9 and minor radius rb = 0.7.

The whole problem is written in elliptic coordinates, which could defined by a conformal

mapping x + iy = a0cosh(ρ + iΘ), i.e. x = a0cosh ρ cos Θ and y = a0sinh ρ sin Θ. Here

a0 = pra2− r2b ≈ 0.565685 is the focus length of the ellipse. The interface Γ can be

represented by Γ = {ρ = ρ0 ≈ 1.039721, 0 ≤ Θ ≤ 2π} in elliptic coordinates, which

divides the domain into inner region Ω0 and outer part Ω1. The exact solution Ψ is

Ψ(ρ, Θ) =      a3

0(cosh2ρ sinh ρ cos2Θ sin Θ + sinh3ρ sin3Θ) if ρ < ρ0

c exp(−3ρ) sin(3Θ) + d exp(−ρ) sin Θ if ρ > ρ0

, (78)

where c ≈ 1.267135 and d ≈ 1.128542 are given to provide the continuity of solution

across interface. The constant coefficient inside the interface b0 is set as b0 = 0.2, while

the one of outer region b1 is equal to 100. Thus, the piecewise constant coefficient b can

be defined as b = b1+ I(x)(b0− b1). The right hand side function f to this problem is

f (ρ, Θ) =      8b0a0sinh ρ sin Θ if ρ < ρ0 0 if ρ > ρ0 , (79)

and the source term on the interface g can be derived by g = (b1∇Ψ1− b0∇Ψ0) · n, where

Ψ1 and Ψ0 are the limit of Ψ from outer region and that from inner area along normal

direction, respectively.

Figure 4.13 shows the comparison between analytic solution and numerical solution on different grids, where we use the snapshot of the values along x = 0. One can see the maximum error occurs inside the interface, and the scale of error is decreasing as the mesh size is getting smaller. Table 4.22 lists the convergence rates of numerical solutions in different norms, which is about first-order.

M × N ku − Uhk∞ ratio ku − Uhk2 ratio ku − Uhk1 ratio

32×32 2.1097E-01 – 1.6389E-01 – 2.0212E-01 –

64×64 1.2118E-01 0.7999 9.5354E-02 0.7813 1.1720E-01 0.7863

128×128 6.4839E-02 0.9021 5.1056E-02 0.9012 6.2913E-02 0.8974 256×256 3.3470E-02 0.9539 2.6293E-02 0.9574 3.2434E-02 0.9558 512×512 1.7063E-02 0.9719 1.3393E-02 0.9731 1.6537E-02 0.9718

Table 4.22: Convergent test for Example 4.3.2.1.

−1 −0.5 0 0.5 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 2D case comparison exact sol numer sol, N=32 numer sol, N=128 numer sol, N=512

Figure 4.11: Comparison between numerical and analytic solutions for 2D case 4.3.2.1.

Example 4.3.2.2: The final example of this section is modified from Stokes problem

in [14] and [23]. By using the pressure we have solved in example 4.1.2.2, we try to solve the velocity field by Eq. (1). The domain to the problem is Ω = [−2, 2] × [−2, 2] with an interface inside, which is a unit circle centered at (0, 0), i.e. X(θ) = (cos θ, sin θ). For