在流體體積法中使用高解析離散法計算兩相流

國 立 交 通 大 學

機 械 工 程 學 系

碩 士 論 文

在流體體積法中使用高解析離散法計算兩相流

Use of high-resolution scheme in volume-of-fluid method

for two-fluid flow

研 究 生：林 仕 文

指 導 教 授：崔 燕 勇 博士

在流體體積法中使用高解析離散法計算兩相流

Use of high-resolution scheme in volume-of-fluid method for two-fluid flow

研 究 生：林 仕 文 Student：Shi-Wen Lin

指導教授：崔 燕 勇 Advisor：Yeng-Yung Tsui

國立交通大學

機械工程研究所

碩士論文

A Thesis

Submitted to Institute of Mechanical Engineering

Collage of Engineering

National Chiao Tung University

In Partial Fulfillment of the Requirements

For the degree of

Master of Science

In

Mechanical Engineering

July 2008

Acknowledgments

I would like to thank Prof. Yeng-Youg Tsui for his constructive criticism, guidance provided throughout my research, patience and continuous support.

I would like to thank all individuals of CFD Laboratory.

Thanks to the seniors, Yu-Chang Hu, Tong-Ting Chen, and Tian-Cherng Wu for their teaching and helping when I have difficulties in the investigate.

Thanks to Xin-En Wu and Jun-Yan Lee for their friendships and companions. They give me great help during the two years.

Then, I would like to thank all members of my family for their attention during the two years.

Finally, I would like to thank my parents, dad Lin and mom Zhong for their support and encouragement. They let me know how to deal with problems and overcome the difficulties.

Use of high-resolution scheme in volume-of-fluid method for two-fluid flow

Student：Shi-Wen

Lin

Advisor：Prof. Yeng-Yung Tsui

Institute of Mechanical Engineering National Chiao Tung University

ABSTRACT

A numerical method for direct simulations of two-fluid flows is established in this study. The motion of the interface is captured by the solution of a transport equation for the volume fraction. Some numerical schemes, such as high resolution and compressive schemes are discussed in this study. The high resolution schemes preserve the shape of the interface but can not reduce the numerical diffusion. The compressive schemes are able to reach less numerical diffusion but let the interface deformed. Most composite schemes switch the compressive scheme and high resolution scheme with a switching function about the slope of the interface in order to overcome the above problem. The aim of this study is to develop a composite of the modified MUSCL and modified bounded downwind scheme with a switching function. This scheme presents the high accurate results in test cases and can be used on the quadrilateral and triangular mesh.

在流體體積法中使用高解析離散法計算兩相流

研究生 林仕文 指導教授 崔燕勇 博士 國立交通大學機械工程學系 摘要 在本研究中將建立ㄧ種直接模擬兩相流的數值方法。在此方法中，藉由求解流體體 積分率的傳輸方程式，來抓取兩流體間的介面運動。在本研究中將討論許多數值的離散 法，像是高解析離散法和壓縮式離散法。高解析離散法可以維持兩流體間界面的外形， 卻沒有辦法減少在介面上的數值擴散。而壓縮式離散法可以達到較少的數值擴散，卻會 造成介面形狀的變形。為了處裡上述的問題，大部分混合離散法都是採用一個與介面斜 率相關的轉換函數，在高解析離散法和壓縮式離散法之間做切換。本研究主要的目的是 在於發展一種在 modified MUSCL 和 modified bounded downwind scheme 之間轉換的 混合離散法。此方法在本文中所測試的案例中都有相當高的準確性，而且本方法可以在 三角形和四邊形的網格系統中使用。

Content

ABSTRACT ...iii

Content ... v

List of Table... vii

List of Figure...viii

Nomenclature... xii

Abbreviation ... xiv

Chapter1

Introduction... 1

1.1 Background ... 1 1.2 Related studies ... 2 1.2.1 Lagranian schemes ... 2 1.2.2 Eulerian schemes ... 2

1.3 Outline of this thesis ... 7

Chapter2

Mathematical Model... 9

2.1 Introduction ... 9

2.2 General transport equation ... 9

2.3 Mass and momentum conservation equation... 10

2.4 VOF equation ... 10

2.5 Surface tension... 12

2.6 Boundary conditions ... 13

Chapter3

Numerical Method... 14

3.1 Introduction ... 14

3.2 Discretization of the VOF equation... 14

3.3 Discretization of the momentum equation ... 15

3.3.1 Unsteady term ... 15

3.3.2 Convection term... 16

3.3.3 Diffusion term ... 16

3.3.4 Source term ... 17

3.3.5 Arrangement of the difference transport equation ... 18

3.4 Pressure-velocity coupling of the PISO algorithm ... 18

3.5 Solution procedure... 24

Chapter4

High Resolution Differencing Schemes ... 25

4.1 Introduction ... 25

4.2 Convective flux of volume fraction ... 25

4.3 CBC and TVD constraints ... 27

4.4 Linear and non-linear schemes... 27

4.5 Composite scheme with switching function ... 28

5.2 Uniform density flow ... 32

5.3 Shear flow ... 34

5.4 Broken dam ... 36

5.5 Filling process in an open tank ... 38

Chapter6

Conclusion and Future Work ... 40

List of Table

Table 4.1 The NVD equation and flux limiter function of linear and non-linear

difference schemes ... 47

Table 5.1 Errors of different schemes in

Co=0.25

(quadrilateral mesh) ... 48

Table 5.2 Errors of different schemes in

Co=0.75

(quadrilateral mesh) ... 48

Table 5.3 Errors of different schemes in

Co=0.25

(triangular mesh)... 48

Table 5.4 Errors of different schemes in

Co=0.75

(triangular mesh)... 49

Table 5.5Errors of different composite schemes in

Co=0.25

(quadrilateral mesh)

... 49

Table 5.6Errors of different composite schemes in

Co=0.75

(quadrilateral mesh)

... 49

Table 5.7Errors of different composite schemes in

Co=0.25

(triangular mesh)49

List of Figure

Figure 1.1 The method of two-fluid flow (a) Lagrangian (b) Eulerian scheme . 50

Figure 1.2 Front tracking method... 51

Figure 1.3 Maker and cell method ... 51

Figure 1.4 Line techniques... 52

Figure 1.5 Donor and accepter cell configuration... 53

Figure 2.1General form of the conservation law ... 54

Figure 2.2 VOF method on the Eulerian grids... 54

Figure 2.3 Continuity of the velocity and discontinuity of the momentum... 55

Figure 2.4 Fluid arrangements and the sign of the curvature ... 55

Figure 3.1 Illustration of the primary cell P and the neighbor cell nb with a

considering face... 56

Figure 4.1 The relationship of a control volume and its neighbor cells ... 57

Figure 4.2 The CBC constraint in NVD... 57

Figure 4.3 The CBC constraint in NVD... 58

Figure 4.4 The TVD condition in TVD diagram ... 58

Figure4.5 Linear schemes in Normalized Variable Diagram... 59

