A computational approach for two-dimensional, gas-liquid flows with surface tension and viscosity

Material Information

A computational approach for two-dimensional, gas-liquid flows with surface tension and viscosity
Khunatorn, Yottana
Publication Date:
Physical Description:
viii, 74 leaves : illustrations ; 29 cm

Thesis/Dissertation Information

Master's ( Master of Science)
Degree Grantor:
University of Colorado Denver
Degree Divisions:
Department of Mechanical Engineering, CU Denver
Degree Disciplines:
Mechanical Engineering
Committee Chair:
Welch, Sam W.
Committee Co-Chair:
Trapp, John A.
Committee Members:
Clohessy, William H.


Subjects / Keywords:
Surface tension ( lcsh )
Two-phase flow ( lcsh )
Surface tension ( fast )
Two-phase flow ( fast )
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )


Includes bibliographical references (leaves 73-74).
General Note:
Submitted in partial fulfillment of the requirements for the degree, Master of Science, Department of Mechanical Engineering.
Statement of Responsibility:
by Yottana Khunatorn.

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:
37757065 ( OCLC )
LD1190.E55 1996m .K48 ( lcc )

Full Text
Yottana Khunatom
B.Eng., Chiang Mai University, 1992
A thesis submitted to the Faculty of
the Graduate School of the
University of Colorado at Denver
in partial fulfillment of the
requirements for the degree of
Master of Science
Department of Mechanical Engineering

This thesis for the Master of Science
degree by
Yottana Khunatom
has been approved
William H. Clohessy
V3 U

Khunatom, Yottana (M.S., Mechanical Engineering)
A Computational Approach for Two-dimensional, Gas-Liquid Flows with
Surface Tension and Viscosity
Thesis directed by Assistant Professor Sam W. Welch
Two-phase fluid flows including surface tension effects are studied
using a numerical method. An existing finite difference method with
staggered-grid and interface tracking are used as our basic numerical method.
The surface tension force is expressed in terms of the curvature of the
interface by considering the surface tension as a force that acts on the
interface. The curvature is found by geometrical method. The discrete form of
the surface tension force is found by interpolation to the computational cell.
The viscous terms are calculated with consideration of the discontinuity in
viscosity across the interface. Simulations of two-dimensional buoyancy-
driven bubbles are presented as a numerical result with surface tension and
viscosity. In addition, an oscillating bubble and the Rayleigh-Taylor
instability, which are inviscid problems, are attempted to test the applicability
of our numerical approach. The numerical results from the buoyancy-driven
bubbles concur very well with the literature. The Rayleigh-Taylor instability
simulation gave reasonable results in both unstable and stable wave number
regions. The results from oscillating bubble are accurate only for the first half
cycle and tend to deteriorate afterward. This, we believe, is due to the
numerical dissipation attendant with the donoring convection scheme. This
work not only shows the accuracy of our approach, but also illustrates
interesting fluid mechanical behavior in the simulations presented.

This abstract accurately represents the content of candidates thesis. I
recommend its publication.
Sam W. Weld

List of
1. Introduction.......................................................1
1.1 Motivation.....................................................1
1.2 Scope of Study.................................................2
2. Numerical Method...................................................4
2.1 The Governing Equations........................................4
2.2 The Numerical Procedure........................................6
2.3 The Finite Difference Method...................................6
2.4 The Gauss-Seidel Iteration.....................................9
2.5 Tracking the Interface........................................12
2.6 Density Calculation.......................................... 16
2.7 Viscous Stress Calculation....................................17
2.8 Surface Tension and Curvature.................................22
2.8.1 Surface Phenomena..........................................22
2.8.2 The Approximation of Surface Tension Effects...............26
2.8.3 The Calculation of Curvature...............................30
2.8.4 The Numerical Method...................................... 32
3. Numerical Results.................................................34
3.1 Buoyancy Driven Bubbles...................................... 34
3.1.1 Theoretical Background.....................................34
3.1.2 Numerical Result...........................................38
3.2 Oscillating Bubble............................................45
3.2.1 Theoretical Background.....................................45
3.2.2 Numerical Result...........................................49
3.3 Rayleigh-T aylor Instability..................................51
3.3.1 Theoretical Background.....................................51
3.3.2 Numerical Result...........................................59
3.4 Conclusions and Recommendations...............................72
Bibliography ........................................................73

Figure Page
2.1. Location of all parameters on a staggered cell..................7
2.2. Geometry in interface cell.....................................15
2.3. Mass fluxed to the right cell by velocity u....................17
2.4. Void and viscosity index.......................................19
2.5. Interface displaced through distance dQ, with
pressure pi and p2 in fluids 1 and 2 respectively.............23
2.6. The surface tension force at the interface on each side
of the control volume..........................................26
2.7. Geometry of curvature calculation..............................31
3.1 Graces graphical correlation..................................37
3.2 Bubble shape, velocity and pressure field.
Buoyancy-driven bubble case 1..................................40
3.3 Bubble shape, velocity and pressure field.
Buoyancy-driven bubble case 2..................................41
3.4 Bubble shape, velocity and pressure field.
Buoyancy-driven bubble case 3..................................42
3.5 Bubble shape, velocity and pressure field.
Buoyancy-driven bubble case 4..................................43
3.6 Bubble shape, velocity and pressure field.
Buoyancy-driven bubble case 5..................................44

3.7 Oscillating bubble at start, half period, and full period.....50
3.8 Interface shape, velocity and pressure field
Rayleigh-Taylor instability case 1 at time 0.15 sec...........61
3.9 Interface shape, velocity and pressure field
Rayleigh-Taylor instability case 2 at time 0.15 sec...........62
3.10 Interface shape, velocity and pressure field
Rayleigh-Taylor instability case 3 at time 0.15 sec...........63
3.11 Interface shape, velocity and pressure field
Rayleigh-Taylor instability case 4 at time 0.15 sec...........64
3.12 Interface shape, velocity and pressure field
Rayleigh-Taylor instability case 1 at time 0.25 sec...........65
3.13 Interface shape, velocity and pressure field
Rayleigh-Taylor instability case 2 at time 0.25 sec...........66
3.14 Interface shape, velocity and pressure field
Rayleigh-Taylor instability case 3 at time 0.25 sec...........67
3.15 Interface shape, velocity and pressure field
Rayleigh-Taylor instability case 4 at time 0.25 sec...........68
3.16 Interface shape, velocity field
Rayleigh-Taylor instability case 1 at time 0.45 sec...........69
3.17 Interface shape, velocity field
Rayleigh-Taylor instability case 2 at time 0.45 sec...........70
3.18 Interface shape, velocity field
Rayleigh-Taylor instability case 3 at time 0.45 sec...........71

