Citation
A level set approach for two-dimensional, two-fluid flows with surface tension and viscosity

Material Information

Title:
A level set approach for two-dimensional, two-fluid flows with surface tension and viscosity
Creator:
Wilson, John
Publication Date:
Language:
English
Physical Description:
viii, 92 leaves : illustrations ; 28 cm

Subjects

Subjects / Keywords:
Two-phase flow -- Computer simulation ( lcsh )
Fluid dynamics -- Mathematical models ( lcsh )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 91-92).
General Note:
Department of Mechanical Engineering
Statement of Responsibility:
by John Wilson.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
44105082 ( OCLC )
ocm44105082
Classification:
LD1190.E55 1999m .W35 ( lcc )

Full Text
A LEVEL SET APPROACH FOR TWO-DIMENSIONAL, TWO-FLUID
FLOWS WITH SURFACE TENSION AND VISCOSITY
by
John Wilson
B.S.M.E., University of Colorado at Denver, 1997
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Mechanical Engineering
1999


This thesis for the Master of Science
degree by
John Wilson
has been approved
by
Ken Ortega
Date


Wilson, John (M.S., Mechanical Engineering)
A Level Set Approach for Two-dimensional, Two-Fluid Flows with Surface Tension
and Viscosity
Thesis directed by Assistant Professor Sam W. Welch
ABSTRACT
A level set approach is applied to two-fluid flows with surface tension and
viscosity. The level set method defines a smooth function, , as the signed normal
distance from the interface. The zero level set of (f> locates the interface. The
interface is defined with a set thickness that must be maintained throughout the
evolution of the flow. Flow properties are transitioned from the interior to the
exterior across the interface thickness. The curvature calculation used in determining
the surface tension is formulated using the level set function. The level set function
is advected with the flow field velocity. The flow field velocities and pressures are
computed using a semi-implicit version of the Mac method. The equations of motion
are discretized on a staggered grid. The fluids are modeled with constant densities
and viscosities. Simulations are conducted to ascertain the level set methods
capability of modeling flow situations in which there are analytical or experimental
results with which to compare. An oscillating bubble is simulated to verify the
correctness of the surface tension model. To further examine the surface tension
accuracy the Raleigh-Taylor instability is simulated. Buoyancy driven flows were
considered in flow regimes where there are experimental and computational results
available with which to compare. In addition simulations are conducted using the
in


level set method to determine its ability to accurately model two-fluid flows where
merging and breaking of the interface occurs.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication .
Signed
IV


ACKNOWLEDGMENT
I would like to thank Sam Welch for his guidance in everything that is
computational fluids. I would like to thank my fiance, Antoinette, for everything that
isnt computational fluids. I thank her for all of her support through the times that I
wasnt so pleasant to be around. I would like to thank her for hanging my bubbles on
the refrigerator as a sign of support and encouragement. I would also like to thank
my family that supported me through all of the lean years.


CONTENTS
Figures.....................................................................viii
Chapter
1. Introduction.............................................................1
2. Basic Equations..........................................................3
2.1 Control Volume Containing Interface......................................3
2.2 Level Set Formulation....................................................5
2.3 Free Boundary Condition at the Interface.................................6
2.4 Navier Stokes Equations..................................................9
3. Numerical Method........................................................12
3.1 Equations of Motion.....................................................12
3.2 Mac Method..............................................................13
3.2.1 Pressure Iteration......................................................17
3.3 Spatial Discretization..................................................17
3.3.1 Convection Terms........................................................18
3.3.2 Viscous Stress Terms....................................................20
3.3.3 Surface Tension Terms...................................................22
3.4 Fluid Properties........................................................24
3.5 Maintaining Level Set as a Distance Function............................24
3.6 Convergence of Level Set Method.........................................27
4. Numerical Simulations and Results.......................................28
4.1 Oscillating Bubble......................................................28
4.1.1 Derivation of Analytical Frequency......................................29
4.1.2 Numerical Results.......................................................33
4.2 Raleigh-Taylor Instability..............................................36
4.2.1 Derivation of Critical Wave Number......................................37
4.2.2 Numerical Results......................................................42
4.3 Buoyancy Driven Bubbles.................................................52
4.3.1 Necessity of Reinitialization Procedure.................................54
4.3.2 Spherical Bubble........................................................59
4.3.3 Ellipsoidal Bubble......................................................61
4.3.4 Dimpled Ellipsoidal-Cap Bubble.........................................63
4.3.5 Skirted Bubble.........................................................65
4.4 Merging and Breaking Interfaces.........................................68
vi


4.4.1 Dropping Bubble....................................................68
4.4.2 Two Bubbles........................................................76
4.4.3 Unstable Interface.................................................82
5. Conclusions and Recommendations....................................89
Bibliography............................................................91
vii


FIGURES
Figure
2.1 Control volume containing interface......................................4
2.2 Unit normals on interface................................................4
3.1 Momentum cells..........................................................14
3.2.2.1 Viscous stress on x-momentum cell.......................................21
3.2.3.1 Curvature stencil.......................................................23
4.1.1.1 Oscillating bubble, initial position....................................35
4.1.1.2 Oscillating bubble, half period.........................................35
4.1.1.3 Oscillating bubble, full period.........................................36
4.2.1.1 Perturbed interface.....................................................37
4.2.2.1 Enlarged interfacial disturbance........................................43
4.2.2.2 Velocity field, k = 50, t = 0.125.......................................43
4.2.2.3 Velocity field, k = 50, t = 0.25........................................44
4.2.2.4 Velocity field, k = 50, t = 0.375.......................................44
4.2.2.5 Pressure field, k = 50, t = 0.375.......................................45
4.2.2.6 Velocity field, k = 95, t = 0.25........................................46
4.2.2.7 Velocity field, k = 95, t = 0.375.......................................46
4.2.2.8 Pressure field, k = 95, t = 0.375.......................................47
4.2.2.9 Velocity field, k = 105, t = 0.25......................................48
4.2.2.10 Velocity field, k = 105, t = 0.375.....................................48
4.2.2.11 Velocity field, k = 105, t = 0.5.......................................49
4.2.2.12 Velocity field, k = 150, t = 0.25......................................50
4.2.2.13 Velocity field, k = 105, t = 0.375.....................................50
4.2.2.14 Velocity field, k = 150, t = 0.5.......................................51
4.3.1 Graphical correlation for buoyancy driven flows.........................53
4.3.1.1 No reinitialization, E = 1.0, M = 104..................................56
4.3.1.2 Reinitialization, E = 1.0, M = 10"4.....................................57
4.3.1.3 Reynolds number, E = 1.0, M = 10"4.....................................58
4.3.1.4 Interface thickness, no reinitialization, E = 1.0, M = 10"4.............58
4.3.1.5 Interface thickness, reinitialization, E = 1.0, M = 10-4................58
4.3.2.1 Spherical bubble shape and velocity field...............................59
4.3.2.2 Spherical bubble Reynolds number.......................................60
4.3.2.3 Spherical bubble pressure field.........................................60
4.3.3.1 Ellipsoidal bubble shape and velocity field.....................61
4.3.3.2 Ellipsoidal bubble Reynolds number.....................................62
viii


4.33.3 Ellipsoidal bubble pressure field........................................62
4.3.4.1 Dimpled ellipsoidal-cap bubble shape and velocity field..................63
43.4.2 Dimpled ellipsoidal-cap bubble pressure field............................64
43.5.1 Skirted bubble shape and velocity field..................................65
43.5.2 Skirted bubble pressure field............................................66
4.3.53 Skirted bubble shape and velocity field after break......................67
4.4.1.1 Dropping bubble velocity field, t = 0.1..................................70
4.4.1.2 Dropping bubble velocity field, t = 0.15.................................71
4.4.1.3 Dropping bubble velocity field, t = 0.17.................................72
4.4.1.4 Dropping bubble velocity field, t = 0.175................................73
4.4.1.5 Dropping bubble velocity field, t = 0.18.................................74
4.4.1.6 Dropping bubble velocity field, t = 0.2..................................75
4.4.2.1 Two bubbles , velocity field, t = 0.02.................................77
4.4.2.2 Two bubbles , velocity field, t = 0.05.................................78
4.4.23 Two bubbles , velocity field, t = 0.1..................................79
4.4.2.4 Two bubbles , velocity field, t = 0.15.................................80
4.4.2.5 Two bubbles , velocity field, t = 0.165................................81
4.43.1 Perturbed interface......................................................82
4.43.2 Unstable interface, velocity field, t = 0.1..............................83
4.4.33 Unstable interface, velocity field, t = 0.2..............................84
4.43.4 Unstable interface, velocity field, t = 03...............................85
4.43.5 Unstable interface, velocity field, t = 0.35.............................86
4.43.6 Unstable interface, velocity field, t = 0.37.............................87
IX


1.
Introduction
There are many numerical methods used for computing solutions for
incompressible two-fluid flow problems. One such method is the level set method.
In this work we investigate the level set method's ability to accurately model
incompressible two-fluid flows where surface tension plays an integral part in the
overall physics.
There are two general numerical approaches for solving incompressible fluid
flow problems with a moving boundary. One approach is based on front tracking
methods where the moving boundary is explicitly tracked. One popular front
tracking method is the Volume-of-Fluid (VOF) method. In this method a fluid
fraction variable,/ is defined and used to identify the position of the interface [1].
The fluid fraction variable is assigned values of 1 and 0 in the two phases. The
location of the interface is identified by fractional values of/ The VOF method
reconstructs the interface for each time step based on volume considerations. The
reconstruction step requires difficult geometrical calculations. The fluid fraction
variable is then tracked as it is advected with the local flow velocity.
The second method is based on front capturing. The front capturing method
considered here is the level set method introduced by Osher and Sethian [2], The
level set method involves defining a smooth function, (j>, as the signed normal
distance from the front. The fluid interface is the zero level set of $. The level set
function is advected with the local flow velocity. The advection of the level set has
the quality that it is uncoupled from the equations of motion for the two fluids and
thus doesnt introduce any greater complexity to the solution of the velocity and
pressure fields.
1


The level set method involves replacing the singular interface with an
interface of a finite thickness. Maintaining the prescribed interface thickness as the
interface evolves is paramount to obtaining accurate solutions. This is enforced by
preserving as a distance function. The benefits of defining ^ as a distance
function include
capturing the interface location as the zero level set of
the ability to transition smoothly from interior fluid
properties to exterior fluid properties
accurately calculating the curvature of the interface based
on the level set function
The implementation of the level set formulation presented here closely
follows the work of Sussman, Smereka, and Osher [3], The solution of the equations
of motion in [3] were obtained using the projection method. The solution of the
equations of motion in the work presented here is obtained using a semi-implicit
version of the Mac method.
In Chapter 2 the equations that describe the motion of the interface are
derived based on two immiscible fluids with constant fluid properties. In Chapter 3
these equations are discretized and the computational solution technique is outlined.
In Chapter 4 simulations are presented that probe the level set method's ability to
accurately model surface tension forces and deforming interfaces. The simulations
considered are an oscillating bubble, Rayleigh-Taylor instability, buoyancy driven
flows, and flows with merging or breaking interfaces. In Chapter 5
recommendations and areas of current improvements to the level set method are
discussed.
2


2.
Basic Equations
In this chapter we summarize the basic principles to fluid flows involving two
immiscible fluids. In addition, only flows where temperature remains constant are
considered, thus density, viscosity, and surface tension are assumed to be constant.
Considering these assumptions the only principles that need to be taken into account
are the conservation of mass and linear momentum. These equations will provide the
theory to solve for the velocity and pressure field for any time in question. Included
with these principles is the level set, which is discussed in section 2.2.
2.1 Control Volume Containing Interface
The conservation of mass and linear momentum will be applied to an arbitrary
control volume shown in figure 2.1.1. The subscripts / and g denote the liquid and
gas regions respectively. The liquid region has a volume of Q, bounded by the
surface dd, + 3Q, where 3Q, is the interface separating the two fluids. The same
nomenclature is used for the gas region. C/ represents the perimeter of 3Q;. Q
without a subscript is used to denote the union of the gas and liquid regions. Figure
2.1.2 shows the unit normal vectors directed out of their respective region. Na
represents the unit vector directed outward from Q in the tangent plane of 3Q,
3


I
o
figure 2.1.2 Unit normals on interface
4


2.2 Level Set Formulation
To locate the interface position a smooth level set function, f), is
introduced such that
dQx = {x| (f>{x,t) = 0}
(2.2-1)
Equation (2.2-1) states that if /) is the level set of the interface the zero level set
of distance from the interface. The level set function is taken as positive outside the
region bounded by the interface and as negative inside the region. The evolution of
the level set in time is given by
a
+ u V (f> = 0
(2.2-2)
Defining 4>(Jc, /) as a distance function provides a convenient means to determine the
location of the interface and as a consequence the fluid properties. The continuous
level set function provides a mechanism of transitioning from the interior fluid
properties to those of the exterior.
p{) =Pg+(p, ~ Pg )# M <2-2-3)
p{)= Pg + (pi-Pg)H{0)
(2.2-4)
5


Where H{ W=
0 if <0
1 if (f> > 0
1
[2
if = 0
(2.2-5)
An additional benefit of defining <}>(£, t) as the signed normal distance from
the interface is that it provides a simple means of calculating the curvature of the
interface which is discussed in Chapter 3.
2.3 Free Boundary Condition at the Interface
By applying the conservation of linear momentum a boundary condition at
the interface that relates surface tension to the jump in normal stress can be derived.
Below is the conservation of linear momentum for the entire control volume
presented in figure 2.1.1 neglecting the contribution from the interfacial momentum.
4 \PludV + ^ {pgudV = jPigdV + lpggdV+ (2.3-1)
a a, a og n, ng
\ T,dS + \ fgdS + | Tads
da, ent c,
g is the gravitational acceleration and represents the only body force considered.
For the gas and liquid regions:
F = T n
(2.3-2)
6


where T is the stress tensor.
T = -p\ + 2 pD (2.3-3)
D is the rate of strain tensor, I is the identity tensor, and p is the pressure.
For the interface:
fa = Ta K (2-3-4)
where T0 is the surface stress tensor given by the following expression when surface
viscous stresses are neglected.
T. = fla (2.3-5)
Where y is the surface tension and I0 is the surface projection tensor. The surface
projection tensor is like the identity tensor in three-dimensions.
Conservation of linear momentum for the liquid region may be expressed as:
4 \ p,udV = J p,gdV + \f,dS (2.3-6)
n, n, |
Conservation of linear momentum for the gas region may be expressed as:
7


(2.3-7)
4 J = \pgdV +
/2g nt xi'+an,
Substituting equations (2.3-2) ,(2.3-4), (2.3-6), and (2.3-7) into equation (2.3-1)
gives
|T, n,dS + jTg = jTa (2.3-8)
an, an, c,
Applying the surface divergence theorem [4] equation (2.3-8) becomes
j|T(-T,J-,dS=J an, <2(2,
Applying the surface divergence operator [4] to the surface stress tensor assuming
the surface tension to be constant results in the following boundary condition at the
interface.
PH = 2yKmnl (2.3-10)
Where |T| = T, -T2 and represents the jump in T across the interface. Km is the
mean curvature of the interface. Equation (2.3-10) shows that there is a discontinuity
in the normal stress proportional to the mean curvature of the interface. The mean
curvature is given by [5]
K=^(ici + k2),
(2.3-11)
8


where k, and k2 are the greatest and least curvatures. For a curve in two
dimensions the least curvature is zero. This results in
K = 2Km (2.3-12)
For two-dimensional problems the boundary condition at the interface becomes
|T| fi, = yml (2.3-13)
2.4 Navier Stokes Equations
Restating the total conservation of linear momentum with the surface
divergence theorem and equation (2.3-2) results in the following equation.
^\p(^dV=iMgdV+ (2.4-1)
a n n
JTrn,dS+ JTg-ngdS+ \yKn,dS
an, ans an,
Applying the Reynold's transport theorem and Greens theorem to equation (2.4-1)
results in
\p( n v 7 Dt n n an,
The last term on the right hand side of equation (2.4-2) may be shown to be [6]
\yK-h,dS = lyic(t)5{tptdV (2.4-3)
an, n
9


The curvature is given by
-W=v.
IwL
(2.4-4)
and Sif) is a one-dimensional Dirac delta function. Substituting equation (2.4-3)
into equation (2.4-2) gives
J[ - pMi V. T - = 0
n\ Ut )
(2.4-5)
Since this is true for an arbitrary volume,
pM pWs -'v t - = o
(2.4-6)
Substituting equation (2.3-3) results in the statement of the conservation of linear
momentum.
= pMg - +v
(2.4-7)
The densities of the liquid and gas are assumed to be constant which results in the
following expression for the conservation of mass.
V-u =0
(2.4-8)
10


Equations (2.4-7) and (2.4-8) were used to solve for the pressure and velocity fields
with suitable boundary conditions on 3D.
11


3.
Numerical Method
In this chapter we develop the numerical method that is implemented for the
two-dimensional flows considered. The equations of motion are solved using the
semi-implicit primitive version of the Mac method.
3.1 Equations of Motion
The equations of motion are listed below. Note that the dependence of the
fluid properties on the level set function is not shown for clarity.
Conservation of linear momentum in the coordinate directions x and z, respectively:
f du du du'
p + u + w--
\dt dx dz)
dx dx dz dx
(3.1-1)
(3.1-2)
where the stress tensor in terms of velocity gradients is given by


Conservation of mass for the incompressible fluids:
du dw n
+ = 0
dx dz
(3.1-4)
Advection equation for the level set function:
d d dd>
+ U + W
dt dx dz
0
(3.1-5)
The advection of the level set function is followed by a reinitialization procedure that
is discussed in section 3.6.
3.2 Mac Method
The Mac Method uses a staggered grid when discretizing the equations of
motion spatially. When considering incompressible flow the centered differences
used to solve the continuity equation allow for a checkerboard distribution of
velocities as a possible solution if the discretizations are not performed on a
staggered grid [7], Similarly, the centered differences used in calculating the
pressure gradients in the momentum equations could yield a checkerboard
distribution of pressures. In other words, centered differences performed on a grid
where pressures and velocities are discretized at the same locations could lead to
alternating values of pressures and velocities that satisfy the equations of motion but
make no physical sense. To eliminate the possibility of calculating erroneous
pressure and velocity fields all scalar fields (pressure, density, viscosity, and the level
set) are defined at cell centers and the velocities are defined at the centers of
momentum cells. Figure 3.2.1 shows a typical x-momentum and z-momentum cell.
13


Uj + 1
f
i-l-j T j i i'+ky -.7 f j

uj-i
figure 3.2.1 Momentum cells
To advance in time when given initial pressure, p, and velocity fields, un and wn,
the Mac Method consists of the following steps:
1. Advance in time by explicitly updating the velocities, un+l/2 and wn*112,
with the momentum equations based on the previous pressure and
velocity fields. This is an intermediate step and does not return the
correct velocity field. The intermediate velocities become
n"+//2 = un + At
i,j i,j
(
- conux conuz
+ (viscx + stenx)
V
PX
Y
i
'i.J
(3.2-1)
At
{P* )t.,Ax

14


I J
w"+j1/2 = w" + At g conwx conwz +(viscz + stem)
\ Pz
(3.2-2)
(3.2-3)
(3.2-4)
con refers to convection terms, vise refers to viscous stress terms, and sten
refers to surface tension terms. The average dynamic viscosity is used in
the viscous stress terms.
2. Implicitly solve for the change in pressure that satisfies the continuity
equation. Because the new velocity field from step 1 will not necessarily
satisfy the continuity equation a pressure equation is constructed from the
continuity equation. This results in a system of linear equations that can
be solved iteratively. Create the pressure equation in the following
manner.
(3.2-6)
15


The discretized continuity equation centered about pt j is
4'j
-u
Ax
n + l
<-.J
+

Az
(3.2-7)
Substituting equations (3.2-5) and (3.2-6) into (3.2-7) gives the pressure
equation
c,4>u + cm'j+c+wx
n+l
n+l
n+l
n+l
i +C6 0
(3.2-8)
^ uddp. j + uddPj_j j wddpt J + uddpt J_
Cy - +'
Ac
Az
(3.2-9)
C,
uddPlj
Ax
C}
uddp,_j j
Ax
C4
wddpt .
Az
c,
wddPij-!
Az
(3.2-10)
(3.2-11)
(3.2-12)
(3.2-13)
16


(3.2-14)
where uddp =----- and wddp =-----
pxAx pzAz
3. Update the velocity fields using equations (3.2-5) and (3.2-6). The
pressure field is updated using equation (3.2-15).
(3.2-15)
3.2.1 Pressure iteration
To calculate the change in pressure necessary to satisfy the continuity
equation the Biconjugate Gradient Stabilized method [8] is used. The linear system
to be solved is sparse, nonsymmetric and nonpositive definite. The Biconjugate
Gradient Stabilized method has the advantage of minimizing the build-up of
rounding errors associated with the Conjugate Gradient Squared algorithm and also
has faster convergence characteristics compared to the Biconjugate Gradient
algorithm [8].
3.3 Spatial Discretization
In this section we describe the discretization of the convection, viscous stress,
and surface tension terms found in the equations of motion. Boundary conditions are
implemented by including ghost cells around the perimeter of the solution cells.
17


3.3.1 Convection Terms
The convection terms in the momentum equations uses the Leonard [9]
method, which is weighted in the upstream direction as follows.
If ut J > 0
du
ij

dx
>.j
2 Ax
6 Ax
(3.3.1-1)
otherwise
du
i.j
' UMJ u,_, j u,+2 j 3uj+, j + 3uu ui_, ]
j
dx
i.j
2 Ax
6 Ax
(3.3.1-2)
The terms w, u, and w are implemented in a similar fashion. The
dz dx dz
Leonard method has the advantage of being second order and simple to implement.
The disadvantage of the Leonard method is it doesn't manage discontinuities well.
The level set has discontinuities in slope at the center of a circular bubble as the level
set of (/> reaches a minimum. As the interface evolves and becomes deformed many
discontinuities in slope may arise. The Leonard scheme would tend to smear out
these slopes removing the discontinuities. The Leonard scheme may also exhibit
oscillatory behavior near the discontinuities. In some instances the discontinuity may
lie within the prescribed boundary layer thickness causing the level set function to
deviate from a distance function. The smearing of the interface or oscillating values
of (f> within the prescribed interface thickness would result in erroneous values of
18


surface tension due to the incorrect calculation of curvature. The Leonard scheme
would also lead to the improper distribution of fluid properties across the interface.
For these reasons a high-resolution scheme is used to convect the level set function.
For the convection terms in the level set equation (3.1-5) a second order ENO
scheme is used. The ENO scheme more accurately maintains the discontinuities in
the gradients of the level set and essentially exhibits no oscillatory behavior, which is
critical for accurately modeling the physics near the interface. The ENO scheme that
is used is described below.
u w u
d<(>
dx
ave
i-l,i . ,/
- + minmod(a, b)
(3.3.1-3)
where
minmo
(3.3.1-4)
a =
h+i.j + 0,-ij
Ax2
(3.3.1-5)

Ax2
(3.3.1-6)
u
ave
(3.3.1-7)
19


The sgn in equation (3.3.1-4) represents the sign function.
If u
ave
<0
dtp
u
dx
Ax
Ax .
----minmo
2
\
d(a,c)
7
(3.3.1-8)
where
ti+ij 2>+i,j+hj
Ax2
(3.3.1-9)
i
A similar approach is taken for the convection term in the z direction, w.
dz
3.3.2 Viscous Stress Terms
Centered differences were used to calculate the viscous stresses on x and z-
momentum cells. The shear and normal stresses on a typical jc-momentum cell are
shown in Figure 3.3.2.1.
20


X
high
--------
w
W
------
r
** low
figure 3.3.2.1 Viscous stress on x-momentum cell
Below are the viscous stress terms used in the x-momentum equation.
"rif*
2m.
Ui + I.I Ui.J
right
Ax

2 Ml
left
U,.,-Ui-I.j
Ax
w,+u-i-wt.j-i u>j
~ MlO* I Ay ZlZ
Tzxhsh Mhigh
W,+lj-Wi.J Ui.J+!~U>-L
Ax
Az
gT t ?
L XX X*Tight____XXftJl
dx Ax
ZX ^ ^high ZXltrit
dz Az
(3.3.2-1)
(3.3.2-2)
(3.3.2-3)
(3.3.2-4)
(3.3.2-5)
(3.3.2-6)
21


Note that the viscosity is a function of the level set and varies through the
thickness of the interface. On the top and bottom of the x-momentum cell the
viscosity is taken to be the average of the viscosities of the four neighboring cells. A
similar approach is taken for the z-momentum viscous terms.
3.3.3 Surface Tension Terms
To calculate the curvature in the surface tension terms centered differences
were used. The curvature as a function of the level set is given by
Figure 3.3.3.1 shows the stencil used to discretize the curvature. For the x-
momentum equation the average curvature of the i,j and the i+l.j cell is used.
Similarly for the z-momentum cell the average curvature of the i,j and the i,j+l cell
is used.
The curvature at a typical cell ij is calculated as follows:
(3.3.3-1)
^ 01+1,j 01-1, j
2 Ax
(3.3.3-2)
2Az
(3.3.3-3)
(3.3.3-4)
0ij+i ~ 20,.j + 0ij-i
Az2
(3.3.3-5)
22


xz ~
2Az
01+i.j+i fii-ij+i fii+ij-i &-1.J-1
2 Ax
2 Ax
(3.3.3-6)
The Dirac delta function is defined as [10]:
s ((f) = j(7 + C05(^/f))/(2f) if M
< £
0 otherwise
(3.3.3-7)
where e is half the overall interface thickness. Defining the Dirac delta as a function
of the level set ensures that surface tension will only be included in the momentum
equations within the interface thickness.
0i,j+i
ti-IJ A. 1 > , ' 1J
+u-i
Figure 3.3.3.1 Curvature stencil
23


3.4 Fluid Properties
In order to solve the equations of motion numerically the discontinuities in
density and viscosity are modeled by a region where the properties transition
smoothly from the interior fluid properties to the exterior fluid properties. This is
accomplished by modeling the interface as a region of finite thickness. In our
simulations the interface thickness, e, is set at 1.5 Ax on either side of the zero level
set. This provided a smooth transition from the gas properties to those of the fluid
over a distance of 3 Ax To implement this numerically a regularized Heaviside
function is used [10]

0 if <-e
0.5(] + /e + sin{n(j)Ie)ln) if |^| < e
1 if (f> > s
(3.4-1)
Using the regularized Heaviside function the continuous density and viscosity
distributions become
Pc M = Pg + (Pi ~ Pg ]Hs U>) (3.4-2)
pA 3.5 Maintaining Level Set as a Distance Function
While equation (3.1-5) convects the zero level set and consequently the
interface at the correct velocity the level set will not necessarily remain a distance
function (|V^| f)[3]. It is essential to maintain the level set as a distance function
in order for the interface to remain an invariant thickness. If the defined interface
thickness is not maintained the smooth transition between fluid properties will be
24


compromised. Also the surface tension force will become inaccurate due to the
incorrect calculation of the curvature as well as not being included in the proper
momentum equations. If the level set deviates from a distance function w'ithin the
interface thickness the level set doesn't reflect the actual discretized interface
thickness, 3Ax. Because of this the Dirac delta function would distribute the surface
tension force over more or less grid cells than the prescribed number. To help
remedy this problem a reinitialization procedure introduced in [3] is implemented. It
requires evolving the following equation to steady state
(t>,+cb-V given the initial level set distribution . S is the sign function and
(3.5-2)
A high-resolution scheme for evolving equation (3.5-1) to steady state is used [3]. A
high-resolution scheme is used for the reasons discussed in section 3.3.1. The
discretization is as follows
a
Ax
h+i.j -0'j
Ax
A j
Az
,j+i ~ Aj
Az
(3.5-3)
25


(3.5-4)
sAA, -r^
+ Ax
G(t),
J max((a+ f, (b_ f)+ max((c+ J, (d~ f) -1
if > ^J(Ax)2 +(Az)2
-J max ((a" f, (b+ f)+ max((c~ f, (d+ f|-1
if C < {hxf + (Az)2
0 otherwise
(3.5-5)
The + and superscripts denote the inclusion of that difference only if it is greater
than or less than zero respectively.
Equation (3.5-1) is discretized temporally in the following manner
r=p-a,s,(tiM,) <3'5-6)
This algorithm has the property that the zero level set of remains unchanged.
Equation (3.5-6) doesn't represent any physics of the solution; it is simply a
convenient method of maintaining a distance function. Equation (3.5-6) is iterated
with an artificial time step of At = Ax/10. Since it is only necessary for the level set
to be a distance function near the interface the solution is considered to be converged
when the error within the interface thickness is within our tolerance. For our
simulations we used an iteration tolerance of AtAx2 [3].
26


3.6 Convergence of Level Set Method
The numerical method used to solve the basic equations is first order in time
and second order in space. The advection of the level set is first order in time and
second order in space using an ENO method. The reinitialization procedure used is
first order. Convergence studies have been performed for the level set method [3]
and have shown the method to converge. For this reason a convergence study of the
level set method is not performed.
27


4.
Numerical Simulations and Results
In this chapter we consider different two-dimensional flow simulations to
demonstrate the effectiveness as well as the weakness of the level set method. The
flows considered are:
oscillating bubble
Raleigh-Taylor instability
buoyancy driven flows
merging and breaking interfaces
With these simulations we will be able to compare our results with
experimental data or analytical theory.
4.1 Oscillating Bubble
In this simulation we consider a gas bubble with a circular interface in
equilibrium with the surrounding fluid. The circular profile is then slightly perturbed
creating an unsteady flow state. Gravity is assumed to be out-of-plane and the flow
is assumed to be free of viscous stresses. For this problem surface tension will
dominate the physics and will serve to affirm the level set method's ability to model
surface tension forces. More specifically, the level set method's ability for accurately
calculating the curvature for a deforming interface. For these flow conditions the
perturbed circular interface will tend to oscillate from the original profile where the
major axis is directed along the horizontal to a position where the major axis is
directed in the vertical direction. In the absence of viscous stresses this oscillation
would never cease.
28


4.1.1 Derivation of Analytical Frequency
Sir Horace Lamb [11] considered the oscillations of a slightly perturbed
column of liquid. Below is a derivation using a similar approach applied to a slightly
perturbed column of gas surrounded by liquid. The flow field is considered to be
two-dimensional at an arbitrary plane that bisects the column of gas. The reasoning
by Lamb was that the uniform flow directed along the axis of the column will not
contribute to the dynamics of the interface oscillations.
For this incompressible, inviscid, two-dimensional flow the equations of
motion can be described in terms of velocity potentials. Applying the conservation
of mass to the interior /, and the exterior, e, of the interface yields
V20, = 0 (4.1.1-1)
V20e = 0 (4.1.1-2)
Where in cylindrical coordinates
d20 1 d0 1 d20 n dr2 r 8r r2 dB2 (4.1.1-3)
Note that the z direction is directed along the axis of the column of gas and is not
included in equation (4.1.1-3).
The conditions imposed on equations (4.1.1-1) and (4.1.1-2) are that the
velocity potential, 0, must vanish as r 0 and r qo as well as satisfy the
boundary condition at the interface. The unperturbed bubble radius is given by a.
29


r-a
r=a
(4.1.1-4)
80,
dr
30.
dr
The momentum equation for this flow reduces to the unsteady Bernoulli
equation for the interior and the exterior of the interface.
M + i + i(v dt p, 2
^ + + -j(V<0j =c,(/)
dt pe 2
Along with the jump in pressure across the interface
r
p.-p,=-R
Where y is the surface tension and R is the curvature of the interface. The curvature
is given by [11]
(4.1.1-6)
(4.1.1-7)
11 1 d2r
R ~ r r2 de2
The bubble radius is perturbed a distance £ where
r = a + £
(4.1.1-8)
(4.1.1-9)
Linearizing equation (4.1.1 -8) with the use of equation (4.1.1 -9) gives
30


(4.1.1-10)
1 1 if d2C
R~ a a2\^ + d0\
The solution for equation (4.1.1-1) follows directly from Lamb's analysis and
is given by
0' = A
/ \n
r '
\a)
cos(n6>)cos(cot + s)
(4.1.1-11)
Using equation (4.1.1-11) along with the boundary condition at the interface (4.1.1-
4) an expression for the exterior potential field can be derived.
01 =-A
a '
\r )
cos(n&)cos(a)t + e)
(4.1.1-12)
Assuming that particles on the interface must remain on the interface provides a
relationship between the velocity potential and the perturbation displacement.
80' _d£
& dt
(4.1.1-13)
Applying the above relationship to equation (4.1.1-12) leads to an analytical
expression for the perturbation displacement.
31


(4.1.1-14)
£ = Acos(n0)sin(aX + e)
coa
Linearizing the unsteady Bernoulli equation by assuming

by
0 = 0+ 0'
p = p + p'
(4.1.1-15)
gives
^ ^ + £~ + |(v( + 0')J = C(/) (
Linearizing equation (4.1.1-16) about the reference state gives the following
relationships at the interface
ri=_d0l
P, St f=a
P\ 90'e
Pe ^ ,mg
Substituting equations (4.1.1-10), (4.1.1-17), and (4.1.1-18) into equation (4.1.1-7)
gives
(4.1.1-17)
(4.1.1-18)
32


Equation (4.1.1-19) can be simplified by noting that
(4.1.1-20)
Substituting equations (4.1.1-11), (4.1.1-12), (4.1.1-14), and (4.1.1-20) into equation
(4.1.1-19) gives
which can be rearranged to give an expression for the frequency of the various
modes
4.1.2 Numerical Results
For this simulation the interface was initialized to be the shape of an ellipse
with the major axis directed along the horizontal. This interface shape corresponds
to the first mode of oscillation (n = 2). Body forces and viscous effects are
neglected. This simulation used a 64x64 grid with the interface perturbed one cell
width. This resulted in a major axis of 0.055 m and minor axis of 0.045 m. The
numerical values used in this simulation are
(4.1.1-21)
(4.1.1-22)
33


a 0.05 m
pe = 8.0 x 106 Pa
y = 1.582 x 10'2 N/m
p, = 42.417 kg/m3
pe = 722.02 kg/m3
The angular frequency using equation (4.1.1-22) is calculated to be
co = 0.99667 rad/sec
which results in a period of
t = 6.30416 sec
Figure 4.1.1.1 shows the initial shape of the bubble. Figure 4.1.1.2 shows the
velocity field and interface position at one half period. Figure 4.1.1.3 shows the
velocity field and interface position at one full period. The figures show that the
interface position has the shape and orientation predicted by the linearized theory.
What isn't evident from the figures is that there was an appreciable loss of area as the
bubble evolved in time. A calculation of the area at the first half period and whole
period resulted in the following loss.
A_
i = t/2
0.92
A_
Ao
T
*0.88
34


This loss of area can be attributed to the fact that nothing in the level set formulation
used here acts to preserve the bubble area. The level set function is advected at the
local flow velocity, the zero level set of formulation restricts the new interface to surrounding the same amount of area.
\
J
i
figure 4.1.1.1 Initial position
figure 4.1.1.2 Half period, ncyc = 500
35


figure 4.1.1.3 Full period, ncyc = 1000
4.2 Raleigh-Taylor Instability
In this simulation we consider the stability of two homogenous fluids of
different densities superposed one over the other. When the two fluids have no
relative horizontal motion the instability of the interface separating the two fluids is
known as the Raleigh-Taylor instability [12]. The equilibrium of the interface will
be investigated by perturbing the interface with disturbances of various wave
numbers. The theory will show that if the denser fluid resides underneath the less
dense fluid the interface is stable for all wave numbers. In this simulation we
consider two inviscid fluids of different densities with the denser fluid above the less
dense fluid Surface tension will be included in the analysis and will serve to
reaffirm our confidence in our calculation of the curvature of the interface.
36


The theory will show that the presence of surface tension has a stabilizing effect on
the interface. For the case of the denser fluid above the less dense fluid the theory
shows that the interface will remain stable for wave numbers above the critical wave
number which is derived in section 4.2.1.
4.2.1 Derivation of Critical Wave Number
Consider two incompressible, inviscid fluids shown in figure 4.2.1.1 denoted
by the subscripts 1 and 2. Initially the two fluids are at rest. The interface is
perturbed from the equilibrium position by rj.

z = 00
figure 4.2.1.1 Perturbed interface
The equations that describe the velocities of the two fluids are
V2, = 0
(4.2.1-1)
V22 = 0
(4.2.1-2)
with the far field boundary conditions
37


0,\ =0
lz=-oo
(4.2.1-3)
02\ =0 (4.2.1-4)
Assuming that particles on the interface will remain on the interface results in the
following condition at the interface
d0 = a1a0 + a1 (42|5)
dz dx dx dt
Equation (4.2.1-5) is linearized about z = 0 where the reference velocity
potential, 0O, is perturbed by 0'. Note that the fluid was initially at rest, which
results in a zero reference velocity potential.
d0\ _ drj
dz 2=0 dt
d0'2 _ drj
dz 2=0 dt
(4.2.1-6)
(4.2.1-7)
The unsteady Bernoulli equation is given by
?f+L(v0y+£-+gz = c(t)
dt 2 p
(4.2.1-8)
The pressure becomes
P = Po+P'
(4.2.1-9)
38


Where p0 is the reference pressure and p' is the perturbation in pressure. Linearizing
equation (4.2.1-8) at z = rj gives
d0' p0 + p /x + y y +gT] = C(t) ot p The reference condition at z = 0 is (4.2.1-10)
^ = C(t) P (4.2.1-11)
Combining equations (4.2.1-10) and (4.2.1-11) gives a relationship between the
perturbed pressure and velocity potential for each fluid.
, (80 ) Pi= P, J+gP \ dt ; (4.2.1-12)
' (d02 \ p,=-p{ s, p) and p'2 are related to each other by the surface tension (4.2.1-13)
p\-p\= Linearizing the curvature, k, considering only small amplitude interfacial disturbances gives [5] 82t] P2 Pi=<* -2 ox (4.2.1-15)
39


Solutions to equations (4.2.1-1) and (4.2.1-2) which satisfy the far field boundary
conditions (4.2.1-3 ) and (4.2.1-4) are assumed to be
0\ = QeV**""0
0'2 = C2e-lae,{kx-*)
With the disturbance of the interface given by
77 =
ae
(kx-ox)
(4.2.1-16)
(4.2.1-17)
(4.2.1-18)
Combining equations (4.2.1-16) through (4.2.1-18) with equations (4.2.1-3) and
(4.2.1-4) gives
a co .
C'=~T
~ aco .
C'=T*
(4.2.1-19)
(4.2.1-20)
Combining equations (4.2.1-12) and (4.2.1-13) with equation (4.2.1-15) gives
P2
d02
~dT
+ grj
+ Pi
d0\
~dT
+ gV
= a
l
dx2
(4.2.1-21)
Combining equations (4.2.1-16) through (4.2.1-18) with equation (4.2.1-21) gives
- p2(-icoC2 + ga)+ p^-icoC, + ga) = -oak2 (4.2.1-22)
40


Combining equations (4.2.1-19) and (4.2.1-20) with equation (4.2.1-22) gives
{ico)2 + (p, -p2)g = -ok2
k
(4.2.1-23)
Solving for (ico)2 gives
(ico)2 = gk
P2-P1
P2+ Pi
ok2
g(p2+Pi)
(4.2.1-24)
From equation (4.2.1-18) the interface will remain stable if
p2 ~Pl~r~~-------\<0 (4.2.1-25)
P2+P, g(P2+Pl)
which is always the case if p2 < p,. For the case p2 > p, stability is only
maintained for wave numbers greater than the critical wave number, kc.
K =
g(p2 ~Pj)
IP
(4.2.1-26)
From equation (4.2.1-25) it can be seen that without the effect of surface
tension the interface would be unstable for all wave numbers if p2 > pr The
presence of surface tension will stabilize the interface for sufficiently short
wavelength disturbances. The presence of surface tension also creates a mode of
41


maximum instability, k,, which can be shown by differentiating equation (4.2.1-24)
4.2.2 Numerical Results
For this simulation the interface was initialized as a sine wave as shown in
figure 4.2.2.1. For clarity the amplitude of the wave in figure 4.2.2.1 is greatly
enlarged. The amplitude of the disturbance for all of the simulations is 10^ m. The
denser fluid is superposed over the less dense fluid with gravity directed down.
Viscous stresses were neglected. The simulations used a 64x64 grid with periodic
boundary conditions. The numerical values used in these simulations are
which results in a critical wave number of kc = 100. Simulations for the wave
numbers 50, 95, 105, and 150 were conducted. The grid size (Ax, Az) was varied so
that exactly one wavelength was contained in the computational domain. From the
figures that follow it can be seen that the results are in complete agreement with the
analytical theory. The interface grows more rapidly for a wave number of 50 as
compared to 95, which agrees with equation (4.2.1-27). A wave number of 50 is
closer to the mode of maximum instability than a wave number of 95. Also at wave
numbers of 105 and 150 the interface appears to be completely stable.
[12].
(4.2.1-27)
P2 = 11.0 kg/m3
Pi = 1.0 kg/m3
g = 10.0 m/sec2
y = 0.01 N/m
42


J
K
. i-
figure 4.2.2.1 Enlarged interfacial disturbance
figure 4.2.2.2 Velocity field, k = 50, t = 0.125 sec
43


30 tO SO (C
figure 4.2.2.3 Velocity field, k = 50, t = 0.25 sec
figure 4.2.2.4 Velocity field, k = 50, t = 0.375 sec
44


figure 4.2.2.5 Pressure field, k = 50, t = 0.375 sec
45


figure 4.2.2.6 Velocity field, k = 95, t = 0.25 sec
r
" i
I
t
t
- jr
to
X
so
figure 4.2.2.7 Velocity field, k = 95, t = 0.375 sec
46


figure 4.2.2.8 Pressure field, k = 95, t = 0.375 sec
47


figure 4.2.2.9 Velocity field, k = 105, t = 0.25 sec
(0
-j
I
1
i
0
figure 4.2.2.10 Velocity field, k = 105, t = 0.375 sec
48


I___________________________I___________________________I____________________________I_____________________________I_____________________________I_____________________________L
10 JO JO 40 JO K
figure 4.2.2.11 Velocity field, k = 105, t = 0.5 sec
49


so r
P
- 1
figure 4.2.2.12 Velocity field, k = 150, t = 0.25 sec
so h
i
figure 4.2.2.13 Velocity field, k = 150, t = 0.375 sec
50


to
figure 4.2.2.14 Velocity field, k = 150, t = 0.5 sec
51


4.3 Buoyancy Driven Bubbles
In this simulation we consider a gas bubble surrounded by an infinite liquid.
Gravity is in the plane of motion, which causes the bubble to rise due to buoyancy
forces. Someone who is not acquainted with this event might consider this to be
uninteresting. However, this simulation proves to be the most exciting and
beneficial of the four that were considered. The motion caused by the combination
of buoyant, viscous, and surface tension forces cause the bubble to undergo
deformations as well as rise in the liquid. Experimental results of this phenomenon
were documented by Grace [13] and provide a valuable resource for comparison.
Grace condensed the results into a graphical representation reproduced in figure
4.3.1. The shape of the bubble and the steady state Reynold's number (Re) are
presented as functions of the Morton (M) and Eotvos (Eo) numbers.
u sa'-4g
pfr3
Eo = ^~
(4.3-1)
Y
Re =
P,deU
Pi
(4.3-2)
(4.3-3)
Ap represents the difference in density between the surrounding liquid and
the gas. The subscript, /, denotes the liquid region. The equivalent diameter, de,
appearing in the Reynold's and Eotvos numbers represents the diameter of a volume
equivalent sphere. In our simulations the initial bubble diameter is used for the
equivalent diameter.
52


figure 4.3.1 Graphical correlation for buoyancy driven flows
53


For given Eotvos and Morton numbers figure 4.3.1 predicts the bubble shape
and terminal bubble Reynold's number. These graphical results are based on three-
dimensional bubbles with infinite boundaries but still serve as a very useful
comparison. As an additional comparison the results presented here are compared
with the two-dimensional computational results by Tryggvason and Unverdi [14].
A total of four simulations were conducted for different combinations of the
Eotvos and Morton numbers. All of the simulations had free slip boundary
conditions with a grid resolution of 65x130. They also had a density ratio, pi/pg- 40
and surface tension, y= 0.01 N/m. In order to make an accurate comparison with
Tryggvason and Unverdi the same viscosity ratio and density ratios are used.
All of the simulations agreed favorably with figure 4.3.1 as well as the results
by Tryggvason and Unverdi. The breaking of the skirted bubble can be attributed to
a lack of grid resolution. When the skirt thickness approaches the order of the
prescribed interface thickness it is difficult to assess the validity of the solution. All
of the simulations had a significant loss of bubble area. The loss of area can be
attributed to the use of only a first order reinitialization scheme. A second and
perhaps more significant, reason is that nothing in the level set formulation specifies
that the bubble area should remain constant as the interface evolves.
4.3.1 Necessity of Reinitialization Procedure
It has been stated that while the level set is convected at the correct velocity it
will not remain a distance function. To illustrate the need for the reinitialization
procedure a comparison is made between two flow simulations. Figure 4.3.1.1
represents the flow field and bubble position at 0.2 seconds without the
reinitialization procedure. Figure 4.3.1.2 shows the flow field and bubble position at
0.2 seconds with the reinitialization procedure. A comparison between Reynold's
numbers is shown in figure 4.3.1.3. Figures 4.3.1.4 and 4.3.1.5 show the interface
position and thickness for the two cases.
54


For Eo = 1.0 and M = 10-4 a spherical steady state bubble shape should be
attained with a terminal Reynold's number of about three. Without the
reinitialization procedure at a time of 0.2 seconds the bubble interface has already
deformed significantly. The Reynold's number increases faster for the case without
reinitialization. For a time beyond 0.2 seconds the Reynold's number increases at an
ever-increasing rate well beyond the predicted terminal Reynold's number. From
figure 4.3.1.5 it is evident that the reinitialization procedure does indeed maintain an
invariant interface thickness. It is apparent from figure 4.3.1.4 that the interface
thickness increased on the trailing edge of the bubble and decreased on the leading
edge. It is difficult to explain the cause of such massive deterioration of the uniform
interface thickness. However, from the contours shown in figure 4.3.1.4 it can be
seen that the fluid properties on the leading edge would almost be discontinuous as
they are transitioned across the minute interface thickness. Also, on the trailing edge
the surface tension force is distributed across many more grid cells by way of the
Dirac delta function. This results in more surface tension force on the trailing edge
of the bubble than the leading edge. These problems combine to give an inaccurate
solution for the flow regime being modeled and illustrate the need for a
reinitialization procedure.
55


120
io 20 x o so ao
figure 4.3.1.1 no reinitialization, E = 1.0, M = 10'4
56


130
100
0
I
i
-i
10
0
so
figure 4.3.1.2 reinitialization, E = 1.0, M = 10"4
57


time (sec)
figure 4.3.1.3 Eo = 1.0, M = 10-4
i * i
[
i
l
figure 4.3.1.4 Eo = 1.0, M = 10"4
no reinitialization
figure 4.3.1.5 Eo = 1.0, M = 1CT4
with reinitialization
58


4.3.2 Spherical Bubble
Eo = 1.0 M = 10~4 = 493 de = 0.05 m

figure 4.3.2.1 Bubble shape and velocity field


:s >-
Re ..

/
/
os t- /
0 ^--------------------------------------------------------------
o oi o? os o< os oe o?
time (sec)
figure 4.3.2.2 Bubble Reynolds number
figure 4.3.2.3 Pressure field
60


4.3.3 Ellipsoidal Bubble
Eo = 10.0 M = 10'1 = 493 de = 0.05 m
~ 0.85
A0
to
I
figure 4.3.3.1 Bubble shape and velocity field
61



i
/
t
r /
j
time (sec)
figure 4.3.3.2 Bubble Reynolds number
figure 4.3.3.3 Pressure field
62


4.3.4 Dimpled Ellipsoidal-Cap Bubble
Eo = 100.0 M = 100.0 = 479 de = 0.1 m
0.85
Ao
10 JC C K
figure 4.3.4.1 Bubble shape and velocity field
63


64


4.3.5 Skirted Bubble
Eo = 104.0 M = 10_1
iii = 85 de = 0.1 m
r
figure 4.3.5.1 Bubble shape and velocity field
65


figure 4.3.5.2 Pressure field
66


I
100 j-
10
figure 4.3.5.3
20
30 40 SO
ao
Bubble shape and velocity field after break
67


4.4 Merging and Breaking Interfaces
For these simulations we consider interfaces that will merge or break to
determine the level set methods ability to track the event. We consider a liquid
bubble dropping into a pool of liquid in section 4.4.1. In section 4.4.2 we consider
two bubbles of equal diameters one above the other. The two bubbles will merge as
the lower bubble catches up to the higher bubble while traveling in its wake. In
section 4.4.3 a flat interface is perturbed such that it will become unstable causing
some of the fluid to break away.
4.4.1 Dropping Bubble
To show the level set methods ability to handle merging interfaces a liquid
bubble was dropped into a pool of liquid. Initially the bubble is surrounded by gas
with gravity directed down. The simulation was conducted on a 65x130 grid with
free slip boundary conditions. Viscous stresses were included in this simulation.
The numerical values used in this simulation are
Pi = 42.417 kg/m3
Pi = 722.02 kg/m3
Pi = 0.10098534 Nsec/m:
p
Pi = 479.0
y = 0.01 N/m
g = 10.0 m2
d = 0.05 m
68


The level set method handled the merging of the interfaces rather well
considering the grid resolution. Figure 4.4.1.4 shows that some of the gas was
trapped during the merger of the two interfaces. At a time of 0.18 sec (figure 4.4.1.5)
the gas has disappeared. It is difficult to say whether the gas should have been
trapped inside the liquid. However, it shouldnt have disappeared once it was
trapped. This can be attributed to the grid resolution and the inefficiency of the
reinitialization procedure at preserving area.
69


h
i
i
i

A
I
I
20
I
1
to
M
40
SO
0
figure 4.4.1.1 Velocity field, t = 0.1 sec
70


I-------------------------1----------------------------1-----------------------------1----------------------------1----------------------------:-------------------------------1------1
!
i
t- -i
100
1
10
M
40
SO
0
figure 4.4.1.2 Velocity field, t = 0.15 sec
71


10 30 30
SO
ao
figure 4.4.1.3 Velocity field, t = 0.17 sec
72


T
120
100
I
20 i
1 i
10 20 90 40 M *0
figure 4.4.1.4 Velocity field, t = 0.175 sec
73


120
10 20 JO 40 90 ao
figure 4.4.1.5 Velocity field, t = 0.18 sec
74


figure 4.4.1.6 Velocity field, t = 0.2 sec
75


4.4.2 Two Bubbles
In this simulation two gas bubbles of equal diameters were placed so that one
was one half of a bubble radius above the other. Both bubbles are surrounded by
liquid with gravity directed down. Both bubbles begin to rise and deform as the
solution evolves. This simulation was conducted on a 65x130 grid with free slip
boundary conditions. Viscous stresses were included in this simulation. The
numerical values used in this simulation are
pg = 42.417 kg/m3
Pi = 722.02 kg/m3
pg = 0.10098534 Nsec/m:
p
pg = 479.0
Y = 0.01 N/m
g = 10.0 m2
d = 0.05 m
The merging of the bubbles was depicted well by the level set method as
shown in the figures. A similar simulation was conducted by Chang et al. [6] with
two bubbles of different diameters. The deformation of the bubbles in this
simulation closely match their results.
76


1I
i

10 JO
SO BO
figure 4.4.2.1 Velocity field, t = 0.02 sec
77


IN
100
M
10 N N
40
SO
ao
figure 4.4.2.2 Velocity field, t = 0.05 sec
78


I-
tO SO 40 SO
figure 4.4.2.3 Velocity field, t = 0.1 sec
79


I
I
j
10
figure 4.4.2.4 Velocity field, t = 0.15 sec
80


120
h
i
i
100
0
10 20 K 40 0C
figure 4.4.2.5 Velocity field, t = 0.165 sec
81


4.4.3 Unstable Interface
In this simulation a flat interface was perturbed as shown in figure 4.4.3.1. A
grid size of 65x130 is used with free slip boundary conditions. The perturbation has
a wave number of k = 50 with an amplitude of 102 m. Viscous stresses are neglected
in this simulation. The purpose of this simulation is to inspect the level set methods
ability of handling a breaking interface. The numerical values used in this simulation
are
pg =0.5 kg/m3
p{ = 10.0 kg/m3
Y = 0.01 N/m
g = 10.0 m2
Pi
Pg
figure 4.4.3.1 Perturbed interface
82


T
T
T
I
i
j-
10 20 X
40 M <0
figure 4.4.3.2 Velocity field, t = 0.1 sec
83


1
10
figure 4.4.3.3 Velocity field, t = 0.2 sec
84


T
12D
!
100
10
ao


ao
figure 4.4.3.4 Velocity field, t = 0.3 sec
85


120
20 X
40 SO
0
figure 4.4.3.5 Velocity field, t = 0.35 sec
86


T
87


In this simulation there are no results with which to compare. However, from
the figures it can be seen that the level set method had no problem managing the
breaking interface. There was no significant change in the area considering the
degree of deformation the interface underwent.
88


5.
Conclusions and Recommendations
The level set method's advantages far outweigh the disadvantages. Every
simulation that was conducted had favorable results. It has an immense ability of
handling the breaking and merging of interfaces. The oscillating bubble and
Raleigh-Taylor simulations matched the analytical theory proving it's accuracy in the
modeling of surface tension forces. The buoyancy driven flow simulations correctly
represented the deformation of the interface and steady state bubble Reynold's
number. These simulations confirmed the level set method's ability to match physical
buoyant flow regimes.
From the simulations that were conducted the most apparent area of
improvement for the level set method is retaining bubble area. The buoyancy driven
bubbles lost appreciable area as the solutions evolved. Some of this can be attributed
to the first order reinitialization scheme that was used. Second order [3] and third
order [6] schemes have been used in the reinitilialization step with increasing
success. However, nothing in the level set method considered here acts to directly
enforce the preservation of area. A new algorithm presented by [15] addresses this
weakness by reformulating the reinitialization iteration based on preserving volume.
This procedure is formulated around the theory that since the interface shouldn't
move during the reinitialization the volume should not change either. They
conducted a drop oscillation simulation (n = 4) with less than a 1.3% fluctuation in
mass as the interface evolved.
Another explanation for the loss of bubble area is due to the grid resolution
used in the simulations. In the skirted bubble simulation the skirts broke off which is
most likely a nonphysical result due to low grid resolution. In the dropping bubble
simulation some of the gas was trapped in the liquid and eventually disappeared. It
89


is evident from these simulations that when the geometry of the bubble is depicted
on only a few grid cells the solution becomes questionable.
The level set is a powerful method for solving two-fluid problems. The
success of the method is contingent upon maintaining the level set a distance
function near the interface. If the interface is maintained at a constant and uniform
thickness as the solution evolves the calculation of surface tension is greatly
simplified. If the level set is not a correct distance function the calculated surface
tension force will be inaccurate. Also, the fluid properties will not be distributed
correctly across the interface. An inaccurate distance function also leads to the loss
of bubble area. The greatest challenge faced by the level set method is how to
efficiently and accurately maintain an invariant interface thickness as the solution
evolves.
90


Bibliography
1. W. Shyy, H. S. Udaykumar, M. M. Rao, R. W. Smith, Computational Fluid
Dynamics with Moving Boundaries, Taylor and Francis, Bristol, Pa., 1996.
2. S. Osher, J. A. Sethian, Journal of Computational Physics, 79(1), 12 (1988).
3. M. Sussman, P. Smereka, S. Osher, A Level Set Approach for Computing
Solutions to Incompressible Two-Phase Flow, Journal of Computational Physics,
114, 146-159(1994).
4. J. C. Slattery, Interfacial Transport Phenomena, Springer-Verlag, New York,
1990.
5. R. Aris, Vectors, Tensors, and the Basic Equations of Fluid Mechanics, Dover,
New York, 1989.
6. Y. C. Chang, T. Y. Hou, B. Merriman, S. Osher, A Level Set Formulation of
Eularian Interface Capturing Methods for Incompressible Fluid Flows, Journal
of Computational Physics, 124, 449-464 (1996).
7. J. D. Anderson, Computational Fluid Dynamics The Basics with Applications,
Mc-Graw Hill, New York, 1995.
8. Y. Saad, Iterative Methods for Sparse Linear Systems, PWS Publishing
Company, Boston, Ma., 1996.
9. B. P. Leonard, A Survey of Finite Differences of Opinion on Numerical Muddling
of the Incomprehensible Defective Confusion Equation, Presented at ASME
winter annual meeting, New York, New York, Dec 2-7, 1979.
10. C. Peskin, Journal of Computational Physics, 25, 220 (1977).
11. H. Lamb, Hydrodynamics, 6th ed., Cambridge at the University Press, 1932.
12. S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Dover, New
York, 1961.
91