.
... .1 IS. ,
L, 1.. IA:
1.2.34.2

121$“.

. I ..
:Ih..?. c
. put".

I...
i J

5:»an

I“ . I.
4-,..." .. .
Unvmm
. n

V .
: .I (A)... it?!
1* ,

E it. 14

 

 

7 gnaww Qwfifié

 

THESS

2 we

This is to certify that the
dissertation entitled

LEVEL SET METHOD ON ADAPTIVE CARTESIAN GRID
FOR MULTIPHASE FLOW COMPUTATION

presented by

ZHU WANG

has been accepted towards fulfillment
of the requirements for the

Doctoral degree in Mechanical Engineering

 

 

 

 

 

Mamofessor’s Signature
/ Z / i l/ c 3"

Date

MSU is an Affirmative Action/Equal Opportunity Institution

 

LIBRARY
Michigan State
University

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

JLIIDIL2U11-ZUUZI

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 p:/ClRC/DateDue.indd-p.1

LEVEL SET METHOD ON ADAPTIVE CARTESIAN GRIDS
FOR MULTIPHASE FLOW COMPUTATION

By

Zhu Wang

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

2005

 

the

0rd

pic

de

512

 

 

ABSTRACT

LEVEL SET METHOD ON ADAPTIVE CARTESIAN GRIDS
FOR MULTIPHASE FLOW COMPUTATION

By

Zhu Wang

The Level Set Method, a powerful numerical approach for analyzing and computing
the motion of interfaces, has been studied and extended to 2D adaptive Cartesian grids in
order to achieve higher accuracy and better efficiency. First, a fast quadtree-based grid
generator is developed to produce 2D adaptive Cartesian mesh using an object-oriented
programming model. After that, both finite volume and semi-Lagrangian methods are
developed to solve the level set equation. Second-order accuracy in smooth regions, good
stability and convergence to viscosity solutions are demonstrated by the finite volume
level set method. Versatile motions of interfaces under passive transport velocity field,
first and second order geometric velocity fields are simulated with the semi-Lagrangian
method to show its ability. A particle correction algorithm is developed to further
improve the accuracy. The results are excellent and have a very good mass-conserving
property. Finally, the level set solver is coupled with a characteristics upwind finite
volume incompressible flow solver to tackle multi-phase flow problems. By using
unstructured adaptive Cartesian grids with automatic refinement near interfaces, the
solution accuracy and efficiency is dramatically improved. Several representative

multi-phase flow problems are tackled to evaluate the effectiveness of the method. The

computational results using the adaptive grids are compared with the results using

uniform grids to demonstrate the efficiency and accuracy of the method.

Copyright by
Zhu Wang
2005

 

 

 

 

 

 

 

CDC

6E.

it

an.

SU

in

 

ACKNOWLEDGMENTS

I would like to thank my advisor, Dr. 2.]. Wang for his invaluable guidance and for
his continuous encouragement and support throughout this work. Without his patience and
encouragement, it would not have been possible to finish this work. It has been a very

enlightening and fruitful experience to work with Professor Wang.

I gratefully acknowledge the financial support from the National Science

Foundation under contracts DCC-0305504 and DSC-O325760.

I would like to thank my committee members, Dr. Andre Benard, Dr. Charles Petty
and Dr. Farhad J aberi for their valuable time and for their constructive comments and
suggestions. Meanwhile I would like to thank Dr. John Strain and Dr. Yong Zhao for

many helpful discussions.

I am also grateful for the invaluable helps granted to me by the faculties and staffs
of Department of Mechanical Engineering, colleagues in CFD lab and close fiiends at

MSU during the course of my research.

At last, I should also thank my family for their support and continuous
encouragement. Most of all I would like to thank my wife, Hui. Her perpetual love and

constant support greatly encourage me throughout my study.

US

US

Ch

Ch

CI-

 

TABLE OF CONTENTS

 

 

 

 

 

 

LIST OF TABLES ....... viii
LIST OF FIGURES ix
Chapter 1. Introduction - - - 1
1.1 Level Set Method .............................................................................................. 2

1.2 Multiphase Flow Computation ......................................................................... 5

1.3 Objectives ............................................................................................ i ............ 7

1.4 Outline of the Dissertation ................................................................................ 8
Chapter 2. Level Set Method - - -- - -- 10
2.1 Level Set Equation .......................................................................................... 10

2.2 Properties of Signed Distance Function .......................................................... 13

2.3 Classic Level Set Numerical Schemes ............................................................ 15
2.3.1 First Order Upwind Schemes ................................................................. 15

2.3.2 Hyperbolic Schemes .............................................................................. 16

2.3.3 ENO and WENO Schemes .................................................................... 19
Chapter 3. Adaptive Cartesian Grid - - - - 21
3.1 Introduction ..................................................................................................... 21

3.2 Adaptive Cartesian Grid Generation ............................................................... 22
3.2.1 Interface Presentation ............................................................................. 22

3.2.2 Quadtree Adaptive Cartesian Mesh ....................................................... 26

3.2.3 Splitting Rules ...................................................................................... 27

3.3 Numerical Examples ....................................................................................... 28
Chapter 4. Finite Volume Level Set Method -- - -- - 32
4.1 Introduction ..................................................................................................... 32

4.2 Numerical Algorithm ...................................................................................... 33
4.2.1 Hamilton-Jacobi Equation ...................................................................... 33

4.2.2 Finite Volume Method ........................................................................... 35

4.2.3 Flux Computation .................................................................................. 37

4.2.4 Linear and Quagdratic Reconstruction .................................................. 38

4.2.5 Time Integration ..................................................................................... 43

4.3 Numerical Results ........................................................................................... 45
4.3.1 Accuracy Study with the 2D Linear Wave Equation ............................. 45

4.3.2 2D Non-Linear Burgers Equation .......................................................... 48

4.3.3 Expansion of Initially Circular Interface ............................................... 50

4.3.4 Moving Circular Interface under Passive Velocity Field ....................... 53

vi

Cl

 

 

 

 

 

 

4.3.5 Propagating Surface ............................................................................... 57

 

 

 

 

4.3.6 Star-Shaped Interface ............................................................................. 59

4.4 Summary ......................................................................................................... 62
Chapter 5. Semi-Lagrangian Level Set Method 63
5.1 Introduction ..................................................................................................... 63

5.2 Numerical Algorithm ...................................................................................... 64
5.2.1 Numerical Schemes of Semi-Lagrangian Level Set Method ................. 64

5.2.2 Fast Redistancing Algorithm ................................................................. 67

5.2.3 Whitney Velocity Extension .................................................................. 71

5.2.4 Contouring Algorithm ............................................................................ 73

5.2.5 Particle Correction ................................................................................. 75

5.3 Numerical Results ........................................................................................... 77
5.3.1 Zalesak’s Disk Rotation Problem .......................................................... 77

5.3.2 Single Vortex ......................................................................................... 84

5.3.3 Motion under First-Order Geometric Velocity ...................................... 90

5.3.4 Motion under Second-Order Geometric Velocity .................................. 93

5.4 Summary ......................................................................................................... 96
Chapter 6. Multiphase Flow Computation 97
6.1 Introduction ..................................................................................................... 97

6.2 Mathematical Formulation .............................................................................. 99

6.3 Numerical Algorithm .................................................................................... 103
6.3.1 Characteristic Upwind Finite Volume Method .................................... 104

6.3.2 Finite Volume Method ......................................................................... 105

6.3.3 Least Square Linear Reconsruction ..................................................... 107

6.3.4 Upwind Characteristic Method ............................................................ 109

6.3.5 Time Integration ................................................................................... 114

6.4 Numerical results .......................................................................................... 115
6.4.1 Broken Dam ......................................................................................... 1 15

6.4.2 Relaxation Bubble ................................................................................ 121

6.4.3 Rising Air Bubble in Fully Filled Container with Shedding ............... 126

6.5 Summary ....................................................................................................... 135
Chapter 7. Conclusions and Future Work 137
7.1 Conclusions ................................................................................................... 137

7.2 Future Work .................................................................................................. 138
APPENDICES 139
BIBLOGRAPHY 142

 

vii

 

 

 

 

la

la

la

la

la

Ta

 

 

LIST OF TABLES

Table 4-1. Accuracy Study for the 2D Linear Wave Equation, T =1.0 ............................. 48
Table 4-2. Accuracy and CPU Time for Circular Interface Expansion, T=4.0 ................. 51
Table 5-1. Zalesak’s Disk: Results with Linear and Adaptive Contouring ....................... 83
Table 5-2. Zalesak’s Disk : Semi—Lagrangian + Particle Correction ................................. 83
Table 6-1. Relaxation of Elliptical Bubble in Viscous Flow: Volume, T=O.3 ................ 122
Table 6-2. List of Test Cases for Bubble Rising Problem ............................................... 128

viii

LIST OF FIGURES

Figure 2-1. Schematic of Level Set Method ..................................................................... 11
Figure 2-2. Signed Distance Function on a 64 by 64 Uniform Cartesian Mesh ............... 14
Figure 2-3. Signed Distance Function on a Level 3 to 7 Adaptive Cartesian Mesh for
Mutli-Component Interfaces .................................................................................. 14
Figure 2-4. Expansion of a Seven-Point Star-Shaped Interface, Second- Order Scheme,
300 by 300 Uniform Cartesian Method, T = O, 0.03 (0.01) ................................... 18
Figure 2-5. The Evolution History of Star-shaped Interface Under Curvature Flow300 by
300 Uniform Cartesian Mesh, Second-Order Scheme ........................................... 19
Figure 3-1. Data Structure for 2D Interface ...................................................................... 24
Figure 3-2. Distance from any Point to a Linear Element ................................................ 24
Figure 3-3. Level 3 Quadtree Mesh and Tree Structure ................................................... 27
Figure 3-4. Seven-Point Star-shaped Interface with 1000 Elements ................................ 29

Figure 3-5. Level 3 to 10 Adaptive Cartesian Grids after Interface-Based Grid
Refinement, Number of Cells = 56764, Time of Generation =0.76 Second .......... 30

Figure 3-6. Adaptive Cartesian Grids near the Interface .................................................. 30
Figure 3-7. Level 3 to 7 Adaptive Cartesian Grids for a Multi-Component Interfaces,

Number of Cells =3751, Time of Generation =0.22 Second .................................. 31
Figure 4-1. Schematic of a Polygon Control Volume (CV) ............................................. 36

Figure 4—2. Quadratic Reconstruction Stencil of Cell i for Adaptive Cartesian Grids ..... 41

Figure 43. Initial ¢ = sin(7r(x + y)) at T =0.0 ................................................................ 47

Figure 4-4. Numerical Solution of ¢ , T =1.0, Level 6 Uniform (64 by 64) Cartesian
Mesh, Second-Order Scheme ................................................................................. 47

Figure 4-5. Numerical Solution of 2D Non-linear Burgers Equation, Level 6 Uniform (64
by 64) Cartesian Mesh, Second-Order Scheme ...................................................... 49

ix

Figure 4-6. Expansion of Circular Interface, V" =1, T= 0, 4 (0.5), Level 6 (64 by 64)

Uniform Cartesian Mesh ........................................................................................ 51
Figure 4-7. Expansion of Circular Interface: Initial Interface, T = 0, 1228 Cells ............ 52
Figure 4-8. Expansion of Circular Interface: Interface at T = 2, 2392 Cells .................... 52
Figure 4-9. Expansion of Circular Interface: Interface at T = 4, 3433 Cells .................... 53
Figure 4-10. Convection of Circular Interface: T=0, 4(1.0), Level 6 (64 by 64) Uniform

Cartesian Mesh, CPU Time/per step = 2.35e-2 second .......................................... 54
Figure 4-11. Convection of Circular Interface: T=0, 4(0.8$), Level 6 to 9 Adaptive

Cartesian Mesh CPU Time/Step is 1.38e-1 second ................................................ 54
Figure 4-12. 6 to 9 Adaptive Cartesian Mesh at T = 4, 23338 cells ................................. 55

Figure 4-13. Rotation of Circle: First-Order T = 0, 7r (7: / 6) , Level 7 Uniform Cartesian
Mesh ....................................................................................................................... 56

Figure 4-14. Rotation of Circle: Second-Order, T = 0, 7r (7: / 6) , Level 7 Uniform
Cartesian Mesh ....................................................................................................... 56

Figure 4-15. Rotation of Circle: Second Order, T = 0, 7: (7r/ 10) , Level 5 to 9 Adaptive
Cartesian Mesh ....................................................................................................... 57

Figure 4-16. Propagating Surface, 6 =0, T =0, 0.6 (0.3) (Note: 45 — 0.35 is displayed for
visual effects at T= 0) ............................................................................................. 58

Figure 4-17. Propagating Surface, 8 =0.1, T =0, 0.6 (0.3) (Note: ¢ — 0.35 is displayed
for visual effects at T= 0) ....................................................................................... 59

Figure 4-18. Expansion of Star-Shaped Interface, V,l =1, T =0, 0.07 (0.01),Level 4 to 8
Adaptive Cartesian Mesh ....................................................................................... 61

Figure 4-19. Motion of Star-Shaped interface, VN = -K , T =0, 0.005 (0.001), Level 8

Uniform Cartesian Mesh ........................................................................................ 61
Figure 5-1. Bilinear Interpolation of Off-grid Value ........................................................ 65
Figure 5-2. Fast Semi-Lagrangian Level Set Method ....................................................... 66

Figure 5-3. Fixing the Sign of Level Set Function: (a) Convex Interface; (b) Concave

Interface .................................................................................................................. 68
Figure 5-4. Concentric Triple of a Quadtree Leaf ............................................................ 69
Figure 5-5. An Example of Level 7 Distance Tree for Multiple Circlers ......................... 70
Figure 5-6. Fast Searching Strategy .................................................................................. 70
Figure 5-7. Level 6 Contour Tree and the Triangulation .................................................. 74
Figure 5-8. Initial Zalesak’s Disk (Notched Interface) ..................................................... 78

Figure 5-9. Zalesak’s Disk Problem: T=628, Level 6, 7, 8 and 9 Adaptive Cartesian Grid,
Linear Contouring .................................................................................................. 79

Figure 5-10. Zalesak’s Disk Problem: Level 7 Adaptive Cartesian Grid, Level 3 Adaptive

Contouring, T=0, 628 (157) ................................................................................... 80
Figure 5-11. Zalesak’s Disk Problem: Diffusive Property of Level Set Method,Results at
Under-Resolved Area (near the sharp comers) ...................................................... 80
Figure 5-12. Zalesak’s Disk Problem: Level 7 Adaptive Cartesian Grid, Particle Level Set
Method, T=1256 ..................................................................................................... 81
Figure 5-13. Zalesak’s Disk Problem: Level 8 Adaptive Cartesian Grid, Particle Level Set
Method, T=628, 1256 (157) ................................................................................... 81
Figure 5-14. Zalesak’s Disk Problem: Level 8 Adaptive Cartesian Grid, Particle Level Set
Method, T=1256 ..................................................................................................... 82
Figure 5-15. Initial Interface and Single Vortex Velocity Field ....................................... 85

Figure 5-16. Single Vortex Problem: T = 0,1 (0.25), Level 6 Adaptive Cartesian Grid,
Level 3 Adaptive Contouring, 2nd order Semi-Lagrangian Scheme ..................... 86

Figure 5-17. Single Vortex Problem: T=3, Level 8 Adaptive Cartesian Grid, dt=
628/2400, level 3 Adaptive Contouring, Second Order SL Scheme ...................... 87

Figure 5-18. Single Vortex Problem: T = 3, Level 7 Adaptive Cartesian Grid, dt= 0.001,
Particle Level Set Method ...................................................................................... 88

Figure 5-19. Single Vortex Problem: T = 3, Level 8 Adaptive Cartesian Grid, dt= 0.001 ,
Particle Level Set Method ...................................................................................... 88

xi

 

 

 

ii.

Pi,

 

 

 

Figure 5-20. Level 8 Adaptive Cartesian Grid, T = 3 ....................................................... 89

Figure 5-21. Single Vortex Problem: T = 5, Level 9 Adaptive Cartesian Grid, dt = 0.001,

Particle Level Set Method ...................................................................................... 89
Figure 5-22. Level 9 Adaptive Cartesian Grid, T = 5 ....................................................... 90
Figure 5-23. Expansion of Circular Interface: T=0, 2 (0.2), dt =0.02, Second-Order SL,

Level 7 Adaptive Cartesian Grid ............................................................................ 91
Figure 5-24. Expansion of A Triangle: T=0, 2.2 (0.1), dt = 0.0125, Second- Order SL,

Level 7 Adaptive Cartesian Grid ............................................................................ 92
Figure 5-25. Shrinking of Triangles: T=0, 0.8 (0.08), dt = 2.667e-3, Second-Order SL,

Level 8 Adaptive Cartesian Grid ............................................................................ 92
Figure 5-26. Expansion of Star-shaped Interface: T=0, 0.07 (0.007), dt=0.0001,

Second-Order SL, Level 8 Adaptive Cartesian Grid .............................................. 93
Figure 5-27. Shrinking of Multiple Circles under Curvature Flow, T=0, 2 Second— Order

SL. Level 8 Adaptive Cartesian Grid ..................................................................... 94
Figure 5-28. Motion of Star-shaped Interface under Curvature Flow, T=0, 0.002 (0.0005),

Second-Order SL, Level 7 Adaptive Cartesian Grid .............................................. 95
Figure 6-1. Linear Reconstruction Stencil for Cell i on Adaptive Cartesian Grid ......... 109
Figure 6-2. The Schematic of t— 5 Characteristics ..................................................... 110

Figure 6-3. (a) Illustration of 2D Broken Dam Problem at T=0; (b) Initial Level 5 to 7
Adaptive Cartesian Grids and Zero-Level Set ...................................................... 116

Figure 6-4. Free Surface Profile at T=0.24 (Red solid line is on level 6 uniform mesh,
blue dashed line is on level 5 to 7 adaptive mesh, and green dotted line is on level
5 to 8 adaptive mesh) ............................................................................................ 117

Figure 6-5. Triangulated Level 5 to 7 Adaptive Cartesian Grids, T=0.24 ...................... 118

Figure 6-6. Free Surface Profile at T=0.42 (Red solid line is on level 6 uniform mesh,
blue dashed line is on level 5 to 7 adaptive mesh, and green dotted line is on level

5 to 8 adaptive mesh) ............................................................................................ 118

Figure 6-7. Triangulated Level 5 to 7 Adaptive Cartesian Grids, T=0.42 ...................... 119

xii

 

 

 

Fl

 

 

Figure 6-8. Comparisons of Numerical and Experimental Results: (a) Non- Dimensional
Length of Surge Front against Time; (h) Non-Dimensional Column Height against
Time ...................................................................................................................... 120

Figure 6-9. Relaxation of Elliptical Bubble in Viscous Flow: T =0, 0.3 (At = 0.01), Level
5 32 by 32 Uniform Cartesian Grids .................................................................... 122

Figure 6-10. Relaxation of Elliptical Bubble in Viscous Flow: Velocity Field at T =O.3,
Level 5 Uniform Cartesian Grids ......................................................................... 123

Figure 6-11. Relaxation of Elliptical Bubble in Viscous Flow: Interface at T =0.3 on
Level 4, 5, and 6 Uniform Cartesian Grids and level 5 to 7 Adaptive Cartesian
Grids ..................................................................................................................... 123

Figure 6-12. Relaxation of Stare-Shaped Bubble in Viscous Flow: Initial Interface, Level
5 to 8 Adaptive Cartesian Grids ........................................................................... 124

Figure 6-13. Relaxation of Stare-shaped Bubble in Viscous Flow: Evolution History, T
=0, 0.5(0.05), Level 5 to 8 Adaptive Cartesian Grids .......................................... 125

Figure 6-14. Relaxation of Stare-shaped Bubble in Viscous Flow: Velocity Field, T =0.l,
Level 5 to 8 Adaptive Cartesian Grids ................................................................. 125

Figure 6-15. Relaxation of Elliptical Bubble in Viscous Flow: Volume Preserving History
versus Time Step .................................................................................................. 126

Figure 6-16. (a) Initial Setup of 2D Rising Bubble Problem; (b) Initial Level 5 to 7
Adaptive Cartesian Grids and Zero-Level Set ...................................................... 127

Figure 6-17. The Evolution of Rising Bubble on Level 6 Uniform Cartesian Grids, T=0.2,
1.0 (At =0.l). (The figures are ordered fiom top to bottom and left to right). ..... 129

Figure 6-18. The Evolution of Rising Bubble on Level 7 Uniform Cartesian Grids, T=0.2,
1.0 (At =0.1) ......................................................................................................... 130

Figure 6-19. The Evolution of Rising Bubble on Level 5 to 7 Adaptive Cartesian Grids,
T=0.2, 1.0 (At =0.l) .............................................................................................. 131

Figure 6-20. The Evolution of Rising bubble on Level 5 to 8 Adaptive Cartesian Grids,
T=0.2, 1.0 (At =0.1) .............................................................................................. 132

Figure 6-21. The Process of Pinching of Shedding Bubble on Level 5 to 8 Adaptive
Cartesian Grids, T=0.74, 0.76 and 0.78 ................................................................ 132

xiii

 

 

ii.

 

Figure 6-22. Level 5 to 8 Adaptive Cartesian Grids, t = 0.66: (a) Adaptive Cartesian
Grids at Time t = 0.66; (b) Triangulated Adaptive Cartesian Grid; (c) Velocity
Field at This Time ................................................................................................ 134

Figure 6-23. The Evolution of Rising Bubble on Level 6 to 9 Adaptive Cartesian Grids,
T=0.2, 1.0 (At =0.1) .............................................................................................. 135

xiv

Chapter 1. Introduction

Propagating interfaces occur and can be seen everywhere, such as ocean waves,
burning flames and material boundaries between different fluids etc. Such an interface is
always of great interest because the interplay between the interface and the surrounding
fluid plays a critical role in physics. It has a long history to numerically compute these
interface problems, which involve external physics and arise in various areas of science
including fluid dynamics, combustion, material science, image processing etc. There are
three numerical approaches to track the motion of the interface: front-tracking method
which is based on Lagrangian view, volume of fluid (V OF) method which discretizes the
underlying domain and fill cell fractions with values that represent the characteristic
function in those cells (usually a mass fraction) and level set method (LSM). In this

