Citation
Direct numerical simulation of a translating vapor bubble with phase change

Material Information

Title:
Direct numerical simulation of a translating vapor bubble with phase change
Creator:
Ye, Tao
Publication Date:
Language:
English
Physical Description:
x, 151 leaves : ill. ; 29 cm.

Subjects

Subjects / Keywords:
Boundary conditions ( jstor )
Cartesianism ( jstor )
Curvature ( jstor )
Fairings ( jstor )
Heat transfer ( jstor )
Liquids ( jstor )
Simulations ( jstor )
Vapor phases ( jstor )
Vapors ( jstor )
Velocity ( jstor )
Dissertations, Academic -- Mechanical Engineering -- UF ( lcsh )
Mechanical Engineering thesis, Ph. D ( lcsh )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 2001.
Bibliography:
Includes bibliographical references (leaves 145-150).
General Note:
Printout.
General Note:
Vita.
Statement of Responsibility:
by Tao Ye.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
027030110 ( ALEPH )
47365878 ( OCLC )

Downloads

This item has the following downloads:


Full Text










DIRECT NUMERICAL SIMULATION OF A TRANSLATING VAPOR BUBBLE
WITH PHASE CHANGE











By

TAO YE











A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY


UNIVERSITY OF FLORIDA


2001





















To My Family














ACKNOWLEDGEMENTS


I have been fortunate to have the privilege of working with Professor Jacob

Chung and Professor Wei Shyy. Professor Chung, through his constant support and

encouragement, has made it possible for me to finish this research. Many insightful

discussions with Professor Shyy have made my Ph.D. study a very rewarding learning

experience. I am also grateful for their tremendous guidance in my personal growth

and development.

I sincerely thank Professor William Tiederman, Professor James Klausner and

Professor Ren-Wei Mei, for serving as members of dissertation committee and gra-

ciously offering thoughtful advice to help improve this work.

I would also like to acknowledge the collaborations with Dr. Udaykumar and

Dr. Mittal in conducting the work of Chapter 3 of this dissertation.

Most of all I would like to express my deepest appreciation to my family for their

patience, love and support throughout all my professional endeavors.

The financial support was provided by the Department of Mechanical Engineer-

ing, Andrew H. Hines, Jr./Florida Progress Endowment Fund and a grant from Eglin

Air Force Base.














TABLE OF CONTENTS




ACKNOWLEDGMENTS ............................. iii

NOMENCLATURE ................................ vi

ABSTRACT .................................... ix

CHAPTERS

1 INTRODUCTION .............................. 1

1.1 Problem Statement ......................... 1
1.1.1 Research Motivation ........................ 1
1.1.2 Physics of a Boiling Vapor Bubble in Two-Phase Flows 2
1.1.3 Schematic of Computational Model ............ 3
1.1.4 Assumptions ......................... 4
1.1.5 Governing Equations ........................ 5
1.1.6 Boundary and Interfacial Conditions ............... 6
1.1.7 Nondimensionalization........................ 9
1.1.8 Dimensionless Parameters .................... 11
1.2 Literature Review .......................... 12
1.2.1 Vapor Bubble Growth or Collapse in a Two-Phase Flow 12
1.2.2 Numerical Methods for Moving Interface Problems .. 22
1.3 Research Objectives ......................... 28
1.3.1 Research Objectives . . . . ... .. 28
1.3.2 Scales and Dimensionless Parameters . . ... ..30

2 NUMERICAL METHOD . . . . . ... .. 32

2.1 Fixed-Grid, Sharp Interface Method . . . ... .. 32
2.2 Interface Representation Algorithms . . . ... .. 33
2.2.1 C2 Cubic B-spline Algorithm . . . .... .. 39
2.2.2 Fairing Algorithm . . . . ... .. 44
2.2.3 Validation of Geometric Algorithms . . ... ..48
2.3 Cartesian Grid Method for Sharp Interfaces . . ... ..51
2.3.1 Fractional Step Method for Unsteady Transient Solution. 59
2.3.2 Discretization on Regular and Cut Cells . ... ..63
2.3.3 Imposition of Boundary Conditions on the Interface 74
2.3.4 Inversion of Discrete Operators . . . ... .. 76
2.4 Moving Interface Algorithms . . . . ... .. 80
2.4.1 Change of Phase . . . . . ... .. 81
2.4.2 Interfacial Conditions . . . . ... .. 82
2.4.3 Iterative Process for Determining the Interface Shape .. 84









2.4.4 Determining Interface Shape with Phase Change . 86
2.4.5 Discontinuity of Material Properties . . ... ..87
2.5 Discretization on Axisymmetric Geometry . . .... ..88

3 VALIDATION OF CARTESIAN GRID METHOD .. ......... 91

3.1 Validation of Global Accuracy of Cut Cell Approach ...... .. .91
3.2 Fidelity of The Numerical Methodology . . .... .. 93
3.3 Summ ary . . . . . . . ... .. 98

4 BUOYANCY-DRIVEN BUBBLE MOTION . . ... ..100

4.1 Summary of Previous Study . . . . ... ..101
4.2 Grid Resolution Study and Convergence Criteria . ... ..102
4.3 Results For Buoyancy Driven Bubble Motion . . ... ..103
4.4 Sum m ary . . . . . . . ... .. 115

5 BUBBLE RISE AND GROWTH . . . . ... ..116

5.1 Stationary Bubble Growth . . . . ... ..116
5.1.1 Case with Thermal Properties of Water . .... ..117
5.1.2 Case with Thermal Properties of Ethanol . ... ..119
5.2 Translating Bubble Growth . . . . ... ..119

6 BUBBLE RISE AND COLLAPSE . . . . ... ..132

7 SUMMARY AND FUTURE WORK . . . . ... ..138

7.1 Summary of Present Work . . . . ... ..138
7.2 Contributions of Present Work . . . . ... ..139
7.3 Future W ork . . . . . . ... .. 140
7.3.1 Applications . . . . . ... .. 140
7.3.2 Numerical Enhancements . . . ... .. 141

APPENDIX

A HEAT TRANSFER-CONTROLLED LIMITS OF SPHERICAL BUB-
BLE GROW TH ................................... 143

A.1 Conduction Dominated Growth ..................... 143
A.2 Convection Dominated Growth ..................... 144

REFERENCES ................................... 145

BIOGRAPHICAL SKETCH ............................ 151















NOMENCLATURE


CPL



CD

Fr

Fo

g

h

Ja

k,

k,,

L

n

Nu

Pd

Pi

Pv

Pe


heat capacity of the liquid phase (J/kg K)

heat capacity of the vapor phase (J/kg K)

drag coefficients

Froude number (= u2/gL)

Fourier number (= ait/(2Ro)2)

normal gravity