I would like to express my sincerest thanks Sam W. Welch for the
thesis topic and his guidance. I would also like to thank to my family for their
love and encouragement. Finally, I would also like to thank Piradee Sograi
and all my friends for their support and encouragement.

1. Introduction
1.1 Motivation
Fluid interfacial motion induced by surface tension plays a
fundamental role in many natural and industrial phenomena. Studies of
capillarity, low-gravity fluid flow, hydrodynamic stability, surfactant behavior,
cavitation, and droplet dynamics in clouds and in fuel sprays used in internal
combustion engines are all examples of phenomenon associated with surface
Surface tension is caused by the following mechanism. The molecules
of a substance usually attract one another. These attraction forces are
different for different substances. Hence at the interface surface that separates
the two media, there is a surface force because of the unsymmetrical
attraction force.
The motion of fluids with interfaces is a moving boundary problem.
Generally, moving boundary problems are complicated and difficult to solve
analytically. These problems not only deal with the transfer of mass,
momentum and heat, but they must also deal with the development of the
interfacial shape and dynamics. Due to the complexity and nonlinearity of
these problems, there are a number of approaches using analytical,

experimental, and numerical methods to understand the physics of .these
flows. This thesis describes a numerical method with which to study a class of
these flows.
1.2 Scope of Study
In this work, we focus on using numerical models in studying the
effect of surface tension on two-phase fluid flows in two dimensions. The
calculations are done under the assumption of an isothermal system. Previous
to this work, a finite difference method with a staggered-grid (MAC type
grid) was developed as the basic numerical scheme. An interface tracking
method based on the Volume-of-Fluid (VOF) method was implemented [1].
In this work we have added numerical models for surface tension and
viscosity for gas-liquid flows. Surface tension is treated as a surface force that
acts on the interface. In our method it is necessary to calculate the curvature
of the interface. This curvature calculation is accomplished by fitting the
osculating circle to the interface at the point where the curvature is needed.
Then the discrete form of the surface tension force is obtained by simple
interpolation onto the computational grid. The sample simulations presented
in this work are

The buoyancy driven bubble
Oscillating bubble
Rayleigh-Taylor instability
The first of these simulations is for viscous fluids. Therefore, a simple
approach to treating the spatial dependent of viscosity is implemented. The
second and third simulations consider inviscid fluids. All three simulations
have the surface tension as an important parameter that affects the physics.

2. Numerical Method
In this chapter we present a numerical method which may be used to
study two-phase flows in which viscosity and surface tension are important.
The method is based on the volume of fluid (VOF) method that is widely used
for treating free surface flows. The existing interface tracking method of
Mortensen, Trapp, and Welch [1] is extended to include surface tension and
discontinuous viscosity.
First, we present the governing equations in differential form including
a surface tension force. Then we review the pre-existing flow solver and
interface tracking algorithm. Lastly, we present the addition of the viscous
stress terms and the surface tension force.
2.1 The Governing Equations
In this work, we consider two dimensional, transient, incompressible,
viscous, and immiscible fluids. The governing equations of mass and
momentum conservation with the surface tension force added are
Conservation of Mass
du dv
dx + dy

Conservation of Momentum
da du da
-----1- u------1- v
dt dx dy
1 dP 1 d da\ Id
p dx p dx\ dx) p dy
da chi'
--- + --
dy dx)
dv dv dv
at + + v 5^
1 8P 1 8
P dy p dy
( dv) 1 8 f du Sv'j
2P V dy) + P dy ^ dy dx)
+ 9y
where w and v are the components of the velocity in x and y direction
respectively. The density is denoted by p, and P refers to the pressure. The
body force components in each direction are gx and gy. The second and the
third terms on the right hand side of the momentum equations are the viscous
stress terms. Note that the viscosity is a function of position. Note also that
the surface tension force components (Fa and Fsy ) are to be modeled as
volume forces [2],

2.2 The Numerical Procedure
The basic finite difference method is a variant of the traditional MAC
method of Harlow and Welch [3], The velocity components and pressure are
defined on an Eulerian, staggered grid. At each time step, finite difference
approximations of the governing equations are constructed. The resulting
linear system of equations is solved by employing a Gauss-Siedel iteration. An
interface tracking procedure has been added and this thesis describes the
addition of surface tension and discontinuous viscosity. A description of the
numerical method is as follows:
Approximate new time velocities are determined by taking an
explicit step in the momentum equation. The surface tension force
and the viscous stresses are treated explicitly during this step.
Using Gauss-Seidel iterations on the mass conservation and
momentum equations solve for the new time pressure and the
exact new time velocity.
Determine the new time density from the mass flux which flows
into or out of each cell with the new time velocity.
Track the new interface by the subgrid counting procedure.
2.3 The Finite Difference Method
As described above, the finite difference discretizations are based on
an Eulerian mesh. A staggered (MAC) grid is employed. In the staggered
grid, the velocity components are calculated for the points that lie on the faces

of the cell, while all other variables are calculated at the cell center. This is
shown in Fig. 2.1.
Location of all parameters on a stagger cell
Figure 2.1
The finite difference approximations of the governing equations are
- Uj-


= -(fuxx) (FUXY) + gx -
Pi+l,j Pi,j

+ (fsx).'. + (viscxx + viscxy)
Paver Pa

v"+1 V. .
= -(FUYX) (FUYY) + gy -
Pj,j+i Pi,j
1 1
+ ----(fsy) + ----(VISCYX + VISCYY)
Pavg Pavg
+ U"w) + - ui.j)
+ u) + )_
Fray = + "*0+ h*i.i|Kj ui.j)
+ i.l) + W|"lH i.i)
2Ax, L *3+s
+ + kj.!^ V*Ul)
- vul)
2A>V *
+ rM*>) + kj^lKj -
+ ^.j) + V)

