等位函數處理兩相流的數值研究

國

立

交

通

大

學

應用數學系

碩

士

論

文

等位函數處理兩相流的數值研究

Numerical Study of Level Set Method

for Two Phase Flows

研 究 生：郭乂維

指導教授：賴明治 教授

等位函數處理兩相流的數值研究

Numerical Study of Level Set Method

for Two Phase Flows

研 究 生：郭乂維 Student：Yi-Wei Guo

指導教授：賴明治 Advisor：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

in

Applied Mathematics June 2008

Hsinchu, Taiwan, Republic of China

等位函數處理兩相流的

數值研究

學生：郭乂維

指導教授

：

賴明治 教授

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

摘

要

本論文介紹使用等位函數結合 Navier-Stokes 方程處理兩相流問題的方法，

尤其是針對幾何結構有較大變動的流體，展現了良好的掌握能力。我們使用 WENO

五階精度來處理空間項，以及 TVD Runge-Kutta 三階精度來處理時間項。此外，

我們也嘗試使用等位函數來處理 Moving Contact Line 的問題。在最後，我們給

了許多流體的例子，來觀察使用等位函數的優點。

Numerical Study of Level Set Method

for Two Phase Flows

Student : Yi-Wei Guo

Advisor : Dr. Ming-Chih Lai

Department (Institute) of Applied Mathematics

National Chiao Tung University

Abstract

In this paper, we introduce that how to use level set function and Navier-Stokes

equations to solve two phase flows problems. Here, we will show the efficiency of

level set methods describe interface with large topological change. We use

Hamilton-Jacobi WENO, which allows us to discretize the spatial terms to fifth-order

accuracy. And we use TVD Runge-Kutta which is third-order accurate in time.

Additionally, we attempt to use level set function to solve moving contact line

problem. Finally, we show many numerical examples of two phase flow and verify the

excellence of level set function.

誌

謝

本篇論文的完成，首先要感謝我的指導老師－賴明治教授。在老師的指點之下，讓我逐漸 熟悉數值偏微分方程與科學計算等領域，並且引導我學習與研究的方向。除了課業上的知識以 外，老師也讓我學到許多求知與做研究的方法，讓我在數值方法與科學計算這片領域中，找到 我的興趣跟做出研究的成就感。更重要的，老師也教了我們做人處事的正確態度。除了指導老 師外，科學計算實驗室的夥伴們也給我不少幫助。感謝曾昱豪學長指點我 Navier-Stokes 方程 的理論和數值計算，讓我在研究的過程中能夠持續獲得進展，也在尋找研究資料時給我許多幫 助。感謝曾孝捷學長和陳冠羽學長指點我程式和學業，讓我在撰寫程式時更加得心應手，課業 上迷惑時能有適時的引導。感謝蔡修齊同學和我一起討論 Level Set Method，讓我研究起來 更加順遂。感謝李信儀學長以他雄厚的分析能力給我在理論方面的幫助，也感謝許多已經畢業 的學長姐留下來的資料。學生在此致上最誠懇的謝意。 在論文口試期間，承蒙王偉仲教授、洪子倫教授與吳金典助理教授費心審閱並提供許多寶 貴意見，使本論文得以更加完備，學生永銘在心。 除了研究生活外，我也要感謝蔡修齊同學、曾鈺傑同學、謝先皓同學、陳建明同學以及胡 偉帆同學。除了課業上的切磋討論外，閒暇之餘能跟他們閒話家常，交換各種意見與心得，讓 我的生活充滿樂趣，在互相鼓勵之下度過研究所的生涯。 最後，我要感謝我的家人，沒有他們的關心與祝福，我無法順利的完成學業。願與他們， 以及所有在我周圍關心我的人，一同分享此篇論文完成之喜悅與榮耀。

目

錄

中文提要 ……… i

英文提要 ……… ii

誌謝 ……… iii

目錄 ……… iv

一、 Introduction ……… 1

二、 Level Set Method ……… 2

2.1 Introduction ……… 2

2.2 Level Set Function ……… 3

2.3 Re-initialization of Level Set Functions …… 3

2.3.1 Introduction ……… 3

2.3.2 Hamilton-Jacobi ENO ……… 5

2.3.3 Hamilton-Jacobi WENO ……… 6

2.3.4 TVD Runge-Kutta ……… 7

2.4

Heaviside Function and Delta Function 8

2.5 Navier-Stoke Equations ……… 9

2.5.1 Boundary Conditions ……… 9

2.5.2 Marker and Cell Method ……… 9

2.5.3 Introduce Navier-Stokes Equations ……… 10

2.5.4 Governing equations ……… 11

2.5.6 The Stability Condition ……… 12

2.5.7 Second-Order Implicit Scheme ……… 12

2.5.8 Second-Order Projection Method ……… 13

2.5.9 Introduce Navier-Stokes Equations ……… 14

2.5.10 Numerical Scheme ……… 14

2.6 Convection Equation ……… 15

2.7 Summary ……… 16

三、 Moving Contact Line Problems ……… 17

3.1 Introduction ……… 17

3.2 Calculate Contact Angle ……… 17

3.3 Distribute the Force at the Contact Line 20

四、 Numerical Results ……… 23

4.1 A Rising Bubble ……… 23

4.2 Merging Two Bubbles with the Same Density … 24

4.3 Rayleigh-Taylor Instability ( I ) ……… 25

4.4 Rayleigh-Taylor Instability ( II ) ……… 26

4.5 Emit ( I ) ……… 27

4.6 Emit ( II ) ……… 28

4.7 Drip ( I ) ……… 29

4.8 Drip ( II ) ……… 30

4.10 Moving Contact Line (obtuse) ……… 32

五、 Concluding Remarks ……… 33

References ……… 34

1 Introduction

The fluid flows are in our daily life, such as the change of climate, and bubbles in the water. Taking the weather forecast as an example, accurate numerical scheme is very important. Therefore, hydrodynamics is popular and is important in computational science.

