Numerical Scheme

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

Step1 Prediction:

We solve first step by Crank-Nicholson method in time

(∇ · (µⁿ⁺¹∇u^∗) − 2ρⁿ⁺¹²

Note that µⁿ⁺¹ and ρⁿ⁺¹² are computed by extrapolation as follows µⁿ⁺¹ = 2µⁿ− µⁿ⁻¹

ρⁿ⁺¹² = 3

2ρⁿ−1 2ρⁿ⁻¹

and Eq.(2) becomes Poisson type equation and can be solved by linear system for u^∗.

Step2 Projection:

Since the new u^∗ is not diverge free, we apply this step to correct u^∗ to be diverge free which means ∇ · u = 0 and then compute new pressure term ∇pⁿ⁺¹².By Helmholtz-Hodge decomposition, we have

u^∗ = uⁿ⁺¹+ ∆t∇ψⁿ⁺¹

∇ · uⁿ⁺¹ = 0 So we have

ρⁿ⁺¹²uⁿ⁺¹− u^∗

∆t = −∇ψⁿ⁺¹ and then take divergence we have

∇ · (∇ψⁿ⁺¹

ρⁿ⁺¹² ) = ∇ · u^∗

∆t After getting ∇ψⁿ⁺¹ we have

uⁿ⁺¹= u^∗− ∆t

ρⁿ⁺¹²∇ψⁿ⁺¹

∇pⁿ⁺¹² = ∇pⁿ⁻¹² + ∇ψ + ∇ · (µⁿ⁺¹(uⁿ⁺¹− u^∗)) Acknowledgment

The whole scheme of solving Navier-Stokes equation as above is provided by Yu-Hau Tseng and again we thank Tseng for providing the code of solving Navier-Stokes equations.

3 Solving PDE on Interface by Level Set Method

3.1 Introduction

In the numerical simulation of partial differential equations, we usually en-counter the problem of handling interface, for example:two-phase flows of Navier-Stokes equation or heat equation on the interface which is irregular or high dimensional. To describe the interface is the main issue when solving PDE with interface problem. We use level set function which we will have details later. There are many advantages that level set method provides, for example the interface can be easy presented and when normal vector on the interface are needed, it also can be calculated easily by level set method.

Our goal in this chapter is to solve heat equation on the circle which means that the heat quantities just diffuse on the circle and do not diffuse out of circle.As we just use Cartesian grid computational domain, how do we handle this problem? We will give details step by step.

3.2 Level Set Function and Signed Distance Function

3.2.1 Level Set Function

The concept of level set function is very simple. We let the symbol φ to represent level set function.On the interface, φ is defined to be zero which is zero level set, outside the interface we define φ > 0 and inside the interface we let φ < 0. With the help of level set function, we can easily distinguish points which are inside or outside the interface.(See figure1)

3.2.2 Signed Distance Function

Next we will explain what the signed distance function is. Let a point P (x₀, y₀) which is arbitrary and we have level function φ .If the value of φ(P ) is equal to the distance from P to the interface,we call that φ is signed distance function. Maybe this description is not clear enough, let us take a example

Example.1 We take a level function of circle, (See figure1) φ(x, y) =p

x²+ y²− 1

with center O(0, 0) and radius r = 1 and we claim that φ is signed dis-tance function. Note that O is inside interface and disdis-tance from O to inter-face is 1, then we find that φ(0, 0) = −1 where − represent inside interinter-face

Figure 1: (a)Definition of level function φ (b)Figure of example.1 and 1 denote the distance. A(1, 0) is on the interface and we find that φ(1, 0) = 0 and B(1.5, 0) is outside interface with distance 0.5 to interface, so φ(1.5, 0) = 0.5.Since φ is signed distance function, we have such pretty property.

From the above example we can figure out the meaning of signed distance function. So for any point plugging into φ, we can know the point is out-side,inside or on the interface and the distance from point to interface is also known when φ is signed distance function.

3.2.3 How To Distinguish Signed Distance Function?

In few cases we can just write down the signed distance function as math-ematical formulation exactly like example.1. Since circle have very good properties, so the above example of circle is special case. Note that the level function of unit circle is not only one case. Consider that φ(x, y) = x²+y²−1 is also a level function for unit circle, but unfortunately,it is no longer a signed distance function. Let us take a test,plug B into φ we get φ(1.5, 0) = 1.25 is not equal to 0.5, so this level function is not signed distance function.

Actually, signed distance function has a very good property which is

|∇φ| = 1 and it is also a condition to check whether a level function is a signed distance function or not. For computing normal n = _|∇φ|^∇φ, if φ is signed distance function, we can calculate n more accurately and n = ∇φ where ∇φ can be calculated by central difference and get second order for n.

If φ is not signed distance function then n = _|∇φ|^∇φ will be first order accuracy and may not accurate enough, this is why using signed distance function to

calculate n will be better.

Most of level set functions are not signed distance function. For example, the level function of ellipse:φ(x, y) =p(^x_a)²+ (^y_b)²−1 is not a signed distance function.It is easy to check that |∇φ| 6= 1. In order to let readers understand more about signed distance function and distinguish whether level function is signed distance function or not, we contour three level functions in the following : function for ellipse. In Eq.(1),the distance between any two adjacent level curves are the same and equal to 5∆x, that is why we call Eq.(1) is signed distance function. Eq.(2) and Eq.(3) do not have such nice property and the distance between any two adjacent level curves are not the same. So when calculating normal on interface, signed distance function will be more accurate. Now we have introduced another method to distinguish signed distance function which is to observe the contour of φ.

Remark

The fifth level curve of each φ in figure2 is the interface, that is φ = 0.We will find that the location of φ = 0 of Eq.(1) and Eq.(2) are at the same location since they are both level functions for unit circle and φ = 0 denotes the interface.

Since most level functions are not signed distance function and the exact mathematical form of signed distance function is hard to find, so how do we handle this problem? In the next section, we will introduce the Re-initialization process which can overcome the problem.

3.3 Re-initialization

In this section, we introduce how to modify a level function which is not a signed distance function into a signed distance function. The process is not complicated and it involves solving a PDE to steady state as follow:

φ_t+ S(φ₀)(|∇φ| − 1) = 0 (6)

Figure 2: Contour of Eq.(3),(4) and (5)

where

Let us explain the meaning of Eq.(6) first. Consider a PDE which similar to Eq.(6) :

φ_t+ |∇φ| = f (x) (7)

If we solve Eq.(7) to steady state, we will get that |∇φ| = f (x). Since signed distance function has a property |∇φ| = 1, so if the level function φ is not signed distance function(i.e. |∇φ| 6= 1), we use the concept of Eq.(7) to enforce |∇φ| = 1, that is to solve :

φ_t+ |∇φ| = 1 (8)

From Eq.(7) we know that the information is propagated in the normal di-rection. So the information is carried from small φ-value to large φ-value.

So we hope the information propagated from interface(i.e. φ = 0). But if we solve Eq.(8) for whole domain Ω, since the level function is negative inside interface,so the information is propagated from the Ω⁻ and obviously it is wrong.

So if we solve Eq.(8) in Ω⁺ ∪ ∂Ω then it means that the information is propagated form interface to Ω⁺and it make sense. Similar concept, we want the information propagated from interface to Ω⁻, so we solve the equation :

φ_t− |∇φ| = −1 (9)

in the domain Ω⁻∪ ∂Ω.

Finally we couple Eq.(8) and Eq.(9) to get Eq.(6). Solving Eq.(6) to steady state, we have |∇φ| − 1 = 0(i.e. |∇φ| = 1 ) and that is the goal of reinitialization. Note that if φ = 0 denotes the interface, then after reinitial-ization process,φ = 0 is static and do not change it’s value since S(φ0) = 0 on ∂Ω. So the reinitialization process ensure that the location of interface will not be changed.

3.4 Spatial Discretization for Reinitialization Process

In this section, we will describe how to handle |∇φ| in Eq.(6). If we just use central difference to compute |∇φ| , then we will get terrible result.Before handling this term, let us look at 1D example :

Example.2 (Burger’s equation)

u_t+ uu_x = 0 (10)

We discretize u_x in Eq.(10) by upwind difference.We denote : D_i⁺u = ui+1− ui

∆x and D_i⁻u = ui− ui−1

∆x

Note that the coefficient u of u_x decides the information from which direc-tion. That is, if u > 0 we know that the characteristics come from left side so we use D⁻_i u to compute (u_i)_x. On the other hand, if u < 0,it means that the information and characteristics come from right side so we use D⁺_i u to compute (u_i)_x.

The above equations explain what is upwind difference for 1D case.Note that Eq.(6) and Eq.(10) are similar for some sense,but Eq.(6) is a 2D PDE and is not as simple as Eq.(10) when handling first derivative term.

3.4.1 Godunov’s Method

Since Eq.(6) is a 2D PDE which is not as simple as Eq.(10), we introduce Godunov’s Method to handle |∇φ| in Eq.(6). Actually, Godunov’s method uses the same concept which is characteristic coming from which direction.

But in 2D domain you may doubt that the characteristic can come from infinite directions and how to decide? To overcome this problem, we just simplify and separate x-direction and y-direction to compute respectively that is just compute φ_x and φ_y respectively. Since we know the information propagated from interface, so when computing (φ_i,j)_x, we choose the point of φ_i+1,j and φ_i−1,j by the direction of information propagating and then compute (φ_i,j)_x. For example,if (i, j) is outside interface and if φ_i−1,j <

φ_i,j < φ_i+1,j which means (i − 1, j) is nearest interface within these three points and the information is propagated from φ_i−1,j to φ_i+1,j, so we use D⁻_i,jφ = ^(φi,j−φ_∆x^i−1,j). On the other hand, if φi−1,j > φi,j > φi+1,j which means information comes from φ_i+1,j to φ_i−1,j, then we use D⁺_i,jφ = ^(φi+1,j_∆x^−φi,j).

Having the basic concept of Godunov’s mothod, we find that in some sense this method is as upwind method. Here we describe Godunov’s method for solving Eq.(6) in the following :

Details of Godunov’s Method

1. If φ_i,j > 0 which means the point is outside interface, then

(a) If D_i,j⁺φ > 0 and D_i,j⁻φ > 0 then it means φ_i−1,j < φ_i,j < φ_i+1,j and (φ_i,j)_x = D_i,j⁻φ

(b) If D_i,j⁺φ < 0 and D_i,j⁻φ < 0 then it means φ_i−1,j > φ_i,j > φ_i+1,j and (φi,j)x = D_i,j⁺φ

(c) If D_i,j⁺φ ≤ 0 and D_i,j⁻φ ≥ 0 then it means φ_i−1,j ≤ φ_i,j and φ_i,j ≥ φ_i+1,j. In this case, the information comes form both two sides and we treat the shock by the larger influence between D⁻_i,j and D_i,j⁺, that is (φ_i,j)_x = max(|D_i,j⁻|, |D_i,j⁺|)

(d) If D_i,j⁺φ ≥ 0 and D_i,j⁻φ ≤ 0 then it means φ_i−1,j ≥ φ_i,j and φ_i,j ≤ φi+1,j and the information flow to both sides from (i, j). This does not make sense since interface do not pass grid (i, j). So in this case we just set (φ_i,j)_x = 0.

2. If φ_i,j < 0 which means the point is inside interface.

(Note that inside interface the φ-value is negative and if φ is larger, it is much near the interface.)

(a) If D_i,j⁺φ > 0 and D_i,j⁻φ > 0 then it means φ_i−1,j < φ_i,j < φ_i+1,j and (φ_i,j)_x = D_i,j⁺φ.

(b) If D_i,j⁺φ < 0 and D_i,j⁻φ < 0 then it means φi−1,j > φi,j > φi+1,j and (φ_i,j)_x = D_i,j⁻φ

(c) If D_i,j⁺φ ≥ 0 and D_i,j⁻φ ≤ 0 then it means φ_i−1,j ≥ φ_i,j and φ_i,j ≤ φi+1,j. In this case, the information comes form both two sides and we treat the shock by the larger influence between D⁻_i,j and D_i,j⁺, that is (φ_i,j)_x = max(|D_i,j⁻|, |D_i,j⁺|)

(d) If D_i,j⁺φ ≤ 0 and D_i,j⁻φ ≥ 0 then it means φi−1,j ≤ φi,j and φi,j ≥ φ_i+1,j and the information flow to both sides from (i, j). This does not make sense since interface do not pass grid (i, j). So in this case we just set (φi,j)x = 0.

After simplifying the Godunov’s method as above, the scheme is obvious and easy to implement by numerical computing. Here we give three remarks as follows:

Remark

1. In Godunov’s method, we discetize φ_x by one-sided difference, so the scheme is first order accurate.

2. The cases φ > 0 and φ < 0 are opposite since the signed function S(φ₀) is 1 outside interface and −1 inside interface.

3. The y-direction case φ_y is as same as computing φ_x. So after solving φ_x and φ_y, the spatial discretization of |∇φ| in Eq.(6) is complete.

3.4.2 Second-Order Hamilton-Jacobi ENO Scheme

Since the Godunov’s method is first order, if we want the scheme more pre-cise, we introduce second-order ENO (essential non-oscillatory) scheme to discretize. The ENO scheme is based on Godunov’s method for first-order term which described above when solving Eq.(6). We write the Taylor ex-pansion of φ:

φ(x) = φ(˜x) + φ⁰(˜x)(x − ˜x) + φ⁰⁰(˜x)

2! (x − ˜x)²+ · · · (11) If we use D⁻φ_i, then it means we plug x_i and x_i−1 into Eq.(11) and get (φ_i)_x. Now the ENO scheme is to add a adjacent point plug into Eq.(11).

The point is chosen to be much smooth between two sides. Let us take a example to explain :

Example.3 If Godunov’s method tell us that we choose the point x_i−1 when computing (φi)x then we need to choose another point between xi+1

and x_i−2 to plug into Eq.(11) and get second-order result. We want choose much smooth side point so we need some criterion for choosing point. Let us denote: and we choose the point xi−2to plug into Eq.(11). Otherwise we choose xi+1

to plug into Eq.(11). No matter what point is chosen, the scheme is always second-order.

Now we describe the whole process for 1D case for ENO scheme based on Godunov’s method :

1. If we use D⁻φ for first-order discretization to compute (φ_i)_x

Figure 3: (a).left-biased (b).right-biased

To handle 2D case, it is as the same method and we split the domain into x-direction and y-direction to handle respectively. Here we have intro-duced the second-order ENO scheme for solving Hamilton-Jacobi equation like Eq.(6). Next we will introduce more higher order scheme.

3.4.3 Third-Order ENO Scheme

The third-order ENO scheme for spatial is based on second-order ENO scheme which is based on Godunov’s method, so we realize that the concept of Godunov method is very important when solving Eq.(6). In second-order ENO, we use three points to calculate φ_x and get second-order result, the third-order is as the same concept as ENO-2. While deciding which three points are needed, we still demand one point to reach third order ENO(ENO-3) which is either left side or right side adjacent point. The fourth point is also chosen in smoother side and the criterion is similar to above ENO-2.

Since the concept is similar, we do not take example again and we just give the details of algorithm as below:

1. If we use D⁻φ for first-order discretization to compute (φ_i)_x

2. If we use D⁺φ for first-order discretization to compute (φ_i)_x

It seems that there are eight cases for ENO-3, but actually it can be reduced to just three cases. That is using {x_k|k = i − 3, i − 2, i − 1, i} or {x_k|k = i − 2, i − 1, i, i + 1} {x_k|k = i − 1, i, i + 1, i + 2} for D⁻ case to plug into Eq.(11). Similarly the right-biased stencil case for D⁺ is {x_k|k = i, i+1, i+2, i+3} or {x_k|k = i−1, i, i+1, i+2} or {x_k|k = i−2, i−1, i, i+1}.

Finally we get just three discretization forms.

3.4.4 Third-Order WENO

In this section we will discuss the third-order WENO(weighted essential non-oscillatory) scheme to discretize φ_x. Again W3 is based on ENO-2.ENO-2 has two possibilities to compute φ_x under D⁻ (See Eq.(13) and

Eq.(15))and we combine this two by weights instead using only one to com-pute φ_x. That is why we call this method Weighted ENO.

The scheme of WENO-3 is more easy than ENO-3 since we only need to judge D⁻ or D⁺. If we need D⁻ for first-order then we use left-biased points {x_k|k = i − 2, i − 1, i, i + 1} to compute φ_x. On the other case we use right-biased points {x_k|k = i − 1, i, i + 1, i + 2} to compute φ_x. The process of WENO is more simple since we do not need to judge which points are needed to add.

We write down the details of WENO-3 which we refer [12]:

1. If D⁻φ is needed then

Note that the form is more simple and is easy to implement, so we rec-ommend the WENO-3 scheme to discretize the φ_x and φ_y of reinitialization PDE.

Remark

Although the sign function of Eq.(6) is just simple 0,1 or −1, [12] suggested that

S(φ) = φ

pφ²+ |∇φ|²(∆x)² (20)

for a better result and reduce the time steps when solving Eq.(6) to steady state.

3.5 Temperal Discretization for Reinitialization Pro-cess

Since we prefer to use WENO-3 to discretize |∇φ| for spatial, if we want the scheme is higher order, we also need third-order method in time. Here we use TVD Runge-Kutta method which we refer [1].

3.5.1 Second-Order TVD Runge-Kutta Method

We let current time step be n and first we move time to n + 2,then take average as below

φⁿ⁺¹= 1

2φⁿ+ 1 2φⁿ⁺² 3.5.2 Third-Order TVD Runge-Kutta Method

For third-order TVD-RK scheme we first need φⁿ⁺¹² which is weighted by φⁿ and φⁿ⁺² and so that

φⁿ⁺¹² = 3

4φⁿ+ 1 4φⁿ⁺²

Next we move φⁿ⁺¹² to φⁿ⁺³² by iteration and finally we get φⁿ⁺¹ as follows:

φⁿ⁺¹ = 1

3φⁿ+2 3φⁿ⁺³²

Note that if ENO-2 is used for spatial, then TVD-RK2 have better been used for time and both ENO-3 and WENO-3 match TVD-RK3.

3.6 Numerical Test for Reinitialization Process

In this section, we will test some examples for reinitialization. We choose above example φ(x, y) = p(_0.6^x )²+ (_0.3^y )² − 1 and a concave case which we refer [12] with polar coordinate φ(r, θ) = (r − 0.5 + 0.1r sin(7θ))³.We contour [−10∆x : 2∆x : 10∆x] for each case. Note that we use WENO-3 and TVD-RK3 method for these examples.

Figure 4: Reinitialization of ellipse

Figure 5: Reinitialization of plum blossom shape

Remark

1. From the figures of two examples, we find that the level function be-comes signed distance function from interface to Ω⁺ and Ω⁻ within a neighborhood of interface. As the time steps increase, the width of neighborhood also increase. The reason is that the information flows from interface to Ω⁺ and Ω⁻.

2. After reinitialization process, the level function becomes signed distance function in a neighborhood of interface and we can not write the level function by mathematical formulation exactly.

3. There is no clear condition of stopping criterion for reinitialization pro-cess. The iteration steps is determined by the width of signed distance function neighborhood of interface. If you want the width wider, the iteration steps needed more.

4. Eq.(6) is solved in whole domain and the above two examples are lo-cated at center region of computation domain, so if we just want a signed distance function of just a neighborhood of interface, we do not need the boundary condition.But if the interface contacts with the boundary, how do we set the boundary condition for Eq.(6)? We will give the idea and algorithm in next section.

3.7 Special Case and Boundary Condition

To handle Eq.(6) with interface contacts with boundary,we must be very carefully on boundary condition. Since the ENO and WENO scheme may need more than two points to compute, so on the boundary we have better not use ENO or WENO scheme, instead we use first-order Godunov’s method to handle since this method just need two points to compute the first derivative term.Here we give the idea of how to set boundary condition for 1D case : Concept of Boundary

If the left end point x₁ is in Ω⁺ and we want compute φ_x. A ghost point x₀ on the left side of x₁ is needed.Since x₁ is in Ω⁺ which means interface is on the right side of x₁ and x₀ is on the left side of x₁, so φ(x₀) ≥ φ(x₁).

Similarly, if x1 is in Ω⁻ then φ(x0) ≤ φ(x1)

Although the point x₀ does not exist in the domain Ω, but the role of this point is important. We will give code of algorithm as below

Figure 6: Contour of special case example Algorithm

For x-direction and left side boundary do j=1,n

if (phi(1,j)>0d0) then

if (phi(1,j)>=phi(2,j)) then DX(1,j)=(phi(2,j)-phi(1,j))/h else

DX(1,j)=0d0 end if

else

if (phi(1,j)<=phi(2,j)) then DX(1,j)=(phi(2,j)-phi(1,j))/h else

DX(1,j)=0d0 end if

end if end do

where DX denotes φ_x and the three other boundary are as same method to handle. We also suggest that the points adjacent boundary also use

Go-dunov’s method, for 1D case we mean x₂. Here we contour a ellipse

for [−10∆x : 2∆x : 10∆x] before and after reinitialization with boundary condition.(See figure6)

Until here, we have presented all process of reinitialization.No matter the interface is concave or convex, the reinitialization process all works that is if we can write down the level function of interface, we can reinitialize it to signed distance function even in the whole domain. The reinitialization process is very useful since it can calculate normal or curvature accurately, so it is really a very good and valuable scheme.

3.8 Normal Extension

In this section we will discuss the normal extension process. If u is some quantity which is defined on the interface, we can extend u in the normal direction.That is along the normal direction, the quantity u is constant. To

In this section we will discuss the normal extension process. If u is some quantity which is defined on the interface, we can extend u in the normal direction.That is along the normal direction, the quantity u is constant. To

在文檔中 等位集與流體體積混合應用在二相流 (頁 9-0)