Figure4.6 Linear schemes in TVD Diagram ... 59

Figure 4.7 SMART... 60

Figure 4.8 MUSCL... 60

Figure 4.9 SUPERBEE ... 60

Figure 4.10 STOIC ... 61

Figure 4.11 OSHER ... 61

Figure 4.12 BDS... 61

Figure 4.13 Van Leer ... 62

Figure 4.14 CHARM... 62

Figure 4.15 Modified BDS... 62

Figure 4.16 Modified MUSCL... 63

Figure 4.17 The NVD of (a) HYPER-C and (b) UQ (

Co

= 0.5)... 63

Figure 4.18 The switching function of CICSAM scheme ... 64

Figure 4.19 The switching function of HRIC scheme ... 64

Figure 4.20 The switching function of the composite of M-MUSCL and M-BDS

... 64

Figure 5.1 (a) The constant velocity field and initial position of hollow circle . 65

Figure 5.1 (b) The constant velocity field and initial position of hollow square 65

Figure 5.2 Triangular computational mesh in uniform density flow (22478 cells)

... 66

Figure 5.3 (a) The exact distribution of the hollow circle ... 67

Figure 5.3 (b) The exact distribution of the hollow square... 67

Figure 5.4 The final shape of hollow circle from different schemes in quadrilateral

mesh... 68

Figure 5.5 The final shape of hollow square from different schemes in

quadrilateral mesh ... 69

Figure 5.6 The final shape of hollow circle from different schemes in triangular

mesh... 70

Figure 5.7 The final shape of hollow square from different schemes in triangular

mesh... 71

Figure 5.8 The final shape of hollow circle from different composite schemes in

quadrilateral mesh ... 72

Figure 5.9 The final shape of hollow square from different composite scheme in

quadrilateral mesh ... 73

Figure 5.10 The final shape of hollow circle from different composite scheme in

triangular mesh... 74

Figure 5.11 The final shape of hollow square from different composite scheme in

triangular mesh... 75

Figure 5.12 The volume fraction distribution in a shear flow field ... 76

Figure 5.13 Triangular mesh in uniform density flow (22494 cells) ... 76

Figure 5.14 The volume fraction distribution in shear flow with Co=0.25

(quadrilateral mesh)... 77

Figure 5.15 The volume fraction distribution in shear flow with Co=0.25

(quadrilateral mesh)... 78

Figure 5.16 The volume fraction distribution in shear flow with Co=0.25

(quadrilateral mesh)... 79

Figure 5.17 The volume fraction distribution in shear flow with Co=0.25

(quadrilateral mesh)... 80

Figure 5.18 The volume fraction distribution in shear flow with Co=0.75

(quadrilateral mesh)... 81

Figure 5.19 The volume fraction distribution in shear flow with Co=0.75

(quadrilateral mesh)... 82

Figure 5.20 The volume fraction distribution in shear flow with Co=0.75

(quadrilateral mesh)... 83

Figure 5.21 The volume fraction distribution in shear flow with Co=0.75

(quadrilateral mesh)... 84

Figure 5.22 The volume fraction distribution in shear flow with Co=0.25

Figure 5.23 The volume fraction distribution in shear flow with Co=0.25

(triangular mesh) ... 86

Figure 5.24 The volume fraction distribution in shear flow with Co=0.25

(triangular mesh) ... 87

Figure 5.25 The volume fraction distribution in shear flow with Co=0.25

(triangular mesh) ... 88

Figure 5.26 The volume fraction distribution in shear flow with Co=0.75

(triangular mesh) ... 89

Figure 5.27 The volume fraction distribution in shear flow with Co=0.75

(triangular mesh) ... 90

Figure 5.28 The volume fraction distribution in shear flow with Co=0.75

(triangular mesh) ... 91

Figure 5.29 The volume fraction distribution in shear flow with Co=0.75

(triangular mesh) ... 92

Figure 5.30 Comparison of errors in the shear flow with Co=0.25 (quadrilateral

mesh) ... 93

Figure 5.31 Comparison of errors in the shear flow with Co=0.75 (quadrilateral

mesh) ... 93

Figure 5.32 Comparison of errors in the shear flow with Co=0.25 (triangular mesh)

... 94

Figure 5.33 Comparison of errors in the shear flow with Co=0.75 (triangular mesh)

... 94

Figure 5.34 Schematic of the broken dam ... 95

Figure 5.35 Experimental results of a collapsing water column by Koshizuka . 96

Figure 5.36 Schematic representation of the non-uniform and quadrilateral mesh

with

56 36×

grids in broken dam ... 97

Figure 5.37 Schematic representation of the triangular mesh with 4506 cells in

broken dam ... 97

Figure 5.38 Schematic representation of the triangular mesh with 12354 cells in

broken dam ... 97

Figure 5.39 Numerical results of the broken dam on the uniform and quadrilateral

mesh with

48 28×

grids... 98

Figure 5.40 Numerical results of the broken dam on the uniform and quadrilateral

mesh with

120 70×

grids ... 99

Figure 5.41 Numerical results of the broken dam on the non-uniform and

quadrilateral mesh with

56 36×

grids ... 100

Figure 5.42 Numerical results of the broken dam on the triangular mesh with 4506

cells... 101

Figure 5.43 Numerical results of the broken dam on the triangular mesh with

12354 cells... 102

Figure 5.44 The position of the leading edge in broken dam ... 103

Figure 5.45 The height of the collapsing water in broken dam ... 103

Figure 5.46 Schematic representation of the filling process in an open tank ... 104

Figure 5.47 Schematic representation of the triangular mesh with 902 cells in

filling process ... 104

Figure 5.48 Schematic representation of the triangular mesh with 2024 cells in

filling process ... 105

Figure 5.49 Schematic representation of the triangular mesh with 3584 cells in

filling process ... 105

Figure 5.50 The volume fraction distribution and velocity field of the filling

process on the uniform and quadrilateral mesh with

28 28×

... 106

Figure 5.51 The volume fraction distribution and velocity field of the filling

process on the uniform and quadrilateral mesh with

40 40×

grids ... 107

Figure 5.52 The volume fraction distribution and velocity field of the filling

process on the uniform and quadrilateral mesh with

80 80×

grids... 108

Figure 5.53 The volume fraction distribution and velocity field of the filling

process on the triangular mesh with 902 cells ... 109

Figure 5.54 The volume fraction distribution and velocity field of the filling

process on the triangular mesh with 2024 cells ... 110

Figure 5.55 The volume fraction distribution and velocity field of the filling

process on the triangular mesh with 3584 cells ... 111

Figure 5.56 The position of leading of the filling process in the open tank ... 112

Figure 5.57 The water volume inside the tank... 112

Nomenclature

A Coefficient of algebraic equation

Co Courant number

C

F Boundary flux due to convection

D

F Boundary flux due to diffusion

f _{Face, point in the centre of the face}

σ

f Surface tension force

g _{Gravitational acceleration}

f