In this paper, we consider the incompressible, immiscible Navier-Stokes equations separated by a free surface. How to deal with the interface efficiently and accurately is the cynosure. There are many well know methods to describe interface, such as front-tracking method, volume of fluid method, level set method, etc. Here we focus on level set method, since it can describe the interface with large topological change conveniently. For example, Chang, Hou, Merriman and Osher [1] consider two phase flow with different viscosity and density, and the incompressible fluid flow is under the action of gravity with no slip boundary condition. We have considered the same problem as [1], but they use vorticity-stream function to solve Navier-Stokes equations, here we use MAC grid and the projection method to replace vorticity-stream function. Additionally, we will apply our numerical method to diversiform problems.

If the interface contacts with boundary, we must consider the moving contact line problems. Huang, Liang and Wetton [2] use Front-tracking method to dispose of the moving contact line problem. Additionally, M. Renardy, Y. Renardy, and Li [3] consider the moving contact line problem using a volume of fluid method. In this paper, we use level set method to solve the problems. There are two challenges needed to be overcome. First, how to calculate the dynamic equilibrium contact angle using level set method? Second, how to distribute the friction force at the contact line using level set method?

2 Level Set Method

2.1 Introduction

Refer to [4]. We define the level set function is a way to describe the interface, which is defined on all of computational domain Ω. If point A is farther from the interface than point B is, the value of point A must be larger than the value of point B after taking absolute value.

We define the distance function is a kind of level set function. The value of each point is to describe the proper distance from interface to each point.

We define the signed distance function is a distance function, but we allow that the value is negative. For example, to consider a distance function with the interface is a circle, and we let the value of the points which are contained in the circle multiply by -1. In general, in the case of two phase flows, positive sign and negative sign represent different parts of fluid as Figure 1.

Figure 1

2.2 Level Set Function

Consider the domain Ω = 0, 1 × [0, 2], we must define the level set function on Ω first.

Example:

a. Let ϕ x, y = (x − 0.5)2+ (y − 1)2− 0.2, where the interface is defined by

the ϕ x, y = 0, which is a circle with radius 0.2 and the center of a circle is at (0.5, 1).

b. In general, it is easy to define a level set function with ϕ x, y = 0 represents a

we want to set two curves that a circle with radius 0.2 and center of the circle is at (0.5, 0.5), and another circle with radius 0.1 and center of the circle is at (0.5, 1). What we must do is to combine, let

ϕ₁ = (x − 0.5)2_{+ (y − 0.5)}2_{− 0.2 and ϕ}

2 = (x − 0.5)2+ (y − 1)2− 0.1 ,

then compare the value of each point of ϕ₁with that of ϕ₂, and take the small one, we will obtain ϕ. That is exactly what we know that the value of each point is to describe the relative distance from interface to each point after taking modulus. In other words, ϕ x, y = min{ϕ₁ x, y , ϕ₂(x, y)}.

c. If the interface is a horizontal at an altitude of y = 1, and two kinds of fluid are

distributed over upside and downside of the interface respectively. If we want to set the level set function, it is easy to calculate the proper distance from interface to each point on Ω. And let the value of points which are contained in upside or downside multiply by -1, therefore we can compartmentalize the different parts of fluids.

d. If we want to set a level set function such that it describes the fluid distribution as

the following diagram, what must do is to combine the approach used in b and c.

Figure 2

2.3 Re-initialization of Level Set Functions

2.3.1 Introduction

distance function. In general, the level set function is defined at the beginning, we cannot promise which is a signed distance function. Even if the initial level set function is a signed distance function, it may not remain a signed distance function after the interface moves. So we must re-initialize the level set function at each time step.

Why we need a signed distance function? To calculate the curvature κ of interface, by definition,

κ = ∇ ∙ 𝑵 = ∇ ∙ (_|∇ϕ|∇ϕ) (1) where N is the interface normal and ϕ is our level set function. We should have signed distance function, otherwise the curvature would not be accurate, and certainly obtains imprecise result. Even if we don’t consider the curvature, in order to describe the topological change of interface, we use the simple level set equation

ϕt + u ∙∇ϕ = 0. (2) By the same reason that we need accurate ∇ϕ, it’s necessary to re-initialize the level set function at each time step.

In fact, on computational domain Ω, most points of level set function don’t have any effects except the points near the interface. Therefore, while doing re-initialization, we only consider the neighborhood of interface.

In order to obtain signed distance function, we solve the steady state solution of the following equation