heat transfer coefficient (= q"/AT)

Jakob number (= pjcpAT/pvA)

thermal conductivity of the liquid phase (W/m K)

thermal conductivity of the vapor phase (W/m K)

length scale (m)

normal direction

Nusselt number (= hL/k)

dynamic pressure (kg/m s2)

pressure in the liquid phase (kg/m s2)

pressure in the vapor phase (kg/m s2)

Peclet number (= piCpjUrL/k1)








Pr Prandtl number (= ii/alt)

q" heat flux (W/m2)

r radial coordinate

Ro initial bubble radius (m)

Re Reynolds number (= piurL/pM)

t time (s)

T temperature (K)

Tint temperature at the bubble interface (K)

T, temperature in the liquid phase (K)

T" temperature in the vapor phase (K)

Ur velocity component in the r direction, or reference velocity scale (m/s)

uz velocity component in the z direction (m/s)

U terminal velocity for the rising bubble (m/s)

(U,)int interface velocity in the normal direction (m/s)

We Weber number (= piu2L/a)

z axisymmetric coordinate


Greek Symbols

a, thermal diffusivity of liquid phase (m2/s)

S under-relaxation factor for interface updating procedure

7 dimensionless radius for collapsing bubble (= R(t)/Ro)

K, curvature of the interfacial curve








AT degree of superheat or subcooling (K)

A latent heat of phase change (J/kg)

Pt dynamic viscosity of the liquid phase (kg/m s)

AV dynamic viscosity of the vapor phase (kg/m s)

p density of the liquid phase (kg/m3)

Pv density of the vapor phase (kg/m3)

a surface tension coefficient (N/m)

T stress tensor

TH dimensionless time (= 4Ja2 ait/7rRo2)

vi kinematic viscosity of the liquid phase (m2/s)

vV kinematic viscosity of the vapor phase (rnm2/s)

0 generic variable for illustration of finite volume integration














Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

DIRECT NUMERICAL SIMULATION OF A TRANSLATING VAPOR BUBBLE
WITH PHASE CHANGE

By

Tao Ye

May 2001

Chairman: Dr. Jacob N. Chung
Major Department: Mechanical Engineering

A fixed-grid direct numerical simulation method has been developed and verified

for single bubble dynamics, heat transfer and phase change phenomena. The original

contribution of this research is the first introduction of a sharp interface between the

vapor and liquid phases, which allows the simulation of density ratios up to the order

of 1000. The interface location is obtained as part of the solution which accounts for

all coupled dynamic and thermal effects.

The current method is believed to be the only method available at present which

resolves the curvature calculation issue, treats the interface as a sharp discontinuity,

honors the mass conservation principle, is capable of handling the large property

ratios and phase change.

For the isothermal cases, the predictions of drag coefficients and the shape defor-

mation of bubbles for Reynolds numbers ranging from 2 to 100 and Weber numbers

from 2 to 20 were shown independently to agree with previous published results.









It was also demonstrated that the current method is capable of predicting accu-

rate results for a wide range of Re, We, Ja, Pe numbers and property jumps at the

interface.

For a stationary and growing bubble, the predicted growth rate was found to

agree with the theoretical limit of oc t1/2. However, for a rising and growing bub-

ble, the predicted growth rate was approaching the theoretical limit of c t2/3 for a

spherical bubble but did not reach it exactly owing to bubble deformation.

The empirical correlation proposed for bubble collapse is assessed by the simu-

lations in the present work. The present numerical model predicts that the dynamics

of a collapsing bubble is qualitatively, but not quantitatively similar to the empirical

correlation proposed by other investigators.














CHAPTER 1
INTRODUCTION

1.1 Problem Statement
1.1.1 Research Motivation

Liquid-vapor phase change is a phenomenon that is frequently encountered in

daily life and that plays an important role in the power, chemical, petroleum and

electronics industries among others. For example, the power generation industry

takes advantage of the high heat transfer rates associated with the liquid-vapor phase

change in the heterogeneous nucleate boiling process to efficiently extract energy

stored in fossil, nuclear or solar fuels. The central mechanism of heat transfer during

nucleate boiling is the so-called ebullition cycle: a complete process of liquid heating,

nucleation, bubble growth and departure (see Carey, 1992). Hence understanding

and prediction of the microscopic vapor bubble behavior have direct importance for

better applications of nucleate boiling heat transfer.

Despite the wide applications of liquid-vapor phase change, a fundamental un-

derstanding of such phase change process is far from being complete. Analytical

analysis of the liquid-vapor phase change has faced mathematical difficulties because

of the nature of nonlinear, partial differential equations which govern the process.

In scientific and engineering research, experimentation and numerical simulation

have been mainly employed as basic methods of approach. Both means have respec-

tive merits and they need to complement each other. Meaningful experiments are of

singular significance for phenomena that are beyond the capability of current compu-

tational technology, and they are worthwhile in providing verification for numerical









computation and guidance for mathematical modeling. On the other hand, numeri-

cal computations can be used to predict microscopic transport mechanisms that are

difficult to measure and visualize in an experiment. Because of limited resources,

most of the experimental projects may be restricted to only one working fluid, simple

geometry, and limited thermodynamic conditions. Numerical computations might

therefore be useful for expanding the data base to a wider spectrum and providing

the information that is otherwise hard to obtain.

Although direct simulations of industrial scale phase change processes are cur-

rently out of reach, a fundamental understanding of the complex physics at small

scale can provide much insight to quantitative predictions of large scale aspects of

processes.

Despite the fact that computational fluid dynamics has advanced to a point

where solutions of relatively complex flows without phase change can be routinely

obtained, numerical techniques for the combination of fluid flow, moving interface

and phase change remain as one of the unresolved challenges since it is necessary to

model the interaction of fluid flow with heat and mass transfer across moving phase

interfaces.

In the present research, our major goal is thus to develop a numerical capability

for simulating liquid-vapor phase change at the small scale, namely, the vapor bubble

behavior in a pool of liquid.

1.1.2 Physics of a Boiling Vapor Bubble in Two-Phase Flows

Nucleate boiling is one of the most efficient heat transfer modes. This mode of

thermal energy transport is accomplished by two integral stages, which transport the

thermal energy from a heater surface to the bulk of a working fluid. The first stage is

the nucleation and growth of bubble embryos on the heater surface. The second stage

is the migration of the bubbles into the bulk of the fluid. The present research is

devoted to the fundamental understanding of the second stage.









In the second stage, if the bulk of the liquid is superheated, the bubble will

continue to grow while it is translating. The growth rate of the bubble in this stage

is primarily limited by the rate of arrival of the heat necessary for vaporization at

the bubble surface. This stage of growth is called heat transfer-controlled growth.

The pressure at the bubble surface inside the bubble approaches the liquid pressure

adjacent to the bubble surface; hence the temperature at the bubble surface is ap-

proximately equal to the equilibrium saturation temperature of the system pressure.

Translating vapor bubbles in a superheated or subcooled liquid environment can

be found in both pool boiling and convective boiling. In pool boiling, bubbles rise

into the bulk of the pool due to buoyancy. In convective boiling, bubbles are carried

into the bulk by flow currents. The interaction between the vapor bubble and its

liquid surroundings takes place after the departure of this bubble from the heater

surface. The interaction between the dispersed phase (bubbles) and the continuous

phase (liquid) involves exchange of momentum, thermal energy and mass. The in-

teractions are two-way coupled which means that the bubble behavior is affected by

the liquid flow and the liquid flow is also modulated by the presence of the bubbles.

These two-way coupled interactions include momentum exchange through the drag

and vaporization/condensation, energy exchange due to convection and latent heat

transport of phase change, and mass exchange during phase change.

The present research will focus on numerically simulating the single vapor bubble

growing or collapsing due to phase change while translating in a liquid.

1.1.3 Schematic of Computational Model

The schematic of computational model for a growing or collapsing vapor bubble

rising in a quiescent liquid is illustrated in Figure 1.1. The motion of the vapor

bubble is simulated in a cylindrical domain filled with liquid of the same substance.

The liquid is superheated in one case or subcooled in the other case. The major
simplification is that the system is assumed axisymmetric.






4




4z










Liquid





u
Bubble








r



Figure 1.1: Schematic of computational model.

1.1.4 Assumptions

The two-fluids model is adopted as the mathematical basis which governs the
liquid-vapor two-phase flow of a single translating vapor bubble; i.e., a separate set
of mass, momentum and energy conservation equations are solved for each phase








(dispersed vapor phase and continuous liquid phase surrounding the bubble) while
sharp discontinuities of material properties are maintained across the interface of zero
thickness.
The assumptions introduced throughout this study are summarized as follows:

Axisymmetric, cylindrical coordinate system

Single component system

Newtonian, constant property fluids in each phase

Thermal equilibrium at the phase interface which means temperatures at the
interface in the liquid and vapor phases are equal

1.1.5 Governing Equations

With the above assumptions, the mass, momentum and energy conservation
equations for both phases can be written as

V- u= 0 (1.1)


P + V (uu) = -Vp+ V (pVu) (1.2)

5 9TI

pcp [ + V (uT) = V.(kVT) (1.3)

where for numerical efficiency, pg = -A(pgh) has been combined with p in place
of the pressure. If the actual pressure is needed, hydrostatic pressure pgh has to be
subtracted from the pressure computed in the equations. Thereby the hydrostatic
pressure (body force term) contribution will only appear in the normal stress balance
at the bubble interface (see Eq. (1.12)).
For this research, the axisymmetric cylindrical coordinates are employed. Equa-
tions (1.1)-(1.3) thus become the following with 9/0aO = 0 and the swirl velocity










10 0
(ru'r) + (u) = 0
r Or Oz
[_ r 10 0 1 Opl10 0
at'[O + 1or(r"r.) + O! (UzU,) = -- r + (rr-) + (rz)
-ru + -(9 r5rT z


[Ou_ I!k 0 1
9 + r0 (ru,,) + (UU9)J



p[E + (ruT) + (u)]


Op+ 1 09 0
= -- (rT"z)+~()
8Z r Or Z


a (rkr) 0 T(kT)
= TOTr- 9 z 8z


(1.6)


(1.7)


where the stress tensor components are


9Ur
rr, = 2p-r

Trz = Tzr = I

r,2 = 2M--z
9z
ree = 2pU-
r


OUzr + Ouz r
(r '9~Z)


(1.8)


1.1.6 Boundary and Interfacial Conditions

Apart from the governing equations, boundary conditions are needed to make
the problem well posed. The interfacial conditions, which are simply the statement
of those conservation principles at the interface, must also be satisfied (see, e.g.,
Delhaye, 1974; Ishii, 1975; Kataoka, 1986; Sadhal, Ayyaswamy & Chung, 1996).

Boundary conditions:

The outflow (zero gradient of the velocity) conditions are specified on the top,
left and right sides of the domain. At the bottom boundary, the symmetric condition
for all variables is used.


Uo = 0:


(1.4)


(1.5)








Interfacial mass conservation conditions:

Due to phase change, for an interface moving at a velocity ujnt, the mass con-
servation conditions on the interface stipulates that

rh = p (ui,,t uj) -n = p, (uint u,) n

:=:: T = pt [(un)int (Un)j] = P [(un)int (Un),] (1.9)

where rh is the interracial mass flux (positive for evaporation), P1 and p, are the liquid
and vapor densities respectively, while the tangential velocity is always continuous

urt=u u-t = (ut)1 = (Ut), (1.10)

Interfacial momentum balance conditions:

The forces balance across the interface is expressed as

p. [(Un), (Un)int.] u + pT -rn = P, [(u)l (Un)int] U1

+pin Tr.n + own

(1.11)

The normal component of the equation (1.11) is

P1 p, + U-K = (-r.n)n (-r,-n)n {pI [(Run) (Un),nt] (Un)i -

Pv [(Un)v (un)int] (Un),} (1.12)

where r represents the viscous stress tensor, a is the surface tension coefficient which
is assumed to be constant, n = 1 + with R 1 and R2 the principal radii of
is aRb oR2) Rb and
interface curvature. The pressure in Eq. (1.12) is the total pressure including the
hydrostatic contribution; i.e., the actual pressure is obtained by subtracting pgh (or
1 h in dimensionless form) from p solved in Eq. (1.1), (1.2).

The normal component of stress balance is reduced to Young-Laplace equation
when there is no motion in both phases.








Interfacial energy conservation conditions:

The energy conservation relation across the interface involves both the sensible
and latent heat portions

(k,,VT,, k1VT1) -n = -Ap, (u,, Uint) n

= k, (') k,, (') = Ap, [(un)int (Un)v1 = A rh (1.13)

where A is the latent heat of evaporation.
Equation (1.13) is used to estimate the normal velocity component of the inter-
face driven by the phase change.

Interface temperature conditions:

By thermal equilibrium, the temperature continuity condition reads as

T1 = Tv = Ti, (1.14)

where Tint is the temperature at the phase change interface. A condition on the
interface temperature must be specified to complete the formulation.
In the case where liquid boiling occurs on a horizontal flat surface, i.e., curvature
n = 0, the interface temperature Tint is equal to saturation temperature Tsat of
the two phase mixture at the corresponding ambient pressure po. For a curved
interface between the liquid and vapor bubbles, when surface tension effect is taken
into account, the thermal and mechanical equilibrium at the interface reveals that the
local saturation temperature Tit along the interface will deviate slightly from that
constant saturation temperature Ts,,t for the fiat interface (see Alexiades & Solomon,
1993), which is the so-called Gibbs-Thomson effect. For the heat transfer controlled
bubble growth, vapor pressure pv at the bubble interface is equal to the system
pressure p, which is the equilibrium pressure of Tsat, i.e.,


Pv =Poo


(1.15)








In this case, a sufficiently accurate formulation for the interface phase change tem-
perature reads

Tint = Tsat 1 + (1.16)

1.1.7 Nondimensionalization

To nondimensionalize the governing equations in the liquid and vapor phases
as well as boundary and interfacial conditions, the reference scales are length L,
velocity Ur, time tr = L/ur. The characteristic temperature scale is AT = Th Tc =
{o T8~t growth
T, T, collapse Too and Tsat are, respectively, ambient temperature and
Tsat- Tlo collapse
saturation temperature.
With the basic scales chosen as above, the dimensionless variables are defined
as x* = x/L, u* = u/u, t* = t/tr, p* = p/(piu72), T* = (T- T,) /AT,
9* = g/g, P* = PIPref, = -l/Iref, k* = k/kref, Cp* = Cp/Cpref, where the

liquid properties are taken to be the reference properties, i.e., pref = p1, Pref = Al,
kref = ki, CPref = Cp.
Finally, the dimensionless equations in the liquid and vapor phases, omitting the
superscript *, can be written as

Liquid Phase:


V.u = 0 (1.17)

r\1 2
+V- (uu) = -Vp + V2u (1.18)
Ot Re
OT1
S+ V. (uT) = V2T (1.19)
Ot Pe

Vapor Phase:


V.u = 0 (1.20)









(P) [ VU.(UU)] =(,V)


(Pu) +V- (UT)]




Reynolds number

Peclet number


= I^Jv2T
(I) (~i ) 1Fe


pte= urL
Re pluL

Pe pcurL
k,


Interfacial Condition:

Mass continuity condition

(U(), ( ).
(U"n~int A --- 7 ~ --
1- P

Momentum balance condition


1
P1 Pv + We',-
IWe


Re [-' On L \On /v

[(u), (un),.J] (Un)j + () [(un), (Un),t] (u),
(1.24)


where p, = Pd + 1/Fr g(zut z), Fr = u'/gL is the total pressure including
hydrostatic and dynamics pressure. Zout is the level of the liquid pool.
Energy conservation condition


(uint (U)v = pJa


,n (k )On


(1.25)


Interface temperature condition


Tint (= F


1 V2U
-RV
Re


where


(1.21)


(1.22)


(1.23)


(1.26)






11


where

pu2 L
Weber number We =- r
or
2
Froude number Fr = -
gL

Jakob number Ja = C ,,AT
Pv A
Tsat '
Parameter F =
AT pIAL

1.1.8 Dimensionless Parameters

After examining the dimensionless governing equations and the associated bound-

ary and interfacial conditions, the major dimensionless parameters, namely, Re, Fr,

Pe, Pr, We and Ja, and their respective physical meanings are listed as follows:

Re = puL Reynolds number
itt'
flow inertia force
viscous force
2
Fr = Froude number
gL'
flow inertia force
gravity force

Pe = PtcurL Peclet number
kt
convection
conduction
piUrL /IItcp,
MIP k,
= Re- Pr

Pr = -, v=-, a Prandtl number
al PA PICP
momentum diffusion capacity
thermal diffusion capacity

We = puL Weber number
a
flow hydrodynamic pressure force
surface tension force









Ja = cpAT Jakob number
Pv A
sensible heat transfer
latent heat transfer

1.2 Literature Review
1.2.1 Vapor Bubble Growth or Collapse in a Two-Phase Flow

Since we are concerned with the heat transfer-controlled growth or collapse of

a single vapor bubble translating in a liquid pool, the review of previous works will

be focused only on those directly relevant to thermally controlled bubble growth or

collapse after the bubble has departed from the heater surface. A number of excellent

publications in the literature are dedicated entirely to the fundamental understand-

ing of the bubble growth dynamics on the heater surface (first stage of boiling), for

example, Klausner, Mei & Bernhard (1993); Zeng, Klausner & Mei (1993b); Zeng,

Klausner, Bernhard & Mei (1993a); Chen, Klausner & Mei (1995) and Mei, Chen &

Klausner (1995). These works are not directly relevant to the scope of the proposed

research, and will not be discussed in detail further here.

Vapor bubble growth or collapse has drawn extensive study over many years

because of its theoretical implication to the nucleate boiling and other two phase

physical phenomena. Due to its intrinsic difficulties, mainly the nonlinearity of con-

vection and the coupling of the bubble evolution with the surrounding liquid motion

as well as the large number of variables involved, it would seem that this problem

has been resolved to a tenable degree of understanding.

The heat transfer controlled bubble growth or collapse can be generally divided
into four basic situations: stationary bubble growth, translating bubble growth, sta-

tionary bubble collapse and translating bubble collapse. The researchers in this field

seem to follow this line to report results for one or more cases in any particular paper.

In our discussion, we will also follow this convention to look at the prior contributions

to all these cases.








Stationary bubble growth

The majority of the research on the heat transfer controlled growth of station-
ary vapor bubbles is based on the simplified model that a spherical bubble in an
unbounded fluid undergoes spherical growth. The simplified model does not solve for
the hydrodynamics of the fluid flow, nor does it locate the moving phase interface.
With the assumption of spherical symmetry, the continuity and momentum equa-
tions for an incompressible fluid lead to the famous Rayleigh equation (see Plesset &
Prosperetti, 1977)
d2R 3 (dR\2 1 f p(T) 1dR\ .2
R -jt _2 + 2 -t) _PV) ~o -411 ~ ~ (1.27)
2 A Rk R -d
in which T, (t) is the liquid temperature at the bubble surface and it is assumed that
the vapor pressure is the thermodynamic equilibrium pressure at the temperature T,.
To analyze the heat transfer controlled growth, the liquid temperature T, at the
bubble surface is obtained by solving the energy equation
RT RdR OT a, 0 (1.28)
-t + r2 dt r r2 Or \ Or ]
r 2dU ar~
in which a, = kl/(plcp,) is the liquid thermal diffusivity.
Plesset & Zwick (1954), Forster & Zuber (1954), Birkhoff, Margulies & Horning
(1958) and Scriven (1959) considered the heat transferred controlled spherical vapor
bubble growth without translation by solving the energy equation (1.28) subject
to certain approximations in conjunction with the dynamical equation (1.27) and
thermodynamic equilibrium relation p, = p,(T,). Their results of bubble radius
varying with time are in good agreement with the experimental data of Dergarabedian
(1953) and, later, Florschuetz, Henry & Khan (1969); basically, what they concluded
is that the bubble radius of the heat transfer controlled growth is proportional to
(time)1/2.
Prosperetti & Plesset (1978) re-examined two assumptions introduced in those
earlier analyses of bubble growth: the validity of a thin thermal boundary layer near















E
o o, .umc-u-U I"
4.0E-03 -

S3.0E-03 -

2.0E-03 -

1.OE-03

O.OE+O00 0.1 0.2 0.3 0.4

Time (sec)

Figure 1.2: Bubble radius vs. time for a vapor bubble growing in superheated water
at atmospheric pressure with 5C superheat [Prosperetti & Plesset (1978)].


the bubble wall and the empirical linear relation for the vapor pressure with vapor
temperature. It was showed that, upon comparison with numerical results of Donne
& Ferranti (1975), both assumptions do not lead to serious errors in the radius-time
behavior for superheat ranges of practical interest. The relation of the bubble radius
with time for heat transfer controlled growth is given as

R = 2k, --) (27)t (1.29)

in which k, is the liquid thermal conductivity, a, = ki/(pip,) is the liquid thermal
diffusivity, (Too T8,t) is the liquid superheat, A is the latent heat, Pv is the vapor
density. The radius increasing with time based on this formulation is plotted in
Figure 1.2.









Translating bubble growth

The stationary vapor bubble is an ideal exception. In practice, the vapor bubble

usually exhibits more or less translator motion. Although the translator motion is

not expected to play a significant role if liquid inertia is the dominant mechanism,

as stated by several investigators, it has a large influence on the growth or collapse

if heat transfer is controlling the phase change. Under such circumstances, the heat

transfer would be enhanced by translation, thus resulting in higher growth or collapse

rates.

The experimental investigation of Darby (1964) studied simultaneous growth and

translation characteristics of vapor bubbles issuing continuously from a nucleation

site. The bubble radius was found to vary as (time)2/3 during the entire recorded

history. It indicates that, as expected, the translational motion increases the growth

rate over that for stationary bubbles.

Ruckenstein & Davis (1971) assumed a potential flow in the region surrounding

the bubble to analyze the effect of bubble translation on vapor bubble growth in a

superheated liquid. He found that for a moderate Jakob number, the translation can

substantially increase bubble growth rate over the rates predicted for a stationary

growing bubble. For sufficiently large Jakob numbers (Ja > 50), bubble growth is not

greatly affected by the relative motion between the bubble center and the surrounding

liquid. The results predicted by Ruckenstein & Davis (1971) are in good agreement

with the experimental data of Florschuetz, Henry & Khan (1969).

In Table 1.1 we summarize the above mentioned important research publications

for heat transfer controlled growth of a vapor bubble.

Stationary bubble collapse

In bubble collapse, the radius variation produced by a given interfacial heat flux

increases in time in contrast to the decrease for the bubble growth case. This interface










Table 1.1: A list of representative investigations on stationary and
bubble growth


translating vapor


Authors Approach Assumptions & Conditions Major Results
Plesset & Zwick Analytical Spherical symmetry The bubble
(1954) Stationary growing vapor bubble radius is pro-
Forster & Zuber Heat transfer controlled vapor bub- portional to
(1954) ble growth the square root
Birkhoff, Mar- The inertia and viscosity of the liq- of time, i.e.
gulies & Horning uid are neglected (time)1/2
(1958) Surface tension is neglected
Scriven (1959) Heat transfer by conduction only
Prosperetti & Constant thermal and physical
Plesset (1978) _______properties_______
Darby (1964) Experiment l Translating vapor bubble growth Bubble ra-
dius varies as
_________________(time)2/3
Ruckenstein & Analytical Spherical symmetry The growth
Davis (1971) Translating vapor bubble growth rate is substan-
Potential flow in the liquid tially increased
Heat transfer controlled growth at moderate
Jakob number
____________________________________(Ja < 50)
Legendre, Boree Numerical Spherical symmetry Confirmed all
& Magnaudet Full Navier-Stokes equations are above prior an-
(1998) solved alytical results
Heat transfer controlled growth and studied the
Stationary and translating bubble hydrodynamic
Decoupled bubble evolution and forces
liquid motion


behavior makes condensation basically more difficult to describe than the evapora-

tion. For the stationary bubble collapse, by virtue of spherical symmetry assumption,

Florschuetz & Chao (1965) examined the relative importance of the liquid inertia and

heat transfer effects on the vapor bubble collapse. A dimensionless parameter is iden-

tified to characterize the collapse mode: either liquid inertia-controlled collapse, heat

transfer dominated collapse or an intermediate case where both liquid inertia and

heat transfer are of comparable importance. For heat transfer-controlled collapse,









they found that the collapse rate is relatively slow and decreases as the collapse pro-

ceeds, the radius-time relation is bounded by two limiting heat transfer-controlled

collapse curves that resulted from two approximations for temperature solutions:

-Y =1--rH (1.30)


and

=H + ( 2-3) (1.31)

in which dimensionless bubble radius 7y = R/Ro, dimensionless time

TH = 4 Ja2 lt (1.32)

with Ja being the Jakob number, a, the liquid thermal diffusivity, R0 the initial
radius of the bubble. The results are visualized in Figure 1.3 which is from the paper

of Florschuetz & Chao (1965).
The analytical results Eq. (1.31) of Florschuetz & Chao (1965) are valid at high

Jakob numbers. In order to obtain results valid at finite Jakob numbers, Okhotsimskii
(1988) solved numerically the heat transfer controlled collapse based on the same
theoretical model of Plesset & Zwick (1954) and tabulated the time evolution of the

bubble radius over a wide range of Jakob numbers ( 10-2 < Ja < 103 ). In his

analysis, he isolated the Jakob number to be the only dominant factor of the collapse

process.

Translating bubble collapse

For more practical cases of translating vapor bubble collapse, Wittke & Chao
(1967) obtained numerical and experimental information of the influence of transla-
tory velocity on the collapse of a vapor bubble. Noting that the translator motion is

not expected to be a significant factor if liquid inertia is the dominant mechanism for
collapse, they examined the bubble traversing at a constant velocity collapsing under




























Figure 1.3: Bubble radius vs. time for heat transfer controlled vapor bubble collapse
[Florschuetz & Chao (1965)].


heat transfer controlled conditions. They found that, again as expected, because the
relative contribution to heat transfer at the bubble surface due to translator mo-

tion is greater when the Jakob number is smaller, the translator motion produces a
greater effect for the smaller Jakob number. Overall, the faster the bubble translates,

the more rapidly it collapses.
Tokada, Yang & Clark (1968) and Dimic (1977) also studied theoretically the ef-
fect of translational velocity of a bubble on the collapse rate. They basically obtained

similar results on the collapse rates.
In these studies, the bubble is relatively large (2mm < Ro < 4mm) and because
of reasons which are not completely clear, it was observed in some experiments on
bubble collapse that bubbles exhibit almost a constant rising velocity (see Wittke
& Chao, 1967; Chen & Mayinger, 1992); hence a constant translation velocity is
assumed in all these studies.








Moalem & Sideman (1973) analyzed the effect of variable rising velocities as-

sociated with relatively small bubbles (R0 < 1 mm) as well as the effect of a cross
flow; their results were in good agreement with the experimental data of Abdelmes-

sih, Hooper & Nanngia (1972). What they concluded is that the effect of variable

velocity as compared to a constant velocity motion is relatively small.

All of the above mentioned results assumed the potential flow for outer liquid

surrounding the bubble.

Chen & Mayinger (1992) performed very sophisticated measurements to deter-
mine the heat transfer at the phase interface of a collapsing vapor bubble rising in

a subcooled liquid. The temporal decrease in the bubble diameter is given in the

correlation using experimental data as

7 = (1 O.56ReO7PrO5JaFo)0'9 (1.33)

in which 7-y = R/Ro, Reynolds number is based on initial bubble diameter. The

Fourier number is
ailt
Fo = al (1.34)
(2Ro)2
Recently, Gumerov (1996) numerically obtained the solutions of the same theoretical
model as that of Wittke & Chao (1967) for the collapse of a vapor bubble with

constant translator velocity. His results are shown to agree substantially better

with the experiments in a number of cases.
In Table 1.2 some of these studies for stationary and translating vapor bubble
collapse are summarized.
Most recently, Legendre, Boree & Magnaudet (1998) made a comprehensive
study of heat transfer-controlled spherical bubble growth and collapse without and
with constant translation velocity. Basically, they reexamined all four cases, and by
comparing with those representative prior results, confirmed previous analytical anal-
ysis. In their study, a numerical approach is employed to solve the full Navier-Stokes










Table 1.2: A list of representative investigations on stationary and
bubble collapse


translating vapor


Authors Approach Assumptions & Conditions Major Results
Florschuetz & Analytical Spherical symmetry Three collapse
Chao (1965) and Stationary collapsing bubble modes are
Experimental Surface tension and viscosity are ig- classified and a
nored dimensionless
No relative velocity between liquid parameter is
and bubble wall identified to
Heat conduction in vapor is ignored characterize
Pressure in the vapor phase is uni- the collapse
form mode.
Perfect gas for vapor
Okhotsimskii Numerical Spherical symmetry Time evolu-
(1988) Stationary collapsing vapor bubble tion of bubble
Heat transfer controlled vapor bub- radius for a
ble collapse wide range
The inertia and viscosity of the liq- of Jakob
uid are neglected numbers is
Surface tension is neglected tabulated.
Constant thermal and physical
properties
Wittke & Numerical Spherical symmetry The transla-
Chao (1967) and Translating vapor bubble tion motion
Experimental Heat transfer controlled collapse generally in-
Constant translational velocity creases the
Potential outer flow collapse rates.
Incompressible for Liquid
Chen & Experimental Translating vapor bubble Correlations of
Mayinger Heat transfer controlled collapse Nusselt num-
(1992) ber and time
evolution of
bubble radius
_________ ~_________are given.
Legendre, Numerical Spherical symmetry Confirmed all
Boree & Full Navier-Stokes equations are above prior
Magnaudet solved analytical
(1998) Heat transfer controlled collapse results and
Stationary and translating bubble studied the
Decoupled bubble evolution and hydrodynamic
liquid motion forces

























10R (b)

I F
I80R



nsuymadMo mI.dy corAdon


Figure 1.4: Computational configuration of investigation of Legendre, Boree & Mag-
naudet (1998).


equations and temperature equation in the liquid to determine the temperature field
around the bubble surface so as to obtain the growth or collapse rates. Despite the
spherical symmetry assumption and disregarding of the coupling between the bubble
evolution and the liquid motion, their solutions, to our knowledge, are the most ad-
vanced ones for the spherical bubble growth and collapse. The second contribution
of their research is regarding the hydrodynamic force experienced by the bubble with

a time-dependent radius.
The computational configuration of Legendre, Boree & Magnaudet (1998) is
shown in Figure 1.4.









1.2.2 Numerical Methods for Moving Interface Problems

The liquid-vapor two-phase problem we outlined in chapter 1 is a special case

of more general situations involving multi-fluids or multi-phases in a computational

domain as shown in Figure 1.5, where multiple immersed interfaces are embedded

in the flow field. The regions enclosed by the interfaces designate certain phases

which could be solid, liquid or gas phases and might be different from that of the

surrounding fluid. The interfaces can be stationary, or moving, deforming through

interaction with the flow around it. Such a scenario may be encountered in several

types of flow problems. For instance, the interface may represent a solid geometry

over which fluid flows. This solid boundary may be rigid as in pure external flow

problems or be flexible exhibiting motions under the influence of the external flow as

in many fluid-structure interaction problems. In other situations, the interface may

represent a phase boundary at which a phase transition occurs. One of the examples

is the liquid-vapor phase change problem involving vapor bubbles as in the present

research, the enclosed phase is the dispersed vapor phase, and the surrounding one

is the continuous liquid phase. In an inverse situation involving liquid drops, the

enclosed phase is the liquid phase, and continuous phase is the vapor phase. In

these cases, sharp discontinuities of the fluid properties such as density, viscosity and

thermal conductivity arise across the interface. Other examples of phase boundary

include solid-liquid phase of the solidification in material processing problems. The

interface may also simply represent a boundary separating two immiscible fluids as

in the case of bubbles or drops immersed in an ambient fluid.

Whether the immersed interface represents a phase interface or a boundary be-
tween two fluids with different transport properties, the numerical simulation of these

processes has to deal with moving interfaces. The numerical investigation of moving

interface problems has posed a difficult challenge for a long time, especially when a

moving phase change interface and extremely large ratios of fluid properties across the









12
D
I pPhase 0
Floo

1'10
all








Figure 1.5: Illustration of multiple immersed interfaces in a computational domain

interface exist simultaneously, as in the present study of the liquid-vapor phase change
problem. Despite its difficulties, in the past several decades, significant advances have
been made in numerical techniques for obtaining the solutions of immersed, moving
interface problems.
In this section, we will briefly examine the available numerical approaches for
immersed, moving interface problems with fluid flow.
Numerical techniques for immersed, moving interface problems could be classi-
fied into three broad categories (see Shyy, Udaykumar, Rao & Smith, 1996): body-
fitted grid approach, fixed grid approach and mixed Eulerian-Lagrangian approach.
The natural approach is to use body-fitted type mesh which could be either struc-
tured or unstructured (e.g., finite element mesh). Relatively complex two-dimensional
free boundary problems were simulated using this approach (see Ryskin & Leal, 1984;









Kang & Leal, 1987; Dandy & Leal, 1989). A body-conforming mesh makes the bound-

ary conditions specification on the interface relatively easy and accurate because the

interfacial conditions are applied at the exact location of the interface without any

smearing operation. The sharp discontinuity of material properties across the inter-

face is always maintained. However, in the body-fitted approach, the computational
expense is an issue because the grid has to be modified or regenerated continuously

to account for the movement of the interface. The body-fitted approaches encounter

major difficulty in dealing with the topological complexity resulted from arbitrary

interface movement with desired mesh quality. Besides, the governing equations are

not in a simple form for body-fitted mesh which causes extra difficulty for the over-

all solution procedure. Multiple connectivity makes designing an efficient algebraic

solver a difficult task.

Another completely different approach is to use the fixed Cartesian grid. In

this purely Eulerian approach, the exact location of the interface is not explicitly

known, in contrast, the interface details are determined from solving a field variable.

Within this approach, there are a few variants, for example, the Volume of Fluid

(VOF) method (see Hirt & Nichols, 1981; Brackbill, Kothe & Zemach, 1992; Kothe

& Mjolsness, 1992; Rider & Kothe, 1998; Scardovelli & Zaleski, 1999) the level-set

method (see Osher & Sethian, 1988; Sethian, 1990; Sussman, Smereka & Osher, 1994;

Sussman, Fatemi, Smereka & Osher, 1998) and phase-field method (see Caginalp,
1984; Kobayashi, 1993; Antanovskii, 1995; Anderson & McFadden, 1997; Jacqmin,

1999) amongst others.
When surface tension on the interface can not be ignored, the accurate estimation
of the interface curvature is crucial to successfully predicting the interface evolution.

Under these circumstances where explicit interface tracking is needed, the mixed
Eulerian-Lagrangian approach seems attractive because in this approach, apart from

the fixed structured grids (Eulerian part) for the flow field, interfaces are explicitly








described by curves or surfaces using computational geometry tools (Lagrangian part)

and the movement of those curves and surfaces are determined by the computational

process.

In the fixed Cartesian based approach, the grid generation is greatly simplified

and the governing equations retain the relatively simple form in Cartesian coordi-

nates. Because the grid is Cartesian, many efficient solvers are easily obtained to

reduce the computational cost.

In the mixed approach, a major issue is how to handle the interaction between

a moving interface and a flow field. The interface would cut through the underlying

Cartesian grid in an arbitrary manner which resulted in fragmented control volumes

(cut cells) around the interface. Fragmented, small cut cells would cause numerical

instability. To circumvent this difficulty, the so called immersed boundary technique

has been widely used to simulate multi-fluid flows (see Peskin, 1977; Unverdi &

Tryggvason, 1992; Goldstein, Handler & Sirovich, 1995; Juric & Tryggvason, 1996,

1998; Kan, Udaykumar, Shyy & Tran-Son-Tay, 1998). In this technique, the interface

is not kept sharp but is represented on the fixed underlying mesh by a thin layer of

finite thickness around the interface. A scalar indicator function is then constructed

from the known position of the interface to smooth out the jumps in the material

properties across this thin layer. The information transfer between the interface and

a flow field is done by introducing additional source terms in the governing equations.

This technique is quite robust, however its fuzzy interface treatment does not reflect

closely the physical reality of the zero thickness interface, and thus it is not suitable

for accurately resolving any boundary layer type of behavior in the thermal and flow

fields when convection is important.

A strictly conservative treatment of small cut cells is the so-called cut-cell tech-

nique. In this treatment, typically, small cut cells are absorbed into neighboring cells









to form mostly quadrilateral cells in the vicinity of the interface. Thereby, the inter-

face is treated as a sharp discontinuity and the conservation property of the solution

procedure is always retained. The interface effect is taken into account by the flow

solver through the reshaped control volumes near the interface.

The finite volume cut cell techniques have been used in dealing with the complex

solid boundaries in the flow field in conjunction with Cartesian mesh (see Quirk,

1994; Udaykumar & Shyy, 1995; Pember, Bell, Colella, Crutchfield & Welcome, 1995;

Udaykumar, Kan, Shyy & Tran-Son-Tay, 1997; Ye, Mittal, Udaykumar & Shyy,

1999). In those applications, the relative locations between solid boundaries and

grid cells are fixed throughout the solution process. There is no need to determine

the boundary shape in the solution process. There is no interaction between rigid

boundaries and flow field which makes it much easier to implement the cut cell

technique.

In the applications where free, moving interfaces with phase change, e.g., in-

terfaces separating bubbles or drops from the surrounding fluids, are involved. It

is much more challenging to implement the cut cell technique because the interface

shape is not fixed any more, rather depending on the dynamic conditions on the

interface. In other words, the instantaneous interface shape must be determined as

part of overall solution process based on the interaction of the interface with the flow

fields in both sides of the interface. Owing to this requirement, no reports has been

found in the literature on the application of finite volume cut cell techniques to free

interfaces like those of bubbles or drops. In the present research, a fixed grid, sharp
interface moving boundary method is developed and applied to simulate the moving

bubble interfaces.

To summarize, pure body-fitted grid approaches are seldom used nowadays, they

are more often combined with fixed grid methods to improve the efficiency. In such

configurations, most of the domain is covered by the fixed grid, body-fitted grid is









only used around the curved interface to capture the movement of the interface. A

communication mechanism has to be designed between global fixed grid and local

body-fitted grid.

Fixed grid methods are mainly aimed at industrial simulations. In these meth-

ods, normally a homogeneous system is applied to the entire domain. As a result, a

single set of governing equations are solved in a domain including multiphase with

different properties. The implementation of these methods is relatively easy. Despite

its simplicity, a sharp discontinuity for the transport properties at the phase front

should be preserved in conformity with the physical reality. In fixed grid methods,

the interface is not explicitly known and has a finite thickness, or sometimes, al-

though the interface itself can be implicitly captured by some field function to be

sharp, large property jumps across the two phases still have to be smoothed over a

few grid points around the interface to overcome the numerical stiffness caused by

the property jumps. In order to maintain a sharp front, the fuzzy interface should be

restricted to within a few grid cells: the less the grid cell, the better the resolution

of the interface.

Our research is intended to capture the local physics near the phase interface for

liquid-vapor phase change process. We developed a numerical method that is based

on fixed grid with explicitly tracked interfaces. At the same time, the explicitly

tracked interface is a sharp discontinuity of fluid properties. The current method is

a mixed Eulerian-Lagrangian approach. However, the major difference between this
method and other mixed Eulerian-Lagrangian methods is the interface treatment.

The sharp interface is maintained in the present method.

In the present method, first, an explicit interface representing and tracking pro-

cedure is developed to locate the interface exactly, and a underlying fixed grid is

employed to cover the entire domain for field solution. Secondly, the interface is kept

sharp in terms of property jumps. Thirdly, two set of governing equations are solved









in two phases, respectively in conjunction with corresponding transport properties.

As a result, the fixed grid over the entire domain is virtually divided into two sub-

domains using cut cell techniques. The coupling of the two phases is treated with

the implicit updating of the interface location under the constraints of interfacial

conditions.

1.3 Research Objectives
1.3.1 Research Objectives

The interests in the vapor bubble phenomena has stimulated continuous, inten-

sive efforts from numerous investigators on this subject for half a century, as described

in the previous sections. And liquid-vapor phase change (boiling or condensation)

process is one of the most difficult heat transfer modes to investigate. Because the

vapor bubble behavior is described by a set of coupled, nonlinear partial differential

equations and because there are not mathematical tools to solve those equations an-

alytically, the majority of the theoretical works on the vapor bubble dynamics in the

literature have to resort to more or less simplified theoretical models so as to make

the problem tractable. Although those simplified theoretical models sometimes do

yield quite good results in comparison with experimental data, it is certainly desired

and has been ceaseless attempted to have the theoretical model as close as possible

to that set of coupled, nonlinear partial differential equations. The present research

thus serves to bring that closeness to a level which has not been achieved before,

however, the gap between the theoretical model and the ultimate description is not

completed removed.

On the other hand, owing to experimental difficulties, only a few investigations
have focused on the micro-phenomena at the phase change interface. In order to

achieve a comprehensive understanding of the physics of boiling and condensation,

an effort is needed in the development of a direct numerical simulation technique that








can be used to study the microscopic transport mechanisms at the moving liquid-

vapor phase change interface.

The objective of the present research is thereby to develop a numerical method to

simulate the liquid-vapor phase change transport mechanisms of a translating vapor

bubble in a superheated or subcooled liquid. Based on the numerical simulation

results, insights are offered into the physical and numerical characteristics of the

bubble dynamics. An assessment of some conclusions or correlations proposed in the

literature regarding the bubble growth or collapse dynamics is also presented.

The development of the present fixed grid, sharp interface moving boundary

method consists of the implementations of the following major components:

Interface representation algorithms using B-spline curve fitting,

Fairing algorithm for continuous curvature calculations,

Intersection of the curved interface with grid lines,

Assemble and reshaping of the cut cells,

Discretization of governing equations on regular and cut cells,

Fractional step procedure for solving the coupled governing equations,

Imposition of the boundary conditions for velocity, pressure and temperature

at the interface,

Iterative algorithm for updating the interface location based on normal stress

balance conditions at the interface,

Algorithm for the interface evolution owing to phase change effect.

This challenging task is accomplished owing to the latest development of numeri-

cal techniques for solving time-dependent free boundary problems in fluid mechanics.








In recent years, it has become possible to solve coupled, nonlinear mass, momen-

tum and energy conservation equations (partial differential equations) numerically

because advanced numerical techniques and fast computer hardware were available.

For problems with even more physical complexity such as those involving multi-phase

fluid flow, heat and mass transfer and/or moving interfaces etc., the reported simu-

lation works on Navier-Stokes equations start to appear in the literature only in the

last ten years yet there are many numerically relevant issues which are not resolved.

In our case, the method for solving the free boundary problems is developed and

large scale computations are carried out to solve coupled Navier-Stokes and energy

equations in liquid and vapor phase respectively at the same time. On the phase

transition interface, the mass, momentum and energy conservation principles are

matched from the liquid and vapor sides of the interface.

The method developed in the current research is not restricted to only liquid-

vapor phase change problems, it is generally applicable to other phase change pro-

cesses such as solid-liquid phase change of melting or solidification because the

method is developed to handle property discontinuity and moving interfaces which

are common to the liquid-vapor and solid-liquid phase change.

1.3.2 Scales and Dimensionless Parameters

Several major non-dimensional parameters including Reynolds, Weber, Jakob

and Fourier numbers are varied to gain a better understanding of physical and nu-

merical aspects of the solution procedure. To determine the above dimensionless pa-

rameters, characteristic scales are chosen based on the problem features. The length
and velocity scales are determined as well as the maximum temperature difference,
AT, in the system. Naturally, we choose the initial bubble diameter, Do, as the

length scale L. The velocity scale can be defined based on the diffusion mechanism:

aj/L for stationary bubble growth cases, buoyancy effect: (gL)112 for rising bubble









Table 1.3:


Saturation properties of water at Psat = 101.3 kPa, Tsat = 373.15 K

Saturation Properties Liquid Vapor
p (kg/m3) 958.3 0.597
A (N s/rM2) 277.53 x 10-6 12.55 x 10-6
c (J/kgK) 4,220 2,030
k (W/m K) 0.679 0.025
Pr 1.72 1.02
A (J/kg) 2,256,700
a (N/m) 0.0589


Table 1.4: Saturation properties of Ethanol at Psat = 101.3 kPa, Tsat = 351.45 K

Saturation Properties Liquid Vapor
p (kg/m3) 757.0 1.435
i (N s/m2) 428.7 x 10-6 10.4 x 10-6
cp (J/kgK) 3000 1830
k (W/m K) 0.1536 0.0199
Pr 8.37 0.96
A (J/kg) 963,000
ao (N/m) 0.0177

with phase change cases, or the bubble terminal velocity U for bubble rising with no
phase change. The choice in each case will be specified individually.
The rest of the information required to specify the dimensionless parameters are
the thermo-physical properties of the liquid. An example of working fluids is water
whose saturation properties at the atmospheric pressure are listed in Table 1.3. Then
the typical magnitude of those dimensionless property ratios of the liquid and vapor
phases of water are Pi/P, = 1605.2, /t/ltv = 22.114, ki/k, = 27.16, Cpj/cp,, = 2.0788.
The range of Reynolds number we consider in this research is Re = 1 1000.
The bubble diameter is 0.5 mm 1.5 mm, and the maximum temperature difference
AT = 5C.
Besides water, we also consider another working fluid, namely, Ethanol (CH3CH2OH)
whose thermo-physical properties are listed in Table 1.4.














CHAPTER 2
NUMERICAL METHOD

In this chapter, the fixed-grid, sharp interface moving boundary method is de-

scribed. In this method, a underlying fixed Cartesian grid serves as the Eulerian

framework of the algorithm. Within this Eulerian framework, separate marker points

which are connected by piece-wise polynomials are adopted to represent the phase

interface. Those regions separated by the phase interface represent different phases.

The deformation and movement of the interface is realized through the translation of

these markers over the underlying stationary grid, which constitutes the Lagrangian

part of the algorithm. In each phase region, the fractional step method is employed

to solve the coupled governing equations of either liquid or vapor phase.

2.1 Fixed-Grid, Sharp Interface Method

The development of this numerical method consists of three major tasks:

Implementation of curve fitting and curvature calculation algorithms.

Development of Cartesian grid techniques for solving fluid flows over fixed,

immersed interfaces with complex geometries.

Development of moving interface algorithms, such as the interface updating

algorithm.

The content of this chapter is organized, as shown in Table 2.1, so as to describe
these three major components of the numerical method in a logic manner. In the

section 2.2 of this chapter, the implementation of interface representing, tracking and

curvature calculation algorithms are detailed. The C2 cubic B-spline curve fitting









Table 2.1: Organization of Chapter 2


Interface representing and tracking algorithms Section 2.2
C2 cubic B-spline algorithm Section 2.2.1
Fairing algorithm Section 2.2.2
Validation of geometric algorithms Section 2.2.3
Cartesian grid method for fixed, immersed interfaces Section 2.3
Fractional step procedure for solving coupled governing equations Section 2.3.1
Discretization of governing equations on regular and cut cells Section 2.3.2
Imposition of boundary conditions at the interface Section 2.3.3
Algebraic equation solver Section 2.3.4
Moving interface algorithms Section 2.4
Reinitialize the flow field due to phase change Section 2.4.1
Inclusion of interfacial conditions Section 2.4.2
Iterative Process for Determining the Interface Shape Section 2.4.3
Determining Interface Shape with Phase Change Section 2.4.4
Discontinuity of Material Properties Section 2.4.5


is employed in conjunction with fairing algorithm. In the section 2.3, numerical

procedures related to solving fluid flows over fixed, immersed interfaces on Cartesian

grid are described. In the section 2.4, the numerical techniques for moving interface

computations are presented.

The overall solution procedure for numerical simulation of moving interface prob-

lems is illustrated in the flow chart of the code in Figure 2.1 in which the major

function modules interface tracking procedure is detailed in Figure 2.2, fractional step

process is detailed in Figure 2.3, moving and updating interface is detailed in Figure

2.4. In the rest of the chapter, each section is dedicated to certain numerical issues

appeared in these flow charts as also shown in Table 2.1.

Currently, the source code of this solver has nearly 40,000 lines of Fortran90

statement implemented on Compaq Tru64 Unix work stations.

2.2 Interface Representation Algorithms

A representative schematic of the present fixed grid method is depicted in Figure
2.5. The grid is Cartesian and does not conform to the body surface, and the interface










Begin



Input computational parameters/
ut initial interface shape /



Interface tracking procedure (Section 2.2)



Time step n=O



Time step loop
Time step n+l



Moving and updating interface (Section 2.4)
Interface tracking procedure (Section 2.2)
Cartesian Grid solver (Section 2.3)



Until time step > specified number



(End)


Figure 2.1: The flow chart showing the overall solution procedure of moving interface
algorithm.































































bI

Figure 2.2: The flow chart showing the interface tracking procedure as detailed in
Section 2.2.














Initialize the flow field


Solve for intermediate flow field u* (Eqn. 2.12)
* Imposition of boundary conditions (Section 2.3.3)
* Calculate face fluxes using interpolation procedure (Section 2.3.2)



Obtain face center velocity U



Solve for pressure Poisson equation (Eqn. 2.15)
* Imposition of boundary conditions (Section 2.3.3)
* Calculate face fluxes using interpolation procedure (Section 2.3.2)


Correct cell center velocity field u*-* u"*'
Correct cell center velocity field U"-4 U*1
(Eqn. 2.16)


SSolve for energy equation


Figure 2.3: The flow chart showing the fractional step procedure as detailed in Section
2.3.














Iteration k=O



Interface update iteration
Iteration k=k+1


Move interface to new location (Section 2.4.3)



Modify the interface shape to satisfy mass
conservation in each phase (Section 2.4.3)



Interface representation procedure (Section 2.2)


Reinitialize the flow field due to phase change (Section 2.4.1)
Fractional step procedure (Section 2.3.1)


Terminated if all conditions are satisfied


Figure 2.4: The flow chart showing the interface moving and updating procedure as
detailed in Section 2.4.






























Figure 2.5: Schematic of Cartesian grid and interfacial marker points.


is defined explicitly by geometric curves in the computational domain. Basic elements

involved in defining the interface are interface marker points and interfacial curves

connecting the markers. The markers define the terminal points of the interfacial

curves. Given a set of markers, finding a curve which fits all markers is a geometric

interpolation process. Depending on geometric conditions imposed at the marker

points, there are various ways to define numerically the interface characteristics.

In the present work, the C2 cubic interpolatory B-spline curve is employed to

fit a set of marker points. The B-spline curve refers to a set of B6zier curves glued

together at the marker points to form a composite curve representing the interface.

Each piecewise Bezier curve is actually a polynomial curve expressed mathematically

for convenience in special form: Bernstein polynomials. There are efficient algorithms

to do curve manipulation and geometric calculations. Cr denotes the smoothness con-

ditions at the junction points of piecewise Bezier curves, requiring that a composite









curve is r times continuously differentiable at the junction points. C2 thereby means

the composite curve has continuous second derivatives at any marker point. In order

to achieve this level of smoothness, the individual Bezier curve must be at least degree

3 (cubic) polynomial to have continuous second derivatives, and second derivatives

computed from both sides of the junction point have to be equal. This mathematical

requirement of the two polynomials yields a smooth global curve to represent the

entire interface.

Once the B-spline fitting curve of the interface is constructed, the geometric

informations such as location and curvature of any point along the entire interface

can be easily obtained.

The translating, deforming, expanding or shrinking of the interfaces are realized

through the motion of each individual marker point which in turn is determined from

the flow quantities on underlying fixed grid using, e.g., normal stress balance condi-

tions. With the movement of the markers, the instantaneous B-spline representation

of the interface is reconstructed accordingly to keep track of the interface. For a

interface in motion, the curvature calculated based on the continuously constructed

B-spline may exhibit numerical oscillations (see Farmin, 1997). To extract the correct

curvature values along the interface, a so-called fairing algorithm is adopted. As will

be presented, the combination of the cubic B-spline and the fairing algorithm results

in a robust and accurate method for tracking highly distorted interfaces in terms of

location as well as curvature.

2.2.1 C2 Cubic B-spline Algorithm

The interface representation using B-spline is based on the marker locations com-

puted at every time instant. The interfacial marker points are indexed sequentially

and can represent any number of open or closed interfaces based on the assigned

connectivity between them. There is no restriction on the number of interfaces or

on the end conditions that can be imposed on these interface curves. As shown in






























Figure 2.6: Illustration of immersed interfaces, marker points and normal convention.


Figure 2.6, the interfaces can be open or closed, stationary or moving and can enclose

different phases or materials with different properties. The locations of the marker

points are defined by coordinates X(s) which is parameterized in terms of arc length

s and a distance ratio based on the grid spacing. In our investigation, the marker

spacing is initially assigned to be the same as the grid spacing, h. In the course of

computation, markers are redistributed after each time step according to the initial

criterion.

The convention adopted by the present algorithm in indexing the marker points

is that, by choosing the normal of the interface to point from phase 1 to phase 0, the

phase 1 always lies to the right as one traverses the interface along the sequence of

the marker points.

In order to obtain geometric information, we need to construct a smooth curve

fitting this set of marker points. In the computer-aided geometric design applications,

the B-spline offers the greatest flexibility and effectiveness to serve this purpose. The


-Phase 0
s=O


Phase 1









B-spline curves are composed of piece-wise Bezier curves between any two consecutive
markers. When we connect each Bezier curve piece by piece to form the global B-

spline curve, the original data points (here marker points) are junction points of two

consecutive Bezier curves. Now the question arises: how is the smoothness of this

overall B-spline curve? This question is answered by the definition of C' continuity

condition of the B-spline curves. The C' continuity means that at any junction point,

the nth derivative computed from the left Bezier curve equals to that computed from
right Bezier curve. The most popular B-spline curve is the C2 spline. To have C2

continuity, each Bezier curve must be at least degree three; therefore the cubic BMzier

curves can be used as piece-wise polynomial in this context. The curve fitting scheme

we adopted in this research is C2 cubic B-spline curve.

In the following, we describe the construction of C2 cubic B-spline curves and

so-called 'fairing' algorithm which is used to removed numerical noise to recover the

accurate information of the curvature of the interface.

Based on the spline theory (see Ahlberg et al., 1967; Farmin, 1997), to construct

a C2 cubic B-spline curve interpolating a set of data points Xo, ..., XL with corre-

sponding parameter values (or knots) Uo, ...UL, the vertices d, of the B-spline control

polygon have to be determined first, as shown in Figure 2.7. In this example, five
points xo, x5 are assigned initially. The corresponding knot sequence u0,.., 5
is chosen as the chord length at each point. The example here illustrates a closed

(periodic) curve, so x5 = x0. The relationship between the data points axi and the

control vertices dc is (see Farmin, 1997)


(Ai-1 + Ai) xi = aoidi-1 + idi + 7,idi+l


(2.1)












































Figure 2.7: Cubic interpolatory spline curve
control polygon, Bezier control polygon.


of five points along with its B-spline









where we have Aj = Au, = Ui+1 ui and


(Ai)2
a. = --
Ai-2 + Ai-1 + Aj
A A (=2 + A1) + A-1 (A + A +)
A*-2 + A,-1 + Ai Ai-1 + Ai + Ai+I


S A-1 + A + Ai+1

Using the periodic condition, we obtain a linear system of the form

3o 7o ao do ro
a,1 i1 71 di ri

OaL-2 13L-2 UL-2 dL-2 rL-2
7L-1 aL-l J -1 d L- rL-1

where the right-hand sides are of the form


(2.2)






(2.3)


ri = (A-1_ +- Ai) xi

The equation Eq. (2.3) can be solved using the procedure described in, e.g.,

Ahlberg et al. (1967), P. 15, yielding the vertices do, ..- d5 as shown in Figure 2.7.

Then the vertices of control polygon for piece-wise cubic B6zier curves can be

obtained based on this B-spline control polygon, namely, using the parameter values

to determine two vertices on each leg of the B-spline control polygon as shown in

Figure 2.7. Each of the piece-wise Bezier curve is defined by four control vertices,

including the two marker points, denoted by solid circles, and, between them, two

vertices to be decided from the B-spline control polygon. Therefore all four Bezier

curve control vertices are

b3i = Xi

b3i-2 = A-1 +- Ai di-1 + A-2d,
A 'A
b3i-1 A di-1 + Ai-2+ A- d (2.4)









where A = Ai_2 + Ai-1 + Ai. The index i denotes the ith vertice of B-spline control

polygon, and the total number of Bezier control vertices is 3i corresponding to vertice

i of B-spline control polygon for a closed curve.

Once we have the control vertices for each local Bezier curve, by defining a local

parameter 0 < t < 1 for the interval [ui,ui+i] as t = u-ui- we can express the

piecewise cubic Bezier curve in the form of
3
b(t) = b1Bf(t) (2.5)
i=O
where Bernstein polynomial B3(t) is defined by

Bi(t)=( t)-i3 (2.6)

with binomial coefficient

3i\ f 3 i f 0 < i < 3
3 = J !(3-i)! u-: -:
i 0 else

Using Eq. (2.5), the geometric information such as normal and curvature etc.

can be evaluated using analytical formulas as

Yt
nx (X t (2- )1/2
Xt
=y =-(X2

= (x + y)3/2
Xttt-yt

For axisymmetric geometries, the total curvature is the sum of n in Eq. (2.7)

and the other principle curvature, i.e., xt/y(x + y )3/2.

2.2.2 Fairing Algorithm

From the numerical simulation point of view, the smoothness of the interfacial

curve and the smoothness of curvatures are two most important characteristics of

the interface. Using the B-spline fitting, while a curve looks apparently smooth, the

curvatures obtained by Eq. (2.7) can be contaminated by numerical noises. The




























Figure 2.8: A sample interfacial curve from one simulation case


curvature formula involves the second derivatives as well as nonlinear products of
the first derivatives, and is prone to suffer from numerical noises even if the interface

shape and slope are found to be smooth. This usually occurs when the interface

is evolving in the course of computation. As an example to illustrate the problem,
Figure 2.8 shows an interface curve taken from a simulation conducted in this study.

The corresponding curvature plot based on B-spline fitting process is shown in Figure

2.9. There are substantial oscillations in the curvature profile, indicating that noises

associated with the numerical procedures are substantial. Such phenomena are well

known in computer-aided geometric design (see Farmin, 1997).
To treat this difficulty, the so-called curve fairing (smoothing) algorithms have

been developed in the literature (see Sapidis & Farmin, 1990; Farmin, 1997). One such

algorithm, in the C2 cubic spline case, is to make local spline curve segment three

times differentiable which is one order higher than that required by the original

cubic spline. The basic idea is to adjust the vertex locations to make the local






























Figure 2.9: Curvature plot obtained from B-spline fitting for the interfacial curve in
Figure 2.8.


curve segments around them become three times differentiable. Since the curvature

involves only first and second derivatives, this extra differentiable requirement will

ensure that the curvature be differentiable, and prevents discontinuity in the slope

of the curvature, as shown in Figure 2.10. The formula for obtaining new B-spline

control vertex dj are briefly listed as

S= (uj+2 Uj)lj + (uj Uj-2)rj(2.8)
am (2.8)


-j+2 j-2
where the auxiliary points Ij, rj are given by

I = (UJ+l -3)dj-l (Uj+, Uj)dj-2
Uj Uj_3
(uj+3 u l-)dj+l (uj j-l)dj+2


3 ------ \^ ".f)
.7 Uj+3 -Uj k
The geometry interpretation is illustrated in Figure 2.10 which is from page 365 of

Farin (1997). The detailed discussions can be found in Sapidis & Farin (1990); Farin

(1997).


10
8
6
4





-4

-6
-8
-1 0 *, ,i, i 1
0 40 80 120 160
Marker Point


/o n






















Figure 2.10: Knot removal: if dj is moved to dj, the new curve, denoted by dotted
line, is three times differentiable at uj [page 365, Farmin (1997)].


In practice, typically one fairing operation is not sufficient, and the fairing pro-

cedure will be repeated multiple times. The criterion adopted in the present study for

determining the convergence of the curvature profile is that the maximum difference

of the curvature values among all marker points between two consecutive fairing op-
erations must be less than 10-. For the curve in Figure 2.8, the resulting curvature

plot after 100 iterations of fairing operation is shown in Figure 2.11. The maximum

correction in the curvature is within 10-.

It is noted that the fairing algorithm is a geometric operation. It obtains correct

curvatures by modifying the geometry of the interface itself, not by manipulating

the formulas for computing the curvature of a given geometry. On the other hand,

attention needs to be paid to ensure that the interfacial marker locations, as well as

the volume/surface enclosed by the interface, be satisfactorily preserved. A critical

criterion of developing a satisfactory fairing algorithm is that with arbitrary number

of fairing iterations, the geometric information can be maintained at an asymptot-
ically constant state without being continuously smeared. The corrections in the
interfacial locations are usually very small, so the corrected marker locations and the






48


10 curvature before fairing
Curvature after fairing
8

6
4


Cc00
-2

-4
-6
-8
-10 I I I
0 40 80 120 160
Marker Point


Figure 2.11: Curvature plot obtained before and after fairing operation for the curve
in Figure 2.8.


resulted interface shape will not have perceptible changes after fairing operation, as

will be demonstrated in the testing of geometric algorithms.

2.2.3 Validation of Geometric Algorithms

The geometric algorithms for interface representing and tracking described in

the section 2.2.1 and 2.2.2 are capable of handling closed and open interfaces. They

are validated here for accuracy on some cases used in Shyy (1994).

First we examine the expansion of a circle on a 162 x 162 grid as in Figure 2.12.

Initially, the circle has a radius of 1, and expands at a constant speed. The curvature

of the circle at any given time is a constant. Figure 2.12(a) shows the interfacial

shapes at equal intervals of time. Figure 2.12(b) shows the curvature computed

by the aforementioned formulas, at the corresponding time instants as in Figure

2.12(a). As shown there, the interfacial curvature is a constant while, by regulating

the spacing between two neighboring markers, the number of markers increases as the








circle expands. The results shown in Figure 2.12 are obtained from B-spline fitting

algorithm only. No fairing operation is needed for this case; by applying fairing, no

impact is observed, either.
As another example for variable curvatures, Figure 2.13 shows an initial interfa-
cial curve shape expressed in the form of Y, = 2.0+(1-cos(27riXj/10)). The curvature
plot is found to be smooth using B-spline fitting for the initial shape as shown in
Figure 2.13(b). Given the normal velocity of each marker point as Vni = Yi Y1, the

interface moves as depicted in Figure 2.16(a). Figure 2.14 shows the interface shape
and corresponding curvatures at the early stage of the interface evolution. If the

fairing is not applied, the curvature obtained from B-spin representation, in Figure
2.14(b), exhibits oscillations near both ends of the curve after the interface evolves

to a new shape. After 10 fairing operations, the shape and curvatures at the same
time instant are shown in Figure 2.15. The curvature profile in Figure 2.15(b) is free

from the artificial spikes at the ends.
The curvature plot in Figure 2.15(b) becomes smooth and the shape in Figure
2.15(a) is well preserved and can't be distinguished from that in Figure 2.14(a) as if

the shape is identical to previous one. This is further confirmed by plots in Figure
2.18. Figure 2.18 shows the corrections on coordinates of markers by the fairing
operation at various time instants from the beginning to the end of the development
of the interface. All corrections of coordinates are small, which is the reason why
the two shapes before and after the fairing operations are virtually identical to each

other. In Figure 2.16 and Figure 2.17, the shape and the curvatures at the end of the
development of the interface are plotted while the former shows the situation before
fairing is applied, the latter illustrates the effect after the fairing. Furthermore, with a
periodic curve, the curvature should be periodic as well. As shown in Figure 2.13(b),
initially, the computed curvature is periodic; however, without fairing, Figure 2.16(b)









(a)
10


8


6


4


2


0
0 2 4 6 8 10
X
(b)
0.8

0.7 -

0.6 -
2 ______
S0.5 -

0.4 -__-
0 .3-
0.3 ___________________

0.2

0.1 '- I I I
0 100 200 300 400 500

Marker Point

Figure 2.12: (a) Shapes of expanding circle at equal intervals of time. (b) Corre-
sponding curvatures along the curve at the same time instants as in (a).









shows that the computed curvature for a moved interface is no longer periodic. After

fairing, as shown in Figure 2.17(b), the periodicity is restored again.

The shape is well preserved and the oscillation in curvature is removed. The

majority of the curvature plot after fairing is the same as that before fairing. The

number of marker points increases as a result of reorganizing process to follow the

increasing arc length of the bulging out interface.

The B-spline curve fitting in conjunction with the fairing operation is a robust

interface representing and tracking method which can yield smooth interface and

smooth curvatures while preserve the interface shape. The programming task, how-

ever, is quite challenging. For closed and open curves, different program modules and

algorithms have to be implemented for B-spline fitting mainly because the end con-

ditions are different. For closed curves, the end condition is always periodic (cyclic).

In the cases of open curves, there might be other end conditions.

2.3 Cartesian Grid Method for Sharp Interfaces

Once the interface is defined, one needs to identify in which phase each compu-

tational cell lies so that correct transport properties can be assigned. Furthermore,

for any two neighboring cells between which a interface passes through, there may be

a discontinuity in transport properties such as density and viscosity. Those boundary

cells are reshaped to maintain flux conservation around the interface. A major goal

of the present approach is to adopt a finite-volume formulation for all computational

cells so that mass, momentum, and energy conservation is honored in all resolvable

scale, including the computational cells intersecting with phase boundaries. Once a
cell intersects with an interface, it is split into two parts; with the partial cell con-

taining the center of the original Cartesian cell maintaining the initial cell index, and

the other merged into a neighboring cell belong to the same phase. As illustrated

in Fig. 2.19, this procedure results in irregular, trapezoidal shaped cells around the

interfaces. Away from the interfaces, most cells are still structured Cartesian cells.









(a)


0 2 4 6 8 10
X


(b)


Marker Point


Figure 2.13: (a) The initial shape
two fingers using B-spline fitting.


of two fingers. (b) The corresponding curvature of


0.5

0

-0.5

-1

-1.5

-2










10


8


6


4


2



0




1.5

1

0.5

0

-0.5

-1


-1.5

-2

-2.5


(a)


2 4 6 8 10

X

(b)


80 160

Marker Point


240


Figure 2.14: (a) The shape of two fingers at the early stage of the development when
fairing is not applied. (b) The curvature corresponding to the shape denoted by solid
line in (a).


I I I





54
(a)

10


8


6


4
>.




2


0 2 4 6 8 10
X
(b)
1.5

1

0.5



( -0.5


-1.5

-2 10 iterations of fairing

-2.5 1 I I
0 80 160 240
Marker Point

Figure 2.15: (a) The shape of two fingers at the early stage of the development when
fairing is applied. (b) The curvature corresponding to the shape denoted by solid line
in (a).




























0 2 4 6 8 10


(b)


0 80 160 240 320

Marker Point


400


480


Figure 2.16: (a) The shape evolution of two fingers until the later stage of the de-
velopment when fairing is not applied. (b) The curvature corresponding to the final
shape denoted by solid line in (a).


I I I I I


m









(a)


0 2 4 6 8 10


(b)


0 80 160 240 320


400 480


Marker Point

Figure 2.17: (a) The shape evolution of two fingers until the later stage of the devel-
opment when fairing is applied. (b) The curvature corresponding to the final shape
denoted by solid line in (a).


10 iterations of fairing
l I l l l


I I I I l I"


II









(a)


I I I I I I I
0 80 160 240 320 400 480


Marker Point


_____ ~(b)______


L.jdf


L 6


L -.A 4


0 80 160 240 320 400 480
Marker Point
Figure 2.18: (a) The corrections of the interface location by fairing algorithm at
time instants corresponding to shapes shown in 2.17(a). (b) The corrections of the
interface location by fairing algorithm at the same markers as in (a).


4.0E-03

3.0E-03

2.0E-03

1.0E-03

0.OE+00


5.0E-03

4.0E-03

3.0E-03

2.0E-03

1.0E-03

0.OE+00


, j J J





























Figure 2.19: Illustration of the resulting situation when cut-cell approach is applied,
i.e. fragments of cells which are cut by the interface are absorbed into neighboring
cells. The newly formed cells are shown in dashed lines on both sides of the interface


In such an algorithm, while a nominal structured grid index system is maintained,

the flux calculation across the surface is conducted based on interpolation of varying

degrees of polynomials. The resulting flow solver needs to account for both Cartesian

and trapezoid cells, using a finite volume, fractional step method.

After the domain is partitioned into several sub-domains representing different
phases. The so-called Cartesian grid method is developed to deal with the discretiza-

tion and solution of governing equations in such a domain left with mostly structured

rectangular cells and some unstructured trapezoidal shape cells in the vicinity of the
interface. Essentially, what the Cartesian grid method does is to solve the governing

equations on a domain with fixed, complex geometry interfaces without resorting to

the body-fitted grid. Although the method is developed in this research as one of

the major components in the entire moving interface simulation package, the method









alone is an good alternative to the body-fitted grid approach to solve the flow over

solid bodies with complex geometry, as demonstrated in Ye et al. (1999).

This globally second-order accurate Cartesian grid flow solver for fixed, immersed

interfaces, (see Ye et al., 1999), is summarized in the following:

1. The interface is represented numerically using C2 cubic B-spline interpolatory

curve;

2. A fractional step procedure (see Chorin, 1968; Kim & Moin, 1985) is employed

to solve the coupled governing equations, which results in a efficient solution of

time accurate, unsteady flows;

3. Colocated variable arrangement (see Ferziger & Peric, 1996) is adopted, which

leads to simplicity in coping with cut cells near the interface,

4. Second-order accurate central difference scheme is used for spatial discretiza-

tion, Crank-Nicolson scheme for temporal discretization;

5. A compact interpolation scheme near the immersed interfaces is devised that

allows us to retain second-order accuracy and conservation property of the

solver.

6. A preconditioned conjugate gradient iteration method is used for solving the

Poisson equation for pressure, which takes advantage of the underlying struc-

tured nature of the mesh and also substantially accelerates the convergence of

the Poisson equation for pressure.

2.3.1 Fractional Step Method for Unsteady Transient Solution

To use finite volume discretization, the dimensionless governing equations listed

in the section 1.1.7 on page 9 have to be written in the integral forms by integrating

each term over the control volume cv and wherever possible, using Gauss theorem








to convert volume integrals into surface integrals. To solve the momentum equations
(1.18) and (1.21) under the constraints (1.17) and (1.20) in the liquid and vapor
phases respectively, a so called fractional step procedure is applied. After the mo-
mentum equations are solved, the energy equations (1.19) and (1.22) can then be
solved in the updated flow field.
All governing equations are discretized on a Cartesian mesh using a colocated
arrangement (see Ferziger & Peric, 1996) of the primitive variables which are located
at the cell center. For convenience of describing the numerical procedure, only the
integral form of dimensionless governing equations in the liquid phase are considered
throughout this chapter as examples, the vapor phase transport equations are solved
in the exactly same procedure except additional dimensionless fluid property ratios
are present on those cells belong to the vapor phase.
The integral form of Navier-Stokes equations in the liquid phase are:
Mass conservation:

Su. ndS=0 (2.10)
Ca
Momentum conservation:
f f f i
SdV+ fu(u .n)dS=- pndS+ -] Vu-ndS (2.11)
Ci, CS CS CS
In the above equations cv and cs denote the control volume and control surface
respectively and n is a unit vector normal to the control surface. The above equations
are to be solved with u(x, t) = v(x, t) on the boundary of the flow domain where v
is the prescribed boundary velocity.
A second order accurate, two step fractional step method (see Chorin, 1968; Kim
& Moin, 1985; Zang et at, 1994) is used for advancing the solution in time. In this
time stepping scheme, the solution is advanced from time step 'n' to 'n+1' through
an intermediate advection diffusion step where the momentum equations without the








pressure gradient terms are first advanced in time. The semi-discrete form of this
advection diffusion equation for each cell shown in Figure 2.20 can be written as:

fdV = -J[3Un (Un.n)- Un1(Un'1-n)]dS
CV CS
+1 J (Vu* + Vu') .n dS (2.12)
Mne
CS
where u* is the intermediate cell-center velocity. A second-order Adams-Bashforth
scheme is employed for the convective terms and the diffusion terms are discretized
using the implicit Crank-Nicolson scheme. This eliminates the viscous stability con-
straint which can be quite restrictive in the simulation of viscous flows.
The velocity boundary condition imposed at this intermediate step corresponds
to that at the end of the full time step, i.e., u* = v"+'1.
At this stage, in addition to the cell center velocities which are denoted by u,
we also introduce face center velocities U. In a manner similar to a fully staggered
arrangement, only the component normal to the cell face is computed and stored
(see Figure 2.20). The face center velocity is used for computing the surface flux
from each cell in our finite volume discretization scheme. The reason of separately
computing the face center velocities will be addressed later in this section.
Following the advection diffusion step, the intermediate face center velocity U*
is computed by interpolating from the cell center intermediate velocities.
The advection diffusion step is followed by the pressure correction step
Un+l -- U*
= = VP (2.13)
At

where we require that the final velocity field satisfy the integral mass conservation
equation given by

f (U+1n) dS = 0 (2.14)
CS






62


I I
I I I


'' [/




4j- -- --
I I
I N














0 5-point stencil for regular cell


Figure 2.20: Schematic of underlying Cartesian mesh and arrangement of cell-center
and face-center velocities


This results in the following equation for pressure


I (Vp) -n dS= f (U n)dS (2.15)
C3 C

which is the integral version of the Poisson equation for pressure. Note that the

pressure correction step is represented by the inviscid equation (2.13) and is well

posed only if the velocity component normal to the boundary is specified. Therefore

the velocity boundary condition consistent with Eq. (2.13) is ul+.n = V1f1l where

n is the unit normal to the boundary of the flow domain. It can be easily shown

that this implies that (ApnIl) -n = 0 be used as the boundary condition for equation

(2.15). Once the pressure is obtained by solving this equation, both the cell center

and face center velocities are updated separately as follows:


u 1= At (Vpe)" ; = At (Vep+')fc (2.16)









It is a well known fact that in a colocated mesh, when the compact stencil is used

for the Poisson equation for pressure, the pressure field usually does not exhibit non-

physical oscillation, however, the final velocity field does not satisfy the divergence

free constraint exactly. To overcome this problem, two types of approaches could be

used, i.e., the momentum interpolation (see Rhie & Chow, 1983) scheme in which a
pressure gradient term is introduced into the interpolation process for obtaining cell

face velocities. As it is shown by Ferziger & Peric (1996), this technique produces a
third-order dissipation in estimating the cell face velocities. The momentum interpo-

lation technique is not straightforward to implement in cut cells encountered in the

mixed Eulerian-Lagrangian algorithm for flow with immersed interfaces.

The other approach for dealing with this problem was first used by Zang et al.

(1994). In their method, both cell center and face center velocities are stored while

momentum equations and Poisson equation for pressure are solved at the cell centers

only. After Poisson equation for pressure is solved, both cell center and face center

velocities are updated. The pressure gradient at the cell face is computed in a compact

stencil manner. They have used this procedure in conjunction with a curvilinear mesh

solver to simulate the turbulent flows with LES. In our algorithm, we implemented
this method and test results show that this approach is well suited for high Reynolds

number, unsteady flows.

2.3.2 Discretization on Regular and Cut Cells

To discretize governing equations in the configuration as shown in Figure 2.19,
we have to deal with two types of situations: discretization on cells which are away

from the interface thus not cut by the interface and on those cells which are cut by

the interface and reshaped to become trapezoidal cells.
The numerical procedure for solving Navier-Stokes equations just outlined can

be easily implemented on a regular Cartesian grid cell. In this section, detailed

finite volume integration of each term is described. For demonstration purpose, basic





























Figure 2.21: Schematic of cells around the interface. (a) boundary cells with interface
located south of cell center. (b) boundary cells with interface located west of cell
center. (c) typical reshaped quadrilateral cells corresponding to case (a). (d) typical
reshaped quadrilateral cells corresponding to case (b).


procedure for evaluating volume and surface integrals on a regular cell is explained
first, followed by the special treatment on those cells cut by the immersed interfaces
and reshaped, as shown in Figure 2.19. The emphasis is placed on the latter.
In general, the so called cut cell approach is used to treat those reshaped cells:
the small fragments of cut cells are absorbed into neighboring cells in the same phase
to form mainly trapezoidal cells as shown in Figure 2.21, after that, the governing
equations are integrated on those reshaped cells.
It is straightforward to obtain second order accurate spatial discretization on
regular cells. The key for retaining global second-order accuracy of the solver is to
evaluate mass, convective and diffusive fluxes and pressure gradients to second order
accuracy on the cell face of those trapezoidal cells resulted from reshaping process.


(a)











(c)








As shown in equation (2.12), a finite volume discretization of Navier-Stokes

equations requires the estimation of surface integrals on the faces of each cell. The

integrand (denoted here by f can either involve the value or the normal derivative of

a variable. An example of the former is the convective flux denoted by (pev-n) and
of the latter, the diffusive flux given by (LFVq.n) where q0 is a generic scalar variable

and F is diffusive coefficient. In addition to this, the Poisson equation for pressure
also requires evaluation of the normal pressure gradient. In order to estimate these

surface integrals with second order accuracy, the midpoint rule can be used and this

requires an accurate evaluation of the integrand at the center of the face. For regular

cells which are away from the interface the integrand can be evaluated at the face

center to second order accuracy in a straightforward manner by assuming a linear

profile between nodes on either side of the face.

However, this is not the case for the trapezoidal cells. Unlike regular Cartesian

cells, for those trapezoidal cells near the interface, simple linear approximation be-
tween neighboring cell centers can not be made to calculate the fluxes and pressure

gradients on the face centers to second order accuracy. This is because the centers

of some of the faces of such a cell (marked by a shaded arrow in Figure 2.21c and

d) may not lie in a location which puts it in the middle of neighboring cell centers
where a linear approximation would produce second order accurate estimate of the

gradients. Furthermore, some of the neighboring cell centers do not even lie on the

same side of the immersed interface and therefore cannot be used in the differencing
procedure. Thus, not only do we need a procedure for computing these face center
quantities accurately, but we also require that the procedure adopted be capable of

systematically handling reshaped boundary cells with a wide range of shapes.

Our solution has been to use a compact two dimensional polynomial interpo-
lating function which allows us to obtain a second order accurate approximation of
the fluxes and gradients on the faces of the trapezoidal boundary cells from available























Figure 2.22: Schematic of interpolation for cell face values and derivatives at trape-
zoidal cells (a) various fluxes required for trapezoidal cell (b) trapezoidal region and
stencil used in computing faw.


neighboring cell center values. The current interpolation scheme coupled with the
finite volume formulation guarantees that the accuracy and conservation property of
the underlying algorithm are retained even in the presence of curved immersed inter-
faces. In the following, we describe the interpolation function for a typical trapezoidal
boundary cell.
Consider the trapezoidal cell ABCD& in Figure 2.22a. The face ABC of the
trapezoidal cell is composed of two pieces: AB coming from cell P and BC coming
from cell S. The integral on this face can be decomposed as

I f dy = Jf dy + Jf dy (2.17)
AC AB BC
A second order approximation to this integral can then be obtained as

f f dy f (yA yB) + fsW (YB yc) (2.18)
AC
where fand f,, are computed at the centers of segments AB and BC respectively. On
the other hand, if the face is cut by the immersed interface such that it is smaller than


(b)

I I







0 Points in the six-point stencil for fsw








a nominal cell face, as in the case of face TD then the integral can be approximated
as

f fdy fe(Ye -yv) (2.19)
Ve
where fe is the flux computed at the center of the segment VS. For non-boundary
cells, these face center values can be evaluated to second order accuracy quite easily
by a linear approximation and we would therefore like to evaluate f", fs'., fe to within
second order accuracy also. Approximation of fw to second order accuracy is quite
straightforward and is done in the same way as for the face of a non-boundary cell.
For instance, if fw, requires the value of 0, this can be evaluated to second order
accuracy as

,w = OwAw + Op (1 Aw) (2.20)

where the linear interpolation factor Aw is defined as

w = X (2.21)
Xp XW

Alternatively, if f, requires the normal gradient of q as it would for the diffusion or
pressure gradient terms, this can be approximated by a central difference scheme as
follows:
(0 _p-Ow (2.22)
aXIw Xp -XW

This approximation is second order accurate when the cell face is midway between
P and W, i.e., when the mesh is uniform. However, expressions such as Eq. (2.20)
or Eq. (2.22) cannot be used in evaluation of fsw or fe to second order accuracy
since in many instances some of the neighboring nodes can be inside the immersed
interface. For instance, for the situation shown in Figure 2.22a, the south node is
inside the immersed interface and cannot be used in the evaluation of fs.. Even
if neighboring nodes are available, as they are for the east face, it is not clear how








a second order accurate scheme can be constructed since fe is not located on the
line joining the neighboring cell centers and consequently, expressions such as Eq.

(2.20) or Eq. (2.22) cannot approximate this flux to second order accuracy. Thus, a
different approach is needed here for evaluating these fluxes.

Our approach is to express the flow variables in terms of a two dimensional
polynomial interpolating function in an appropriate region and evaluate the fluxes
such as fl. or fe based on this interpolating function. For instance, in order to

approximate fsw,, we express q in the shaded trapezoidal region shown in Figure
2.22b in terms of a function that is linear in x and and quadratic in y

0 = Clxy2 + c2y2 + c3xy + c4y + Csx + C6 (2.23)

where Cl to c6 are six unknown coefficients. If fsw involves the normal derivative of

, this can be obtained by differentiating the interpolating function, i.e.

ao = cly2 + c3y + CS (2.24)
y-=CiY +C3y/+C5 (.4

The rationale for choosing Eq. (2.23) as the interpolating function for evaluating

fs,,, is as follows: the objective here is to evaluate (94//ax) at the center of BC to
within at least second order accuracy. Furthermore, we would like to do this with

the most compact interpolating function so as to minimize the size of the stencil
required for the boundary cell. Clearly, a biquadratic interpolating function in the

trapezoid shown in Figure 2.22b would lead to a second order accurate evaluation of
the derivative anywhere inside the trapezoid. However, a biquadratic function has
nine unknown coefficients and therefore requires a large nine point stencil. It turns
out however that for the trapezoid shown in Figure 2.22b, second order accurate
evaluation of the derivative on the cell face can be achieved by using an interpolating

function that is quadratic in y but only linear in x. This is because BC is midway
between the two parallel sides of the trapezoidal and in a manner analogous to central
differencing, linear interpolation in the x-direction leads to a second order accurate








evaluation of derivative at this location. On the other hand, this situation does not

exist in the y-direction for the cell shown in Figure 2.22b and therefore a quadratic
interpolation is necessary in this direction in order to obtain a second-order accurate

approximation to (,9/9x) at the center of BC. In Appendix 1, we have demonstrated

numerically that the linear quadratic interpolating function shown in Eq. (2.23) does
indeed result in second order accurate evaluation of values and derivatives on a line
which is located midway between the two parallel sides of a trapezoid.

It can be seen in Figure 2.22b that the sides of the trapezoid in which the
interpolation is performed pass through four nodal points and two boundary points.
Thus, the six unknown coefficients in 2.23 can be expressed in terms of the values of

4 at these six locations. To solve for c, we obtain the following system of equation

by expressing the 0 at the six location in terms of the linear quadratic interpolating

function:
^ 'Xlyl y1 x-^y- y1 Xi1 c
{ 0 X2Y2 Y2 y x2Y2 Y2 X2 1 c2 (2.25)
06 X6yj y2 x6y6 Y6 x6 1 c6

The coefficients can now be expressed in terms of values of 0 at the six points by

inverting Eq. (2.25), i.e.
6
cn = bnjj, n = 1,6 (2.26)
j=1
where bnj are the elements of the inverse of the matrix in Eq. (2.25).

After cq is obtained, the value of 0 at the center of BC is expressed in the form
of

0sw = CiXswyS + y2 + C3XawYsw + C4Ysw + C5Xsw + C6 (2.27)

and using Eq. (2.26) this can be rewritten as
6
03w =Ea.jj0 (2.28)
j=1








where

aj = bljxsw,,yw + bj3yw + b3jxswysw + b4jysw + b5sjx,w + b6j (2.29)

The value of (9qp/9x) at the center of BC is expressed as

c( ) + C3yw + C5 (2.30)
lox = W

and using Eq. (2.26) this can be written as

=0 (2.31)
sw j=1
where

O = biy' + b3,Ysw + b5j (2.32)

Note that a and / are coefficients that depend only on the mesh, the location and
orientation of the immersed interface. Therefore these can be computed once at the

beginning of the solution procedure. Subsequently, relationships such as (2.28) and
(2.31) can be used in the spatial discretization of the governing equations (2.12)-

(2.16).
A similar interpolation procedure is also used for approximating fe. For this,

a linear quadratic interpolating function is used in the trapezoidal region shown in
Figure 2.23 and a relationship similar to Eq. (2.28) and (2.31) is developed for ap-
proximating fe. The six points contained in this stencil are also shown in Figure
2.23. It should be pointed out that the north face of the particular cell being con-

sidered here does not need special treatment since face center values and derivatives
can be computed to second order accuracy using a linear approximation. However,
in general there are also boundary cells which have their north or south faces cut by
the immersed interface (as shown in Figure 2.21b). For these boundary cells too, the
same approach is used to evaluate the fluxes on these cut faces. The only difference
here is that the interpolating function is linear in y and quadratic in x.























0 Points in the six-point stencil for fe


Figure 2.23: Trapezoidal region and stencil used in computing fe.


Now we turn to the calculation of the flux on cell face CD which lies on the
immersed interface as shown in Figure 2.22a. The integrated flux on this face can
again be evaluated to second order accuracy using the midpoint rule and as before
we would like to evaluate the integrand at the center of face CD (denoted here by

lt) to second order accuracy. In general both convective and diffusive fluxes are
needed on this face and this requires approximation of variable values as well normal
derivatives at the center of CD. The value is usually available from a Dirichlet type
boundary condition and hence no interpolation is required for this. Here we describe
the approximation procedure for the normal derivative. The normal derivative on
face CD can be decomposed as
o oq
= xn, + (2.33)

where nx and ny are the two components of the unit vector normal to face CD.
Since we know the shape of the immersed interface, n, and ny are known. Therefore
computation of the normal flux requires estimation of 9Ob/9x and 9O/9y at the center
























0 Points in the stencil.


Figure 2.24: Stencil for calculating interface flux. (a) stencil for calculating 00/19y
(b) stencil for calculating 9ck/Ox


of the line segment CD. For the cell being considered here, Q90/9y is computed to
second order accuracy with relative ease by expressing the variation along the

vertical line in terms of a quadratic in y as follows

S= aiy2 + a2y + a3 (2.34)

The coefficients in the quadratic can be expressed in terms of the values of 0 at the
three points indicated in Figure 2.24a. Subsequently, the normal derivative at the
center of face CD is evaluated as:

( = 2aiyint + a2 = r (2.35)
int j=1

where again rI are coefficients which depend solely on the geometry of the boundary
cell.
Unlike the calculation of 9o/Oy for this cell, the calculation of O9/Ox is not
straightforward. However, an approach consistent with the computation of fsw and

fe can be used to estimate the value of this derivative to desired accuracy. Consider









the trapezoid shown in Figure 2.24b. Again, because the x-coordinate of the center

of CD) is located midway between the two parallel sides of this trapezoid, expressing

the variable in this trapezoid in terms of an interpolating function which is linear in
x and quadratic in y allows us to obtain a second order accurate approximation to

(09/Ox)M at the center of the line segment C). The procedure for this follows along
lines similar to that shown earlier for (q/O/x)8 and we get the following expression

for the x-derivative on the interface:
6

(fi') mX0 (2.36)
ax int j=1

where rj depend on the location and orientation of the immersed interface in the

neighborhood of the cell under consideration. Finally we obtain an expression of the
form


(ar j j (2.37)
ic'7nt j=1

for the normal gradient where Trj can be obtained from Eq. (2.33), (2.35) and (2.36).
Thus we obtain a nine-point stencil for the flux on the interface and the points in this
stencil are shown in Figures 2.24a and b. As can be seen from these figures, of these
nine points, three points lie on the immersed interface and their values are available

from the prescribed boundary condition.

It should be pointed out that although most cells are four-sided trapezoidal cells,
some five-sided cells and three-sided triangular cells may also be encountered. How-
ever, the discretization of the governing equations for these cells can also be handled
within the framework of the current interpolation scheme. With all of these features,
the current solver can in principle, handle arbitrarily complex geometries. Further-
more, multiple immersed interfaces can be handled as easily as a single interface.
This is in contrast to body fitted grid where the grid topology can get quite compli-
cated in the presence of multiple interfaces. Finally, since both sides of the immersed








interface are covered by the grid, we have the capability to solve a different set of
equations on each side of the immersed interface. For instance, the multi-phase flow
with separate governing equations in different phases across the interface could be
simulated in this manner.

2.3.3 Imposition of Boundary Conditions on the Interface

This section describes the imposition of boundary conditions on curved, immersed
interfaces on the underlying Cartesian grid
In the discretization process on cut cells as discussed in the last section, the
interpolation scheme used to obtain the face flux such as fe, fsw and fint as shown in
Figure 2.22 involves the values of generic variable 0 at point (5), (6) on the interface
as illustrated in Figure 2.22, 2.23 and 2.24. The values of on the interface are
simply known if Dirichlet boundary conditions are specified on the interface. In the
case of Neumann boundary conditions on the interface, the explicit values of 0 on
the interface need to be extrapolated from values in the interior cells.
The extrapolation process is explained in Figure 2.25. For implementing Neu-
mann boundary conditions on the interface, it is necessary to obtain the expression
for gradients 0 on the interface. Taking one side of the interface as example as shown
in Figure 2.25, the gradient can be evaluated by using a normal probe, i.e., drawing a
normal into one phase from the interfacial marker point location. Since second order
accuracy is desired, two nodes on the normal probe are needed which are located at
distance h and 2h from the interface. Values at each node, Oni, 0.2 are obtained by
bilinear interpolation from the neighboring grid points which are denoted by shaded
circles. Then based on the values at two points on the normal probe and the known
gradient values on the interface, an 0(h2) estimate of the normal gradient at the
marker point k location can be written as:

(0) 4ni n2 3 (2.38)
k9n 2h

























Figure 2.25: Illustration of implementing Neumann boundary conditions on the in-
terface.


Suppose the Neumann boundary condition on the interface gives (o) k = 0, the value
of 0 at marker point k can be extrapolated based on this expression for gradient as
40,1 0,2 0.2h (2.39)
3

At each marker point, the value is extrapolated in the same manner. All values 4
at sequential marker points are parameterized in terms of arc length s as

4(s) = aos2 + bos + co (2.40)

In the cut cell procedure described earlier, however, the values of at the inter-
sections of the interface with grid lines are needed, e.g., 05,6 in the Figure 2.25.
Let (xk,yk) be the coordinate of this interfacial marker point. Previously the
functions x(s) and y(s) on the interface for each marker point has been obtained.
Thus, the arclength value S5,6 along the interface where the curve intersects the line
x = xi is obtained by solving


(ax)ks2,6 + (b.)kS5,6 + (C.)k = Xi(


2nd normal probe
node On2



1st normal probe
node ni


(2.41)








Once s5,6 is obtained, based on the parameterization of 0 in terms of s, the value of

05,6 is

05,6 = (ao)ks,6 + (b)kss,6 + (CO)k (2.42)

These values can then be incorporated into the interpolation scheme to calculate

the flux on the cut cell faces in the discretization process.

2.3.4 Inversion of Discrete Operators

This section discusses our choice of algebraic equations solver.
The finite volume discretization of Eq. (2.12) or (2.15) in a given cell P can be

written as
M
S k = bp (2.43)
k=1
where X's denote the coefficients of the nodal values within a stencil of size M and

bp is the source term that contains the explicit terms as well as the terms involv-
ing boundary conditions. The term on the left-hand side represents a discretized
Helmholtz operator in the case of the advection-diffusion equation and a Laplacian

operator in the case of the Poisson equation for pressure. In the present method, M
is equal to five for non-boundary cells and this five-point stencil is shown in Figure
2.20. For the trapezoidal boundary cells, the linear quadratic interpolation scheme
described earlier results in a six-point stencil and for the particular boundary cell
considered in Figure 2.22, the resultant six point stencil is shown in Figure 2.26.

The discretized advection-diffusion equation and Poisson equation for pressure
result in a coupled system of linear algebraic equations which requires the inversion of
a large, sparse, banded matrix. The structure of this matrix is for the most part simi-
lar to that obtained on a Cartesian mesh without any immersed interfaces. However,
the presence of the immersed interface alters this matrix since the coefficients in Eq.
(2.43) are different for the trapezoidal boundary cells. Furthermore, the rows in the























0 Points in the six-point stencil


Figure 2.26: Six points stencil of finite volume integration on a typical boundary cell


matrix corresponding to the trapezoidal boundary cells also have additional elements
since the stencil of the trapezoidal boundary cell is different from regular cells. The

alternative direction line successive-overrelaxation (SOR) method (see Press et al.,
1992) is one of the most widely used iterative methods for solving equations resulting

from discretization on structured grids. A typical implementation of this technique,

which we have adopted here, involves alternating sweeps along the two families of
grid-lines. We find that even in the presence of the immersed interfaces, this method

is extremely effective for the numerical solution of the discretized advection-diffusion

equation and the residual can be reduced to an acceptable level within a few itera-

tions.
The discretized Poisson equation for pressure however exhibits a slower conver-

gence than the advection-diffusion equation. This is because the time-derivative term
in the advection-diffusion equation tends to improve the diagonal dominance of the
corresponding discretized operator. In fact due to its slow convergence, the solution









of the discretized Poisson equation for pressure is usually the most time consum-
ing part of a fractional-step algorithm. In the presence of immersed interfaces this

behavior of the Poisson equation for pressure can be further exacerbated since the

stencil for the trapezoidal cells contains dependence on some neighboring cells which

are not included in the line-SOR sweeps. For example, in Figure 2.26 the coupling of
cell P with cells (1) and (4) in the stencil is not included directly in any of the line

sweeps. Furthermore, depending on the aspect ratio of the trapezoidal boundary cell

and the angle at which the immersed interface cuts the cell, diagonal dominance in

the pressure operator can be severely weakened. In the various simulations that have
been performed using the current method, it has been found that the simple line-SOR

procedure can result in an extremely slow convergence for the pressure equation.

One remedy in such a situation is to resort to a multigrid method (see Brandt,

1977; Briggs, 1987). However, the presence of the immersed interface can complicate

the implementation of a multigrid procedure since operations such as prolongation
and restriction are difficult to perform near the immersed interface. In contrast,

Krylov subspace methods (see Golub & Loan, 1989) are an attractive alternative

since these are designed for general sparse matrices and therefore do not assume any

structure in the matrix operator. Thus, the presence of the immersed interface does

not pose any additional complication for the implementation of these methods. A

particularly suitable method in this class is the bi-conjugate gradient stabilized (Bi-

CGSTAB) method (see van der Vorst, 1992; Barrett et al., 1993) which is applicable

to non-symmetric matrices and provides relatively uniform convergence. However, as
with all conjugate gradient methods, the convergence rate of the Bi-CGSTAB proce-

dure depends critically on the choice of the preconditioner. The Jacobi preconditioner
which has a trivial construction phase is used routinely in conjunction with unstruc-

tured grids. However, this preconditioner only produces marginal improvement in

the convergence rate of conjugate gradient type algorithms. Other preconditioners









such as those based on incomplete factorization can substantially increase the con-

vergence rate but these usually have a non-trivial and expensive construction phase

(see Barrett et al., 1993).

The structure of the matrix that results from the current Cartesian grid ap-

proach however presents us with another alternative choice of a preconditioner. As

pointed out before, the presence of an immersed interface only alters the underlying

matrix operator in the rows corresponding to the boundary cells. Although this alter-

ation slows the convergence rate of the line-SOR procedure, the convergence rate is

still significantly faster than what is obtained with a simple point Jacobi method. It

follows that the line-SOR procedure would serve as a better preconditioner than the

Jacobi preconditioner. A further advantage of using the line-SOR as a preconditioner

is that this procedure only requires the solution of tridiagonal systems and this can be

accomplished with ease using the Thomas algorithm. Thus, in the solver developed

here, the line-SOR procedure was used as a preconditioner in the Bi-CGSTAB algo-

rithm and a significant improvement in the convergence rate over a simple line-SOR

iterative procedure was found.

This completes the discretization and solution procedure of Navier-Stokes equa-

tions in a domain with immersed, stationary interfaces. The same solution procedure

can be applied to solve energy equations considering it is an advection-diffusion equa-

tion similar to momentum equations.

When the immersed interfaces are no longer stationary due to various physical

mechanisms such as phase change or hydrodynamic forces effect, the situation be-
comes much more complicated. Many issues associated with solving moving interface

problems arise. If the discontinuity of fluid properties become appreciable and even

severe, the overall numerical procedure will face additional challenges.

In the following sections, key issues regarding solving transport equations for

problems with moving interfaces are addressed.









2.4 Moving Interface Algorithms

In the previous chapters, the numerical procedure used to solve the coupled,

nonlinear transportation equations in the domain with immersed fixed interfaces is

described. The method is fixed-grid approach in which the grid lines don't conform

to the interface with complex geometry.

When the immersed interfaces in the flow domain translate over the underlying

fixed Cartesian grid due to physical mechanisms, various computational issues related

to the moving interface can arise.

Overall, the numerical procedure involving immersed moving interfaces on a fixed

Cartesian grid is to solve a fixed interface problem at any time instant provided that

the instantaneous interface location is known. However, since the interface movement

is not independent of the flow field, this coupling between the interface evolution and

the flow field makes the moving interface problem not a simple addition of many

individual fixed interface problems because determining the interface location at any

time step is also part of the solution process. For fixed interface problems, the global

iteration process is only carried out among the coupled mass, momentum and energy

equations to solve for the flow field. With the introduction of moving interfaces, the

interface location is coupled with the flow field, another outer iteration is needed to

deal with this higher level of coupling. The major difficulty of the entire simulation

methodology arises in this iteration process to solve for the flow field and to adjust

the interface shape together.

Physically, the interface location is to be determined based on, in the case of

liquid-vapor two phase flow considered here, the thermal effects such as evaporation

or condensation of liquid or vapor as well as hydrodynamic forces such as pressure,

viscous stresses in the two phases. The interface must be advected to a position where

multiple interfacial conditions such as mass, momentum and energy conservation

conditions (see Equations 1.23, 1.24, 1.25), and mass conservation constraint within









each phase are simultaneously satisfied. In order to achieve this, a major numerical

technique is required to determine the interface location in a iterative manner within

any particular time step.

2.4.1 Change of Phase

As it is mentioned earlier in this chapter, when a interface which represents

the phase front of the liquid and vapor state moves across the underlying stationary

Cartesian grid, some grid points may undergo discontinuous change of the phase due
to the interface motion as illustrated in the Figure 2.27. At the time level n, the

interface is located such that the grid point (i,j) lies in the vapor phase as shown in

Figure 2.27a. At the next time level n+1, the interface moves to the updated location

shown in Figure 2.27b, as can be seen, some grid points which previously lie in the

vapor phase now lie in the liquid phase. In the time advancement solution procedure

for the transportation equations, the value of any variable 0 at previous time level n

is needed at current time level n + 1 by the temporal discretization scheme adopted.

However, for grid point that have just changed the phase since last time level n such

as (i, j) in the Figure 2.27, the current computation at time level n + 1 is performed
in the liquid phase while previous history for this point is in the vapor phase which

couldn't be used, i.e., there is no physically meaningful information regarding the

previous value of on at (i, j). This is a purely computational issue encountered in the

numerical approach that treats the interface as sharp discontinuity in terms of the

fluid properties. In the purely Eulerian methods, the interface is not explicitly tracked

but deduced from certain field variables which are smoothed over the fixed grid, thus
the interface is spread over a band of grid points rather a sharp discontinuity. In that

case, the problem of change of the phase is not a significant issue if it is a issue at
all. However, if the sharp discontinuity of the interface is desired to be maintained

like in our approach, those grid points have to be dealt with.








A solution to this problem is to estimate the value of any variable at the grid
point (i, j) in the Figure 2.27 by interpolation. At time level n + 1, the grid point
(i, j) falls into liquid phase, but this point was in the vapor phase at last time level n.
Therefore the value of Oi, doesn't exist in terms of the liquid phase and is obtained
through interpolation between the surrounding grid points in the liquid phase and
points on the interface.
The procedure of the interpolation is as following. First a normal to the interface,
which passes through the grid point (i, j) and has a length of the grid spacing h, is
drawn. The end of the normal is called /3, the intersection of the normal with the
interface is denoted as a. The value of is estimated by interpolation through
values at a and 03 as

=' ata + bop

where a, b are interpolation coefficients based on length ratios between a, # and (i, j).
0,, is the known boundary value on the interface, Op could be obtained by constructing
a bilinear interpolation within point 1, 2,3,4 illustrated in the Figure 2.27:

linj + 622 + 6303 + 6404 = -00

where the 61-4 are geometric coefficients in the bilinear interpolating function. Com-
bining these two expressions gives the estimated value at the grid point (i, j) as
Saq, + b (6202 + 6303 + 644)
23 =1 -- b-l

Note that if two adjacent grid points change the phase at the same time, the
above process is still applicable except that additional iteration process between those
adjacent grid points is needed to find their values simultaneously.

2.4.2 Interfacial Conditions

The interface shape in our simulation is not known as a priori, instead, it is
determined as a part of the entire solution process. At each time step, in order to
























Figure 2.27: Illustration of change of the phase on some grid points following the
interface movement. (a) Location of the interface at the time level n. The shaded
region designates the vapor phase. At this instant, the grid point (i, j) lies in the
vapor phase. (b) Location of the interface at the time level n + 1 after moving. The
grid point (i, j) now lies in the liquid phase. Those grid points which undergo the
change of the phase are marked by open circles. Also shown in (b) is the schematic
of finding the value at the (i, j) by drawing a normal probe passing through (i, j),
the value at the end of the normal probe is obtained by using a bilinear interpolating
function over points numbered 1 to 4.


determine the correct interface location which satisfies all interfacial conditions (1.23),
(1.24), (1.25) and (1.26), an iterative process for updating the interface location is
needed as described in section 2.4.3 and 2.4.4.
The interfacial temperature is specified using Eq. (1.26), the vapor pressure at
the interface is poo, then the liquid pressure at the interface can be determined by
Eq. (1.24).
The energy conservation condition Eq. (1.25) is used to estimate the interface
velocity relative to vapor velocity at the interface each time, thus it is naturally
satisfied. The mass conservation condition Eq. (1.23) serves for determining the
liquid phase velocity at the interface when the interface velocity and vapor phase
velocity at the interface are known. Both liquid and vapor phase velocities at the


Marker points (1


(a)


' a








interface are assigned as the boundary conditions respectively for solving the flow

fields inside and outside the vapor bubble.

2.4.3 Iterative Process for Determining the Interface Shape

In this section, we will describe the process for determining the interface shape

using normal stress balance coupled with flow solver via an iterative process. This

is the most important technique and difficult numerical procedure in this work. The

method depicted in the following is used to determine the bubble deformation caused

by the momentum balance with no phase change occurs. The technique implemented

in the present research is the one similar to that proposed by Ryskin & Leal (1984)

in their method using body-fitted moving grid. Here we implemented it in the fixed

grid which results in a more efficient overall solution procedure.

If phase change is to be considered and contributes, in conjunction with the mo-

mentum balance, to the shape change of the bubble, additional algorithm is required

and developed in the present study. They are given in the next section.

Given an initial interface shape and flow field, material properties are assigned

in the liquid and vapor phases, the flow field is then solved until a converged flow field

is obtained when the interface is fixed. This flow field serves as the starting point

for subsequent moving interface problems. Then the solution algorithm proceeds

iteratively through the following steps starting from iteration number k=0 for next

time step n + 1, the initial guess for flow variables at time step n + 1 is 0,+1,0 = On,

the initial interface velocity is uin1'0 = un.

The interface shape is determined by the normal-stress balance Eq. (1.24). The

technique proposed by Ryskin & Leal (1984) is adopted to obtain the estimate of

the interface shape for next iteration. The idea behind this technique is very simple.
This approach essentially takes the local imbalance between the total normal stress

Tn, which includes both static and dynamic pressure and viscous contributions, and








capillary forces
1
I(s) = Tn- We (2.44)

as a "driving force" which causes a local displacement of the interface in the normal
direction, the magnitude of the local displacement is proportional to H(s). Thus
n+l,k+l n+ l,k+ n 00 k
int = int+/H s)n
Yin ,k+ nX,k + OPI1(8) l
n+1,k+1 = n +n (2.45)
,-nt = +1t, /3Hk(s)

where the under-relaxation coefficient /3 has to be determined by numerical experi-
ment, its typical values being 0(10-4 to 10-3) in our computations.
In the cases where interface completely encloses one of the two fluids, e.g., bubble
in our study, the local increment in the location of interface marker points must be
done in such a way as to satisfy the mass conservation constraint, i.e. preserve the
bubble volume if density is constant. We know the change in the volume between
the kth and (k + 1)th iterations is

/ 1n+lk+1 _- tlk)2 + (n+lk+l lk)2]1/2ds (2.46)
[[ IXint -Xint ) Yint Yint ) ] (2.46)

where the integral is taken along the interface. Since
[ n+lk+l n+lk2 yn+lk+l n1lk-)2]1/2 k(s) (2.47)
(Xint --Xint ) Ynt -- nt (2.47

we have the mass constraint

foI k(s)ds = 0 (2.48)

This constraint determines a free constant of the pressure field in the total normal
stress Tn at each iteration k so Eq. (2.48) is satisfied.
Even after this constraint has been satisfied, the bubble may still change volume
slightly at each iteration due to higher-order numerical errors. To prevent these
small changes from accumulating and becoming significant, a simple scaling of the








interface can be done following the Eq. (2.45). The uniform scaling magnitude A of

the interface location in the normal direction can be determined from

f Alds= AV (2.49)


where AV is the error in the volume. This simple scaling is very effective, normally

one or two iterations of this process is sufficient to bring the percentage error in the

volume down to 10-5.
The global iterative process within each time step is therefore as follows.

1. For a given shape of the bubble, the flow field is computed by solving the
Navier-Stokes equations with a small number of iterations on pressure Poisson

equation.

2. Knowing the flow field, the normal stress balance at the interface is calculated
and checked. If it is not satisfied, the interface shape is modified according
to Eq. (2.45) so as to reduce the imbalance between the total stress and the
surface tension force.

3. After each interface update, the mass conservation is enforced by rescaling the

interface using Eq. (2.49).

4. The interface normal velocity is calculated from the kinematic condition u+l =
t int -

(_nt xnt)/At, which is the boundary conditions for solving the flow field.

5. Return to step 1 and repeat until all equations and boundary conditions are
satisfied to a predetermined level of accuracy.

2.4.4 Determining Interface Shape with Phase Change

The interface movement when phase change is present can be decomposed into
two components: local propagation velocity due to evaporation or condensation, and
the global movement resulted from buoyancy when there is no phase change.








The interface velocity component due to phase change effect in the normal di-
rection, which is also the velocity component relative to the vapor phase, is as Eq.
(1.25), i.e.,

)n+1_Ja [ r _(k.)(IT)] (2.50)
Un )iMtp = Pe [-n \k-- a\n

where the superscript 'n+1' on the left hand side denotes the n + 1 time step while
the right hand side uses the field values at time step n. The subscript 'p' on the left
hand side means this is the interface velocity component due to phase change effect.
So the interface movement in the normal direction due to phase change is

n & n + (U n) At (2.51)
int,p -- tnt ,int~p

After we have this new location of the interface, we can use the process in the
section 2.4.3 to further determine the shape satisfying the momentum balance condi-
tion. The final interface velocity at time step n + 1 obtained from kinematic condition
is the combination of the components owing to phase change and momentum balance.
Hence, the vapor phase velocity at the interface in the normal direction is given
by:

)n1+ nl (u )n+1 (2.52)
(n)v "- Unlint -- Unlint~p

which is the boundary condition for solving the flow field in the vapor bubble.
The liquid phase velocity at the interface in the normal direction, according to
Eq. (1.23), is
), n+1 ) n+1 I) v\. v n+lr co
(Un = (Un]nt *[ (~)] +(p ) +(Un)n (2.53)

which is the boundary condition for liquid phase.

2.4.5 Discontinuity of Material Properties

In the currently available numerical methods reported in the literature for sim-
ulating two-phase flow, most of the single field approaches encounter the numerical








instability caused by large property jumps of the two phases, the maximum density

ratio for phase change problems which is reported so far in the prior publications, to

our best knowledge, is 100. However, for liquid-vapor phase phenomena, the density

ratio often tops 1000. In a typical situation of water at atmospheric pressure, the

density ratio of liquid and vapor is roughly 1600. For some air bubble problems with

no phase change, we saw publications dealing with density ratio as large as 1000 but

with a finite thickness interface treatment to prevent instability.

In the single field approaches, for any cell around the interface, different proper-

ties will appear in the same discretized equation on that cell which cause the stiffness
in that equation so as to deteriorate the numerical stability although easy implemen-

tation of the single field algorithm is quite attractive.

In our approach, because we employ the explicit interface tracking and cut cell

technique around the interface, and our carefully designed computational procedure

solves separately the two fields corresponding to two phases while only on the inter-

face, the two fields communicate with each other to adjust there own flow field to

satisfy the interfacial conditions, the fundamental difference from other approach is

that the discretized equation for any cell including those around the interface will
involve identical material properties corresponding to the phase that particular cell

lies in. Therefore, there is no stiffness in the discretized equation on any cell, the

large density jump will not be expected to exert severe effect on the numerical stabil-

ity in our approach although the interface has zero thickness. This benefit is at the

cost of complex algorithm design and coding involved in the development of present
techniques.

2.5 Discretization on Axisymmetric Geometry

In the simulations, the computational domain is axisymmetric where it is as-

sumed that 0/80 = 0 and the swirl velocity uO is zero.








In a differential formulation, the 2-D dimensional conservation equations, Eqs.
(1.4)-(1.7) written in an axisymmetric cylindrical coordinate system, can be further
expressed as


18 1 a
1- (rur) + 1- (rue)
rOr r9z
9r[_ 19 1a
p)9 +[ (ruru) + (ru-uT)
[Ot r Or rz ZJruU
Qu, 1 0 1 9
1 -9- (ru,.u) + r- (ruzuz)]
T (ruT) + (ruT)
r r 5r-r9zJ


= 0
9P + /I V2u
= -^+t (vl*u-.
Op
= --Ikv2T
- kV2T


(2.54)


or equivalently as


V-u




pcp [ +V.(uT)
Ot I..r


= 0
Op Ur
Or + V, r2
Op
= -z+ V- Vuz
SVzVT
-kV-VT


If the coordinates z and r of the cylindrical coordinate system are replaced by
x and y, the analogy between Eq. (2.54) and the equations in Cartesian coordinates
becomes apparent. If r is also set to 1, Eq. (2.54) becomes identical to those in
Cartesian coordinates except with an additional term -IUr/r2 with u, = u. and
ur. = uy for differential forms.
In terms of the finite volume discretization procedure, the conservation equations
written in integral form can be obtained from Eq. (2.55) by integrating the both sides
over an axisymmetric control volume. The integral equations are the same as those for
planar 2D cases in Cartesian coordinate except for an additional term -kuu,/r2 AV
in the u, momentum equation, where AV is the volume of the control volume. Any


(2.55)









method described in previous chapters are equally applicable for discretizing the 2D

axisymmetric equations in cylindrical coordinates. However, some differences need

to be addressed in discretizing the equations in an axisymmetric control volume.

In addition to the surface integrals of pressure force over cell faces zi = const,

zi-1 = const, rj = const, rj-1 = const, as was the case in planar 2D problems, the

radial component of pressure forces onto the two faces 0 = const has to be considered

also.

The other difference compared to planar 2D problems is the calculation of cell

face areas and cell volume. Since the control volume size in the 0-direction is one

radian, the areas of cell faces are computed in the same way as in the planar 2D case

with the length in 0-direction taken into account, i.e., inclusion of a factor rfc (where

rfc denotes the cell face center). The areas of the front and back faces are calculated

in the same way as the volume in the planar geometry. The volume of axisymmetric

CVs with any number of faces is (see Ferziger & Peric, 1996) :


AV -= 6 (zi-1 zi) (rl1 + r riri-i) (2.56)
i=1

where N, denotes the number of vertices, counted counterclockwise, with i = 0

corresponding to i = Nv.

With the similarity of the governing equations for the planar 2D problems in

Cartesian coordinate and axisymmetric 2D problems in cylindrical coordinate, the

same computer code can be used for both planar and axisymmetric 2D flows; for ax-

isymmetric problems, one sets r = y to account for the length in the third dimension
of the axisymmetric control volume, for planar problems, one sets r = 1 because the

length of the third dimension is unity for planar 2D control volume.




Full Text

PAGE 1

DIRECT UMERICAL SIMULATION OF A TRA SLATI G VAPOR BUBBLE WITH PHASE CHA GE By TAO YE A DISSERTATIO PRESENTED TO THE GRADUATE SCHOOL OF THE U IVERSITY OF FLORIDA IN PARTIAL FULFILLME T OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2001

PAGE 2

To My Family

PAGE 3

ACK OWLEDGEME TS I have been fortunate to have the privilege of working with Professor Jacob Chung and Professor Wei Shyy. Professor Chung, through his constant support and encouragement has made it possible for me to finish this research. Many insightful discussions with Professor Shyy have made my Ph.D. study a very rewarding learning experience. I am also grateful for their tremendous guidance in my personal growth and development. I sincerely thank Professor William Tiederman, Professor James Klausner and Professor Ren-Wei Mei, for serving as members of dissertation committee and gra ciously offering thoughtful advice to help improve this work. I would also like to acknowledge the collaborations with Dr. Udaykumar and Dr. Mittal in conducting the work of Chapter 3 of this dissertation. Most of all I would like to express my deepest appreciation to my family for their patience, love and support throughout all my professional endeavors. The financial support was provided by the Department of Mechanical Engineer ing, Andrew H. Hines, Jr./Florida Progress Endowment Fund and a grant from Eglin Air Force Base. lll

PAGE 4

ACKNOWLEDGMENTS NOMENCLATURE ABSTRACT CHAPTERS 1 INTRODUCTION. TABLE OF CONTENTS lll Vl lX 1 1.1 Problem Statement . . 1 1.1.1 Research Motivation 1 1.1.2 Physics of a Boiling Vapor Bubble in Two-Phase Flows 2 1.1.3 Schematic of Computational Model 3 1.1.4 Assumptions . . . . . . . 4 1.1.5 Governing Equations . . . . 5 1.1.6 Boundary and Interfacial Conditions 6 1.1 7 N ondimensionalization . . . . 9 1.1.8 Dimensionless Parameters . . . 11 1.2 Literature Review . . . . . . . 12 1.2.1 Vapor Bubble Growth or Collapse in a Two-Phase Flow 12 1.2.2 Numerical Methods for Moving Interface Problems 22 1.3 Research Objectives . . . . . . . 28 1.3 1 Research Objectives . . . . . 28 1.3.2 Scales and Dimensionless Parameters 30 2 NUMERICAL METHOD . . . . 32 2.1 Fixed-Grid Sharp Interface Method 32 2.2 Interface Representation Algorithms 33 2.2.1 C 2 Cubic B-spline Algorithm 39 2.2.2 Fairing Algorithm . . . . 44 2.2 3 Validation of Geometric Algorithms 48 2.3 Cartesian Grid Method for Sharp Interfaces . . . . . . 51 2.3.1 Fractional Step Method for Unsteady Transient Solution 59 2 3.2 Discretization on Regular and Cut Cells . . . 63 2.3.3 Imposition of Boundary Conditions on the Interface 7 4 2.3.4 Inversion of Discrete Operators . . . . . 76 2.4 Moving Interface Algorithms . . . . . . . . 80 2.4 1 Change of Phase . . . . . . . . . . 81 2.4 2 Interfacial Conditions . . . . . . . . 82 2.4.3 Iterative Process for Determining the Interface Shape 84 lV

PAGE 5

2.4.4 Determining Interface Shape with Phase Change 86 2.4 5 Discontinuity of Material Properties . 87 2.5 Discretization on Axisymmetric Geometry . 88 3 VALIDATION OF CARTESIAN GRID METHOD 91 3.1 Validation of Global Accuracy of Cut Cell Approach 91 3.2 Fidelit y of The Numerical Methodology 93 3.3 Summary . . . . . . . . 98 4 BUOYANCY-DRIVEN BUBBLE MOTION 100 4.1 Summary of Previous Study . . . . . 101 4.2 Grid Resolution Study and Convergence Criteria 102 4.3 Results For Buoyancy Driven Bubble Motion 103 4.4 Summary . . . . . 115 5 BUBBLE RISE AND GROWTH. . . . . 116 5.1 Stationary Bubble Growth . . . . . 116 5 1.1 Case with Thermal Properties of Water 117 5.1.2 Case with Thermal Properties of Ethanol 119 5 2 Translating Bubble Growth . 119 6 BUBBLE RISE AND COLLAPSE . 132 7 SUMMARY AND FUTURE WORK 138 7.1 Summary of Present Work . 138 7 2 Contributions of Present Work 139 7 3 Future Work . . . . . 140 7.3.1 Applications . . . 140 7 3 2 Numerical Enhancements 141 APPENDIX A HEAT TRANSFER-CONTROLLED LIMITS OF SPHERICAL B U BBLE GROWTH . . . . . . . . 143 A 1 Conduction Dominated Growth. A.2 Convection Dominated Growth REFERENCES . . . . BIOGRAPHICAL SKETCH. V 143 144 145 151

PAGE 6

NOMENCLATURE Cpz heat capacity of the liquid phase (J/kg K) cP v heat capacity of the vapor phase (J/kg K) CD drag coefficients Fr Froude number(= u~/gL) Fo Fourier number ( = azt/(2R 0 )2) g normal gravity h heat transfer coefficient ( = q" / b.T) Ja Jakob number(= pzcp 1 b.T/pv>..) k 1 thermal conductivity of the liquid phase (W/m K) k v thermal conductivity of the vapor phase (W/m K) L length scale ( m) n normal direction Nu Nusselt number(= hL/k) Pd dynamic pressure (kg/m s 2 ) p 1 pressure in the liquid phase (kg/m s 2 ) P v pressure in the vapor phase (kg/m s 2 ) Pe Peclet number (= pzcp 1 urL/k1) Vl

PAGE 7

Pr II q r Re t T u z Prandtl number ( = zl a1) heat flux (W/m 2 ) radial coordinate initial bubble radius ( m) Reynolds number (= p1urL/ ) time (s) temperature (K) temperature at the bubble interface (K) temperature in the liquid phase (K) temperature in the vapor phase (K) velocity component in the r direction or reference velocity scale ( m/ s) velocity component in the z direction ( m/ s) terminal velocity for the rising bubble (m/s) interface velocity in the normal direction ( m/ s) Weber number (= p1u~L/a) axisymmetric coordinate Greek Symbols a 1 thermal diffusivity of liquid phase (m 2 /s) /3 under-relaxation factor for interface updating procedure 'Y dimensionless radius for collapsing bubble ( = R( t) / Ro) K curvature of the interfacial curve vu

PAGE 8

D.T degree of superheat or subcooling (K) A latent heat of phase change (J/kg) 1 dynamic viscosity of the liquid phase (kg/ m s) v dynamic viscosity of the vapor phase (kg/ms) Pt density of the liquid phase (kg/m 3 ) Pv density of the vapor phase (kg/m 3 ) CY surface tension coefficient ( N / m) r stress tensor 'TH dimensionless time ( = 4J a 2 a 1 t/1r R,e/) 111 kinematic viscosity of the liquid phase (m 2 / s) llv kinematic viscosity of the vapor phase ( m 2 / s) generic variable for illustration of finite volume integration Vlll

PAGE 9

Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy DIRECT NUMERICAL SIMULATION OF A TRANSLATING VAPOR BUBBLE WITH PHASE CHANGE By Tao Ye May 2001 Chairman: Dr. Jacob N. Chung Major Department: Mechanical Engineering A fixed-grid direct numerical simulation method has been developed and verified for single bubble dynamics, heat transfer and phase change phenomena. The original contribution of this research is the first introduction of a sharp interface between the vapor and liquid phases, which allows the simulation of density ratios up to the order of 1000. The interface location is obtained as part of the solution which accounts for all coupled dynamic and thermal effects. The current method is believed to be the only method available at present which resolves the curvature calculation issue, treats the interface as a sharp discontinuity honors the mass conservation principle, is capable of handling the large property ratios and phase change For the isothermal cases, the predictions of drag coefficients and the shape defor mation of bubbles for Reynolds numbers ranging from 2 to 100 and Weber numbers from 2 to 20 were shown independently to agree with previous published results. IX

PAGE 10

It was also demonstrated that the current method is capable of predicting accu rate results for a wide range of Re, We, Ja, Pe numbers and property jumps at the interface For a stationary and growing bubble, the predicted growth rate was found to agree with the theoretical limit of ex t 1 1 2 However, for a rising and growing bub ble the predicted growth rate was approaching the theoretical limit of ex t 2 3 for a spherical bubble but did not reach it exactly owing to bubble deformation. The empirical correlation proposed for bubble co ll apse is assessed by the simu lations in the present work. The present numerical model predicts that the dynamics of a collapsing bubble is qualitatively, but not quantitatively similar to the empirical correlation proposed by other investigators. X

PAGE 11

CHAPTER 1 INTRODUCTION 1.1 Problem Statement 1.1.1 Research Motivation Liquid-vapor phase change is a phenomenon that is frequently encountered in daily life and that plays an important role in the power, chemical, petroleum and electronics industries among others. For example, the power generation industry takes advantage of the high heat transfer rates associated with the liquid-vapor phase change in the heterogeneous nucleate boiling process to efficiently extract energy stored in fossil, nuclear or solar fuels. The central mechanism of heat transfer during nucleate boiling is the so-called ebullition cycle: a complete process of liquid heating, nucleation, bubble growth and departure (see Carey, 1992). Hence understanding and prediction of the microscopic vapor bubble behavior have direct importance for better applications of nucleate boiling heat transfer. Despite the wide applications of liquid-vapor phase change, a fundamental un derstanding of such phase change process is far from being complete. Analytical analysis of the liquid-vapor phase change has faced mathematical difficulties because of the nature of nonlinear, partial differential equations which govern the process. In scientific and engineering research, experimentation and numerical simulation have been mainly employed as basic methods of approach. Both means have respec tive merits and they need to complement each other. Meaningful experiments are of singular significance for phenomena that are beyond the capability of current compu tational technology, and they are worthwhile in providing verification for numerical 1

PAGE 12

2 computation and guidance for mathematical modeling. On the other hand, numeri cal computations can be used to predict microscopic transport mechanisms that are difficult to measure and visualize in an experiment. Because of limited resources, most of the experimental projects may be restricted to only one working fluid, simple geometry, and limited thermodynamic conditions. Numerical computations might therefore be useful for expanding the data base to a wider spectrum and providing the information that is otherwise hard to obtain. Although direct simulations of industrial scale phase change processes are cur rently out of reach, a fundamental understanding of the complex physics at small scale can provide much insight to quantitative predictions of large scale aspects of processes. Despite the fact that computational fluid dynamics has advanced to a point where solutions of relatively complex flows without phase change can be routinely obtained, numerical techniques for the combination of fluid flow, moving interface and phase change remain as one of the unresolved challenges since it is necessary to model the interaction of fluid flow with heat and mass transfer across moving phase interfaces In the present research, our major goal is thus to develop a numerical capability for simulating liquid-vapor phase change at the small scale, namely, the vapor bubble behavior in a pool of liquid. 1.1.2 Physics of a Boiling Vapor Bubble in Two-Phase Flows Nucleate boiling is one of the most efficient heat transfer modes. This mode of thermal energy transport is accomplished by two integral stages, which transport the thermal energy from a heater surface to the bulk of a working fluid. The first stage is the nucleation and growth of bubble embryos on the heater surface. The second stage is the migration of the bubbles into the bulk of the fluid The present research is devoted to the fundamental understanding of the second stage.

PAGE 13

3 In the second stage, if the bulk of the liquid is superheated the bubble will continue to grow while it is translating. The growth rate of the bubble in this stage is primarily limited by the rate of arrival of the heat necessary for vaporization at the bubble surface. This stage of growth is called heat transfer-controlled growth. The pressure at the bubble surface inside the bubble approaches the liquid pressure adjacent to the bubble surface; hence the temperature at the bubble surface is ap proximately equal to the equilibrium saturation temperature of the system pressure. Translating vapor bubbles in a superheated or subcooled liquid environment can be found in both pool boiling and convective boiling. In pool boiling, bubbles rise into the bulk of the pool due to buoyancy. In convective boiling bubbles are carried into the bulk by flow currents. The interaction between the vapor bubble and its liquid surroundings takes place after the departure of this bubble from the heater surface. The interaction between the dispersed phase (bubbles) and the continuous phase (liquid) involves exchange of momentum, thermal energy and mass. The in teractions are two-way coupled which means that the bubble behavior is affected by the liquid flow and the liquid flow is also modulated by the presence of the bubbles. These two-way coupled interactions include momentum exchange through the drag and vaporization/ condensation, energy exchange due to convection and latent heat transport of phase change and mass exchange during phase change The present research will focus on numerically simulating the single vapor bubble growing or collapsing due to phase change while translating in a liquid. 1.1.3 Schematic of Computational Model The schematic of computational model for a growing or collapsing vapor bubble rising in a quiescent liquid is illustrated in Figure 1.1. The motion of the vapor bubble is simulated in a cylindrical domain filled with liquid of the same substance. The liquid is superheated in one case or subcooled in the other case. The major simplification is that the system is assumed axisymmetric.

PAGE 14

Liquid z ' __ ..,. r Figure 1.1: Schematic of computationa l model. 1.1.4 Assumptions 4 The two-fluids model is adopted as the mathematical basis which governs the liquid vapor two-phase flow of a sing l e trans l ating vapor bubble; i e. a separate set of mass, momentum and energy conservation equations are solved for each phase

PAGE 15

5 (dispersed vapor phase and continuous liquid phase surrounding the bubble ) while sharp discontinuities of material properties are maintained across the interface of zero thickness. The assumptions introduced throughout this study are summarized as follows: Axisymmetric cylindrical coordinate system Single component system Newtonian constant property fluids in each phase Thermal equilibrium at the phase interface which means temperatures at the interface in the liquid and vapor phases are equal 1.1.5 Governing Equations With the above assumptions, the mass, momentum and energy conservation equations for both phases can be written as V-u = 0 (1.1) (1.2) pep [8:: + V (uT)l = V (kVT) (1.3) where for numerical efficiency pg = -i1(pgh) has been combined with p in place of the pressure. If the actual pressure is needed hydrostatic pressure pgh has to be subtracted from the pressure computed in the equations. Thereby the hydrostatic pressure (body force term) contribution will only appear in the normal stress balance at the bubble interface (see Eq. (1.12)). For this research the axisymmetric cylindrical coordinates are employed. Equa tions (1.1)-(1.3) thus become the following with 8/80 = 0 and the swirl v e locity

PAGE 16

6 Uo = 0: 1 o o -l (rur) + l ( Uz) = 0 r ur u z (1.4) [ OUr l O O l op l O 0 p fJt +-;. Or (rurUr) + O Z (UzUr) = or +-;. Or (rTrr) + OZ (Trz) Too r (1.5) op 1 o o -l + -l (rTzr) + l (Tzz) u z r ur u z (1.6) !~ (rka.I') + (ka.I') r or or oz oz (1. 7) where the stress tensor components are Tzz (1.8) 1.1.6 Boundary and lnterfacial Conditions Apart from the governing equations, boundary conditions are needed to make the problem well posed. The interfacial conditions, which are simply the statement of those conservation principles at the interface, must also be satisfied (see, e.g., Delhaye 1974 ; Ishii, 1975; Kataoka, 1986; Sadhal, Ayyaswamy & Chung, 1996). Boundary conditions: The outflow (zero gradient of the velocity) conditions are specified on the top, left and right sides of the domain. At the bottom boundary, the symmetric condition for all variables is used.

PAGE 17

7 Interfacial mass conservation conditions: Due to phase change, for an interface moving at a velocity Uint, the mass con servation conditions on the interface stipulates that rh = Pl (uint u1) n Pv (uint U v) n (1.9) where mis the interfacial mass flux (positive for evaporation), P1 and Pv are the liquid and vapor densities respectively, while the tangential velocity is always continuous Interfacial momentum balance conditions: The forces balance across the interface is expressed as The normal component of the equation (1.11) is Pl Pv + aK, = (rrn)n (rvn)n {Pi [(un)1 (un)int] (un)1 Pv [ ( Un)v ( Un)int] ( Un)v} (1.10) (1.11) (1.12) where T represents the viscous stress tensor, a is the surface tension coefficient which is assumed to be constant, "' = ( J 1 + Ji 2 ) with R 1 and R 2 the principal radii of interface curvature. The pressure in Eq. (1.12) is the total pressure including the hydrostatic contribution; i.e., the actual pressure is obtained by subtracting pgh ( or J.rh in dimensionless form) from p solved in Eq. (1.1), (1.2). The normal component of stress balance is reduced to Young-Laplace equation when there is no motion in both phases.

PAGE 18

8 Interfacial energy conservation conditions: The energy conservation relation across the interface involves both the sensible and latent heat portions (1.13) where A is the latent heat of evaporation. Equation (1.13) is used to estimate the normal velocity component of the inter face driven by the phase change. Interface temperature conditions: By thermal equilibrium, the temperature continuity condition reads as Tz = Tv = ~nt (1.14) where Tint is the temperature at the phase change interface. A condition on the interface temperature must be specified to complete the formulation. In the case where liquid boiling occurs on a horizontal flat surface, i.e., curvature K, = 0, the interface temperature Tint is equal to saturation temperature Tsat of the two phase mixture at the corresponding ambient pressure p 00 For a curved interface between the liquid and vapor bubbles, when surface tension effect is taken into account, the thermal and mechanical equ ilibrium at the interface reveals that the local saturation temperature Tint along the interface will deviate slightly from that constant saturation temperature Tsat for the flat interface (see Alexiades & Solomon, 1993), which is the so-called Gibbs-Thomson effect. For the heat transfer contro ll ed bubble growth, vapor pressure Pv at the bubble interface is equal to the system pressure p 00 which is the equilibrium pressure of Tsat, i .e., Pv = Poo (1.15)

PAGE 19

9 In this case, a sufficiently accurate formulation for the interface phase change tem perature reads {1.16) 1.1. 7 N ondimensionalization To nondimensionalize the governing equations in the liquid and vapor phases as well as boundary and interfacial conditions, the reference scales are length L, velocity Ur, time tr = L/ur. The characteristic temperature scale is 11T =Th-Tc= { T 00 -Tsat growth T. d T 1 b 00 an sat are, respective y, am ient temperature and Tsat Tc:,o collapse saturation temperature. With the basic scales chosen as above, the dimensionless variables are defined as x = x/ L u = u/u r, t* = t/tr, p* = p/(piu/), T* = (T Tc)/ 11T g = g/g, p* = Pf Pref, *=/ref, k* = k/kref, c;,* = c;,/c;,ref' where the liquid properties are taken to be the reference properties, i.e., Pref = Pl, ref = ,, kref = kt, Cvref = Cvt Finally, the dimensionless equations in the liquid and vapor phases, omitting the superscript *, can be written as Liquid Phase: Vapor Phase: au -+V(uu) at 8T -+ V-(uT) at V-u 0 1 2 -Vp+-V u Re 0 (1.17) (1.18) (1.19) {1.20)

PAGE 20

where Reynolds number Peclet number lnterfacial Condition: Mass continuity condition ( v ) 1 2 -v1p+ -v7 u Re Re Pe PtUrL PtCp,urL k1 (un)1 (Ti) (un\ 1-(Ti) Momentum balance condition 10 (1.21) (1.22) (1.23) where Pl = Pd + l/ Fr g(zout z), Fr = u~/ gL is the total pressure including hydrostatic and dynamics pressure. Zout is the level of the liquid pool. Energy conservation condition (1.25) Interface temperature condition 1int = I'K, (1.26)

PAGE 21

where Weber number Froude number Jakob number Parameter 1. 1. 8 Dimen s ionles s Parame t ers W e Fr Ja I' gL Pl ep,t1T P v >. Tsat <7 i1T p1>.L 11 After examining the dimensionless governing equations and the associated bound ary and interfacial conditions the major dimensionless parameters namely R e, Fr, Pe, Pr We and Ja and their respective physical meanings are listed as follows: R e = Fr Pe Pr W e PtUrL Reynolds number flow inertia force viscous force u2 r gL Fro u de number flow i nertia force grav i ty force P1Cp 1 UrL k1 convection conduction PtUrL Cp 1 k1 RePr Peclet number k1 lit = a1 = -, a1 Pi P1Cp, Prandtl number momentum diffusion capacity thermal diffusion capacity P1U~L Weber number (j flow hydrodynamic pressure force surface tension force

PAGE 22

12 Ja Pt CpliJ.T Pv. >. Jakob number sensible heat transfer latent heat transfer 1.2 Literature Review 1.2.1 Vapor Bubble Growth or Collapse in a Two-Phase Flow Since we are concerned with the heat transfer-controlled growth or collapse of a single vapor bubble translating in a liquid pool, the review of previous works will be focused only on those directly relevant to thermally controlled bubble growth or collapse after the bubble has departed from the heater surface. A number of exce llent publications in the literature are dedicated entirely to the fundamental understand ing of the bubble growth dynamics on the heater surface {first stage of boiling), for example, Klausner, Mei & Bernhard (1993); Zeng, Klausner & Mei (1993b ); Zeng, Klausner, Bernhard & Mei (1993a); Chen, Klausner & Mei (1995) and Mei, Chen & Klausner (1995). These works are not directly relevant to the scope of the proposed research, and will not be discussed in detail further here. Vapor bubble growth or collapse has drawn extensive study over many years because of its theoretical implication to the nucleate boiling and other two phase physical phenomena. Due to its intrinsic difficulties, mainly the nonlinearit y of con vection and the coup ling of the bubble evolution with the surrounding liquid motion as well as the large number of variables involved it would seem that this problem has been resolved to a tenable degree of understanding. The heat transfer controlled bubble growth or collapse can be generally divided into four basic situations: stationary bubble growth, translating bubble growth, sta tionary bubble collapse and translating bubble collapse. The researchers in this field seem to follow this line to report results for one or more cases in any particular paper. In our discussion, we will also follow this convention to look at the prior contributions to all these cases.

PAGE 23

13 Stationary bubble growth The majority of the research on the heat transfer controlled growth of station ary vapor bubbles is based on the simplified model that a spherical bubble in an unbounded fluid undergoes spherical growth. The simplified model does not solve for the hydrodynamics of the fluid flow, nor does it locate the moving phase interface. With the assumption of spherical symmetry, the continuity and momentum equa tions for an incompressible fluid lead to the famous Rayleigh equation (see Plesset & Prosperetti, 1977) d 2 R 3 (dR) 2 1 { 2a 1 (dR)} R dt2 + 2 dt = Pl Pv (Ts) Poo R 4 R dt (1.27) in which T 8 (t) is the liquid temperature at the bubble surface and it is assumed that the vapor pressure is the thermodynamic equilibrium pressure at the temperature T 8 To analyze the heat transfer controlled growth, the liquid temperature Ts at the bubble surface is obtained by solving the energy equation 8T R 2 dR 8T a1 8 ( 2 8T) at + i2 dt 8r = r 2 8r r 8r in which a 1 = kz/ (pzcp 1 ) is the liquid thermal diffusivity. (1.28) Plesset & Zwick (1954), Forster & Zuber (1954), Birkhoff, Margulies & Horning (1958) and Scriven (1959) considered the heat transferred controlled spherical vapor bubble growth without translation by solving the energy equation (1.28) subject to certain approximations in conjunction with the dynamical equation (1.27) and thermodynamic equilibrium relation Pv = Pv(Ts)Their results of bubble radius varying with time are in good agreement with the experimental data of Dergarabedian (1953) and, later, Florschuetz, Henry & Khan (1969); basically, what they concluded is that the bubble radius of the heat transfer controlled growth is proportional to (time) 1 1 2 Prosperetti & Plesset (1978) re-examined two assumptions introduced in those earlier analyses of bubble growth: the validity of a thin thermal boundary layer near

PAGE 24

14 8.0E-03 7.0E-03 6.0E-03 5.0E 03 E u, 4.0E 03 .2 "C m a: 3 0E-03 2.0E -0 3 1.0E 03 0 0E + 00 0 0 1 0.2 0.3 0 .4 Time (sec) Figure 1.2: Bubble radius vs. time for a vapor bubble growing in superheated water at atmospheric pressure with 5C superheat [Prosperetti & Plesset (1978)]. the bubble wall and the empirical linear relation for the vapor pressure with vapor temperature. It was showed that, upon comparison with numerical results of Donne & Ferranti (1975), both assumptions do not lead to serious errors in the radius-time behavior for superheat ranges of practical interest. The relation of the bubble radius with time for heat transfer controlled growth is given as (1.29) in which k1 is the liquid thermal conductivity, a 1 = kt/(p 1 epi) is the liquid thermal diffusivity, (T 00 Tsat) is the liquid superheat, .X is the latent heat, Pv is the vapor density. The radius increasing with time based on this formulation is plotted in Figure 1.2.

PAGE 25

15 Translating bubble growth The stationary vapor bubble is an ideal exception. In practice, the vapor bubble usually exhibits more or less translatory motion. Although the translatory motion is not expected to play a significant role if liquid inertia is the dominant mechanism, as stated by several investigators it has a large influence on the growth or collapse if heat transfer is controlling the phase change Under such circumstances, the heat transfer would be enhanced by translation, thus resulting in higher growth or collapse rates. The experimental investigation of Darby (1964) studied simultaneous growth and translation characteristics of vapor bubbles issuing continuously from a nucleation site. The bubble radius was found to vary as (time) 2 1 3 during the entire recorded history. It indicates that, as expected, the translational motion increases the growth rate over that for stationary bubbles. Ruckenstein & Davis (1971) assumed a potential flow in the region surrounding the bubble to analyze the effect of bubble translation on vapor bubble growth in a superheated liquid. He found that for a moderate Jakob number, the translation can substantially increase bubble growth rate over the rates predicted for a stationary growing bubble. For sufficiently large Jakob numbers (Ja > 50), bubble growth is not greatly affected by the relative motion between the bubble center and the surrounding liquid. The results predicted by Ruckenstein & Davis (1971) are in good agreement with the experimental data of Florschuetz, Henry & Khan (1969). In Table 1.1 we summarize the above mentioned important research publications for heat transfer controlled growth of a vapor bubble. Stationary bubble collapse In bubble collapse, the radius variation produced by a given interfacial heat flux increases in time in contrast to the decrease for the bubble growth case. This interface

PAGE 26

16 Table 1.1: A list of representative investigations on stationary and translating vapor bubble growth Authors Plesset & Zwick (1954) Forster & Zuber (1954) Birkhoff, Mar gulies & Horning (1958) Scriven (1959) Prosperetti & Plesset (1978) Darby (1964) I Approach I Assumptions & Conditions Analytical Spherical symmetry Stationary growing vapor bubble Heat transfer controlled vapor bubble growth The inertia and viscosity of the liq uid are neglected Surface tension is neglected Heat transfer by conduction only Constant thermal and physical properties Experimental Translating vapor bubble growth Ruckenstein & Analytical Davis (1971) Spherical symmetry Translating vapor bubble growth Legendre, Boree & Magnaudet (1998) Numerical Potential flow in the liquid Heat transfer controlled growth Spherical symmetry Full Navier-Stokes equations are solved Heat transfer controlled growth Stationary and translating bubble Decoupled bubble evolution and liquid motion I Major Results The bubble radius is pro portional to the square root of time, i.e. (time) 1 1 2 Bubble ra dius varies as (time) 2 1 3 The growth rate is substan tially increased at moderate Jakob number (Ja < 50) Confirmed all above prior an alytical results and studied the hydrodynamic forces behavior makes condensation basically more difficult to describe than the evapora tion. For the stationary bubble collapse, by virtue of spherical symmetry assumption, Florschuetz & Chao (1965) examined the relative importance of the liquid inertia and heat transfer effects on the vapor bubble collapse. A dimensionless parameter is iden tified to characterize the collapse mode: either liquid inertia-controlled collapse heat transfer dominated collapse or an intermediate case where both liquid inertia and heat transfer are of comparable importance For heat transfer-controlled collapse

PAGE 27

17 they found that the collapse rate is relatively slow and decreases as the collapse pro ceeds, the radius-time relation is bounded by two limiting heat transfer-controlled collapse curves that resulted from two approximations for temperature solutions: (1.30) and 1 (2 2 ) TH= 3 7 +1 -3 (1.31) in which dimensionless bubble radius 1 = R/ Ilo, dimensionless time (1.32) with Ja being the Jakob number, a 1 the liquid thermal diffusivity, Ilo the initial radius of the bubble. The results are visualized in Figure 1.3 which is from the paper of Florschuetz & Chao (1965) The analytical results Eq. (1.31) of Florschuetz & Chao (1965) are valid at high Jakob numbers. In order to obtain results valid at finite Jakob numbers, Okhotsimskii (1988) solved numerically the heat transfer controlled collapse based on the same theoretical model of Plesset & Zwick (1954) and tabulated the time evolution of the bubble radius over a wide range of Jakob numbers ( 102 < Ja < 10 3 ). In his analysis, he isolated the Jakob number to be the only dominant factor of the collapse process. Translating bubble collapse For more practical cases of translating vapor bubble collapse, Wittke & Chao (1967) obtained numerical and experimental information of the influence of transla tory velocity on the collapse of a vapor bubble. oting that the translatory motion is not expected to be a significant factor if liquid inertia is the dominant mechanism for collapse, they examined the bubble traversing at a constant velocity collapsing under

PAGE 28

1.0 -----r---~---,---,---r----r--~---, ~o.P >--091-\.,---l--~~--1f----f--~---t--,---., ; a: O.SJ---:--.-+-~--+---TBASED CN PLESSET-ZW10< w 'ZERO OR:>ER' TaHRATI..RE SCl..l.JTCN t!IJ-.--1f--___;~---+---'~..,t:: -~=---~ I 0 0 2 04 0.6 OB 1.0 DMENSICNLESS TIME, Tt. 11 : l 0 12 L4 1.6 18 Figure 1.3: Bubble radius vs. time for heat transfer controlled vapor bubble collapse [Florschuetz & Chao (1965)]. heat transfer controlled conditions. They found that, again as expected, because the relative contribution to heat transfer at the bubble surface due to translatory mo tion is greater when the Jakob number is smaller, the translatory motion produces a greater effect for the smaller Jakob number. Overall, the faster the bubble translates, the more rapidly it collapses. Tokada, Yang & Clark (1968) and Dimic (1977) also studied theoretically the ef fect of translational velocity of a bubble on the collapse rate. They basically obtained similar results on the collapse rates. In these studies the bubble is relatively large (2mm < Flo < 4mm) and because of reasons which are not completely clear, it was observed in some experiments on bubble collapse that bubbles exhibit almost a constant rising velocity (see Wittke & Chao, 1967; Chen & Mayinger, 1992); hence a constant translation velocity is assumed in all these studies.

PAGE 29

19 Moalem & Sideman (1973) analyzed the effect of variable rising velocities as sociated with relatively small bubbles (~ < 1 mm) as well as the effect of a cross flow; their results were in good agreement with the experimental data of Abdelmes sih Hooper & Nanngia (1972) What they concluded is that the effect of variable velocity as compared to a constant velocity motion is relatively small. All of the above mentioned results assumed the potential flow for outer liquid surrounding the bubble. Chen & Mayinger (1992) performed very sophisticated measurements to deter mine the heat transfer at the phase interface of a collapsing vapor bubble rising in a subcooled liquid The temporal decrease in the bubble diameter is given in the correlation using experimental data as ( 1.33) in which 'Y = R/ ~ Reynolds number is based on initial bubble diameter. The Fourier number is 0:1t Fo = (2~)2 (1.34) Recently, Gumerov (1996) numerically obtained the solutions of the same theoretical model as that of Wittke & Chao (1967) for the collapse of a vapor bubble with constant translatory velocity. His results are shown to agree substantially better with the experiments in a number of cases. In Table 1.2 some of these studies for stationary and translating vapor bubble collapse are summarized. Most recently Legendre, Boree & Magnaudet (1998) made a comprehensive study of heat transfer-controlled spherical bubble growth and collapse without and with constant translation velocity Basically, they reexamined all four cases, and by comparing with those representative prior results, confirmed previous analytical anal ysis. In their study a numerical approach is employed to solve the full Navier-Stokes

PAGE 30

20 Table 1.2: A list of representative investigations on stationary and translating vapor bubble collapse Authors Florschuetz & Chao (1965) Okhotsimskii (1988) Wittke & Chao (1967) Chen & Mayinger (1992) Legendre, Boree & Magnaudet (1998) I Approach I Assumptions & Conditions Analytical and Spherical symmetry Stationary collapsing bubble Experimental Surface tension and viscosity are ig nored Numerical Numerical and Experimental No relative velocity between liquid and bubble wall Heat conduction in vapor is ignored Pressure in the vapor phase is uniform Perfect gas for vapor Spherical symmetry Stationary collapsing vapor bubble Heat transfer controlled vapor bubble collapse The inertia and viscosity of the liq uid are neglected Surface tension is neglected Constant thermal and physical properties Spherical symmetry Translating vapor bubble Heat transfer controlled collapse Constant translational velocity Potential outer flow Incompressible for Liquid Experimental Translating vapor bubble Heat transfer controlled collapse Numerical Spherical symmetry Full Navier-Stokes equations are solved Heat transfer controlled collapse Stationary and translating bubble Decoupled bubble evolution and liquid motion I Major Results Three collapse modes are classified and a dimensionless parameter is identified to characterize the collapse mode Time evolu tion of bubble radius for a wide range of Jakob numbers tabulated. is The translation motion generally increases the collapse rates. Correlations of Nusselt num ber and time evolution of bubble radius are given Confirmed all above prior analytical results and studied the hydrodynamic forces

PAGE 31

. ~ 0 C .5 free-stream condition 160R 21 (a) (b) 80R Figure 1.4: Computational configuration of investigation of Legendre, Boree & Mag naudet (1998). equations and temperature equation in the liquid to determine the temperature field around the bubble surface so as to obtain the growth or collapse rates. Despite the spherical symmetry assumption and disregarding of the coupling between the bubble evolution and the liquid motion, their solutions to our knowledge, are the most ad vanced ones for the spherical bubble growth and collapse The second contribution of their research is regarding the hydrodynamic force experienced by the bubble with a time-dependent radius. The computational configuration of Legendre, Boree & Magnaudet (1998) is shown in Figure 1.4.

PAGE 32

22 1.2.2 Numerical Methods for Moving Interface Problems The liquid-vapor two-phase problem we outlined in chapter 1 is a special case of more general situations involving multi-fluids or multi-phases in a computational domain as shown in Figure 1.5, where multiple immersed interfaces are embedded in the flow field The regions enclosed by the interfaces designate certain phases which could be solid, liquid or gas phases and might be different from that of the surrounding fluid. The interfaces can be stationary or moving, deforming through interaction with the flow around it. Such a scenario may be encountered in several types of flow problems. For instance, the interface may represent a solid geometry over which fluid flows. This solid boundary may be rigid as in pure external flow problems or be flexible exhibiting motions under the influence of the external flow as in many fluid-structure interaction problems In other situations, the interface may represent a phase boundary at which a phase transition occurs. One of the examples is the liquid-vapor phase change problem involving vapor bubbles as in the present research the enclosed phase is the dispersed vapor phase and the surrounding one is the continuous liquid phase. In an inverse situation involving liquid drops, the enclosed phase is the liquid phase, and continuous phase is the vapor phase. In these cases, sharp discontinuities of the fluid properties such as density, viscosity and thermal conductivity arise across the interface. Other examples of phase boundary include solid-liquid phase of the solidification in material processing problems The interface may also simply represent a boundary separating two immiscible fluids as in the case of bubbles or drops immersed in an ambient fluid. Whether the immersed interface represents a phase interface or a boundary be tween two fluids with different transport properties the numerical simulation of these processes has to deal with moving interfaces. The numerical investigation of moving interface problems has posed a difficult challenge for a long time especially when a moving phase change interface and extremely large ratios of fluid properties across the

PAGE 33

{} Phase 0 an 23 Figure 1.5: Illustration of multiple immersed interfaces in a computational domain interface exist simultaneously, as in the present study of the liquid-vapor phase change problem. Despite its difficulties in the past several decades significant advances have been made in numerical techniques for obtaining the solutions of immersed moving interface problems In this section we will briefly examine the available numerical approaches for immersed, moving interface problems with fluid flow. Numerical techniques for immersed, moving interface problems could be classi fied into three broad categories (see Shyy, Udaykumar Rao & Smith, 1996) : body fitted grid approach fixed grid approach and mixed Eulerian-Lagrangian approach. The natural approach is to use body-fitted type mesh which could be either struc tured or unstructured (e.g. finite element mesh). Relatively complex two-dimensional free boundary problems were simulated using this approach (see Ryskin & Leal 1984;

PAGE 34

24 Kang & Leal, 1987; Dandy & Leal, 1989). A body-conforming mesh makes the bound ary conditions specification on the interface relatively easy and accurate because the interfacial conditions are applied at the exact location of the interface without any smearing operation The sharp discontinuity of material properties across the inter face is always maintained. However in the body-fitted approach the computational expense is an issue because the grid has to be modified or regenerated continuously to account for the movement of the interface. The body-fitted approaches encounter major difficulty in dealing with the topological complexity resulted from arbitrary interface movement with desired mesh quality. Besides, the governing equations are not in a simple form for body-fitted mesh which causes extra difficulty for the over all solution procedure. Multiple connectivity makes designing an efficient algebraic solver a difficult task. Another completely different approach is to use the fixed Cartesian grid. In this purely Eulerian approach, the exact location of the interface is not explicitly known, in contrast the interface details are determined from solving a field variable. Within this approach, there are a few variants, for example, the Volume of Fluid (VOF) method (see Hirt & ichols 1981; Brackbill, Kothe & Zemach 1992; Kothe & Mjolsness, 1992; Rider & Kothe, 1998; Scardovelli & Zaleski, 1999) the level-set method (see Osher & Sethian 1988; Sethian, 1990; Sussman, Smereka & Osher 1994; Sussman, Fatemi, Smereka & Osher, 1998) and phase-field method (see Caginalp, 1984; Kobayashi 1993; Antanovskii, 1995; Anderson & McFadden 1997; Jacqmin, 1999) amongst others. When surface tension on the interface can not be ignored, the accurate estimation of the interface curvature is crucial to successfully predicting the interface evolution. Under these circumstances where explicit interface tracking is needed the mixed Eulerian-Lagrangian approach seems attractive because in this approach, apart from the fixed structured grids (Eulerian part) for the flow field, interfaces are explicitly

PAGE 35

25 described by curves or surfaces using computational geometry tools (Lagrangian part) and the movement of those curves and surfaces are determined by the computational process. In the fixed Cartesian based approach, the grid generation is greatly simplified and the governing equations retain the relatively simple form in Cartesian coordi nates. Because the grid is Cartesian, many efficient solvers are easily obtained to reduce the computational cost. In the mixed approach, a major issue is how to handle the interaction between a moving interface and a flow field. The interface would cut through the underlying Cartesian grid in an arbitrary manner which resulted in fragmented control volumes ( cut cells) around the interface. Fragmented, small cut cells would cause numerical instability. To circumvent this difficulty, the so called immersed boundary technique has been widely used to simulate multi-fluid flows (see Peskin, 1977; Unverdi & Tryggvason, 1992; Goldstein, Handler & Sirovich, 1995; Juric & Tryggvason, 1996, 1998; Kan, Udaykumar, Shyy & Tran-Son-Tay, 1998). In this technique, the interface is not kept sharp but is represented on the fixed underlying mesh by a thin layer of finite thickness around the interface. A scalar indicator function is then constructed from the known position of the interface to smooth out the jumps in the material properties across this thin layer. The information transfer between the interface and a flow field is done by introducing additional source terms in the governing equations. This technique is quite robust, however its fuzzy interface treatment does not reflect closely the physical reality of the zero thickness interface, and thus it is not suitable for accurately resolving any boundary layer type of behavior in the thermal and flow fields when convection is important. A strictly conservative treatment of small cut cells is the so-called cut-cell tech nique. In this treatment, typically, small cut cells are absorbed into neighboring cells

PAGE 36

26 to form mostly quadrilateral cells in the vicinity of the interface. Thereby, the inter face is treated as a sharp discontinuity and the conservation property of the solution procedure is always retained. The interface effect is taken into account by the flow solver through the reshaped control volumes near the interface. The finite volume cut cell techniques have been used in dealing with the complex solid boundaries in the flow field in conjunction with Cartesian mesh (see Quirk 1994; Udaykumar & Shyy, 1995; Pember, Bell Colella, Crutchfield & Welcome 1995 ; Udaykumar, Kan Shyy & 'Iran-Son-Tay, 1997; Ye, Mittal, Udaykumar & Shyy 1999). In those applications, the relative locations between solid boundaries and grid cells are fixed throughout the solution process. There is no need to determine the boundary shape in the solution process. There is no interaction between rigid boundaries and flow field which makes it much easier to implement the cut cell technique. In the applications where free, moving interfaces with phase change e.g. in terfaces separating bubbles or drops from the surrounding fluids are involved. It is much more challenging to implement the cut cell technique because the interface shape is not fixed any more, rather depending on the dynamic conditions on the interface. In other words, the instantaneous interface shape must be determined as part of overall solution process based on the interaction of the interface with the flow fields in both sides of the interface. Owing to this requirement no reports has been found in the literature on the application of finite volume cut cell techniques t o free interfaces like those of bubbles or drops. In the present research, a fixed grid sharp interface moving boundary method is developed and applied to simulate the moving bubble interfaces To summarize, pure body-fitted grid approaches are seldom used nowadays they are more often combined with fixed grid methods to improve the efficiency In such configurations most of the domain is covered by the fixed grid body-fitted grid is

PAGE 37

27 only used around the curved interface to capture the movement of the interface. A communication mechanism has to be designed between global fixed grid and local body-fitted grid. Fixed grid methods are mainly aimed at industrial simulations. In these meth ods, normally a homogeneous system is applied to the entire domain. As a result a single set of governing equations are solved in a domain including multiphase with different properties. The implementation of these methods is relatively easy. Despite its simplicity, a sharp discontinuity for the transport properties at the phase front should be preserved in conformity with the physical reality. In fixed grid methods, the interface is not explicitly known and has a finite thickness, or sometimes, al though the interface itself can be implicitly captured by some field function to be sharp, large property jumps across the two phases still have to be smoothed over a few grid points around the interface to overcome the numerical stiffness caused by the property jumps. In order to maintain a sharp front the fuzzy interface should be restricted to within a few grid cells: the less the grid cell, the better the resolution of the interface. Our research is intended to capture the local physics near the phase interface for liquid-vapor phase change process. We developed a numerical method that is based on fixed grid with explicitly tracked interfaces. At the same time, the explicitly tracked interface is a sharp discontinuity of fluid properties. The current method is a mixed Eulerian-Lagrangian approach. However, the major difference between this method and other mixed Eulerian-Lagrangian methods is the interface treatment The sharp interface is maintained in the present method In the present method first, an explicit interface representing and tracking pro cedure is developed to locate the interface exactly, and a underlying fixed grid is employed to cover the entire domain for field solution Secondly the interface is kept sharp in terms of property jumps. Thirdly two set of governing equations are solved

PAGE 38

28 in two phases, respectively in conjunction with corresponding transport properties. As a result, the fixed grid over the entire domain is virtually divided into two sub domains using cut cell techniques. The coupling of the two phases is treated with the implicit updating of the interface location under the constraints of interfacial conditions. 1.3 Research Objectives 1.3.1 Research Objectives The interests in the vapor bubble phenomena has stimulated continuous, inten sive efforts from numerous investigators on this subject for half a century, as described in the previous sections. And liquid -va por phase change (boiling or condensation) process is one of the most difficult heat transfer modes to investigate. Because the vapor bubble behavior is described by a set of coupled, nonlinear partial differential equations and because there are not mathematical tools to solve those equations an alytically, the majority of the theoretical works on the vapor bubble dynamics in the literature have to resort to more or less simplified theoretical models so as to make the problem tractable. Although those simplified theoretical models sometimes do yield quite good results in comparison with experimental data, it is certainly desired and has been ceaseless attempted to have the theoretical model as close as possible to that set of coupled, nonlinear partial differential equations The present research thus serves to bring that closeness to a level which has not been achieved before, however, the gap between the theoretical model and the ultimate description is not completed removed. On the other hand, owing to experimental difficulties, only a few investigations have focused on the micro-phenomena at the phase change interface. In order to achieve a comprehensive understanding of the physics of boiling and condensation, an effort is needed in the development of a direct numerical simulation technique that

PAGE 39

29 can be used to study the microscopic transport mechanisms at the moving liquid vapor phase change interface. The objective of the present research is thereby to develop a numerical method to simulate the liquid-vapor phase change transport mechanisms of a translating vapor bubble in a superheated or subcooled liquid. Based on the numerical simulation results, insights are offered into the physical and numerical characteristics of the bubble dynamics An assessment of some conclusions or correlations proposed in the literature regarding the bubble growth or collapse dynamics is also presented The development of the present fixed grid, sharp interface moving boundary method consists of the implementations of the following major components: Interface representation algorithms using B-spline curve fitting Fairing algorithm for continuous curvature calculations, Intersection of the curved interface with grid lines Assemble and reshaping of the cut cells, Discretization of governing equations on regular and cut cells Fractional step procedure for solving the coupled governing equations Imposition of the boundary conditions for velocity, pressure and temperature at the interface Iterative algorithm for updating the interface location based on normal stress balance conditions at the interface, Algorithm for the interface evolution owing to phase change effect. This challenging task is accomplished owing to the latest development of numeri cal techniques for solving time-dependent free boundary problems in fluid mechanics

PAGE 40

30 In recent years, it has become possible to solve coupled, nonlinear mass momen tum and energy conservation equations (partial differential equations) numerically because advanced numerical techniques and fast computer hardware were available. For problems with even more physical complexity such as those involving multi-phase fluid flow, heat and mass transfer and/or moving interfaces etc., the reported simu lation works on avier-Stokes equations start to appear in the literature only in the last ten years yet there are many numerically relevant issues which are not resolved In our case, the method for solving the free boundary problems is developed and large scale computations are carried out to solve coupled Navier-Stokes and energy equations in liquid and vapor phase respectively at the same time. On the phase transition interface the mass, momentum and energy conservation principles are matched from the liquid and vapor sides of the interface The method developed in the current research is not restricted to only liquid vapor phase change problems, it is generally applicable to other phase change pro cesses such as solid-liquid phase change of melting or solidification because the method is developed to handle property discontinuity and moving interfaces which are common to the liquid-vapor and solid-liquid phase change. 1.3.2 Scales and Dimensionless Parameters Several major non-dimensional parameters including Reynolds, Weber, Jakob and Fourier numbers are varied to gain a better understanding of physical and nu merical aspects of the solution procedure To determine the above dimensionless pa rameters, characteristic scales are chosen based on the problem features The length and velocity scales are determined as well as the maximum temperature difference 6.T, in the system. Naturally we choose the initial bubble diameter Do as the length scale L. The velocity scale can be defined based on the diffusion mechanism: aif L for stationary bubble growth cases buoyancy effect: (gL ) 1 1 2 for rising bubble

PAGE 41

31 Table 1.3: Saturation properties of water at Psat = 101.3 kPa Tsat = 373 15 K II Saturation Properties II Liquid Vapor p (kg/m 3 ) 958.3 0.597 (N s/m'l) 277.53 X 106 12.55 X 10 6 Cp (J /kg K) 4 220 2,030 k (W/mK) 0 679 0 025 Pr 1.72 1.02 -X (J/kg) 2,256,700 a (N/m) 0.0589 Table 1.4 : Saturation properties of Ethanol at Psat = 101.3 kPa, Tsat = 351.45 K II Saturation Properties II Liquid Vapor II p (kg/m 3 ) 757.0 1.435 (N s/m'l) 428.7 X 10 6 10.4 X 106 Cp (J/kg K) 3000 1830 k (W/mK) 0 1536 0.0199 Pr 8.37 0.96 ,X (J /kg) 963 000 a (N/m) 0.0177 with phase change cases, or the bubble terminal velocity U for bubble rising with no phase change. The choice in each case will be specified individually The rest of the information required to specify the dimensionless parameters are the thermo-physical properties of the liquid An example of working fluids is water whose saturation properties at the atmospheric pressure are listed in Table 1.3. Then the typical magnitude of those dimensionless property ratios of the liquid and vapor phases of water are pt/ Pv = 1605.2 t/ v = 22.114 kt/kv = 27.16, cpzf cP v = 2.0788. The range of Reynolds number we consider in this research is Re = 1 v-. 1000 The bubble diameter is 0 5 mm v-. 1.5 mm, and the maximum temperature difference l:lT = 5C Besides water, we also consider another working fluid, namely, Ethanol (CH 3 CH 2 0H) whose thermo-physical properties are listed in Table 1.4.

PAGE 42

CHAPTER 2 UMERICAL METHOD In this chapter, the fixed-grid, sharp interface moving boundary method is de scribed. In this method, a underlying fixed Cartesian grid serves as the Eulerian framework of the algorithm. Within this Eulerian framework, separate marker points which are connected by piece-wise polynomials are adopted to represent the phase interface. Those regions separated by the phase interface represent different phases. The deformation and movement of the interface is realized through the translation of these markers over the underlying stationary grid, which constitutes the Lagrangian part of the algorithm. In each phase region the fractional step method is employed to solve the coupled governing equations of either liquid or vapor phase. 2.1 Fixed-Grid, Sharp Interface Method The development of this numerical method consists of three major tasks: Implementation of curve fitting and curvature calculation algorithms. Development of Cartesian grid techniques for solving fluid flows over fixed, immersed interfaces with complex geometries. Development of moving interface algorithms, such as the interface updating algorithm The content of this chapter is organized, as shown in Table 2.1, so as to describe t hese three major components of the numerical method in a logic manner. In the section 2.2 of this chapter, the implementation of interface representing tracking and curvature calculation algorithms are detailed. The C 2 cubic B-spline curve fitting 32

PAGE 43

33 Table 2.1: Organization of Chapter 2 Interface representing and tracking algorithms Section 2.2 C 2 cubic B-spline algorithm Section 2.2.1 Fairing algorithm Section 2.2.2 Validation of geometric algorithms Section 2.2.3 Cartesian grid method for fixed, immersed interfaces Section 2.3 Fractional step procedure for solving coupled governing equations Section 2.3.1 Discretization of governing equations on regular and cut cells Section 2.3.2 Imposition of boundary conditions at the interface Section 2.3.3 Algebraic equation solver Section 2.3.4 Moving interface algorithms Section 2.4 Reinitialize the flow field due to phase change Section 2.4.1 Inclusion of interfacial conditions Section 2.4.2 Iterative Process for Determining the Interface Shape Section 2.4.3 Determining Interface Shape with Phase Change Section 2.4.4 Discontinuity of Material Properties Section 2.4.5 is employed in conjunction with fairing algorithm. In the section 2.3, numerical procedures related to solving fluid flows over fixed, immersed interfaces on Cartesian grid are described. In the section 2.4, the numerical techniques for moving interface computations are presented. The overall solution procedure for numerical simulation of moving interface prob lems is illustrated in the flow chart of the code in Figure 2.1 in which the major function modules interface tracking procedure is detailed in Figure 2.2, fractional step process is detailed in Figure 2.3, moving and updating interface is detailed in Figure 2.4. In the rest of the chapter, each section is dedicated to certain numerical issues appeared in these flow charts as also shown in Table 2.1. Currently, the source code of this solver has nearly 40,000 lines of Fortran90 statement implemented on Compaq Tru64 Unix work stations. 2.2 Interface Representation Algorithms A representative schematic of the present fixed grid method is depicted in Figure 2.5 The grid is Cartesian and does not conform to the body surface and the interface

PAGE 44

Begin Input computational parameters Input initial interface shape Interface tracking procedure (Section 2.2) Time step n=O Time step loop Time step n+ 1 Moving and updating interface (Section 2.4) Interface tracking procedure (Section 2.2) Cartesian Grid solver (Section 2.3) Until time step> specified number End 34 Figure 2.1: The fl.ow chart showing the overall solution procedure of moving interface algorithm.

PAGE 45

35 Redistribute the marker points to have equal spacing Parameterization of the interface through marker points Find the intersecitons of the interface with grid lines Identify the phase which each cell lies in and flag cells Creat and assemble the cut cells Assign the material properties Figure 2 2: The flow chart showing the interface tracking procedure as detailed in Section 2.2.

PAGE 46

36 I Initialize the flow field I I Solve for intermediate flow field u (Eqn. 2.12) Imposition of boundary conditions (Section 2.3.3) Calculate face fluxes using interpolation procedure (Section 2.3.2) I Obtain face center velocity u I t Solve for pressure Poisson equation (Eqn. 2.15) Imposition of boundary conditions (Section 2.3.3) Calculate face fluxes using interpolation procedure (Section 2.3.2) Correct cell center velocity field u un+J Correct cell center velocity field u lJ"+ 1 (Eqn. 2.16) I Solve for energy equation 1"+ 1 I Figure 2.3: The flow chart showing the fractional step procedure as detailed in Section 2.3

PAGE 47

37 I Iteration k=O I If Interface update iteration Iteration k=k+ 1 1 I Move interface to new location (Section 2.4.3) I 1 Modify the interface shape to satisfy mass conservation in each phase (Section 2.4.3) I Interface representation procedure {Section 2.2) I Reinitialize the flow field due to phase change (Section 2.4.1) Fractional step procedure (Section 2.3.1) I Terminated if all conditions are satisfied j Figure 2.4: The flow chart showing the interface moving and updating procedure as detailed in Section 2.4.

PAGE 48

38 Vapor Ma ker points Figure 2.5: Schematic of Cartesian grid and interfacial marker points. is defined explicitly by geometric curves in the computational domain. Basic elements involved in defining the interface are interface marker points and interfacial curves connecting the markers. The markers define the terminal points of the interfacial curves. Given a set of markers, finding a curve which fits all markers is a geometric interpolation process. Depending on geometric conditions imposed at the marker points, there are various ways to define numerically the interface characteristics. In the present work, the C 2 cubic interpolatory B-spline curve is employed to fit a set of marker points. The B-spline curve refers to a set of Bezier curves glued together at the marker points to form a composite curve representing the interface. Each piecewise Bezier curve is actually a polynomial curve expressed mathematically for convenience in special form: Bernstein polynomials There are efficient algorithms to do curve manipulation and geometric calculations. er denotes the smoothness con ditions at the junction points of piecewise Bezier curves, requiring that a composite

PAGE 49

39 curve is r times continuously differentiable at the junction points. 0 2 thereb y means the composite curve has continuous second derivatives at any marker point In order to achieve this level of smoothness the individual Bezier curve must be at least degre e 3 (cubic) polynomial to have continuous second derivatives, and second deriv a tives computed from both sides of the junction point have to be equal. This mathematical requirement of the two polynomials yields a smooth global curve to represent the entire interface. Once the B-spline fitting curve of the interface is constructed the geometric informations such as location and curvature of any point along the entire int e rfac e can be easily obtained The translating deforming expanding or shrinking of the interfaces are realized through the motion of each individual marker point which in turn is determined from the flow quantities on underlying fixed grid using, e.g., normal stress balance c ondi tions. With the movement of the markers the instantaneous B-spline representation of the interface is reconstructed accordingly to keep track of the interface. For a interface in motion the curvature calculated based on the continuously constructed B-spline may exhibit numerical oscillations (see Farin, 1997). To extract the correct curvature values along the interface a so-called fa i ring algorithm is adopted. As will be presented the combination of the cubic B-spline and the fairing algorithm r e sults in a robust and accurate method for tracking highly distorted interfaces in terms of location as well as curvature. 2.2.1 0 2 Cubic B-spline Algorithm Th e interfa ce representation using B-spline is based on the marker locations c om puted at every time instant. The interfacial marker points are indexed sequentiall y and can represent any number of open or closed interfaces based on the assign e d connectivity between them. There is no restriction on the number of interfa ce s or on the end conditions that can be imposed on these interface curves. As sho w n in

PAGE 50

40 Phase 1 Phase 0 s=O Phase 1 Figure 2.6: Illustration of immersed interfaces, marker points and normal convention. Figure 2.6, the interfaces can be open or closed, stationary or moving and can enclose different phases or materials with different properties. The locations of the marker points are defined by coordinates X ( s) which is parameterized in terms of arc l e ngth s and a distance ratio based on the grid spacing. In our investigation the marker spacing is initially assigned to be the same as the grid spacing h. In the course of computation, markers are redistributed after each time step according to the initial criterion. The convention adopted by the present algorithm in indexing the marker points is that by choosing the normal of the interface to point from phase 1 to phase 0 th e phase 1 always lies to the right as one traverses the interface along the sequen c e of the marker points. In order to obtain geometric information we need to construct a smooth c urve fitting this set of marker points. In the computer-aided geometric design applications the B-spline offers the greatest flexibility and effectiveness to serve this purpose. The

PAGE 51

41 B-spline curves are composed of piece-wise Bezier curves between any two consecutive markers. When we connect each Bezier curve piece by piece to form the global spline curve, the original data points (here marker points) are junction points of two consecutive Bezier curves. ow the question arises : how is the smoothness of this overall B-spline curve? This question is answered by the definition of e n continuity condition of the B-spline curves. The en continuity means that at any junction point the nth derivativ e computed from the left Bezier curve equals to that computed from right Bezier curve The most popular B-spline curve is the e 2 spline. To have e 2 continuity, each Bezier curve must be at least degree three ; therefore the cubic Bezier curves can be used as piece-wise polynomial in this context. The curve fitting scheme we adopted in this research is e 2 cubic B-spline curve. In the following we describe the construction of e 2 cubic B-spline curves and so-called fairing algorithm which is used to removed numerical noise to recover the accurate information of the curvature of the interface. Based on the spline theory (see Ahlberg et al. 1967 ; Farin, 1997) to construct a e 2 cubic B-spline curve interpolating a set of data points x 0 . XL with c orre sponding parameter values (or knots) u 0 ... uL, the vertices d i of the B-spline control polygon have to be determined first as shown in Figure 2.7. In this example five points Xo , X s are assigned initially. The corresponding knot sequence u 0 , u 5 is chosen as the chord length at each point. The example here illustrates a closed (periodic) curve so x 5 = x 0 The relationship between the data points x i and the control vertices d i is (see Farin 1997) (2.1)

PAGE 52

42 Figure 2.7: Cubic interpolatory spline curve of five points along with its B-spline control polygon Bezier control polygon.

PAGE 53

43 where we have L'.li = L'.lui = U i+i ui and L'.li-2 + L'.li-1 + L'.li L'.li (L'.li-2 + L'.li-1) L'.li-1 (L'.li + L'.li+i) ------+ ----'-----L'.li-2 + L'.li-1 + L'.l i L'.li-1 + L'.li + L'.li+l 'Yi = (L'.li-1) 2 (2.2) Using the periodic condition, we obtain a linear system of the form do ro d1 r1 (2.3) 0:L-2 f3L-2 'YL-2 dL-2 rL-2 'YL-1 0:-1 f3L-1 dL-1 TL-1 where the right-hand sides are of the form The equation Eq. (2.3) can be solved using the procedure described in, e.g., Ahlberg et al. (1967), P. 15 yielding the vertices d 0 , d 5 as shown in Figure 2.7. Then the vertices of control polygon for piece-wise cubic Bezier curves can be obtained based on this B-spline control polygon namely, using the parameter values to determine two vertices on each leg of the B-spline control polygon as shown in Figure 2.7. Each of the piece-wise Bezier curve is defined by four control vertices, including the two marker points, denoted by solid circles, and, between them, two vertices to be decided from the B-spline control polygon. Therefore all four Bezier curve control vertices are (2.4)

PAGE 54

44 where L1 = L1 i_ 2 + L1 i -l + L1 i The index i denotes the ith vertice of B-spline control polygon and the total number of Bezier control vertices is 3i corresponding to vertice i of B-sp l ine control polygon for a closed curve. Once we have the control vertices for each local Bezier curve by defining a local parameter O ::; t ::; 1 for the interval [ui, Ui+1] as t = u;:~~~ ; we can express the piecewise cubic Bezier curve in the form of 3 b (t) = L b i Bf (t) i = O where Bernstein polynomial Bf (t) is defined by with binomial coefficient (~)= { 3! i!(3 i)! 0 if O ::; i ::; 3 else (2.5) (2 6) Using Eq. (2 5), the geometric information such as normal and curvature etc can be eva l uated using analytical formulas as K, = Yt (x; + y;)l/2 Xt (x; + y;)l/2 XttYt YttXt (x; + y;)3 / 2 (2.7) For axisymmetric geometries the total curvature is the sum of,.,, in Eq (2.7) and the other principle curvature i.e xtf y(x; + y;) 3 1 2 2 2.2 Fairin g A lgorith m From the numerical simulation point of view, the smoothness of the interfacial curve and the smoothness of curvatures are two most important characteristics of the interface. Using the B-spline fitting, while a curve looks apparently smooth the curvatures obtained by Eq. (2. 7) can be contaminated by numerical noises. The

PAGE 55

45 Figure 2.8: A sample interfacial curve from one simulation case curvature formula involves the second derivatives as well as nonlinear products of the first derivatives, and is prone to suffer from numerical noises even if the interface shape and slope are found to be smooth. This usually occurs when the interface is evolving in the course of computation. As an example to illustrate the problem Figure 2.8 shows an interface curve taken from a simulation conducted in this study. The corresponding curvature plot based on B-spline fitting process is shown in Figure 2.9. There are substantial oscillations in the curvature profile, indicating that noises associated with the numerical procedures are substantial. Such phenomena are well known in computer-aided geometric design (see Farin, 1997). To treat this difficulty, the so-called curve fairing (smoothing) algorithms have been developed in the literature (see Sapidis & Farin, 1990; Farin, 1997). One such algorithm, in the C 2 cubic spline case, is to make local spline curve segment three times differentiable which is one order higher than that required by the original cubic spline. The basic idea is to adjust the vertex locations to make the local

PAGE 56

46 1 0 8 6 4 2 ::::, cu 0 ::::, 2 0 4 6 8 -1 0 0 40 80 120 160 Marker Point Figure 2.9: Curvature plot obtained from B-spline fitting for the interfacial curve in Figure 2.8 curve segments around them become three times differentiable. Since the curvature involves only first and second derivatives, this extra differentiable requirement will ensure that the curvature be differentiable, and prevents discontinuity in the slope of the curvature, as shown in Figure 2.10. The formula for obtaining new B-spline control vertex d i are briefly listed as d ~ (uj+2 ui) l i + (ui Uj_ 2 )ri 1Uj+2 Uj-2 where the auxiliary points lj ri are given by (ui+1 Uj_3) d j-1 (ui+l ui)di-2 Uj Uj-3 ( Uj+3 Uj-1) d j+1 ( Uj Uj-1)dj+2 Uj+3 Uj (2.8) (2.9) The geometry interpretation is illustrated in Figure 2.10 which is from page 365 of Farin (1997). The detailed discussions can be found in Sapidis & Farin (1990); Farin (1997).

PAGE 57

47 ..... d J Figure 2.10: Knot removal: if di is moved to dj, the new curve denoted by dotted line, is three times differentiable at Uj [page 365, Farin (1997)). In practice typically one fairing operation is not sufficient, and the fairing pro cedure will be repeated multiple times. The criterion adopted in the present study for determining the convergence of the curvature profile is that the maximum difference of the curvature values among all marker points between two consecutive fairing op erations must be less than 10 3 For the curve in Figure 2.8, the resulting curvature plot after 100 iterations of fairing operation is shown in Figure 2.11. The maximum correction in the curvature is within 10 3 It is noted that the fairing algorithm is a geometric operation. It obtains correct curvatures by modifying the geometry of the interface itself, not by manipulating the formulas for computing the curvature of a given geometry. On the other hand, attention needs to be paid to ensure that the interfacial marker locations as well as the volume/surface enclosed by the interface be satisfactorily preserved. A critical criterion of developing a satisfactory fairing algorithm is that with arbitrary number of fairing iterations the geometric information can be maintained at an asymptot ically constant state without being continuously smeared. The corrections in the interfacial locations are usually very small so the corrected marker locations and the

PAGE 58

10 8 6 4 2 :::, cu 0 :::, -2 (.) 4 6 -8 -10 0 40 80 Marker Point curvature before fairing curvature after fairing 120 160 48 Figure 2 11: Curvature plot obtained before and after fa i ring operation for the curve in Figure 2.8. resu l ted interface shape will not have perceptib l e changes after fairing operation as will be demonstrated in the testing of geometric algorithms. 2.2.3 Validation o f Geome t ric Algorithms The geometric algorithms for interface representing and tracking describ e d in the section 2.2.1 and 2.2.2 are capab l e of handling closed and open interfaces. They are validated here for accuracy on some cases used in Shyy (1994) First we examine the expansion of a circle on a 162 x 162 grid as in Figure 2.12. Initially, the circle has a radius of 1 and expands at a constant speed. The curvature of the circle at an y given time is a constant. Figure 2 12(a) shows the interfa c ial shapes at equal intervals of time. Figure 2 12(b) shows the curvature computed by the aforementioned formulas, at the corresponding time instants as in Figure 2.12(a). As shown there the interfacial curvature is a constant while, by regula t ing the spacing between two neighboring markers, the number of markers increases as the

PAGE 59

49 circle expands The results shown in Figure 2 12 are obtained from B-spline fitting algorithm only. o fairing operation is needed for this case; by applying fairing, no impact is observed, either. As another example for variable curvatures, Figure 2.13 shows an initial interfa cial curve shape expressed in the form of~= 2.0+(1-cos(21rXi/lO)) The curvature plot is found to be smooth using B-spline fitting for the initial shape as shown in Figure 2 13(b). Given the normal velocity of each marker point as Vn i = Y 1 the interface moves as depicted in Figure 2.16(a). Figure 2.14 shows the interface shape and corresponding curvatures at the early stage of the interface evolution. If the fairing is not applied, the curvature obtained from B-spin representation in Figure 2.14(b), exhibits oscillations near both ends of the curve after the interface evolves to a new shape. After 10 fairing operations, the shape and curvatures at the same time instant are shown in Figure 2 15. The curvature profile in Figure 2.15(b) is free from the artificial spikes at the ends. The curvature plot in Figure 2.15(b) becomes smooth and the shape in Figure 2.15(a) is well preserved and can't be distinguished from that in Figure 2.14(a) as if the shape is iden t ical to previous one. This is further confirmed by plots in Figure 2.18 Figure 2.18 shows the corrections on coordinates of markers by the fairing operation at various time instants from the beginning to the end of the development of the interface All corrections of coordinates are small, which is the reason why the two shapes before and after the fairing operations are virtually identical to each other. In Figure 2.16 and Figure 2.17, the shape and the curvatures at the end of the development of the interface are plotted while the former shows the situation before fairing is applied the latter illustrates the effect after the fairing. Furthermore with a periodic curve, the curvature should be periodic as well As shown in Figure 2.13(b) initially, the computed curvature is periodic; however, without fairing, Figure 2 16(b)

PAGE 60

50 (a) > 0 0 2 4 6 8 10 X (b) 0.8 0.7 0.6 .. e :::::, .... 0 5 ca 0.4"'" :::::, CJ 0.3 0.2 0.1 I I I I I 0 100 200 300 400 500 Marker Point Figure 2.12: (a) Shapes of expand i ng circle at equal interva l s of time. (b) Corresponding curvatures along the curve at the same time instants as in (a).

PAGE 61

51 shows that the computed curvature for a moved interface is no longer periodic. After fairing, as shown in Figure 2.17(b) the periodicity is restored again The shape is well preserved and the oscillation in curvature is removed. The majority of the curvature plot after fairing is the same as that before fairing The number of marker points increases as a result of reorganizing process to follow the increasing arc length of the bulging out interface The B-splin e curve fitting in conjunction with the fairing operation is a robust interface representing and tracking method which can yield smooth interfac e and smooth curvatures while preserve the interface shape. The programming t ask how ever, is quite challenging. For closed and open curves, different program modules and algorithms have t o be implemented for B-spline fitting mainly because the end con ditions are different. For closed curves, the end condition is always periodic (c y clic). In the cases of open curves there might be other end conditions 2.3 Cartesian Grid Method for Sharp Interfaces Once the interface is defined, one needs to identify in which phase each compu tational cell lies so that correct transport properties can be assigned. Furthermore for any two neighboring cells between which a interface passes through, there may be a discontinuity in transport properties such as density and viscosity. Those boundary cells are reshaped to maintain flux conservation around the interface A major goal of the present approach is to adopt a finite-volume formulation for all computational cells so that mass momentum and energy conservation is honored in all resolvable scale, including the computational cells intersecting with phase boundaries Once a cell intersects with an interface, it is split into two parts; with the partial cell con taining the center of the original Cartesian cell maintaining the initial cell index and the other merged into a neighboring cell belong to the same phase. As illus t rated in Fig. 2.19, this procedure results in irregular trapezoidal shaped cells around the interfaces. Awa y from the interfaces, most cells are still structured Cartesian cells

PAGE 62

52 (a) 8 6 > 4 2 0 0 2 4 6 8 10 X (b) 2 1.5 1 0.5 ::J ... (U 0 ::J 0.5 0 -1 -1.5 -2 0 70 140 210 Marker Point Figure 2.13: (a) The initial shape of two fingers. (b) The corresponding curvature of two fingers using B-sp l ine fitting

PAGE 63

53 (a) 8 6 > 4 2 0 0 2 4 6 8 10 X (b) 1.5 1 0.5 0 ::l ... ca -0.5 C: ::l -1 0 -1.5 -2 -2.5 0 80 160 240 Marker Point Figure 2 14: (a) The shape of two fingers at the early stage of the development when fairing is not applied. (b) The curvature corresponding to the shape denoted by solid line in (a).

PAGE 64

54 (a) 10 --------------, 8 6 > 4 ... .. ... 2 0 0 2 4 6 8 10 X (b) 1.5 1 0.5 0 ::, .., ro -0.5 ::, -1 0 -1.5 -2 10 iterations of fairing -2.5 0 80 160 240 Marker Point Figure 2.15: (a) The shape of two fingers at the early stage of the development when fairing is applied. (b) The curvature corresponding to the shape denoted by solid line in (a).

PAGE 65

55 (a) 10 ----------------, 8 6 > 4 ... .. ... \ . . ii i~t)ii~1i~\/ ".;,;j:-~ 2 0 .,__ __ .__ __ .__ ___ ________ 0 2 4 6 8 10 X (b) 0 ::::, 2 ... ::::, -4 (.) -6 -8 --------------------0 80 160 240 320 400 480 Marker Point Figure 2.16: (a) The shape evolution of two fingers until the later stage of the de velopment when fairing is not applied. (b) The curvature corresponding to the final shape denoted by solid line in (a).

PAGE 66

> (a) 8 6 4 2 0 -----------------0 2 4 6 8 10 X (b) 1----------------, 0 -1 -2 -3 -4 10 iterations of fairing 5 ---------------0 80 160 240 320 400 480 Marker Point 56 Figure 2 17: (a) The shape evolution of two fingers until the later stage of the devel opment when fairing is applied. (b) The curvature corresponding to the final shape denoted by solid line in (a).

PAGE 67

57 (a) 4.0E-03 Q) t,s C: 3.0E 03 "E 0 0 (.) >< 2.0E-03 0 Q) 0 1.0E-03 C: e 0 0E+00 C l Nh M 11 J ) ) / / 0 80 160 240 320 400 480 Marker Point (b) 5.0E-03 Q) t,s C: 4 0E-03 "E 0 0 3 0E-03 (.) > 2.0E-03 0 Q) 0 C: 1 0E-03 e j j Q) Vj I~) :E 0 0E+00 '" /),._ /J J C 0 80 160 240 320 400 480 Marker Point F i gure 2.18: (a) The corrections of the i nterface location by fa i ring algorithm at time instants corresponding to shapes shown in 2 17(a) (b) The corrections of the interface location by fairing algorithm at the same markers as in (a).

PAGE 68

58 Figure 2.19: Illustration of the resulting situation when cut-cell approach is applied i.e. fragments of cells which are cut by the interface are absorbed into neighboring cells. The newly formed cells are shown in dashed lines on both sides of the int e rface In such an algorithm while a nominal structured grid index system is maintained the flux calculation across the surface is conducted based on interpolation of varying degrees of polynomials. The resulting flow solver needs to account for both Car t esian and trapezoid cells using a finite volume fractional step method. After the domain is partitioned into several sub-domains representing diff e ren t phases. The so-called Cartesian grid method is developed to deal with the discr e tiza tion and solution of governing equations in such a domain left with mostly structured rectangular cells and some unstructured trapezoidal shape cells in the vicinit y o f the interface Essentially what the Cartesian grid method does is to solve the gov e rning equations on a domain with fixed, complex geometry interfaces without resorting to the body-fitted grid. Although the method is developed in this research as one of the major components in the entire moving interface simulation package the m e thod

PAGE 69

59 alone is an good alternative to the body-fitted grid approach to solve the flow over solid bodies with complex geometry, as demonstrated in Ye et al (1999). This globally second-order accurate Cartesian grid flow solver for fixed immersed interfaces (see Ye et al., 1999) is summarized in the following : 1. The interface is represented numerically using C 2 cubic B-spline interpolatory curve ; 2. A fractional step procedure (see Charin, 1968; Kim & Moin 1985) is employed to solve th e coupled governing equations, which results in a efficient solution of time accurate unsteady flows; 3. Colocated variable arrangement (see Ferziger & Perie 1996) is adopted which leads to simplicity in coping with cut cells near the interface, 4. Second-order accurate central difference scheme is used for spatial discretiza tion, Crankicolson scheme for temporal discretization; 5 A compact interpolation scheme near the immersed interfaces is devised that allows us to retain second-order accuracy and conservation property of the solver. 6. A precondi t ioned conjugate gradient iteration method is used for solving the Poisson equation for pressure, which takes advantage of the underlying struc tured nature of the mesh and also substantially accelerates the convergence of the Poisson equation for pressure. 2.3.1 Fractional Step Method for Unsteady Transient Solution To use finite volume discretization the dimensionless governing equations listed in the section 1.1.7 on page 9 have to be written in the integral forms by integrating each term over the control volume c v and wherever possible using Gauss th e orem

PAGE 70

60 to convert volume integrals into surface integrals. To solve the momentum equations (1.18) and (1.21) under the constraints (1.17) and (1.20) in the liquid and vapor phases respectively, a so called fractional step procedure is applied. After the mo mentum equations are solved, the energy equations (1.19) and (1.22) can then be solved in the updated flow field. All governing equations are discretized on a Cartesian mesh using a colocated arrangement (see Ferziger & Perie, 1996) of the primitive variables which are located at the cell center. For convenience of describing the numerical procedure, only the integral form of dimensionless governing equations in the liquid phase are considered throughout this chapter as examples, the vapor phase transport equations are solved in the exactly same procedure except additional dimensionless fluid property ratios are present on those cells belong to the vapor phase. The integral form of Na vier-Stokes equations in the liquid phase are: Mass conservation: ju ndS = 0 (2.10) cs Momentum conservation: j ~: dV + ju ( u n) dS = j pn dS + ~e j \l u n dS (2.11) CV cs cs cs In the above equations cv and cs denote the control volume and control surface respectively and n is a unit vector normal to the control surface. The above equations are to be solved with u(x, t) = v(x t) on the boundary of the flow domain where v is the prescribed boundary velocity. A second order accurate, two step fractional step method (see Chorin, 1968; Kim & Moin, 1985; Zang et al., 1994) is used for advancing the solution in time. In this time stepping scheme, the solution is advanced from time step 'n' to n+ 1' through an intermediate advection diffusion step where the momentum equations without the

PAGE 71

61 pressure gradient terms are first advanced in time. The semi-discrete form of this advection diffusion equation for each cell shown in Figure 2.20 can be written as: cs +1 -/ (\7u* + \7u n) n dS 2Re cs (2. 12) where u* is the intermediate cell-center velocity. A second-order Adams-Bashforth scheme is employed for the convective terms and the diffusion terms are discretized using the implicit Crankicolson scheme. This eliminates the viscous stability con straint which can be quite restrictive in the simulation of viscous flows. The velocity boundary condition imposed at this intermediate step corresponds to that at the end of the full time step, i.e., u* = vn +i. At this stage, in addition to the cell center velocities which are denoted by u, we also introduce face center velocities U In a manner similar to a fully staggered arrangement, only the component normal to the cell face is computed and stored (see Figure 2.20). The face center velocity is used for computing the surface flux from each cell in our finite volume discretization scheme. The reason of separately computing the face center velocities will be addressed later in this section. Following the advection diffusion step, the intermediate face center velocity U is computed by interpolating from the cell center intermediate velocities. The advection diffusion step is followed by the pressure correction step u n+l u* ____ = -\7pn+l ~t (2 .13) where we require that the final velocity field satisfy the integral mass conservation equation given by J (un + 1 n) dS = 0 (2.14) cs

PAGE 72

: : N -!.------,------(!)--..... _,_ : I Uy, I I I I I ~WW I -w I I -i--------1 I u_,, I pU.r Ux E EE I I s I I --(f)-- 'I I 8 I 0 5 -point stencil for regular cell 62 Figure 2.20: Schematic of underlying Cartesian mesh and arrangement of cell-center and face-center velocities This results in the following equation for pressure J (Vp) ndS = ~t J (U*n) dS (2.15) cs cs which is the integral version of the Poisson equation for pressure. ote that the pressure correction step is represented by the inviscid equation (2.13) and is well posed only if the velocity component normal to the boundary is specified. Therefore the velocity boundary condition consistent with Eq. (2.13) is u n+l. n = vn+ 1 .n where n is the unit normal to the boundary of the flow domain. It can be easily shown that this implies that (~pn+I) n = 0 be used as the boundary condition for equation (2.15). Once the pressure is obtained by solving this equation, both the cell center and face center velocities are updated separately as follows: un+i = U ~t (Vpn+i) Jc (2.16)

PAGE 73

63 It is a well known fact that in a colocated mesh, when the compact stencil is used for the Poisson equation for pressure, the pressure field usually does not exhibit non physical oscillation however, the final velocity field does not satisfy the divergence free constraint exactly To overcome this problem two types of approaches could be used i.e the momentum interpolation (see Rhie & Chow, 1983) scheme in which a pressure gradient term is introduced into the interpolation process for obtaining cell face velocities. As it is shown by Ferziger & Perie (1996), this technique produces a third-order dissipation in estimating the cell face velocities. The momentum interpo lation technique is not straightforward to implement in cut cells encountered in the mixed Eulerian-Lagrangian algorithm for flow with immersed interfaces. The other approach for dealing with this problem was first used by Zang et al (1994). In their method, both cell center and face center velocities are stored while momentum equations and Poisson equation for pressure are solved at the cell centers only. After Poisson equation for pressure is solved, both cell center and face center velocities are updated The pressure gradient at the cell face is computed in a compact stencil manner. They have used this procedure in conjunction with a curvilinear mesh solver to simulate the turbulent flows with LES. In our algorithm, we implemented this method and test results show that this approach is we ll suited for high Reyno l ds number, unsteady flows. 2 3.2 Dis c r etiz a ti o n o n R eg ular an d C u t Ce ll s To discretize governing equations in the configuration as shown in Figure 2.19, we have to deal with two types of situations: discretization on cells which are away from the interface thus not cut by the interface and on those cells which are cut by the interface and reshaped to become trapezoidal cells. The numerical procedure for solving avier-Stokes equations just outlined can be easily implemented on a regular Cartesian grid cell. In this section, detailed finite volume integration of each term is described. For demonstration purpose basic

PAGE 74

64 (a) (b) liui -.-r ... ... I (c) (d) Figure 2.21: Schematic of cells around the interface. (a) boundary cells with interface located south of cell center. (b) boundary cells with interface located west of cell center. ( c) typical reshaped quadrilateral cells corresponding to case (a). ( d) typical reshaped quadrilateral cells corresponding to case (b). procedure for evaluating volume and surf ace integrals on a regular cell is explained first, followed by the special treatment on those cells cut by the immersed interfaces and reshaped, as shown in Figure 2.19. The emphasis is placed on the latter. In general, the so called cut cell approach is used to treat those reshaped cells: the small fragments of cut cells are absorbed into neighboring cells in the same phase to form mainly trapezoidal cells as shown in Figure 2 21, after that the governing equations are integrated on those reshaped cells. It is straightforward to obtain second order accurate spatial discretization on regular cells. The key for retaining global second-order accuracy of the solver is to evaluate mass convective and diffusive fluxes and pressure gradients to second order accuracy on the cell face of those trapezoidal cells resulted from reshaping process.

PAGE 75

65 As shown in equation (2.12), a finite volume discretization of Navier-Stokes equations requires the estimation of surface integrals on the faces of each cell. The integrand ( denoted here by f can either involve the value or the normal derivative of a variable. An example of the former is the convective flux denoted by (p
PAGE 76

66 (a) (b) ltqui O Points in the six-point stencil for fsw Figure 2.22: Schematic of interpolation for cell face values and derivatives at trape zoidal cells (a) various fluxes required for trapezoidal cell (b) trapezoidal region and stencil used in computing fsw neighboring cell center values. The current interpolation scheme coupled with the finite volume formulation guarantees that the accuracy and conservation property of the underlying algorithm are retained even in the presence of curved immersed inter faces. In the following, we describe the interpolation function for a typical trapezoidal boundary cell. Consider the trapezoidal cell ABCV in Figure 2.22a. The face ABC of the trapezoidal cell is composed of two pieces: AB coming from cell P and BC coming from cell S. The integral on this face can be decomposed as J f dy = J f dy + J f dy AC AB BC A second order approximation to this integral can then be obtained as J f dy fw (YA YB)+ f s w (YB Ye) AC ( 2.17) ( 2.18) where !wand sw are computed at the centers of segments AB and BC respectivel y On the other hand, if the face is cut by the immersed interface such that it is smaller than

PAGE 77

67 a nominal cell face, as in the case of face VE then the integral can be approximated as J l dy -;::::, le (Y Yv) V (2.19) where le is the flux computed at the center of the segment VE. For non-boundary cells, these face center values can be evaluated to second order accuracy quite easily by a linear approximation and we would therefore like to evaluate lw, lsw, le to within second order accuracy also. Approximation of lw to second order accuracy is quite straightforward and is done in the same way as for the face of a non-boundary cell. For instance, if lw requires the value of this can be evaluated to second order accuracy as (2.20) where the linear interpolation factor >-w is defined as Xp Xw AWXp-Xw (2.21) Alternatively, if lw requires the normal gradient of as it would for the diffusion or pressure gradient terms, this can be approximated by a central difference scheme as follows:
PAGE 78

68 a second order accurate scheme can be constructed since le is not located on the line joining the neighboring cell centers and consequently, expressions such as Eq. (2.20) or Eq (2.22) cannot approximate this flux to second order accuracy. Thus, a different approach is needed here for evaluating these fluxes. Our approach is to express the flow variables in terms of a two dimensional polynomial interpolating function in an appropriate region and evaluate the fluxes such as lsw or le based on this interpolating function. For instance, in order to approximate lsw, we express
PAGE 79

69 evaluation of derivative at this location On the other hand, this situation does not exist in the y-direction for the cell shown in Figure 2 22b and therefore a quadratic interpolation is necessary in this direction in order to obtain a second-order accurate approximation to (8/8x) at the center of BC. In Appendix 1, we have demonstrated numerically that the linear quadratic interpolating function shown in Eq. (2.23) does indeed result in second order accurate evaluation of values and derivatives on a line which is located midway between the two parallel sides of a trapezoid. It can be seen in Figure 2.22b that the sides of the trapezoid in which the interpolation is performed pass through four nodal points and two boundary points. Thus, the six unknown coefficients in 2.23 can be expressed in terms of the values of at these six locations. To solve for en, we obtain the following system of equation by expressing the at the six location in terms of the linear quadratic interpolating function: [ x,yl Y? X1Y1 Y1 X1 )-] C1 X2Y2 y~ X2Y2 Y2 X2 C2 (2.25) 2 Yl X6Y6 Y6 X6 X6Y6 The coefficients can now be expressed in terms of values of at the six points by i nverting Eq. (2.25) i.e. 6 Cn = L bnjj j=l (2.27) (2.28)

PAGE 80

70 where (2.29) The value of (8/ox) at the center of BC is expressed as ( 8) 2 OX sw = C1Ysw + CJYsw + Cs (2.30) and using Eq. (2.26) this can be written as ( 8) 6 ax = I: 13 jj SW j=l (2.31) where (2.32) Note that a and /3 are coefficients that depend only on the mesh, the location and orientation of the immersed interface. Therefore these can be computed once at the beginning of the solution procedure. Subsequently, relationships such as (2.28) and (2.31) can be used in the spatial discretization of the governing equations (2.12) (2.16). A similar interpolation procedure is also used for approximating leFor this, a linear quadratic interpolating function is used in the trapezoidal region shown in Figure 2.23 and a relationship similar to Eq. (2.28) and (2.31) is developed for ap proximating leThe six points contained in this stencil are also shown in Figure 2.23. It should be pointed out that the north face of the particular cell being con sidered here does not need special treatment since face center values and derivatives can be computed to second order accuracy using a linear approximation. However, in general there are also boundary cells which have their north or south faces cut by the immersed interface (as shown in Figure 2.21b). For these boundary cells too, the same approach is used to evaluate the fluxes on these cut faces. The only difference here is that the interpolating function is linear in y and quadratic in x.

PAGE 81

71 I ,vapo I 0 Points in the six-point stencil for fe Figure 2.23: Trapezoidal region and stencil used in computing feNow we turn to the calculation of the flux on cell face CV which lies on the immersed interface as shown in Figure 2.22a. The integrated flux on this face can again be evaluated to second order accuracy using the midpoint rule and as before we would like to evaluate the integrand at the center of face CV ( denoted here by !int) to second order accuracy. In general both convective and diffusive fluxes are needed on this face and this requires approximation of variable values as well normal derivatives at the center of CV. The value is usually available from a Dirichlet type boundary condition and hence no interpolation is required for this. Here we describe the approximation procedure for the normal derivative. The normal derivative on face CV can be decomposed as (2.33) where nx and ny are the two components of the unit vector normal to face CV. Since we know the shape of the immersed interface, nx and ny are known. Therefore computation of the normal flux requires estimation of 8/ ax and 8/ 8y at the center

PAGE 82

72 (a) ( b ) - va or va or 0 P o ints in the stencil. Figure 2.24: Stencil for calculating interface flux. (a) stencil for calculating 8/oy (b) stencil for calculating 8 / ox of the line segment CV For the cell being considered here, 8/oy is computed to second order accuracy with relative ease by expressing the


PAGE 83

73 the trapezoid shown in Figure 2.24b Again because the x-coordinate of the c enter of CV is located midway between the two parallel sides of this trapezoid expressing the variable in this trapezoid in terms of an interpolating function which is linear in x and quadratic in y allows us to obtain a second order accurate approximation to ( 8 / 8x \w at the center of the line segment CV. The procedure for this follows along lines similar to that shown earlier for ( 8 / 8x \w and we get the following expr e ssion for the x-derivative on the interface: ( 2.36) where r f depend on the location and orientation of the immersed interface in the neighborhood of the cell under consideration. Finally we obtain an expression of the form ( 8) 9 = LTj
PAGE 84

74 interface are covered by the grid, we have the capability to solve a different set of equations on each side of the immersed interface. For instance the multi-phase flow with separate governing equations in different phases across the interface could be simulated in this manner. 2.3.3 Imposition of Boundary Conditions on the Interface This section describes the imposition of boundary conditions on curved immersed interfaces on the underlying Cartesian grid In the discretization process on cut cells as discussed in the last section the interpolation scheme used to obtain the face flux such as fe, !sw and !int as shown in Figure 2.22 involves the values of generic variable at point (5), (6) on the int erface as illustrated in Figure 2.22, 2.23 and 2.24. The values of on the interface are simply known if Dirichlet boundary conditions are specified on the interface. In the case of Neumann boundary conditions on the interface, the explicit values of on the interface need to be extrapolated from values in the interior cells. The extrapolation process is explained in Figure 2.25. For implementing Neu mann boundary cond itions on the interface, it is necessary to obtain the expression for gradients on the interface. Taking one side of the interface as example as shown in Figure 2.25, the gradient can be evaluated by using a normal probe, i.e., drawing a normal into one phase from the interfacial marker point location. Since second order accuracy is desired, two nodes on the normal probe are needed which are located at distance h and 2h from the interface. Values at each node, ni n 2 are obtained by bilinear interpolation from the neighboring grid points which are denoted by shaded circles. Then based on the values at two points on the normal probe and the known gradient values on the interface, an O(h 2 ) estimate of the normal gradient at the marker point k location can be written as: 4nl c/Jn2 3c/Ja 2h (2.38)

PAGE 85

2nd normal probe node
PAGE 86

76 Once s 5 6 is obtained based on the parameterization of in terms of s the value of
PAGE 87

77 liquid ' .. ~3 -@@ ~ I vapor 0 Points in the six-point stencil Figure 2.26: Six points stencil of finite volume integration on a typical boundary cell matrix corresponding to the trapezoidal boundary cells also have additional elements since the stencil of the trapezoidal boundary cell is different from regular cells The alternative direction line successive-overrelaxation (SOR) method (see Press et al. 1992) is one of the most widely used iterative methods for solving equations resulting from discretization on structured grids. A typical implementation of this technique, which we have adopted here, involves alternating sweeps along the two families of grid-lines. We find that even in the presence of the immersed interfaces, this method is extremely effective for the numerical solution of the discretized advection-diffusion equation and the residual can be reduced to an acceptable level within a few itera tions. The discretized Poisson equation for pressure however exhibits a slower conver gence than the advection-diffusion equation This is because the time-derivative term in the advection-diffusion equation tends to improve the diagonal dominance of the corresponding discretized operator. In fact due to its slow convergence, the solution

PAGE 88

78 of the discretized Poisson equation for pressure is usually the most time consum ing part of a fractional-step algorithm. In the presence of immersed interfaces this behavior of the Poisson equation for pressure can be further exacerbated since the stencil for the trapezoidal cells contains dependence on some neighboring cells which are not included in the line-SOR sweeps. For example, in Figure 2.26 the coupling of cell P with cells (1) and (4) in the stencil is not included directly in any of the line sweeps. Furthermore, depending on the aspect ratio of the trapezoidal boundary cell and the angle at which the immersed interface cuts the cell, diagonal dominance in the pressure operator can be severely weakened. In the various simulations that have been performed using the current method, it has been found that the simple line-SOR procedure can result in an extremely slow convergence for the pressure equation. One remedy in such a situation is to resort to a multigrid method (see Brandt 1977; Briggs, 1987). However, the presence of the immersed interface can complicate the implementation of a multigrid procedure since operations such as prolongation and restriction are difficult to perform near the immersed interface. In contrast Krylov subspace methods (see Golub & Loan, 1989) are an attractive alternative since these are designed for general sparse matrices and therefore do not assume any structure in the matrix operator. Thus, the presence of the immersed interface does not pose any additional complication for the implementation of these methods. A particularly suitable method in this class is the bi-conjugate gradient stabilized (Bi CGSTAB) method (see van der Vorst, 1992; Barrett et al., 1993) which is applicable to non-symmetric matrices and provides relatively uniform convergence. However, as with all conjugate gradient methods, the convergence rate of the Bi-CGSTAB proce dure depends critically on the choice of the preconditioner. The Jacobi preconditioner which has a trivial construction phase is used routinely in conjunction with unstruc t ured grids. However, this preconditioner only produces marginal improvement in t he convergence rate of conjugate gradient type algorithms. Other preconditioners

PAGE 89

79 such as those based on incomplete factorization can substantially increase the con vergence rate but these usually have a non-trivial and expensive construction phase (see Barrett et al ., 1993). The structure of the matrix that results from the current Cartesian grid ap proach however presents us with another alternative choice of a preconditioner. As pointed out before the presence of an immersed interface only alters the underlying matrix operator in the rows corresponding to the boundary cells. Although this alter ation slows the convergence rate of the line-SOR procedure, the convergence rate is still significantly faster than what is obtained with a simple point Jacobi method. It follows that the line-SOR procedure would serve as a better preconditioner than the Jacobi preconditioner. A further advantage of using the line-SOR as a preconditioner is that this procedure only requires the solution of tridiagonal systems and this can be accomplished with ease using the Thomas algorithm. Thus in the solver developed here, the line-SOR procedure was used as a preconditioner in the Bi-CGSTAB algo rithm and a significant improvement in the convergence rate over a simple line-SOR iterative procedure was found. This completes the discretization and solution procedure of Navier-Stokes equa tions in a domain with immersed stationary interfaces The same solution proc e dure can be applied to solve energy equations considering it is an advection-diffusion e qua tion similar to momentum equations. When the immersed interfaces are no longer stationary due to various ph y sical mechanisms such as phase change or hydrodynamic forces effect the situation be comes much mor e complicated. Many issues associated with solving moving interface problems arise If the discontinuity of fluid properties become appreciable and e ven severe the overall numerical procedure will face additional challenges In the following sections, key issues regarding solving transport equations for problems with moving interfaces are addressed.

PAGE 90

80 2.4 Moving In t erface Algorithms In the previous chapters, the numerical procedure used to solve the coupled nonlinear transportation equations in the domain with immersed fixed interfaces is described The method is fixed-grid approach in which the grid lines don t conform to the interface with complex geometry. When the immersed interfaces in the flow domain translate over the underlying fixed Cartesian grid due to physical mechanisms, various computational issues related to the moving interface can arise. Overall, the numerical procedure involving immersed moving interfaces on a fixed Cartesian grid is to solve a fixed interface problem at any time instant provided that the instantaneous interface location is known However since the interface movement is not independent of the flow field this coupling between the interface evolution and the flow field makes the moving interface problem not a simple addition of many individual fixed interface problems because determining the interface location at any time step is also part of the solution process. For fixed interface problems the global iteration process is only carried out among the coupled mass momentum and energy equations to solve for the flow field. With the introduction of moving interfaces the interface location is coupled with the flow field, another outer iteration is needed to deal with this higher level of coupling. The major difficulty of the entire simulation methodology arises in this iteration process to solve for the flow field and to adjust the interface shape together. Physically the interface location is to be determined based on, in the case of liquid-vapor two phase flow considered here the thermal effects such as evaporation or condensation of liquid or vapor as well as hydrodynamic forces such as pressure viscous stresses in the two phases The interface must be advected to a position where multiple interfacial conditions such as mass, momentum and energ y conservation conditions (see Equations 1.23, 1.24 1.25) and mass conservation constraint within

PAGE 91

81 each phase are simultaneously satisfied. In order to achieve this a major num e rical technique is required to determine the interface location in a iterative manner within any particular time step. 2.4.1 Change of Phase As it is mentioned earlier in this chapter when a interface which repr e sents the phase front of the liquid and vapor state moves across the underlying stationary Cartesian grid some grid points may undergo discontinuous change of the phas e due to the interface motion as illustrated in the Figure 2 27 At the time level n, the interface is located such that the grid point ( i, j) lies in the vapor phase as shown in Figure 2.27a. At the next time level n+ l, the interface moves to the updated location shown in Figure 2.27b as can be seen, some grid points which previously lie in the vapor phase now lie in the liquid phase. In the time advancement solution proc e dure for the transportation equations, the value of any variable at previous time l e vel n is needed at current time level n + l by the temporal discretization scheme adopted. However for grid point that have just changed the phase since last time level n such as ( i j) in the Figure 2.27 the current computation at time level n + l is performed in the liquid phase while previous history for this point is in the vapor phas e which couldn t be used i.e there is no physically meaningful information regarding the previous value of q> n at ( i, j) This is a purely computational issue encountered in the numerical approach that treats the interface as sharp discontinuity in terms of the fluid properties In the purely Eulerian methods the interface is not explicitly tracked but deduced from certain field variables which are smoothed over the fixed grid thus the interface is spread over a band of grid points rather a sharp discontinuity. In that case, the problem of change of the phase is not a significant issue if it is a issu e at all. However, if the sharp discontinuity of the interface is desired to be maintained like in our approach those grid points have to be dealt with

PAGE 92

82 A solution to this problem is to estimate the value of any variable at the grid point ( i, j) in the Figure 2.27 by interpolation. At time level n + l, the grid point ( i, j) falls into liquid phase, but this point was in the vapor phase at last time level n. Therefore the value of ".!,j doesn't exist in terms of the liquid phase and is obtained through interpolation between the surrounding grid points in the liquid phase and points on the interface. The procedure of the interpolation is as following. First a normal to the interface, which passes through the grid point ( i, j) and has a length of the grid spacing h, is drawn. The end of the normal is called /3, the intersection of the normal with the interface is denoted as a. The value of ".!,j is estimated by interpolation through values at a and /3 as : = aa + b (6 + 6 + ~4
PAGE 93

83 (a) Figure 2.27: Illustration of change of the phase on some grid points following the interface movement. (a) Location of the interface at the time level n. The shaded region designates the vapor phase. At this instant the grid point ( i, j) lies in the vapor phase. (b) Location of the interface at the time level n + l after moving. The grid point ( i, j) now lies in the liquid phase Those grid points which undergo the change of the phase are marked by open circles. Also shown in (b) is the schematic of finding the value at the ( i, j) by drawing a normal probe passing through ( i, j), the value at the end of the normal probe is obtained by using a bilinear interpolating function over points numbered 1 to 4. determine the correct interface location which satisfies all interfacial conditions (1.23), (1.24), (1.25) and (1.26), an iterative process for updating the interface location is needed as described in section 2.4.3 and 2.4.4. The interfacial temperature is specified using Eq. (1.26), the vapor pressure at the interface is p 00 then the liquid pressure at the interface can be determined by Eq. (1.24). The energy conservation condition Eq. (1.25) is used to estimate the interface velocity relative to vapor velocity at the interfa ce each time, thus it is naturally satisfied. The mass conservation condition Eq. (1.23) serves for determining the liquid phase velocity at the interface when the interface velocity and vapor phase velocity at the interface are known. Both liquid and vapor phase velocities at the

PAGE 94

84 interface are assigned as the boundary conditions respectively for solving the flow fields inside and outside the vapor bubble. 2.4.3 Iterative Process for Determining the Interface Shape In this section, we will describe the process for determining the interface shape using normal stress balance coupled with flow solver via an iterative process. This is the most important technique and difficult numerical procedure in this work. The method depicted in the following is used to determine the bubble deformation caused by the momentum balance with no phase change occurs. The technique implemented in the present research is the one similar to that proposed by Ryskin & Leal (1984) in their method using body-fitted moving grid Here we implemented it in the fixed grid which results in a more efficient overall solution procedure. If phase change is to be considered and contributes, in conjunction with the mo mentum balance, to the shape change of the bubble, additional algorithm is required and developed in the present study. They are given in the next section. Given an initial interface shape and flow field, material properties are assigned in the liquid and vapor phases, the flow field is then solved until a converged flow field is obtained when the interface is fixed. This flow field serves as the starting point for subsequent moving interface problems. Then the solution algorithm proceeds iteratively through the following steps starting from iteration number k=O for next time step n + l the initial guess for flow variables at time step n + l is 4p+i o = n the initial interface velocity is u~n~l,O = uint The interface shape is determined by the normal-stress balance Eq. (1.24) The technique proposed by Ryskin & Leal (1984) is adopted to obtain the estimate of the interface shape for next iteration. The idea behind this technique is very simple. This approach essentially takes the local imbalance between the total normal stress Tn, which includes both static and dynamic pressure and viscous contributions and

PAGE 95

capillary forces 1 II(s) = r. -r;, n We 85 (2.44) as a "driving force" which causes a local displacement of the interface in the normal direction, the magnitude of the local displacement is proportional to II ( s). Thus n+l,k+l xint n+l,k+l Yint (2.45) where the under-relaxation coefficient (3 has to be determined by numerical experi ment, its typical values being 0(104 to 10 3 ) in our computations In the cases where interface completely encloses one of the two fluids, e.g., bubble in our study the local increment in the location of interface marker points must be done in such a way as to satisfy the mass conservation constraint, i.e. preserve the bubble volume if density is constant. We know the change in the volume between the kth and (k + l)th iterations is 1 s[( n+l,k+l n+l,k)2 + ( n+l,k+l n+l,k)2]1/2d xint xint Yint Yint 8 0 (2.46) where the integral is taken along the interface. Since [( ~+l,k+l n+l,k)2 + ( ~+l,k+l ~+l,k)2]1/2 ~ Ilk( ) Xmt Xmt Ymt Ymt S (2.47) we have the mass constraint (2.48) This constraint determines a free constant of the pressure field in the total normal stress Tn at each iteration k so Eq. (2.48) is satisfied. Even after this constraint has been satisfied, the bubble may still change volume slightly at each iteration due to higher-order numerical errors. To prevent these small changes from accumulating and becoming significant, a simple scaling of the

PAGE 96

86 interface can be done following the Eq. (2.45). The uniform scaling magnitude L1 of the interface location in the normal direction can be determined from (2.49) where L1 V is the error in the volume. This simple scaling is very effective normally one or two iterations of this process is sufficient to bring the percentage error in the volume down to 105 The global iterative process within each time step is therefore as follows. 1. For a given shape of the bubble, the flow field is computed by solving the Navier-Stokes equations with a small number of iterations on pressure Poisson equation. 2. Knowing the flow field, the normal stress balance at the interface is calculated and checked. If it is not satisfied, the interface shape is modified according to Eq. (2.45) so as to reduce the imbalance between the total stress and the surface tension force. 3. After each interface update, the mass conservation is enforced by rescaling the interface using Eq. (2.49). 4. The interface normal velocity is calculated from the kinematic condition ufn1 1 = ( xfn1 1 xfut) / L1t, which is the boundary conditions for solving the flow field. 5. Return to step 1 and repeat until all equations and boundary conditions are satisfied t~ a predetermined level of accuracy. 2.4.4 Determining Interface Shape with Phase Change The interface movement when phase change is present can be decomposed into two components: local propagation velocity due to evaporation or condensation, and the global movement resulted from buoyancy when there is no phase change.

PAGE 97

87 The interface velocity component due to phase change effect in the normal di rection, which is also the velocity component relative to the vapor phase, is as Eq (1.25), i.e., (unt+l = Ja. [[)Tz (kv) (8Tv)] mt,p Pe 8n k1 8n (2.50) where the superscript 'n+ 1' on the left hand side denotes the n + 1 time step while the right hand side uses the field values at time step n. The subscript 'p' on the left hand side means this is the interface velocity component due to phase change effect. So the interface movement in the normal direction due to phase change is n+l n + ( n+l ) "t xint,p n = xint n uint,p n L.l (2.51) After we have this new location of the interface, we can use the process in the section 2.4.3 to further determine the shape satisfying the momentum balance condi tion. The final interface velocity at time step n+ 1 obtained from kinematic condition is the combination of the components owing to phase change and momentum balance. Hence, the vapor phase velocity at the interface in the normal direction is given by: ( ) n+l ( )n+l ( )n+l Un v = Un int Un int,p (2.52) which is the boundary condition for solving the flow field in the vapor bubble. The liquid phase velocity at the interface in the normal direction, according to Eq. (1.23), is which is the boundary condition for liquid phase. 2.4.5 Discontinuity of Material Properties (2.53) In the currently available numerical methods reported in the literature for sim ulating two-phase flow, most of the single field approaches encounter the numerical

PAGE 98

88 instability caused by large property jumps of the two phases the maximum d e nsity ratio for phase change problems which is reported so far in the prior publications to our best knowledge is 100. However, for liquid-vapor phase phenomena the d e nsity ratio often tops 1000. In a typical situation of water at atmospheric pressur e, the density ratio of liquid and vapor is roughly 1600 For some air bubble problems with no phase change we saw publications dealing with density ratio as large as 1000 but with a finite thickness interface treatment to prevent instability In the single field approaches for any cell around the interface different proper ties will appear in the same discretized equation on that cell which cause the stiffness in that equation so as to deteriorate the numerical stability although easy impl e men tation of the single field algorithm is quite attractive. In our approach because we employ the explicit interface tracking and cu t cell technique around the interface and our carefully designed computational proc e dure solves separatel y the two fields corresponding to two phases while only on the inter face, the two fields communicate with each other to adjust there own flow fi e ld to satisfy the interfacial conditions the fundamental difference from other approa c h is that the discretized equation for any cell including those around the interfac e will involve identical material properties corresponding to the phase that particular cell lies in. Therefor e, there is no stiffness in the discretized equation on an y cell the large density jump will not be expected to exert severe effect on the numerical stabil ity in our approach although the interface has zero thickness. This benefit is at the cost of complex algorithm design and coding involved in the development of pr e sent techniques. 2.5 Discretization on Axisymmetric Geometry In the simulations, the computational domain is axisymmetric where it is as sumed that 8 / 80 = 0 and the swirl velocity u 0 is zero.

PAGE 99

89 In a differential formulation, the 2-D dimensional conservation equations Eqs. (1.4)-(1.7) written in an axisymmetric cylindrical coordinate system can be further expressed as or equivalently as Vu p [ 0 ;; + V (uur)] p [ 8 ;z + V (uuz)] PC-p [8;; + V (uT)l 0 Op ( 2 Ur) --+ V Ur-fJr r 2 fJp 2 fJz + V Uz (2.54) 0 Op Ur -+ V Vur fJr r 2 fJp f) z + V Vuz kV-VT ( 2.55) If the coordinates z and r of the cylindrical coordinate system are replaced by x and y, the analogy between Eq. (2.54) and the equations in Cartesian coordinates becomes apparent. If r is also set to 1, Eq. (2.54) becomes identical to those in Cartesian coordinates except with an additional term -ur/r 2 with U z = U x and Ur = Uy for differential forms. In terms of the finite volume discretization procedure, the conservation equations written in integral form can be obtained from Eq. (2 55) by integrating the both sides over an axisymmetric control volume. The integral equations are the same as those for planar 2D cases in Cartesian coordinate except for an additional term -ur / r 2 .6. V in the Ur momentum equation, where .6. V is the volume of the control volume. Any

PAGE 100

90 method described in previous chapters are equally applicable for discretizing the 2D axisymmetric equations in cylindrical coordinates. However, some differences need to be addressed in discretizing the equations in an axisymmetric control volume. In addition to the surface integrals of pressure force over cell faces Z i = canst zi-l = canst, ri = canst, ri-l = canst, as was the case in planar 2D problems, the radial component of pressure forces onto the two faces 0 = canst has to be considered also. The other difference compared to planar 2D problems is the calculation of cell face areas and cell volume. Since the control volume size in the 0-direction is one radian, the areas of cell faces are computed in the same way as in the planar 2D case with the length in 0-direction taken into account i.e., inclusion of a factor r f c ( where r Jc denotes the cell face center). The areas of the front and back faces are calculated in the same way as the volume in the planar geometry. The volume of axisymmetric CVs with any number of faces is (see Ferziger & Perie, 1996) : l N v 6. V = 6 L (zi-1 zi) (r;_ 1 + r; + riri-1) i=l ( 2.56) where Nv denotes the number of vertices, counted counterclockwise with i = 0 corresponding to i = Nv. With the similarity of the governing equations for the planar 2D problems in Cartesian coordinate and axisymmetric 2D problems in cylindrical coordinate the same computer code can be used for both planar and axisymmetric 2D flows; for ax isymmetric problems, one sets r = y to account for the length in the third dimension of the axisymmetric control volume, for planar problems, one sets r = l because the length of the third dimension is unity for planar 2D control volume.

PAGE 101

CHAPTER 3 VALIDATION OF CARTESIAN GRID METHOD In Chapter 2, we detailed the numerical method for dealing with fixed, moving, deforming geometry in the underlying Cartesian grid system. Also the phase change effect is incorporated into the moving interface updating procedure. The moving interface algorithm is based upon the fundamental Cartesian grid method for fixed complex geometries using cut cell discretization. Therefore, validating our numerical method is a two-step process, first, the spatial discretization scheme using Cartesian grid method for fixed geometries is to be checked. Secondly, on top of the fixed geometry Cartesian grid method, the algorithm for moving and updating dynamic interfaces has to be checked. It is the task of this chapter to fulfill the first step validation of the overall numerical method, i.e., to examine the accuracy and fidelity of the fixed interface Cartesian grid method described in 2 3. In the next Chapter, we will make some checks on the overall moving interface algorithm involving techniques described in 2.4, which completes our second step validation on the numerical techniques. 3.1 Validation of Global Accuracy of Cut Cell Approach As we have shown in section 2.3, a second-order accurate spatial discretization scheme adopting cut cell approach is employed in the vicinity of the interface. To confirm that we obtain the nominal second-order accuracy when solving the flow field over complex geometries, the most straightforward way is to compute a flow which has a curved immersed interface and one for which an analytical solution exists. 91

PAGE 102

92 The flow chosen here corresponds to the two-dimensional Stokes flow past a circular cylinder placed next to a moving wall (see Figure 3.1). The exact solution to this flow was given by Wannier (1950) and is reproduced here. Here we have simulated this flow using our solver on four different uniform meshes. The meshes have equal spacing in the x and y directions and have N x and Ny points in these two directions, respectively. The global error in the numerical solution is computed as N,Ny = N lN L I /;um e rical ef>';xact I x y j=l (3.1) In order to simulate Stokes flow, the convection terms have been turned off in our simulation Computations have been carried out in the domain shown in Figure 3 1 with the exact solution imposed on the boundaries. In Figure 3 2 we show a log-log plot of the error for both velocity components u and v versus Nx. Also shown is a line with a slope of -2 which corresponds to second-order accurate convergence Moving wall u .. Figure 3 1: Computational domain and computed streamline pattern of Wannier flow

PAGE 103

93 10-2 u V 103 w 1 o4 1 oS ..,__....._....._....._....,......, ______ 10 1 10 2 10 3 N Figure 3.2: Global error of u and v as a function of mesh points for Wannier flow. The plot clearly shows that the global error in our computed solution decreases in a manner consistent with a second-order accurate scheme This test therefore proves that the current approach of treating the fluxes in the boundary cells does indeed result in a solver which is globally second-order accurate. 3.2 Fidelity of The Numerical Methodology The exact solution of the Wannier flow allows us to confirm the accuracy of the solver in the Stokes flow regime Here we validate the solver in the finite Reynolds number regime by simulating steady and unsteady flow past a circular cylinder im mersed in an unbounded uniform flow over a wide range of Reynolds numbers where the Reynolds number is defined as Red = U 00 d/v with d being the cylinder diam eter and U 00 the freestream velocity. This flow had been studied quite extensively in the past and a number of numerical and experimental datasets exist for this flow which are useful for the purpose of validation Simulations have been performed at Red= 20, 40 80, & 300 and results were compared with established experimental and

PAGE 104

11 111111111111111111111111111 11111111111111111111 111111111111111111111111111 111111111111111111111111111111111 11111111111111 1111111111111111111111111 1111 1 111111111111111111111111111111 11111111111111111111111 11111111111111111 1111 1111 111111111111111r .. ... 111111111111111111 11111 1 111111 1111111111111111 ll111111111111111111111111 1 111111111111111111111111 11111111111111111111111111 1 1111 111111111111111111111111 111111111111111111111111111 1111111111111111111111111 1111111111111 1 1111111111111111111111 11111111111111 ,,,111 111111111111111111111 11111111111111 ) Tj0.8431 0.8471 0.8627 rg4.1602 0 0 10.3 330.19 646.31 Tm(111111111111 1 ) Tj0.8431 0.8471 0.8627 rg12.2686 0 0 12.9 337.35 641 Tm<3131319595B7B7B7B7B7B7B7B7B7B7B7B7B7B7B72D2D> Tj/Suspect <> BDC0.6275 0.6235 0.6353 rg0 Tc24.5 0 0 24.5 302.58 590.74 Tm<95> Tj0.8431 0.8471 0.8627 rg24.5 0 0 24.5 304.27 590.74 Tm<9520> TjEMC/Fh0 1 Tf5.4001 0 0 5.4001 209.28 558.85 Tm<7F> Tj/Fh1 1 Tf-0.035 Tc11.6674 0 0 12.1 212.91 558.85 Tm<2D2D2DB7B7B7B7B7B7B7B7B7B7B7B7B7B7B7B7B7B7B7B7B7B7B7B7B7B795> Tj0.7059 0.7059 0.7176 rg3.165 0 0 12.1 301.5 558.85 Tm(,,,111111111111 11111111 111 1111 111111111111111111111111 1111111111111111111 111111111111111111111 ) Tj0.7059 0.7059 0.7176 rg-0.035 Tc2.4758 0 0 4.6 304.38 542.42 Tm(1111111111 1 111 1111111111111111111111 1 11111111111111111111111 1111111111111111111 11111111111111111111 1 1111111111111111111111111111 1111111111111111111 1111111111111111111111 11 1 111111111111111111111111111 1111111111111111111 1111111111111111111111111 1 111111111111111111111111111 11111111111111111111 111111111111111111111111 1 111111111111111111111111111 111111111111111111 1111111111111111111111111 IIIIHlllllllllllllllllllllllll lllllllllllllllllll llllllllllllllllllllll IIII 11 11111111111111111111111111111 111111111111111111 111111111111111111111111111 1 111111111111111111111111111111 111111111111 11111111111111111111111111 111 -10 0 10 20 x/d 94 Figure 3.3: Typical non-uniform mesh for simulation of flow past a circular cylinder. numerical results. All these simulations have been performed in a large 30d x 30d domain so as to minimize the effect of the outer boundary on the development of the wake. Figure 3.3 shows the 152 x 156 non-uniform mesh used in the low Reynolds number simulations. At the inlet and along the top and bottom boundaries we spec ify velocities corresponding to the potential flow past a cylinder and a homogeneous Neumann boundary condition is applied at the exit boundary. We have also tested larger domain sizes in order to ensure that the results presented here are independent of the domain size. For all these simulations we first impose a small asymmetric dis turbance at the inflow boundary for a short period of time and then allow the flow to evolve naturally after this. For Red = 20 and 40 the wake eventually attains a steady symmetric state and this is consistent with the well established result that the cylinder wake is stable to perturbations below Red= 46 1 (see Jackson, 1987; Provensal et al., 1987; Williamson, 1996). Once the flow has reached a steady state

PAGE 105

95 8.5 "C >8 7.5 4 5 6 7 8 x/d 8.5 "C 8 >7.5 4 5 6 7 8 x/d 4., Figure 3.4 : Streamline p l ot of flow past a circular cylinder. (a) Red= 20 (b) Red= 40 we compute the drag coefficient defined by Cn = (~%1~';!:d and the length of the recirculation zon e and then compare these with established results. The streamline p l ots in Figure 3.4a and b show the mean recirculation regions beh i nd the circular cylinder at Red = 20 and 40 respect i ve l y In Table 3 1 results in this steady flow reg i me obtained by the current method are compared with nu merical simulation by D ennis & Chang (1970) as well as ex p erimental measurements of Tritton (1959). It is fo u nd that our resu l ts compare well with other numerical sim u lations and experiments. It i s generally accepted that the wake of a cylinder immersed in a freestream first becomes unstable to perturbations at a critical Reyno l ds number of about R e d = 46 1 (see Jackson 1987; Provensal et al 1987). Above this Reynolds number a small asymmetric perturbation in the near wake will grow in time and lead to an unsteady wake and Karman vortex s h edding. This is indeed what we find for

PAGE 106

96 Table 3.1: Comparison of mean drag coefficient, length of wake bubble Lw (measured from rear end of the cylinder) and Strouhal number with established results. ReynoldsN um b e r --+ 20 4 0 80 30 0 Mes h --+ 152 X 156 152 X 156 217 X 183 217 X 183 St u d y .J.. CD Lw/d CD Lw/d CD St CD St Tri tt o n (1 9 5 9 ) 2.22 -1.48 -1.29 ---We i se ls be rger (1 922 ) 2.05 -1.70 -1.45 -1.22 -D ennis & C hang (1 9 7 0 ) 2.05 0.94 1.52 2.35 ----Fornberg (1 9 8 0 ) 2.00 0.91 1.50 2.24 ----Will i amso n ( 1996 ) -----0.15 -0.20 Cu r r ent 2.03 0.92 1.52 2.27 1.37 0.15 1.38 0.21 our simulation at Red = 80 Figure 3.5 shows the variation of the lift and drag coefficients with time and it shows how vortex shedding develops to a periodic state in time. The computed mean drag coefficient from the current simulation is about 1.37 which lies between the two experiments results. The Strouhal number for vortex shedding is defined as St = f d/U 00 where f is the shedding frequency and is one of the key quantities that characterizes the vortex shedding process Here we have estimated the Strouhal number from the periodic variation of the lift coefficient and the value comes out to be 0.15 which compares very well with the value obtained from experiments (see Williamson, 1996). Figure 3.6 shows a plot of the streamline pattern and spanwise vorticity contour at one instant and both plots show the classical Karman vortex street. In addition to the low Reynolds number simulations, we have also carried out a simulation at a moderately high Reynolds number of 300 This simulation serves to demonstrate that the current methodology is capable of resolving thin boundary layers that develop in flows at these Reynolds numbers. The mesh used for this simulations is the same as that used for Red = 80. This is a relatively coarse mesh and this coarse resolution severely tests the discretization scheme used in our solver for the boundary cells. The mean drag and Strouhal numbers have been computed from time

PAGE 107

2.5 2 1.5 Cd 0 1 ,:, 0 0.5 0 c, -0.50 50 100 150 200 tU J d Fig ur e 3.5: Variation of lift a n d drag coeffic i ents w i th t i me for Red = 80 8.5 8 "C >, 7 5 7 9 ,, 8 >7 6 4 5 6 6 8 7 x/d 10 x/d 8 12 9 10 14 16 97 Fig ur e 3 6: Instantaneous streamline p l ot and vorticity contour plot in the near wake of the circu l ar cylinder for Red = 80 ( a) Streamline p l ot (b) Vorticity contours

PAGE 108

98 variation oflift and drag coefficients and they are included in Table 3.1. It can be seen that the Strauhal number matches well with the experiments of Williamson (1996). It should be pointed out that at this Reynolds number the cylinder wake is intrinsically three-dimensional whereas our simulation is two-dimensional and therefore does not allow spanwise variations. It is well known that the 2-D simulation of a flow which is really 3-D over-predicts the drag. This is indeed what we observe for the current simulations The drag coefficient from our 2-D simulations is about 12% higher than the experimentally determined value of Wieselsberger (1922). In Figure 3.7 we have shown contour plots of spanwise vorticity at one instant. Figure 3.7a shows a view of of the wake that extends to about 10 diameters down stream from the cylinder and as expected, this plot shows the formation and evolution of compact Karman vortices in the wake. Figure 3. 7b is a closeup view of the flow around the cylinder and the mesh superposed on the greyscale contour plot clearly shows that there are fewer than five points in the attached boundary layer. It can be seen that even with a relatively low resolution provided here, the boundary lay ers on the cylinder surface are smooth indicating that the current treatment of the boundary cells adequately resolves thin boundary layers. 3 3 Su mmary In this chapter, the numerical algorithm employing Cartesian grid method to simulate flow with fixe d complex geometries are validated through extensive test cases The results have confirmed the accuracy and fidelity of the current compu tational method. The first step in our validation process is completed. In the next chapter, we will conduct the second step of the validation process, i.e. validate nu merical techniques for m o ving interfaces.

PAGE 109

9 8 7 C 6---.....---------------~------4 8.5 8 ,, 7.5 7 6 8 10 12 x/d 4.5 5 5.5 6 6.5 7 7.5 x/d 1 4 99 Figure 3. 7: Spanwise vorticity contour plots in the wake of the circular cylinder for Red= 300 (a) View extending to 9d downstream of the cylinder. (b) Closeup view showing the resolution provided to the attached boundary layers and separated shear layers

PAGE 110

CHAPTER 4 BUOYANCY-DRIVEN BUBBLE MOTION The three major components of the present numerical methodology ar e int e rfac e representation s c heme, Cartesian grid flow solver and moving interface algorithms In previous chapters the capability of two of the three components are examined. The accuracy of the interface tracking algorithms is tested and demonstrated The accuracy and fidelity of the Cartesian grid method for solving fluid flows over fixed sharp interfaces with complex geometries are also validated through extensive c ases. Since the main objective of this research is to study vapor bubble heat transfer by direct numerical simulations we present in this chapter the development of pre vious numerical methods on an axisymmetric cylindrical coordinates system All the vapor bubble simulations are performed on this system. Specifically the solution pro cedure dealing with the axisymmetric geometry, the B-spline curve fitting and fairing algorithms for interface representation and curvature evaluation, the interfac e up dating procedure based on normal stress balance and the global iterative process for solving flow fields and interface motion simultaneously are validated. The problems simulated in this c hapter involve no phase change. The flow investigated is unsteady axisymmetric and laminar. For this type of problems, there have been man y reported investigations. The results abound and physical aspects are well understood so that we may validate the entire numerical method by comparing our solutions with o t hers On the other hand, the buoyancy-driven motion of a bubble is a typical moving boundary problem in fluid mechanics The physical process is highly complex subject to multiple effects By simulating such flows with the dynamic interface, w e can obtain a valuable understanding and insight of the solution characteristics both th e 100

PAGE 111

101 flow field and bubble deformation of our method. These simulations thus serve as the essential precursor to numerical simulations of rising bubble with phase change that will be detailed in the following chapters. 4.1 Summary of Previous Study Although a single bubble rising in a liquid has been studied for a number of years most theoretical results are limited to very small deformation at either low or high Reynolds numbers For example, at very low Reynolds numbers there e xists the theoretical model by Taylor & Acrivos (1964) based on an asymptotic theory. At high Reynolds numbers, only boundary-layer approximations (see Moore 1963 ; Harper & Moore 1968) and semi-empirical models (see Parlange, 1969) are available. All of the above investigations assume that the bubble maintains a spherical shape which is rather unrealistic at high Reynolds or Weber numbers. Experimental investigation of Saffman (1956) and Bhaga & Weber (1981 ) pro vided a fairly detailed picture of the gas bubble upward motions in a liquid. Ryskin & Leal (1984) have reported the first successful theoretical solution for motion of bubbles with finite degree of deformation using body-fitted moving grid techniques. In their model the interface is treated as zero thickness ; the shape is explicitly determined by the stress balance at the interface However, the problems they considered involve one viscous fluid surrounding the bubble and a void bubble with Pv = 0, v = 0 There is no flow field inside the bubble. So the interfacial conditions involve forces on only one side of the bubble interface. Later, Dandy & Leal (1989) have extended the numerical method of Ryskin & Leal (1984) to consider the deformable drop problems involving two viscous fluids both inside and outside of the drop In both Ryskin & Leal (1984) and Dandy & Leal (1989) only steady state problems are considered. Furthermore, orthogonal body-fitted coordinates are adopted to generate the grid system.

PAGE 112

102 In recent years, many numerical simulations of the bubble motion have been reported in the literature. However in virtually all cases, the interface is not sharply defined, and the stress balance is enforced via several cell spacings, instead of at precisely defined locations. Examples include the immersed boundary method ,and volume of fluids method as listed in Chapter 1. The nature of those methods does not render possibly the accurate prediction of the exact location of deformed interface. While our method solves the two separate fields with stress conditions matched at the interface to obtain the exact, instantaneous interface shape. The comparison between our results and those obtained recently (e.g., see Chen, Garimella, Reizes & Leonardi, 1999) is therefore not meaningful owing to the difference in dealing with interface stress conditions. To help assess the performance of the present method, we will make direct com parison with the results reported by Ryskin & Leal (1984) for the bubble motion at different Reynolds and Weber numbers. It is noted that the results of Ryskin & Leal ( 1984) agree well with experimental study of Saffman (1956) and Bhaga & Weber ( 1981). 4.2 Grid Resolution Study and Convergence Criteria Referring to Figure 4.1, unless otherwise mentioned, all the cases reported in the above are computed on a 10 x 3 domain with a 250 x 75 uniform mesh. The boundary conditions for the cases presented below are: at three far sides of the boundary the outflow (zero velocity gradient) condition is specified for velocity. At the line of symmetry, the mirror condition for all variables is used. The effect of grid resolution on solution accuracy is examined first. Unless otherwise mentioned, all simulations in this study employ the same grid resolution around the interface that is, the number of cells across the initial bubble diameter is 25. In order to assess the grid dependency of the solution, we have conducted computations for one case: Re = 100, We= 4, Pv/ p 1 = 0.001, v/ 1 = l.0, using

PAGE 113

103 L H rl_ z Figure 4.1: Schematic of the computational domain and boundary conditions For cases involving stationary bubble, L/ R.o = 48, H / R.o = 24; For rising bubble cases without phase change, L/ R.o = 20, H/ R.o = 6; For rising bubble cases with phase change, L/ R.o = 20, H/ R.o = 8 500 x 150 grid with 50 cells across the initial bubble diameter. Figure 4.2 shows the time history of the aspect ratio of the bubble. The aspect ratio is defined as the length along the major axis divided by that along the minor axis. As can be seen in Figure 4.2, the difference between the results from the two grids is small. The criteria for determining the convergence are: the maximum norm of the absolute error must be less than 106 for momentum equations (2.12); 104 for pres sure Poisson equation (2 .15 ); 106 for energy equation in phase change cases; 103 for evaluating the normal stress balance (2.45); When these convergence criteria are satisfied, the maximum norm of the residual of mass conservation for the entire fl.ow field is less than 103 4.3 Results For Buoyancy Driven Bubble Motion Representative cases of Ryskin & Leal (1984) are simulated in the present study. We compute the solutions for Re in the range of 1 :S Re :S 100 and We from 0 up to 20 for Re :S 20 and up to 10 for Re 50. The parameter ranges are in line

PAGE 114

2.2 0 2 :;:: co a: 1.8 u G) C. 1.6 U) <( G) 1.4 :c .c :::, 1.2 m 1 0 1 2 3 Time 250 X 75 500 X 150 4 104 5 Figure 4.2: Comparison of the aspect ratio history for the rising bubble for Re = 100, We= 4, Pv/ Pl = 0.001, v/ = l.O on the 250 x 75 grid and the 500 x 150 grid with those of Ryskin & Leal (1984). All computations are done in time dependent manner. To facilitate the direct comparison, the condition of balance between the drag force and buoyancy force used in the work of Ryskin & Leal (1984) is employed in the present study to determine the Froude number for the hydrostatic pressure term of (1.24). The relation used by Ryskin & Leal (1984) is 2Rg ~C U2 4 D (4.1) where R is the bubble radius, U is the terminal velocity of the bubble, Cn is the drag coefficient. The left hand side of ( 4.1) is actually 1 /Fr. For each case, we use (4.1) to find the Froude number from a given Cn value. This procedure ensures that the scaling processes between the current and that used in Ryskin & Leal (1984) are consistent. Of course, the drag coefficients are computed from the solution obtained. The drag coefficients obtained in the present study given in Table 4.1, are those when the bubble reaches constant rising velocity.

PAGE 115

105 The drag coefficients obtained by our simulations and by Ryskin & Leal (1984) are summarized in Table 4.1. Also included in the table is the error estimate reported in Ryskin & Leal (1984) based on energy dissipation analysis of their numerical simulation. This information helps one gain a sense of the accuracy level in that work. Table 4.1: Comparison of drag coefficients from present simulations (first row) with those obtained by Ryskin & Leal (second row) using integration of the forces at the surface. The third row shows the relative deviation of drag coefficients computed via energy dissipation in Ryskin & Leal's computations. In Ryskin & Leal's computations, the bubble is consider to be a void, while in the present case, Pi/ Pv = 1605, if v = 22. I Re.J,. \We I 2 3 4 6 8 10 12 15 16 20 2 10.9 11.4 11.6 11.8 11.9 10.6 11.0 11.1 11.2 11.3 2% 2% 3% 4% 6% 5 5.0 5.5 6.1 5.00 5.48 5.99 0.2% 0.5% 4% 10 2.9 3.4 4.1 4.4 2.92 3.41 4.00 4.25 0% 0.5% 3% 8% 20 1.7 2.2 2.6 3.0 3.3 3.7 3.7 1.74 2.16 2.56 2.94 3.22 3.55 3.60 0.5% 1% 2% 5% 10% 5% 4% 50 0.9 1.2 2.3 0.88 1.23 2.18 0.5% 0% 12% 100 0.5 0.8 0.54 0.81 0.5% 1% The favorable overall agreement in drag coefficients between the two simulations shows that our method is capable of correctly predicting the dynamic behavior of a coupled system involving the liquid flow field and vapor bubble.

PAGE 116

106 The computed bubble shapes for selected cases in Table 4.1 are shown in Figure 4.3. The shapes are the ones when the unsteady bubble motion reaches the terminal velocity.

PAGE 117

2 3 4 6 8 10 12 15 16 20 Re 2 0 0 C) C) 5 0 0 10 0 0 C") C") 20 0 0 C) 50 0 0 C=-=> 100 0 0 Figure 4 3: Computed terminal axisymmetric shapes of rising bubbles as a function of R e and W e with ptf P v = 1605 tf v = 22

PAGE 118

108 Overall, the trend of bubble shapes changing with increased We is in agreement with common experimental observation: spherical to oblate-ellipsoidal and then to oblate-ellipsoidal/spherical cap Bhaga & Weber (1981). Figure 4.4 shows development of the bubble shapes for: (a) Re = 5, We = 10, pz/ Pv = 1605, z/ v = 22, (b) Re = 2, We = 16, pz/ Pv = 1605, z/ v = 22. The corresponding streamlines for the two cases, plotted based on the coordinate fixed at the middle of the lower surface of the moving bubble, are shown in Figure 4.5. Two recirculating structures are observed in each case, one inside the bubble, and the other caused by the interaction between the bubble and the liquid. It is interesting to observe with Re= 2 and We= 16, the recirculating flow in the liquid phase is, as expected, attached to the bubble; for Re = 5 and We = 10, it tends to detach from the bubble. Figure 4.6 compares the steady bubble shapes at three density ratios. The differences observed are small. The reason for this phenomenon is that by fixing Re and We, the only impact from the density ratio is via the unsteady and convection terms in the momentum equation in the vapor domain. For a rising bubble, since the fluid dynamics inside the bubble is induced by the interface movement, for the present Re and We ( defined based on the properties of the liquid phase), the impact from the vapor portion of fluid dynamics is limited. Accordingly, only minor differences are observed in Figure 4.6. Nevertheless, it is reasonable to observe that as the density ratio increases, the bubble becomes less deformed. Figure 4. 7 shows the bubble shapes at selected time instants for the three density ratios. These observations have been reported previously in other studies, including Dandy & Leal (1989). The drag coefficients for these three cases from our simulations are 1.29, 1.32, 1.34 for density ratios of 0.1, 0.01, 0.001, respectively. The drag coefficient reported in Dandy & Leal (1989) for varying density ratios under the corresponding Reynolds and Weber numbers is 1.29.

PAGE 119

109 Re=2, We=16 Re=5, We=10 w w 0 0 0 (a) (b) F i gure 4.4: D evelopment of bubb l e s h a p es: (a) Re = 2, We = 16 pi/ Pv 1605, if v = 22 and (b) Re= 5, We= 10, Pi/Pv = 1605, if v = 22.

PAGE 120

110 Re=2, We:16 .. u Re=S, We:10 .. u F i gure 4.5: Flow structures at the term i nal state for cases corresponding to Figure 4.4, (a) Re= 2, We= 16, PdPv = 1 605, if v = 22, and (b) Re= 5, We= 10, pif Pv = 1605, if v = 22. The stream l ines are observed on the reference frame attached to the moving b u bble.

PAGE 121

-::--.;;"7 // i ft I I "',,, ,,,. ,,, ,,, .::.:,,........ _,;. -Figure 4.6: The terminal state shapes for cases with Re 1, Pv/ Pl = 0.1, 0.01, 0.001, v/ l = 1 Pv / p 1 = 0.001 Pv / P1 = 0.01 Pv I P1 = 0.1 111 100, We = 4, Fr To further illustrate the effect of density ratio on the computational performance, Figure 4.8 compares the convergence histories between two cases with different density ratios. The residues of both the Young-Laplace equation, (1.24), and the pressure Poisson equation, at a given time step, are shown. The residues are based on the sum of the absolute value of the residue computed in each cell. The levels shown in Figure 4.8 are not normalized. The figure demonstrates that the present method is robust in terms of handling the large property variations across the phase boundary. Figure 4.9 shows the flow structures corresponding to three density ratios, each with Re = 100, We = 4, v/ = 1.0. This figure corresponds to the same pa rameters as those shown in Figure 4.6 and 4. 7. In all cases, the recirculating wake is detached from the bubble. Again, there is no significant difference for different density ratios.

PAGE 122

C) C) C) 8 C) C) C) 8 C) C) C) 8 112 Figure 4.7: The shape evolution for cases with Re= 100 We= 4 Fr= 1 v/ = 1 at equal time intervals.

PAGE 123

10 4 10 3 10 2 10 1 cu ::::J 10 "C "' 10 1 Q) a: 10 2 10 3 10 4 10 5 0 10 Pressure, pjp 1 =0.001 Normal Stress Balance, pjp 1 =0.001 Pressure, pjp 1 = 0.1 Normal Stress Balance, pjp,=0.1 20 Iteration 30 113 40 Figure 4.8: Convergence paths of the Young-Laplace equation at the interface and the pressure equation in the entire domain within a time step. Here two different density ratio cases are shown with Re= 100, We= 4, Fr= 1, v/ 1 = 1.

PAGE 124

114 Re:100, We:4, Pv / P, = 0.1 ... u Re:100, We:4, Pv / p, = 0.01 u Re:100, We:4, Pv / p 1 = 0.001 ... u Figure 4 9: Flow structures at the terminal state for cases with Re = 100, We = 4, Fr= 1, v/ = 1 Pv/ Pl = 0.1, 0.0 1, 0.001; The stream li nes are observed on the reference frame attached to the mov i ng b u bble.

PAGE 125

115 4.4 Summary By applying the present method to simulating the buoyancy-driven motion of a bubble, the method is proved to be accurate and capable of resolving the interfacial dynamics and shapes. It should be noted that for the bubble motion c onsidered here to my best knowledge comparison of drag coefficients with experiments and established theoretical results were not reported in other fixed grid based methods because those methods usually are unable to resolve the exact sharp interface While the drag coefficient is strongly influenced by interfacial dynamics, correct calculation of drag coefficients requires the method to resolve the dynamic interface To this extent the present method is suitable for problems when interfacial dynamics plays an important role in the system.

PAGE 126

CHAPTER 5 BUBBLE RISE AND GROWTH DUE TO EVAPORATION In Chapter 4 the simulation results based on the current method for buoyancy driven bubble motion without phase change are compared with other s work. In this chapter, the results for the heat transfer-controlled bubble growth due to evaporation in a superheated liquid under either zero gravity or normal gravity conditions are presented. 5.1 Stationary Bubble Growth Under the zero gravity condition, the bubble motion resembles the ideal case of a stationary bubble growth studied in many early works such as Forster & Zu ber (1955), Scriven (1959) and Prosperetti & Plesset (1978). In those studies the bubble is assumed in a spherical shape, thus a 1-D problem with the bubble radius as the dependent variable is solved in conjunction with the thermal boundary layer approximation. The momentum effect in the liquid and vapor phases on the bubble growth and its shape is totally discarded In doing so, the bubble growth rate i.e ., the radial velocity depends only on the Jakob number. All those studies concluded that the time evolution of the bubble growth radius follows the t 1 1 2 law. For the heat transfer-controlled stationary bubble growth, the major transport mechanism is the heat diffusion, the appropriate scaling for the velocity will be based on the diffusion mechanism : Ur = aif L, where a 1 is the liquid thermal diffusivity and L is the initial bubble diameter. Using the diffusion scale for the velocity the Peclet number is always 1.0. 116

PAGE 127

117 5.1.1 Case with Thermal Properties of Water First we computed a stationary bubble growth based on the thermal properties of water under the atmospheric conditions, as listed in Table 1.3. A !::1T = 1 C superheat is considered. Under this superheat, the bubble departure diameter is about 0.3 mm according to empirical calculations (see Carey, 1992), which is our initial bubble diameter, and thereby the length scale L Using the thermal properties of the water listed in Table 1.3, the thermal diffu sivity of water at atmospheric pressure is a 1 = 1.679 x 107 m 2 / s. Hence the velocity scale is Ur= at/ L = (1.679 X 107 )/(0.3 X 103 ) = 5. 5967 X 104 m/ s From the velocity and length scales, all required dimensionless parameters for this case are estimated as R = P1UrL = {958.3)(5.5967 x 104 )(o.3 x 103 ) = 0 58 e 1 277 53 x 106 W = PlU~L = {958 3)(5.5967 x 104 ) 2 (0.3 x 103 ) = 1 53 10-6 e u 0.0589 X J El.. c,, 1 t:.T (958 3) ( 4220)(1) 3 O a P v >. (0 597) 2256700 The Prandtl number is 1.72. The Weber number at this magnitude essentially means no driving force for the deformation because everything involved is radially symmetric. The above dimensionless parameters are thereby used in the simulation. At time t = 0, the liquid phase is uniformly superheated, i.e., an initial temperature step exists: T(lxl >Ro)= T 00 T(lxl :S Ro)= Tv. A non-uniform mesh on a domain of 48Ro x 24Ro is used in the computations. Ro is the initial bubble radius. The resulted transient dimensionless bubble radius R( t) /Ro as a function of time is plotted in Figure 5.l(a). The stationary bubble growth exhibits an asymptotic slope of 1/2 in the log-log plot. This agrees with the well accepted theory that the growth radius is proportional to the square root of the time. The 1/2 law is valid when the thermal boundary layer around the bubble is well developed after the initial stage.

PAGE 128

(a) 10 2 -------.-------.---..,.....---, 10 1 .... 0 er: er: 10 1 ff 1 ..__.i....i..,1,,,,1,,1,,1,1.1.1.... ...................... __._ ................ ____. ............. ___ 1 ff 5 1 ff 4 1 ff 3 1 ff 2 1 ff 1 10 Time (b) 9 10 l,-----+-----+---1-"'-'--,,.,,C---+------l ,... I 0 er: cc 1 ff 1 ________ ._.._ __ ------i ___ ...,..._ __ 103 ........... -------------1 ff 5 1 ff 4 1 ff 3 1 ff 2 1 ff 1 10 Time 118 Figure 5.1: Growth rate of the bubble radius for a stationary bubble obtained in the present simulation for Re= 0.58, We= 1.53 x 106 Ja = 3.0, Pr= 1.72, Pe= 1.0 pi/ Pv = 1605, if v = 22. (a) Dimensionless bubble radius (b) Growth portion of the bubble radius.

PAGE 129

119 In Figure 5.l(b), we plotted only the growth portion (R/R 0 1.0) of the bubble radius versus time The plot also shows that the bubble growth rate approaches the diffusion controlled limit t 1 1 2 The above case with water required approximately 20, 000 time steps for the bub ble growth rate to approach the diffusion controlled limit because of the magnitude of the Jakob number Ja = 3.0, which requires very small time steps. 5.1.2 Case with Thermal Properties of Ethanol In the second case, we use the thermal properties of Ethanol as listed in Table 1.4. To be efficient, a smaller Jakob number Ja = 0.l is chosen so that a larger time step could be used. All the other dimensionless parameters are determined in a way similar to that of the water case described earlier. The Prandtl number is Pr = 8.37. All dimensionless parameters are listed as Re= 0.12 Pe= 1.0 We= 3.91 x 107 Pr= 8.37 Ja = 0.l Figure 5.2 shows the growth portion of the bubble radius versus time for this case. The growth radius is also following the relation R(t) ex: t 1 1 2 The two cases confirm the analytical prediction for the growth of the stationary bubble, namely, R(t) ex: t 1 1 2 5.2 Translating Bubble Growth Under the influence of gravity the bubble would be rising and growing simulta neously. Darby (1964) and Ruckenstein & Davis (1971) suggest that the growth rate

PAGE 130

120 10 1 ---------------------1 0 :-----+----+-----+-----+----t C? ,'o 10-1 1,-----+----.+----+,,,C...-~.J-----I a: a: 1 03 ............................................................................... _._.. ....... ...._ ................ .... 1 04 1 0 3 1 02 1 01 1 0 1 0 1 Time Figure 5.2: The dimensionless growth radius for a stationary bubble obtained in the present simulation versus time for Re= 0.12, We= 3 91 x 101 Ja = O.l, Pr= 8.37, Pe = l.0, pi/ Pv = 527, if v = 41. is significantly higher when the relative motion between the bubble and the surround ing liquid is present based on their analysis. A growth rate R(t) ex: t 2 1 3 is speculated under such conditions. For the rising and growing bubble due to buoyancy an appropriate scale for the velocity is Ur = ,/gl where g is gravity and L is the initial bubble diameter. Using t h is scale, the Froude number is always 1.0. Two cases were computed for the rising bubble case based on the thermal prop e r ties of water and ethanol, respectively The dimensionless parameters for the two cases are summarized in the following The simulations were conducted in the time dependent manner as in all other cases in this work. The simulations were stopped when the thermal boundary layer around the bubble became fully developed. The d i mensionless times when the two cases were stopped respectively were not the same.

PAGE 131

121 1. The first case we simulated is based on the thermal properties of water under atmospheri c conditions. However for convenience the Prandtl number is kept as 1.0. Re y nolds number is 10 Peclet number is 10 because of the value of Prandtl number Pr = 1.0 The Weber number is set as 0.2 which corresponds to a slight deformation case. The Jakob number is 1.0. The system dimensionless parameters are listed in the following. Re= 10 0 Fr= 1.0 Pe= 10.0 We= 0 2 Pr= 1.0 Ja = 1.0 2. The second case we simulated is based on the thermal properties of Ethanol under the atmospheric conditions as listed in Table 1.4. A !:1T = 10 G su perheat is considered. Under this superheat the bubble departure diameter is about 1.0 mm according to empirical calculations (see Carey 1992) which is our initial bubble diameter and thereby the length scale L. The velocity scale is thus Ur = .fgL = 9.9 x 10 2 m/ s. With these reference scales the dimensionless parameters for this case are estimated using the procedure given in Section 5 1.1 and listed as follows: Re= 174.81 Fr= 1.0 P e = 1463. 7 We= 0.42 Pr= 8.37

PAGE 132

122 Ja = 16.43 The bubble growth radius history for these two cases is shown in Figure 5.3. The dotted lines have slopes 1/2 and 2/3 respectively. The slopes of the bubble growth radius from the present simulation fall in between the two dotted lines. This figure shows that the growth rate is higher when bubble is rising as compared to a stationary bubble However, the growth rate is approaching but does not reach the t 2 1 3 limit, which is probably due to the deformation of the bubble from a sphere. At the end of the simulations, the streamlines for the two cases, plotted based on the coordinate fixed at the middle of the lower surface of the moving bubble are shown in Figure 5.4. The recirculation inside the bubble, which starts a short distance from the interface, is observed even for the bubble undergoing phase change In the vicinity of the interface, the streamlines emanating from the interface are a result of evaporation. The detached recirculating wake is also observed for the case with higher Reynolds number. The corresponding temperature profiles around the bubble for these two cases are shown in Figure 5 5. The thermal boundary layer is formed where the boundary layer is thinner around the upper surface, and thicker around the lower surface of t he rising bubble owing to the relative motion between the rising bubble and the surrounding liquid. It is also clear that the higher the Reynolds number the thinner t he thermal boundary layer. The tail-shaped structure in the wake results from the separation of the boundary layer. The local Nusselt numbers, i.e., dimensionless heat flux, along the bubble surface at the end of the simulations are depicted in Figure 5.6. The marching direction in the figure is from the middle of the lower surface to the middle of the upper surface of the bubble for the Nusselt numbers. The highest heat flux occurs at certain point where the relative velocity between the liquid and the bubble is the highest The

PAGE 133

(a) C? ,I 0 a: a: (b) 10 2 ae-----,------,------r-----, 10 1 ... 3 10 2 .. . . . ... 1 .... 1 ff 1 2 1 ff 2 .__..i,_1...1,,,1.....i,..._ ............ i..i....i. ......................................................... .... 103 102 10 1 10 1 o 1 Time 10 1 ..-----,----"'T"'""---""T""---~ 3 .... 2 .... .. .. 10 -------+---------r'7"/_~-------t .... :: : ... r 1 !? , ,, ::: :::: : : : .... 2 a: 1 ff 1 J----+":.....,---+-----+-----1 10 2 _..._. ................. ...._ ........................ ...._ ........................ ...._ ....................... 1 ff 2 1 ff 1 1 o 0 10 1 10 2 Time 123 Figure 5.3: Growth rate of the b u bb l e radius for a r i sing bubble obtained in the present sim u lation for cases (a) Re= 10.0, Fr= 1.0, We= 0.2, Ja = 1.0, Pr= 1.0, Pe = 10.0, pi/ Pv = 1605, z/ v = 22, {b) Re = 174 81, Fr = 1.0, We = 0.42, J a = 16.43, P r = 8 37, P e = 1463. 7, pi/ Pv = 527, z/ v = 41.

PAGE 134

12 4 (a) u (b) .. u Figure 5.4: Flow structures at the end of t h e s i m ul at i ons for cases (a) R e = 10.0 Fr = 1.0 W e = 0.2 J a = 1.0 Pr = l.0 Pe = 10 0 Pd P v = 1605 d v = 22 ( b ) R e = 17 4 .81 F r = 1.0 W e = 0.42 Ja = 16.43 Pr = 8.37 P e = 1 4 63.7 Pd P v = 527 d v = 41. The stream li nes are observed on the reference frame attached to the mov i ng bubble.

PAGE 135

125 spots of the maximum heat flux are not at the middle of the upper surface but somewhere downstream. The overall Nusselt number versus time for the two cases are plotted in Figure 5.7. The average heat flux almost became constant which means the thermal bound ary layers are fully developed at the time the simulations ceased. The case with the higher Reynolds number and Jakob number has higher average heat flux than the other. The time evolutions of the rising velocity of the bubble, which is the mean of the velocity values at the middle of the lower and upper surfaces of the moving bubble, are shown in Figure 5.8. The rising velocity of the bubble is increasing owing to its increased volume. It is noted that the present simulation computed the solution under the varying bubble rising velocity resulted from the buoyancy effect. Other i nvestigations on the bubble growth often employ an uniform mean flow to simulate t he effect of the relative velocity on the bubble growth. The present simulation is closer to the realistic situations. In Chapter 4 we studied the effect of density ratio on the bubble dynamics for Re = 100, andW e = 4. Here we compute the bubble growth under the same Reynolds and Weber numbers. Figure 5.9 shows the development of bubble shapes, under the influence of phase change and buoyancy, at selected time instants for three cases In these cases, the Re, We and viscosity ratio are fixed, while the density ratio is varied from 0.1 to 0.001. It is noted that the density ratio directly influences the Ja and the interface speed, as indicated in (1.25) Between the density ratios of 0.1 and 0.01, while the bubble size grows slightly faster with a smaller density ratio, the shapes are similar between the two cases. However, as the density ratio is reduced to 0 001, which is closer to a normal boiling heat transfer case, significant differences in bubble growth rate and

PAGE 136

126 (a) Level T 10 0 85 9 0.76 8 0 66 7 0 57 6 0.48 0 38 0 29 0 20 (b) Level T 10 0 85 9 0 76 8 0 66 7 0 57 6 0.48 5 0.38 4 0.29 3 0 20 2 0 10 0 01 Figure 5.5: The temperature profile at the end of the simulations for cases (a) Re= 10.0, Fr= 1.0, We= 0.2 Ja = 1.0, Pr= 1.0, P e= 10.0, pif Pv = 1605 if v = 22, ( b) Re = 174.81, Fr = 1.0, We = 0.42, Ja = 16.43, Pr = 8.37, Pe= 1463.7, pif Pv = 527, if v = 41.

PAGE 137

127 (a) 2.8 a.. 2 6 Q) .c 2.4 E 2.2 ::::, z 2 Q) 1. 8 en en ::::, 1. 6 z 1 .4 cu 0 1. 2 0 ..J 1 0 8 0 0 2 0 4 0 6 0 8.0 Location along the surface (b) 1 6 ... .8 1 2 E ::::, z = 8 ffl Cl) ::::, z 4 cu (,) 0 -I 0 0 0 2 0 4 0 6 0 8 0 Location along the surface Figure 5.6: The usselt number along the interface at the end of the simulations for cases (a) Re = 10.0, Fr = l.O, We = 0 2, Ja = l.O, Pr = l.O, Pe = 10 0, pi/ Pv = 1 6 0 5, if v = 22, (b) Re = 174.8 1 Fr = l.O, We = 0.42, Ja = 16.43, Pr= 8.37, Pe= 1463.7, pi/ Pv = 527, if v = 41.

PAGE 138

128 (a) 20 16 12 ::::J z 8 4 0 0 1 2 3 Time (b) 50 40 30 ::::J z 20 10 o----------------0 1 2 3 4 5 Time F i gure 5.7: The overall N u sse l t number versus time for cases (a) Re= 10.0, Fr= l.0, We = 0 2, J a = l.0, Pr = l.0, Pe = 10 .0, pt/ Pv = 1605, if v = 22, ( b ) Re = 1 74.81, Fr = l.0, We = 0.42, Ja = 16.43, Pr = 8.37, P e = 1463.7, pt/ Pv = 527, iJv = 41.

PAGE 139

129 (a) 1.2 c:; 0 8 .2 Q) > C) C: .!!! 0.4 a: o.__ ____ ....._ _________ 0 1 2 3 Time (b) 1. 2 c:; 0 8 .2 Q) > C) C: :~ 0.4 a: o.__ __ ....._ ___________ 0 1 2 3 4 5 Time Figure 5.8: The bubble rising velocity versus t i me for cases (a) Re= 10.0, Fr= l.0, We = 0 2, J a = l.0, Pr = l.0, Pe = 10.0, ptf Pv = 1 605, tf v = 22, ( b ) Re = 174 81, Fr= l.0, We= 0.42, Ja = 16.43, Pr= 8.37, Pe= 1463 7, PdPv = 527, if v = 41.

PAGE 140

130 pjp 1 =0.001 pjp 1 =0.01 P)P1 =0.1 C) C) C) C) C) C) 8 8 Figure 5.9: Development of phase change shapes for Re = 100, We = 4, Fr = l, v/ 1 = l with different density ratios. interface shape are observed. Figure 5 10 shows the flow structures at the final stage of each case.