m& _{Mass flow rate cross face}

P Pressure

Q_{φ Source term in discretized momentum equation}

α

S _{Source term in discretized the indicator equation}

t Time

f

w _{Weighting factor}

y

x, _{Components of Cartesian coordinate}

d

δr Vector pointing from center of primary cell to the center of neighbor one

φ Property

∀

Δ Volume

) (r

γ Flux limiter equation

α _{Volume fraction}

σ _{Surface tension coefficient}

ρ _{Fluid density}

Δ Difference operator

∇ Gradient

f

θ Angle between the normal to interface and orthogonal component of the

face area vector ( )_f

ψ θ _{Switching function}

n _{New time level}

Abbreviation

BDS Bounded Downwind Scheme

CBC Convective Boundedness Criterion

CDS Central Difference Scheme

CICSAM Compressive Interface Capturing Scheme for Arbitrary Meshes

CUS Cubic upwind difference scheme

CV Control Volume

DDS Downwind Difference Scheme

HRS High Resolution Scheme

LUS Linear Upwind Scheme

MAC Marker and cell

MUSCL Monotonic Upwind Scheme for Conservation Law

NVD Normalized Variable Diagram

PISO Pressure Implicit with Splitting of Operator

QUICK Quadratic Upwind Interpolation for Convection Kinematics

SLIC Simple Line Interface Calculation

SMART Sharp and Monotonic Algorithm for Realistic Transport STOIC Second- and Third-Order Interpolation for Convection

TVD Total Variation Diminishing

UDS Upwind Difference Scheme

UQ ULTIMATE-QUICKEST

M-MUSCL Modified Monotonic Upwind Scheme for Conservation Law

Chapter1 Introduction

1.1 Background

The flow involving two immiscible fluids has been of interest to many investigators during the last decades. How to predict the position and the movement of the interfaces in the two-fluid flow accurately is very important in many scientific and technical applications. The objective of this study is the development of a numerical method which can cope with the above problem. In many engineering problems and industrial processes, such as marine engineering, biochemical engineering, tube/channel flows, and casting, welding, molding, injection or extrusion processes, the simulation of two-phase flow with discrete interface is a rather popular issue. This simulation has also played an important part in IC package process and the production of LED screen. In marine engineering, some numerical results of free surface flow with wave breaking are applied to the motion of sea water. Muzaferija et al. have presented the flow around ship hulls or submerged hydrofoils [1]. In biochemical engineering, this technique is used to simulate the transportation of biochemical fluids in the capillary channels, such as the blood in veins. On the other areas, mold-filling process with heat transfer is an important application in casting. K.A. Percleous et al. investigate the collapse of a liquid column in a sealed cavity and simulate cooling process of a step-like model which is a common test one [2]. In the tube/channel flows area, the numerical results of two-phase flow are widely employed. In [3], Yang et al. present the flow boiling of refrigerant R134B in a horizontal coiled tube. They predict the temperature profile and the phenomenon of the flow boiling in the straight and bending parts of the coiled tube. In [4], the simulation of two-phase flow can also be used in the fuel cell. Yun Wang et al. [4] have developed a two-phase model for the flow in mini-channels of a proton exchange membrane (PEM) fuel cell. The results of their investigations can be applicable to many common two-phase flow behaviors across the micro- or mini- channels. Although the application of the two-phase flow is so extensive, the representative problems of the flow with two immiscible fluids can be classified into three

categories [5]. The first one is dispersed flow that the two fluids in the flow field are considered as suspensions without a defined interface. The second is that the two fluids are separated by a sharp interface without breaking. The third is the transitional flow that the interface of two-phase flow may or may not be broken. In this study, the numerical method mainly deals with the last two problems.

1.2 Related studies

Several numerical methods of two-fluid flow with a moving interface have been posed. The methods which are used to predict the phenomenon of two-phase flow can be divided into two main categories: Lagrangian and Eulerian schemes (Fig. 1.1) [6].

1.2.1 Lagranian schemes

The first method can keep free surface sharp between the two fluids and present the exact position of the free surface with re-meshing as the calculation proceeds. The mesh of this scheme is deformed and changed all the time. The main procedure of this method is that the position of the free surface at next time step is calculated by using the velocity field which is known. When the new free surface boundary is defined, we can reconstruct grids and update new properties for the new flow field [7]. The position of the interface can be predicted precisely, because the boundary mesh matches the free surface. Although the accurate prediction of the free surface can be carried out by the Lagrangian scheme, this method can not be employed to the quite complex flow field. Many deformations and stretches which result from the breaking, overturning, or gravity wave may cause numerical errors as well as reducing the precision. In [8], there is a Langranian scheme which is presented by Peric et al. This scheme can be employed in some simple and no large interface deforming cases, but its drawback is that the scheme cannot be used while the interface is unfolded.

1.2.2 Eulerian schemes

The second method can reduce errors which result from the deforming grids by using fixed grids that are generated before calculating the movement of the interface. The main

disadvantage of this method is the fact that it is prone to result in numerical diffusion. The diffusion will make the interface to spread over several mesh cells and the interface between the two immiscible fluids is going to be no longer sharp. In reality, the interface remains sharp due to the surface tension and the action of gravity, which separates immiscible fluids of different densities [2]. The interface of the two-phase flow must be tracked by employing some special treatments because its motion can not match the mesh any more as the calculation proceeds. A lot of techniques have been developed to cope with the problems of the multi-fluid flow systems in the past decade. These techniques can be classified to three main categories: 1. tracking the interface by using a set of mass-less particles; 2. using several mass-less maker particles to point out the only one kind of fluids and interface; 3. capturing the interface by a indicator function, such as a level set function or a volume fraction function. There are several Eulerian schemes will be introduced in the following part.

(A) Front tracking method

This front tracking method [9] (Fig. 1.2) is applied to construct the interface between the liquid and gas by a simple trajectory technique. A lot of mass-less particles are uniformly distributed over the interface in the first instance, but the Navier-Stokes equations are solved in a fixed and Eulerian grid system. The numbers of particles on the interface may be increased or reduced as the calculation proceeds. The new positions of these particles can be obtained by integrating the Eulerian fluid velocity field near the particles for each time step. This method has been used to deal with the motion of rising bubble, the breaking of water waves, and the collapse of an unsupported water column. This method is quite accurate but rather complex. Its first drawback is that the re-meshing of the Lagranian mesh is needed. Another difficult is that transforming the mesh data of the Lagranian system into the Eulerian is quite complicated. In the three-dimensional problems, the front tracking method has another problem because the particles on the interface are not a string one any more. This problem will cause the calculational time and the computer storage to increase significantly.

(B) Level set method

A level set method for moving interfaces was proposed in [10]. The interface is identified as the zero level set of a smooth distance function from the front of the interface. This method not only eliminates the problems of the numerical diffusion which will smear the sharp front, but also avoids adding or reducing points to the moving grid. This method presents the interface by solving a scalar convection equation of the level set function. This method is easy to code due to the use of Eulerian grid and can result in more accurate results when the flow motion of the interface coincides with one of the coordinate axis. This method can also be easily generalized to three dimensions. However, the main drawback of this method is that level set methods loses its accuracy because the mass is not conserved when the interface is significantly deformed. Sussman et al. used it to simulate the flows of bubbles and droplets [10], and Li presented the results of Rayleigh-Taylor instability [11].

(C) Marker and cell method

Marker and cell method (MAC) (Fig 1.3) was proposed by Harlow and Welch in [12]. Several mass-less marker particles are distributed over a space which is filled with one particular fluid with a free surface, and these marker particles are used to calculate the motion of the flow field including the free surface. This method is quite accurate and can be used accurately to deal with many complex problems, such as an interface subjected to shearing and vorticity, and wave breaking in two-dimensional system, but it may become expensive of operating in three-dimensional one. More marker particles will be added when treating problems with interface stretching, shrinking, breaking, or merging in three dimensions. The above process results in increase of computational time and computer storage. There have been many studies about this method [13, 14].

(D) Volume-of-fluid method

In volume-of-fluid method, the fluids of two-fluid flow are represented by one scalar indicator function called volume fraction. The value of the volume fraction is bounded

between zero and unity. The value of unity denotes one of the fluids. The volume of zero denotes the other fluid. The volume fraction value between zero and unity indicate the interface. This method is quite popular and easy to code in the finite-volume method. The scalar indicator function is convected through the computational domain by solving a scalar convective equation like other transport equations. The scalar indicator function can not maintain a step function on the interface because most convective schemes result in numerical diffusion and dispersion. There are three categories of this volume-of-fluid method as follows.

Line techniques

This method has been implemented in two-dimensional problems, but the reconstruction of the interface in three-dimensional flows is difficult. The methods are used for interface reconstruction can be classified into three categories as follows.

The first method is SLIC method (Simple Line Interface Calculation) which was proposed by Noh and Woodward in 1976 [15]. The interface is approximated by using lines parallel to one of the coordinate axes. The volume fractions of the left and right cells of the prime cell are used to reconstruct the interface in the prime cell approximately when the sweeping direction coincides with the x-axis. On the other hand, the volume fractions in the cells above and under the prime cell are used in the y-axial sweeping. The second method is the one with improvement on the SLIC method by Chorin [16]. All direct neighbors of the prime cell will be used for interface reconstruction in the prime cell. The third method (PLIC or Youngs’ VOF) which is posed by Youngs [17] is more accurate than the SLIC method. In this method, the interface is approximated by using oblique lines. Unlike the SLIC method, all neighboring cells are used to approximate the slope of the interface in Youngs’ VOF. The above methods are illustrated in figure 1.4.

Donor-acceptor techniques

In this method, the value of the volume fraction transported through a cell face between two cells can be approximated by the volume fraction value of the downwind cell (Figure 1.5).

This method will cause the volume fraction values unbounded, i.e. the values of the volume fraction may become greater than one or less than zero. In order to ensure the boundedness, the method which improves the level of volume fraction value on the face by using the value of the donor cell is proposed by Ramshaw and Trapp [18] , but it will cause the incorrect steeping on the interface due to change any finite gradient into a step. As mention above, the volume fraction on the face correlates closely with the flow and interface direction. Another method was proposed to cope with the above problem in [19]. Hirt and Nichols calculate the volume fraction value on the face by including some information on the slope of the interface into fluxing algorithm.

High-order differencing schemes

In this method, the convective scalar transport equation is discretised by using a high order scheme or a blending scheme to predict the interface of two-fluid system. The main errors of this method are numerical diffusion which smears the front of the fluids and numerical dispersion which causes non-physical oscillation. The first-order upwind scheme has the numerical diffusion. This diffusion becomes significantly strong when the flow direction is normal to the interface direction. In order to reduce the numerical diffusion, the linear upwind scheme (LUDS) [20] and the quadratic upstream interpolation for convective kinematics (QUICK) scheme [21] were proposed. The former is second-order accurate and interpolated by the two upwind values. The latter is third-order accurate and interpolated by the two upwind and one downwind values. These high order schemes can reduce numerical diffusion, but they may cause numerical dispersion, such as oscillations, in the strong gradient regions.

In order to cope with the dispersion problem, the flux-blending and flux-limiter technique have been proposed. The former can be classified into two classes. The first class is based on adding an anti-diffusion flux to a first order upwind scheme [22] and used to resolve sharp gradient without over-/under-shoots. The second method is based on introducing some

smoothing diffusive fluxes into an unbounded high-order scheme, and it can prevent oscillations. The flux-blending technique will become expansive due to their multi-step nature and balancing the two fluxes. In [22], the flux corrected transport (FCT) method has posed. FCT schemes are non-diffusive in nature, but create unphysical flotsam and jetsam.

The flux-limiter technique can remove non-physical oscillations and is based on the numerical flux on the interface of a cell which can be adjusted by using the flux-limiter function that enforce the boundedness. High resolution schemes (HRs) [23] are the schemes which obey the above criterion. The methods, such as Normalized Variable (NV) and Normalized Variable Diagram (NVD) [24], can be used to employ the flux-limiter technique. The flux limiter function is presented by Van Leer [25]. Sweby developed the Total Variation Diminishing (TVD) [26] approach for high resolution schemes. In the past decades, many high resolution schemes have been proposed, such as SMART of Gaskell and Lau [27], GAMMA of Jasak [28], SUPERBEE of Roe [29], STOIC of Darwish [30], MUSCL and Van Leer of Van Leer [31].

Numerical diffusion can be classified into two main components, namely cross-stream and stream-wise. These two numerical diffusions can be associated with the angle between the flow and interface direction. The blending strategy was proposed in order to improve the accuracy and less numerical diffusion including cross-stream and stream-wise. The key issue in the composite scheme is not just when to switch, but how to switch [32]. Hence, the best approach must have a continuous switching function whereby the values of compressive and high resolution schemes are blended together with a blending factor. This method has been used in utilized in the HRIC of Muzaferija [33], STACS of Darwish [34], CICSAM of Ubbink [35] and the composite of MUSCL and SUPERBEE [36].

1.3 Outline of this thesis

In this study, a high resolution scheme with switching function in the volume-of-fluid method used to solve the two-fluid flow is employed. This composite scheme can enhance the

accuracy of numerical results, preserve the sharpness on the interface and prevent the smearing.

In Chapter 2, the governing equations used to simulate the two-fluid flow, such as the continuity, momentum and volume-of-fluid equations will be introduced. The surface tension term in the momentum equation and the boundary conditions used in our simulation will be addressed.

In Chapter 3, governing equations will be discretized by using the finite volume method. The coupling between the velocity and pressure is treated by the PISO algorithm. The solution procedure of our numerical method will be introduced.

In Chapter 4, the method used to calculate the face value of the volume fraction will be introduced, and several schemes will be formulated. The former switching function will be introduced and a new composite scheme will be developed.

In Chapter 5, four cases will be tested. The first two cases, such as uniform density flow and shear flow will be used to evaluate the accuracy and availability of our numerical method by comparing with the result in the previous papers. The last two cases will be simulated by the new composite scheme. The results of them will show the superiority and accuracy of this method against the other composite schemes.

Chapter2 Mathematical Model

2.1 Introduction

In our calculation, the different fluids can be defined as a single and continuous fluid, and the fluid properties have a jump at the interface. The volume fraction is denoted as a step function on the interface. This volume fraction will be used to affect the properties of the fluid and separate the two immiscible fluids by a well defined interface. The main subject of this chapter is the description of the mathematical model which is used to solve the two-fluid system. In this study, there are several basic assumptions in our mathematical model. The model simulates the unsteady, incompressible, viscous, two-dimensional, and two-fluid systems, and the body force term includes both gravity and surface tension force. The surface tension term is a rather important effect on the interface in the two-fluid flow. The surface tension is modeled by the continuum surface force (CSF) model proposed by Brackbill et al. [37].

2.2 General transport equation

The conservation laws for mass, momentum, and energy are used to describe the physical behavior of the fluid flow. The general form of the conservation equation for a flow property φ in the control volume (C.V.) system shown in figure 2.1, is

C D S V S S V S d F dS F dS Q d Q dS t ϕ ∀ ∂ ∀ + ⋅ = ⋅ + ∀ + ⋅ ∂

∫

∫

∫

∫

∫

r r r (2.1)

where t is the time, Q the internal source, _V Q the source at the boundary, _S F_C =ρVrϕ the

flux over the boundary due to convection, Vr the fluid velocity, FD the flux over the

boundary due to diffusion, ∀ the control volume and S the control surface. We can use the Gauss’s theorem to rewrite equation (2.1):

C D S V V V V V d F d F d Q d Q d t ϕ ∀ ∂ ∀ + ∇ ⋅ ∀ = ∇ ⋅ ∀ + ∀ + ∇ ⋅ ∀ ∂

∫

∫

∫

∫

∫

(2.2)

The above equation can be rewritten to the general conservative differential form when the control volume is contracted to a single point:

S V D C F Q Q F t +∇⋅ =∇⋅ + +∇⋅ ∂ ∂ϕ (2.3)

2.3 Mass and momentum conservation equation

The conservation equations for mass and momentum can be obtained by substituting φ=ρ by neglecting the diffusion term and the source terms, and φ=Vr with the assumption of a

laminar Newtonian working fluid under unsteady and incompressible conditions with body force and surface tension force. The mass and momentum conservation equations can be written as following: 0 = ⋅ ∇ + ∂ ∂ V t r ρ ρ (2.4)

( )

ρ

( )

μ ρ σ ρ f g V p V V t V + + ∇ ⋅ ∇ + −∇ = ⋅ ∇ + ∂ ∂ r rr r (2.5)

where ρ is the fluid density, Vr the velocity, P the pressure, μ the viscosity coefficient,

g the gravitational acceleration and f the surface tension. _σ

The density and viscosity of the effective fluid in the equation (2.4) and (2.5) can be calculated by the volume fraction, as

) 1 ( 2 1α ρ α ρ ρ = + − (2.6) ) 1 ( 2 1α μ α μ μ = + − (2.7)

where the subscripts 1 and 2 denote the two fluids, α is the volume fraction. The density and viscosity of the different fluids are considered as variables through the full domain but constants in each kind of fluids. As mentioned above, all properties are piecewise continuous due to the volume fraction.

2.4 VOF equation

In the following, we will define different fluids by using the volume fraction in the volume-of-fluid (VOF) method on the Eulerian grid system (see figure 2.2). The value of volume fraction can be defined as

Volume Contral of Volume Total fluid of Volume 　 　 　 　 　　 　 　 1 = α (2.8)

categories by the volume fraction as ⎪ ⎩ ⎪ ⎨ ⎧ < < = region al transition the inside s po the for fluid the inside s po the for fluid the inside s po the for 　 　 　 　 　 　 　 　 　 　 　 　 　 ０　 　 　 　 　 　 　 　 int 1 0 2 int 1 int 1 α α (2.9)

The two-fluid system is propagated as the Lagrangian invariant and thus has a zero material derivative [19]: 0 = ∇ ⋅ + ∂ ∂ = α α α V t Dt D r (2.10) The above equations (2.4) to (2.10) can describe the fluid flow of the two-fluid system. However, the form of the volume fraction, shown as (2.10), is not a conservative one and is not suitable for numerical solution. Because of the reason, it must to be reformulated [33]. The mass conservation equation (2.4) is a conservation form. It can be rewritten as

0 1 1 (ln ) V V t D D V V t Dt Dt ρ _{ρ ρ} ρ _ρ ρ ρ ρ ρ ∂ + ⋅∇ + ∇ ⋅ = ∂ − ∂⎛ ⎞ − ⎛ ⎞ ⇒ ∇ ⋅ = _⎜ + ⋅∇ _⎟= _{⎜ ⎟}= − ∂ ⎝ ⎠ ⎝ ⎠ r r r r (2.11)

This non-conservation form of the mass conservation equation is much suitable for the two-fluid system with high density ratio, because the Vr on the interface is defined as

continuous. Figure 2.3 shows the densities at the inlet and outlet are not the same in the closed domain. The velocity Vr of the fluid of entering and leaving the domain is the same, but the

momentum Vρr of the fluid entering and leaving the domain is different. In addition, the fluids of this study are the assumption that they are incompressible. By substituting equation (2.6) into equation (2.11) the non-conservative equation becomes

(

1 2

)

2 2 1 1 = = = 0 D V Dt D Dt α ρ ρ ρ ρ ρ ρ α ρ − ∇ ⋅ ⎡_⎣ − + ⎤_⎦ − ⎛ ⎞ ⎜ ⎟ ⎝ ⎠ r (2.12)

condition by recognizing that ∇⋅αVr =α⋅∇Vr+Vr⋅∇α as: 0 = ⋅ ∇ + ∂ ∂ V t r α α (2.13) In the present study, the continuity equation (2.12), the momentum equation (2.5) and the VOF equation (2.13) together with the equation (2.6) and (2.7) will be employed to model the two-fluid flow.

2.5 Surface tension

As mentioned above, the surface tension will be modeled by the continuum surface force (CSF) model [37]. Surface tension creates a pressure jump which supplies the mean interface curvature with its necessary work on the interface. The surface tension coefficient σ exists for any pair of fluids and its magnitude is determined by the nature of the fluids. The value of σ is always positive for immiscible fluids and negative for miscible fluids [38]. The pressure jump is a function of the mean interface curvature, and it can be shown as [39]:

1 2 1 1 i o P P P R R σ⎛ ⎞ σκ Δ = − = _⎜ + _⎟= ⎝ ⎠ (2.14)

where R and ₁ R are the principal radii of curvature of the surface, ₂ P is the pressure on i

the concave side of the curved surface, P the pressure on the convex side, o σ is the surface

tension coefficient and κ is the mean interface curvature. For κ>0 fluid 1 lies on the concave side of the interface and for κ<0 fluid 2 lies on the concave side (figure 2.4). The gradient of α which is zero everywhere except at transient region, gives the normal vector, which always point from fluid 2 toward fluid 1 (figure 2.4):

α

∇ =

nr (2.15)

Thus, the mean interface curvature κ can be rewritten in terms of divergence of the unit normal vector as:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∇ ∇ ⋅ −∇ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⋅ −∇ = α α κ n n r r (2.16)

momentum equation can be expressed as α α α σ α σκ σ _⎟⎟∇ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∇ ∇ ⋅ ∇ − = ∇ = Δ = ∇ = P Pn f r (2.17)

2.6 Boundary conditions

Inlet: A velocity distribution is specified at the inlet.

Outlet: The outlet boundary condition uses the fixed pressure boundary condition. The

boundary values are obtained from convective boundary condition [40] 0 = ∇ ⋅ + ∂ ∂φ _φ c V t v (2.18)

where φ represent the transported property and Vc

v

is the convective velocity.

Rigid boundary (walls): A rigid boundary is generally defined as a non-slip boundary

Chapter3 Numerical Method

3.1 Introduction

In the chapter 2, the mathematical model of the two-fluid flow has been described in detail. It is a necessary to choose a suitable discretization method, such as the finite difference (FD), finite volume (FV) or the finite element (FE) methods. These methods approximate the differential equations by a system of algebraic equations. Finite volume method is the method which uses integral form of conservation equations. The calculated domain can be divided into many several control volumes. In our computation, the VOF equation and the momentum will be discretized by using the finite volume method. The coupling between pressure and velocity will be treated by the PISO algorithm [41].

3.2 Discretization of the VOF equation

The finite volume method for the VOF equation of equation (2.13) is first integrated over a control volume, and then can be transformed by the Gauss divergence theory as:

0 V V d Vd t α _α ∂ ∀ + ∇ ⋅ ∀ = ∂

∫

∫

r (3.1) 0 C S d V dS t α _α ∂ ∀ + ⋅ = ∂

∫

∫

r r (3.2)

The unsteady term can be discretized as:

(

n o

)

D D V d t t α _{α α} ∂ _{∀ ≈} Δ∀ ₋ ∂ Δ

∫

(3.3)

where Δ∀ is the volume of the cell, the superscripts n and o denote respectively the new and old time steps, D is the donor cell.

The second term can be discretized as:

* ( ) ( )_{f f} f S V

α

⋅dS ≈

∑

V

α

⋅S

∫

v v v v (3.4) where

(

)

* 1 2 n o f f f α = α +α (3.5)

In the above equation, the value of *

f

α is obtained by using the second-order Crank-Nicolson

scheme.

The volume fraction on the considering face can be determined by a function of neighbor cells and a flux limiter function γ

( )

r , which will be introduced in next chapter, shown as:

( )

2 A D f D r α α α =α +γ ⎛_⎜ − ⎞_⎟ ⎝ ⎠ (3.6)

Substituting the above equations into equation (3.1) yields

∑

+ = nb n nb nb n P P A S A α α _α (3.7) t A A nb nb P Δ ∀ Δ + =

∑

(3.8) ) 0 , max( 2 1 v nb F A = − & (3.9)

( )

₍

₎

1 max( ,0)( ) 2 2 o o o V nb P A D V P nb r S F F t α γ α α φ φ α ⎡ ⎤ Δ∀ = _⎢ − − − − _{⎥ Δ}+ ⎣ ⎦

∑

& & _(3.10) where FV =Vf ⋅Sf uuv uuv

& is the volume flux, the subscript P denotes the primary node and nb

denotes the neighbor node.

3.3 Discretization of the momentum equation

The momentum equations can be expressed by

( )

V

(

)

Q t φ ρφ _{ρ φ μ φ} ∂ + ∇ ⋅ = ∇ ⋅ ∇ + ∂ r (3.11)

where Q_φ is the source term of momentum equation, and φ represents the velocity components. Then, take a volume integral of the above equation and make use of Gauss theorem to yield:

( )

(

)

V S S V d V d S d S Q d t φ φ ρ∂ ∀ + ρ φ ⋅ = μ φ∇ ⋅ + ∀ ∂

∫

∫

r uv

∫

uv

∫

(3.12)

3.3.1 Unsteady term

(

)

o n o P P P V d t t ρ ρφ Δ∀ _{φ φ} ∂ ∀ ≈ − ∂ Δ

∫

(3.13)

3.3.2 Convection term

The surface integral for the convection term can be approximated by

( )

( ) C f f f f f f f f S V d S S F m ρ φ ⋅ ≈

∑

ρ φ∇ ⋅ =

∑

=

∑

φ

∫

r uv v & _(3.14)

where m&_f is the mass flux through the considering face.

The convective flux on the considering face can be estimated as a function of neighbor cells by a flux limiter function γ

( )

r , which will be introduced in the next Chapter.

( )

_⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ − + = 2 D A D f r φ φ γ φ φ (3.15)

where the subscript D denotes the donor node and A the acceptor node. In our calculation, the Van Leer scheme will be employed in the momentum equations.

The form of flux limiter function of Van Leer will be described in Chapter 4. The convetive fluxes at donor node and acceptor node can be calculated as:

⎪⎩ ⎪ ⎨ ⎧ < = = > = = 0 , 0 , f P A nb D f nb A P D m for m for & & 　 　 　 　 　 　 φ φ φ φ φ φ φ φ (3.16)

where the subsript P denotes the primary node and nb stands for the neighboring node (see Figure 3.1).

3.3.3 Diffusion term

The surface integral of the diffusion term is approximated by

(

)

( o ) D f f f f f S d S S F μ φ∇ ⋅ ≈

∑

μ φ∇ ⋅ =

∑

∫

uv v (3.17) Let,

(

S d

)

d Sr_f = r+ r_f − r (3.18)

neighboring volume center (Figure 3.1). The length dr was considered to be the factor affecting the diffusion dominancy and numerical stability. Hence, the over-relaxed approach for dr was introduced

2 f d d f S d S δ δ = ⋅ r r r r r (3.19)

The diffusion flux can then be expressed as:

(

)

(

S d

)

S S F _nbn _Pn o_f o_{f f} f d f o f D f r r r r r − ⋅ ∇ + − ⋅ = φ φ μ φ δ μ 2 (3.20)

3.3.4 Source term

The volume integral of the diffusion term is approximated by

( )P

V

Q dφ ∀ ≈ qφΔ∀

∫

(3.21)

The source terms in the momentum equation are pressure, gravitational acceleration and surface tension terms. In the following, each term will be introduced.

Pressure term

The surface integral of the pressure term is approximated as

o o f f f S Pd S =

∑

P S = ∇ Δ∀P

∫

uv r (3.22)

Gravitational acceleration term

The gravitational acceleration term can be approximated as

o V

gd g ρ ∀ = ρ Δ∀

∫

v v (3.23)

Surface tension term

The volume integral of the surface term can be obtained from the approximation

( )

_∑

∀ Δ = ∇ f f f P S α

α 1 v , where αf is the face value obtained by the interpolation from the

( )

( )

1 P P f f f V _P f f f f f _f f d S S S α σκ α σκ α σ α α α σ _α α ⎡ ⎛ _∇ ⎞⎤ ∇ ∀ = ∇ Δ∀ = − ⎢∇ ⋅⎜_{⎜ ∇} ⎟_⎟⎥ _Δ∀ Δ∀ ⎢ ⎝ ⎠⎥ ⎣ ⎦ ⎡ ∇ ⎤ − = ⋅ ⎢ ⎥ Δ∀ _⎢⎣ ∇ _⎥⎦

∑

∫

∑

∑

v v v (3.24)

3.3.5 Arrangement of the difference transport equation

The discretized form of the momentum equation is approximated by the following form.

n n o P P nb nb nb Aφ =

∑

A φ +Q_φ − ∇ Δ∀P _(3.25) where

∑

+ Δ_Δ∀ = nb P nb P t A A ρ

(

,0

)

max 2 f f d f o f nb m S S A r r & r − + ⋅ = δ μ

( )

₍

₎

₍

₎

(

)

2 o o f A D f f f f o o o o P P P r Q m S d g t φ γ φ φ μ φ ρ φ ρ σκ α ⎧ ⎫ = _⎨− − + ∇ − _⎬ ⎩ ⎭ Δ∀ + + Δ∀ + ∇ Δ∀ Δ

∑

& r r v 　　 (3.26)

3.4 Pressure-velocity coupling of the PISO algorithm

The method of Pressure-Implicit with Splitting of Operators, which is proposed by Issa [41], is called PISO algorithm. In this study, the PISO algorithm will be used to deal with the unsteady problems. In the following, the procedure of PISO algorithm is addressed.

Predictor step

The predictor step is to solve the momentum equation using the prevailing pressure field.

(

)

* * o

P P nb nb P nb

A Vr =

∑

A Vr + S− ∇ Δ∀P _(3.27)

The above equation solves the velocity field but the mass conservation law has not been satisfied yet. Dividing the above equation by A yields _P

* * nb nb nb o P P P P P P A V S V P A A A ⎛Δ∀⎞ = + −_{⎜ ⎟} ∇ ⎝ ⎠

∑

r r _(3.28) Let

P nb nb nb P A S V A H + =

∑

* * r (3.29) Then, * P Vr is written as * * o P P P P P V H P A ⎛Δ∀⎞ = −_{⎜ ⎟} ∇ ⎝ ⎠ r (3.30)

In the equation (3.29), the superscript “─” stands for the value interpolated from the primary cell P and the neighbor cell C with a weighting factor w . f

First corrector step

The corrector steps are taking care of the mass conservation law of the flow field by updating the corresponding pressure. The new velocities and the corresponding new pressure assumed to be obtained from the first corrector step are denoted with superscript ** and *.

(

)

** * *

P P nb nb P nb

A Vr =

∑

A Vr + S− ∇ Δ∀P _(3.31)

Dividing the above equation by A and using the definition of (3.29) P

* * ** nb nb nb * * P P P P P P P P P P A V S V P H P A A A A ⎛Δ∀⎞ ⎛Δ∀⎞ = + −_{⎜ ⎟} ∇ = −_{⎜ ⎟} ∇ ⎝ ⎠ ⎝ ⎠

∑

r r _(3.32)

The first velocity correction equation (V ′P

v ) is obtained by ** * V V V′≡ − , P′≡P* −Po, and shown as P P P P P A V _⎟⎟ ∇ ′ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − = ′ v (3.33)

At the cell face, the velocity correction equation can be obtained in a similar manner.

f f p f P A V _⎟⎟ ∇ ′ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ∀Δ − = ′ r (3.34)

f f f P f f f P S A S Vr r r & ⋅ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − = ⋅ ′ = ∀′ (3.35)

Over-relaxed approach is employed to let Sf = +d

(

Sf −d

)

v v

v v

for better numerical diffusion control.

(

)

(

)

⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − ⋅ ′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ + ⋅ ′ ∇ ⋅ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − = ∀′ P S d A P S S A P _f f f Pnb f f Pnb f f P f r r r r r & δ δ 2 (3.36) where _nb f nb f S S d δ δ v v v v v ⋅ = 2

Then, replacing the term Pf δPnb

r ⋅ ′ ∇ by Pnb′ −PP′ yields

(

)

(

)

⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − ⋅ ′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ + ′ − ′ ⋅ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − = ∀′ P S d A P P S S A P _f f f P nb f Pnb f f P f r r r r & δ 2 (3.37)

The first volume flux correction equation can be presented as

(

)

(

)

⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ − ⋅ ′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ + ′ − ′ ⋅ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ∀ = ∀ P S d A P P S S A _{Pnb f} nb P _{P f} f f f f P f f r r r r & & δ 2 * * * (3.38)

The volume flux at the face ( *

f

∀& ) of the above equation can be obtained by the following.

The relationship between velocity and the pressure at the face can be written as the form of (3.30) similarly. o f f P f f P A H V _⎟⎟ ∇ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − = * * r (3.39) where o f f P f f P A V H _⎟⎟ ∇ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ + = * * r (3.40) * f

Vr and ∇P_fo are written by a weighting factor w and shown as _f o P f o nb f o f w P w P P = ∇ + − ∇ ∇ (1 ) (3.41)

* * * _{(1 )} P f nb f f w V w V Vr = r + − r (3.42)

Substituting equation (3.40) into equation (3.39) to obtain:

o f f P f f P f f P A P A V V _⎟⎟ ∇ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ + = * * r r (3.43) where f P A ⎟⎟_⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ

is an average value of primary cell and neighbor cell.

⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ + ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ = ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ nb P P P f P A A A 2 1 (3.44)

Then, the volume flux at the face ( *

f

∀& ) can be obtained.

The continuity equation is discretized as

The pressure correction equation is obtained by substituting (3.38) into (3.46).

where d P P A S V S P P A S V S V o f f f P f f f o f f f P f f f f f ⋅ ∇ − ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⋅ ≈ ⋅ ∇ − ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⋅ = ⋅ = ∀ ) ( ) ( * * * * * * r r r r r r r & (3.45) 0 * * * _{= ∀ + ∀′ =} ∀

∑

∑

∑

f f f f f f & & & _(3.46) 2 1 * P P f f f nb P nb P P PP A P S S A ′ =

∑

′ −

∑

∀& + + _(3.47)

(

)

∑

∑

∑

− ⋅ ′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ = ∀ − = ⋅ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ = = f f f f P P f f P f nb f f P P nb f P nb P P d S P A S S S S A A A A v v & v v v 2 * 1 2 δ

Second corrector step

To enhance the SIMPLE procedure PISO performs a second corrector step. Similarly, the momentum equation is taken as

(

)

∑

+ −∇ Δ∀ = nb P nb nb P PV A V S P A r*** r** ** _(3.48)

where the new velocities and the corresponding new pressure are denoted with superscript *** and **.

Similarly, the second velocity corrector can be deduced as

f f P f P f nb nb P P A A V A V _⎟⎟ ∇ ′′ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ′ = ′′

∑

r r (3.49) where *** **_, ** *_, ** * P P P V V V V V V′′≡ − ′≡ − ′′≡ −

Again, the volume flow rate corrector is

f f f P f P f nb nb f P S A A V A r r & ⋅ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ′′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ − ⎟⎟ ⎟ ⎟ ⎠ ⎞ ⎜⎜ ⎜ ⎜ ⎝ ⎛ ′ = ∀ ′′

∑

(3.50)

The second corrector step is

2 2 2 1 P P f nb P nb p P PP A P S S A ′ =

∑

′ + + _(3.51) where

(

)

∑

∑

∑

− ⋅ ′′ ∇ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ = ⋅ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ ′ − = ⋅ ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ ∀Δ = = f f f f P P f f P f nb nb P f Pnb f f P P nb f P nb P P d S P A S S A V A S S S A A A A r r r r v r 2 2 2 1 2 δ (3.52)

Although more corrector steps are needed completed satisfy for the conservation law, two corrector steps are sufficient to have the accuracy of solution within temporal truncation error.

Solution procedure of PISO

Step 1. read the velocities and pressure of the flow field from the old time level. Step 2. solve the momentum equation (3.27) to get V *

Step 3. compute P′ by solving the first pressure correction equation (3.47) to update

velocities and pressure to get **

V and P . *

Step 4. compute P ′′ by solving the second pressure correction equation (3.51) to further

update velocities and pressure to get ***

V and P . **

Step 5. if the required time step is achieved, then stop the calculation and output the data

otherwise proceed to the next time step and repeat all over the way from step 1 to step 4.

3.5 Solution procedure

The VOF equation and momentum equation have been discretized. In this section, the solution procedure of two-fluid flow system will be described.

1. Initialize all variables at initial time t ₀

2. Solve the VOF equation for volume fraction α by using the old time volumetric fluxes. 3. Update the coefficients of the momentum equations. Use the new α values to obtain an

estimate for new viscosity and density.

4. Solve the momentum equation and continue with PISO algorithm.

5. If the final time step has not been reached, advance to the next time step and return to step 2.

Chapter4 High Resolution Differencing Schemes

4.1 Introduction

As mentioned above, an effective scheme adopted to solve the two-phase flow must have some features, such as small diffusion, boundedness, and maintenance of sharp interface. The schemes which are introduced in the first chapter can not include above features at the same time. For example, the first-order upwind difference scheme is bounded but too diffusive. The problem which results from numerical diffusion is very important in the two-fluid system with the interface between two fluids. The strong numerical diffusion will smear the characteristic of the step function on the interface. In this chapter, the linear and non-linear scheme will be presented, and a composite scheme with switching function which can maintain the sharpness and shape of the interface will be proposed.

4.2 Convective flux of volume fraction

As mentioned above, the discretization of the VOF equation has been established. In the calculation process, the value on the considering face of a cell must be approximated. The approximation of the face value is necessary to ensure the accuracy and stable. The face value

f

α can be estimation about a function of neighbor cells. Only two neighbor cells should be

considered in the unstructured grid system. Figure 4.1 shows a control volume and its neighbor cells including the upwind and accepter. The subscript U, D and A denote upwind cell, donor cell and accepter cell.

In general, the methods, such as the upwind difference scheme (4.1) and the central difference scheme (4.2), are adopted to approximate this value.

f D α =α (4.1) 1 ( ) 2 f A D α = α +α (4.2)

where αf is the value of volume fraction on the face, αD the value of volume fraction of

These two schemes can be associated by the variable of γ shown as: ( ) 2 f D A D γ α =α + α −α (4.3)

The term which has γ in the equation (4.3) is called the anti-diffusion correction to the upwind differencing. When γ=1, the equation (4.3) becomes central difference scheme. Then, this approximation will result in oscillations in the regions where the gradients are large. Because of the above, the variable γ must be limited. The schemes with limitation will present high accuracy and resolution results which guarantee boundedness. The schemes with total variation diminishing (TVD) flux limiters were proposed in [23] in order to ensure the bounded solution. These schemes are implemented in the context of the normalized variables formulation (NVF) [24] for the development of normalized variables diagram (NVD) schemes originally. The limiter γ is defined as a function of the gradient r [26], shown as:

D U A D r α α α α − = − (4.4)

where α_U is the value of volume fraction of the upwind cell. The value of α can be normalized as [24]:

U A U α α α α α − = − % (4.5)

with the normalization, we can get following equations:

0 1 D U D A U f U f A U U A α α α α α α α α α α α α − = − − = − = = % % % % (4.6)

Introduce this normalized variable into equation (4.3), and it can be rewritten as following equation, and α~ is just the function of f α~D.

1 ( )( ) 2 where 1 D f D D D r r α α α γ α α − = + = − % % % % % (4.7)

4.3 CBC and TVD constraints

In order to ensure a bounded value, the high resolution schemes (HRs) [23] must satisfy the Convective Boundedness Criterion (CBC) or total variation diminishing (TVD) condition. The high resolution schemes will prevent the oscillation or wiggles and get more accurate results around shocks and discontinuities in the two-fluid flow simulation. The Convective Boundedness Criterion (CBC) which was proposed by Gaskell and Lau [27] can be shown in the NVD (Fig 4.2): 0 1 1 0 1 f D D D D f D for or for α α α α α α α = < > ⎧⎪ ⎨ _{≤ ≤ ≤ ≤} ⎪⎩ % % % % % % % 　 　 　 　 　 　　 (4.8)

Sweby [26] has proposed another constraint which makes the scheme satisfy the TVD condition and it can be shown as:

( ) 0 ( f , ( )) 2 f r r r γ γ ≤ ≤ (4.9)

The above constraints can be illustrated from Fig. 4.2 to Fig. 4.4. Fig. 4.2 and Fig 4.3 present the comparison with TVD constraint and the CBC constraint in the NVD. TVD constraint in the TVD diagram is showed in Fig. 4.4, and the hatched region is known as the second-order regime.

4.4 Linear and non-linear schemes

In this section, several high order schemes will be introduced. Generally speaking, these schemes can be divided into linear and non-linear schemes. First, linear schemes are built and can be explained by using of the combination of UDS and an anti-diffusion term or by the flux limiter function γ( )r . The normalized variable and the flux limiter function of the linear

schemes can be found in table 4.1. Furthermore, these schemes will be plotted in the normalized variables diagram (NVD) and total variation diminishing (TVD) diagram (see Fig.