ϕ_t+ S ϕ₀ (|∇ϕ −1 = 0 (3) where S ϕ₀ is a sign function taken as 1 in Ω+_{, -1 in Ω}−_{, and we want ϕ to stay}

identically equal to zero.

An important characteristic of signed distance function ϕ is |∇ϕ| = 1 , so it’s standard while doing re-initialization. We calculate |∇ϕ| of the points whose distance from the interface is less than 2.5∆x, where ∆x is the numerical grid size. In the case that ∆x = 1/128, if total error is less than 3 × 10−6, we accept the signed distance function is accurate enough.

An example is as follow, the interface is an ellipse with major axis 0.35 and minor axis 0.2 initially. Here, ϕ x, y = ( x − 0.5 2/0.352) + ( y − 0.5 2/0.22) − 1 with computational domain Ω = 0,1 × [0,1].

Figure 3 : Grid 128 × 128, ∆t = 0.25∆x. Contour taken over [−10∆x: 2∆x: 10∆x].

2.3.2 Hamilton-Jacobi ENO (Essentially Nonoscillatory)

In order to understand HJ WENO, we must introduce HJ ENO first. If readers interest in this, you can refer to [4], [5], [6], and [7]. Giving some definitions as follows, 𝐷+_{ϕ = ϕ} x + ₌ ϕi+1−ϕi ∆x 𝐷−_{ϕ = ϕ} x − ₌ ϕi−ϕi−1 ∆x

ϕ and then differentiate to get ϕ_x. Define 𝐷_i0 _{ϕ = ϕ} i 𝐷_i+1/21 _{ϕ =}𝐷i+10 ϕ−𝐷i0 ϕ ∆x 𝐷_i2 _{ϕ =}𝐷i+1/21 ϕ−𝐷i−1/21 ϕ 2∆x 𝐷_i+1/23 _{ϕ =}𝐷i+12 ϕ−𝐷i2 ϕ 3∆x ϕ x = 𝑄₀ x + 𝑄₁ x + 𝑄₂ x + 𝑄₃ x And we have ϕ_x x_i = 𝑄₁′ _x i + 𝑄2′ xi + 𝑄3′ xi First, let 𝑄₁ x = (𝐷_k+1/21 _{ϕ)(x − x} i) 𝑄₁′ _x i = 𝐷k+1/21 ϕ

Here, we start with k = i − 1 to find ϕ_x−, and start with k = i to find ϕ_x+. And ϕ_x x_i = 𝑄₁′ x_i is the first-order accuracy.

Second, let

𝑄₂ x = c x − x_k (x − x_k+1) 𝑄₂′ _x

i = c(2 i − k − 1)∆x

Here, k is the same as above. If 𝐷_k2 _{ϕ ≤ 𝐷}

k+12 ϕ , we set c = 𝐷k2 ϕ and

k∗ _{= k − 1, otherwise, we set c = 𝐷}

k+12 ϕ and k∗ = k.

ϕ_x x_i = 𝑄₁′ _x

i + 𝑄2′ xi is the second-order accuracy.

Third, let

𝑄₃ x = c∗ _{x − x}

k∗ (x − xk∗+1)(x − xk∗+2)

𝑄₃′ _x

i = c∗(3 i − k∗ 2− 6 i − k∗ + 2)(∆x)2

If 𝐷_k3∗_+1/2 ϕ ≤ 𝐷_k3∗_+3/2 ϕ , we set c∗ = 𝐷_k3∗_+1/2 ϕ . Otherwise we set

c∗ _{= 𝐷}

k3∗+3/2 ϕ .

Now, ϕ_x x_i = 𝑄₁′ x_i + 𝑄₂′ x_i + 𝑄₃′ x_i is the third-order accuracy comes from HJ ENO.

2.3.3 Hamilton-Jacobi Weighted ENO (WENO)

Here, we use Hamilton-Jacobi WENO, which allows us to discretize the spatial terms in equation (3) to fifth-order accuracy. We take (ϕ_x−₎

i as an example to introduce HJ

WENO, and consider the case that first-order approximation of (ϕ_x−)_i is ϕi−ϕi−1

∆x . To

calculate (ϕ_x−₎

definitions as follows, 𝑣₁ = ϕi−2−ϕi−3 ∆x , 𝑣2 = ϕi−1−ϕi−2 ∆x , 𝑣3 = ϕi−ϕi−2 ∆x , 𝑣4 = ϕi+1−ϕi ∆x , 𝑣1 = ϕi+2−ϕi+1 ∆x .

According to these definitions, allows us to write

ϕ1_{x =} 𝑣1 3 − 7𝑣2 6 + 11𝑣3 6 , ϕ_x2 _{= −}𝑣2 6 + 5𝑣3 6 + 𝑣4 3, ϕ_x3 ₌𝑣3 3 + 5𝑣4 6 − 𝑣5 6,

as the three potential HJ ENO approximations to (ϕ_x−₎

i with third-order accuracy. To

choose the single approximation with the smallest error by choosing the smoothest possible polynomial interpolation of ϕ is the goal of HJ ENO. But in the case that interface is smooth, it is overkill to pick exactly one of three candidate stencils. So, the weighted ENO (WENO) method is to take a convex combination of the three ENO approximations. Giving some definitions as follows again,

𝑆₁ =13₁₂(𝑣₁− 2𝑣₂+ 𝑣₃)2 ₊1 4(𝑣1− 4𝑣2+ 3𝑣3)2, 𝑆₂ =13₁₂(𝑣₂− 2𝑣₃+ 𝑣₄)2₊1 4(𝑣2− 𝑣4)2, 𝑆₃ =13₁₂(𝑣₃− 2𝑣₄+ 𝑣₅)2₊1 4(3𝑣3− 4𝑣4+ 𝑣5)2, α1 =_(𝑆0.1 1+ε)2 α2 =_(𝑆0.6 2+ε)2 α3 =_(𝑆0.3 3+ε)2 ε = 10−6_{max 𝑣}

12, 𝑣22, 𝑣32, 𝑣42, 𝑣52 + 10−99, where the 10−99 term is set to avoid

division by zero in the definition of the α_𝜅. Let ω_n = αn

α₁+α₂+α₃, n = 1, 2, 3.

Then ϕ_x = ω₁ϕ_x1+ ω₂ϕ_x2 + ω₃ϕ_x3 is fifth-order accurate by WENO.

2.3.4 TVD Runge-Kutta

In order to introduce TVD Runge-Kutta, we take the equation ϕ_t+ u ∙▽ϕ = 0

as an example. The forward Euler time discretization in equation (4) is only first-order accurate in time. Therefore, we use TVD Runge-Kutta which is third-order accurate in time. The third-order accurate TVD RK scheme proposed in [6] is as follows.

First, we take an Euler step to advance the solution to time tn _{+ ∆t,}

ϕn +1_−ϕn

∆t + u

n_{∙ ∇ϕ}n _{= 0.}

Taking an Euler step again to advance the solution to time tn_{+ 2∆t,}

ϕn +2_−ϕn +1

∆t + u

n+1_{∙ ∇ϕ}n+1 _{= 0.}

Using ϕn and ϕn+2, we can take an averaging step,

ϕn+12 =3

4ϕ n₊1

4ϕ n+2_.

That produces an approximation to ϕ at time tn₊1 2∆t.

Now, another Euler step is taken to advance the solution to time tn₊3 2∆t, ϕn +32−ϕn +12 ∆t + u n+1_{2∙ ∇ϕ}n+1_{2 = 0.} Using ϕn and ϕn+ 3

2_{, we can take an averaging step again}

ϕn+1 ₌1 3ϕ

n ₊2 3ϕ

n+3_2.

That produces a third-order accurate approximation to ϕ at time tn_{+ ∆t.}

2.4 Heaviside Function and Delta Function

In fact, the densities of two kinds of fluids change immediately as passing the interface. But, in the cause of numerical calculation, we need to maintain the continuity. Here, we refer to [1].

We define the regularized delta function δ_ε as

δε ϕ = 1 2ε 1 + cos πϕ ε if |ϕ| < 𝜀 0 otherwise (5)

𝐻ε ϕ =

0 if ϕ < −ε ϕ + ε 2ε+ sin πϕ_ε 2π if ϕ ≤ ε 1 if ϕ > ε

(6)

Where ϕ is a signed distance function comes from the re-initialization.

It’s easy to check the Heaviside function satisfies the relation 𝑑𝐻_ε(ϕ) 𝑑 ϕ = δ_ε(ϕ). Using the regularized Heaviside function 𝐻_ε, we can define the corresponding continual density function ρ and the continual viscosity function μ as

ρ_ε 𝑥 = ρ₁+ ρ₂− ρ₁ 𝐻_ε(ϕ(𝑥)) μ_ε 𝑥 = μ₁+ μ₂− μ₁ 𝐻_ε(ϕ(𝑥)).

We can consider ρ_ε and μ_ε as smooth variable density and variable viscosity.

2.5 Navier-Stokes Equations

2.5.1 Boundary Conditions

Figure 4

1. No slip condition : The continuous velocities should vanish at the boundary. 2. Free slip condition : The velocity component normal to the boundary should

vanish along with the normal derivative of the velocity component tangent to the boundary, as Figure 4, i.e. 𝜕𝑣(A or B,y)

∂𝑛 =

𝜕𝑢 (C or D,x) ∂𝑛 = 0

In this paper, we use the no slip boundary condition.

2.5.2 Marker and Cell Method

Refer to [4], [9] and [12]. We use marker and cell (MAC) method and let the computation domain Ω be squared with staggered grid in Cartesian coordinates. This specially defined grid decomposes the computational domain into cells with velocities and forces defined on the cell faces and scalars defined at cell centers such as pressure, density and viscosity. We take the mesh size h = ∆x = ∆y and let 𝑁_𝑥, 𝑁_𝑦 be the

number of mesh in the x and y directions respectively. In general, we take the same number of grid points in each direction, and denote that 𝑁 = 𝑁_𝑥 = 𝑁_𝑦. MAC grid is as the following diagram comes from [9],

Figure 5 : The computational domain Ω using staggered grid with mesh h.

First Part : Density and Viscosity Are Constant

2.5.3 Introduce Navier-Stokes Equations

The Navier-Stokes equations are

∂w

∂t

+ w

∙ ∇ w

+ ∇𝑝 =

1

𝑅𝑒

∆w

+ 𝑓

∇ ∙ w

= 0

(

7

)

where w = _{𝑣(𝑡,𝑥,𝑦)}u(t,x, y) = the velocity in the 𝑥 − direction_{the velocity in the 𝑦 − direction}

and p : p(t,x,y) denotes the pressure. f(t,x,y) is the external force defined on the computational domain, which contains the forces resulted from the following, a. Surface tension. The form is τκ(ϕ)∇ϕδ(ϕ), where τ is the surface

and κ ϕ = ϕy

2_{ϕxx−2ϕxϕyϕxy+ϕx}2_ϕyy

(ϕ_x2_+ϕ y

2₎3/2 ,

δ is a one-dimension Dirac delta function and ϕ is our signed distance function.

b. Gravity.

c. Frictional force of moving contact line.

2.5.4 Governing equations

For simplicity, consider the case ρ ≡ 1, low Re and zero external force : 𝑢t+ 𝑢𝑢x + 𝑣𝑢y + 𝑝x =_𝑅𝑒1 (𝑢xx + 𝑢yy)

𝑣_t+ 𝑢𝑣_x + 𝑣𝑣_y + 𝑝_x = _𝑅𝑒1 (𝑣_xx + 𝑣_yy) 𝑢_x + 𝑣_y = 0

(8)

𝑢 0, 𝑥, 𝑦 , 𝑣 0, 𝑥, 𝑦 , 𝑢|_𝜕Ω, and 𝑣|_𝜕Ω are know.

2.5.5 Pressure Poisson’s Equation

Define the time derivative : [𝑢_t]n+1 = 𝑢n +1−𝑢n

∆t + 𝑂(∆t)

Therefore we need to solve the coupling system : 𝑢n+1 _{− 𝑢}n ∆t + (𝑢𝑢x+ 𝑣𝑢y)n + 𝑝xn+1 = 1 𝑅𝑒∆𝑢n 𝑣n+1_{− 𝑣}n ∆t + (𝑢𝑣x + 𝑣𝑣y)n+ 𝑝yn+1 = 1 𝑅𝑒∆𝑣n 𝑢_xn+1_{+ 𝑣} yn+1 = 0 [𝑢_x]_i+1 2,j = 𝑢_i+3 2,j − 𝑢i−12,j 2ℎ + 𝑂(ℎ2) [𝑢y]_i+1 2,j = 𝑢_i+1 2,j+1 − 𝑢i+12,j−1 2ℎ + 𝑂(ℎ2) [𝑢_xx]_i+1 2,j = 𝑢_i+3 2,j − 2𝑢i+12,j + 𝑢i−12,j 2ℎ + 𝑂(ℎ2) [𝑢_yy]_i+1 2,j = 𝑢_i+1 2,j+1− 2𝑢i+12,j + 𝑢i+12,j−1 2ℎ + 𝑂(ℎ2) [𝑣]_i+1 2,j = 𝑣_i,j−1

2+ 𝑣i+1,j−12 + 𝑣i,j+12+ 𝑣i+1,j+12

4 + 𝑂(ℎ2)

Similar arguments for v can be done. Discretizing (8), we obtain

𝑢 i+1_2,j n+1 _{= 𝐹} i+1_2,j n ₋∆t ℎ (𝑝i+1,jn+1 − 𝑝i,jn+1) Where 𝐹 i+1_2,j n _{= 𝑢} i+1_2,j n _{+ ∆t(} 1 𝑅𝑒∆𝑢 − 𝑢𝑢x − 𝑣𝑢y)i+1_2,j n And 𝑣 i,j+1₂ n +1 _{= 𝐺} i,j+1₂ n ₋∆t ℎ (𝑝i,j+1n+1 − 𝑝i,jn+1) Where 𝐺 i,j+1₂ n _{= 𝑣} i,j+1₂ n _{+ ∆t(} 1 𝑅𝑒∆𝑣 − 𝑢𝑣x − 𝑣𝑣y)i,j+1₂ n

From the continuity equation, we have (∇ ∙ w )_i,jn+1 ₌ 𝑢 i+1_2,j n+1 _{− 𝑢} i−1_2,j n+1 ℎ + 𝑣 i,j+1₂ n+1 _{− 𝑣} i,j−1₂ n+1 ℎ = 0

And then the resulting equation is (𝑝i+1,j− 2𝑝i,j+ 𝑝i−1,j

ℎ2 +

𝑝_i,j+1− 2𝑝_i,j + 𝑝_i,j−1

ℎ2 )n+1 = 1 ∆t( 𝐹_i+1 2,j − 𝐹i−12,j ℎ + 𝐺_i,j+1 2− 𝐺i,j−12 ℎ )n

This implies that (∆𝑝)_i,jn+1 = _∆t1 ∇ ∙ (𝐹_G)_i,jn

2.5.6 The Stability Condition

Courant-Friedrich-Lewy (CFL) conditions :

|𝑢|_max∆t < ∆𝑥, 𝑎𝑛𝑑 |𝑣|_max∆t < ∆𝑦

which states that no fluid particle can travel a distance greater then mesh spacing ∆x or ∆y, in the time ∆t. The other constraint is _𝑅𝑒∆t(_∆x1₂+_∆y1₂) <1₂, coming from the diffusion term. Therefore, we choose

∆t = τ ∙ min {𝑅𝑒 2 1 ∆x2+ 1 ∆y2 −1 , ∆x 𝑢 _max , ∆y 𝑣 _max} for a stable computation and τ is a safety factor.

2.5.7 Second-Order Implicit Scheme

[𝑢_t]n+1 ₌3𝑢n +1−4𝑢n+𝑢n −1 2∆t + 𝑂(∆t 2_{), and} 𝑢𝑢_x+ 𝑣𝑢_y n+1 = 𝑢𝑢_x+ 𝑣𝑢_y |_n−1n = 2[𝑢𝑢_x + 𝑣𝑢_y]n _{− [𝑢𝑢} x+ 𝑣𝑢y]n−1 + 𝑂(∆t2).

Therefore we need to solve the coupling system :

3𝑢n +1−4𝑢n+𝑢n −1 2∆t + [𝑢𝑢x + 𝑣𝑢y]|n−1 n _{+ 𝑝} xn+1 = _𝑅𝑒1 ∆𝑢n+1 3𝑣n +1−4𝑣n+𝑣n −1 2∆t + [𝑢𝑣x + 𝑣𝑣y]|n−1n + 𝑝yn+1 = 1 𝑅𝑒∆𝑣n+1 𝑢xn+1+ 𝑣yn+1 = 0

2.5.8 Second-Order Projection Method

We use second-order projection method to solve Navier-Stokes equations.

Step1. Using Helmholtz-type solver to solve

3w ∗_−4w n_+w n −1

2∆t + w ∙ ∇ w |n−1

n _{+ ∇𝑝}n ₌ 1 Re∆w

∗ ₍₉₎

Where w ∗ is the intermediate velocity field between w n and w n+1. Here w ∗ is solved implicitly, the constraint of ∆t will be diminished.

Step2. From Hodge decomposition, since w ∗ is a vector field, there are a potential function φn+1 and a divergence-free vector field w n+1 such that

w n +1−w ∗

2∆t/3 = −∇φn+1 (10)

And if we take the divergence operator into (10), we obtain ∆φn+1 ₌ 3

2∆t∇ ∙

w

∗

(11) A Poisson solver is used to obtain φn+1.

Step3. Projection from w ∗

onto

w n+1. w

n+1 _{= w} ∗₋2∆t

3 ∇φn+1 (12)

Step4. It’s clearly that we need to get ∇𝑝n+1. Write w ∗ as the form of the combination of w n+1 and ∇φn+1, and substituting it into (9). Finally, compare the new schemes and the coupling system term by term, we obtain

∇𝑝n+1 _{= ∇𝑝}n _{+ ∇φ}n+1₋ 1

Second Part : Density and Viscosity Are Not Constant

2.5.9 Introduce Navier-Stokes Equations

ρ 𝐮t+ ∇ ∙ 𝐮𝐮 = −∇𝑝 + ∇ ∙ 2μ𝐃 + 𝑓

∇ ∙ 𝐮 = 𝟎

(14)

where u is velocity, ρ is density, p is pressure, μ is viscosity, f is the external force

and D is the rate of deformation tensor, whose components are 𝐷_𝑖,𝑗 = 0.5(𝑢_𝑖,𝑗 + 𝑢_{𝑗 ,𝑖}).

2.5.10 Numerical Scheme

In order to achieve our scheme, we refer to [12] and deeply appreciate Y. H. Tseng provides code of solving Navier-Stokes equations.

We solve Navier-Stokes equations by two steps as follows.

Step1. Prediction:

First, we solve it by Crank-Nicholson method in time ∇ ∙ 𝜇n+1_∇𝐮∗ _{− 2ρ}n+1₂𝐮∗ ∆t = 2∇𝑝 n−1_{2− 2ρ}n+1₂𝐮n ∆t − ∇ ∙ (𝜇∇𝐮)n + (3𝐅n − 𝐅n−1) (15) where 𝐅 = ρ 𝑢𝑢x+ 𝑣𝑢y − (μ𝑣x)y − (μ𝑢x)x − 𝑓1 ρ 𝑢𝑣_x + 𝑣𝑣_y − (μ𝑣_y)_x − (μ𝑢_y)_y − 𝑓₂ and 𝑓 = (𝑓₁, 𝑓₂) means the external force term.

Note that 𝜇n+1 and ρn+

1

2 _{are computed by extrapolation as follows,}

𝜇n+1 _{= 2𝜇}n _{− 𝜇}n−1

ρn+12 =3

2ρn − 1 2ρn−1

and (15) becomes Poisson type equation and it can be solved by linear system for 𝐮∗.

Step2. Projection:

Since the new 𝐮∗ is not diverge free, we apply this step to correct 𝐮∗ to be diverge free which means ∇ ∙ 𝐮 = 0 and then compute new pressure term ∇𝑝n+

1 2_{. By}

Helmholtz-Hodge decomposition, we have

𝐮∗ _{= 𝐮}n+1 _{+ ∆t∇𝜓}n+1

∇ ∙ 𝐮n+1 _{= 0}

ρn+12𝐮

n+1_{− 𝐮}∗

∆t = −∇𝜓n+1 and then take divergence we have

∇ ∙ ∇𝜓n+1 ρn+12

= ∇ ∙ 𝐮∗ ∆t After obtaining ∇𝜓n+1 we have

𝐮n+1 _{= 𝐮}∗₋ ∆t

ρn+12

∇𝜓n+1

∇𝑝n+12 = ∇𝑝n−12+ ∇𝜓 + ∇ ∙ (𝜇n+1(𝐮n+1− 𝐮∗))

2.6 Level set Equation

In order to define the evolution of our implicit function ϕ, we use the simple level set equation ϕ_t + u ∙▽ϕ = 0, where ϕ is a signed distance function, and u is the velocity. By solving the Navier-Stokes equations, we can obtain the velocity u of each point on computational domain at this time step. Therefore, according to the equation ϕ_t+ u ∙▽ϕ = 0, we can get ϕ at next time step. Here, we use HJ WENO to discretize the spatial terms to fifth-order accuracy and TVD Runge-Kutta to get third-order accuracy in time.

2.7 Summary

Refer to figure 6, our procedure is as follows.

1. Construct the level set function on computational domain where ϕ 𝑥 = 0

stands for the interface.

2. In order to adjust the level set function to a signed distance function, we must

re-initialize the level set function. It’s important that re-initialization doesn’t change the shape of interface.

3. Define Heaviside function and use it to define the corresponding regularized

density function and the regularized viscosity function.

4. Solve the Navier-Stokes equations to obtain the velocity of each point on

computational domain at this time step.

5. Using the velocity obtained from Navier-Stokes equations and combining it with

the level set equation, we can get the level set function at next time step.

6. Remark that the level set function obtained from process 5 is not a signed distance

function, therefore we should return to process 2.

3 Moving Contact Line Problems

3.1 Introduction

The moving contact line problem is very practical and interesting, because it is in our daily life. For example, the capillarity, lubricated pipeline transport, fibers, and the bead on a leaf. We know that a bead on the floor and a drop of mercury on the floor whose shapes are different clearly, this phenomenon is mainly caused by surface tension.

Figure 7

Refer to [2]. The surface tension occurs on interface, therefore as the above diagram, let ς_sg, ς_sl and ς be the surface tension coefficients between solid-gas, solid-liquid, and gas-liquid phases, respectively. Because the different surface tensions, there are three forces at points A as Figure 6. These forces would attain to balance, at this time, ς_sg, ς_sl, ς and the contact angle satisfy the well-known Young’s equation

ς_sg = ς_sl + ς ∙ cosθ₀ (16) here θ₀ is the static equilibrium contact angle.

When the shape of interface changes, the contact angle θ will differ from θ₀ until an equilibrium is established. Let θ_𝑑 represent the dynamic equilibrium contact angle at the contact line, the friction force F_s at the contact line satisfies the following equation

F_s = ς_sg − ς_sl− ς ∙ cosθ_𝑑 (17) By (16) and (17), we can obtain that F_s = ς ∙ (cosθ₀− cosθ_𝑑).

It’s clearly there are two challenges needed to be overcome as following if we use level set methods to solve moving contact line problems.

1. How to calculate the dynamic equilibrium contact angle using the level set method?

3.2 Calculate Contact Angle

In order to obtain the friction force at the contact line, from equation (17), we need to calculate the dynamic equilibrium contact angle at the contact line using level set method.

First, we know that interface which is near the boundary approximates tangent line of the interface. It’s reasonable if the mesh size is small enough.

Here, we assume that the level set function is a signed distance function. As the following diagram,

Figure 8

where θ is the contact angle, and the signed distance function is defined at cell center. What we know are y, z and mesh size ∆x and unknown is h, by the property of similar triangles, we have

𝐲 𝐡 = 𝐳 ∆x − 𝐡 𝐲 ∙ ∆x − 𝐲𝐡 = 𝐳𝐡 𝐲 ∙ ∆x = (𝐳 + 𝐲)𝐡 𝐡 = 𝐲 ∙ ∆x 𝐳 + 𝐲 sinθ = 𝐲 𝐡= 𝐲 + 𝐳 ∆x

Now, we obtain the sine of contact angle, but since F_s = ς ∙ (cosθ₀− cosθ_𝑑), what we need is the cosine of contact angle. Although we have sin2_{θ + cos}2_{θ = 1, but we}

can’t ascertain that cosθ is positive or negative. In other words, we should determine that contact angle is acute or obtuse.

Refer to figure 7, we define the part of signed distance function which represents liquid is negative, the other part which represents gas is positive. Remark that signed distance function is defined at cell centers. Let ϕ 𝑖, 𝑗 be our signed distance function, if ϕ 𝑖, 1 > 0 and ϕ 𝑖 + 1,1 < 0 for some 𝑖, it shows that interface is between the two points. For this 𝑖, if ϕ 𝑖, 2 > ϕ 𝑖, 1 , we can say that contact angle is acute, otherwise, the contact angle is obtuse. As the following diagrams,

Figure 9 : ϕ 𝑖, 2 > 0, ϕ 𝑖, 1 > 0, ϕ 𝑖 + 1,1 < 0 and ϕ 𝑖, 2 > ϕ 𝑖, 1 , θ is acute.

Figure 11 : ϕ 𝑖, 2 > 0, ϕ 𝑖, 1 > 0, ϕ 𝑖 + 1,1 < 0 and ϕ 𝑖, 2 < ϕ 𝑖, 1 , θ is obtuse.

By the same way, we can calculate the dynamic equilibrium contact angle using level set method in the other case. Therefore, we have F_s .

3.3 Distribute the Friction Force at the Contact Line

We have calculated the friction force at the contact line, but the force terms are defined at cell faces. Therefore we need to distribute the friction force at the contact line using discrete delta function as (5)

δ_ε ϕ = 2ε1 1 + cos πϕ

ε if |ϕ| < 𝜀

0 otherwise

Here, taking ε = 2.5∆x, ∆x is the mesh size.

Now, we should calculate distance between the source of force and the point of MAC grid where the force terms are defined at. As the following diagram

Figure 12

According to our MAC grid, force is defined at cell faces like point A and B. Signed distance function is defined at cell center like point a and b, therefore y and z we have known. Now, we need to calculate Ao and Bo .

cb = 𝐳 𝐲 + 𝐳∙ ∆x dc = do ∙ cotθ = 1 2∆x ∙ cotθ dB = dc + cb − Bb = 1 2∆x ∙ cotθ + 𝐳 𝐲 + 𝐳∙ ∆x − 1 2∆x dA = ∆x − dB = 3 2∆x − 1 2∆x ∙ cotθ − 𝐳 𝐲 + 𝐳∙ ∆x Ao = dA 2_{+ do} 2 Bo = dB 2_{+ do} 2

By the same way, we can calculate distance between the source of force and the point of MAC grid where the force terms are defined at. To combine discrete delta function (5), it’s easy to distribute the friction force at the contact line.

Here, we let the incidence of the friction force at the contact line be a circle with radius 2.5∆x . Refer to the following diagram,

Figure 13

By section 3.2, we can obtain the friction force at the contact line F_s . By section 3.3, we can obtain the lengths of F , … , F_sF1 respectively. _sF9

We know that the incidence of discrete delta function (5) is a circle, but we let the force outside of the boundary be zero. In order to maintain the integral of delta function is 1, we let the discrete delta function inside of the boundary multiply by 2. Taking the force F1 which is resulted from the friction force at the contact line as an example,

F1 = Fs × 2 × δε F sF1

4 Numerical Results

4.1 A Rising Bubble

In this easy case, we try to simulate a rising bubble in the water under the influence of gravity. The fluid is at a standstill initially. Our computational domain Ω = 0,1 × [0,2], using 128× 256 grid points with ∆t = 0.25∆x. Viscosity for the fluid inside the bubble is equal to μ = 0.001873. Viscosity for the fluid outside the bubble is equal to μ = 0.00089. We take the density inside the bubble to be 0.8 and the density outside the bubble to be 1. The surface tension is 0.05. The initial position of the bubble is a circle. The circle is centered at (0.5, 0.65) with radius 0.15. In Figure 4.1, we plot the evolution of the bubble at time T = 0,100₁₂₈,200₁₂₈,300₁₂₈,₁₂₈400,500₁₂₈,600₁₂₈,700₁₂₈ .

4.2 Merging Two Bubbles with the Same Density

Here, we compute the interaction of two bubbles of the same density under the influence of gravity. Our computational domain Ω = 0,1 × [0,2], using 128× 256 grid points with ∆t = 0.25∆x. The fluid is at a standstill initially. Viscosity for the fluid inside the two bubbles is equal to μ = 0.001873. Viscosity for the fluid outside the two bubbles is equal to μ = 0.00089. We take the density inside the two bubbles to be 0.8 and the density outside the bubbles to be 1. The surface tension is 0.05. The initial position of the two bubbles is two circles. The lower one is centered at (0.5, 0.35) with radius 0.1. The upper one is centered at (0.5, 0.65) with radius 0.15. In Figure 4.2, we plot the evolution of the two bubbles at time T = 0,250₅₁₂,500₅₁₂,750₅₁₂,1000₅₁₂,

1250 512 , 1500 512 , 1750 512 .

4.3 Rayleigh-Taylor Instability ( I )

Here, we compute the Rayleigh-Taylor instability under the influence of gravity. Our computational domain Ω = 0,1 × [0,2], using 128 × 256 grid points with ∆t = 0.025∆x. The fluid is at a standstill initially. Viscosity for the upper fluid is equal to μ = 0.05. Viscosity for the lower fluid is equal to μ = 0.01. We take the density of the upper fluid to be 10 and the density of the lower fluid to be 1. The surface tension is 0.05. Initially, the two kinds of fluids are separated by the interface. The interface is a part of circle. The circle is centered at (0.5, 2.3) with radius 1.2. In Figure 4.3, we plot the evolution at time T = 0,₂₅₆9 ,₂₅₆18 ,₂₅₆27 ,₂₅₆36 ,₂₅₆45 ,₂₅₆54 ,₂₅₆63 .

4.4 Rayleigh-Taylor Instability ( II )

In this case, the two kinds of fluids are separated by the interface initially. The interface is a part of circle. The circle is centered at (0.5, -0.3) with radius 1.5. All other parameters are as case 4.3. In Figure 4.4, we plot the evolution at time T = 0,₂₅₆9 ,₂₅₆18 ,₂₅₆27 ,₂₅₆36 ,₂₅₆45 ,₂₅₆54 ,₂₅₆63 .

4.5 Emit ( I )

Here, we try to simulate the air bubble emits from water surface under the influence of gravity. Our computational domain Ω = 0,1 × [0,2], using 128× 256 grid points with ∆t = 0.025∆x. The fluid is at a standstill initially. Viscosity for the fluid A is equal to μ = 0.001. Viscosity for the fluid B is equal to μ = 0.5. We take the density of the fluid A to be 1 and the density of the fluid B to be 10. The surface tension is 0.05. The bubble is a circle which is centered at (0.5, 0.6) with radius 0.1. In Figure 4.5, we plot the evolution of the bubble at time T = 0,120₅₁₂,276₅₁₂,288₅₁₂,300₅₁₂,312₅₁₂,

324 512,

360 512 .

4.6 Emit ( II )

Most parameters are as 4.5. But we take the density of the fluid A to be 1 and the density of the fluid B to be 100 here. In Figure 4.6, we plot the evolution of the bubble at time T = 0,120₅₁₂,240₅₁₂,246₅₁₂,252₅₁₂,258₅₁₂,264₅₁₂,360₅₁₂ .

4.7 Drip ( I )

Here, we try to simulate the water drop drips into the water under the influence of gravity. Our computational domain Ω = 0,1 × [0,2], using 128× 256 grid points with ∆t = 0.025∆x. The fluid is at a standstill initially. Viscosity for the fluid A is equal to μ = 0.089. Viscosity for the fluid B is equal to μ = 0.1873. We take the density of the fluid A to be 10 and the density of the fluid B to be 8. The surface tension is 0.1. The bubble is a circle which is centered at (0.5, 1.2) with radius 0.1. In Figure 4.7, we plot the evolution of the bubble at time T = 0,330₅₁₂,340₅₁₂,350₅₁₂,360₅₁₂,

370 512, 380 512, 500 512 .

4.8 Drip ( II )

Most parameters are as 4.7. But we take the density of the fluid A to be 10 and the density of the fluid B to be 5 here. Viscosity for the fluid A is equal to μ = 0.2. Viscosity for the fluid B is equal to μ = 0.09. The surface tension is 0.05. In Figure 4.8, we plot the evolution of the bubble at time T = 0,100₅₁₂,170₅₁₂,200₅₁₂,240₅₁₂,270₅₁₂,

280 512,

320 512 .

4.9 Moving Contact Line (acute)

In this case, we compute the behavior of a bubble on the floor under the influence of the friction force at the contact line without gravity. Our computational domain Ω = 0,1 × [0,1], using 128× 128 grid points with ∆t = 0.025∆x. We set the static equilibrium contact angle θ₀ = 90. The fluid is at a standstill initially. The interface is a part of circle. The circle is centered at (0.5, -0.15) with radius 0.3, therefore the initial contact angle θ = 60. The viscosity is 0.001 in both fluids, the interfacial tension is 0.03. In Figure 4.9, we plot the evolution of the bubble at time T = 0,2₅,4₅,6₅,8₅, 2 .

4.10 Moving Contact Line (obtuse)

As 4.9, we compute the behavior of a bubble on the floor under the influence of the friction force at the contact line without gravity. The fluid is at a standstill initially. To set the static equilibrium contact angle θ₀ = 90. The interface is a part of circle. Here, we set the circle is centered at (0.5, 0.1) with radius 0.2, therefore the initial contact angle θ = 120. All other parameters are as 4.9. In Figure 4.10, we plot the evolution of the bubble at time T = 0,2₅,4₅,6₅,8₅, 2 .

5

Concluding Remarks

We have introduced the numerical scheme of how to use level set methods to describe interface and displayed many numerical experiments of two phase flows. These results show that the efficiency of level set methods describe interface with large topological change.

But, there are many details need to be improved and there are many problems what we can challenge. For example, the re-initialization of level set function is essential if we use level set methods. Because it uses iteration, we must cost much time. Therefore, it’s a good problem that how to increase our work efficiency.

The problem of surfactant is also very practical, many authors have researched this problem carefully. But most of them are based on front tracking where the free interfaces are explicitly tracked. It’s a challenge that how to solve the problem of surfactant using level set methods [10].

We need to solve the equation defined on the interface Γ 𝑓_t + u ∙ ∇𝑓 − 𝐧 ∙ (∇u𝐧)𝑓 = 𝐷_𝑠∇_𝑠2_𝑓

Where f is the consistency of surfactant, n is the normal vector to Γ directed towards outside, 𝐷_𝑠 is the surfactant diffusivity, ∇_𝑠= (𝐼 − 𝐧⨂𝐧)∇ is the surface gradient, u is velocity.

Since the level set method doesn’t like the front-tracking method which can control the information on interface completely. Therefore we should combine level set method with the other numerical methods or extend the equation to a band near the interface instead of restraining on the interface.

It’s a challenge to extend this method to three-dimension. But in order to describe physical phenomenon accurately, it’s a meaningful work we must do. As it should be, to maintain the efficiency and accuracy when extending this method into three-dimension is very important.

References

[1] Y. C. Chang, T.Y. Hou, B. Merriman, S. Osher, “A Level Set Formulation of Eulerian Interface Capturing Methods for Incompressible Fluid Flows.”Journal of Computational Physics 124, 449-464 (1996).

[2] Hauxiong Huang, Dong Liang, and Brian Wetton, “Computation of a Moving Drop/ Bubble on a Solid Surface Using a Front-Tracking Method.”August 23, 2004.

[3] Michael Renardy, Yuriko Renardy, and Jie Li, “Numerical Simulation of Moving Contact Line Problems Using a Volume-of-Fluid Method. ” Journal of Computational Physics 171, 243-263 (2001).

[4] S. Osher, R. Fedkiw, “Level Set Methods and Dynamic Implicit Surfaces.” Springer, 2003.

[5] A. Harten, B. Engquist, S. Osher, and S. Chakravarthy, “Uniformly High-order Accurate Essentially Non-Oscillatory Schemes III. ” J. Comput. Phys. 71, 231-303 (1987).

[6] C. W. Shu, and S. Osher, “ Efficient Implementation of Essentially Non Oscillatory Shock Capturing Schemes.”J. Comput. Phys. 77, 439-471 (1988). [7] C. W. Shu, and S. Osher, “ Efficient Implementation of Essentially Non

Oscillatory Shock Capturing Schemes II (two).”J. Comput. Phys. 83, 32-78 (1989).

[8] S. Osher, and J. Sethian, “Fronts Propagating with Curvature Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulations.” J. Comput. Phys. 79, 12-49 (1988).

[9] H. C. Tseng, “Numerical Methods and Applications for Immersed Interface Problems.”Master Thesis.

[10]J. J. Xu, Z. L. Li, J. Lowengrub, H. Zhao, “A Level-set Method for Interfacial Flows with Surfactant.”J. Comput. Phys. 212, 590-616 (2006).

[11]D. P. Peng, B. Merriman, S. Osher, H. K. Zhao, and M. J. Kang, “A PDE-Based Fast Local Level Set Method.” J. Comput. Phys. 155, 410-438 (1999).

[12]M. C. Lai, Y. H. Tseng, H. X. Huang, “An immersed boundary method for interfacial flows with insoluble surfactant.”submitted.