dissertation, the third approach LSM is concerned.

The level set method, a much more powerful numerical technique for tracking the
motion of the interface, has been developed greatly during past two decades. It has been
applied to some of the most complex problems in fluid-interface motion. Level set
method relies on an implicit representation of the interface whose equation of motion is
numerically approximated. The resulting techniques are capable of solving problems in
which the speed of the evolving interface (speed function) depends on local properties

such as curvature and local unit normal.

 

 

 

[A

F.)

fl

 

In this chapter, level set method is reviewed firstly. Then its application on two-phase
interfacial flow is described. After that, the objectives are presented. The last section gives

the outline of this dissertation.

1.1 Level Set Method

The level set method, introduced by Osher and Sethian [46] in 1988, is a powerfirl
numerical technique for analyzing and computing motion of interfaces in two or three
dimensional space (i.e., curves in 2D and 3D or surfaces in 3D). This method relies on an
implicit representation of the interface whose equation of motion is numerically
approximated. During the past two decades, numerical algorithms for solving level set
equation have been tremendously advanced as summarized by Osher and Fedkiw [48, 49]
and Sethian [58, 59]. At the same time, it has been applied greatly to various areas of
science, such as physics, chemistry, fluid mechanics, combustion, image processing,
material sciences, control science, computational geometry, computer-aided-design etc.
The examples of interface problems involving external physics include compressible flow
[26, 44], incompressible flow [67, 68], Stefan problems [11], kinetic crystal growth,
epitaxial growth of thin films, vortex dominated flows, body fitted grid generation, image

processing etc.

The level set methods can be mainly categorized into three groups according to the

numerical algorithm. The first group of level set solvers was focused on solving

 

 

Hamé
COIlSi
by C I
[46].
L310:
[33].
mm
high-
Zhar.
AIIIlt
cases
is ur.
inter
alter
The
geon
IllIIl.
athigi
adap

Casi?

 

Hamilton-Jacobi equation and they were extended from the schemes built for hyperbolic
conservation laws. They include the monotone scheme with finite difference (FD) method
by Crandall and Lions [19, 20], non-oscillatory hyperbolic schemes by Osher and Sethian
[46], high-order essentially non-oscillatory (ENO) schemes by Osher and Shu [47],
Lafon and Osher [39], discontinuous Galerkin (DG) finite element method by Hu and Shu
[33], Petrov-Galerkin method by Barth and Sethian [4], high order weighted essentially
non-oscillatory schemes (WENO) on uniform structured mesh by Jiang and Peng [35],
high-order weighted essentially non-oscillatory (WENO) schemes for triangular mesh by
Zhang and Shu [79], central WENO schemes by Levy, Nayak, Shu and Zhang [40].
Although these algorithms can accurately simulate the motion of the interfaces in some
cases, most of them were carried out on structured Cartesian meshes only. Since the mesh
is uniform everywhere, it is not possible to have higher resolution for the grids near the
interface. This may be very costly for some applications. There also have been a few
attempts to use unstructured triangular meshes for the level set equation such as [4, 79].
The unstructured grids are much more flexible than uniform mesh in handling complex
geometries and performing grid adaptation locally. However, the numerical results may be
influenced by the quality of the computational mesh. Adaptive Cartesian grids are
advocated in this study. With the help of tree-based data structure, the generation of the

adaptive Cartesian grids can be very fast and interface-based adaptation can be performed

easily.

 

 

 

uasr
on ur
volur
recor

adap'

set f
numb
dfilu

loss.

senn

and

 

accu
Ioth

IEVc}

lfih.

 

The second group of level set method is finite volume (FV) level set method which
was newly developed. A flux-based finite volume method was proposed by Frolkovic [27]
on uniform mesh in 2003. And in 2004 Zhu Wang and Z]. Wang [75] developed a finite
volume level set method on 2D adaptive Cartesian mesh realized by a least-square data
reconstruction algorithm. Through the automatic interface-based grid refrnement and
adaptation, the efficiency and accuracy are fitrther improved near the interface. The two
groups of level set methods described above are Eulerian methods, which solve the level
set function on some fixed background grids. Although Eulerian level set methods
mathematically defines the merging and breaking of interface well, they are numerically
diffusive. In another word, Eulerian level set method can not conserve mass well due to the

loss of information in under-resolved area.

The third group of level set method is Lagrangian level set method, which is from a
Lagrangian view. In 1999, a semi-Lagrangian scheme was developed by Strain in [63] to
simulate the interface on 2D Quadtree-mesh, and it was modulated into a fast
semi-Lagrangian level set method by Strain [65] in 2000. Through a fast redistancing [64]
and adaptive contouring algorithm [66], the motion of interface can be simulated
accurately and efficiently. One drawback of this method is that it is too intricate to extend
to three dimensional problems. In order to conserve the mass or volume better, a hybrid
level set method was developed by Enright et al. [24] in 2002. This method coupled the

fifth-order WENO schemes for level set equation with third-order Runge-Kutta scheme for

 

 

IIIO‘

v01

der

V8]

C0

sir

de

Ir‘

thi

be
lI'lt
Ilo
of.
10;)

film

moving the mass less particles. This method dramatically improved the property of mass or
volume preserving through a particle correction algorithm. This method was only
developed for passive transported velocity field. Meanwhile the quantity of particles can be
very large by distributing the particles randomly in a band of area around the interface.
Thus the computation for moving particles could be time-consuming. In 2005 Enright et al.
[25] employed first order semi-Lagrangian level set method coupled with particle
correction algorithm. Similar results as in [24] were achieved while it is more efficient
since only first order scheme was used to solve the level set equation. Even newer level set
method such as hybrid particle level set on Octree adaptive Cartesian grids is under

developing as stated in [43].

1.2 Multiphase Flow Computation

The efforts to compute the motion of multiphase flow have been a long history and
they could be as old as computational fluid dynamics (CFD). In this research, two-phase
interfacial incompressible flow is under consideration. For such problems, the interplay
between the interface and the surrounding fluids are subtle and can play a critical role on
the flow. All of these three interface tracking approaches have been used to multiphase
flow computation. Compared to front-tracking method [3, 52, 60-62, 70, 71] and volume
of fluid (V OF) [32, 8, 82], level set method is more powerful and can naturally handle the
topological changes such as breaking up and merging. Meanwhile, the terms in the speed

function of the interface involving geometric properties such as local unit normal and

 

 

 

CUIV;

appli

incur
Sim;
1995
“16>
air I
unil‘
ever
awar

be er

 

adl‘c

refin

0i S(

then

curvature can be easily approximated with the derivatives of level set function. The

application of level set method to multiphase flow is reviewed next.

Since Sussman et a1. [67] applied the level set approach for solving 2D
incompressible two-phase flow in 1994, there has been a great deal of interest in the
simulation of multiphase flows involving interface capturing with the level set method. In
1998, Sussman et al. [68] modified the numerical scheme and improved the accuracy.
These methods were accurate and the topological changes such as merging and pinching of
air bubbles were captured accurately. But these methods were designed for structured
uniform mesh only. In order to improve the accuracy, one has to refine the mesh
everywhere. This could be too costly and waste a lot of time on computing the area far
away from the interface, where velocity fields are smooth and lower grid resolution would

be enough.

In order to improve the accuracy and efficiency, the spatially adaptive techniques are
advocated to reduce memory and computation time (CPU time) through adaptively
refining the grids in under-resolved areas, where the interface usually lies, steep gradients
of solution variables exist and sharp comer occurs. As summarized in Losasso et al. [43],
there are three approaches to speed up the level set method. The first approach, Narrow
Band method, simply allocates the full amount of memory everywhere but only carries out
the computation near the interface as described in Adalsteinsson and Sethian [2] and Peng

et al. [50]. The second approach is called adaptive mesh refinement (AMR) method, in

 

 

 

CV61?

ICCIll

[Dori

 

 

mell

adl'c
dlna

adap

which a number of uniform grids of different resolution are used for different regions. It
could be wasteful in the sense that too many fine mesh cells are introduced in each region
even when only a few cells are sufficient there. Sussman et al. [69] employed such an AMR
technique on solving incompressible two-phase flows. The third approach is to using much
more eflicient quadtree-based adaptive Cartesian grids in 2D [75] and octree-based
adaptive Cartesian grids [28, 42] in 3D. The octree/quadtree based level set method is
much more general and efficient than AMR techniques and it is based on more complicated
data structures through which many operations related to searching and sorting can be done

efficiently.

In this research, adaptive Cartesian grids are used to improve the accuracy and
efficiency. A fast semi-Lagrangian level set method is coupled with a second-order
characteristic upwind finite volume method [74, 81, 82] for 2D incompressible flow to

solve multiphase problems.

1 .3 Obiectives

The first objective of this study is applying spatially adaptive techniques to level set
method to improve accuracy and efficiency. In this research, the adaptive Cartesian mesh is
advocated. The adaptive Cartesian grids have been widely used in computational fluid
dynamics (CFD). Its ability in handling complex geometries and in solution based grid

adaptations is demonstrated in many applications [72. 16]. One unique feature of the

 

 

adap
local

for a

field
accu
beer.

mei‘:

inco
this
inco

50h

bad.

 

adaptive Cartesian grid is that diverse length scales can be handled very efficiently with
local feature-based grid adaptations. For this purpose, a finite volume level set method [75]
for adaptive Cartesian mesh was developed in 2004. Following the work by Strain [65], a
second-order semi-Lagrangian level set solver is developed on 2D adaptive Cartesian mesh.
This solver is more general and it can maintain the signed distance function accurately
through a fast redistancing algorithm. The semi-Lagrangian method is capable of solving
2D interface moving problems under different velocity fields, including passive transport
fields, first-order and second-order geometric velocity. In order to further improve the
accuracy and properties of conserving volume or mass, a particle correction algorithm has
been developed for the semi-Lagrangian level set method. Thecomparison between these

methods is presented in Chapter 5.

Another objective of this research is to couple the level set method with 2D
incompressible flow solver to multiphase flow computation on adaptive Cartesian grids. In
this dissertation, the semi—Lagrangian level set method has been coupled with a 2D
incompressible flow solver to solve multiphase problems. Different validation cases are

solved and the numerical results are presented and discussed in Chapter 6.

1 .4 Outline of the Dissertation

The remainder of this thesis is organized as follows. Following the introduction, the

backgrounds and basic knowledge of level set method, including level set equation,

 

 

propi
Char
that l
5, th
The
com]

extei

 

 

properties of signed distance function and classic numerical schemes, are presented in
Chapter 2. Next the generation of adaptive Cartesian grid is described in Chapter 3. After
that the finite volume level set method is described and discussed in Chapter 4. In Chapter
5, the semi-Lagrangian level set method is presented with different numerical examples.
The level set method is coupled with a 2D incompressible flow solver to multiphase flow
computation in Chapter 6. In the last chapter, conclusions are drawn and possible future

extension is discussed based on the current study.

 

 

 

 

the I
num:

Can

Velo
by n

the(

Gen.
Into
30 i.
“fie:

Cat.

 

Chapter 2. Level Set Method

This chapter reviews basic background material of the level set approach, including
the basics and concepts of level set method, signed distance function and classic

numerical schemes of Hamilton-Jacobi solvers, which are developed on uniform

Cartesian grids.
2.1 Level Set Equation

 

Let F denote the interface bounding a region 0 (possibly multiply connected) and SE
is the position vector in Rn space. For example, the interface can be represented as curve

in two spatial dimensions and surface in three spatial dimensions. Assume that the
velocity at each point 55 on interface F is 17(3) , the motion of interface can be simulated

by moving all the points on the interface with this velocity. The simplest way is to solve

the ordinary differential equation (ODE)
(if _
— = V ‘ . 2.1
dt (x) ( )

Equation (2.1) is the Lagrangian formulation of the interface evolution equation.
Generally there are infinites points on the interface. It is usually to discretize the interface
into a finite number of pieces or elements (linear segments in 2D and small triangles in
3D). Then the endpoints are moved by solving equation (2.1). This type of methods is
referred as front-tracking method. But when topological changes exist or velocity field

causes large distortion even connectivity change of the elements, the accuracy of

10

 

 

lllllll;

spec:
detat

AII e

fune

 

numerical method can deteriorate very quickly and the method becomes unstable. Hence
special treatments such as smoothing and regularizing techniques, the surgical methods of
detaching and reattaching interfacial elements, and even re-discretization are required.

All of these can make the method rather complicated and hard to implement.

In 1988, Osher and Sethian [46] used a smooth and Lipschitz continuous signed

distance function ¢(3c',t) (at least), which represents the interface as the zero set where ¢
equals to zero. With this definition, the level set function 95 has the following properties:
¢(‘,t)< 0 for if 6 0(0’ ),
¢(",t)> 0 for 2 is mot),
¢(5c‘,t)= 0 for 56 e 60 = F(t).
Thus, the interface is to be captured implicitly by finding zero set of level set

function, i.e., ¢(Jr‘,t) is equal to zero. This is illustrated in Figure 2-1.

 

Figure 2-1. Schematic of Level Set Method

11

 

 

C011\

20%

the

prop

bash

l'elo

files

a:
I!

 

With all of these definitions, the motion of interface F (t) can be analyzed by
convecting the level set firnction ¢()‘c,t) with external velocity field 17(5c‘,t). The
governing equation for ¢ is level set equation (2.2)

:32 ‘. =
at+V W o, (2.2)

where V¢ is the gradient of ¢. 17 is an external velocity field, which is extended from

the interface to the whole computational domain. One can compute I7 through local
properties (curvature and local unit normal of interface) or external physics such as the

basic laws of fluid mechanics. The signed distance function is selected as ¢ due to its

nice properties, which are discussed in section 2.2. Because topological changes such as
breaking up and merging of interface are defined naturally by finding the zero set of the
level set function, the level set equation is significant and used to simulate the movement

of interface in 2D and 3D.

Suppose ii is out-warding normal of the interface, then the normal component of the

velocity is V" = l7 - ii . The level set equation can be rewritten in the following form,
52—? + V,,(5i,v¢,x]v¢| = 0. (2.3)

Here V" could be a function of if , V¢ and curvature K. If signed distance function

are selected as level set function, unit normal and curvature can be simply calculated by

fi=V¢ and K=V~fi.

12

 

 

 

 

56'

de

sh.

Illl

2.2 Properties of Signed Distance Function

As mentioned in previous section 2.1, the signed distance function is chosen as level
set function ¢(5c‘,t) in the numerical computation. A signed distance function d(5c’,t) is

defined as:

d(f,t)=sign(Jr‘)-min(]5c‘—-5c‘1|) for all i, 6F(t), (2.4)

where sign(3c‘) is positive unit on the exterior region (2+, negative unit on the interior

region Q" and zero on the interface F. The minimum distance between a point 55 and an
interface is the Euclidean distance between current point and the closest point on the

interface. Furthermore, the signed distance function satisfies
|w| =1. (2.5)

This property of distance function will greatly simplify the computation of
outwarding normal and curvature, which will be discussed in later chapters. In the
interface tracking and multi-phase computation, it is very important to maintain ¢ as a
signed distance function. Now it is better to make clear that he distance ftmction and level

set function will have same notation ¢(f, t) in later chapters.

Several examples of signed distance function are given to illustrate. Figure 2-2
shows a signed distance function for a single seven-point star shaped interface on

uniform Cartesian mesh. Figure 2-3 shows the signed distance function for an interface

13

 

 

 

 

ch

with multiple components on adaptive Cartesian mesh, which will be discussed in next

chapter.

4'"
¢

 

80

Figure 2-3. Signed Distance Function on a Level 3 to 7 Adaptive Cartesian
Mesh for Mutli-Component interface

 

 

 

 

IIQ

2.3 Classic Level Set Schemes

Classic numerical schemes for level set equations are developed on uniform
Cartesian grids. They include first order upwind differencing scheme, hyperbolic
schemes for Hamilton-Jacobi equation, high order essentially non-oscillatory (ENO)

schemes and weight ENO schemes.
2.3.1 First Order Upwind Scheme

If an interface is moved by externally generated velocity field, upwind differential
scheme is preferable. In such a problem, ¢ and velocity field I7 are well defined on all

grid points (or points close to the interface) in the computational domain. In order to

evolve (1) forward in time, upwind differencing schemes are applied to move the interface

across the grid. Generally, the upwind schemes choose the forward differencing or
backward differencing according to sign of velocity components u, v and w, which

decide the direction where the characteristic information is coming from.

Starting from the level set equation (2.2) directly, the first-order accurate forward

Euler differential scheme for the time discretization is
¢n+l __ ¢n _.
———At + V" -v¢" = o. (2.6)

Equation (2.6) can be firrther rewritten as

¢n+l _ ¢n

At + un¢£ + v"¢; + w"¢;’ = 0. (2.7)

15

The first order upwind scheme for equation (2.7 ) is

W“ =¢" —At(+ ax(v",0)¢; + in(v",0 ¢; ). (2.8)

 

 

2.3.2 Hyperbolic Schemes

Recall that the level set equation is the initial value formation of ¢. Ignoring the

second order curvature terms in equation (2.3), it can be rewritten as
‘1‘” + F( x ,v¢]v¢|=o (2.9)

where F is the velocity component normal to the interface and called speed function.

Equation (2.9) can be further written into the form of Hamilton-Jacobi equation

%—¢+H(¢,,¢y,¢z,x, y,z 2:) o, (2.10)

where the function H is known as “Hamiltonian”. Basing the techniques from hyperbolic
conservation law, Osher and Sethian [46] designed a class of high order upwind type
schemes, which take advantage of high order ENO schemes developed by Harten et al
[31]. Due to properties of parabolic equation, the second order curvature terms are
resembled at right hand side of H-J equation and discretized by central differential

schemes in these schemes. The first order scheme is

(1551:]: ¢tjk "(Himx (Fty'ktOIV tminiF ijk’ (,OIVJ, (2-11)

16

where

1/2

—x - +x
max ijk ,0)2 +mrn(Dijk ,0)2 +

v” = max(D;{,o)z +min(D;ky,0)2 + (2.12)

max 5,? ,0)2 +min(D,TJ'-"kz ,0)2
L .

 

 

~ *1/2
+x - —x
max(ng ,0)2 + mrn(Dl.j ,0)2 +

V" = max(1);ky,o)2 +min(1)ij‘.ky,o)2 + . (2.13)

ax(D,TJI]f ,0)2 + min(D,-J—-,f ,0)2

 

 

The forward difference and backward difference in x direction are defined as

05% = (¢,~+1 — ¢i)/ Ax and D1}; ¢ = (46,- — ¢,-_1)/ Ax respectively. The formulations

of forward and backward difference operators in y and z direction are similar.

The second order hyperbolic scheme is extended based on first order scheme. The

scheme is the same as first order scheme, however the differential operator V“ and V"
are given by higher order differential approximation. The details will not be discussed

here, and interested readers can refer to [46, 58].

The above schemes were implemented as a test at the beginning of this study. Figure
2-4 and Figure 2—5 are two numerical examples of using hyperbolic schemes on uniform

Cartesian grids. The first example is expansion of an initially star-shaped interface and

17

second example is the motion under curvature flow for an initially seven-point star

shaped curve.

 

 

 

 

0.2:-

0.1:-

>-0-

0.1—

-0.2L
l l rlrrirlr er
-02 -01 0x 0.1 02

Figure 2-4. Expansion of a Seven-Point Star-Shaped Interface, Second- Order
Scheme, 300 by 300 Uniform Cartesian Method, T = 0, 0.03 (0.01)

18

 

 

   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

C a :

0-2': 0.2: b

0.1_- 0.1}
i .

>-0- >0-

-0.1_~ 01?
. t

-0.2:- -02-

l 1 [iii 1 ilii lLr_r l 1 I l

-02 -01 0x 01 0.2 -0.2 -01 0x 01 02
‘ i

0.2} c 0-2.’

01- 0.1-

>-0E- >0}

01} -011