PAGE 141

131 Re=100, We=4, pjp 1 = 0.001 u Re=100, We=4, pjp 1 = 0.01 u Re=100, We=4, pjp 1 = 0.1 u Figure 5.10: Flow structures at the end of the simulaions for cases with Re = 100, We = 4 Fr = 1 v/ = 1, Pv/ Pt = 0.1, 0.01, 0.001; The streamlines are observed on the reference frame attached to the moving bubble.

PAGE 142

CHAPTER 6 BUBBLE RISE AND COLLAPSE DUE TO CONDENSATION The collapse process is more challenging numerically than the growth case be cause of the fact that the ratio of bubble surface area and its volume is increasing with time in contrast to decreasing in the growth case. This trend in geometrical relation indicates that the radius variation induced by interfacial heat flux would be increasing as the bubble size is reduced. In the growth situation, the ratio between the thermal boundary layer thickness and bubble radius will eventually reach a constant value because they both grow in the same direction. The thin boundary layer theory could thus be used to find the heat flux on the bubble surface. However that is not the case for bubble collapse. Since the thermal boundary layer surrounding a collapsing bubble becomes thicker as the bubble size is reduced, the thin boundary layer approximation is not valid for describing the heat transfer in this expanding thermal region. Unlike bubble growth, there is no analytical solution valid for the entire collaps ing process. Some used the thin boundary layer assumption, such as Florschuetz & Chao (1965), to obtain an analytical result which is valid only at the beginning of the collapse for large Jakob numbers. Okhotsimskii (1988) solved numerically the heat transfer-controlled collapse based on the same theoretical model of Plesset & Zwick (1954) and tabulated the time evolution of the bubble radius over a wide range of Jakob numbers ( 102 < Ja < 10 3 ). In his analysis, he isolated Jakob number to be the only dominant factor of the collapse process. 132