The cells are numbered by indices i and j. These indices indicate cell center
position along the horizontal and vertical directions, respectively. Cell
boundary positions are labeled with half-integer values of the indices. The
dimensions of the rectangular cells are Ax and Ay. In contrast to the other
variables, which use the same index as the cell numbering, the velocity
components, u and v, are indexed by the same index number but at the
boundary of the cell. The indexing scheme is shown in Fig. 2.1.
In addition to the space index subscripts, we use a superscript to
number the time cycle. For example, u?*1 designates the horizontal velocity
at the time t = (n+J)At, in which At is the time increment per cycle. When
the cycle number superscript is omitted, it is assumed that its value is n,
corresponding to the time t = (n)At. The average density, pavg, is used in the
difference equations due to the presence of two-phase cells. The calculation
of this average density is discussed in the following section. Because of
discontinuity of the viscosity across the interface, our differential
approximation of the viscous terms requires special treatment. This will be
discussed in detail in a later section.
2.4 The Gauss-Seidel Iteration
The Gauss-Seidel iteration is a simple but efficient point-iterative
procedure for a large system of equations. Convergence is related to the
condition of diagonal dominance of the matrix of coefficients [4], The
procedure for a general system of algebraic equations is

Make an initial guess for all unknowns.
Solve each equation for the unknown whose coefficient is on the
diagonal, using guessed values initially and the most recently
computed values thereafter for the other unknowns in each
Repeat iteratively the solution of the equations in this manner until
changes in the unknowns become small.
The exact new time velocities must satisfy the momentum equations
with the new time pressures. The discrete momentum equations in the x-
direction, Eq.(2.5), may be written
u;1 = ultj At(FUXX) At(FUXY) + gx -
pn+l ____ pn+l
ci+1.3 *1.1
+ (Fsx). + (VISCXX + VISCXY)
Introducing the pressure increment, AP, we have
Tpn+1 pn+ll _
[*1+1.1 *i,1 J ~
i + 1
Pavg^1 + i
Pav^X. ,
[Pi+i,3 Pi,j]
[APiVi Api"j+1]
combining Eq. (2.9) and (2.10) we obtain

Un+1 = uexp
i/J hi
At [a- AP"1]

u*xp = Uifj At(FUXX) At(FUXY) + gx -
Ax. _
+ (psx). + (VTSCXX + VISCXY)
Using the same procedure as above in the y-direction, we obtain
vn+1 = vexp -

[ap,:;.1, ap^;1]
~ At(FUYX) At(FUYY) + gy -
P avg

Ay. j
J j+I
+ (psy). (VJSCifX + viscyy)
^ ''"10 p
And the continuity equation is
^ + ytf Ax,
= 0

Substituting Eq.(2.11) and Eq.(2.13) into Eq.(2.15) we obtain
<7 - ~
In the procedure above, the initial guess of the new time velocities are
found by solving the momentum equations with old time pressure values.
These velocities are denoted by uexp in Eq.(2.12) and ve!tp in Eq.(2.14). The
Gauss-Seidel iteration is used to solve Eq. (2.16) for the pressure increments.
This pressure increment then used in Eq.(2.12) and (2.14) to update uexp and
vexp, respectively. The iteration is repeated until the pressure increments
converge to a sufficiently small value. After convergence we have the new
time velocities and pressures.
2.5 Tracking the Interface
This section describes the simple subgrid counting procedure
presented by Mortensen, Trapp, and Welch [1], The method uses the void

fraction to indicate the existence of interface cells. The cell is pure liquid
when the void fraction is 1 and pure gas when the void fraction is 0. Two-
phases cells exist when the void fraction is between 0 and 1. The void
fraction, a, is defined in term of the densities of the two-phases as
P~ P*
Pi ~ Pg
where p is the cell density and pg and pi are the density of the gas and
the liquid phases, respectively. Assuming the void fraction is continuous, the
gradient of the void fraction is a vector perpendicular to the surfaces of
constant void fraction pointing in the direction of maximum void fraction
increase. The normal vector at interface in any two-phase cell pointing into
the liquid region can be calculated by using an averaged second order gradient
of void fraction
a* = I-
lV A a
+ i A
AJKAxJ^ V2AAx>i(j
+ |i¥^
av =
^ A
+ ,i

The vector is normalized as follows

Since the unit normal vector at the interface in a cell and the volume
place the interface at the correct location within the cell. These geometric
calculations are done using a simple subgrid counting scheme based on a
subgrid mesh in the cell. At the interface cell, the subgrid is constructed. If the
dimension on each side of the cell is Ac, the dimension of the subgrid cells is
Axsub=Ax-H7sub, where nsiib is the number of subgrid cells in each direction.
As shown in Fig. 2.2 the possible location of the interface is along the
direction of the normal vector and can be found by comparing the true void
fraction with the void fraction calculated from this counting. The true void
fraction calculation is discussed in section 2.6.
of the cell occupied by gas are known, geometric calculations must be done to

Figure 2.2
The distance along the unit normal vector measured from the cell center is
defined as S. Dividing the number of subgrid cells on the liquid side of the
plane by the total number of subgrid cells in each cell is approximately the
void fraction as a function of the interface location S. S can be solved for
consistent with the known void fraction by mapping S into a function of cell
void fraction and applying the Secant iteration method with this function. As
shown in Fig. 2.2 the vector f is a vector from the intersection of the
normal vector and the interface to any subgrid cell. The liquid cell is counted
if the inner product f h is positive and the cell is counted as vapor cell if
the inner product is negative. If a subgrid cell center is within the distance
AxSUb+2 of the interface, the unit counter for this subgrid cell is weighted

linearly between zero and one. It is zero if the subgrid cell center is exactly a
distance Axmb^2 from the interface on the vapor side. It is one if the subgrid
cell center is exactly a distance 4x^+2 from the interface on the liquid side.
2.6 Density Calculation
The new time density of each cell is calculated by the mass flux across
the cell consistent with the new time velocity. The mass flux is determined by
calculating the volume of fluid that is fluxed out of the cell during each time
step. Shown in Fig. 2.3 is an example mass flux calculation in the x-
direction. The positive x-direction velocity will flux the mass on the right face
of the cell within the distance uAt of the right face of the cell. The mass flux is
calculated as the density in this cell multiplied by the flux volume. In a two-
phase cell the mass flux is calculated using the same method but the subgrid
counting procedure is employed to find the volume of gas and the volume of
liquid. The gas mass flux is then the gas density multiplied by the gas volume
and the liquid mass flux is the liquid density multiplied by the liquid volume.

Mass in liquid and vapor region fluxed to the right cell by velocity u
Figure 2.3
2.7 Viscous Stress Calculation
The Navier-Stokes equations with constant viscosity and density are
applicable to each single phase. Since these variables are both discontinuous
across the interface, we may expect either excessive numerical diffusion or
problems with oscillations around the interface if no special treatment is used
at the interface. As described in an earlier section, the density at the interface
is calculated by the subgrid counting method. The momentum equations,
Eq.(2.2) and Eq.(2.3), can be expressed as

da du du 1 (do dx
+ u + v = - +
dt dx dy p \ dx dy )
+ ~ Fsx + &x (2-21),
dv dv dv 1
+ u + v =
dt dx dy p
V dx
We are considering isotropic, incompressible, Newtonian fluids, hence the
stresses are
ox = -P + 2p
_ dv
a = -P + 2p
* xy ~ Tyx ~ M
' dv_ du_
\dx dy>
In constructing finite difference approximations for the viscous stress
terms in Eq.(2.21) and (2.22), we must take into account the spatial
distribution of viscosity. Consider the x-direction momentum control volume
shown in Fig.2.4,

Hi-ij+i Hij+i Vj j Ui+lj+1
l ^ a Otij+1/2 Jij -1/2 j ;J a V:; 1 +1/2J +lj
Ui-ij- IJ 1 Ctij-1/2 Uij-1 Ui+Ij-
Void and viscosity index
Figure 2.4
the shear stress, txy, acts on the top and bottom faces, therefore we need the
viscosity at the top and bottom faces. Similarly the viscous stress, Txx, acts on
the left and right faces hence we need the viscosity at those faces. The finite
difference approximations for the viscous terms in the momentum equations


f f \ ( Y\
l Ayi+1 J Pi.l l AVj JJ
Where the indices on the viscosity indicate its special location, as shown in
Fig. 2.4. The viscosity at the cell center is calculated by exact void weighting,
the viscosity at the right and left junctions in the x-direction momentum cells
is calculated by exact void weighting. At the top and bottom faces of the x-
direction momentum cells, the viscosity is calculated by average void
weighting. A similar procedure is used for the y-direction momentum cells.
This may stated mathematically as,

htj = + (l ~ aitj)fi9
1 i+i.:K (2.27),
1 ai,j+l)Mg (2.28),
i)/' + - Km),)"* (2.29),
>)/J + i1 - (2.30),
i).M + i1 - K=.hLK (2.31),
i)f+ t1 - (2.32),
[a i i) = fa. a + a. i)
\ 1+i'J+\)y 2 \ 1,3+*