-02} -02}
‘iiliiiiliiiiliiiiliiiilii I'itlirrrliiiiliiiiliiiilii
-0.2 -0.1 0)( 0.1 0.2 0.2 01 0X 0.1 0.2

Figure 2-5. The Evolution History of Star-shaped Interface Under Curvature
Flow300 by 300 Uniform Cartesian Mesh, Second-Order Scheme

2.3.3 ENO and WENO Schemes

Based on the first order upwind differential scheme in section 2.3.1, the high order

differential schemes can be constructed with more accurate approximation of ¢; and (ti; .

The sign of velocity component 11 (v and w for y and z direction respectively) is still used

19

to select g5; or g}; , but 45; and 45; can be greatly improved through essentially non-

oscillatory (ENO) polynomial interpolation.

A general framework for constructing high order ENO schemes for H-J equation
was given by Osher and Shu [47], in which a variety of monotone fluxes were presented
including global and local Lax-Friedrichs flux, Godunov flux and Roe flux with entropy

fix etc. The more details and discussion of ENO schemes can be found in [47, 39].

Another simple way of constructing third order accurate approximation of (if; and

¢; is to use the smoothest possible polynomial interpolation (Newton polynomial
interpolation) to find (0 and then differentiate to get ¢x. The details of construction

second order and third order schemes can be found in [49].

A weighted ENO (WENO) method that takes a combination of all different stencils
was proposed by Liu et a1 [41]. And their work was improved further to obtain the
optimal fifth-order accuracy in smooth regions by Jiang and Shu [34]. Following the
work in [47], Jiang and Peng [35] extended WENO to H-J framework. This H-J WENO
scheme is very useful for solving level set equation and it can reduce the errors by more
than one order of magnitude over the third order H-J ENO schemes. Interested reader can

be referred to the papers mentioned above.

20

Chapter 3. Adaptive Cartesian Grids

This chapter reviews the development of adaptive Cartesian grid and application to
level set method first. And the generation of 2D adaptive Cartesian grid is presented next.
At last several examples of adaptive Cartesian grids for difierent interfaces are given to

demonstrate the robustness and efficiency of our quadtree-based grid generator.

3.1 Introduction

Adaptive Cartesian grids have been used in grid generation and solving fluid flow
problems for decades. In order to obtain more accuracy in regions interested or required,
adaptive grid methods can be employed to save on memory and CPU time by adapting
the areas that are under-resolved. These methods typically rely on recursive subdivision
to generate a computational mesh about complex geometries. With the help of these
techniques, the grid generation can be automatic and extremely fast. It has been also
demonstrated to be very viable for solving inviscid and viscous flow with complex

geometries [16, 17, 36, 53, 72].

In order to optimize the computation for level set method, a number of works [1, 50,

13] have been done by only locally updating ¢ in a band near the interface, which is so
called narrow band level set method. While these methods can reduce computation effort

and CPU time, they still require and inordinate amount of memory. Truly adaptive

21

techniques, referred as patch-based adaptive mesh refinement (AMR) techniques, were
developed in [6, 7] and applied to two-phase incompressible flow in [69]. Much more
efficient tree-based adaptive grid techniques (quadtree in 2D and octree in 3D) were

applied to level set method in [63-66, 24-25, 75] and two-phase flow computation in [42].

To make the level set method more efficient and accurate, adaptive mesh refinement
strategy is ad0pted in this study. Generally grid refinement is desired in regions where the
interface locates, large gradients and high curvature occurs, and the speed function changes
rapidly. To accurately simulate the motion of interfaces, the meshes are adaptively refined
near the zero level set (the interface). In this research, the quadtree-based adaptive

Cartesian grid is employed to implement this strategy.

3.2 Adaptive Cartesian Grid Generation

A quadtree-based adaptive Cartesian grid generator and a grid adaptor have been
developed for the lever set solver based on the Object-Oriented Programming (OOP)
techniques. The OOP techniques make the grid generation, operation of searching
neighboring cells and computation of signed distance very efficient through the recursive
strategy. The detailed definition and background of quadtree-based Cartesian mesh can be

found in [63].

3.2.1 Interface Presentation

22

In this study, the interfaces in two dimensions (open curves intersecting the boundary
or close curves interior) are under the consideration. The interfaces can have more than one
curve, and each separate curve is called a component of the interfaces. Although level set
methods implicitly capture the moving interfaces by find the zero set of level set function,
it is a better idea to create a front data structure to store the interface explicitly. By
borrowing this from front-tracking method, the calculation of sign distance to the interface,
the Lagrangian particle correction and post processing such as volume conservation can be

done easily and efficiently.

A data structure that consists of points connected by elements is employed in this
study. Both particles (the ending points) and elements (linear 2D segments) are stored in
linked lists, which contain the pointers to the previous element and next element. The order
of all the elements is completely arbitrary hence it is not necessary to store the elements
following the actual order of the interface. The advantage of using linked list is that it
makes the addition and removal of particles and elements very simple. The data structure is
shown in Figure 3-1. This data structure can be extended to 3D with similar mechanism
while the elements are substituted by linear triangle, which have three ending points and

three adjacent elements.

23

 

Pointers to particles or end points

 

    
  

Element i

Pointers to adjacent elements

 

 

 

Figure 3-1. Data Structure for ZD Interfaces

The resolution of the interface depends on the discretization described above.
Obviously the interface is more accurate when the elements are finer (typically 3 order
higher resolution than background grids). In this study, with OOP techniques, an object
class named interface is developed and local properties of the interface such as outwarding
normal, curvature can be calculated efliciently and accurately by coupling with the

adaptive grid strategies.

 

Figure 32. Distance from any Point to a Linear Element

24

 

 

 

PTO.

cal

The distance from any point to a linear segment can be easily calculated with a

projection method and dot production, as illustrated in Figure 3-2. The distance L can be

calculated in the following way:

(1) Calculate the square of length BE

2
BE = (xe “bexe 'xb)+(ye — beye —}’b)'
(2) Find the projection of vector BA on vector BE by calculating their dot

production

—-> —>
dot = 3435 = (Xe —bex-xb)+(ye ‘Ybe—yb)-
(3) Calculate the parameter 0 which is belong to [0, 1] by
e = maxIO, minil, dot ”9152 )J.
(4) Find the closest point 0 on vector BE by

1‘0 = Xb +90% 'er
yo = yb +901; we).

(5) Find the distance L to the closest point 0 on vector BE

 

L=\/(x-xo)2 +(y-y0)2-

(3.1)

(3.2)

(3.3)

(3.4)

(3-5)

The signed distance function for any point (x, y) can be directly calculated by finding

the distance between that point and the closest linear segment, which is described above.

The sign can be decided by checking the outwarding normal of adjacent elements. But it is

25

 

 

th

CO

 

 

noted that directly computing the signed distance is extremely time-consuming if one loops
all over the segments. But it can be efficient computed in the (NlogN) efforts on adaptive

Cartesian grids with quadtree-based data structure. This is discussed in Chapter 5.

3.2.2 Quadtree Adaptive Cartesian Mesh

Quadtree-based adaptive Cartesian grids have been developed for the lever set solver
based on the Obj cot-Oriented Programming (OOP) techniques. The OOP techniques make
the grid generation, operation of searching neighboring cells of all the sides and

computation of signed distance very efficient.

A special data structure, called cell-based tree, is used for generation of adaptive
Cartesian mesh. The basic unit of this data structure is cell (or called node) and all the cells
are arranged in a tree. In 2D, regularly each node or cell at level [has one parent at level l-l
(except root) and four children cells at level 1+1 by splitting along both 2: and y direction,
hence it leads to a name of Quadtree mesh. Similarly if each cell has one parent and eight
children cells in 3D, then it is called Octree mesh (Note: root is the cell which has no parent
and its level I is defined to 0, and leaf cell is a cell which has no children). A simple level 3
quadtree mesh and its tree structure are depicted in Figure 3-2 and it has 23 points and 13

leaf cells.

26

 

 

 

 

 

 

 

 

 

LEVEL ROOT

Figure 3-3. Level 3 Quadtree Mesh and Tree Structure

The adaptation including refinement and coarsen can be easily realized by splitting a
cell or simply delete its children on a Quadtree-based adaptive Cartesian grids. Following

some adaptation criterions, desired mesh with high quality can be generated to satisfy

different requirements.

3.2.3 Splitting Rules

27

 

 

dill

COT

 

A graded (also called balanced) adaptive Cartesian mesh is generated under several
different criterions of splitting from a root cell, which represents the complete

computational domain. These criterions are:

1. Split any cell uniformly until size of each leaf cell less than the maximum size
(minimal level).

2. Split any cell whose size is larger than the distance between cell centroid and the
interface.

3. Split any cell whose size is large than the minimum size (maximum level) in a
band nearby.

4. Smoothing whole mesh to ensure that neighbor cells differ in size by no more than

a factor of 2.

Through recursive quad-tree based subdivisions, one can efficiently and
automatically generate fine cells near particular flow features such as a material interface
which can have multiple components. If one is interested in only the interface, the level set
equation can be solved locally and the Narrow Band method is preferred because of its
efficiency. With the Cartesian adaptive mesh, the minimum grid resolution is achieved near

the interface within the narrow band thus high accuracy is guaranteed there.

3.3 Numerical Example;

In this section, several examples of interfaces and adaptive Cartesian grids are given

28

(All numerical experiments here and in later chapters are implemented on a 2.8 GHz

Pentium IV personal computer).

The first example is the generation of adaptive Cartesian grids around a seven-point

star-shaped interface which is defined parameter equation:
y(s) = 20 - [0.1 + (0.065 sin(7 - 2713) (cos 2m, cos 275))1, s 6 [0,1]. (3.6)

For this case, the interface is resolved by 1000 linear elements. The time for generating

such an interface is trivial (only 0.0069 second). The interface is shown in Figure 3-4

 

ITUU

c':
2..
r

 

-0.15;
-0.15 -0.1 -0.05 0)( 0.05 0.1 0.15

 

Figure 3-4. Seven-Point Star-shaped Interface with 1000 Elements

Under such a discretization for the interface, a level 3 to 10 adaptive Cartesian mesh is
generated. Figure 3-5 shows the global adaptive Cartesian grids and Figure 3-6 shows the

local grid after zooming near the interface.

29

 

Figure 3-5. Level 3 to 10 Adaptive Cartesian Grids after Interface-Based Grid
Refinement, Number of Cells = 56764, Time of Generation =0.76 Second

0.1
0.08
0.06
0.04

0.02

 

Figure 3-6. Adaptive Cartesian Grids near the Interface

Another example is adaptive Cartesian mesh generated for interface with multiple

components. This is shown in Figure 3-7. An interface composed of three circular

30

components and each of them is discretized by 300 elements.

.ANOO-P-UIO‘JVm

   

Figure 3-7. Level 3 to 7 Adaptive Cartesian Grids for a Multi-Component interfaces,

Number of Cells =3751, Time of Generation =0.22 Second

31

 

 

ar

1:“

Ci

 

Chapter 4. Finite Volume Level Set Method

A finite volume (FV) formulated algorithm is presented in this chapter for solving
level set equation on two dimensional unstructured adaptive Cartesian grids. First the
introduction of this method is given. The detailed algorithm of finite volume level set
method is described in section 4.2. After that, the finite volume level set method is tested
and applied to several two-dimensional problems in section 4.3 to demonstrate its ability

and accuracy. Summary are presented finally in section 4.4.

4.1 Introduction

Most of the simulations with the level set method described in the introduction
chapter were carried out on structured Cartesian meshes [19, 20, 47, 39, 35]. These
algorithms have been applied successfully to a variety of applications. Since the mesh is
uniform everywhere, it is not possible to have higher resolution for the interface than
elsewhere. This may be very costly for some applications. In order to optimize the
computation for level set method, Narrow Band method has been developed by only
locally updating ¢ in a band near the interface. There also have been a few attempts to use
unstructured triangular meshes for the level set equation. Barth and Sethian [4] developed
such a level set method in 1998. The unstructured grids are much more flexible in handling

complex geometries and adaptation. However, the numerical results may be influenced by

32

 

 

ad

in

he

 

the quality of the computational mesh.

In this study, we extend the level set method to two-dimensional unstructured
adaptive Cartesian grid to improve the solution accuracy and efficiency. A quadtree- based
grid generator is developed for 2D adaptive Cartesian mesh generation. Following ideas
from the finite volume method for hyperbolic conservation laws and data reconstruction
for unstructured grids, a new level set numerical algorithm for the adaptive Cartesian grid
has been developed. A variety of classic two-dimensional numerical examples, including
an accuracy study of the linear wave equation, non-linear Burgers equation, motion of
circular interface under different external velocity fields, propagating surface and motion
under curvature flow etc, are solved to illustrate the capability of the developed algorithm.
Second-order accuracy in smooth regions, good stability and convergence to viscosity

solutions are demonstrated.

4.2 Numerical Algorithm

From sub section 4.2.1 to 4.2.5, the detailed numerical algorithm are described
respectively, including basic formulations for finite volume method, flux computation,

least square method for distance function, data reconstruction, and time integration.

4.2.1 Hamilton-Jacobi Equation

If second and higher order derivative terms in Equation (2.3) are ignored, then it can

33

 

 

 

b

to

\l

 

 

be re-written into the first order Hamilton-Jacobi (H-J) equation:

%+H(5€,V¢)=O. (4.1)

Furthermore, if the curvature term is under consideration, a more general formulation

of a H-J liked equation can be written as

(M - _ - n .

?37 + H(x,V¢,K) _ 0, (x,t) e R x R (4.2)
with the initial condition ¢(5E,0) = ¢0(x°).

In this study, we only study the motion of interface in two-dimension space. In order
to apply the finite volume method to Equation (4.2), one can take the gradient of Equation

(4.2) and the following two equations are obtained:

u +H x, ,u,v,x' =0
{I ( y )x (4.3)

v, +H(x,y,u,v,lr)y =0

. . . . . . “(xaya0):| [u0(x’y)]
where u, v = , , With the rmtral condition = .
( ) I¢x ¢y) [V(x, yaO) V0 (x, y)

Equations (4.3) can be viewed as a hyperbolic conservation law, for which many
numerical methods have been developed in the last three decades, including the finite
volume method. Here the outward normal ii and curvature x in 2D space can be

expressed as

34

_V_¢_ = (they)

'W' Jiiettief

tutti — 22M... + ¢yy¢3
leer)“

ii:

 

(4.4)

KzV-fiz

 

It is noticed that because equations (4.3) are not a strictly hyperbolic system, they may
not be equivalent to equation (4) if u and v are treated as independent variables. Thus at is
still selected as the solution variable, which can be a piecewise polynomial over each cell

of the 2D unstructured mesh. Then the derivatives of the polynomial are taken to obtain

a¢ a¢

u=— and v=—.
6x

4.2.2 Finite Volume Method

In order to find the numerical solution of the hyperbolic system (4.3) on a 2D arbitrary
unstructured mesh, the finite volume method is selected due to its efficiency and simplicity.
Assume Q is divided into general polygons of size h. Each cell is called a control volume
(CV). The solution unknowns are the cell-averages of the level set function ¢. Alternatively,
the cell average ¢, can also be taken to be the point value of (15 at the cell centroid since the
numerical scheme is designed to be second-order accurate. By integrating Equations (4.3)

over CV i, the following equations are obtained:

35

fluidxdy +_”iH(x, y,u, v, k)x dxdy = 0

Hivtdxdy +ILH(x,y,u,v,k)ydxdy : 0 (4.5)

Applying the divergence theorem, the volume integral terms of equation (4.5) can be

transformed into surface integrals

. 6 _
vol(z)--5-t-u,- + an (H -nx)ds = 0 (4.6)
and
vol(i)-2Tfi + (H-n I15 = 0 (4 7)
6t ' aVi y ° °

where the vector (nx, ny) is the normal of each face of cell i, vol(i) is the volume, and 17,.
and V, are the cell averages of u and v respectively. These variables and their arrangement

are illustrated in Figure 4-1.

 

 

 

Inpnyj

Cell i

Figure 4-1. Schematic of a Polygon Control Volume (CV)

36

Once the solution unknowns ¢,, :7, and 17,- at time step n are given, they are used to
build piecewise polynomials in each CV to achieve second-order accuracy. No solution
continuity is enforced at the cell interfaces, therefore the fluxes through these interfaces are
not uniquely defined because different solutions exist on both sides of the interfaces.
Similar to the 2D Euler equation, Riemann numerical fluxes are employed for equations

(4.6) and (4.7). Thus (4.6) and (4.7) can be rewritten in the following flux forms:

ne
v01(i)-[L-‘i+ Zfilksik = 0 (4.8)
dt k=1
and
ne
voz(i).-d—‘—"'—+ grins”, = 0. (4.9)
dt k=l

For the 1"” cell, Si, denotes length of the k’h face, NE is the number of faces, and the
numerical fluxes [ilk and I?“ are defined at the k’h face. The numerical fluxes [in
and Flu are Lipschitz continuous monotone fluxes which are consistent with function
H (:e,v¢, K) in Equation (4). Details of the Riemann flux can be found in the paper of

Osher and Shu [47].
4.2.3 Flux Computation

Because the polynomial for each cell is constructed locally within a stencil, the left

37

solution (u‘,v’) and the right solution (u+ ,v”) of a cell face are different. A simple local
Lax-Friedrichs (LLF) flux is used for the numerical tests presented in this paper. Given the
left and right solutions, LLF flux is computed and has less numerical dissipation than the

global Lax-Friedrichs flux. The LLF flux is defined as
‘ + - + — _1 + + — — 1 + -
Hlku ,u ,v ,v —5Hu ,v +Hu ,v ~nx—Ea-u —u (4.10)

and

H2k(u+,u‘,v+,v‘)=%[H(u+,v+)+H(u‘,v’)]-ny ---;-,B-(v+ —v‘), (4.11)

where parameters a and B are defined as

6H(u,v) .

av (4.12)

and ,6 = maxu’v

 

If the velocity is imported from an external known velocity field i7 , then we have

6H(u,) V =an|Vx| (1611. (,u av,)=V| IVI

4.2.4 Linear and Quadratic Data Reconstruction

In the finite-volume level set solver, solution variables ¢ , u and v are defined only at
the centroid of a control volume. In order to evaluate the fluxes through cell faces, the
solutions variables at both sides of the face are needed to be calculated. This can be

fulfilled by data reconstruction. In this study, two reconstruction approaches are developed

38

and implemented.

In the first approach, a simple least square approach is employed to construct a

piecewise linear polynomial over each cell. Let Q represent one of the solution variables ¢ ,

u and v, then the gradient of Q over each cell can be calculated through the least square

approach. A linear polynomial is constructed for each control volume (CV) as below:

Qi(x’y)=§i+Qx'(x—xci)+Qy'(y—yci) (4-13)
Here (xci , ya) is the coordinate of centroid of the ith cell. The square of error over all the
cells in the stencil is defined as

5:22}: [0,-+0..-(x—x..)+Q.-(y—y..)—0,-f (444)
J. .

J

=0 and-fl—=0, a set of

an aQy

 

By differentiating Equation (4.14) and letting
equations are obtained as

Qy'ny'I'Qx'Jxx—JQI =0

. (4.15)
Q, Jyy +Qx .ny 4Q, =0

where

39

Jxx = 2(ij —xci)2

J

J”, =Z(ycj 'yciy
j

ny =Z(xcj ‘xcixycj _}’ci)
1

JQx =Z(_Q-j Taxxcj ”xci)

J
JQy = Z(§j T 51'ij ’J’ci)
1'
By solving equations (4.15), the gradients are calculated by

Qx = Jyy-JQx—ny 4Q,

 

 

JxxtJW—Jfiy
Jn-ng—nythx'
y Jn-Jyy—ny

(4.16)

(4.17)

For the first order scheme, one-time linear reconstruction is sufficient. In order to

(face neighbors) are used in a relatively small reconstruction stencil.

construct second order scheme, a linear reconstruction is performed firstly for the signed
distance function ¢ to find the gradient V¢, i.e., u and v. Then given u and v, another

linear reconstruction is carried out to obtain the second-order derivatives ¢xx , (15,0, and
(15”,. Both reconstructions are performed with the least-squares approach, and can be

implemented efficiently through using a face-based data structure for 2D arbitrary

unstructured mesh. In the linear reconstruction, only neighbor cells sharing the same face

In the second approach a cell-wise quadratic reconstruction designed by Delanaye etc

40

[21] is performed to construct second-order scheme directly. A larger reconstruction stencil
is necessary to build a quadratic polynomial over a cell. In addition to face neighbors used
in linear reconstructions, the cells sharing the same node with the current cell (node
neighbors) are also chosen in the quadratic stencil. Figure 4-2 shows a typical stencil used

to perform quadratic reconstruction in 2D adaptive Cartesian mesh.

 

 

 

 

 

 

Cx /\ HR {ix
/ \J K/ \J
G p o
C’\ KN {‘N K‘\
J \J \J \J
n
r C
I
m
C'\ f\ f'\ K'L
/ \J \J \T
.I k
S I
c o c e

 

 

 

 

 

 

Figure 4-2. Quadratic Reconstruction Stencil of Cell i for Adaptive Cartesian Grids

In order to construct a second order polynomial over the cell i, let us consider the

Taylor expansion of the level set function 05 j (x j, y j) (j is the index of cells involving in

the stencil) about the cell center (xi, yi) by omitting third and higher order terms:

¢j =¢(x,-.yj)=¢.+¢. ~4x+¢y -Ay+-§-I¢..-(At)2 +2¢.yAx-Ay+¢yyi(4yfi
(4.18)

where Ax = x j -— x,- and Ay = y j — yi . Assume N j is the number of cells involved in the

41

stencil supporting the reconstruction and regularly bigger than number of unknowns.

Substituting the coordinates of the centroid of cell j into Equation (4.18), the following

system is obtained

where A is a two dimensional matrix

and

 

 

Aj1=xj
Aj2=)’j—)’i
Aj3=(xj—xiXxj—xi)/2
Aj4 =ij wily,- —y,-)/2
Aj5=(Xj—xiX)’j—yi)

A¢j=¢j’¢i-

(4.19)

(4.20)

(4.21)

Because Equation (4.19) is an over-determined equation system, it can be solved

through least-squares method. Singular value decomposition (SVD) method is employed

to compute the pseudo inverse of matrix A, which is denoted by A+. Then the first and

second order derivatives of 05 can be solved directly through the simple matrix operation

42

A+ -A¢.

Both approaches are tested and desired order of accuracy is achieved. The first
approach is preferred because a small stencil is sufficient. In the process of reconstruction,
large error occurs near the area with discontinuity. If the solution variables are smooth
enough, data limiting is not necessary since it is second-order accurate already. However,
for non-smooth problems with discontinuities, a limiter similar to that developed by

Barth21 is employed to avoid or reduce numerical oscillations.

4.2.5 Time Integration

Equations (4.8) and (4.9) can be rewritten in the form of ordinary differential equation

(ODE)
fig. = 1.64:5.) (4.22)
dt 1 l ’
where
j=nek
21H t'jStr
". = __J_=__
res(Q, ) vol (1) (4.23)
and
Q- = [If]. (4.24)
Vi

43

For time discretization, a simple two-stage Runge-Kutta (R-K) time integration

scheme is applied to solve Equations (4.22), i.e.

Q0) =5," +At-res(§,-")
— 2 1— 1—1 1 —1
Qi( )=§Qin+§Qi()+§At’reS(Qi( ))° (4-25)

Qn+l = Q0.)

The two-stage Runge-Kutta time integration scheme can achieve second order time
accuracy. Since Equations (4.25) will only update the gradients V¢ , we need to update the

solution variable ¢ by solving Equation (4.2) and requiring that

[an (96? + H(5e, V¢))dxdy = 0 . (4.26)

This updating procedure can be coupled with solving Equations (4.25).

Since the scheme for V¢ only determines 96 for each cell up to a constant, another
way to update 45 in whole computational domain is to calculate the path integral by using

the following equation starting from a particular point:
- - f2
¢(x2 , z) = ¢(xl,t)+ I521 (¢xdx + ¢ydy). (4.27)

The integral path should be taken to avoid crossing a derivative discontinuity;
otherwise another starting point needs to be selected. Both approaches are implemented in

our numerical experiments. For smooth problems, the results of both approaches are

44

similar and same order of accuracy is achieved while the first approach gives slight better
results. However, for unsmooth problems with discontinuity and singularity in and its
derivatives, the first approach often produce dents. By using the second approach one
could choose the integral path to avoid crossing derivative singularities. Therefore, the
smooth problems including linear wave equations and circle expansion in Sect. 4 are
computed with first approach while the others with discontinuity are solved through the

second approach.

4.3 Numerical Results

In this section, the finite volume level set method described in the previous sections is
tested and applied to several two-dimensional problems. Firstly, an accuracy study is
carried out for the linear wave equation. Then 2D Burgers equation is solved to test the
stability of the algorithm for non-smooth problems with discontinuities. After that, the
motions of an initially circular interface are computed. A grid adaptation procedure is
tested on the adaptive Cartesian grids. Finally, more sophisticated problems evolved with
outward normal ii and second-order terms (mean curvature K) are solved to demonstrate
the ability of developed numerical method. All numerical experiments are implemented on

a 2.8 GHz Pentium IV personal computer.

4.3.1 Accuracy Study with 20 Linear Wave Equation

As the first test case, the 2D linear wave equation is chosen to study the numerical

45

order of accuracy of the finite volume method. The 2D linear wave equation takes the

following form
¢t +2. +¢y =0. (tape-[0.21402] (4.28)

with initial condition ¢(x, y,0)= sin(7r(x+ y)) and periodic boundary conditions. The

Hamiltonian satisfies thatH(u,v) = u + v. The exact solution of this example is
¢(x, y,t) = sin[:2r- (x + y —— 20). (4.29)

For this accuracy study, the computations were carried out on uniform Cartesian

meshes. Both the first order and second-order schemes were implemented. The initial (15 at
t = 0 are shown in Figure 4-3 and ntunerical solution of (b are shown in Figure 4-4. The
grid refinement studies are carried out for this case and the L1, Loo errors and orders of

accuracy are presented in Table 4-1. Note that the designed first order and second order of

accuracy are achieved.

46

 

Figure 4-3. initial 45 = sin(7r(x+ y)), T :00

 

Figure 4-4. Numerical Solution of ¢, T =1 .0, Level 6 Uniform (64 by 64)

47

Cartesian Mesh, Second-Order Scheme

Table 4-1. Accuracy Study for the ZD Linear Wave Equation, T =1.0

 

 

 

 

 

 

 

 

 

 

 

 

3:31:22: Grid Level L,o Error Lao Order L1 Error L1 Order
16* l 6 4 7.08e-01 --- 4.56e-01 ---
32*32 5 4.60e-01 0.62 2.94e-01 0.63
1st order 64*64 6 2.65e-01 0.80 1.69e-01 0.80
128*128 7 1.43e-01 0.89 9.10e-02 0.89
256*256 8 7.42e—02 0.95 4.72e-02 0.95
16* 1 6 4 8.19e-02 --- 5.34e-02 «-
32*32 5 2.03e-02 2.01 1.30e-02 2.04
2nd order 64*64 6 5.06e-03 2.00 3.23e-03 2.01
128*128 7 1.26e—03 2.01 8.04e-04 2.01
256*256 8 3.16e-04 2.00 2.01e-04 2.00

 

 

 

 

 

 

 

 

4.3.2 20 Non-Linear Burgers Equation

To demonstrate the ability of the finite volume level set method in solving
non-smooth problems, the non-linear two-dimensional Burgers equation is solved. The

same initial condition as in Hu and Shu [33] is selected

 

+ (¢x + ¢y + 1)2 = O, (x, y) 6 [—2,2] x [—2,2]
2 (4.30)

¢<x,y.o)=-cos[-’,i(x+y))

¢z

with periodic boundary conditions again. In this case, the second-order scheme is

48

employed. Since the discontinuity exists, a Barth limiter is applied in the data
reconstruction to avoid numerical oscillations. The smooth solution at t=l/7r2 and

non-smooth solution at t = 2 / 7:2 with discontinuous derivatives are shown in Figure 4-5.

Note both solutions are non-oscillatory at least to the naked eye.

 

(b) t=2/7r2

Figure 4-5. Numerical Solution of 20 Non-linear Burgers Equation, Level 6 Uniform

(64 by 64) Cartesian Mesh, Second-Order Scheme

49

 

4.3.3 Expansion of Initially Circular interface

In this case, expansion of an initially circular interface is simulated. Equation (4.1) is

solved with unit normal velocity V" =1 . Then the level set equation becomes

a
£+J¢f+¢§ :0. (431)

For this case, the initial interface is a circle with a radius of 2. To analyze the accuracy,
the initialization of the distance fimction ¢ is exact to avoid introducing extra numerical
errors by calculating the value ¢ = '55 —J‘cO|-2 for each grid point. The exact boundary
conditions are imposed and a second-order scheme was employed. Table 4-2 presents the
numerical order of accuracy of the solution for -1 S ¢ 51 on uniform Cartesian meshes. It
is noticed that the desired order of accuracy is achieved while CPU time per time step is
proportional to number of total cells. Figure 4—6 displays the evolution history of the

interface from t = 0 to t = 4 with interval of 0.5 second.

The computations with adaptive refinement and a narrow band procedure are also
carried out, and the numerical results are shown in Figure 4-7, 8 and 9. The maximum level
of the cells is 7 (assuming the root cell being level 0) while the minimum level is 3. The
smallest and largest cells in the computation differ by a factor of 16 (4 levels). Note that
high resolution of the interface was achieved near the interface through interface based grid

adaptations with much fewer number of grid cells. Similar accuracy is achieved while the

50

 

total of CPU time is only one seventh of level 7 uniform Cartesian meshes.

Table 4-2. Accuracy and CPU Time for Circular Interface Expansion, T=4.0

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Order of Grid L.o Order L, Order CPU

accuracy Time/Step

16*16 4 5.41e-02 -- 1.39e-02 --- 1.80e-3 s

32*32 5 8656-03 2.64 2.56e-03 2.44 6.756—3 s

2"" order 64*64 6 2.81e—O3 1.62 5.818-04 2.14 2.668-2 s

128*128 7 9.11e-04 1.63 1.35e-04 2.11 1.14e—l s

256*256 8 2.90e-04 1.65 3.18e-05 2.09 4.456-1 s

L l l l l l I l I l I l l I l
-4 0 4 8

Figure 4-6. Expansion of Circular interface, V" =1, T= 0,4 (0.5), Level 6 (64 by 64)

 

 

 

 

Uniform Cartesian Mesh

51

 

 

 

‘88 4 0 4

Figure 4-7. Expansion of Circular Interface: Initial Interface, T = 0, 1228 Cells

 

‘8-8 .4 0 4 8

Figure 4-8. Expansion of Circular Interface: interface at T = 2, 2392 Cells

52

 

 

 

 

Figure 4-9. Expansion of Circular interface: Interface at T = 4, 3433 Cells

4.3.4 Moving Circular Interface under Passive Velocity Field

Two kinds of motions under different passive velocity fields for an initially circular
interface are simulated in this section. The first one is the convection of interface with a
constant velocity V = (1, l) . Thus the solution is just a circle translating from the
bottom-left to the top-right and the shape of circle should remain the same. The exact
boundary conditions are imposed at x=0 and y=0. Simple extrapolations are applied at the

other two boundaries. The second-order scheme is applied in the first simulation.

The numerical experimental over both level 6 uniform and level 6 to 9 adaptive
Cartesian meshes are implemented and the results are shown in Figure 4-10 and 4-11,

respectively. Figure 4-12 shows the adaptive Cartesian mesh at T = 4. Note that the solution

53

 

on the adaptive mesh is as good as that on level 9 (512 by 512) uniform Cartesian mesh,

and was obtained with much fewer cells. The CPU time of level 6 to 9 adaptive Cartesian

meshes is only one tenth of that of uniform meshes with same minimum level.

 

8

k
I U—Y I I

 

”A

“7

 

 

00

l 8

Figure 4-10. Convection of Circular interface: T=0, 4(1.0), Level 6 (64 by 64)

Uniform Cartesian Mesh, CPU Time/per step = 2.35e-2 second

 

8

TV I

 

2‘3

 

 

I 8

Figure 4-11 . Convection of Circular Interface: T=0, 4(0.8$), Level 6 to 9 Adaptive
Cartesian Mesh CPU Time/Step is 1.38e-1 sec

54

 

 

 

Figure 4-12. Level 6 to 9 Adaptive Cartesian Mesh at T = 4, 23338 cells

In the second case the velocity field is set to I7 =(y,—x). This velocity field is
time-independent and leads to a solid body rotation of a circle of radius 0.5 initially
centered at (-l, 0). Both first order and second order schemes are tested in this case. The
computation is carried out on both uniform and adaptive Cartesian meshes. Figure 4-13
shows the numerical result on level 7 (128 by 128) uniform Cartesian mesh with the first
order reconstruction scheme. Base on the result, the numerical solution is pretty diffusive
and conservation of volume is not satisfied well. By applying the second order scheme to
same grids, the results are much better (seen in Figure 4-14). After that, adaptation
technique and second order scheme are employed on level 5 to 9 adaptive Cartesian
meshes to make comparison. With even less number of total cells (around 15000) than

level 7 uniform mesh (16384), the interface is captured with higher resolution, higher

 

 

accuracy (similar as that on level 9 uniform mesh) and ahnost same CPU time as level 7

uniform mesh. The results are shown in Figure 4-15.

2

 

 

 

 

2 rlurgulgrrr
-

l I 1 1 I I I l
-2 -1 0 1 2
Figure 4-13. Rotation of Circle: First-Order T = 0, 7: (7r/ 6), Level 7 Uniform

 

Cartesian Mesh

 

 

 

 

 

2
I.

1-

i

0..

-L

1L
hilljlllllljlll'llll

'2-2 -1 0 1 2

Figure 4-14. Rotation of Circle: Second-Order, T = 0, 7: (7r/ 6), Level 7 Uniform

Cartesian Mesh

56

 

 

 

 

 

2
1-
0-
l'
-1-
-2 IJIILIILIIIIIIIILL
-2 -1 0 1 2

Figure 4-15. Rotation of Circle: Second Order, T = 0, II (n/ 10), Level 5 to 9
Adaptive Cartesian Mesh

4.3.5 Propagating Surface

The motion of a propagating surface is computed in this section to demonstrate

correct “viscosity solution”. The governing equation for the propagating surface is:

Z—‘t — (1 — ext/1+ ¢f + as; = 0, x 6 [0,1], y 6 [0,1] (4.32)

t

with a smooth initial condition
¢(x,y,0)=1.0 —:11-(cos270c —1)(cos27zy — 1), (4.33)

where a is a small constant and K is the mean curvature defined by

57

K: ¢..(1+¢§)—2¢.¢y¢., +¢yyit+¢ii (4.34)

/2
(l+¢3+ 5

 

In order to make comparison with results in Hu and Shu [33], the computation is
implemented on uniform Cartesian mesh only. A level 6 (64 by 64) uniform Cartesian mesh
is selected and periodical boundary conditions are applied. In the case of 3 =0,
discontinuous derivatives appear, data limiting is applied in the reconstruction procedure,
and the convergence to a viscosity solution is shown in Figure 4-16. The smooth solution
for e=0.1 is shown in Figure 4-17. This case demonstrates the ability of the numerical

method to solve curvature-dependent problems.

 

 

Figure 4-16. Propagating Surface, 3:0, T =0, 0.6 (0.3) (Note: ¢—0.35 is

displayed for visual effects at T: 0)

58

 

 

Figure 4-17. Propagating Surface, e=0.1, T=0, 0.6 (0.3) (Note: ¢—0.35 is

displayed for visual effects at T: 0)

4.3.6 Star-Shaped Interface

In this section, the motions of an initially seven-point star-shaped interface defined by

parameter equation
7(5): 20 - [0.1 + (0.065 sin(7 - 271's) (cos 2715', cos 2m))], 5 6 [0,1]. (4.35)

are simulated. These cases were considered by S. Osher and J.A. Sethian [46], and are
selected here to demonstrate the capability of our method to handle more complex

interfaces and make comparison.

59

Two types of motions are simulated in this section. First, it is assumed that the normal
velocity component VN equals to unit. Since discontinuity occurs and sharp comers exist
in this case, data limiting is necessary. In order to reduce the CPU time while keep the
similar accuracy as in Osher and Sethian [46], level 4 to 8 adaptive Cartesian mesh are
employed to do simulation. To avoid dent and bump, four comer cell and one center cell
are used to update the distance function 45. The non-smooth solution of the interface,
which is converged to a viscosity solution, is shown in Figure 4-18. The numerical results

agree well with simulations presented in the literature.

Next the speed function VN of the interface is set to - K' , where K is curvature
defrne by Equation (4.4). This curvature-dependent motion of the interface is simulated
over level 8 (256 by 256) uniform mesh in order to compare with literature. During the
computation, the absolute value of curvature is limited to 1/Ax, where Ax is the size of
minimum cell. Figure 4-19 shows the smooth numerical solution of the interface. The

result matches that in [46] well.

60

 

 

 

 

 

Figure 4-18. Expansion of Star-Shaped Interface, V, =1, T =0, 0.07 (0.01 ),Level 4
to 8 Adaptive Cartesian Mesh '

 

 

 

Figure 4-19. Motion of Star-Shaped interface, VN = —x, T =0, 0.005 (0.001), Level

8 Uniform Cartesian Mesh

61

4.4 Summary

A finite volume method for the level set equation has been developed for arbitrary 2D
polygon grids including adaptive Cartesian grids. Two different least squares
reconstruction approaches were employed. One is linear and the other is quadratic. Both
approaches gave satisfactory numerical results Accuracy studies have been performed with
the linear wave equation, and the present finite volume method was capable of achieving
the designed numerical order of accuracy. Interface based grid adaptations have been
demonstrated to capture the interfaces with high resolution. A variety of test cases have
been used to test the finite volume level set method, and results have been satisfactory. It
can be extend to three dimensional through a Octree-based adaptive Cartesian mesh to
tackle more general real interface problem and couple the level set solver with a flow

solver to solve multi-phase flow problems.

62

Chapter 5. Semi-Lagrangian Level Set Method

In this study, semi-Lagrangian (SL) level set method is implemented and extended
with a particle correction algorithm on 2D unstructured adaptive Cartesian grids. Firstly
the semi-Lagrangian level set method is introduced in section 5.1. Then the detailed
algorithm of semi-Lagrangian level set method is described in section 5.2. After that in
section 5.3, the semi-Lagrangian level set method is tested and applied to simulate the
various motions of the interface under different velocity fields. Summary are drawn last in

section 5.4.

5.1 Introduction

Since Osher and Sethian [46] proposed the level set method in 1988, there have been a
number of schemes developed, which has been mentioned in introduction chapter. Among
them, the fast semi-Lagrangian level set method is quite simple yet very accurate, efficient
and specially designed for quadtree-based adaptive Cartesian grid by Strain [65]. In this
study, a 2D semi-Lagrangian level set solver is developed based on the tracing-back
algorithm presented in [63]. Quadtree-based adaptive Cartesian technologies are used. A

particle correction algorithm is also employed in order to achieve better mass conservation.

63

5.2 Numerical Algorithm

From sub section 5.2.1 to 5.2.7, the detailed numerical algorithm of semi-Lagrangian
level set method is described, including first order and second order numerical scheme,
Whitney extension of velocity, redistancing algorithm for maintaining the signed distance

function, contouring algorithm for zero level set, and coupling of particle level set method.

5.2.1 Numerical Schemes of Semi-Lagrangian Level Set Method

Starting from a Lagrangian view, one can observe that level set equation (2.2)

propagates level set function ¢ along the characteristic curve 56 = 8' (t) which is defined

by

2(2): §(t)= V6.3) (5.1)

and initial condition given by x0 = 8(t0). Hence the value of level set function ¢ at time
step n+1 (t + dt) can be computed by finding the characteristic curve passing through point
(55, t + dt) and following it backwards to some previous point (Y, t), where the value of ¢
is already known at the time t, then setting ¢"+1(x‘,t+ dt) = (2" (it). This procedure can
be numerically realized by solving the characteristic ordinary difference equation (ODE)
(5.1). Once the velocity field l7(x,t) is extended from the interface or computed by
external physics, it is easy to locate the off-grid point (3?, t), on which the signed distance

function ¢(Y,t) can be evaluated by numerical interpolation or directly calculating to the

discretized interface with the help of an auxiliary distance tree. That will be described in
section 5.3. The simplest way of calculating E is to use the first order accurate

approximation

3‘ =55 —dt . 17(r,t), (5.2)

and linear interpolation of ¢(3c', t) can calculated by bilinear interpolation

atrttl= (1 -aX1 —/t)¢oo + all who. (1 -a),6¢m watt“, (5.3)

where (1500, (1501, P11, (1510, a and ,6 are defined and showninFigure 5-1.

 

 

 

 

¢01 ¢11
dx
a:—
Ax
-------- of
4 dx 9? Ax
= -9
dyEIB-Ax
V
¢00 Ax ¢IO

Figure 5-1. Bilinear interpolation of Off-grid Value

The first-order semi-Lagrangian scheme for solving the convection function (2.2) is
called Courant-Issaacson-Rees (CIR) scheme [18], which was developed in [63] and

applied in [64-66]. The formulation of CIR scheme is

65

¢(5c',t+dt)= ¢(52 —d:-i7(r,z), t). (5.4)

In order to obtain second-order accuracy in time direction, a CIR predictor (5.3) and a
second-order trapezoidal corrector are combined together. Firstly a predicted function a

is acquired by using the CIR scheme, i.e.,
$65,: + dt)= ¢(3c', t) = ¢(55 — dt - 17(x,:), t). (5.5)

Then the extensive velocity I7 is evaluated from the predicted function $6131 + dt)

at time t + dt. The second-order semi-Lagrangian scheme is given by

¢"+1(x°,t + dt) = ¢"[i _%{,,7(;,,)_5:_’.i7”(3,, + dtl ti (5-6)

 

 

 

 

 

Figure 5-2. Fast Semi-Lagrangian Level Set Method

The fast semi-Lagrangian level set methods are illustrated in Figure 5-2. The CIR and

second-order semi—Lagrangian schemes are explicit yet unconditionally stable and satisfy

66

the CFL condition. The solutions can converge with large time step ( dt = 0(Ax) ) even
for parabolic problems where the velocity is depended on curvature terms. Compared to
second-order semi-Lagrangian scheme, the CIR scheme is only first-order accurate in time
direction and sometime dissipative. In this study, both CIR scheme and second-order
scheme are tested and numerical results shown in this chapter are computed by second

order semi-Lagrangian scheme if there is no extra explanation.

5.2.2 Fast Redistancing Algorithm

As stated in Chapter 3, the interface is discretized with N linear segments and the
signed distance from any point to the closest linear element can be calculated by a
projection method but it will be extremely time-consuming if one searches through all
linear elements. With the help of distance tree introduced in [64], fast redistancing

algorithm of maintaining signed distance function is efficient and accurate.

Firstly the sign of distance function need to be fixed. The sign of 45 can be fixed
efficiently by defining an outward unit normal vector for every element of interface F and
checking which side of the closest element it lies on. Two special cases of fixing sign of
¢ need to be considered carefully, and they are shown in (a) and (b) of Figure 5-3 for
intersecting elements which form sharp comers. In the first case, 0 is the closest point
hence both elements are nearest segments to point r . The point 55 is out side of element 1'

and inside of element i+l, while the sign of ¢ is positive. For the convex case (a), the

67

sign is negative only when the point It falls inside of both element 1‘ and element 1'+1,
otherwise the sign is positive. For the concave case (b), the sign is positive only when x

falls outside of both element 1' and element i+l, otherwise the sign is negative.

 

Convex Case Concave Case

Figure 5-3. Fixing the Sign of Level Set Function:

(a) Convex interface; (b) Concave interface

Distance tree is a special quadtree mesh and resolves interface F at optimal 0(N) cost.
Meanwhile it allows a fast evaluation 0 (NlogN) of the signed distance firnction for I‘. The
concentric triple of any cell C shares same center and the size is three times of cell C. The

illustration of concentric triple is shown in Figure 5-4.

68

 

 

 

 

 

 

 

 

3b

Figure 5-4. Concentric Triple of a Quadtree Leaf

To generate such a distance tree, one can follows the following two rules:

a) Split any cell whose size is larger than the distance between cell and the
interface;

b) Split any cell whose concentric triple intersects F.

One example of the level 7 distance tree for an interface composed of multiple circles

is shown in Figure 5-5.

69

NQAO‘IONCD

 

Figure 5-5. An Example of Level 7 Distance Tree for Multiple Circlers

 

 

 

 

 

 

 

 

 

 

 

‘ F
ii
[H/ T. Keir/A
T** C***
T***

 

 

 

Figure 5-6. Fast Searching Strategy

70

In order to evaluate the signed distance fi'om a vertex of cell C to the interface, a fast
searching strategy is developed in [64] and employed in this study. This strategy builds the
distance tree in 0 (NlogN) time and space complexity. Any element intersecting T* beats
all the elements outside of T***, which is illustrated in Figure 5-6. This property will
exclude most of the elements and greatly reduce the computation efforts hence the
searching always terminates in a bounded number of steps. Interested readers can refer to
the work done by Strain [64]. The fast redistancing algorithm is accurate and efficient on
tree-based adaptive Cartesian meshes, the evaluation of signed distance is straightforward

and the accuracy is dependent on the resolution of interface.

5.2.3 Whitney Velocity Extension

For the motion of the interface, theoretical only the velocity of the interfaces
themselves contains real information of interface propagation. Since the velocity field off
the interface has nothing to do with the motion of the interface, how to extend velocity on
the background grids is a very important issue for level set method. What kind of extension

velocity field has the desirable properties?

Since the original level set method [46], there have been various approaches
developed to extension velocities. The simplest one is to directly use the velocity field
solved by external physics (for example fluid mechanics) as the extension velocity, as done

in [56, 10, 67]. Another approach is to extend the velocity from the interface, including a

71

boundary integral model extending velocity through jump condition across the interface
[57], the extension algorithm through solving the partial differential equations by Chen,
Merriman, Osher and Smerkeka [1 l ], and the Whitney velocity extension in [65, 66] which

is employed in this study. Interested reader can refer the details in Chapter 11 of [5 8].

By defining the velocity 17 at any point r as the interface velocity of the closest
point on the interface SEC is an extremely smart idea. By doing this, it makes the velocity
off the interface approximately constant in the normal direction near the interface. As
shown by Zhao et al [80], using the closest interface point velocity to convect the interface
can maintain the signed distance well. This method was successfirlly employed by Chen et

al [11], and Adalsteinsson and Sethian in [2].

In this study, the method of numerical Whitney extension is employed. This method
couples the classical Whitney extension [77] and fast searching method of Quadtree-base
adaptive Cartesian grid. The velocity field is only extended near the interface and the level
set function is only updated locally. By solving the level set function locally, the distance
tree speed up the evaluation of extensive velocity greatly because these points lies in the
finer cells with only few nearby linear elements of interface, where the searching operation

is extremely fast. Readers can refer to the worked done by Strain [64].

If the velocity of the interface is dependent on outwarding normal ii and curvature

K (also called local geometric velocity), the higher order ENO or WENO schemes can be

72

implemented locally near the interface since the cells there are finest yet uniform. In this
study, the second and third-order ENO schemes are employed to evaluate the first-order

and second-order derivatives of level set function ¢ locally. Then the normal and

curvature are computed by

ii=—V£— I¢x+¢yI and K=V fi=22£+an—y. (5.7)
IV¢I ¢3+§ 6x 6y

5.2.4 Contouring Algorithm

Since the level set function (15 is an implicit representation of interface, only the zero
isocontour is considered and need to be extracted given (0 on the adaptive Cartesian grid.
The most straightforward way is to use a contour plotting algorithm to locate the zero
isocontour (curves in 2D and surfaces in 3D) and discretize them by linear elements (linear

segments in 2D and linear triangles in 3D).

In order to implement such a fast contouring algorithm on adaptive Cartesian grid, an
auxiliary contour tree is generated then triangulated. The contour tree is a balanced
quadtree, in which the adjacent cells difi‘er in size at most a factor of 2. For such a contour
tree, the vertices and cell center can be triangulated by checking size of neighboring cells.

An example of contour tree and correspondent triangulated tree mesh is shown in Figure

5-7 (a) and (b).

73

Nwhmmflm

wammflm

 

Figure 5-7. Level 6 Contour Tree and the Triangulation

Once ¢ at the vertices and cell centers of the triangulated quadtree mesh is given, it
is straightforward to plot the zero isocontour by extracting the zero set of the continuous
piecewise linear interpolant on the triangulation. Since the triangulated quadtree mesh is

conforming and the interpolant is continuous, the zero line segments form a polygonal

74

interface. Based on above linear contouring algorithm, a fast adaptive contouring method
was developed by Strain [66]. The fast adaptive contouring method uses bisection,
adaptive refinement, and prune algorithm, and it can capture the sub-grid details of the
interface. The details of this method can be found in [66]. In this study, adaptive contouring
algorithm is implemented to plot the interface in order to acquire better resolution of the
interface. Although the fast semi-Lagrangian contouring algorithm is efficient through
taking the advantages of adaptive Cartesian mesh, it can be time-consuming and even fail
to plot the interface correctly for some extremely difficult situations (for instance, the
vortex problem in section 5.3). In this study, a particle correction algorithm is developed on
quadtree-based adaptive Cartesian grids to further improve the accuracy and robustness,

which is discussed in next section.

5.2.5 Particle Correction

Particle level set method was firstly developed by Enright et a1, and in this method a
layer of Lagrangian marker particles near the interface are randomly distributed near the
interface. The particles are then moved with high order Runge-Kutta schemes and used to
correct the interface in under-resolved area and reinitialize the signed distance function.
This hybrid method couples the H-J WENO scheme and Particle correction algorithm, and
it dramatically improves the conservation of the mass. Based on it, Enright, Losasso and
Fedkiw developed a fast and accurate semi-Lagrangian particle level set method, in which

the level set method are solved by CIR scheme and the interface same particle correction

75

algorithm are employed. Since the CIR scheme is much simple than high order WENO H-J

scheme, the method is more efficient. Compared to these two methods, a more efficient

particle correction algorithm is developed in this study. The particles are treated as the

vertices of linear elements, which is stored with a linked list data structure described in

Chapter 3. Much less number of particles are needed to move with high order Runge-Kutta

method compared to the above two method since the particles are only placed on the

interface itself (the particles will represent the interface). By using the linked list, the

removing unnecessary particles and adding new particles are quite simple. With the fast

redistancing algorithm to evaluate the signed distance function, the signed distance can be

maintained accurately and efficiently.

The algorithms of particle level set method in this study are:

(1)
(2)

(3)

(4)

Evaluate the signed distance function with first-order CIR scheme;

Moving the particles with second or three-order Runge-Kutta schemes by using
the extensive velocity field;

Use the particles and elements stored in linked list to correct the level set
function locally (within a band of interface);

Update the particles and elements and check whether adding particles and

removing particles are necessary.

The particle correction algorithms are tested in study. The numerical results are

76

amazingly good and comparison to regular semi-Lagrangian level set method without

particle correction will be shown in next section.

5.3 Numerical Results

Versatile motions of interface under passive transport velocity field, first order and
second order geometric velocity field are simulated with semi-Lagrangian level set method
to show its ability. The particle correction algorithms are also tested and the numerical
results are compared to regular semi-Lagrangian level set method without particle

correction.

5.3.1 Zalesak’s Disk Rotation Problem

The rigid body rotation of Zalesak’s disk (notched circle in 2D) under a constant
vorticity velocity field [78] is employed as the first example to demonstrate the ability of
semi-Lagrangian level set method to conserve mass. This case has been studied by
Sussman [68], Enright and Fedkiw [24], and Enright, Losasso and Fedkiw [25] for the
purpose of improving mass conservation of their numerical schemes. It has been known
that the classic level set methods are diffusive and tend to lose mass while the hybrid

level set method in [24, 25] can achieve much better mass conservation.

The initial interface and velocity field over the computational domain are shown in

Figure 5-8. The notched circle centers at (50, 75) with radius 15 and the slot has a width

77

of 5 and length of 24.79. And the velocity is given by

(7r/3l4X50- y)
v = (z/314)(x-50)°

I:
II

(5.8)

The period of rotation is 628.
1 00

 

80

r‘rrjrr—rr

60

40

20

 

 

rrrrirr—rrirrrrj‘fi

 

OILJIIIIIIIIIIIIIIIIIIIIII

20 40 60 80 1 00

0

Figure 5-8. initial Zaiesak’s Disk (Notched Interface)

In this study, the adaptive Cartesian grid technology and second-order fast
semi-Lagrangian schemes are applied to Zalesak’s disk rotation problem. The numerical
results acquired by linear contouring, adaptive contouring and particle correction under

different levels of adaptive Cartesian grids are shown and compared.

Figure 5-9 shows the numerical results after one period under level 6, 7, 8 and 9
adaptive Cartesian grids with linear contouring. Second-order semi-Lagrangian level set

scheme is used for all simulations. From the figure, it is noticed that by refining the mesh

78

with higher level quadtree mesh the accuracy of the numerical results and mass
conservation are improved, while solution on the under-resolved area near the sharp

comers of the interface is still diffusive.

 

level 9 . .
level 8

level 7

 

 

 

 

Figure 5-9. Zaiesak’s Disk Problem: T=628, Level 6, 7, 8 and 9 Adaptive Cartesian

Grid, Linear Contouring

Figure 5-10 shows the numerical results after one period on the level 7 adaptive
Cartesian grids with level 3 adaptive contouring techniques. The accuracy and mass
conservation are dramatically improved compared to those of linear contouring strategy.
The diffusive property of level set method is greatly reduced and limited within one or

two finest cells at the under resolved area, which can be seen in Figure 5-11.

79

iOO

 

 

 

 

 

 

 

 

 

 

i H
80-
60E-
: i
40:-
i
20;
C
OIIIJJJIJgLLLIJIlIIIIlILIJ L
0 20 40 60 80 100

Figure 5-10. Zalesak’s Disk Problem: Level 7 Adaptive Cartesian Grid, Level 3
Adaptive Contouring, T=0, 628 (157)

80

 

60

40 60

Figure 5-11. Zaiesak’s Disk Problem: Diffusive Property of Level Set

Method,Results at Under-Resolved Area (near the sharp corners)

80

90

80

70

 

60
40 50 60

Figure 5-12. Zalesak’s Disk Problem: Level 7 Adaptive Cartesian Grid, Particle
Level Set Method, T=1256

1 00
90
80
70
60
50

. 40
30
20
1 0

0 1 1 1 1 I 1 1 r L

 

 

l'lIIIIIIII'IIfilIIII

IITITII'IIIIIIIIIIIIII'III

 

 

 

O
01
O

100

Figure 5-13. Zaiesak’s Disk Problem: Level 8 Adaptive Cartesian Grid, Particle
Level Set Method, T=628, 1256 (157)

81

90

80

70

60

 

Figure 5-14. Zaiesak’s Disk Problem: Level 8 Adaptive Cartesian Grid, Particle
Level Set Method, T=1256

Figure 5-12 shows numerical the results after two periods respectively on level 7
adaptive Cartesian grids. In this simulation, the particle correction algorithm is applied.
Figure 5-13 and Figure 5-14 shows the numerical results and grids after two periods on
level 8 adaptive Cartesian grids. For these simulations, the interface is discretized with
3000 linear elements and the particles (vertices of the elements) are moved with
third-order Runge-Kutta scheme with the predefined external velocity field. Based on the
numerical results, the interface is captured amazingly accurate and conservative. The
mass conservation is kept very well (less than 0.01% after two rotation periods on level 8
adaptive Cartesian grids). In order to compare between these methods, the volume loss
and CPU time with linear and adaptive contouring methods are shown in Table 5-1 while

results with particle level set method are shown in Table 5-2.

82

Table 5-1. Zalesak’s Disk: Results with Linear and Adaptive Contouring

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(T=1256)

 

 

 

 

 

 

Adaptive Cartesian Grid Volume Volume lost CPU Time/step
Exact Solution 582.20 -- --
Level 6 + Linear Contour 471.52 19.01% 0.052
Level 7 + Linear Contour 529.78 9.00% 0.176
One Level 8 + Linear Contour 550.89 5.38% 0.378
PCFiOd Level 9 + Linear Contour 569.43 2.19% 0.697
(T=628) Level 6 + Adaptive Contour 579.17 0.52% 0.773
Level 7 + Adaptive Contour 579.86 0.40% 2.434
Level 8 + Adaptive Contour 581.17 0.18% 3.520
Table 5-2. Zaiesak’s Disk : Semi-Lagrangian + Particle Correction
I Adaptive Cartesian Grid Volume Volume lost CPU Time/step
Exact Solution 582.20 -- ~-
One Level 6 581.94 0.045% 0.984
Period
Level 7 582.12 0.014% 1.545
~ (T=628)
Level 8 582.16 0.007% 2.143
TWO Level 7 582.12 0.014% --
Period
Level 8 582.16 0.007% --

 

83

 

 

5.3.2 Single Vortex

Zalesak’s disk rotation problem is used to demonstrate the ability of our method to
conserve mass. While the single vortex problem is used to test the ability of the
semi-Lagrangian level set method to accurately resolve thin filaments in stretching and
tearing flows. An example of such flows was introduced by Bell et a1 [5] and called
“single vortex” problem. This case has been studied with the particle level set methods in

[24, 25].
The velocity field is defined by the stream function
‘1’ = — 71r-sin2(m)sin2(792). (5.9)
Hence the velocity field is computed by

9:

u = = — sin2 (70c)sin(27zy)

(5.10).

v :33: sin2(7oc)sin(27zy)
fix

The initial interface and velocity field are shown in Figure 5-15. The interface is a
circle center at (0.5, 0.75) and the computational domain is l by 1. When time increases,
the interface keeps stretching into a long and slim filament-like interface. Classic level set
method will fail to capture the interface at later time because of diffusive property of level

set method. Adaptive Cartesian grids technologies are applied in this example and the

84

finest cells with highest level are generated around the interface. Higher accuracy can be
achieved than uniform mesh with comparable number of cells and it is much more efficient
than using finest cell everywhere. By coupling the particle correction algorithm, the
interface can be captured very accurately even the width of interface is less than the scale

of finest Cartesian cells.

 

o.-..o----o.-.-.noo-o ~. a
g‘¢.--——-------------—-~—

_ . . , . .. .. a.vn-arav'0—————“—na—-—----------Q.-u§‘hfi

~~~s~~---~. an

5', _.....aoaaa
.rqqrooblIllllll’
.gA-oofllllll Ill
..l.tllll"’
O8 .....;opIIIII
.

0.6

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.4

 

 

 

 

 

 

 

 

 

 

0231;:

 

 

 

 

 

 

 

 

0 51191313313113731 1J. r r”'-':1'i.:-:.'I1 1 1 r I r 111

0 0.2 0.4 0.6 0.8 1

Figure 5-15. initial Interface and Single Vortex Velocity Field

For this problem, the numerical results of linear contouring algorithm are pretty
diffusive and it cannot capture sub-grid information of the interface, interface cannot be
captured accurately after t = 1. Hence the single vortex problem is solved by adaptive

contouring and particle correction algorithms under different levels of adaptive Cartesian

grids.

85

Figure 5-16 shows the numerical results of t from 0 to 1 under level 6 adaptive
Cartesian grids with level 3 adaptive contouring. Second—order semi-Lagrangian level set

scheme is used.

 

0.8 -

0.6 -

0.4

III'IITI

 

0.2

IIII'I

 

 

0 111I1111I1111L1111I1L14

0 0.2 0.4 0.6 0.8 1

 

Figure 516. Single Vortex Problem: T = 0,1 (0.25), Level 6 Adaptive Cartesian

Grid, Level 3 Adaptive Contouring, 2nd order Semi-Lagrangian Scheme

Afier t = 1, the interface stretches longer and slirrnner very fast. And Figure 5-17
shows the numerical results at t = 3 on level 8 adaptive Cartesian mesh with level 3
adaptive Contouring. The end of interface is too slim and can not be resolved thus it is
smeared. The results start getting worse after t = 3 since the filaments are too slim to be
resolved even with the level 10 adaptive Cartesian grids, which can achieve a resolution

of 1024 by 1024 uniform mesh around the interface.

86

 

0.81-

0.61-

0.4

 

 

0.21P

 

 

 

Figure 5-17. Single Vortex Problem: T = 3, Level 8 Adaptive Cartesian Grid, dt=
628/2400, Level 3 Adaptive Contouring, Second-Order SL Scheme

The particle correction algorithm needs to be applied to simulate the motion after
=3. Initial circular interface are discretized by linear elements with comparable size, and
the resolution is at least three orders higher than the resolution of background grids. For
example, if level 8 adaptive Cartesian grids are used, then the size of each element is
comparable to 1/2”. Since the interface becomes longer, a new particle is added at the

midpoint by checking the length of each element during the simulation.

Figure 5-18 and Figure 5-19 shows the numerical results at t=3 solved on level 7 and
8 adaptive Cartesian grids by particle level set method described in section 2.5.
Compared with the interface in Figure 5-17, the interface captured here by particle level

set method is much more accurate even on level 7 adaptive Cartesian grids. For this

87

simulation, the level 8 adaptive Cartesian grid at t=3 is shown in Figure 5—20.

 

Figure 5-18. Single Vortex Problem: T = 3, Level 7 Adaptive Cartesian Grid, dt=

0.001, Particle Level Set Method

 

Figure 5-19. Single Vortex Problem: T = 3, Level 8 Adaptive Cartesian Grid, dt=
0.001, Particle Level Set Method

 

Figure 5-20. Level 8 Adaptive Cartesian Grid, T = 3
Figure 5-21 shows the numerical results at t=5 solved on level 9 adaptive Cartesian
grids and Figure 5-22 shows the adaptive Cartesian on that time. It can be seen that the
semi-Lagrangian level set solver with particle correction algorithm developed in this
study are capable of simulating slim interfaces and it is very robust and powerfirl.

1

0.8
0.6
0.4

0.2

   

0o 0.2 0.4 0.6 0.8 1

Figure 5-21. Single Vortex Problem: T = 5, Level 9 Adaptive Cartesian Grid, dt=
0.001,Particle Level Set Method

89

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. .i..i.l..l.l. iii.
0 0.2 0.4 0.8

 

 

l. I ii. I .
0.6
Figure 5-22. Level 9 Adaptive Cartesian Grid, T = 5

5.3.3 Motion under First-Order Geometric Velocity

In this section, the motions of interface under first order geometric velocity field are
simulated. The velocity 17' is dependent on the first-order derivatives of level set function
41 , i.e, I7 is equal to 17(x, (fix, ¢y,t). In the simulation, ¢x and (by are computed locally
with second-order ENO schemes, then velocities of all vertices of interface are evaluated
by numerical interpolation. After that Whitney numerical extension is employed to extend
the velocity only near the interface. Thus level set function is solved locally within i 5Ah

of the interface.

The first test case is expansion of a single circle initially centered at (0, 0) with a

radius of 1, which has been solved with finite volume level set method in Chapter 4. Here it

90

is solved with semi-Lagrangian method. The numerical results are quite accurate with a
volume of 28.22 at t = 2 (compare to exact solution 28.27 and volume loss is 0.17%). The

evolution history is shown in Figure 5-23.

 

NOD-h

llll'rlI'Illl'IIlI

 

i'\)
IIIFIII'II'I'ITII'I

 

 

4bIIIIlIIIIIIIIIIIIIIIIIIIIILIIDIIIIIIIII I

-4-3-2-101234

 

Figure 5-23. Expansion of Circular Interface: T=0, 2 (0.2), dt = 0.02, Second-Order
SL, Level 7 Adaptive Cartesian Grid

As the second case, expansion and shrinking of a triangle are simulated to
demonstrate that viscosity solutions can be achieved by semi-Lagrangian method for
non-smooth interfaces with sharp comers, where the velocity is non-smooth or even
doesn’t exist. The expansion history is shown in Figure 5-24 and the shrinking case is

shown in Figure 5-25. Based on the results, similar viscosity solution is obtained as in [46].

91

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 5-24. Expansion of ATriangie: T=O, 2.2 (0.1), (it = 0.0125, Second-Order SL,
Level 7 Adaptive Cartesian Grid

 

 

 

 

 

 

 

 

 

 

P
3 IJIIIIIIIIIIIIIIIIIIIIIIJI
-

-3“-2 -1012‘3

 

Figure 5—25. Shrinking of Triangles: T=0, 0.8 (0.08), dt = 2.667e-3, Second-Order
SL, Level 8 Adaptive Cartesian Grid

The third case is expansion of a non-convex seven-point star-shaped interface. This is

92

has been also solved by finite volume level set method in Chapter 3. The results under level

7 adaptive Cartesian mesh and semi-Lagrangian level set method are shown in Figure 5-26.

 

0.2 -

0.1

IT'IWIIIIIUI'

I

.5
N
I

 

_ T

I L

I I I I I - I I I l J I L1
-0.2 -0.1 0 0.1 0.2

b

 

 

 

Figure 5-26. Expansion of A Star-shaped Interface: T=0, 0.07 (0.007), dt =
0.0001 ,Second-Order SL, Level 8 Adaptive Cartesian Grid

5.3.4 Motion under Second-Order Geometric Velocity

In this section, the motions of interface under second order geometric velocity field
are simulated, that is, I7 is equal to 17(3, ”, it). In the simulations, all the derivatives are
computed locally with second-order ENO schemes, then velocities of all vertices of
interface are evaluated. After that Whitney numerical extension is used to extend the

velocity locally and ¢ is only updated locally.

The first test case is motion of interface composed of multiple circlers under curvature

flow, Where the speed function is equal to V = —K. Initially a small circle centers at (-2.1,

93

-2.1) with a radius of 1.0 and a large circle centers at (1.1, 1.1) with a radius of J5 .
Simulation starts from t = 0 to t = 2. For a circle under curvature flow, the speed of

shrinking is equal to V = l/r . It is easy to construct the following ODE

—Vodt=ldt=dr (5.11)
r

and integrate Equation (5.11) further to find that
r2(t)=r02 —2(t—t0). (5.12)
From the equation (5.12), the small circle will vanish at r-=0.5 second and the large

circle shrinks to a circle with radius equals to l at t = 2. The numerical results are shown in

Figure 5-27 and they match analytical solution very well.

 

 

IIII'IrrI'IIIIIIIIIIIIIl'rIrI'lIII'IIII

IILLIIIIIIIIIIIIIIIIIIIIIILIIIIIIIIIII

-3 -2 -1 0 1 2 3 4

 

 

A“:

 

Figure 5-27. Shrinking of Multiple Circles under Curvature Flow, T=0, 2
Second-Order SL. Level 8 Adaptive Cartesian Grid

94

The second case is motion of an initially non-convex seven-point star-shaped
interface under the curvature flow. According the theory of interface moving, any close
curves will shrink into a circle and finally vanish. This case is used to validate the
semi-Lagrangian solver developed in this study. The simulation starts fi'om t =0 to t = 0.002.

The numerical results are shown in Figure 28 a, b, c and d.

 

 

   

 

 

 

 

 

 

 

 

 

 

 

 

 

 

i- " b
0.2” a 0.2»_
04; 01}
0; 0;
I- ..
' C
cr- on;
I
-O.2- -0.2_-

1141—11I1111I111111111I111 "11I1rrhlLJarlrrrrlrnnrlr11

-0.2 -0.1 0 0.1 0.2 '02 -0.1 0 0.1 0.2
i c C

02" 0.2‘" d
0.1} 01:.
- C
0; 0;
)- ..
I- e
cA- an;
L: l'
-0.2- .02;

h'l‘J—lllll‘lllflll"11'4 :11I1411I1111I1111I1111IL1

-02 01 0 01 02 -02 or 0 Q1 02

Figure 5-28. Motion of Star-shaped Interface under Curvature Flow, T=0, 0.002
(0.0005), Second-Order SL, Level 7 Adaptive Cartesian Grid

95

5.4 Snmmagy

A semi-Lagrangian method for the level set equation has been realized and extended
with particle correction algorithm for 2D adaptive Cartesian grids. Both first-order CIR
scheme and second-order trapezoidal scheme are presented and tested. With the help of
fast redistancing algorithm, numerical velocity extension, fast contouring algorithm, the
semi-Lagrangian method is robust and accurate for interfacial simulation. Versatile
motions of interfaces under passive transport velocity field, first and second order
geometric velocity fields are simulated with the semi-Lagrangian method to show its
ability. The particle correction algorithm has applied to further improve the accuracy. The
results are excellent and have a very good mass-conserving property. The semi-Lagrangian

level set method can be extended to 3D adaptive Cartesian grids further based on this study.

96

 

 

Chapter 6. Multiphase Flow Computation

Numerical simulations of 2D incompressible two-phase interfacial flow by coupling
a level set method are presented in this chapter. In this research, different phases of fluids
are modeled as one fluid with variable material properties so that only one set of
governing equations are solved for the whole computational domain. Interfacial terms
such as surface tension are taken into consideration by adding the appropriate sources
terms near the boundary separating the phases. The 2D unsteady incompressible flow is
solved by characteristic finite volume methods on Zadaptive unstructured Cartesian grid,
and the material interface between different phases is simulated with semi-Lagrangian

level set method.
6.1 Introduction

In this chapter the numerical solution of 2D incompressible multiphase flow is
considered. The fluids are immiscible and steep gradients of density and viscosity exist
across the interface. For such problems, it is critical to track the interface accurately since

the material interface plays a critical role in the overall flow computation.

Since Sussman et al. [67] applied the level set approach for solving 2D
incompressible two-phase flow in 1994, there has been a great deal of interests in the
simulation of multi-phase flows involving interface capturing using the level set method.
In 1998, Sussman et al. [68] modified the numerical scheme and improved both the

accuracy and efficiency of the algorithm. These methods are accurate and can naturally

97

handle the topological changes of moving interface such as merging and breaking up, but
they are designed for structured uniform mesh only. In order to improve the accuracy, one

has to refine the mesh everywhere. This could be very costly and time-consuming.

In order to improve both accuracy and efficiency, the spatially adaptive techniques
can be employed to reduce total number of grids hence save a lot of memory and
computation time (CPU time) by adaptively refining the grids in wanted areas. As
summarized by Losasso et al. [43], there are three approaches to speed up the level set
method. The first approach is narrow band method, which simply allocates the full
amount of memory everywhere but only carries out the computation locally near the
interface as in Adalsteinsson and Sethian [2] and Peng et al. [50]. The second approach is
called adaptive mesh refinement (AMR), in which a number of uniform grids of different
resolution are used for different regions. It could be wasteful in the sense that too many
fine mesh cells are introduced in each region even when only a few cells are sufficient
there. Sussman and Almgren [69] employed such an AMR technique on solving
incompressible two-phase flows. The much more general and efficient approach is using
quadtree-based adaptive Cartesian grids in 2D [75, 76] and octree-based adaptive
Carteisan grids in 3D [28, 42]. The octree/quadtree based level set method is much more
general and efficient than AMR techniques because the total number of cells can be
minimized, and many operations related to searching and sorting can be done efficiently

using tree data structure.

In this study, a fast semi-Lagrangian level set method is used to simulate the motion
of interface accurately and efficiently. Through the fast quadtree-based redistancing

algorithm [64], the signed distance function is maintained accurately at each time step.

98

 

Solving the reinitialization equation of ¢ is avoided. Meanwhile, a second-order

characteristics upwind finite volume method [74, 81, 82] is employed to compute the
flow field, which is used as the extension velocity for evolving the interface. By
introducing the Chorin's artificial compressibility [14], the incompressible flows are
solved on 2D adaptive Cartesian grids with a dual-time stepping procedure until a
divergence-free velocity field is achieved at each time step. The flow solver developed in
this study is cell-centered, using a face-based data structure for flux calculation, and

capable of arbitrary 2D unstructured grids.

The remainder of this chapter is organized as follows. Following introduction, the
physical background and concepts are reviewed then the mathematical formulations are
given in section 6.2. Then, the detailed numerical algorithms are described and discussed
in section 6.3. After that, numerical examples are given to verify the accuracy, stability
and convergence properties in section 6.4. Finally concluding remarks are summarized in

section 6.5.
6.2 Mathematical Formulation

In this study, two-phase incompressible interfacial flow is under consideration.
Before writing the governing equation for multiphase flows, it is helpful to discuss the

single phase incompressible flow. For a single phase flow, the fluids are taken to be

1ncompress1ble such that the densrty of flurd partlcle remains constant, that rs, —E- = 0.
t

Hence the equation of mass conservation is reduced to
V - V,- = 0. (6.1)

99

 

The governing equation is Navier—Stokes equation

2%fi+v.p,r§r7, =—vr;- +2v.p,15,+p,-f',. (6.2)

‘

Here i denotes the ith phase, V,- is the velocity vector of the ith phase, P,- is pressure, D,-
is the rate of deformation tensor with the components D jk = (ajuk + aka,- I/ 2, f, is the

body force, constants p,- and y,- are density and viscosity of the ith phase fluid

respectively. The viscous effects, large density ratio and viscosity ratio across the

interface and surface tension at the interface are included.

In this study, only two-phase flows are solved for simplicity consideration. The
viscous effects, large density ratio and viscosity ratio across the interface and surface
tension at the interface are included. A key to two-phase flow computation here is that
only a single set of conservation equations are used for the whole flow field by modeling
density and viscosity as continuous functions over all computation domain. In order to
realize that, it is natural to model the density and viscosity as the function of signed

distance function as p(¢) and p(¢). The distance to the interface is used for modeling
p and p with an Heaviside function H ((15) and the sign of (15 is used to decide which

side of interface it lies such that proper material properties can be chosen. In our
simulation, the Heavside function defined in [67, 68] is employed for modeling density

and viscosity, that is,

p(¢) = pi + (pz — pi)H.(¢) (6.3)

and

100

 

#(¢) = #1 + (#2 - #1 )H e (615) (64)

where p1 and p2 , and ,ul and ,uz are density and viscosity of phase 1 and phase 2
respectively. 8 represents the uniform thickness used to model the material properties

within a band of [- e , + a] near the interface. H 5 ((15) is the Heaviside function defined in

[67, 68] and given by
'0 if ¢ < —a
H,(¢) = < 1[1 + g + lsin[”—¢)] if |¢| s e. (6.5)
2 c 7r 3
.1 if ¢ > e

 

According to our computational experience, the thickness 8 cannot be set too large
otherwise the numerical result is too diffusive. In our computation, a is set to 1.5Ax,

where Ax is size of finest Cartesian cell near the interface.

The jump conditions across the interface F between two phases (phase 1 and phase 2)

are described next. Here [] is used to represent the difference quantity across the

interface (H = []2 —[]1), t' is tangential unit vector and ii is outwarding unit normal to

the interface. Since it is viscous flow, the tangential stress is continuous hence lead to

[TD-ii] =0. (6.6)
F
The velocity at the interface is continuous thus we have

[17] r: 0. (6.7)

101

 

With consideration of surface tension, the normal stress condition satisfied
[fi-(—Pi+pD)-a] T=0'K, (6.8)

where I is unit matrix, K is curvature and 0' is surface tension coefficient.

In addition to accounting for the effects due to different material properties across
the interface, interfacial phenomena such as surface tension are modeled as appropriate
source terms of the governing equations. Since these terms are only concentrated at the
interface and usually discontinuous across the interface, they have to be smoothen and
approximated by smooth functions otherwise infinite gradients exist due to the
discontinuity. One good way to over come this difficulty is to use a smoothed delta
firnctions 6. As a matter of fact, the surface tension has been modeled by such a delta
function near the interface by Unverdi and Tryggvason [71], Brackbill and Kothe et al [8],

Sussman et al [67], and Chang et a1 [10]. The smoothed delta function corresponding to

 

Equation (6.5) is given by
'0 if ¢ < -.t:
_ 21H, _ i. gt ] .
6(¢)— 61¢ _(2£[1+c6s[ 8)] If |¢|s e. (6.9)
[0 if (9 > e

 

Based on above assumptions, the 2D Navier-Stokes equation for two-phase

interfacial flow can be transformed into a level set formulation

2+ -. -__ vp +V-(2p(¢)DI_0K6(¢)ir'+ -
a. (V VI‘ .0) p0) p(¢) f ’ ‘6“)

102

where the density and viscosity are modeled as a continuous function and (15 is a signed
distance firnction here. ii is the local unit normal of the interface, Kis the curvature,
6(¢) is a Dirac delta function defined by Equation (6.9) and a is the surface tension

coefficient.

The level set function ¢ is evolved with the external velocity field Vex, , which is

computed by incompressible flow solver. It leads to level set equation

%?+i7m-v¢=0. (6.11)

6.3 Numerical Algorithm

The governing equations (6.10) and level set equation (6.11) are solved on 2D
adaptive Cartesian mesh. The grid generation and adaptation are very fast and efficient by
employing the quadtree data structure. The detailed algorithm has been described in
Chapter 3. The level set equation is solved by the semi-Lagrangian level set method
presented in Chapter 5 since the fast semi-Lagrangian level set method is quite simple yet

very accurate, efficient and capable of quadtree based adaptive Cartesian grids.

In this section, we only focus on the algorithms specific to the incompressible flow
computation. The introduction to the numerical methods for Navier-Stokes equations is
Peyret and Taylor [51], including artificial compressibility method (ACM) by Chorin
[l4] and the pressure projection method by Chorin [15]. The projection method has been
greatly developed and applied in [38, 5, 67, 22, 54, 68, 9, 37, 30]. The artificial

compressibility method was reviewed by Puckett [55] in 1997 and further developed by

103

 

Zhao et al [81, 82, 74]. In this study, a characteristic upwind finite volume method is used

to solve the two-phase flow field on 2D adaptive Cartesian grids.

From section 6.3.1 to 6.3.5, the detailed numerical algorithm of characteristic finite
volume method is described, including the conservative formulations of expanded
governing equation, integral over a finite control volume, least square reconstruction of

solution variables, upwind characteristic scheme, and time integration.
6.3.1 Characteristic Upwind Finite Volume Method

By introducing artificial compressibility ,8 and a pseudo time I, the governing

Equations (6.10) can be expressed in a vector form with pseudo term 661 , inviscid term
I

V - EC , viscous term V - E, and source term S :

filfigwrfivfivm, (6.12)
at at

where

In Equation (6.13), the surface tension term F ST is given by

0' K 6(¢)fi

, 6.14
p(¢) ( )

FST =—

104

 

and unit normal ii and curvature K of the interface are computed by

5:,322
Will
an (6l5)
K:V.fi:_a_’.l_x. 4
8x 6y

Constant fl is 8 called artificial compressibility coefficient which affects the
convergence rate. It is noted that vector Q is a subset of W. If W is known, Q is also

known. Therefore, W is selected to be the solution variables. Obviously, if the solution

a . . .
reaches a “steady state” where 3— = 0 1n terms of the pseudo-tlme 2', then the physrcal
r

Navier-Stokes equations will be satisfied.

The above equations can also be written in a non-dimensional form based on

predefined length scale L, and velocity scale V, resulting in the non-dimensional

pVZL
0'

parameters such as Reynolds number Re = £111 , Weber number We = and

 

Froude number F, = —‘/V—_Z. In the study, dimensional equation (6.12) is solved directly
g

with characteristics upwind finite volume method in dimensional forms.
6.3.2 Finite Volume Method

The governing Equation (6.12) can be discretized on an arbitrary unstructured grid
with possibly polygon in 2D or polyhedral cells in 3D. A cell-centered scheme is adopted
here, i.e., all solution variables in vector W are stored at the cell centroid of the grid cells.

By integrating Equation (6.12) over a control volume V,, we obtain the following integral

105

 

equation

dideV+§IQdV+IVchV= [v-F,dV+ [Sam (6.16)
TV,- tV, V,- V,- V,-

Applying the divergence theorem to Equation (6.16), we obtain the following form

iIWdV+iIQdV+ [Fe-rm: §Ev-fidA+ISdV, (6.17)
d 1 dt

Vi Vi 6V,- 6V,- Vi
where ii is the unit normal of the control surfaces. The solution unknowns or degrees-of-

freedom (DOFs) are chosen to be the cell-averaged dependent variables,

[ WdV
V.

177:1 , 6.18
1 Vi ( )

 

In a second-order finite volume method, the solution variables are assumed to be

cell-wise linear. In this sense, the cell-averaged solution is identical to the solution at the
cell-centroid, i.e., W, = W, (the solution at the centroid of cell i). Given the DOFs, the

computation of the convective and viscous fluxes becomes the key of the numerical
algorithm. In our simulation, a Godunov-type finite volume method is used to compute
the fluxes [29]. There are several key components in a Godunov finite volume method,
including the solution reconstruction, flux computation, and time integration. An efficient
approach to compute the viscous flux presented in Wang [73] is employed here. The
approach is stable and efficient because it does not require a separate viscous

reconstruction. Interested readers can refer to Wang [73]. Other details of the method are

106

 

 

given next.
6.3.3 Least Square Linear Reconstruction

In a cell-centered finite-volume method, flow variables are lmown in a cell-average
sense. No indication is given as to the distribution of the solution over the control volume.
In order to evaluate the flux through a face using the Godunov approach, flow variables
are required at both sides of the face. This task is fulfilled through data reconstruction. In
this study, a least-squares reconstruction algoritlnn capable of preserving a linear function
on arbitrary grids is employed. This linear reconstruction also makes the finite volume

method second-order accurate in space.

The reconstruction problem reads: Given cell-averaged solution variables for all the
cells of the computational grid, build a linear distribution for each cell (c. g. c) using data
at the cell itself and its neighboring cells sharing a face with c. We use the fact that the
cell-averaged solutions can be taken to be the point solutions at the cell centroids without
sacrificing the second-order accuracy. Therefore, we seek to reconstruct the gradient (Wx,

W,) for cell c, which produces the following linear distribution
W(x,y)=Wc+Wx(x—xc)+Wy(y—yc), (6.19)

where (xc, yc) is the cell-centroid coordinates. The following equations can be easily

derived from a least-squares approach:

107

 

(Wf,c —Wc Xf,c —xc)
. (6.20)

M2 3M2

(Wf,c TWchf,c ”YC)

 

 

1

r
K.‘
II

where N is the number of faces in cell c, Wf, c is the solution variable at a neighboring

cell sharing face f with cell c, x f, c and y f, c are the coordinates of the cell centroid, and

matrix M is defined by

where

1 —1
M =l[ W W], (6.21)
A —1,y 1,,
N
1n=2(xf,c—xc)2
f=l
it it i
I = x , —x , —y
xyf=1fc cfc C. (6.22)
N
IW:Z(yf,c'—yc)2
f=1
A =1,,,,1yy-1,y2

Note that matrix M is symmetric and dependent only on the computational grids. If

one stores three elements of M for every cell, the reconstruction can be performed

efficiently through a loop over all faces. Figure 6-1 shows a typical stencil used to

perform linear reconstruction in an adaptive Cartesian grid.

108

 

 

 

 

 

 

 

 

 

 

 

 

 

<8. r . fix xix
g R! J \7
n
{‘x xx 2 xx
v R/ \J \J

m
0 I
<\ (x rx {*1
J 14/ \J u}
o e o c

Figure 6-1. Linear Reconstruction Stencil for Cell ion Adaptive Cartesian Grid

6.3.4 Upwind Characteristic Method

The Euler equations modified by Chorin's method are rewritten in partial differential

form in a Cartesian coordinate system for the derivation of the method of characteristics:

6a
61) +fl—J—O
6t axj
a“ (6.23)
Q-u—l—+uj%-+ui—J+2‘?—=O
6t axj 0x,- 6x,-

Here subscripts i and j equal 1 and 2, representing the two coordinates. Suppose that if is

a new coordinate normal to the surface of a control volume (outward normal direction).
In order to extend the method of characteristics to the unstructured grid solver, it is

assumed that flow in the 5 direction is approximately 1D and the above equations can

then be transformed into

109

0_p Tfl 8611;:sz

:t a; , (6.24)
u,-
67+ “j auwéxag__§x1.+au§:§éxi_0_

where 5,1. =% and tij- — g—f. In the t- 2; space as shown in Figure 6- 2, flow
1 J

variables W at time level n+1 can be calculated along characteristics k using a Taylor

series expansion and the initial value at time level n (IV):

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

w”+1 = W" + W§§,At + WtAt. (6.25)
And therefore
W W’Wk —W§ (626)
I At 6 I' .
t A
w
n+1
lit
WR WL
I L g ’
L R

Figure 6-2. The Schematic of t— 2,” Characteristics

A wave speed 21" is introduced as cf, = 21" 1in 5x1. and the normal vector

110

6e
components are defined as n . = J . By substitutin com onents of W into
x J \/— g P t
xi xi

 

Equation (6.24), we have

1 P-Pk k
—p2 +,6(u n +vn )=0
Jae, At 5 5" 5”
1 u—uk

45., 5., At

1 v—vk

 

+ u; (20 — 2" )+ u(u§nx + vény)+ pgnx = 0, (6.27)

+ v§(/to — 21k )+ v(u6nx + v§ny)+ pgny = 0

where 210 is the contra-variant velocity [to = 1mJr + vny.

In order to derive the compatibility equations, the spatial derivatives, such as u 5 ,

v5 and p; have to be eliminated from the above equations. Following the approach of

Eberle [23] for compressible flow equations, each of the above four equations is
multiplied by an arbitrary variable and all the resulting equations are summed to form a

new equation:

1
————f—_A+p B+u C+v D+TE=0, (6.28)

where
A=aip-pk)+biu—u")+civ—v"I+diT-T"I

k
B=—a/i +bnx +cny,

lll

C = anxfl + b(/l.0 — 21k + unx)+ cvnx ,

D = anyfl+ bun), + C(llo -/lk + vnx ),

E = d(2° — 2" I
and a, b, c, and d are the arbitrary variables used to multiply the equations. We define the
coefficients of the partial space derivatives to be zero, i.e., A, B, C, and D are zero such
that a linear system expressed as (DX = 0 with X = {a,b,c,d}. Variables a, b, c and d are

generally non-zero, thus the system of equations has non-trivial solution. This means that

det(<l)) = 0 and we can derive the following eigenvalues based on this condition:

a=h=fl,

23 =21=2°+(I(2°)2+,8=2°+c,
2,, =22 =20—1/(20)2+,6=20—c.

For each eigenvalue or characteristic speed, characteristic equations can be derived.

bn +cn
For example, for ,1 = 210 we have a = x—L. Substituting this into Equation (6.28),

’10

we obtain

flifluflohbb —u°)+ c(v—v°)+ d(T—T°)= 0, i.e.,

b[n,(p —p°)+ 2°(u —u0)]+c[ny(p—p0) + 2°(v—v°)]+ d2°(r— T0) = 0.

For any b, c, and d the above equation is always satisfied. Therefore all the terms in the

112

 

 

square brackets are zero. As a result, we have where

(u-u0)ny —(v—v0)nx =0, (6.29)
For 21:21
1_ l l l
p—p ——}I[(u—u )nx+(v—v )ny]. (6.30)
For 2:212
p-pz =—22[(u—u2)n,+(v—v2)ny]. (6.31)

Finally it, v, and p are determined using the above characteristics equations (6.29), (6.30)

and (6.31):
u = fnx +u0n}2, -v0nxny = 0,
v = fny +v0n3 —u0nxny = 0,
p=p1-21[nx(u—u1)+ny(v—vlb,
or
p=p2 —22[nx(u-u2)+ny(v—v21],
where

c = (10Y+,6, f =2i[(p—p°)+n,(;l‘u1-22u2)+ny(21v1—22v2)].
c
Flow quantities at (n+1) time level obtained from the above equations on the

113

characteristic lines are then used to calculate convection fluxes at the control volume
interface, whilst those on different characteristics at (n) time level are approximately

evaluated by an upwind scheme using the signs of the characteristics:
Wj = %[(l + signaj ))WL + (l — sign(/1j))WR] , (6.32)

where WL and WR are obtained using the reconstruction presented earlier.
6.3.5 Time Integration

A fully implicit second-order backward-difference scheme is used to discretize the

physical time (denoted by t), i.e.,

n+1 _ _n p—l
am + 3Q! 4Q! + l =Ri(Wn+1), (633)
at 2At

 

where Ri(W) denotes the steady residual. Let RAW) denote the unsteady residual given

by

Rpm): R.(W..i)_ 30;”1-40." +2?"
I l 2At

9

then Equation (6.12) becomes

 

n+1
an; = 1‘1,(W"+1). (6.34)
T

The equation (6.34) is then driven at each time step to a steady-state, in which

114

B,(W"+l)= O. The time-integration scheme in the pseudo-time is a multi-stage Runge-

Kutta scheme with local time-stepping.
6.4 Numerical Results

In this section, the coupling of the semi-Lagrangian level set approach and the finite
volume characteristics upwind finite volume method described in the previous sections is
performed and validated using several two-dimensional classical multi-phase flow
problems. First the classic broken darn problem is solved in order to validate the present
algorithm. The second case is relaxation of the interface due to the surface tension.
Finally a 2D rising air bubble in a fully filled container with bubble shedding is simulated.
We carried out all the computations on both uniform Cartesian and adaptive Cartesian

meshes to compare the accuracy and efficiency.
6.4.1 Broken Dam Problem

In order to test the accuracy and robustness of the method, the two-dimensional
broken darn problem is selected as a benchmark test. The broken dam case is a classic
problem in two-phase flow with a free interface. The name “broken darn” basically
means the sudden collapse of a rectangular column of fluid onto a horizontal surface and
it is used in modeling the sudden failure of a dam. It is popular because of the relatively

simple geometry and initial condition. Figure 6-3 illustrates the geometry, initial interface

and an initial adaptive Cartesian grid employed for this case.

115

 

1 .5 Air
pm, =1.205kg/m3

   
     

0-5 b=0.9
pwater =1000kgh"

lllillllllllll

0.5 ' 1 1.5 :2

 

 

Figure 6-3. (a) Illustration of 2D Broken Dam Problem at T=0; (b) Initial Level 5 to
7 Adaptive Cartesian Grids and Zero-Level Set

The initial conditions for the broken darn problem are prescribed as follows. The
computational domain is a square with an edge length of 2m. The height and width of the

water column are 0.9m and 0.45m respectively. The density for each phase is shown in
Figure 6—3. Viscosity ,uwater and pair are 1.01 x10—3 Pa -s and 1.81 x10'5Pa-s

respectively. The initial velocity is s zero everywhere and the pressure distribution is
defined to be hydrostatic pressure relative to the top surface of the computational domain.

Thus for any point within the water column (x S 0.45m and y S 0.9m ),
P(xt.V) = (pairhair + pwaterhwater )g

and for any point outside the water column, pressure is given as

POW) = pairhairg -

Symmetric velocity boundary conditions are applied to the four walls. At time t =0, the

dam is removed instantaneously and the water column collapses due to gravity. For

116

comparison purposes, three meshes were tested in for this simulation, a level 6 (64 by64)
uniform Cartesian grids, level 5 to 7 (mesh size between 1/25 and 1/27 of the domain size)
and level 5 to 8 adaptive Cartesian grids. Time refinement study for each grid has been
performed to guarantee that the solution is time step independent. The profiles of the free
surface between air and water at time t = 0.24 and t = 0.42 on different meshes are shown
in Figure 6-4 and Figure 6-6 respectively. The triangulated level 5 to 7 adaptive Cartesian
grid at these times are shown in Figure 6-5 and Figure 6-7. There is no visible difference
between the free surfaces computed on level 5 to 7 and level 5 to 8 adaptive Cartesian
grids. The free surface generated on level 6 uniform mesh is less accurate, and lags that

computed on the two adaptive meshes.

 

Level 6 uniform
_ _ ._ .. 5to7Adaptive

__ 5 to 8 adaptive

1.6

1.2

0.8

0.4

rjI'lI’Ir'IfrT'TIII'lI—II

   

 

 

OIIIIIIIIII IIIIIIIIIIJLI

0 0.4 0.8 L12 1.6 2

Figure 64. Free Surface Profile at T=0.24 (Red solid line is on level 6 uniform

mesh, blue dashed line is on level 5 to 7 adaptive mesh, and green dotted line is

 

on level 5 to 8 adaptive mesh)

117

 

   

 

 

 

 

4
2
0
=
T
S
d -
2w .»..
.4»..». .»..4»..4 ...4....4..»..4..4..»..»..4..»...»..u_.4. .4..4...»4..4. 2 -
.4»..4»..4»..»t. .u..4»..n ..4»..u..4»..»4..4»...4..4»..»..u..u..».. ..m..»6..m..»u. n m e e -
.n..u..u..u. .»..n..6..»4. .4 .u..n.4..4..4..u..4..4.4..»..4..4»....» .B .W .
n.u..»u..»u. ...4...4..M4..»4. ..6..»6..4». ..4»..4 ..4»..4..4»..4..4 ..4 ..4..4»..u..»..M. S O I -
.u..n..»4 ..»4..4.W ..»..»4.. .»4....»..4..»4. ..»4..»..».......4..»..»..n..4...4..n.4..u..n..4». e an D. D.
4....»....4....4».n.»4..»..»..»..»..»....u..u....»..»4..4..u..4..4.n.»4..»..4»..n..u..n. .n n a a
»..u.4..u....»..u..u..n u..».u..».1..».n..u.......4....»....4..4..4». a u d M
C 6 A 8
I
.n.u.. ...6..m. e e 7 0
“New. .»..»4.» .N V 0 t
.4»..H..4»..n..u. .n..»..v..u ..m. 9 I
.4»..4»..4»..»4..u.4 .u..n..n..4». .n..»»..4»..n..6»..»4.. 6..»...4...»..M. a L 5 5
.4»..4»..4»..4»..N» .N..M..4»..H. ...4....N..M..4..4..4..4»..4..»..4..4.. d
.4.»..M..N..N..u. .4»..»4..u..u A
human new. _
flfififififififirdfiflfifififififififlfl4 7 _
afifififimmfifl. ..6..4»..4..»4..4..4..»4...»4..n..u. O
.N..M..N..N..M..M..N..N.. .4»..4 ..4..M..4..u..»..u..n..n..»4..4»..m..m..44..m .I .
.4»....»..»4..u..»4..n..»4..4»..» ...»..4»..»4..4..»..4..n..4»..4»..u..»4.. »..»..n..4..6...n...» 5
.u..»4..u..»4..»..n..4»..4..4».» 4. ......4....4..4...n..n.u.»..4.... _
.H..». ..4..4 ..4..4. 4..4.....4..4..4..4..b..4....n .N. .N..4..4..4.. .424.“ a
.H..b4..m..6>..m..4n..4.. ...4....4....4..4..>..4.4..4.N..4>. . V
fififlfixfifimfiflfinfixfifimmwfl e
4..N..N..N..4..4..H..4 ..W..4..4>.fl..4>..4..4»..4>...4...»4.u.4. . L
.1».. ....n........».4..4....u..n..4»..u..n..4»..u..u.
.»4..4..»4..4.. ..4..4..u..n...4...4»..4»..n...4..u..n..u..u.u. » .w
.m..4..u..u..n..4. ».u. t
.»..»4. .4.. N..»4. a
.».»n..n..n. new. w
.4»..4»."m.» MM...» 9 _ p .
PALV’AVAVAVZVDO n __-—___—__-—___
2 .5. 1 5. o .m 2 5 2 8.
1 0 T 1 1 0
. :
6
6
m
u
.m.
F.

me is

is on level 6 uniform

me

we mesh, and green dotted i

118

0.42 (Red solid I

level 5 to 7 adapt
on level 5 to 8 adaptive mesh)

Ine IS on

Figure 66 Free Surface Profile at T

mesh, blue dashed l

.N..)4..N..H..N..N.N..H¢H..4¢4. .N..N..N.&.N..M.v4 .4>.. >24KAVM..4..N..H..N.M..N¢N..H 2K4..4).
Wfivuhvu. ..MMEW .Emmfifififi.“ .KAEN...M..H¢H¢H..4 ..MNWAEW.
‘
D
‘

V6>'
b

..4. ..N. .N. .N..4v..N..N..H..4..4.vH..N.
.»4..M..4..4 ..4....M..4».
..4>..»4..4..»4..N..4..

.»4..4 2.»..64..4..

»»» .u.
...»...K...4..b..r..>..4. ..6..6..»4..v4..>6..n. ...
4..4..»1..4 6..N..»4..»4..4..4..»4..4.4 .63». .M..N..»4..M..>4..b...
..H..4..4..4..4»..u..4 »..4». ..4..4..4»..4..»..H..N..
".».N. .4..»4. ..4..4..N. .»4..4..4..»..4..4 ..4..4..4». .N.
.»N..N..4..4 .>.>..».2.4.2.4.22 .4
..N..H..N..4...4...»4..4..4..»4..»A ..4..»4..4..m.
..H..N..4»..4..4>..4..4..4..N..M..m..»4.. 2
..u..n..4..4»..n...4...u..»..4»..n..».4
>..»..N..»4..»4..H..41..N..M. 4..4»..H.
"422.452.". n..»4..n..u. . .
Wen. M. 4..4»..M..4. .
2
.H.
.H..H..M.N..H..4»..N. .»4.
fifififififlfiflfififififiwmm .6..4.»4..4. "42.4.3232 ..4..4.
.N..4..4..4..4..4.2.4..4...4..2.4..1..H.2.h..b..>..>.hue». N..4..4..4...6m
..4..M..64..6..14..6..4..4..4.. ..r..>.2.>..4..>..>..>.

024

'.
202024.
>203.

50.0.1

m “

4:2:
224
K4'

#2020202

v

>'4K4K4::
0:020:020V51K02“

K4

4
£020.03

'24 4:
r
26262624

A
A

"W
’vv
PA

A

W.
.N.
N.

m..m..46..»4 .4..M..n..4..4..u....u..p..» ..6..6»..4». .6 wei..6..m...
».u..n. .4..M.4..4..»4..n..4..m..».»4..6..m..»..4..u..4..
».Vm..4»..»44..4..u»..n4..4. 4. ..4.. .4..»4.4. ..4. 4..4..4..4

‘ A‘ «A 0

..».. ...

...»..»..4.. ..»4.. .»4..»4...
..4..4..6. ..4..H

26262:
»'6z4»:6»
K4K4K4

1.6
K4
020.4

‘P AP
">66 6..»H.. 6.W..W..6»..H..N..»4..4»..4»..
N..4..N .4..4..4. .N..4..4..4..4..4..».
.4..4..»4.. 6..6..6..4»..N..4..4»..N..4 ..4).
..m..m 4..4»..64..4..4 ..4>..4»..4>.
.4>..4..N..46..4 .N.. 4..H..M..M...N»..4 ..M. ..4..4..>.2.b..>..m..4..>
»

A

4

 

Figure 6-7. Triangulated Level 5 to 7 Adaptive Cartesian Grids, T=0.42

1 height

1mens1ona

a1 length IJa and non-d'

1mension

The evolution history of the non-d

t/(a/ g)“2 has been plotted

. . '6
ltrmet

1mensrona

H/b of the water column versus non—d

rical solution and experimental surge front

between nume

aI‘ISOIl

6-8. The comp

igure

1nF

1 solution was computed on

this figure. The numerica

1n

and column height can be found

1

It shows an agreement between the numerica

ran mesh

the level 5-7 adaptive Cartes

solution and experimental results in [45].

119

 

— Numerical result

3.5 - .
0 Experiment

 

 

 

 

 

 

 
 
   

 

 

 

 

 

 

L/a
2
1
— Numerical result
0.9 — 0 Experiment
0.8 ~
0.7 —
H/b 0.6 _
0.5 —
0.4 l 1 L l
0 0.5 l 1.5 2
t' = t/(a/g)"2
(b)

Figure 6-8. Comparisons of Numerical and Experimental Results:
(3) Non-Dimensional Length of Surge Front against Time
(b) Non-Dimensional Column Height against Time

120

6.4.2 Relaxation Bubble

The second test case that we tested is bubble relaxation, which has been studied by
Chen et al. [12]. For this case, the motion of non-circular bubble in viscous flow under
nonzero surface tension and zero gravity conditions is simulated. According to [12], any
fluid particle whose shape deviates from circular (spherical in 3D) should recover the

circular (spherical in 3D) static and stable shape.

First, the relaxation of an elliptical bubble was simulated on uniform mesh to

compare. The initial ellipse axes are 5/ 16m and 1/5m.The fluid properties outside and

inside of the bubble are ,01=p2=1.0kg/m3 , p1=p2=0.1Pa-s , g = O and
0’ =1.0 N/m . Under such conditions, the elliptical bubble will restore to a circular

interface within 0.3 second. The simulation has been carried out on level 4, 5 and 6
uniform Cartesian grids and level 5 to 7 adaptive Cartesian grids for comparison purpose.
Figure 6-9 shows the evolution history of the interface and Figure 6-10 shows the
velocity field at t =0.3 second on level 5 uniform Cartesian grids. Figure 6-11 shows the
interface at t =O.3 second under different grids and it is noted that the solutions are
converged with higher accuracy through refining the mesh. Table 6-1 shows the area the

interface enclosed at t = 0.3 second under different grids.

12]

Table 6-1. Relaxation of Elliptical Bubble in Viscous Flow: Volume, T=O.3 second

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Grids Volume(m2) Volume Loss
(%)
Exact Solution 0.1963 --
Level 4 Uniform 0.1819 7.30
=0.3 Level 5 Uniform 0.1905 2.95
Second
Level 6 Uniform 0.1942 1.07
5 to 7 Adaptive 0.1950 0.66
0.4 :-
02 E-
L
0 E
L
-0.2 _-
-o.4 :-
plllllLllLJJlLlllllLllllll

 

 

 

-0.4 -0.2 0 0.2 0.4

Figure 6-9. Relaxation of Elliptical Bubble in Viscous Flow: T =0, 0.3 (At = 0.01),

Level 5 32 by 32 Uniform Cartesian Grids

122

 

‘~“. ‘I
\ ‘ ‘ I I
‘\‘\':'I
\ I I

 
  

     

’ §~-‘-.‘.-
04 / ..--
o \\\\ "

\II/

.V.I.1 I . .1.
L---. Cl 4 4

ova-on.-
‘ "’.a.—o—.—
\\ ‘ I’ M“

0.2

O
J-U- J'lJ-I-L

5':
N
J

.I-

IIIIII
1 Q

 

 

.
-

.
"r
.l
O

I
H

I
I
I-I
I
bl
I

4
-

0
4‘
p4
-
. .
_.
'l

l
-

4
' l
....
I
hi
I
'l

I
n

I
I
"I
0
.I
|

0
-

|
I'—.‘
,_.
r...
h-

 

i-Oéi ' -o.2 o 0.2 0.4
Figure 6-10. Relaxation of Elliptical Bubble in Viscous Flow: Velocity Field at T

=0.3, Level 5 Uniform Cartesian Grids

 

Level 4 UrTIform

- - - ° Level 5 Uniform
- ------ Level 6 Uniform

_ ......... 5 to 7 Adaptive

0.4

0.2

 

-0 .4

O
lIIIIIIIIITTIIIIYIIYIIIII

 

 

lgllljllllllllllllJLJLllll

-0.4 -0.2 0 0.2 0.4

 

Figure 6-11. Relaxation of Elliptical Bubble in Viscous Flow: Interface at T =O.3
on Level 4, 5, and 6 Uniform Cartesian Grids and level 5 to 7 Adaptive Cartesian
Grids

123

Next, the relaxation of a five-point star-shaped bubble was simulated. The equation
of the interface is given by
r = a(1+ a cos(n 0))
where a is mean radius, a is coefficient of perturbation and n is the mode. In this

example, a= 0.5m, a=0.4 and n =5.

The fluid properties outside and inside of the bubble are same as relaxation of
elliptical bubble. And the simulation is carried out on level 5 to 8 adaptive Cartesian
mesh. The initial interface at t = 0 on adaptive Cartesian grids are shown in Figure 6-12.
Figure 6-13 shows the evolution history of the interface from 0 to 0.5 second. The
interface approaches steady state (a circular interface) afier t=0.4 second with small-
amplitude oscillation around it. Figure 6-14 shows the velocity field at t =0.1 second. In
order to check the property of mass conservation for incompressible flow, the history of
volume preserving versus the time step is plotted in Figure 6-15.The final volume loss at
t = 0.5 second is very small (only 0.52%).

1

-1-

 

Figure 6-12. Relaxation of Stare-Shaped Bubble in Viscous Flow: Initial Interface,
Level 5 to 8 Adaptive Cartesian Grids

124

0.5

-1 -O.5

 

ljjllilllll

 

 

r r I r l l I I

IIIIlIIIIlIIIIlIIIr

 

 

Figure 6-13. Relaxation of Stare-shaped Bubble in Viscous Flow: Evolution

History, T =0, 0.

 

5(0.05), Level 5 to 8 Adaptive Cartesian Grids

 

Figure 6-14. Relaxation of Stare-shaped Bubble in Viscous Flow: Velocity Field, T

:01, Level 5 to 8 Adaptive Cartesian Grids

125

..L
tr]

r

 

0.8

IIIU'IIII'ITTIIIFIT'

 

oIIIIlrrrllrrlllrLLrlIlrrl

O 100 2 300 400 500
Wme Step

Figure 6-15. Relaxation of Elliptical Bubble in Viscous Flow: Volume Preserving

History versus Time Step
6.4.3 Rising Air Bubble in Fully Filled Container with Bubble Shedding

To further demonstrate the ability of the present method in solving multi-phase
problems, the evolution of a two-dimensional rising air bubble with surface tension is

simulated. For this problem, the Reynolds and Weber number are computed by

/ 2 2
Re : pwater\/E(2R0)3 and Wb = pwaterg(2R0)
llwater 0'

We choose the case with typical parameters Reynolds number Re=100, Weber
number W, = 200, Froude number F, = 1, density ratio 1000/1 and viscosity ratio 100/ 1.
The initial set up is illustrated in Figure 6-16. The computational domain is a square with

edge length = 0.2m and the initial radius of the bubble R0 is 0.025m. The center of the

126

initial bubble is located at xo=0 and yo=-0.065. The density of water and air bubble are

1000 kg/ m3 and 1kg/ m3 respectively. The viscosity coefficients for water and air are

0.1118 Pa - s and 0.001118 Pa ~s . The static hydrodynamic pressure is chosen as the initial

condition as defined in Zhao [82]. Viscous wall boundary condition is employed for this

 

 

 

 

 

case.
0.1 _
0.05_- pm... =1000kg/m
-0.05 -
I 1kg/m3
at. ‘WETT - a‘ ‘ ' ”obs"

Figure 6-16 (a) Initial Setup of 20 Rising Bubble Problem; (b) Initial Level 5 to 7
Adaptive Cartesian Grids and Zero-Level Set

In order to compare the accuracy and efficiency, both uniform and adaptive
Cartesian meshes with different levels are tested. The computation time starts from 0 to
1.0 second (equivalent non-dimensional time 0 to 4.47 seconds). Table 6-2 shows the
different cases on different meshes. Figure 6-17 to 6-21 shows the numerical results in
the order as Table 6-2 defined. For each case, time-step refinement was carried out to

guarantee that the solution is time step independent.

127

Table 6-2. List of Test Cases for Bubble Rising Problem

 

Case Type of Figure Level Number of total Initial number
number Grid number time steps of cells
1 Uniform 6-17 6 1000 and 2000 4096
2 Uniform 6-18 7 1500 and 3000 16384
3 Adaptive 6-19 5 to 7 1000 and 2000 2125
4 Adaptive 6-20 5 to 8 1500 and3000 3511
5 Adaptive 6-22 6 to 9 1500 and 3000 9262

 

 

 

The bubble begins to deform due to the circulation of velocity field through the
center of the bubble and a mushroom-shaped bubble is formed with time increase.
Meanwhile, the bubble rises due to the buoyancy. The diameter of the deformed bubble
increase as it rises up. Also both inner and outer radii of the mushroom-shaped bubble
increase. At later time, the bubble breaks up and two small bubbles detach from the main
bubble. Two detached small bubbles are captured accurately in the numerical tests by
using higher level adaptive meshes. Although all the cases show a similar evolution
patterns, the accuracy and efficiency are different. The numerical results for these cases

are discussed for comparison.

Figure 6-17 shows a quite low resolution for the interface on coarse level 6 uniform
Cartesian grids. Because grid resolution is too low, the details of pinching are under
resolved and the numerical solutions are very diffusive. By refining the mesh uniformly,
we carried out the simulation on a level 7 (128 by 128) uniform Cartesian mesh. The
resolution and accuracy are improved greatly. The formation of mushroom shaped air

bubble and process of pinching can be seen in Figure 6-18 more clearly. But the

128

computational efforts increase a lot because the number of cells increases by a factor of 4

and time step reduce to half of that in first case.

 

GD 09

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 6-17. The Evolution of Rising Bubble on Level 6 Uniform Cartesian Grids,
T=0.2, 1.0 (At =0.1). (l' he figures are ordered from top to bottom and left to right)

To demonstrate the advantage of adaptive Cartesian mesh, the motion of rising
bubble is simulated on a level 5 to 7 adaptive mesh with 2511 cells initially, which has
less munber of cells than level 6 uniform mesh but the grid resolution is one level higher.
The sequences of evolution are shown in Figure 6-19. The equivalent resolution and

accuracy are achieved as in case 2 using a level 7 uniform mesh, but the number of cells

129

is only around 1/8 of level 7 uniform meshes. This demonstrates that using the adaptive
Cartesian grid significantly reduces CPU time while achieving high resolution and

accuracy.

 

0b

 

 

 

 

 

 

.. .11..

 

 

a A C

Figure 6-18. The Evolution of Rising Bubble on Level 7 Uniform Cartesian Grids,

 

 

 

 

 

T=0.2, 1.0 (At =o.1)

In order to obtain even better numerical results, level 5 to 8 adaptive Cartesian mesh
is employed. Initial total number of cells is only 3511 but the effective resolution is same
as level 8 256 by 256 uniform meshes (total 65536 cells). Figure 6-20 shows much more

accurate numerical results. On level 5 to 8 adaptive grids, the stretching of the bubble

130

skirt and process of breaking up can be observed much more clearly in Figure 6-21. It
appears that the numerical diffusion is greatly reduced meanwhile the movement of
detached small bubbles such as rotation and rising are captured accurately. Figure 6-23
shows the best results simulated on level 6 to 9 adaptive Cartesian grids. They are even

more accurate than those on level 5 to8 adaptive Cartesian grids. The volume loss of the

air bubble at final time is only 2.2%.

 

 

 

I”

 

 

6% (3

 

 

65 A 55

Figure 6-19. The Evolution of Rising Bubble on Level 5 to 7 Adaptive Cartesian
Grids, T=0.2, 1.0 (At =0.1)

 

 

 

 

 

 

 

l3l

 

 

 

 

 

 

 

 

 

u ..

 

 

O 0

00 00

 

 

 

 

 

 

 

Figure 6-20. The Evolution of Rising bubble on Level 5 to 8 Adaptive Cartesian
Grids, T=0.2, 1.0 (At =0.1)

 

 

A R Q

o o "O 0' O O

 

 

 

 

 

 

Figure 6-21. The Process of Pinching of Shedding Bubble on Level 5 to 8
Adaptive Cartesian Grids, T=0.74, 0.76 and 0.78

132

The adaptive Cartesian mesh, triangulated adaptive Cartesian mesh and the velocity

field at a certain time t =0.66 are shown in Figure 6-22 (a) (b) and (c) respectively.

0.025

-0.025

-0.05

 

(a)

0.025

-0.025

-0.05

 

(b)

133

0.025

\ \t 11111 ... I Million, I

I;
1;,

 

0.04

-0.02

0.02

 

-o.o4‘

(C)

Figure 6-22. Level 5 to 8 Adaptive Cartesian Grids, t

0.66:

0.66

(a) Adaptive Cartesian Grids at Time t

(b) Triangulated Adaptive Cartesian Grid

(c) Velocity Field at This Time

134

 

 

 

 

 

 

 

 

 

 

@ @' 5".

 

 

 

 

0 O

O O 0 0

 

 

 

 

 

 

 

Figure 6-23. The Evolution of Rising Bubble on Level 6 to 9 Adaptive Cartesian
Grids, T=0.2, 1.0 (At =0.1)

6.5 Summarv
A semi-Lagrangian level set method is coupled with a characteristics upwind finite
volume incompressible flow solver to tackle multi-phase flow problems. Given an

accurate flow field, the semi-Lagrangian level set solver can evolve the interface

accurately and efficiently. The characteristics upwind finite volume incompressible flow

135

solver is stable and second order accurate. Through the adaptive refinement, the
computational efforts are greatly reduced while achieving higher solution accuracy for
two dimensional multi-phase problems. Several classic validation cases have been used to

test the present method, and results are satisfactory compare to the work by Sussman et al.

[67].

136

Chapter 7. Conclusions and Future Work

7.1 Conclusions

To apply the spatially adaptive techniques to level set method and multi-phase computation,
a fast and general interface-based grid generator for 2D adaptive Cartesian mesh has been
developed in this study. It is efficient and successful for smooth/non-smooth, convex/non-
convex, close/non-close, single/multi-component interfaces. Through interface-based adaptation,

the grids can be automatically refined near the interface through recursive subdivision.

A new finite volume (FV) level set method for 2D arbitrary unstructured grids is developed
in this research. The finite volume level set method is successful implemented on adaptive
Cartesian mesh. Two types of data reconstruction approaches were carried out and designed
accuracy is achieved through a variety of numerical experiments. To further compare the
accuracy and efficiency, a fast and general second-order semi-Lagrangian level set solver is
developed. A particle correction algorithm has been extended to semi-Lagrangian method to
further improve the accuracy and mass or volume conservation. Versatile motions of interfaces
under different velocity fields have been successfully simulated to demonstrate the ability of the

solver we developed.

Finally the level set method is applied to 2D incompressible two-phase interfacial flow. The
level set solver is successfully coupled with a finite volume incompressible flow solver to solve
2D multi-phase problems. The characteristics upwind finite volume incompressible flow solver

is stable and second order accurate. Several classical two-phase problems are simulated

137

successfully. The motions of interfaces are captured accurately. With the help of adaptive

techniques, the computational efforts are greatly reduced while achieving high solution accuracy.

7.2 Future Work

In this research, finite volume level set solver, semi-Lagrangian level set solver coupled
with particle correction algorithm, and a 2D upwind characteristic finite volume solver for
incompressible flow have been developed and used to simulate the motion of interface in 2D
two-phase flow. In the future, current work can be extended to 3D easily based on this study.
Coupled with an accurate 3D incompressible flow solver, the level set method will be able to

tackle more complicated 3D multiphase problems.

138

APPENDICES

139

N h K 38

@IIVIQt‘QN

I

('5

I

<

a?“
~a

ass

«g $1

APPENDIX A

NOMENCLATURE

Level set function or signed distance function

Unit normal of the interface
Curvature of interface

Artificial compressibility coefficient

Physical time
Pseudo time

Density
Viscosity coefficient
Surface tension coefficient

Velocity

Rate of deformation tensor

Inviscid flux vector
Viscous flux vector

Surface tension
Vector of conserved variables

Extended vector of solution variable

Reynolds number
Froude number

Weber number

140

APPENDIX B

PUBLICATIONS ASSOCIATED WITH THIS WORK

Wang, Zhu and Wang, Z.J., “Multiphase Flow Computation with Semi-Lagrangian Level
Set Method on Adaptive Cartesian Grids”, to be submitted to Journal of Computational
Physics.

Wang, Zhu and Wang, Z.J., “Finite Volume Level Set Method on Adaptive Cartesian
Grids for Interface Capturing”, to be submitted to International Journal of Numerical
Method in Engineering.

Wang, Zhu and Wang, Z.J., “A Fast Level Set Method with Particle Correction on
Adaptive Cartesian Grids”, AIAA paper, AIAA-2006-0887, 2006.

Wang, Zhu and Wang, Z.J., “Multi-Phase Flow Computation with Semi-Lagrangian
Level Set Method on Adaptive Cartesian Grids”, AIAA paper, AIAA-2005-1390, 2005.

Wang, Zhu and Wang, Z.J., “The Level Set Method on Adaptive Cartesian Grid for
Interface Capturing”, AIAA paper, AIAA-2004-0082, 2004.

141

[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

BIBLIOGRAPHY

Adalsteinsson, D. and Sethian J ., “A Fast Level Set Method for Propagating
Interfaces”, J. Comput. Phys, Vol. 118, pp.269-277, 1995.

Adalsteinsson, D. arnd Sethian, J ., “The Fast Construction of Extension Velocities
in Level Set Method”, J. Comput. Phys, Vol. 148, pp. 2-22, 1999.

Agresar, G., Linderman, J ., Tryggvason, G., and Powell, K., “An Adaptive
F ront-Tracking Method for the Motion, Deformation and Adhesion of Circulating
Cells”, J. Comput. Phys, Vol. 143, pp.346-380, 1998.

Barth, T]. and Sethian, J .A., “Numerical Schemes for the Hamilton-Jacobi and
Level Set Equations on Triangulated Domains”, J. Comput. Phys, Vol. 145,
pp.l-40, 1998

Bell, J ., Colella, P., and Glaz, H., “A Second-Order Projection Method for the
Incompressible Navier-Stokes Equations”, J. Comput. Phys, Vol. 85, pp. 257-283,
1989.

Berger, M. and Oliger, J ., “Adaptive Mesh Refinement for Hyperbolic Partial
Differential Equations”, J. Comput. Phys, Vol. 53, pp. 484-512, 1984.

Berger, M. and Colella, P., “Local Adaptive Mesh Refinement for Shock
Hydrodyanrnics”, J. Comput. Phys, Vol. 82, pp. 64-84, 1989.

Brackbill, J .U., Kothe, DB, and Zemach, C., “A continuum Method for
Modeling Surface Tension”, J. Comput. Phys, Vol. 100, pp.335-353, 1992

Brown, D., Cortez, R., and Minion, M., “Accurate Projection Methods for the
Incompressible Navier-Stokes Equations”, J. Comput. Phys, Vol. 168,
pp.464-499, 2001. ‘

Chang, Y.C., Hou, T.Y., Merriman, B., and Osher, S., “A Level Set Formulation of
Eulerian Interface Capturing Methods for Incompressible Fluid Flows”, J.
Comput. Phys, Vol. 124, pp. 449-464, 1996.

Chen, S., Merriman, B., Osher, S., and Smereka, P., “A Simple Level Set Method
For Solving Stefan Problems”, J. Comput. Phys, Vol. 135, pp.8-29, 1997.

Chen, T., Minev, RT, and Nandakumar, K., “A Projection Scheme for
Incompressible Multiphase Flow Using Adaptive Eulerian Grid”, Int. J. Numer.
Metho. Fluids, Vol. 49, pp.l-l9, 2004.

Chopp, D., “Computing Minimal Surfaces via Level Set Curvature Flow”, J.
Comput. Phys, Vol. 106, pp 77-91, 1993.

142

[14] Chorin, A., “A Numberical Method for Solving Incompressible Viscous Flow
Problems”, J. Comput. Phys, Vol. 2, pp. 12-26, 1967.

[15] Chorin, A., “Numberical Solution of Navier-Stokes Equations”, Math. Comp., Vol.
22, pp. 745-762, 1968.

[16] Coirier, W. J. and Powell, K. G., “A Cartesian, Cell-Based Approach for
Adaptively—Refined Solutions of the Euler and Navier-Stokes Equations”, AIAA
paper. 0566, 1995.

[17] Coirier, W. J. and Powell, K. G., “Solution-Adaptive Cartesian Cell Approach for
Viscous and Inviscid Flows”, AIAA Journal, Vol. 34, pp. 938-945, 1996.

[18] Courant, R., Isaacson, E., and Rees, M., “On the Solution of Nonlinear
Hyperbolic Differential Equations by Finite Differences”, Comm. Pure Appl.
Math., Vol.5, pp.243-249, 1952.

[19] Crandall, MG. and Lions, P-L., “Viscosity Solutions of Hamilton-Jacobi
Equations”, Tran. AMS, Vol. 277, pp. 1-43, 1983.

[20] Crandall, MG and Lions, P-L., “Two Approximations of Solutions of
Hamilton-Jacobi Equations”, Math. Comp., Vol.43, pp. 1-19, 1984.

[21] Delanaye, M. and Liu, Y., “Quadratic Reconstruction Finite Volume Schemes on
3D Arbitrary Unstructured Polyhedral Grids”, 14 AIAA CFD Conference paper,
AIAApaper. 99-3259, 1999.

[22] E, W. and Liu, J .-G, “Finite Difference Schemes for Incompressible Flows in the
Velocity-Irnpulse Density Formulation, J. Comput. Phys, Vol. 130, pp. 67-76,
1997.

[23] Eberle, A., “3D Euler Calculations Using Characteristic Flux Extrapolation”,
AIAAPaper 1985-0119, 1985.

[24] Enright, D., Fedkiw, R., Ferziger, J ., and Mitchell, 1., “A Hybrid Particle Level
Set Method for Improved Interface Capturing”, J. Comput. Phys, Vol. 183,
pp.83-116, 2002.

[25] Enright, D., Losasso, F., and Fedkiw, R., “A Fast and Accurate Semi- Lagrangian
Particle Level Set Method”, Computers and Structures, Vol. 83, 479-490, 2005.

[26] Fedkiw, R., Aslarn, T., Merriman, B., and Osher, S., “A Non-Oscillatory Eulerian
Approach to Interface in Multimaterial Flow (The Ghost Fluid Method)”, J.
Comput. Phys. Vol. 152(2), pp.457-492, 1999.

[27] Frolkovic, P. and Mikula, K., “Flux Based Level Set Method: A Finite Volume
Method for Evolving interfaces”, Preprint IWR/SFB No. 2003-15,
Interdisciplinary Center for Scientific Computing, University of Heidelberg, 2003

143

[28] Gerris, S. Popinet., “A Tree-Based Adaptive Solver for the Incompressible Euler
Equations in Complex Geometries”, J. Comput. Phys, Vol. 199, pp.465-502,
2004.

[29] Godunov, S.K., “A F mite-Difference Method for the Numerical Computation of
Discontinuous Solutions of the Equations of Fluid Dynamics”, Mat. Sb., Vol. 47,
pp.271, 1959.

[30] Ham, F. and Young, Y.N., “A Cartesian Adaptive Level Set Method for
Two-Phase Flows”, Center for Turbulence Research, Annual Research Briefs
2003. 227, 2003.

[31] Harten, A., Engquist, B., Osher, S., and Charkavarthy, 8., “Uniform High-Order
Accurate Essentially Non-Oscillatory Schemes III”, J. Comput. Phys, Vol. 71,
pp.231-303, 1987.

[32] Hirt, CW. and Nichols, B.D., “Volume of Fluid (V OF) Method for the Dynamics
of Free Boundaries”, J. Comput. Phys, Vol. 39, pp.201-225, 1981.

[33] Hu, C. and Shu, C.W., “A Discontinuous Galerkin Finite Element Method for
Hamilton-Jacobi Equations”, ICASE Report No.98-2, NASA/CR-1998-206903,
1998. '

[34] J iang, GS. and Shu, C.W., “Efficient Implementation of Weighted ENO schemes”,
J. Comput. Phys, Vol. 126, pp.208-228, 1996.

[35] Jiang, GS. and Peng, D., “Weighted ENO Schemes for Hamilton-Jacobi
Equations”, SIAM J. Sci. Comput., Vol. 21, No.6, pp.2126-2143, 2000.

[36] Karman, S.L., “Splitflow: A 3D Unstructured Cartesian/Prismatic Grid CFD Code
for Complex Geometries”, AIAA Paper. 0343, J an., 1995.

[37] Kim, D. and Choi, H., “A Second—Order Time-Accurate Finite Volume Method
for Unsteady Incompressible Flow on Hybrid Unstructured Grids”, J. Comput.
Phys, Vol.162, pp.411-428, 2000.

[3 8] Kim, J. and Moin, P., “Application of a Fractional-Step Method to Incompressible
Navier-Stokes Equations”, J. Comput. Phys, Vol. 59, pp.308-323, 1985.

[39] Lafon, F. and Osher, S., “High Order Two Dimensional Non-Oscillatory Methods
for Solving Harrrilton-Jacobi Scalar Equations”, J. Comput. Phys, Vol. 123,
pp.235-253, 1996.

[40] Levy, D., Nayak, S., Shu, C.W., and Zhang, Y., “Central WENO Schemes for
Hamilton-Jacobi Equations on Triangular Meshes”, Scientific computing report,
BrownSC-2004-08, 2004.

[41] Liu, X., Osher, S., and Chan, T., “Weighted Essentially Non-Oscillatory schemes”,

144

J. Comput. Phys, Vol. 126, pp.202-212, 1996.

[42] Losasso, F., Gibou, F., and Fedkiw, R., “Simulating Water and Smoke with an
Octree Data Structure”, ACM Trans. Graph, pp 457-462, 2004.

[43] Losasso, F., Fedkiw, R., and Osher, S., “Spatially Adaptive Techniques for Level
Set Methods and Incompressible Flow”, Computers and Fluids, (in press).

[44] Mulder, W., Osher, S., and Sethian, J .A., “Computing Interface Motion in
Compressible Gas Dynamics”, J. Comput. Phys, Vol. 100, pp.209-228, 1992.

[45] Martin, J .C. and Moyce, W.J., “An Experimental Study of the Collapse of Liquid
Columns on a Rigid Horizontal Plan”, Philos. Trans. R. Soc. London, Ser. A 244,
312, 1952.

[46] Osher, S. and Sethian, J .A., “Fronts Propagating with Curvature Dependent Speed:
Algorithms Based on Hamilton-Jacobi Formulations”, J. Comput. Phys, Vol. 79,
pp.l2-49, 1988.

[47] Osher S. and Shu C.W., “High-Order Essentially Non-oscillatory Schemes for
Hamilton-Jacobi Equations”, SIAM J. Numer. Anal., Vol. 28, No. 4, 'pp.907-922,
1991.

[48] Osher, S. and Fedkiw, R., “Level Set Methods: An Overview and Some Recent
Results”, J. Comput. Phys, Vol. 169, pp.463-502, 2001.

[49] Osher, S. and Fedkiw, R., “Level Set Methods and Dynamic Implicit Surfaces”,
New York : Springer, 2003.

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

[51] Peyret, R. and Taylor, T., “Computational Methods for Fluid Flow”,
Springer-Verlag, NY, 1983.

[52] Popinet S. and Zaleski S., “A Front-Tracking Algorithm for Accurate
Representation of Surface Tension”, Int. J. Numerical Methods in Fluids, Vol. 30,
pp.775-793, 1999.

[53] Powell, K. G., “Solution of Euler and Magneto-Hydrodynamics Equations on
Solution-Adaptive Cartesian Grids”, von Karman Institute of Fluid Dynamics,
Leture Series 1994-05, Belgium, 1996.

[54] Puckett, E.G, Almgren, A., Bell, J., Marcus, D. and Rider, W., “A High Order
Projection Method for Tracking Fluid Interfaces in Variable Density
Incompressible Flows”, J. Comput. Phys, Vol. 130, pp.269-282, 1997.

[55] Puckett, E.G., “An Introduction to Numerical Method for Solving Incompressible

145

Viscous Flow Problems”, J. Comput. Phys, Vol.135, pp. 134-164, 1997.

[56] Rhee, C., Talbot, L., and Sethian, J.A., “Dynamical Study of a Premixed V
Flame”, J. Fluid Mech., Vol. 300, pp.87-115, 1995.

[57] Sethian, J .A. and Strain, J ., “Crystal Growth and Dendritic Solidafication”, J.
Comput. Phys, Vol. 98, pp. 231-253, 1992.

[58] Sethian, J.A., “Level Set Methods and Fast Marching Methods”, Cambridge
University Press, 1999.

[59] Sethian, J. A. and Smereka, P., “Level Set Methods for Fluid Interface”, Annu.
Rev. Fluid Mech., 352341—72, 2003.

[60] Shin, S. and Juric, D., “Modeling Three-Dimensional Multiphase Flow Using a
Level Contour Reconstruction Method for Front Tracking without Connectivity”,
J. Comput. Phys, Vol. 180. No. 2, pp.427-470, 2002.

[61] Singh, R., Dri, N., Uzgoren, E., Shyy, W., and Garbey M., “Three-dirnensional
Adaptive Cartesian Grid Method for Multiphase Flow Computation”, AIAA
paper No. 2005-1389, 2005.

[62] Sousa, F. S., Mangiavacchi, N., Nonato, L. G, Castelo, A., Tomé, M. F., Ferreira,
V. G, Cuminato, J. A., and McKee, S., A Front-Tracking/ Front-Capturing
Method for the Simulation of 3D Multi-Fluid Flows with Free Surfaces, J.
Comput. Phys, Vol. 198, No. 2, pp.469-499, 2004.

[63] Strain, J ., “Tree Methods for Moving Interfaces”, J. Comput. Phys, Vol. 151,
pp.6l6-648, 1999.

[64] Strain, J ., “Fast Tree-based Redistancing for Level Set Computations”, J. Comput.
Phys, Vol.152, pp.664-686, 1999.

[65] Strain, J ., “A Fast Modular Semi-Lagrangian Method for Moving Interfaces”, J.
Comput. Phys, Vol. 161, pp.512-536, 2000.

[66] Strain, J ., “A Fast Semi-Lagrangian Contouring Method for Moving Interfaces”, J.
Comput. Physl69, pp.l-22, 2000.

[67] Sussman, M., Smereka, P., and Osher, S., “A Level Set Approach for Computing
Solutions to Incompressible Two-Phase Flow”, J. Comput. Phys, Vol. 114,
pp.l46-159, 1994.

[68] Sussman, M., Fatemi, E., Smereka, P., and Osher, S., “An Improved Level Set
Method for Incompressible Two-Phase Flows”, Computers and Fluids, Vol. 27,
No. 5-6, pp. 663-680, 1998.

[69] Sussman, M., Almgren, A., Bell, J., Colella, P., Howell, L., and Welcome, M.,

146

“An Adaptive Level Set Approach for Incompressible Two-phase Flow”, J.
Comput. Phys, Vol. 148, pp.8l-124, 1999.

[70] Tryggvason, G, Bunner, B., Esmaeeli, A., Juric, D., Al-Rawahi, N., Tauber, W.,
Han, J ., Nas, S., and Jan, Y.-J., “A Front Tracking Method for the Computations
of Multiphase Flow”, J. Comput. Phys, Vol. 169 pp.708-759, 2001.

[71] Unverdi, SO. and Tryggvason, G, “A Front-Tracking Method for Viscous,
Incompressible, Multi-Fluid Flows”, J. Comput. Phys, Vol. 100, pp.25-37, 1992.

[72] Wang, Z.J., “A Quadtree Based Adaptive Cartesian Quad Grid Flow Solver for
Navier Stokes Equations”, Computers and Fluids, Vol.27, No.4, 1998.

[73] Wang, Z.J., “A Fast Nested Multi-grid Viscous Flow Solver for Adaptive
Cartesian/quad Grids”, International Journal for Numerical Methods in Fluids,
vol. 33, No. 5, pp. 657-680, 2000.

[74] Wang, Z]. and Zhao, Y., “A Characteristic Upwind Finite Volume Method for
Incompressible Flow and Heat Transfer on Arbitrary Grid”, Computer Methods in
Applied Mechanics and Engineering, (in review).

[75] Wang, Zhu and Wang, Z.J., “The Level Set Method on Adaptive Cartesian Grid
for Interface Capturing" AIAA paper, no.2004-0082, 2004.

[76] Wang, Zhu and Wang, Z.J., “Multi-Phase Flow Computation with
Semi-Lagrangian Level Set Method on Adaptive Cartesian Grids”, AIAA paper,
no.2005-1390, 2005.

[77] Whitney, H., “Analytic Extensions of Diflerentiable Functions Defined in Closed
Sets”, Trans. Am. Math. Soc., Vol. 36, pp.63-89, 1934.

[78] Zalesak, S., “Fully Multidimensional Flux-Corrected Transport Algorithm for
Fluids”, J. Comput. Phys, Vol. 31, pp.335-362., 1979.

[79] Zhang, Y. and Shu, C.W., “High Order WENO Schemes for Hamilton-Jacobi
Equations on Triangular Meshes”, SIAM J. Sci. Comput., Vol. 24, pp.1005-1030,
2003.

[80] Zhao, H.K., Chan, T., Merriman, B., and Osher, S., “A Variational Level Set
Approach to Multiphase Motion”, J. Comput. Phys, Vol. 127, pp.l79-195, 1996.

[81] Zhao, Y. and Zhang, B., “A High Order Characteristics Upwind Finite Volume
Method for Incompressible Flow and Heat Transfer Simulation on Unstructured
Grids”, Comput. Methods Appl. Mech. Engr, Vol. 190, pp.733-756, 2000.

[82] Zhao, Y., Tan, H. and Zhang, B., “A High-Resolution Characteristics-based
Implicit Dual Time-stepping VOF Method for Free Surface Flow Simulation on
Unstructured Grids”, J. Comput. Phys, Vol. 183, pp.233-273, 2002.

147

   

ll[[il[i[l[[iil[[i][i[il[i[iliil