PAGE 143

133 For a translating and collapsing bubble, Tokada et al. (1968) Akiyama (1973) Moalem & Sideman (1973), and Dimic (1977) etc. obtained equations for predicting the temporal decrease of the condensing bubble diameter. They are reviewed in Chen & Mayinger (1992) along with some other formulas. Chen & Mayinger (1992) gave an empirical correlation describing the temporal decrease in the bubble diameter using their experimental data as (6.1) in which 'Y = R/ Ro, Fourier number Fo = a 1 t/(2Ro) 2 and the Reynolds number based on the initial bubble diameter. For the rising and collapsing bubble, the same velocity scale Ur = ..fgL in the translating bubble growth cases is employed. Thereby the Froude number is always 1.0. Two cases were computed for the rising and collapsing bubble based on thermal properties of water and ethanol, respectively. The dimensionless parameters for the two cases are summarized in the following. 1. The first case we simulated is based on the thermal properties of Ethanol un der the atmospheric conditions, as those listed in Table 1.4. A D..T = l0C superheat is considered. Under this superheat, the bubble departure diameter is about 1.0 mm according to empirical calculations (see Carey, 1992), which is the length scale L. The velocity scale is thus Ur = ..fgL = 9.9 x 102 m/ s. With these reference scales and material properties, the resulted dimensionless parameters are Re= 174.81 Fr= 1.0 Pe= 1463.7

