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.