. + a
2.8 Surface Tension and Curvature
Flows involving surface tension are the main focus of this work. The
including this phenomena into the governing equations.
2.8.1 Surface Phenomena
Following [5], we seek to define the surface tension in Laplaces
formula, which gives the surface pressure, as a cause of the pressure
difference of the jump across the interface when the interface surface is
curved. Consider a segment of an interface with area dA cL^dl,,
displaced in the normal direction by an infinitesimal virtual displacement, 5£.
dJ\ and d_L, are the lengths on each side of the element.
following section will describe briefly surface phenomena and our method of

Interface displaced through distance dwith
pressure pi and P2 in fluid 1 and 2 respectively.
Figure 2.5
The virtual work done in displacing the area element shown is
J (- Pi + P2)5CdA (2.34),
where pi and p2 are the pressure of fluid on each side of the interface. If we
define 5W is the total virtual work done in displacing the surface we must add
the virtual work done in changing the area of the surface. This work is
proportional to the change of the area of the surface, 8A. Thus the total
virtual work is
5w -J p2)g£dA + odA

where ct is the fluid surface tension coefficient. If the system is in an
equilibrium condition, the virtual work, 8W, is zero. If Ri and R2 are the
principal radii of the curvature of the surface. According to the infinitesimal
size of the normal virtual displacement, the corresponding virtual increments
of length d\ and dL, are (sQR^dl^ and {hC,/R^)dl2, respectively. Then
the virtual surface area increment can be expressed as
( 5Cl r
1 + dJL
l + £
. V
1 +

The total virtual change in area of the interface is
5A = | sc|
Substituting Eq.(2.37) into Eq.(2.35) and setting the total virtual work to
zero we obtain the equilibrium condition
J Sqta p2) a
>dA = 0

Since S£ is arbitrary, the expression in the bracket must be zero:
(Pi Pi)
^ +
Eq.(2.39) is Laplaces formula, which relates the surface tension to the fluid
pressures and the surface curvature. For the interface between two viscous
fluids in motion, Laplaces formula may be generalized to
(Pi ~ P2h = (hj* ~ T2i*K + CT '
1 1^
Since the surface tension coefficient is not always constant. There is another
force tangential to the surface, grada. The momentum balance at the
interface may be shown to be
Note that above equation is the projection of the normal and tangent
coordinate into the Cartesian coordinate. The last term on both the left hand
side and the right hand side are the surface tension terms in the normal and
tangential directions of the interface, respectively. This work considers
isothermal problems, hence the surface tension coefficient is constant. The

last term on the right hand side is zero. Hence, the effect of the surface
tension on the momentum at the interface is only the normal direction.
2.8.2 The Approximation of Surface Tension Effects
This section describes how the surface tension forces is modeled and
applied to the momentum equation. The effect of surface tension is a surface
force that acts on the interface, as shown in Fig. 2.6.
The control volume form of the momentum equation can be written as:
The surface tension force at the interface on each side of
the control volume.
Figure 2.6
where T is the stress tensor and h is the outward directed unit normal vector
to the surface of the control volume. The body force is represented by g and

Fsv is the surface tension force. Considering the interface surface as
infinitesimally thin, the surface tension force is added into the momentum
equation as a line integral. By employing transport theorem [6]
f 'FdV = f
cJt v
+ 'Fdivv
where Y can be any scalar, vector, or tensor, and the divergence theorem:
j^divvdV = ndS (2.44),
the momentum equation, Eq.(2.42), can be rewritten as
Dt = jvdiv^dv + \vP$dV + jc^dL (245a)>
l ^ - divT pg dV = £FsvdD
Consider the right hand side term, which is the surface tension force term.
Viewing Fig 2.6, we calculate

£rsvdi = £otdi
= J crtjdz + J atzdz
+ J atfdx + J atadx
where t is a unit vector in the tangent plane of the surface directed
outward on the boundary. Since we consider two dimensional problems, the
components of above equation in the third dimension vanish because t£
and ta are equal in magnitude but opposite in direction. We have
tz -£j + At, Eq.(2.46)becomes
i^L = J afc + tr)dz
= j* a(At)dz
Substituting this into Eq.(2.42) we have
£ (pv) divf pg dV = J cr(At^jdz
Our discrete approximation to the control volume equation is
(pv) divf pg AxAyAz = cf(a£)Az (2-49).
Dividing by AxAyAz we get

D x
(yov) = divT + pg +
DtK '
From the definition of curvature, ic = , where t is the tangent vector
of the curve, and dS is an increment along the length of the curve.
Substituting a discrete approximation to the curvature into Eq.(2.50) we have
D , \
(pv) = divT + pg +
The last term on the right hand side is the surface tension term which is added
to the momentum equation in this work. Note that we explicitly show the
direction of the force is in the direction normal to the interface.
Surface Tension Force = -L- (2.52)
The x and y component of the surface tension force was included in Eq.(2.2)
and Eq. (2.3), where Fsx and FSy are the components of the surface tension
force in x and y directions, respectively.
We consider flows with constant surface tension. The arc length, AS,
and the curvature, k, in Eq.(2.52) are found by a simple geometric method.

These parameters influence the accuracy of the surface tension force. Of the
two parameters, the curvature is the most difficult to calculate. The next
section will describe the curvature calculation.
2.8.3 The Calculation of Curvature
This section discusses the calculation of the curvature that is
employed in this work. Since by definition, the curvature vector is the change
in the unit tangent vector along a curve with respect to changes in the arc
length S, along the curve,
k = \ (2.54),
where R is the radius of curvature.
The concept our method uses is to solve for the radius of curvature
and then the reciprocal of this radius is the curvature, as expressed in Eq.
(2.54). Consider Fig 2.7, (x0, yo), (x/, yj), (x2, y2), are points on the interface
which define a curve that we want to find the curvature of.

Geometry of curvature calculation.
Figure 2.7
Let Si be a straight line connecting points (x0, yo) and (xj, yf), and S2 be a
straight line connecting points (xi, yj) and (x2, y2). If Li and L2 are the
perpendicular bisectors to the line Sj and S2 respectively, the point of
intersection of Lj and L2, (xc, yc), is the center of the curvature of the circle
that passes through these three points. In effect, we approximate small
sections of the interface as circular arcs. These circular arcs are known in the
literature as osculating circles [7],

Each cell that contains an interface may have two or more neighboring
cells that contain interfaces. We have found that selecting neighboring cells
which have void fractions closest to 0.5 give us a more accurate curvature
calculation. These cells have a longer length of interface. Generally, the
longer the interface length, the more accurate the unit normal vector is. The
unit normal vector is important in applying the above geometrical method for
calculate the curvature. The curvature vector can be calculated by substituting
this curvature into Eq. (2.54).
2.8.4 The Numerical Method
As described in section 2.3, a staggered grid is used in our finite
difference approximation. The surface tension calculations discussed
previously provide a cell centered quantity. We interpolate this cell centered
quantity to momentum cells by simple averaging. Hence the surface tension
terms in Eq. (2.2) and Eq. (2.3) are expressed as
= \ ((p4 + (Oj (2.55a),
for the x-direction momentum, and
^ = | (W.+ teA)

for the y-direction momentum. Here (f^x) and are the components of
the surface tension forces in the x direction of the left hand side and the right
hand side of the x-momentum control volume, respectively. Similarly, (fsy) ,
and (^,y j are the components of the surface tension force in the y direction
of the top side and the bottom side of the y-momentum control volume,
respectively. According to Eq. (2.5) and Eq. (2.6), the numerical expression
of Eq. (2.55) are
Xf-Lu, + (
where (normx) and (normy) are unit normal vector components of the
interface in x and y directions of two-phase cells.

3. Numerical results
This chapter presents some background theory and computational
results from the three simulations attempted. These simulations were chosen
to test key features of the method. The simulation of the buoyancy driven
bubbles is a case where the fluids can have arbitrary density difference, and
surface tension and viscosity are important parameters that effect the
characteristic of the bubbles. A number of bubbles with various Eo (Eotvos
number) and M (Morton number) were simulated and the results compare
favorably with results in the literature. The second simulation is an oscillating
bubble. This inviscid bubble oscillation is driven by surface tension without
gravity. The last simulations are of Rayleigh-Taylor instabilities of two
inviscid fluids of different density. Details of all simulations will be discussed
in the following sections.
3.1 Buoyancy Driven Bubbles
3.1.1 Theoretical Background
Bubbles can be considered as fluid particles which occur in many
natural phenomena ( e.g. boiling) and industrial processes (e.g. fermentation).
The behavior of bubbles differ from that of solid particles in two main

the bubble is deformed by hydrodynamics forces, and
an internal flow will result from momentum transfer across the
gas-liquid interface.
Due to buoyancy forces, unconstrained gas bubbles rise in liquid. The
motion of the liquid induces a change in pressure. The inertia forces, pressure,
and surface tension lead to deformation of the bubbles. In spite of a number
of previous works, there is no satisfactory theory except for the cases of very
small deformation at either high or low Reynolds number [10]. However, for
engineering purposes, a number of experimental works have been published in
terms of dimensionless groups and empirical correlations. In this case of
bubbles rising freely in an infinite media, as shown in Fig. 3.1, Grace[8] has
summarized available experimental data and presented it as a graphical
correlation in term of three dimensionless parameters;
Eotvos number
Morton number
_ gpAp
_2 _3
and Reynolds number

Pde U
Where Ap is the difference between the density of surrounding fluid
and bubble fluid. The dimension of the bubble is given as an equivalent
diameter, de and the terminal velocity of the bubble is U. The subscript T
means the quantity is the property of the surrounding fluid. As expressed in
Equation (3.1)-(3.3), the Eotvos number represent the ratio of the buoyancy
force to the surface tension. As far as we know, there is no simple
interpretation for the Morton number. The Morton number does include
surface tension, inertia, and viscosity. The Reynolds number is the bubble
Reynolds number.
The effects of these dimensionless parameters on the deformation of
the bubble is represented by Graces graphic correlation At a constant Morton
number, increasing bubble diameter or decreasing surface increased the
Eotvos number. Therefore, the shape of bubble changed from spherical to
ellipsoidal or wobbling and finally to spherical cap or dimpled. At a constant
Eotvos number, increasing the Morton number causes the bubble shape to
change from wobbling or spherical cap to ellipsoidal or dimpled and finally to
spherical. The Morton number can be increased by increasing surrounding
fluid viscosity, by decreasing surrounding fluid density, or by decreasing
surface tension.
Note that the Reynolds number increases with decreasing Morton
number at constant Eotvos number. The Reynolds number increases with an
increase of Eotvos number at constant Morton number. The effect of the
Reynolds number on the bubble motion will be discussed later.

Figure 3.1 Graces graphical correlation

According to the previous studies about the effect of Eotvos number
and Morton number on the bubble motion, Tomiyama, Zun, Sou, and
Sakaguchi [9], they have been presented that the fluctuation motion of the
bubble can be terminated by increasing the Morton number or by decreasing
the Eotvos number. As mention above, the Reynolds number decrease when
the Morton number increase.
3.1.2 Numerical Results
A number of bubbles for various Eo and M are simulated to compare
to results in the literature. The approximate steady state shape of the bubbles
with velocity fields and pressure profiles are shown in Fig. 3.2-3.6. All
simulation are done at 65 by 33 grid resolution except Fig. 3.2 and 3.6 which
are done at a finer grid resolution, 130 by 66. All calculations are done with
the density ratio pj pb =40 and free slip boundaries. The Eotvos
number, Morton number, and the viscosity ratio are noted in the caption of
each figure.
Considering Fig.3.2-3.5, the bubbles in steady state at small Eo ( small
bubble or large surface tension) are close to circular shape. This is because at
small Eo, the surface tension and viscous forces are much more important
than the inertia force. Hence the deformation of the bubble is very small. At
higher Eo number the bubbles obtain a steady state at elliptic shape or nearly
flat bottom and semicircular top. When Eo number becomes large ( large
bubble or small surface tension), the flow separation causes the wake at the
rear of the bubbles to expand. As a result, the bubbles become dimple shapes
or ellipsoidal/spherical cap shapes, as shown in Fig. 3.6.

Considering small Eo and small M bubble (Fig. 3.2), the bubble is in
the ellipsoidal region of Graces chart. Note that the bubble has the top part
slightly flatter than the bottom part which is because the front stagnation
pressure becomes dominant as described in Ryskin and Leals work [10].
Comparing Fig. 3.2 with Fig. 3.3, the pressure change in the liquid around the
bubble is higher in the lower Morton number bubble. The bubble flow in the
less viscous fluid causes a larger change in the pressure field of the outer
We compare our simulations with Graces graphical correlation by
ignoring that Graces chart is based on three-dimensional and infinite
boundary flows. A more precise comparison is made with the two-dimension
results of Tryggvason and Unverdi [11]. For the same Eotvos number,
Morton number, viscosity ratio, and density ratio, our bubble shapes are in

Mo = I0~7 Eo=1.0
Bl = 40 ^-=88
P* Pi
Figure 3.2 Bubble shape, velocity and pressure field. Buoyancy-driven bubble
case 1

Mo =10~ 4 Eo=1.0
= 40 = 493
Figure 3.3 Bubble shape, velocity and pressure field. Buoyancy-driven bubble
case 2.

Pb n*
Figure 3.4 Bubble shape, velocity and pressure field. Buoyancy-driven bubble
case 3.

Mo =10'2 Eo = 10.0
40 = 277
Pi P*
Figure 3.5 Bubble shape, velocity and pressure field. Buoyancy-driven bubble
case 4.

Mo = 100 Eo= 104.0
P;/Pi=40 Ji,/fi6=479
Figure 3.6 Bubble shape, velocity and pressure field. Buoyancy-driven bubble
case 5

3.2 Oscillating Bubble
3.2.1 Theoretical Background
This is a good test simulation as there exists an analytic solution with
which to compare result. We consider a gas bubble surrounded by liquid in an
equilibrium configuration. The flow is invisid and we neglect gravity.
According to these conditions, the surface tension force becomes dominant in
stress balancing. When the bubble is perturbed from its equilibrium shape the
surface tension will tend to restore the bubble to its original shape. An
analysis of the vibration of a circular column of gas in surrounding liquid by
Sir Horace Lamb [12] applies to this flow. The analytic procedure is as
Consider the potential equation
£l + L^L.JL£l = o
dr2 r dr r2 dd2
a solution for the above equation is

We have the conditions that the potential vanishes at the center of the bubble
(r=0), and at infinity (r=o). At the radius of the bubble

It can be shown that the potential inside and outside of the bubble are

*1 =
cos (nff) cos (at + e)

cos (n6) cos (at + s)
respectively. The radius of perturbed bubble can be expressed as
r = a + 4
where a is the mean radius of the bubble. Particles at the interface move with
the interface, then
<%_ _ c%
ft r=a ft
Hence the perturbation displacement is
£ =
-A cos (n&) sin (at + e
The jump of pressure across the interface is
(r -) = <{D

_ i 1 d2r
~ r ~ r2 d02
Substitute this equation into Eq.(3.8) and linearizing we have
£ a a2 l d.Qz j
Consider the Bernoullis equation
a + p + w
dt p 2
for the perturbation
= + (/>'
P = Po + P
p = P0 + P'
4^+ + (po + ^') + (V(^+ ^))2
dt p 2
= C(t)
Linearizing with the reference terms canceling the right hand side, Eq. (3.16)
P_ _
p dt

According to this equation, the perturbation pressure for inside and outside
the bubble are
P;' = -p^A cos (nQ) sin(o)t + e) (3.18a),
P'e = pe(oh cos (n9) sin (cot + s) (3.18b),
Substituting Eq. (3.18), Eq. (3.15), and Eq.(3.13) into Eq. (3.11), we have
P0i P0e ~ (Pi + Pe)>
1 1 n n-3 T
a a2 coa taaj
P01 P0e ~
We obtain the frequency of the oscillation for the various modes expressed as
cr(n3 n)
(Pi + p*y

3.2.2 Numerical Result
A bubble is set up in the computational domain with the gravity and
viscous effects neglected. The calculation is done with 40 by 40 grid
resolution and free slip boundaries. The time step of the calculation is set to
1/100 of the complete period of oscillation. Initially, the bubble is perturbed
from its circular shape to an elliptical shape. We compare the bubbles shape
with the linearized analysis above. The numerical results for the first half
period and one full period are shown in Fig. 3.7. The simulation shows a
good result only for the first half period but deteriorates afterward. This
simulation didnt improve on a finer grid resolution. This problem needs to be

i to a a
Figure 3.7 Oscillating bubble at start, half period, and full period

3.3 Rayleigh-Taylor Instability
3.3.1 Theoretical Background
When the two fluids of different densities are superposed with one
over the other, or accelerated towards each other; the instability of the plane
interface is called the Rayleigh-Taylor instability. This is an interesting
problem because the perturbation of the interface grows exponentially in time.
Another interesting aspect to this problem is that the surface tension will
stabilize small wave length perturbation. This is a useful simulation because
we should be able to observe this behavior in our computational model.
Following is a derivation of the critical wave length.
Consider the two dimensional Navier-Stokes equation for
incompressible, inviscid flow.
Conservation of Mass
du dv_
dx + dy
Momentum Equation
+ pu
+ pv

dv dv
+ PUT~ + Pv ~Z~
dx dy
- + pg
We perturb these equations with
u = u0 + u'
V = v0 + V
P = Po + P'
p = p0 + P'
Where the subscript O means the values of parameters at the reference
configuration, and the reference configuration velocity is the velocity at a
non-disturbed point, u0 = v0 =0. Substitute Eq.(3.24) into Eq.(3.22)-
dP dP
(3.23), and neglect quadratic terms. Since - = 0, and - = p0g,
dx dy
Eq.(3.22) and (3.23) become
Conservation of Mass
du' dv'
--- + ---
dx dy
Momentum Equation
du' dP'
--- + ----
dt dx

P 0
dv' dP'
dt dy
= P'9
We try solution of the form
v' = v(y)eMkx~tot)
p'= p{y)eUkx-^]
P'= P(y)eMkx~eot]
and substitute Eq.(3.27) into Eq.(3.25) and (3.26) to obtain
ikfiM + = 0 dy (3.28),
-icop0u(y) + JcP(y) = 0 (3.29)
dP(y) iap0v(y) + = p(y)g dy (3.30).
Since particles must have constant density
Dp dp dp dp = + u + v = 0 Dt dt dx dy (3.31).
Then employ Eq.(3.31) with Eq.(3.24) to obtain
= 0 dt dy (3.32).

Substitute Eq.(3.32) into Eq.(3.27) to get
-i coply) + vly) = 0
and rearrange as
ico dy
Combining Eq.(3.28) and Eq.(3.29) we get
Ply) =
icop0 dvly)
k2 dy
Eliminate Ply and ply in Eq.(3.30) by employing Eq.(3.35) and Eq.(3.34)
resulting in
. icop cfv -ig dp0 ,
-icop0vly) + f =-------------vly)
k dy co dy
and rearrange to
gkz dpQ k2
eozp0 dy
vly) = 0
The general solution, in each phase, of Eq.(3.37) is
vly) = C,e*y + C2e


In region 1, v = 0 at y > -oo,then
v.iy) = Cxeky + C2e~ky (3.39),
implies C2 = 0 and
^(y) = Cieky (y < 0) (3.40).
For region 2, v = 0 at y > co, then
v2(y) = C3eky + C4e~ky (3.41),
implies C3 = 0 and
v2(y) = C4eky (y > 0) (3.42).
At y = 0, v3 = v2 then
C2 = C4 = C
and Eq.(3.40) and Eq.(3.42) become
v.iy) = Ceky (y < 0) (3.43),
v2(y) = Ce~ky (y > 0) (3.44).
Since at the boundary of the separation between the fluids there is a
discontinuity in the distribution of density, the surface tension effect must be

employed. The effect of surface tension on a fluid interface is included by
employing Laplaces equation.
P, P, = kg a
where k is the curvature of the surface and R is the principal radius of the
curvature. Let the interface displacement be 4(x,t), then the right hand side of
Eq.(3.45) can be expressed as
To include the surface tension effect in the momentum equation, Eq.(3.26),
we must include on the right hand side of Eq.(3.26b) the additional term
s(y 4 where 5(y £(x)) denotes Diracs 8-fimction. Thus the momentum
equation, Eq.(3.26b), becomes
p' g
8{y £(x))
In the above equation, we can relate E, to the velocity at the interface as

£ =
Substitute Eq.(3.50) into Eq.(3.48) and apply with Eq.(3.27) to obtain
v a>Po tfviy)
-aPov(y) +
oo oy
8{y £(*))
As the equation involves 8-functions, we must integrate the equation
across the interface over the point of discontinuity, v is continuous across
the interface. Integrating across the interface
cop0 d2v[y) g dpQ
T ~ O)p0v(y) + v{y)
k oy oo dy
S(y |(x)) dy
where vs is a velocity at the interface, we obtain
aPo 8v{y)£
k2 dy
P*vsy _c + P0(y)^s
c oo
-----vsa = 0

Taking the limit as £ > 0, the second term vanishes and v at y = e is
v8. Eq. (3.53) becomes

Solving for
(Pi ~ Pi)
(P2 + Pi)
k*_ a
Q (p2 + Pi)
According to Eq.(3.55), the left hand side of the equation expresses the
square of the perturbation growth rate, (icu). The arrangement is always
stable if p2 < px. If p2 > px, the arrangement is unstable for all wave
numbers, k, in the range of 0 < k < kc, where
(Pi ~ Pl)7
Further discussion, including the viscous case is found in Chandrasekhar [13].

3.3.2 Numerical Result
We simulate the case of a heavy fluid laid over a light fluid with
downward pointing gravity. The simulation is done with 50 by 50 grid
resolution, (p2/Pi) of 11, surface tension 0.01 N/m, gravity 10N/m2 and free
slip boundaries. The critical wave number, according to Eq.(3.56), is 100.
Four different wave numbers are chosen from the unstable region and the
stable region. Fig. (3.8)-(3.11) are the simulation results at physical time 0.15
sec., and Fig. (3.12)-(3.15) are the simulation results at physical time 0.25
sec.. Note that the computational domain of each wave number are varied to
fit with their wave length. Then the computational domain of the wave
number 200 is four times smaller than the computational domain of the wave
number 50.
Fig. (3.8) and (3.12), show simulation results for the wave number 50.
This wave number is in the unstable region. The figures clearly show a
developing of the unstable wave at the interface. In contrast, the results with
an initial wave number from the stable region show the oscillation of the
interface instead of the growth of the unstable wave, as shown in Fig. (3.12)
and (3.15). According to the size of the computer domain, we can see that the
amplitude of the result of wave number is four time smaller than of wave
number 50.
Considering the results of wave number 90 and 110, the results is very
similar, but the amplitude of the unstable wave at the interface of the wave
number 90 is larger. Fig. (3.16)-(3.18) show the simulation at physical time

0.45 sec. of wave number 50, 90, and 110, respectively. These figures show
the most unstable situation of the simulations for these wave number.
Since the simulations are done with close boundaries, the effect from
the boundaries causes new wave numbers to exist and lead to instability. As
shown in the results of the wave number 110 and 200, which are in the stable
region, the boundary effects cause the unstable wave length to appear for the
simulation with initial wave number 110 and cause an increase in amplitude of
the oscillation for the wave number 200.

k = 50 ncyc = 150 time = 0.15 sec
Figure 3.8 Interface shape, velocity and pressure field. Rayleigh-Taylor
instability case 1 at time 0.15 sec.

k = 90 ncyc = 150 time = 0.15sec
Figure 3.9 Interface shape, velocity and pressure field. Rayleigh-Taylor
instability case 2 at time 0.15 sec.

k = 110 ncyc = 150 lime = 0.15sec
Figure 3.10 Interface shape, velocity and pressure field. Rayleigh-Taylor
instability case 3 at time 0.15 sec.

k = 200 ncyc = 150 time= 0.15sec
Figure 3.11 Interface shape, velocity and pressure field. Rayleigh-Taylor
instability case 4 at time 0.15 sec.

k 50 ncyc = 250 time = 0.25sec
Figure 3.12 Interface shape, velocity and pressure field. Rayleigh-Taylor
instability case 1 at time 0.25 sec.

k = 90 ncyc = 250 time 0.25 sec
Figure 3.13 Interface shape, velocity and pressure field. Rayleigh-Taylor
instability case 2 at time 0.25 sec.

k = 110 ncyc = 250 time = 0.25sec
Figure 3.14 Interface shape, velocity and pressure field. Rayleigh-Taylor
instability case 3 at time 0.25 sec.

k = 200 ncyc = 250 time 0.25 sec
Figure 3.15 Interface shape, velocity and pressure field. Rayleigh-Taylor
instability case 4 at time 0.25 sec.

k =50 ncyc = 450
time = 0.45 sec
Figure 3.16
Interface shape, velocity field. Rayleigh-Taylor instability case
at time 0.45 sec.

k =90 ncyc = 450
time = 0.45 sec
Figure 3.17
Interface shape, velocity field. Rayleigh-Taylor instability case
at time 0.45 sec.

Figure 3.18
£=110 ncyc = 450 time 0.45sec
Interface shape, velocity field. Rayleigh-Taylor instability case 3
at time 0.45 sec.

3.4 Conclusions and Recommendations
For the buoyancy driven bubble, the simulations show good
agreement with the literature. The simulation of dimple, ellipsoidal/spherical
cap bubbles should be attempted at finer grid resolution. This is because the
expansion of the wake behind bubble generate bubble tail which have high
The simulation of an oscillating bubble give us an accurate result in
the early part of the simulation but not accuracy for the whole simulation.
However, when we compare with the analytical solution, this simulation gives
some reasonable results. We believe that the inaccurate results were caused
by the numerical dissipation of the donoring covection scheme. This problem
needs to be improved.
The Rayleigh-Taylor instability simulation provides us some
information on the behavior of the flow in the unstable region. The simulation
in the stable region indicates the accuracy of the numerical method. However,
the boundary condition needs to be improved to simulate infinite boundaries.
This is because the closed boundaries cause the creation of new wave
numbers that lead to instability.

1. Glen Mortensen, John A Trapp, and Sam Welch. Interface Tracking in
Two-Phase Flow Simulations Using a Simple Subgrid Counting
Procedure, Numerical Method in Multiphase Flows 185, ASME, 1994.
2. J. U. Brackbill, D.B. Kothe, and C. Zemach. A Continuum Method for
Modeling Surface Tension, J. of Comp. Phys. 100, 335 (1992).
3. Francis H. Harlow and J. Eddie Welch, Numerical Calculation of Time-
Dependent Viscous Incompressible Flow of Fluid with Free Surface,
Phys. Fluids 8 (1965), 2182.
4. Joe D. Hoffman, Numerical Methods for Engineers and Scientists,
McGraw-Hill Book Company, New York, 1992.
5. L. D. Landau, E. M. Lifshitz, Fluid Mechanics, Pergmon Press, London,
6. John C. Slattery, Momentum, Energy, and Mass Transfer in Continua,
2nd ed., Robert E. Krieger Publishing Company, Huntington, New York,
7. Martin M. Lipschutz. Theory and Problems of Differential Geometry,
Schaums Outline Series, McGraw-Hill Book Company, New York,
8. R. Clift, J. R. Grace, and M. E. Weber, Bubbles, Drops and Particles,
Academic Press., New York London, 1978.
9. Akio Tomiyama, Akira Sou, and Iztok Zun, Numerical Analyses of A
Two-Dimensional Bubble with VOFMethod,
10. G. Ryskin, and L. G. Leal, Numerical solution of free-boundary problems
in fluid mechanics. Part 2., J. Fluid Mech. 148, 19 (1984).
11. Gretar Tryggvason, Salih Ozen Unverdi, .4 front-Tracking Methodfor
Viscous, Incompressible, Multi-fluid Flows, J. Comp. Phys. 100, 25

12. Sir Horace Lamb, Hydrodynamics, 6th ed., Dover Publications, New
York, 1945.
13. S. Chandrasekhar, Hydrodynamic andHydromagnetic Stability, Dover
Publications, Inc., New York 1961.