PAGE 144

We= 0.42 Pr= 8.37 Ja = 16.43 134 2. The second case we simulated is based on the thermal properties of water under the atmospheric conditions, as those listed in Table 1.3. A !:::,.T = 3C superheat is considered. Under this superheat, the bubble departure diameter is about 0.9 mm according to empirical calculations (see Carey, 1992), which is used as the length scale L. The velocity scale is thus Ur = -JgL = 9.39 x 102 m/ s. With these reference scales, the resulted dimensionless parameters are Re= 291.86 Fr= 1.0 Pe= 503.41 We= 0.13 Pr= 1.72 Ja = 9.0 In Figure 6.1, the temporal variation of bubble radius for case 1 obtained by the present simulation, correlation of (6.1) and prediction of Akiyama (1973) are all plotted. In our simulation, we compute the flow fields both inside and outside the bubble. A certain number of grid cells have to be enclosed inside the bubble. Therefore we were unable to compute the entire collapse process but to stop at a certain size to ensure that we have enough cells to resolve the flow fields inside the bubble. Thereby the dimensionless bubble radius didn't reach zero at the end of curves. Figure 6.1 indicates that the bubble collapse speed is faster in the earlier stage than that predicted by Chen & Mayinger (1992). After that, the collapse speed is

PAGE 145

1.0 K.'..:---------:=================:=:::iL, 0.8 0.6 0 ~ .. ... .. .. . . .. 0.0002 .. .. . .. ........ 0.0004 Fo Present Chen & Mayinger (1992) Akiyama (1973) 0.0006 0.0008 135 Figure 6.1: Comparison of the temporal change of bubble radius predicted by (6.1) Akiyama (1973) and the present simulation for Re = 174.81 Fr = l.O W e = 0.42 Ja = 16.43 Pr = 8.37 Pe = 1463.7, pi/ P v = 527, if v = 41; 'Y = R/ Ro Fo = o:,t/(2Ro)2.

PAGE 146

1 0 . \ --.. :::..: . . . .. . . 0 8 ' ' .... .... 0 6 ..... ...... .... ...... ...... ..... ..... 0.4 P r esent C h en & May i nge r (1 992 ) ........... Ak i yama ( 19 7 3 ) Florschuetz & Chao ( 1965 ) 0 2----------------------0 0 0002 0 0004 Fo 0 0006 0 0008 136 Figure 6.2: Comparison of the temporal change of bubble radius predicted b y ( 6.1) Florschuetz & Chao (1965), Akiyama (1973) and the present simulation for R e = 291.86 Fr = 1.0 W e = 0.13 Ja = 9.0 Pr = 1.72 Pe = 503.41 pi/ P v = 1605 if v = 22; 'Y = R/ !lo Fo = a,t/ (2flo) 2 slower during the majority of the collapse process than that of Chen & Mayinger (1992). Moreover present results are closer to those by Chen & Mayinger (1992) which were obtained from measurements Akiyama (1973) used the heat transfer rate similar to that around a moving solid sphere to obtain his analytical solution which probably resulted in a slower collapse rate than that of a deformed bubble. In Figure 6.2 the temporal variation of bubble radius for case 2 obtained by the present simulation, correlation of (6 1) and prediction of Akiyama (1973) as well as Florschuetz & Chao (1965) are plotted.

PAGE 147

137 From these two cases, the observation is that the present simulation is not in favorable agreement with the correlation proposed by Chen & Mayinger (1992) for the time evolution of the bubble radius of a collapsing and rising bubble. Chen & Mayinger (1992) also suggested that their data were better fitted by (6 1) for values of Pe/ Ja equal to several hundred. When the ratio of Pe/ Ja drops below 60, a significant disagreement is observed (see Legendre, Boree & Magnaudet 1998). This is indeed the outcome from our simulations, case 1 has a P e / Ja ratio over 60 so that a better match is observed. Case 2 has a ratio lower than 60 and shows substantial disagreement

PAGE 148

CHAPTER 7 SUMMARY A D FUTURE WORK 7.1 Summary of Present Work The objective of the present work has been to develop a predictive numerical capability for single bubble dynamics with phase change. In order to do this, the fixed-grid, sharp interface methodology has been developed (Ye et al., 2001b,a). The major components developed in this work for the entire numerical simulation method are A C 2 cubic spline curve fitting in conjunction with the fairing algorithm is im plemented to numerically represent, track the interface, and obtain the accurate geometric informations. The critical advantage of this scheme is by eliminat ing the oscillations in the curvature calculation, to enable the success of the computations since the curvature is the most notorious geometric parameter which often leads to the failure of the overall simulation of the interfacial flows in which the capillary force is important. This is the first time the efficient and capable geometric fairing algorithm is exploited and implemented to solve the curvature calculation problem in the area of computational fluid dynamics. A Cartesian grid based flow solver is developed for solving the coupled govern ing equations on a fixed grid with curved boundaries. Although this part is developed as an integrated component of the present overall method, it alone is an good alternative to body-fitted grid approach for solving flows over complex geometries. The major advantage of this method is that it eliminates the grid generation phase and thereby results in a more efficient solution procedure. 138

PAGE 149

139 The interface updating schemes and procedures are implemented and developed to allow the accurate prediction of the exact interface shape which satisfies all dynamical and thermal conditions. All major components of this method are tested or validated through compar isons with the established results and grid resolution study. The capability of the method is demonstrated in simulations of buoyancy-driven bubble motion bubble growth and bubble collapse under the gravity. These are extremely complex physical phenomena where multiple interacting effects are required to be modeled simultaneously in addition to a constantly evolving interface. 7. 2 Co nt ribu t ion s o f P resent Work 1 A fixed-grid direct numerical simulation method has been developed and verified for bubble dynamics, heat transfer and phase change phenomena. The original contribution of this research is the first introduction of a sharp interface between the vapor and liquid phases which allows the simulation of density ratios up to the order of 1000. The interface location is obtained as part of the solution which accounts for all coupled dynamic and thermal effects. 2 The current method is believed to be the only method available at present which resolves the curvature calculation issue treats the interface as a sharp discontinuit y, honors the mass conservation principle is capable of handling the large property ratios and phase change. 3. For the isothermal cases, the predictions of drag coefficients and the shape deformation of bubbles for Reynolds numbers ranging from 2 to 100 and We ber numbers from 2 to 20 were shown independently to agree with previous published results of Ryskin & Leal (1984).

PAGE 150

140 4. It was also demonstrated that the current method is capable of predi c ting accurate results for a wide range of Re We Ja Pe numbers and propert y jumps at the interface. 5. For a stationary and growing bubble, the predicted growth rate was found to agree with the theoretical limit of ex t 1 / 2 (see Appendix A) However for a rising and growing bubble, the predicted growth rate was approaching the theor e tical limit of ex t 2 1 3 (see Appendix A) for a spherical bubble but did not reach i t exactly owing to bubble deformation 6. The empirical correlation proposed for bubble collapse is assessed by the sim ulations in the present work. The present numerical model predicts tha t the dynamics of a collapsing bubble is qualitatively but not quantitativ e ly similar to the empirical correlation proposed by Chen & Mayinger (1992) 7.3 Future Work 7.3.1 Applications There are many important phenomena in our nature can be modeled as multi phase, multicomponent flows with dynamic interfaces such as fluid-fluid fluid-solid interaction problems. Compared to pure single phase materials and fluids th e nu merical simulations of multiphase flows with dynamic interfaces are among the most difficult computational problems because of the major challenge of resolving the in terfacial motion. In this research a numerical capability has been developed to begin probing this frontier. The method is applied to study the vapor bubble dynamics which is just one of examples of flows with dynamic interfaces. A target application of the present method is to investigate the flows in th e bio logical systems. In biomedical engineering investigations of the interaction of blood cells with elastic vessels are underway to explore the mechanical mechanisms of heart and vessel diseas e s. The study of blood flow in the small artificial organs helps d ev elop

PAGE 151

141 better devices with reduced blood damage when blood circulates through them For these types of problems, blood can be modeled as numerous deformable cells (largely red blood cells) suspended in a Newtonian fluid (plasma). The blood cells are bod ies containing fluid gel within a solid membrane. The blood damage occurs when those blood cells are broken due to unsustainable stress and deformation caused by interaction among the blood cells, plasma and the solid walls. The dynamics of the blood cell membranes are central to the cell damage process. In order to quantify the membrane stresses and the deformations, the microscopic computations at each cell level are necessary to resolve the numerous dynamically deforming cellular interfaces simultaneously. 7.3.2 Numerical Enhancements The present method is developed to resolve the interfacial dynamics within a fixed grid frame. The sharp interface treatment renders it well suitable for simulating the problems outlined in the above. The implementation of geometric algorithms can handle numerous interfaces as described in the Chapter 2. However, several enhancements have to be incorporated into the present numerical methodology The large number of blood cells needed to be covered in a computation so as to have practically meaningful results, plus adequate spatial resolution necessitate the use of massively parallel computers where thousands of processors can each compute a sub-portion of the problem. In this regard the present method has a significant advantage, that is the parallelization task is relatively easy owing to the static-grid frame The parallelization issues for fixed grid solver are reasonably well understood. On the other hand, even for single interface problem such as the single vapor bubble dynamics studied in this thesis research current overall solution procedure relies on extensive iteration at multiple levels The computing cycle often takes hours on single processor computer. The leap in modeling capability and quality of the present method is thereby expected by parallelization of the numerical algorithms.

PAGE 152

142 The fixed, Cartesian grid can be readily decomposed into sub-domains which are mapped to multiple processors. Each processor then operates independently on the portion assigned to it communicating only the much smaller shared border data when necessary The message passing approach is the primary model for writing parallel pro grams. The program execution is divided among a number of processes which need not be running on different processors. The process is a more general term referring to an independent instruction and data flow. Each process thereby can work on some local data. The data sharing between processes is achieved by the message passing that is, by explicitly sending and receiving data between processes. On modern parallel computers, a library of functions or subroutines called MPI (Message Passing Interface) is provided for programmer to insert into source code to perform data communication between processes. The programmer must explicitly divide data and work and manage the communications among the processes to achieve parallel execution. The merger /breakup module for handling the interface interaction and topologi cal change is needed for computing multiple interfaces problems. The merger/breakup function is implemented in the present code and yet to be tested. Turbulence might be important at high Reynolds number. The direct numerical simulation of the turbulence is still out of reach even with parallel computers and algorithms because of the extremely high spatial resolution requirement The LES (Large Eddy Simulation) is a good compromise which solves for the flow structures the grid size can resolve, and models those smaller than the grid size. Various modeling approaches have been proposed for those sub-grid eddies. The most reasonable model that is, dynamic sub-grid scale model of LES has been implemented in the research yet to be tested and incorporated into the present method.

PAGE 153

APPENDIX A HEAT TRANSFER-CO TROLLED LIMITS OF SPHERICAL BUBBLE GROWTH For a spherical bubble undergoing heat transfer-controlled growth, the overall energy balance for the evaporation at the bubble surface could be expressed as (A. 1 ) where R(t) is the instantaneous bubble radius, A is the latent heat, Pv is the vapor density, h is the heat transfer coefficient, ~T is the degree of superheat. To estimate h, Ranz & Marshall (1952) 's correlation, considering both the con duction and convection effects, may be used Nu 1 1 2 + 0.6Re2 Pr3 2+0.6 (u:R) Prl hR k (A.2) where the first term 2 represents the contribution from conduction, the second term 1 1 0.6Re2 Pr3 denotes the contribution by the convection. From Eq. (A.2), we have (A.3) A.1 Conduction Dominated Growth If only the first term of Eq. (A.3) is substituted for h, Eq. (A.l) becomes (A.4) 143

PAGE 154

144 Integration of Eq. (A.4) yields l R(t) = CU (A.5) where C includes all parameters that are not function of the time such as thermal properties, Prandtl number and degree of superheat, etc. (A.5) indicates that the conduction dominated growth rate is proportional to the square root of time. A.2 Convection Dominated Growth With only the second term of Eq. (A.3) substituted for h, Eq. (A.l) then becomes (A.6) Integration of Eq. (A.4) gives 2 R(t) = Ct3 (A.7) From (A.7), the convection dominated growth rate is proportional to the two thirds power of time.

PAGE 155

REFERENCES ABDELMESSIH, A.H., HOOPER F. C. & NANNGIA, S. 1972 Flow effects on bubble growth and collapse in surface boiling. Int. J. Heat Mass Transf. 15 115-120. AHLBERG, J. H. NILSON, E. N & WALSH, J. L. 1967 The Theory of Splin es and Their Applications. Academic Press. AKIYAMA M. 1973 Bubble collapse in subcooled boiling Bull. JSME 16 (93), 530575. ALEXIADES, V. & SOLOMON, A. D. 1993 Mathematical Modeling of Melting and Freezing Processes Philadelphia, PA: Taylor & Francis. ANDERSON, D. M. & McFADDEN, G. B. 1997 A diffuse-interface description of internal waves in a near-critical fluid. Phys. Fluids 9, 1870 1879. ANTANOVSKII L. K. 1995 A phase field model of capillarity. Phys. Fluids 7 747 753. BARRETT, R., BERRY, M., CHAN, T., DAMMEL, T., DONATO, J., DONGARRA J., EIJKHOUT v. Pozo, R. ROMINE, C. & VAN DER VORST, H. 1993 Tem plates for the solution of linear systems: Building blocks for iterative methods. Philadelphia, PA: SIAM. BHAGA, D. & WEBER, M. E. 1981 Bubbles in viscous liquids: shapes wakes and velocities. J Fluid Mech. 105, 61 85. BIRKHOFF, G., MARGULIES, R. S. & HORNING W. A. 1958 Spherical bubble growth. Phys. Fluids 1, 201 -2 04. BRACKBILL, J. U., KOTHE, D. B. & ZEMACH, C. 1992 A continuum method for modeling surface tension. J. Comput. Phys. 100 335-354. BRANDT, A. 1977 Multilevel adaptive solutions to boundary value problems. Math. Comput. 31, 333 390. BRIGGS, W. L. 1987 A Multigrid Tutorial. Philadelphia, PA: SIAM. CAGINALP, G. 1984 Applications of Field Theory to Statistical Mechanics Berlin: Springer-Verlag. CAREY, V. P. 1992 LiquidVapor Phase-Change Phenomena. Bristol, PA : Taylor & Francis. CHEN, L., GARIMELLA, S. V. REIZES, J. A. & LEONARDI, E. 1999 The devel opment of a bubble rising in a viscous liquid. J. Fluid Mech. 387, 61 96. 145

PAGE 156

146 CHEN, W. C., KLAUSNER, J. F. & MEI, R. 1995 A simplified model for predicting vapor bubble growth rates in heterogeneous boiling. J Heat Transfer 117, 976-980. CHEN, Y. M. & MAYINGER, F. 1992 Measurement of heat transfer at the phase interface of condensing bubbles. Int. J. Multiphase Flow 18, 877-890. CHORIN, A. J. 1968 Numerical solution of the Navier-Stokes equations. Math. Com put. 22, 745 762. DANDY, D. S. & LEAL, L. G. 1989 Buoyancy-driven motion of a deformable drop through a quiescent liquid at intermediate Reynolds numbers. J. Fluid Mech. 208 161 192. DARBY, R. 1964 The dynamics of vapor bubbles in nucleate boiling. Chem Eng Sci. 19, 39 49. DELHAYE, J. M 1974 Jump conditions and entropy sources in two-phase systems. Local instant formulation. Int. J. Multiphase Flow 1, 395 409. DENNIS, S. C. R. & CHANG, G.-Z. 1970 Numerical solution for steady flow past a circular cylinder at reynolds number up to 100. J. Fluid Mech. 42, 471. DERGARABEDIAN, P. 1953 The rate of growth of vapor bubbles in superheated water. J. Appl. Mech. 20, 537 545. DIMIC, M. 1977 Collapse of one-component vapour bubble with translatory motion. Int. J. Heat Mass Transf. 20, 1325 1332. DONNE, M. D. & FERRANTI, M. P. 1975 The growth of vapor bubble in super heated sodium Int. J. Heat Mass Transf. 18 477 493. FARIN, G. 1997 Curves and Surfaces for Computer-Aided Geometric Design, A Prac tical Guide, 4th edn. San Diego: Academic Press FERZIGER, J. H. & PERIC, M. 1996 Computational Methods for Fluid Dynamics New York: Springer-Verlag. FLORSCHUETZ, L. W. & CHAO, B. T. 1965 On the mechanics of vapor bubble collapse J. Heat Transfer 87, 209 220. FLORSCHUETZ, L. W., HENRY, C. L. & KHAN, A. R. 1969 Growth rates of free vapour bubbles in liquids at uniform superheats under normal and zero gravity conditions. Int. J. Heat Mass Transf. 12, 1465 1489. FORSTER, H. K. & ZUBER, N. 1954 Growth of a vapor bubber in a superheated liquid J. Appl. Phys. 25, 474-478. FORSTER, H. K. & ZUBER, N. 1955 Dynamics of vapor bubbles and boiling heat transfer. AIChE J. 1 531. GOLDSTEIN, D. HANDLER, R. & SIROVICH, L. 1995 Direct numerical simulation of turbulent flow over a modeled riblet covered surface. J. Fluid Mech. 302, 333 376. GOLUB, G. H. & LOAN, C. F. V. 1989 Matrix Computations. Baltimore : Johns Hopkins University Press.

PAGE 157

147 GUMEROV N A 1996 The heat and mass transfer of a vapour bubble with t rans lator y motion at high Nusselt numbers Int J. Multiphase Flow 22 259 27 2. HARPER J. F. & MOORE D. W. 1968 The motion of a spherical liquid drop at high Reynolds number. J Fluid Mech. 32 367 391. HIRT, C. W. & NICHOLS B. D. 1981 Volume of Fluid (VOF) method for the dynamics of fr e e boundaries. J. Comput. Phy s. 39 201 225 ISHII, M. 1975 Thermo-Fluid Dynamic Th e ory of Two-Phas e Flow. Eyrolles Paris JACKSON C. P 1987 A finite element study of the onset of vortex shedding in flow past variously shaped bodies J Fluid Mech. 182, 23. JACQMIN D. 1999 Calculation of two-phase Navier-Stokes flows using phas efield modeling. J. Comput. Phys 155 96 127. JuRIC D & TRYGGVAS0N G. 1996 A front-tracking method for dendritic so lidi fication J. Comput. Phys. 123 127 148 JURIC D & TRYGGVASON, G. 1998 Computations of boiling flows Int J Multi phase Flow 24 (3) 387 410 KAN, H. c., UDAYKUMAR, H s. SHYY, w & TRAN-SON-TAY, R. 1998 Hy drodynamics of a compound drop with application to Leukocyte modeling. Phys. Fluids 10 (4) 760 774 KANG, I. S. & LEAL, L G 1987 Numerical solution of axisymmetric unst e ad y free boundary problems at finite Reynolds number. I. finite-differen c e schem e and its applications to the deformation of a bubble in a uniaxial straining flow. Phy s Fluids 30 1929 1940. KATAOKA I. 1986 Local instant formulation of two-phase flow. Int. J Mult i phas e Flow 12 7 45 758. KIM J & Mo IN P 1985 Application of a fractional step method to incompr es sibl e Navier-Stokes equations J Comput. Phys 59 308 323. KLAUSNER J. F. MEI R. & BERNHARD D. M. 1993 Vapour bubble departur e in forced-convection boiling Int J Heat Mass Transf e r 36 651 662. KOBAYASHI R. 1993 Modeling and numeri c al simulation of dendritic crystal growth. Physica D 63 410 423 KOTHE, D. B. & MJOLSNESS, R. C. 1992 Ripple : A new model for incompressible flows with free surfaces AIAA J 30 2694 2700. LEGENDRE D ., BOREE J. & MAGNAUDET, J. 1998 Thermal and dynami c evo lution of a spherical bubble moving steadily in a superheated or subcooled liquid. Phys. Fluids 10 (6), 1256 1272. MEI R. CHEN W & KLAUSNER J. F. 1995 Vapor bubble growth in heteroge neous boiling I IL Int. J. Heat Mass Transf e r 38 909 934. MOALEM D. & SIDEMAN, S. 1973 The effect of motion on bubble collapse. Int. J Heat Mass Transf. 16 2321 2329.

PAGE 158

148 MOORE, D. W. 1963 The boundary layer on a spherical bubble. J. Fluid Mech. 16, 161-176. OKHOTSIMSKII, A. D. 1988 The thermal regime of vapor bubble collapse at different Jakob numbers. Int. J. Heat Mass Transf. 31, 1569 1576. OSHER, S. & SETHIAN, J. A. 1988 Fronts propagating with curvature dependent speed: Algorithms based in Hamilton-Jacobi formulations. J. Comput. Phys. 79 12 -4 9. PARLANGE, J. Y. 1969 Spherical-cap bubbles with laminar wakes. J. Fluid Mech. 37, 257-263. PEMBER, R. B ., BELL, J.B., COLELLA, P., CRUTCHFIELD, W. Y. & WELCOME, M. L. 1995 An adaptive Cartesian grid method for unsteady compressible flow in complex geometries. J Comput. Phys. 120, 278 304. PESKIN, C. S. 1977 Numerical analysis of blood flow in the heart. J. Comput. Phys. 25 220-252. PLESSET, M. S. & PROSPERETTI, A. 1977 Bubble dynamics and cavitation. Ann. Rev. Fluid. Mech. 9, 145 185. PLESSET, M. S. & ZWICK, S. A. 1954 The growth of vapour bubbles in superheated liquids. J. Appl. Phys. 25, 493-500. PRESS, W. H., TEUKOLSKY, S. A., VETTERLING, W. T. & FLANNERY, B. P. 1992 Numerical Recipes in FORTRAN The Art of Scientific Computing Cam bridge, UK: Cambridge University Press. PROSPERETTI, A. & PLESSET, M. S. 1978 Vapour-bubble growth in a superheated liquid. J. Fluid Mech 85, 349-368. PROVENSAL, M., MATHIS, C. & BOYER, L. 1987 Benard-von Karman instability: Transient and forced regimes. J. Fluid Mech. 182, 1. QUIRK, J. 1994 An alternative to unstructured grids for computing gas dynamics flows around arbitrarily complex two-dimensional bodies Computers Fluids 23(1), 125-140. RANZ, W. E. & MARSHALL, W. R. 1952 Evaporation from drops. Chem Eng. Prag. 48 (3), 141, 173. RHIE, C. M. & CHOW, W. L. 1983 Numerical study of turbulent flow past an airfoil with trailing edge separation. AIAA Journal 21, 1525 1532. RIDER, W. J. & KOTHE, D. B. 1998 Reconstructing volume tracking. J. Comput. Phys. 141, 112 152. RuCKENSTEIN E. & DAVIS E. J. 1971 The effect of bubble translation on vapour bubble growth in superheated liquid. Int. J. Heat Mass Transf. 14, 939 952. RYSKIN, G. & LEAL L. G. 1984 Numerical solution of free boundary problems in fluid mechanics. Part I, II, III. J. Fluid Mech. 148, 1 -4 3. SADHAL, S. S., AYYASWAMY P. S. & CHUNG, J. N. 1996 Transport Phenomena with Drops and Bubbles. Springer-Verlag.

PAGE 159

149 SAFFMAN, P. G. 1956 On the rise of small air bubbles in water. J. Fluid Mech. 1 249-275. SAPIDIS, N. & FARIN, G. 1990 Automatic fairing algorithm for b-splin curves. Computer Aided Design 22 (2), 121 129. SCARDOVELLI, R. & ZALESKI, S. 1999 Direct numerical simulation of free-surface and interfacial flow. Annu. Rev. Fluid. Mech. 31 567 584. SCRIVEN, L. E. 1959 On the dynamics of phase growth. Chem. Eng. Sci. IO 1 13. SETHIAN, J. A. 1990 Numerical algorithms for propagating interfaces: Hamilton Jacobi equations and conservation laws. J. Diff. Geom. 31 131 161. SHYY, W. 1994 Computational Modeling for Fluid Flow and Interfacial Transport. Amsterdam, Netherlands: Elsevier. SHYY, W., UDAYKUMAR, H. S., RAO, M. M. & SMITH, R. W. 1996 Com putational Fluid Dynamics with Moving Boundaries. Philadelphia, PA: Taylor & Francis. SUSSMAN, M., FATEMI, E., SMEREKA, P. & OSHER S. 1998 An improved level-set method for incompressible two-phase flows. Comput. Fluids. 27 663 680. SUSSMAN, M., SMEREKA, P. & OSHER, S. 1994 A level-set approach for comput ing solutions to incompressible two-phase flow J. Comput. Phys. 114, 146-159. TAYLOR, T. D. & AcRivos, A. 1964 On the deformation and drag of a falling viscous drop at low Reynolds number. J. Fluid Mech. 18 466 476. TOKADA, N., YANG, W. J. & CLARK, J. A. 1968 Dynamics of moving gas bubbles in injection cooling. J. Heat Transfer 90 371-378. TRITTON, D. J. 1959 Experiments on the flow past a circular cylinder at low reynolds number. J. Fluid Mech. 6 547. UDAYKUMAR, H. S., KAN, H. C., SHYY, W. & TRAN-SON-TAY, R. 1997 Mu lti phase dynamics in arbitrary geometries on fixed Cartesian grids. J. Comput. Phys. 137 366-405. UDAYKUMAR, H. S. & SHYY, W. 1995 A grid-supported marker particle scheme for interface tracking. Num. Heat Transf. B 27 (2), 127-153. UNVERDI, S. & TRYGGVASON, G. 1992 A front tracking method for viscous incom pressible flows. J. Comput. Phys. 100 25-37. VAN DER VoRST, H. 1992 Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of non-symmetric linear systems. SIAM J. Sci. Statist. Comput. 13, 631-644. WANNIER, G. H. 1950 A contribution to the hydrodynamics of lubrication. Quart. Appl. Math. 8 1 10. WIESELSBERGER, C. 1922 New data on the laws of fluid resistance. Tech. Rep .. NACA TN 84.

PAGE 160

150 WILLIAMSON, C. H. K. 1996 Vortex Dynamics in the Cylinder Wake. Ann. Rev. Fluid Mech. 28 477 539 WITTKE, D. D. & CHAO, B. T 1967 Collapse of vapor bubbles with translatory motion. J. Heat Transfer 89, 17 24. YE, T., MITTAL R., UDAYKUMAR, H S. & SHYY, W. 1999 An accurate Cartesian grid method for viscous incompressible flows with complex immersed boundaries. J. Comput. Phys. 156, 209 240. YE, T., SHYY, W. & CHUNG J. N. 2001a A fixed-grid, moving boundary method for bubble dynamics with phase change 4th International Conference on Multi phase Flow (ICMF-2001}, New Orleans LA May 27-June 1 2001 YE, T., SHYY, W. & CHUNG, J. N. 2001b A fixed-grid, sharp-interface method for bubble dynamics and phase change. Submitted to J. Comput. Phys . ZANG, Y., STREET, R. L. & KOSEFF, J. R. 1994 A non-staggered grid frac tional step method for time-dependent incompressible Navier-Stokes equations in Curvilinear coordinates. J. Comput. Phys. 114, 18 33. ZENG, L. Z., KLAUSNER, J. F. BERNHARD D. M. & MEI, R. 1993a A unified model for the prediction of bubble detachment diameters in boiling systems II. flow boiling. Int. J. Heat Mass Transfer 36, 2271 2279. ZENG, L. Z., KLAUSNER, J. F. & MEI, R. 1993b A unified model for the prediction of bubble detachment diameters in boiling systems I. pool boiling. Int. J. Heat Mass Transfer 36, 2261 2270.

PAGE 161

BIOGRAPHICAL SKETCH Tao Ye was born in Jiangsu China. After completing 11-years primar y and secondary education he moved to Hangzhou where he received his bachelor s and master's degree in engineering mechanics from Zhejiang University At Zhejiang Uni versity, he conducted his thesis research on the numerical modeling of 3-D turbulent flow. He then worked as a lecturer in Shanghai Jiao Tong University (SJTU) Shang hai, and later as technical staff in Siemens (China) Ltd., Shanghai, for a few y ears. In the fall of 1997, he came to University of Florida to pursue his doctorate Tao Ye can be reached at: tao_ye_2000 @ yahoo.com. 151

PAGE 162

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scop e and quality, as a dissertation for the degree of Doctor of Philosophy Jacob N. Chung Chair Professor of Mechanical Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scop e and quality, as a dissertation for the degree of Doctor of Philosophy William G. Tiederman Professor of Mechanical Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scop e and quality, as a dissertation for the degree of Doctor of Philosophy J a s F. Klausner Pr ssor of Mechanical Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate in scop e and quality as a dissertation for the degree of Doctor of Philosophy Wei Shyy Professor of Aerospace Engineering, Mechanics and Engineering Science

PAGE 163

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate in scope and quality as a dissertation for the degree of Doctor of Philosophy l< MJv mcu 1 Ren-Wei Mei Associate Professor of Aerospace Engineering, Mechanics and Engineering Scienc e This dissertation was submitted to the Graduate Faculty of the College of En gineering and to the Graduate School and was acc~~.,---~-artial fulfillment of the requirements for the degree of Doctor of Philo y May 2001 f M. J. Ohanian Dean, College of Engineering Winfred M Phillips Dean, Graduate School

PAGE 164

L UNIVERSITY OF FLORIDA II I II IIIIII Ill Ill lllll lllll II IIIIII IIII IIII II 1111111111111111111 3 1262 08555 3575


xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID E1908GZLX_PNTIKK INGEST_TIME 2014-05-23T23:21:27Z PACKAGE AA00020455_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES