MODELING TWO-PHASE TRANSPORT DURING CRYOGENIC CHILLDOWN IN
A PIPELINE
By
JUN LIAO
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
2005
Copyright 2005
by
Jun Liao
ACKNOWLEDGMENTS
I would like to express my appreciation to all of the individuals who have assisted
me in my educational development and in the completion of my dissertation. My greatest
gratitude is extended to my supervisory committee chair, Dr. Renwei Mei. Dr. Mei's
excellent knowledge, boundless patience, constant encouragement, friendly demeanor,
and professional expertise have been critical to both my research and education. Dr.
James F. Klausner also deserves recognition for his knowledge and technical expertise. I
would like to further thank Dr. Jacob N. Chung for kindly providing his experiment data
of chilldown.
I would like to additionally recognize my fellow graduate associates Christopher
Velat, Jelliffe Jackson, Yusen Qi, and Yi Li for their friendship and technical assistance.
Their diverse cultural background and character have provided an enlightening and
positive environment. Special appreciation is given to Kun Yuan for his kindness
providing his experiment data and insight on chilldown.
I would like to further acknowledge the Hydrogen Research and Education
Program for providing funding to this study. This research was also funded by NASA
Glenn Research Center under contract NAG3-2750.
Finally, I would like to recognize my wife Xiaohong Liao and my parents for their
continual support and encouragement.
TABLE OF CONTENTS
page
A C K N O W L E D G M E N T S ....................................................................... .....................iii
LIST OF TABLES .............. ................ .......... ............... ............. vii
LIST OF FIGURES ............... ............... ......................... ........... viii
N O M E N C L A T U R E ........................................................................ ........................... xiv
A B S T R A C T ...............................x..................................................... .. x ix
CHAPTER
1 IN TR O D U C T IO N ............ ............................................. ........ .. ........... .. 1
1.1 B background ...................... ............... ........... ... ............................... . 1
1.2 L literature R eview ......................................................... ............ .. ............ 4
1 .3 S c o p e ...................... .. ............. .. ..................................................... 1 0
2 TWO-PHASE FLOW MODELING AND FLOW BOILING HEAT TRANSFER OF
CRY O GEN IC FLU ID.............................................................. .......................... 13
2.1 Flow Regim e and Heat Transfer Regim e.................................. .................... 13
2.2 Flow Models in Cryogenic Chilldown.................................................... 18
2.2.1 Hom ogeneous Flow M odel.................................................................... 18
2 .2 .2 T w o-F luid M odel ................................................................. ............... ... 22
2.3 Heat Transfer between Cryogenic Fluid and Solid Pipe Wall............................ 26
2.3.1 Heat Transfer between Liquid and Solid wall ....................................... 27
2.3.1.1 Film boiling ................. .............. .. .. .. ...... ..... .. ..... 27
2.3.1.2 Forced convection boiling and two-phase convective heat transfer 30
2.3.2 Heat Transfer between Vapor and Solid Wall................. .................... 33
3 VAPOR BUBBLE GROWTH IN SATURATED BOILING.............. ................ 34
3 .1 In tro d u ctio n ............ .. .............. .............................. ............... 3 4
3.2 Formulation ................ ......... ......... ........ 39
3.2.1 On the Vapor Bubble ............ .... ......... ...................... 39
3.2.2 M icrolayer............................. .............. 41
3.2.3 Solid Heater .......................................... 42
3.2.4 O n the B ulk Liquid .............. ........................................................... 43
3.2.4.1 V elocity field ... ..................................................... ... .. ....... .. 43
3.2.4.2 Tem perature field .............................. ............. ............... 44
3.2.4.3 Asymptotic analysis of the bulk liquid temperature field during early
stages of growth ...... ........................ .......... .............. 46
3.2.5 Initial Conditions .............. ..... ............. ....... ............ ................. 50
3.2.6 Solution Procedure............................................... .......................... 52
3.3 R results and D iscu ssions ................................................................ ............... ... 52
3.3.1 Asymptotic Structure of Liquid Thermal Field ......................... ............. 52
3.3.2 Constant Heater Temperature Bubble Growth in the Experiment of
Y addanapudi and K im .......................................................... .................. 56
3.3.3 Effect of Bulk Liquid Thermal Boundary Layer Thickness on Bubble
G row th .................. ............................... ....... ......... ...... 5 9
3 .4 C o n clu sio n s ............... .............................................................. 6 3
4 ANALYSIS ON COMPUTATIONAL INSTABILITY IN SOLVING TWO-FLUID
MODEL ........ ........... ........................... 64
4.1 Inviscid T w o-Fluid M odel .............. ................................................................ 65
4 .1.1 Introdu action .......................................... ...... ............ 65
4.1.2 G governing E quations ........................................... .................... ...... 67
4.1.3 Theoretical A analysis ..................................................................... 69
4.1.3.1 Characteristic analysis and ill-posedness ...................................... 69
4.1.3.2 Inviscid Kelvin-Helmholtz (IKH) analysis and linear instability.... 72
4.1.4 Analysis on Computational Instability ............................................... 73
4.1.4.1 Description of numerical methods........................... .... .............. 73
4.1.4.2 Code validation- dam-break flow ........................................ 78
4.1.4.3 Von Neumann stability analysis for various convection schemes .. 81
4.1.4.4 Initial and boundary conditions for numerical solutions............... 86
4.1.5 R results and D discussion .............................................................. ... 87
4.1.5.1 Computational stability assessment based on von Neumann stability
analysis ............................................................. 87
4.1.5.2 Schem e consistency tests........................................................... 94
4.1.5.3 Computational assessment based on the growth of disturbance...... 95
4.1.5.4 Discussion on the growth of short wave................................ .. 101
4.1.5.5 Wave development resulting from disturbance at inlet............... 104
4 .1.6 C onclu sions........................................... .................... .. .............. 106
4.2 V iscous Tw o-Fluid M odel .......................................... ........................... 110
4.2.1 Introduction..................................... 110
4.2.2 Governing Equations ................................... .. ...................... 111
4.2.3 Theoretical Analysis ...................................... .............. 112
4.2.3.1 Characteristics and ill-posedness............................... .............. 112
4.2.3.2 Viscous Kelvin-Helmholtz (VKH) analysis and linear instability 113
4.2.4 Analysis on Computational Intability ............ ............ ............. 115
4.2.4.1 Description of numerical methods.................. .. ................ 115
4.2.4.2 Von Neumann stability analysis for various convection schemes 116
4.2.4.3 Initial and boundary conditions for numerical solution............... 119
4.2.5 Results and Discussion ................... ........ .............. 119
4.2.5.1 Computational stability assessment based on von Neumann stability
a n a ly sis ............... ... ..... ............. ..... .... .. ................... 1 1 9
4.2.5.2 Computational assessment based on the growth of disturbance.... 126
4.2.5.3 Wave development resulting from disturbance at inlet............... 128
4 .2 .6 C on clu sion s........................................... .. ..................... ................. 130
5 MODELING CRYOGENIC CHILLDOWN..................................... 133
5.1 Homogeneous Chilldown Model .......................................................... 133
5.1.1 Analysis ....... ....... ............................ 134
5.1.2 R results and D discussion ...................................... ................. ... ........... 136
5.2 Pseudo-Steady Chilldown M odel... .................... .. .................. 140
5.2 .1 F orm ulation .................. ... ......................... ..... ........... .. ............ 14 1
5.2.1.1 H eat conduction in solid pipe ............. ................ ................... 141
5.2.1.2 Liquid and vapor flow ......................................... .............. 144
5.2.1.3 Film boiling correlation ...... .......... ....................... ............. 145
5.2.1.4 Forced convection boiling correlation................. .................... 151
5.2.1.5 Heat transfer between solid wall and environment........................ 152
5.2.2 Results and Discussion .................................. 155
5.2.2.1 Experim ent of Chung et al .............. ......................... ................ 156
5.2.2.2 Comparison of pipe wall temperature .......................................... 157
5.2.3 D discussion and R em arks .................................... ................................... 163
5.2 .4 C onclu sion s............................................. .. .............. .......... ....... 166
5.3 Separated Flow Chilldown M odel .............. ...... ......................................... 167
5 .3 .1 F o rm u latio n ........................................ ... ............... ....... ................ .. 16 7
5 .3 .1.1 F lu id fl ow .... .................. .. .................................. .............. 16 8
5.3.1.2 H eat conduction in solid pipe ............. ................ ................... 168
5.3.1.3 Heat and mass transfer............................. .............. 169
5.3.1.3 Initial and boundary conditions ................................... ............... 172
5.3.2 Solution Procedure............................. .............. 173
5.3.3 Results and Discussion ........................................ 174
5.3.3.1 Comparison of solid wall temperature................ .......... .... 177
5.3.3.2 Flow field and fluid tem perature .................................................. 181
5.3 .4 C onclu sion s....................................... .. ................. ...... ... .. ..... 186
6 CONCLUSIONS AND DISCUSSION ............................................................. 187
6 .1 C o n clu sio n s .............................................................................. 18 7
6.2 Suggested Future Study ....................................... .............. 188
LIST O F REFEREN CE S ................................................... ................................. 190
B IO G R A PH IC A L SK E T C H .......................................................................................... 198
LIST OF TABLES
Table page
4-1. Analytical solution for dam-break flow.......................................... .............. 80
4-2. A (O) for different discretization schem es............................................ ... ... .............. 85
5-1. Heat and mass transfer relationship used in separated flow chilldown model........ 173
LIST OF FIGURES
Figure page
1-1. Schematic of filling facilities for LH2 transport system from storage tank to space
shuttle external tank............................................. ................ ......... 3
1-2. The schematic of chilldown and heat transfer regime. ........................................... 9
2-1. Schematic of two-phase flow regime in horizontal pipe. ...................................... 14
2-2. Schematic of two-phase flow regime in vertical pipe............................ ............ 14
2-4. Typical wall temperature variation during chilldown............................ ............ 17
2-5. Schem atic for hom ogeneous flow m odel.............................................. ... ................. 19
2-6. Schematic of the two-fluid model ................................................................. 22
2-7. Schem atic of heat transfer in chilldow n.............................................. .... .. .............. 27
3-1. Sketch for the growing bubble, thermal boundary layer, microlayer and the heater
w all ................................. ...................... ........ ......... ...... 39
3-2. Coordinate system for the background bulk liquid................ ............. .............. 43
3-3. A typical grid distribution for the bulk liquid thermal field with S,, = 0.65,
S = 0.73 and R = 10 ............ ....... .................................. ...... .. ..... 46
3-4. Comparison of the asymptotic and the numerical solutions at '=0.001, 0.01, 0.1 and
0.3 for 0 40 and 7 1 ............................................................................... 54
3-5. Effect of parameter A on the liquid temperature profile near bubble. ..................... 55
3-6. The computed isotherms near a growing bubble in saturated liquid at T=0.01,
.1, 0 .3, and 0 .9 .................... ................ ........................ .............. 57
3-7. Comparison of the equivalent bubble diameter d,, for the experimental data of
Yaddanapudi and Kim (2001) and that computed for heat transfer through the
m icrolayer (c, =3.0). ........ ...................... ................ ...................... ........ ........ .... 58
3-8. Comparison of bubble diameter, d(t), between that computed using the present
model and the measured data of Yaddanapudi and Kim (2001) ............. ............. 60
3-9. Comparison between heat transfer to the bubble through the vapor dome and that
through the m icrolayer ..................................................................................... 60
3-10. The computed isotherms in the bulk liquid corresponding to the thermal conditions
reported by Yaddanapudi and Kim (2001). ....................................................... 61
3-11. Effect of bulk liquid thermal boundary layer thickness Son bubble growth........... 62
4-1. Schematic of two-fluid model for pipe flow......... ........................................... 68
4-2. Staggered grid arrangement in two-fluid model ........................................... 74
4-3. Flow chart of pressure correction scheme for two-fluid model................................. 78
4-4. Schematic for dam-break flow model ........................... .......... ............. ................. 79
4-5. Water depth at t=50 seconds after dam break.................................. .............. 80
4-6. Water velocity at t=50 seconds after dam break .................. ............. .............. 81
4-7. Grid index number in staggered grid for von Neumann stability analysis ............ 81
4-8. Comparisons of growth rates of various numerical schemes. N = 200, a, = 0.5,
u, = Im /s, Ug = 17m /sand CFL, = 0.1 ........... ............... ...... ............ ........ 89
4-9. Growth rate of CDS scheme at different AU = Ug -u,. N = 200, a, = 0.5,
u = Im /s and C F L = 0.1 ........................... .................................. ............ ..... 90
4-10. Growth rate of FOU scheme at different AU = Ug u. N = 200, a, = 0.5,
uz = m / s and CFL1 = 0.1 ........................... .............. .................................... 90
4-11. Growth rate of CDS scheme at different u,. N = 200, AU = 16m/s, a, = 0.5, and
A t/A x = O s /m .............. ...... ............................................................... 9 2
4-12. Growth rate of FOU scheme at different u,. N = 200, AU = 16m/s, a, = 0.5 and
A t/A x = O s /m .............. ...... ............................................................... 9 2
4-13. Growth rate of CDS scheme at different At/Ax N = 200, u = Im / s,
A U = 16m /s an d a = 0 .5 .................................................................................. 93
4-14. Growth rate of FOU scheme at different At/Ax N = 200, u, = Im / s,
A U = 16m /s and a = 0.5 ............................................................................... 93
4-15. Comparison of fi growth using CDS scheme on different grids. u, = lm/s,
Ug =17.5m/s, CFL, = 0.1, and a = 0.5 ................................... .............. 95
4-16. ti using CDS scheme in the computational domain. N = 200, u, =m I/s,
AU =14m /s, CFL, = 0.05, a, = 0.5, and t =4s .............................................. 97
4-17. Amplitude of liquid velocity disturbance i, using CDS scheme. N = 200,
u = Im/s, AU = 14m/s, CFL, = 0.05, a, = 0.5, and t = 4s ............................ 97
4-18. fi using CDS scheme after 10399 steps of computation, N = 200, u = Im / s,
AU =16.5m/ s, CFL, = 0.1, a, = 0.5 and t = 5.2s. ......................................... 98
4-19. Growth history of iI solved using CDS scheme, N = 200, u, = Im/s,
AU =16.5m/ s, CFL, = 0.1, a, = 0.5 and t = 5.2s. ......................................... 98
4-20. Growth rate of FOU scheme, N = 200, u, = 0.5m/s, AU = 16m/s, CFL, = 0.02,
and a, = 0.5 ....................... .................................................. 100
4-21. ti using FOU scheme after 12000 steps of computation. N = 200,
u, = 0.5m/s,AU =16m/s, CFL, = 0.02, and a, = 0.5 ................................ 102
4-22. Growth rate of SOU scheme. N= 200, u, = Im/s, AU = 16m/s, CFL, = 0.05,
and a, = 0.5 ....................... .................................................. 103
4-23. il using SOU scheme after 3000 steps of computation. N = 200, u, =m / s,
AU = 16m / s, CFL, = 0.05, and a, = 0.5 ............ .............................. .... 103
4-24. Growth history of fi under different initial amplitude using FOU scheme........... 104
4-25. ti propagates in the pipe with FOU at well-posed condition, quasi-steady state. 107
4-26. fi propagates in the pipe with FOU scheme at ill-posed condition, quasi-steady
state ........................................................................................... 1 0 7
4-27. fi propagates in the pipe with CDS at well-posed condition, quasi-steady state. 108
4-28. ti propagates in the pipe with CDS at ill-posed condition, an instance before the
computation breaks down. ....................... ............................... 108
4-29. Comparison of growth rate between CDS and FOU schemes. N = 200, u, = Im / s,
Ug = 21m /s, CFL, = 0.05, and a = 0.5 ............................................................ 109
4-30. Schematic depiction of viscous two-fluid model........................... .... ................. 111
4-31. Comparisons of growth rate of different schemes. N = 200, u, = 0.3m I s,
Ug, = 66 m /s ,and C F L = 0. 1 ......... .................................................................... 120
4-32. Comparisons of growth rate of different schemes at low k. N = 200, u, = 0.3m / s,
U g, = 6 m /s ,and C F L = 0. 1 ......... .................................................................... 12 1
4-33. Growth rate for CDS scheme with VKH unstable. N = 200, us = 0.3m Is,
U g, = 6 m /s ,and C F L = 0. 1 ......... .................................................................... 122
4-34. Growth rate for CDS scheme with VKH instability. lwate = 10-2Pa s, N= 200,
u, = 0.3m / s, Ug = 6m /s ,and CFL, = 0.1 ....... ... ..................... ............. .. 123
4-35. Growth rates for CDS scheme with VKH instability. ,lwae = 10 -Pa s, N= 200,
u, = 0.1m/s Ug = 2m/ s,and CFL, = 0.01 ......................................... 124
4-36. Growth rates for FOU scheme with VKH instability. N = 200, u,, = 0.3m /s,
U g = 6 m /s and C F L = 0 .1 ................................................................................ 12 5
4-37. Growth rates for FOU scheme with VKH instability. ,l wae = le 1Pa s,
N =200, u, = 0.1m/s Ug =2m/s, and CFL, =0.01. ............................... 125
4-38. Growth history of u, using CDS scheme. N = 200, u, = 2m/s,
Ug = 0.998174m/s, 8 = -0.0617144, a, = 0.98, and CFL, = 0.05................... 127
4-39. Growth history of i, using FOU scheme. N = 200, u, = 2m/s,
Ug = 0.998174m/s, 8 = -0.0617144, a, = 0.98, and CFL, = 0.05 ................. 128
4-40. Disturbance of i, propagates in the pipe with FOU and CDS schemes at VKH
unstable and w ell-posed condition. ............. ................................. .............. 129
4-41. Disturbance of i, propagates in the pipe with FOU and CDS schemes at both VKH
unstable and w ell-posed condition. ............. ................................. .............. 130
5-1. Schematic of homogeneous chilldown model. ............................................... 134
5-2. Schematic for evaluating film boiling wall friction ......................... ................. 135
5-3. Distribution of vapor quality based on the homogenous flow model...................... 137
5-4. Pressure distribution based on the homogenous flow model................................. 138
5-5. Velocity distribution based on the homogenous flow model.................................. 139
5-6. Solid temperature contour based on homogenous flow model............................ 139
5-7. Schematic of cryogenic liquid flow inside a pipe............................. .............. 141
5-8. Coordinate systems: laboratory frame is denoted using z; moving frame is denoted
using Z. ............................................. 142
5-9. Schematic diagram of film boiling at stratified flow........................ ............ 145
5-10. Numerical solution of the vapor thickness and velocity influence functions........ 150
5-11. N um erical solution of G ((0 ) ................................................................. ... 151
5-12. Schematic of vacuum insulation chamber. .................................. .............. 153
5-13. Schematic of Yuan and Chung (2004)'s cryogenic two-phase flow test apparatus. 157
5-14. Experimental visual observation of Chung et al. (2004)'s cryogenic two-phase flow
experiment. ................ ............................ ............... 158
5-15. Computational grid arrangement and positions of thermocouples ...................... 159
5-16. Comparison between measured and predicted transient wall temperatures of
positions 12 and 15 at the bottom of pipe during film boiling chilldown ......... 160
5-17. Comparison between measured and predicted transient wall temperatures of
positions 12 and 15 at the bottom of pipe during convection boiling chilldown. 161
5-18. Comparison between measured and predicted transient wall temperatures of
positions 12 and 15 is at the bottom of pipe during entire chilldown. .................. 161
5-19. Comparison between measured and predicted transient wall temperatures of
positions 11 and 14, which is at the bottom of pipe during entire chilldown........ 162
5-20. Cross section wall temperature distribution at t=0, 50, 100 and 300 seconds....... 162
5-21. Computed wall temperature contour on the inner surface of inner pipe.............. 163
5-22. Comparison between measured and predicted transient wall temperatures of
positions 12 and 15 with Chen correlation (1966). .......................................... 165
5-23. Schematic of separated flow chilldown model ................................ ................. 168
5-24. Schematic of heat and mass transfer in separated flow chilldown model ........... 169
5-25. Flow chart of separated flow chilldown model................................................ 175
5-26.Geometry of the test section and locations of thermocouples .............................. 176
5-27. Comparison between measured and predicted transient wall temperatures of
positions 12 and 15. .... ........... .......................... ..................... .... .......... 178
5-28. Comparison between measured and predicted transient wall temperatures of
positions 11 and 14. ..................................................... ............. 178
5-29. Comparison between measured and predicted transient wall temperatures of
position 6 and 9 (the measured T 9 is obviously incorrect).............................. 179
5-30. Comparison between measured and predicted transient wall temperatures of
positions 5 and 8. ....................................................... .............. 179
5-31. Comparison between measured and predicted transient wall temperatures of
position 3 and the numerical result of temperature at position 4 ....................... 180
5-32. Comparison between measured and predicted transient wall temperatures of
positions 1 and 2. ....................................................... .............. 180
5-33. Liquid nitrogen depth in the pipe during the chilldown. .............. ............... 182
5-34. Vapor nitrogen velocity in the pipe during the chilldown .................................. 183
5-35. Liquid nitrogen velocity in the pipe during the chilldown ............... .............. 184
5-36. Vapor nitrogen temperature in the pipe during the chilldown............................ 185
5-37. Liquid nitrogen temperature in the pipe during the chilldown. .......................... 185
A
Ab
Am
Bo
C,
CFL
D
D, and Dg
d
deq
E
f
flo
G
g
H, and Hg
NOMENCLATURE
dimensionless parameter for bubble growth, cross section area, surface
area
area of vapor bubble dome exposed to bulk liquid
area of wedge shaped interface
Boiling number
ratio of wedge shaped interface radius and vapor bubble radius, wave
speed
microlayer wedge angle parameter; empirically determined
Courant number
diameter of pipe
liquid layer and gas layer hydraulic diameter
bubble diameter
equivalent bubble diameter
common amplitude factor
friction factor
friction factor for liquid phase in homogeneous model
mass flux, amplification factor
gravity
liquid layer and gas layer hydraulic depth
h, and hg liquid layer and gas layer depth
h
h B
pool
h,,c and hgc
hfg
heat transfer coefficient
film boiling heat transfer coefficient
pool boiling heat transfer coefficient
forced convection heat transfer coefficient for liquid and gas
latent heat of vaporization
imaginary unit, V-
enthalpy
Jacob number
thermal conductivity, wavenumber
effective thermal conductivity
Ja
k
kif
L
Nu
1th
n
P
Po
Pc
Pr
q
qrad
q frc
q'f
qW
local microlayer thickness, characteristic length
Nusselt number
mass transfer rate between liquid and gas per unit length
normal direction
pressure
pressure in the liquid-vapor interface
Peclet number
Prandtl number
heat transfer rate per unit length
radiation heat flux
free convection heat flux
Heat flux from wall to fluid
R
R', Wand
R'
Rb
Ro
Ra
Re
r
S
S, and S,
T
Ts,
T,.
Tb
t,
to
Uand V
u and v
u
vapor bubble radius, pipe radius
bubble growth rate
spherical coordinates
dimensionless radial coordinate
radius of wedge shaped interface
initial bubble radius
Rayleigh number
Reynolds number
radial coordinate
suppression factor in flow nucleate boiling, perimeter
stretching factor in computation
temperature
saturated temperature
initial solid temperature
bulk liquid temperature
time
characteristic time
waiting period
initial time
averaged velocities
velocities
mean u velocity
Vb vapor bubble volume
x, y, and z Cartesian coordinates
z, r, and (p cylindrical coordinates
X boundary layer coordinate
Z coordinate in the direction normal to the heating surface
Greek symbols
a thermal diffusivity, volume fraction
3 volumetric thermal expansion coefficient
Xtt Martinelli number
AT solid wall superheat
6 superheated bulk liquid thermal boundary layer thickness, vapor film
thickness
6* dimensionless thickness of unsteady thermal boundary layer
c emissivity, amplitude
P velocity potential function for liquid flow, general variable
Smicrolayer wedge angle, azimuthal coordinate, phase angle
oi friction multiplier
r7 and S computational coordinates
A characteristic root of a matrix
v kinematic viscosity
0 dimensionless temperature, azimuthal coordinate, pipe incline angle
o0 initial dimensionless temperature of liquid
p density
a
'FB
superscripts
in
out
subscripts
b
FB
eva
1, g, and i
i and o
/
ml
NB
w
v
CX3
stretched time in computation, Stefan Boltzmann constant
dimensionless time, shear stress
wall shear stress in film boiling regime
inner solution
outer solution
quantity per unit length
quantity per unit area
bubble
film boiling
evaporation
liquid, gas, and interface
inner and outer pipe
liquid
microlayer
nucleate boiling
wall
vapor
far field condition
xviii
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
MODELING TWO-PHASE TRANSPORT DURING CRYOGENIC CHILLDOWN IN
A PIPELINE
By
Jun Liao
August 2005
Chair: Renwei Mei
Major Department: Mechanical and Aerospace Engineering
Cryogenic chilldown process is a complicated interaction process among liquid,
vapor and solid pipe wall. To model the chilldown process, results from recent
experimental studies on the chilldown and existing cryogenic heat transfer correlations
were reviewed together with the homogeneous flow model and the two-fluid model. A
new physical model on the bubble growth in nucleate boiling was developed to correctly
predict the early stage bubble growth in saturated heterogeneous nucleate boiling. A
pressure correction algorithm for two-fluid model was carefully implemented to solve the
two-fluid model used to model the chilldown process. The connections between the
numerical stability and ill-posedness of the two-fluid model and between the numerical
stability and viscous Kelvin-Helmholtz instability were elucidated using von Neumann
stability analysis. A new film boiling correlation and a modified nucleate boiling
correlation for chilldown inside pipes were developed to provide heat transfer correlation
for chilldown model. Three chilldown models were developed. The homogeneous
chilldown model is for simulating chilldown in a vertical pipe. A pseudo-steady
chilldown model was developed to simulate horizontal chilldown. The pseudo-steady
chilldown model can capture the essential part of chilldown process, provides a good
testing platform for validating cryogenic heat transfer correlations based on experimental
measurement of wall temperature during chilldown and gives a reasonable description of
the chilldown process in a frame moving with the liquid-vapor wave front. A more
comprehensive separated flow chilldown model was developed to predict both the flow
field and solid wall temperature field in horizontal stratified flow during chilldown. The
predicted wall temperature variation matches well with the experimental measurement. It
provides valuable insights into the two-phase flow dynamics, and heat and mass transfer
for a given spatial region in the pipe during the chilldown.
CHAPTER 1
INTRODUCTION
One of the key issues in the efficient utilization of cryogenic fluids is the transport,
handling, and storage of the cryogenic fluids. The complexity of the problems results
from, in general, the intricate interaction of the fluid dynamics and the boiling heat
transfer. Chilldown of the pipeline for transport cryogenic fluid is a typical example. It
involves unsteady two-phase fluid dynamics and highly transitory boiling heat transfer.
There is very little insight into the dynamic process of chilldown. This study will focus
on the understanding and modeling of the unsteady fluid dynamics and heat transfer of
the cryogenic fluids in a pipeline that is exposed to the atmospheric condition.
1.1 Background
Presently there exists considerable interest among U.S. Federal agencies in driving
the U.S. energy infrastructure with hydrogen as the primary energy carrier. The
motivation for doing so is that hydrogen may be produced using all other energy sources,
and thus using hydrogen as an energy carrier medium has the potential to provide a
robust and secure energy supply that is less sensitive to world fluctuations in the supply
of fossil fuels. The vision of building an energy infrastructure that uses hydrogen as an
energy carrier is generally referred to as the "hydrogen economy," and is considered the
most likely path toward widespread commercialization of hydrogen based technologies.
Hydrogen has the distinct advantage as fuel in that it has the highest energy density
of any fuel currently under consideration, 120 MJ/kg. In contrast, the energy density of
gasoline, which is considered relatively high, is approximately 44 MJ/kg. When
launching spacecraft, the energy density is a primary factor in fuel selection. When
considering liquid hydrogen to propel advanced aircraft turbo engines, it is a very
attractive option due to hydrogen's high energy density. One drawback with using liquid
hydrogen as a fuel is that it's volumetric energy capacity, 8.4 MJ/liter is about one
quarter that of gasoline, 33 MJ/liter. Therefore, liquid hydrogen requires more
volumetric storage capacity for a fixed amount of energy. Nevertheless, liquid hydrogen
is a leading contender as a fuel for both ground-based vehicles and for aircraft propulsion
in the hydrogen economy.
When any cryogenic system is initially started, (this includes turbo engines,
reciprocating engines, pumps, valves, and pipelines), it must go through a transient
chilldown period prior to operation. Chilldown is the process of introducing the
cryogenic liquid into the system, and allowing the hardware to cool down to several
hundred degrees below the ambient temperature. The chilldown process is anything but
routine and requires highly skilled technicians to chilldown a cryogenic system in a safe
and efficient manner.
A perfect example of utilization and chilldown cryogenic system exists in NASA's
Kennedy Space Center (KSC). In the preparation for a space shuttle launch, liquid
hydrogen (as fuel) is filled from a storage tank to the main liquid hydrogen (LH2)
external tank (ET) through a complex pipeline system (Figure 1-1). The filling procedure
consists of 5 steps:
* Facility and orbiter chilldown.
* Fill transition and initial fill (fill ET to 2%).
* Fast fill ET (to 98%).
* Fill ET (to 100%).
* Replenish (maintain ET 100%).
- ,, ,- I ,I-F ,7r
n -I. -'
r;K;r
Calegor; III FaI Fill lo 9.8; ;
I I. .' .
4 ;7 -
F' FiL
|I,,: IL,
r ,, r
rF
F-
," r I- .F -
- .IJ
F F
, L
Figure 1-1. Schematic of filling facilities for LH2 transport system from storage tank to
space shuttle external tank.
While the engineers have a general understanding of the process in the initial fill
and rapid fill stages, there has been very little insight about the process of chilldown,
which is the first procedure to be initiated. There is not a single formula or computer
code that can be used to estimate the elapse time during the chilldown stage if certain
operating condition changes. The absence of guidelines stems from our lack of
fundamental knowledge in the area of cryogenic chilldown. Many such engineering
issues are present in the transport, handling, and storage of cryogenic liquid in industry
applications.
1.2 Literature Review
Experimental studies: Studies on cryogenic chilldown started in the 1960s with the
development of rocket launching systems. Early experimental chilldown studies started in
the 1960s by Burke et al. (1960), Graham (1961), Bronson et al. (1962), Chi and Vetere
(1963), Steward (1970) and other researches. Burke et al. (1960) and Graham (1961)
experimentally studied the cryogenic chilldown in a horizontal pipe and in a vertical pipe,
respectively. However, none of these studies provided the flow regime information in
chilldown. Bronson et al. (1962) visually studied the flow regimes in a horizontal pipe
during chilldown with liquid hydrogen as the coolant. The results revealed that the
stratified flow is prevalent during the cryogenic chilldown.
Flow regimes and heat transfer regimes in the horizontal pipe chilldown were also
studied by Chi and Vetere (1963). Information on flow regimes was deduced by studying
the fluid temperature and the volume fraction during chilldown. Several flow regimes
were identified: single-phase vapor, mist flow, slug flow, annular flow, bubbly flow, and
single-phase liquid flow. Heat transfer regimes were identified as single-phase vapor
convection, film boiling, nucleate boiling, and single-phase liquid convection.
Recently, Velat et al. (2004) systematically studied cryogenic chilldown with
nitrogen in a horizontal pipe. Their study included: a visual recording of the chilldown
process in a transparent Pyrex pipe, which is used to identify the flow regime and heat
transfer regime; collecting temperature histories at different positions of the wall in
chilldown; and recording the pressure drop along the pipe. Chung et al. (2004) conducted
a similar study with nitrogen chilldown at relatively low mass flux and provided the data
needed to assess various heat transfer coefficients in the present study.
Modeling efforts: Burke et al. (1960) developed a crude chilldown model based on
1-D heat transfer through the pipe wall and the assumption of infinite heat transfer rate from
the cryogenic fluid to the pipe wall. The effects of flow regimes on the heat transfer rate
were neglected. Graham et al. (1961) correlated the heat transfer coefficient and pressure
drop with the Martinelli number (Martinelli and Nelson, 1948) based on their
experimental data. Chi (1965) developed a one-dimensional model for energy equations
of the liquid and the wall, based on the film boiling heat transfer between the wall and the
fluid. An empirical equation for predicting the chilldown time and the temperature was
proposed.
Steward (1970) developed a homogeneous flow model for cryogenic chilldown.
The model treated the cryogenic fluid as a homogeneous mixture. The continuity,
momentum and energy equations of the mixture were solved to obtain density, pressure
and temperature of mixture. Various heat transfer regimes were considered: film boiling,
nucleate boiling, and single-phase convection heat transfer. Careful treatment of different
heat transfer regimes resulted in a significant improvement in the prediction of the
chilldown time. The homogeneous mixture model was also employed by Cross et al.
(2002) who obtained a correlation for the wall temperature during chilldown with an
oversimplified treatment of the heat transfer between the wall and the fluid.
Similar efforts have been devoted to the study of the re-wetting problem, referred to
as cooling down of a hot object. Thompson (1972) analyzed the re-wetting of a hot dry
rod. The two-dimensional temperature profile inside the solid rod was numerically
calculated. The nucleate boiling heat transfer coefficient between the solid rod and the
liquid was simplified to a power law relation and the heat transfer in the film boiling
stage is neglected. The liquid temperature and velocity outside the rod are assumed to be
constant. Sun et al. (1974), and Tien and Yao (1975) solved similar problems and
obtained an analytical solution for the re-wetting. They considered different heat transfer
coefficients for flow boiling and single-phase convection in order to obtain more accurate
results for re-wetting problems. In those works the thermal field of the liquid is neglected
and the heat transfer coefficients at the boiling and the convection heat transfer stage are
over-simplified, and the results are only valid for the vertical outer surface of a rod or a
tube.
Chilldown in stratified flow regime, which is the prevalent in the horizontal
pipeline, was first studied by Chan and Banerjee (1981 a, b, c). They developed a
comprehensive separated flow model for the cool-down in a hot horizontal pipe. Both
phases were modeled with one-dimensional mass and momentum conservation equations.
The vapor and liquid phase mass and momentum equations were reduced to two wave
equations for the liquid depth and the velocity of the liquid. The energy equation for the
liquid was used to find the liquid temperature and energy equation of vapor phase was
neglected. The wall temperature was computed using a 2-dimensional transient heat
conduction equation and heat transfer in the radial direction was neglected. They also
tried to evaluate the position of onset of re-wetting by studying the instability of film
boiling. Their prediction for the wall temperature agreed well with their experimental
results. Although significant progress was made in handling the momentum equations,
the heat transfer correlations employed were not as advanced.
Following Chan and Banerjee's (1981 a, b, c) separated flow model, Hedayatpour
et al. (1993) studied the cool-down in a vertical pipe with a modified separated flow
model. The flow regime is inverted annular film boiling flow, where the liquid core is
inside and the vapor film separates the cold liquid and the hot wall. This regime
frequently exists in cool-down in a vertical pipe. The modified separated flow model
retains the transient terms in the vapor momentum equation and the vapor phase energy
equation. The procedure is the following: first, the liquid mass conservation equation is
solved to obtain the liquid and vapor volume fractions. Then the vapor mass conservation
equation is used to solve the vapor velocity. The vapor momentum equation is
subsequently solved to obtain the vapor pressure. Finally, the liquid momentum equation
is employed to find the liquid velocity. The iteration stops when the solution is
converged. Although Chan and Banerjee (1981 a, b, c) and Hedayatpour et al. (1993)
were successful in the simulation of chilldown with the separated flow model, their
separated flow model is either incomplete or computationally inefficient.
c) Issues related to two-fluid model
The separated flow model is also called the two-fluid model, which consists of two
sets of conservation equations for the mass, momentum and energy of liquid and gas
phases. It was proposed by Wallis (1969), and further refined by Ishii (1975). Although
the two-fluid model is recognized as a useful computational model to simulate the
stratified multiphase flow in the pipeline, its application to the study of heat transfer in
two-phase flow in the pipeline is still limited.
The numerical scheme for the two-fluid model can be classified into two
categories. One is the compressible two-fluid model, which can be solved by a hyperbolic
equation solver. Examples are the commercial code OLGA (Bendikson et al., 1991),
Pipeline Analysis Code (PLAC) (Black et al., 1990) and Lyczkowski et al. (1978). The
other is the incompressible two-fluid model. Since the hyperbolic equation solver is not
applicable to incompressible two-fluid model, several approaches for incompressible
two-fluid model have emerged. One approach is to reduce the gas and liquid mass and
momentum equations to two wave equations for the liquid depth and velocity, such as in
Barnea and Taitel (1994b) and Chan and Banerjee (1991b). This treatment changed the
properties of two-fluid model. Hedayatpour et al. (1993) approach to two-fluid model is
not widely used due to lack of theoretical analysis on the convergence. Another approach
is to use the pressure correction method, which was initially introduced by Issa and
Woodburn (1998) and Issa and Kempf (2003) for the compressible two-fluid model.
Although their pressure correction scheme is powerful for simulating the multiphase flow
in the pipeline, the accuracy of the scheme is not reported. At the present, application of
pressure correction scheme on the multiphase flow with heat transfer in pipeline, such as
chilldown, does not exist.
d) Heat transfer in chilldown
A typical chilldown process involves several heat transfer regimes as shown in
Figurel-2. Near the liquid front is the film boiling regime. The knowledge of the heat
transfer in the film boiling regime is relatively limited, because i) film boiling has not
been the central interest in industrial applications; and ii) high temperature difference
causes difficulties in experimental investigations. For the film boiling on vertical
surfaces, early work was reported by Bromley (1950), Dougall and Rohsenow (1963) and
Laverty and Rohsenow (1967). Film boiling in a horizontal cylinder was first studied by
Bromley (1950); and the Bromley correlation was widely used. Breen and Westwater
(1962) modified Bromley's equation to account for very small tubes and large tubes. If
the tube is larger than the wavelength associated with Taylor instability, the heat transfer
correlation is reduced to Berenson's correlation (1961) for a horizontal surface.
Wall Liquid Front
---------------- r --- ------------ ;----r------------ -----
Cryogenic Lcuid
S/Valor
Convective Nucleation Film boiling
heat transfer boiling
Figure 1-2. The schematic of chilldown and heat transfer regime.
Empirical correlations for cryogenic film boiling were proposed by Hendrick et al.
(1961, 1966), Ellerbrock et al. (1962), von Glahn (1964), Giarratano and Smith (1965).
These correlations relate a simple or modified Nusselt number ratio to the Martinelli
parameter. Giarratano and Smith (1965) gave detailed assessment of these correlations.
All these correlations are for steady state cryogenic film boiling. Their suitability for
transient chilldown applications is questionable.
When the pipe wall chills down further, film boiling ceases and nucleate boiling
occurs. It is usually assumed that the boiling switches from film boiling to nucleate
boiling right away instead of passing through a transition boiling regime. The position of
the film boiling transitioning to the nucleate boiling is often called re-wetting front,
because from that position the cold liquid starts touching the pipe wall. Usually the
Leidenfrost temperature indicates the transition from film boiling to nucleate boiling.
However, the Leidenfrost temperature is not steady, and varies under different flow and
thermal conditions (Bell, 1967). A recent approach is to check the instability of the vapor
film beneath the liquid core using Kelvin Helmholtz instability analysis (Chan and
Banerjee, 1981c).
Studies on forced convection boiling are extensive (Giarratano and Smith, 1965;
Chen, 1966; Bennett and Chen, 1980; Stephan and Auracher, 1981; Gungor and
Winterton, 1996; Zurcher et al., 2002). A general correlation for saturated boiling was
introduced by Chen (1966). Gungor and Winterton (1996) modified Chen's correlation
and extended it to subcooled boiling. Enhancement and suppression factors for
macro-convective heat transfer were introduced. Gunger and Winterton's correlation can
fit experimental data better than the modified Chen's correlation (Bennett and Chen,
1980) and Stephan and Auracher correlation (1981). Recently, Zurcher et al. (2002)
proposed a flow pattern dependent flow boiling heat transfer correlation. This approach
improves the overall accuracy of heat transfer correlation by incorporating flow pattern.
Kutateladze (1952) and Steiner (1986) also provided correlations for cryogenic fluids in
pool boiling and forced convection boiling. Although they are not widely used, they are
expected to be more applicable for cryogenic fluids since the correlation was directly
obtained from cryogenic conditions. As the wall temperature drops further, boiling is
suppressed and the heat transfer is governed by two-phase convection; this is much easier
to deal with.
1.3 Scope
This dissertation focuses on understanding the unsteady fluid dynamics and heat
transfer of cryogenic fluids in a pipeline that is exposed to the atmospheric condition and
the corresponding solid heat transfer in the pipeline wall. Proper models for chilldown
simulation are developed to predict the flow fields, thermal fields, and residence time
during chilldown.
In Chapter 2, visualized experimental studies on heat transfer regimes and flow
regimes in cryogenic chilldown are reviewed. Based on the experimental observation,
homogeneous and separated flow models for the respectively vertical pipe and horizontal
pipe are discussed. The heat transfer models for the film boiling, flow boiling and forced
convection heat transfer in chilldown are reviewed and qualitatively assessed.
In Chapter 3, a physical model for vapor bubble growth in saturated nucleate
boiling is developed that includes both heat transfer through the liquid microlayer and
that from the bulk superheated liquid surrounding the bubble. Both asymptotic and
numerical solutions reveal the existence of a thin unsteady thermal boundary layer
adjacent to the bubble dome.
In Chapter 4, a pressure correction algorithm for two-fluid model is developed and
carefully implemented. Numerical stability of various convection schemes for both the
inviscid and viscous two-fluid model is analyzed. The connections between ill-posedness
of the two-fluid model and the numerical stability and between the viscous
Kelvin-Helmholtz instability and numerical stability are elucidated. The computational
accuracy of the numerical schemes is assessed.
In Chapter 5, a new film boiling coefficient is developed to accurately predict film
boiling heat flux for flow inside a pipe. The film boiling coefficient with the other
investigated heat transfer models are applied in building chilldown models. A
pseudo-steady chilldown model is developed to predict the chilldown time and the wall
temperature variation in a horizontal pipe in a reference frame that moves with the liquid
wave front. It is of low computational cost and allows for simple validation of the new
film boiling heat transfer correlations. A more comprehensive separated flow chilldown
model for the horizontal pipe is developed to predict the flow field of the liquid and vapor
and the temperature fields of the liquid, vapor and the solid wall in a fixed region of the
pipe flow. The unsteady development of the chilldown process for the vapor volume
fraction, velocities of the two-phases, and the temperatures of fluids and wall are
elucidated.
Chapter 6 concludes the research with a summary of the overall work and
discussion of the future works.
CHAPTER 2
TWO-PHASE FLOW MODELING AND FLOW BOILING HEAT TRANSFER OF
CRYOGENIC FLUID
Information of heat transfer regimes and flow regimes in cryogenic chilldown
obtained from the experimental study provides the foundation for modeling the heat
transfer and multiphase flow in chilldown. Based on the information of flow regimes,
corresponding flow models for simulating chilldown are discussed. For chilldown in the
vertical pipe, homogeneous flow model is preferred due to the prevalence of
homogeneous flow. For chilldown in the horizontal pipe, because the stratified flow is
prevalent, two-fluid model is adopted. The heat transfer models for film boiling, flow
boiling and forced convection heat transfer in chilldown are reviewed and qualitatively
assessed.
2.1 Flow Regime and Heat Transfer Regime
In the study of chilldown, one of the most important aspects of analysis is to
determine the type of flow regime in the given region of the pipe. The flow in cryogenic
chilldown is typically a two-phase flow, because liquid evaporates after a significant
amount of heat is transferred from the wall to the fluid during chilldown. The two-phase
flow regime is determined by many factors, such as fluid velocity, fluid density, vapor
quality, gravity, and pipe size. For horizontal flow, the flow regime is visually classified
as bubbly flow, plug flow, stratified flow, wavy flow, slug flow, and annular flow, as
shown in Figure 2-1. For vertical flow, the flow regimes include bubbly flow, slug flow,
churn flow, annular flow, as shown in Figure 2-2.
r_ -. =
..- -1
'- -
!7
Bubbly Flow
Plug Flow
Stratified Flow
WavyFlow
Slug/Internittent Flow
Annular Flow
Figure 2-1. Schematic of two-phase flow regime in horizontal pipe.
Bubble
Flow
Slug
Flow
Churn
Flow
Figure 2-2. Schematic of two-phase flow regime in vertical pipe.
Annual
Flow
The cryogenic two-phase flow is characterized by low viscosity, small density ratio
of the liquid to the vapor, low latent heat of vaporization, and large wall superheat. For
example, the liquid viscosity, density ratio latent heat of saturated liquid to vapor, and
latent heat of the saturated water at 1 atm are 2.73E-4 Pa*s, 1610, 2256.8kJ/kg,
respectively, while the corresponding data for saturated hydrogen at 1 atm arel.36E-
5Pa*s, 37.9, 444kJ/kg. Furthermore, film boiling, which is prevalent during chilldown,
causes low wall friction. These factors combined with the complex interaction between
the momentum and the thermal transportation make the two-phase flow during the
chilldown to distinguish itself from ordinary two-phase flows.
In the visualized horizontal chilldown experiment by Velat et al. (2004), as shown
in Figure 2-3, the pressure in the liquid nitrogen Dewar drives the fluid. When the liquid
nitrogen first enters the test section, a film boiling front is positioned at the inlet of test
section. This film boiling front produces a significant evaporation accompanied by a high
velocity vapor front traversing down the test section. If the mixture velocity is high
enough due to the large pressure drop between the Dewar and the outlet of the test
section, a very fine mist of liquid is entrained in the vapor flow. Immediately behind the
film boiling front is a liquid layer attached to the wall. The flow regime is either the
stratified flow or annular flow, depending on the flow speed, the pipe size, and the fluid
properties. If the mixture velocity is high, the flow likely appears as annular flow,
otherwise stratified flow or wavy flow is more common. The visual observation shows
that the liquid droplets being entrained in stratified flow and wavy flow is insignificant.
The nucleate boiling front follows the film boiling, indicating the end of film
boiling and the cryogenic liquid starts contacting the wall. The position where the liquid
starts contacting the wall is affected by the wall super heat, the liquid layer velocity and
the thickness of the liquid layer. It is a complex hydrodynamic and heat transfer
phenomenon. Usually Leidenfrost temperature indicates the transition from film boiling
to nucleate boiling. If the wall temperature is lower than the Leidenfrost temperature, the
vapor film cannot sustain the weight of liquid layer and becomes unstable. Therefore, the
liquid starts contacting the wall, and film boiling ceases.
Once the liquid contacts the wall, the nucleate boiling starts. In the nucleate boiling
regime the heat transfer from the wall to the liquid is significantly larger than that in the
film boiling regime, and the wall is chilled down much faster, are shown in Figure 2-4. If
the nucleation sites are not completely suppressed, a region of rapid nucleate boiling is
seen at the quenching front. If most of nucleate sites are suppressed by the subcooled
liquid, the flow directly transforms to the forced convection heat transfer, and nucleate
boiling stage is not visible.
After the nucleate boiling stage, the chilldown process dramatically slows down as
the convection heat transfer dominates. The wall superheat is relatively low at this stage
but the heat leaking from the test section to the environment emerges. These factors lead
to a lower chilldown rate. In the meantime, the liquid gradually builds up in the pipe due
to less vapor generation and the friction between the liquid and the wall. The increase of
the liquid layer thickness eventually leads to the transition of the flow regimes. When the
liquid layer is thick enough, the stratified flow or wavy flow becomes unstable.
Eventually slugs are formed and the flow transforms to the slug flow. In the final stage of
chilldown, the flow is almost a single-phase liquid flow, occasionally with some small
17
slugs. In this stage, the chilldown is almost completed, and the pipe wall temperature
gradually reaches the liquid saturated temperature.
Cryogenic
Liquid I Vapor Flow Increasing
Time
\Film Boiling
Front
Cryogenic Vapor Flow
Liquid
Liquid Film Flow \ Film Boiling
Region
Nucleate Boiling
Front
Cryogenic o oo o C o
Liquid Bubbly Flow Liquid Film Flow--
Figure 2-3. Schematics of observed flow structures in chilldown (Velat et al., 2004).
50
Temp 1 (T)
0 ---- Temp2(T)
-- T mp3(T)
SiTr~i. --,Eli
-50
E -100
H---
-1 50
-200
0 20 40 60 80 100 120 140 160 180 200
Time (sec)
Figure 2-4. Typical wall temperature variation during chilldown. (Velat et al., 2004)
Chilldown in a vertical pipe is practically less important than the chilldown in
horizontal pipe, due to the fact that most of cryogenic transportation pipelines are
horizontal, and only a small part is vertical. The experimental study (Hedayapour et al.
1993; Laverty and Rohsenow, 1967) reveals that the flow regime is mainly bubble flow,
or inverted annular flow if the vapor film of the film boiling is stable, and single-phase
vapor flow and single-phase liquid flow exist at the beginning and the final stage of
chilldown, respectively.
2.2 Flow Models in Cryogenic Chilldown
Based on the experimental investigation, several flow regimes exist in cryogenic
chilldown. At different flow regimes, the models for evaluating velocity and volume
fraction of fluid are different. Two types of flow models are to be discussed in this
section. First is the homogeneous flow model, which is used for modeling the chilldown
in a vertical pipe, where the homogeneous flow is prevalent. Another model is the two-
fluid model, which is mostly used in simulating the stratified flow or wavy flow for the
chilldown in a horizontal pipe.
2.2.1 Homogeneous Flow Model
In the homogeneous flow model, the unsteady mass, momentum, and energy
conservation equations for the mixture are simultaneously solved. The primary
assumptions are: (1) single-phase fluid or two-phase mixture is homogeneous, and each
phase is incompressible; (2) thermal and mechanical equilibrium exists between the
liquid and the vapor flowing together; (3) flow is quasi-one-dimensional; and (4) axial
diffusion of momentum and energy is negligible.
Thus, the continuity equation for the mixture is
S(PA) (P-A)
( + =0, (2.1)
at az
where p is the mixture density of liquid and vapor phase, i is the average fluid
velocity (by the assumption of homogeneous model, both liquid and vapor velocity are
u ), t is time z is the vertical axial coordinate, and A is the cross section area of the
pipeline.
Pipe wall
Vapor bubble
O
O 0 -- Liquid
00
00
-- Mixture front
Figure 2-5. Schematic for homogeneous flow model.
By neglecting the viscous terms, the momentum equation for the mixture becomes
S(PuA) a+ ( A) p + -p A A pgsin (2.2)
at fz dz \z )
where is pressure, is the pressure drop due to wall friction, f is the inclination
( az )f
angle of the pipe. For a vertical pipe, 3 =
2
The energy equation for the homogeneous model is
a(piA) a(puiA)
+ = q"S, (2.3)
at az
where i is the mixture enthalpy, q" is the heat flux from the wall to the fluid, and S is the
perimeter of the pipe.
If the cross section of the circular pipe is constant, the governing equations for
homogenous flow are simplified to the following equations.
a(p) a(p_)
+ = 0, (2.4)
at az
) (P ) -p ( pgsinO, (25)
at az az az
a(Ai) a(puz) 4
+ = q, (2.6)
at az IA
The pressure drop due to the wall friction is evaluated by the correlation for
The pressure dp z
the homogeneous system (Hewitt, 1982). In the correlation (Hewitt, 1982), a friction
multiplier 2o is defined as ratio of two-phase frictional pressure gradient to the
az
frictional pressure gradient for a single-phase flow at the same total mass flux and with
the physical properties of the liquid phase i.e.
w z lo
(dP
d f= 01 ,2 (2.7)
dz
where the friction multiplier 0 can be calculated by
r >-025
P0 Pg 1 -9Pg -0 25
12o 1+x l +x (2.8)
Pg t fg
where subscribes I and g represent the liquid phase and gas phase, respectively. The
single-phase pressure drop -z is evaluated using the standard equation
S)lo
=P 2 (2.9)
dz 1o Dp
where flo is the friction factor and for turbulent flow in a pipe, it is given as
025
fo = 0.079 (2.10)
in which, G is the mixture mass flux.
Compared with the experimentally measured two-phase flow pressure drop, the
homogenous model tends to underestimate the value of two-phase frictional pressure
gradient (Klausner et al., 1990). However, it provides a reasonable lower bound of the
two-phase flow pressure drop.
In the film boiling regime, a layer of vapor film separates the liquid core from the
pipe wall. This vapor film significantly reduces the wall friction, so that the two-phase
flow pressure drop due to the friction is much lower than that in the other heat transfer
regimes. To date, no correlation for the friction coefficient in the film boiling regime
exists. In available chilldown studies, the vapor film is treated as a part of the mixture and
Martinelli type of pressure drop correlation is used, or the wall friction is simply set to
zero.
2.2.2 Two-Fluid Model
In the chilldown inside the horizontal pipe, it is assumed that flow is stratified and
the liquid and the vapor flow at different velocity (Figure 2-6). Two-fluid model (Willis,
1969; Ishii, 1975) is widely used to qualitatively investigate the stratified flow inside
horizontal pipeline with a relatively low computational cost compared with
2-dimensional or 3-dimensional fluid flow models. In the study of the horizontal pipe
chilldown, the fluid volume fractions, velocities, enthalpies are solved with the two-fluid
model.
Pipe wall
S Vapor layer --
1D x
Liquid layer -
^r
Wall heat flux
Figure 2-6. Schematic of the two-fluid model.
The basis of the two-fluid model is a set of one-dimensional conservation equations
for the balance of mass, momentum and energy for each phase. The one-dimensional
conservation equations are obtained by integrating the flow properties over the
cross-sectional area of the flow.
In this study, it is assumed that flow is incompressible as the Mach number of the
gas phase is usually very low for the stratified flow. Hence, continuity equation for the
liquid phase (Chan and Banarjee, 1981c) is
a-(a, )+ (uza)=- (2.11)
at ax APz
where a is volume fraction, p is density, u is the velocity, t is the time, x is the axial
coordinate, and in' is the mass transfer rate between the liquid phase and the gas phase
per unit length; the subscript I denotes liquid.
Similarly, continuity equation for the gas phase is
a(g)+ (ugg= (2.12)
at ax Ap,
where the subscript g denotes gas. It is noted that
a, + a, = 1. (2.13)
The momentum equation for the liquid phase is
(ua,)+ a (u a, a, a=
at ax P, ax (2.14)
Ba, TS, TS mu
gcos OH, a -ag sinO- +
ax AP, AP, Ap,
where p, is the pressure at the liquid-gas interface, g is acceleration of gravity, f is the
angle of inclination of the pipe axis from the horizontal lane, 'is the shear stress, S is the
perimeter over which Tacts, A is the pipe cross section area, H, is the liquid phase
hydraulic depth; the subscript i denotes liquid-gas interface. The second term on the right
hand side of Equation (2.14) represents the effect of gravity on the wavy surface of liquid
layer. The liquid phase hydraulic depth H, is defined as
H1 =- _a (2.15)
aa, lh a'
where h, is the liquid layer depth.
Similarly, the momentum equation for the gas phase is
a (u' a ()+ I (U2 a. ag ap,
a ax P ax (2.16)
SBaa TS TS mi'U
g cos OHg 'a g sin- +-
ax Ap Ap, Ap,
whereHg is the gas phase hydraulic depth. It is defined as
ag ag
H, a (2.17)
aa lahg a,
where hg is the gas layer thickness.
To study heat transfer, appropriate energy equations for both phases are required in
the two-fluid model. Similar to the assumptions made in the homogeneous flow model,
the heat conduction inside the fluid is neglected. Thus the one-dimensional energy
equations for the liquid phase and the gas phase are
3 3 mi) 'i q,
(a + (a, + (2.18)
at ax Ap, Ap,
and
a 1fm qg
(jgag)+-(agUi)= + (2.19)
at ax Apg Apg
where i is enthalpy, and q' is the heat transfer rate to the fluid per unit length.
In the two-fluid model, shear stresses r,, rg and r, must be specified to close the
two fluid model. There are many correlations for shear stresses for separated flow model,
such as those developed by Wallis (1946), Barnea and Taitel (1976), and Andritsos and
Hanratty (1987). No significant difference exists among these models except at the flow
regime transition and at the high-speed flow, which will not be addressed in this study.
Thus, widely accepted shear stress correlations by Barnea and Taitel (1994) are
employed:
2
z, = f, (2.20)
2
Zg = fg Pg2 (2.21)
(U -U, IU g U,
,, = (2.22)
2
where ris shear stress, subscripts 1, g, and i represent interface between the liquid and the
wall, interface between the gas and the wall, interface between the liquid and gas,
respectively. Friction factors fare given by
f, = C1 Re -, and fg = C, Re,", (2.23)
where Re, is defined as
Re UID, (2.24)
fli
where D1 is the liquid hydraulic diameter
4A1
D, = (2.25)
SI
in which A, is liquid phase cross section area, S, is the liquid phase perimeter. In
Equation (2.23) Reg is defined as
PgUDg
Re = ggg (2.26)
flg
where Dg is vapor phase hydraulic diameter
4A
D = g (2.27)
S +S,
in which Ag is vapor phase cross section area, S, is the vapor phase perimeter, and S, is
the liquid-gas interface perimeter.
The coefficients C, and C, are equal to 0.046 for turbulent flow and 16 for
laminar flow, while n and m take the values of 0.2 for turbulent flow and 1.0 for laminar
flow. The interfacial friction factor is assumed to be f = fg or f = 0.014, if
fg < 0.014.
It is supposed that this model works in the flow boiling regime and in the forced
convection heat transfer regime. However, in the film boiling stage, presence of vapor
film dramatically reduces the shear stress between the liquid and the wall. In such a
situation, rT should be evaluated to include the effect of vapor film layer.
2.3 Heat Transfer between Cryogenic Fluid and Solid Pipe Wall
During cryogenic chilldown, the fluid in contact with the pipe wall is either the
liquid or the vapor. The mechanisms of heat transfers between the liquid and the wall and
between the vapor and the wall are different, as shown in Figure 2-7. Based on
experimental measurements and theoretical analysis, liquid-solid heat transfer accounts
for a majority of the total heat transfer. However, the liquid-solid heat transfer is much
more complicated than the heat transfer between the vapor and the wall due to occurrence
of film boiling and nucleate boiling. Thus, the heat transfer between the liquid and the
wall is discussed first.
2.3.1 Heat Transfer between Liquid and Solid wall
The heat transfer mechanism between the liquid and the solid wall includes film
boiling, nucleate boiling, and two-phase convection heat transfer. The transition from one
type of heat transfer to another depends on many parameters, such as the wall
temperature, the wall heat flux, and properties of the fluid. For simplicity, a fixed
temperature approach is adopted to determine the transition point. That is, if the wall
temperature is higher than the Leidenfrost temperature, film boiling is assumed. If the
wall temperature is between the Leidenfrost temperature and a transition temperature, T2,
nucleate boiling is assumed. If the wall temperature is below the transition temperature
T2, two-phase convection heat transfer is assumed. The values of the Leidenfrost
temperature and the transition temperature are determined by matching the model
prediction with the experimental results.
Pipe wall
Vapor layer -Convective heat Vapor
D transfer (vapor)
Liquid layer Liquid
Wall heat flux Thin vapor film
Convective heat Flow boiling Film boiling
transfer (liquid)
Figure 2-7. Schematic of heat transfer in chilldown.
2.3.1.1 Film boiling
Due to the high wall superheat encountered in the cryogenic chilldown, film boiling
plays a major role in the heat transfer process in terms of the time span and in terms of
the total amount of heat removed from the wall, as shown in Figure 2-4. Currently there
exists no specific film boiling correlation for chilldown applications with such high wall
superheat. The research starts from the conventional film boiling correlations.
A cryogenic film boiling heat transfer correlations was provided by Giarratano and
Smith (1965),
Nu 1Bo-04 f(tt), (2.28)
where Nu is Nusselt number
h D
Nu = (2.29)
ki
where hF is the film boiling heat transfer coefficient and k, is the thermal conductivity
of the liquid, Bo is the boiling number
Bo = q (2.30)
hfg *G
where hfg is the evaporative latent heat of the fluid. In Equation (2.28), Nucaic is the
Nusselt number for the two-phase convection heat transfer, which can be obtained using
Nucac = 0.023 ReO8*Pr 4, (2.31)
where Re is Reynolds number of mixture and Pr is Prandtl number of vapor, X, is
Martinelli number
X,,=1 x v I l (2.32)
x P,) ltJv
In Giarratano and Smith (1965) correlation, the heat transfer coefficient is the
averaged value for the whole cross section. Similar correlations for cryogenic film
boiling also exist in the literature. The correlations were obtained from measurements
conducted under steady state. The problem with the use of these steady state film boiling
correlations is that they do not account for information of flow regimes. For example, for
the same quality, the heat transfer rate for annular flow is much different from that for
stratified flow. Available empirical correlations do not make such difference.
Furthermore, in this study, local heat transfer coefficient is needed in order to
incorporate the thermal interaction with the pipe wall. Since the two-phase flow regime
information is available in the present study through the visualized experiment, it is
expected that the modeling effort should take into account the knowledge of the flow
regime. Suppose a liquid-gas stratified flow exists inside a horizontal pipe. Due to
gravity, the upper part of pipe wall is in contact with the gas, and lower part of pipe wall
is in contact with the flowing liquid. Thus, the heat transfer coefficient on upper wall is
significantly different from that on the lower wall. Apparently, the local heat transfer
coefficient strongly depends on the local flow condition instead an overall parameter such
as the flow quality at the given location.
There are several correlations for the film boiling based on the analysis of the vapor
film boundary layer, such as Bromley correlation (1950) and Breen and Westerwater
correlation (1962) for film boiling on the outer surface of a hot tube. Frederking and
Clark (1965) and Carey (1992) correlations, for the film boiling on the surface of a
sphere, are included as well. However, none of these was obtained for cryogenic fluids or
for the film boiling on the inner surface of a pipe or tube.
2.3.1.2 Forced convection boiling and two-phase convective heat transfer
A pool boiling correlation for cryogens was proposed by Kutateladze (1952). The
pool nucleated boiling heat transfer coefficient hpoo0 is
Fk 11 282 P1 750 (c 1 5 1
hP01= 0.487*10-10* (L kp) 0906 6 AT15, (2.33)
pool [hi p, 15 0906 10626
where o7is liquid surface tension, u is viscosity, and ATis wall superheat. Based on this
pool boiling correlation, a convection boiling correlation was proposed (Giarratano and
Smith, 1965). The heat transfer coefficient is contributed by both convection heat transfer
and ebullition:
h= hh,c + hpoo, (2.34)
where h,,c is given by Dittus-Boelter equation which is used in fully developed pipe
flow:
hI,~ = 0.023 Re,8 Pr,04 k,/D, (2.35)
where Re, is defined as
Re, = DG (2.36)
Chen (1966) introduced enhancement factor E and suppression factor S into the
flow boiling correlation. The heat transfer coefficient is given
h = EhlC + Shool (2.37)
Enhancement factor E reflects the much higher velocities and hence forced convection
heat transfer in the two-phase flow compared to the single-phase, liquid only flow. The
suppression factor S reflects the lower effective superheat in the forced convection as
opposed to pool boiling, due to the thinner boundary condition. The value of E and S are
presented as graphs in Chen (1966). The pool boiling heat transfer coefficient in Chen
correlation is
k079 045 049 025
hp 0.100122 : cO, P1 g AT 024AP075 (2.38)
pool 05 P 029hj024 024
Chen correlation (1966) fits best for annular flow since it was developed for vertical
flows. For the stratified flow regime, Chen's correlation may not be applicable.
At the flow boiling heat transfer, Gungor and Winterton correlation (1996) is
widely used due to that it fits much more experimental data. The basic form of Gungor
and Winterton correlation is similar to Chen correlation (1966), Equation (2.37).
However, evaluation of E and S in Gungor and Winterton's correlation takes account for
the influence of heat transfer rate by adding boiling number Bo. Thus, E and S are
presented as
E = 1+ 24000Bo 16 +1.37(1/,tX)86, (2.39)
and
1
S = 1 (2.40)
1+1.15x10-6E2 Re)17
The pool boiling correlation implemented is proposed by Cooper (1984)
hpoo1 = 55pM 12(-log p)-055M-605 q7 (2.41)
The solution of heat transfer correlation in Gungor and Winterton's correlation is
implicitly obtained by iteration.
Although Gungor and Winterton correlation (1996) is widely used due to its good
agreement with a large data set, a closer examination on this correlation shows that it is
based mainly on the following parameters: Pr, Re, and quality x. Similar to the
development of conventional film boiling correlations, these parameters all reflect overall
properties of the flow in the pipe and are not directly related to flow regimes. Thus, it
cannot be used to predict the local heat transfer coefficients required in chilldown
simulation.
Most of existing force convection boiling heat transfer correlations do not
effectively take account the influence of flow regimes and flow patterns. Recently,
Zurcher et al. (2002) proposed a flow pattern dependent heat transfer correlation for the
horizontal pipe. The strategy employed in Zurcher et al. (2002) is that the flow pattern is
obtained using the flow pattern map at the first step. The information of flow pattern
determines the part of wall contacting with the liquid or the vapor, then corresponding
conventional heat transfer correlations is employed to determine the local heat transfer
coefficient. The heat transfer coefficient for the whole pipe is obtained by averaging the
local heat transfer coefficient along the perimeter of the pipe. Although details of the
approach like flow pattern map, and correlations employed are not perfect in study of
Zurcher et al. (2002), their approach to the flow boiling heat transfer is intelligible and
provides insight for studying chilldown.
When wall superheat drops to a certain range all the nucleation sites are suppressed.
The heat transfer is dominated by two-phase forced convection. The heat transfer
coefficient can then be predicated using Equation (2.35), when the flow is turbulent, or
Equation (2.42), when the flow is laminar.
h,xc = 4.36 k /D,. (2.42)
2.3.2 Heat Transfer between Vapor and Solid Wall
The heat transfer between the vapor and wall can be estimated by treating the flow
as a fully developed forced convection flow, neglecting the liquid droplets that are
entrapped in the vapor. The heat transfer coefficient of vapor forced convective flow is
h = 0.023 *Reo8 Pr 04 k /Dg, (turbulent flow) (2.43)
hg, = 4.36 kg /Dg, laminarr flow) (2.44)
CHAPTER 3
VAPOR BUBBLE GROWTH IN SATURATED BOILING
Accurate evaluation of the nucleate boiling coefficient is a critical part of the study
on the chilldown process because it provides the heat transfer rate from the wall to the
cryogenic fluid. During the nucleate boiling the vapor bubble growth rate has a directly
influence on the heat transfer rate. The higher the bubble growth rate, the higher the heat
transfer rate. A physical model for vapor bubble growth in saturated nucleate boiling has
been developed that includes both heat transfer through the liquid microlayer and that
from the bulk superheated liquid surrounding the bubble. Both asymptotic and numerical
solutions for the liquid temperature field surrounding a hemispherical bubble reveal the
existence of a thin unsteady thermal boundary layer adjacent to the bubble dome. During
the early stages of bubble growth, heat transfer to the bubble dome through the unsteady
thermal boundary layer constitutes a substantial contribution to vapor bubble growth. The
model is used to elucidate recent experimental observations of bubble growth and heat
transfer on constant temperature microheaters reported by Yaddanapudi and Kim (2001)
and confirms that the heat transfer through the bubble dome can be a significant portion
of the overall energy supply for the bubble growth.
3.1 Introduction
During the past forty years, the microlayer model has been widely accepted and
used to explain bubble growth and the associated heat transfer in heterogeneous nucleate
boiling. The microlayer concept was introduced by Moore and Mesler (1961), Labunstov
(1963) and Cooper (1969). The microlayer is a thin liquid layer that resides beneath a
growing vapor bubble. Because the layer is quite thin, the temperature gradient and the
corresponding heat flux across the microlayer are high. The vapor generated by strong
evaporation through the liquid microlayer substantially supports the bubble growth.
Popular opinion concerning the microlayer model is that the majority of
evaporation takes place at the microlayer. A number of bubble growth models using
microlayer theory have been proposed based on this assumption such as van Stralen et al.
(1975), Cooper (1970), and Fyodrov and Klimenko (1989). These models were partially
successful in predicting the bubble growth under limited conditions but are not applicable
to a wide range of conditions. Lee and Nydahl (1989) used a finite difference method to
study bubble growth and heat transfer in the microlayer. However their model assumes a
constant wall temperature, which is not valid for heat flux controlled boiling since the
rapidly growing bubble draws a substantial amount of heat from the wall through the
microlayer, which reduces the local wall temperature. Mei et al. (1995a, 1995b)
considered the simultaneous energy transfer among the vapor bubble, liquid microlayer,
and solid heater in modeling bubble growth. For simplicity, the bulk liquid outside the
microlayer was assumed to be at the saturation temperature so that the vapor dome is at
thermal equilibrium with the surrounding bulk liquid. The temperature in the heater was
determined by solving the unsteady heat conduction equation. The predicted bubble
growth rates agreed very well with those measured over a wide range of experimental
conditions that were reported by numerous investigators. Empirical constants to account
for the bubble shape and microlayer angle were introduced.
Recently, Yaddanapudi and Kim (2001) experimentally studied single bubbles
growing on a constant temperature heater. The heater temperature was kept constant by
using electronic feed back loops, and the power required to maintain the temperature was
measured throughout the bubble growth period. Their results show that during the bubble
growth period, the heat flux from the wall through the microlayer is only about 54% of
the total heat required to sustain the measured growth rate. It poses a new challenge to
the microlayer theory since a substantial portion of the energy transferred to the bubble
cannot be accounted for.
Since a growing vapor bubble consists of a thin liquid microlayer, which is in
contact with the solid heater, and a vapor dome, which is in contact with the bulk liquid,
the experimental observations of Yaddanapudi and Kim (2001) leads us to postulate that
the heat transfer through the bubble dome may play an important role in the bubble
growth process, even for saturated boiling. Because the wall is superheated, a thermal
boundary layer exists between the background saturated bulk liquid and the wall; within
this thermal boundary layer the liquid temperature is superheated. During the initial stage
of the bubble growth, because the bubble is very small in size, it is completely immersed
within this superheated bulk liquid thermal boundary layer. As the vapor bubble grows
rapidly, a new unsteady thermal boundary layer develops between the saturated vapor
dome and the surrounding superheated liquid. The thickness of the new unsteady thermal
boundary layer should be inversely related to the bubble growth rate; see the asymptotic
analysis that follows. Hence the initial rapid growth of the bubble, which results in a thin
unsteady thermal boundary layer, is accompanied by a substantial amount of heat transfer
from the surrounding superheated liquid to the bubble through the vapor bubble dome.
This is an entirely different heat transfer mechanism than that associated with
conventional microlayer theory.
In fact, many previous bubble growth models have attempted to include the
evaporation through the bubble dome, such as Han and Griffith (1965) and van Stralen
(1967). However their analyses neglected the convection term in the bulk liquid due to
the bubble expansion, so the unsteady thermal boundary layer was not revealed. This
leads to a much lower heat flux through the bubble dome.
The existence and the analysis on the unsteady thermal boundary layer near the
vapor dome were first discussed in Chen (1995) and Chen et al. (1996), when they
studied the growth and collapse of vapor bubbles in subcooled boiling. For subcooled
boiling, the effect of heat transfer through the dome is much more pronounced due to the
larger temperature difference between the vapor and the bulk liquid. With the presence
of a superheated wall, a subcooled bulk liquid, and a thin unsteady thermal boundary
layer at the bubble dome, the folding of the liquid temperature contour near the bubble
surface was observed in their numerical solutions. The folding phenomenon was
experimentally confirmed by Mayinger (1996) using an interferometric method to
measure the liquid temperature.
Despite those findings, the existence of the thin unsteady thermal boundary layer
near the bubble surface has not received sufficient attention. In the recent computational
studies of bubble growth by Son et al. (1999) and Bai and Fujita (2000), the conservation
equations of mass, momentum, and energy were solved in the Eulerian or
Lagrange-Eulerian mixed grid system for the vapor-liquid two-phase flow. In their direct
numerical simulations of the bubble growth process, the heat transfer from the
surrounding liquid to the vapor dome is automatically included since the integration is
over the entire bubble surface. They observed that there could be a substantial amount of
heat transfer though bubble dome in comparison with that from the microlayer.
However, it is not clear that if these direct numerical simulations have sufficiently
resolved the thin unsteady thermal boundary layer that is attached to the rapidly growing
bubble.
In this study, asymptotic and numerical solutions to the unsteady thermal fields
around the vapor bubble are presented. The structure of the thin, unsteady thermal
boundary layer around the vapor bubble is elucidated using the asymptotic solution for a
rapidly growing bubble. A new computational model for predicting heterogeneous
bubble growth in saturated nucleate boiling is presented. The model accounts for energy
transfer from the solid heater through the liquid microlayer and from the bulk liquid
through the thin unsteady thermal boundary layer on the bubble dome. It is equally valid
for subcooled boiling, although the framework for this case has already been presented by
Chen (1995) and Chen et al. (1996). The temperature field in the heater is simultaneously
solved with the temperature in the bulk liquid. For the microlayer, an instantaneous linear
temperature profile is assumed between the vapor saturation temperature and the heater
surface temperature due to negligible heat capacity in the microlayer. For the bulk liquid,
the energy equation is solved in a body-fitted coordinate system that is attached to the
rapidly growing bubble with pertinent grid stretching near the bubble surface to provide
sufficient numerical resolution for the new unsteady thermal boundary layer. Section 3.2
presents a detailed formulation of the present model and an asymptotic analysis for the
unsteady thermal boundary layer. In Section 3.3, the experimental results of
Yaddanapudi and Kim (2001) are examined using the computational results based on the
present model. A parametric investigation considering the effect of the superheated bulk
liquid thermal boundary layer thickness on bubble growth is also presented.
3.2 Formulation
3.2.1 On the Vapor Bubble
Consideration is given to an isolated vapor bubble growing from a solid heating
surface into a large saturated liquid pool, as shown in Figure 3-1. A rigorous description
of the vapor bubble growth and the heat transfer processes among three phases requires a
complete account for the hydrodynamics around the rapidly growing bubble in addition
to the complex thermal energy transfer. The numerical analysis by Lee and Nydahl
(1989) relied on an assumed shape for the bubble, although the hydrodynamics based on
the assumed bubble shape is properly accounted for. Son et al. (1999) and Bai and Fujita
(2000) employed the Navier-Stokes equations and the interface capture or trace methods
to determine the bubble shape. Nevertheless, the microlayer structure was still assumed
based on existing models.
A\ Z
Bulk liquid
thermal boundary
SR(t) Background
bulk liquid
Liquid
microlayer Solid wall,heat is /
/ supplied from
within or below \
Figure 3-1. Sketch for the growing bubble, thermal boundary layer, microlayer and the
heater wall.
In this study, the liquid microlayer between the vapor bubble and the solid heating
surface is assumed to have a simple wedge shape with an angle 0<<1. The interferometry
measurements of Koffman and Plesset (1983) demonstrate that a wedge shape microlayer
is a good assumption. There exists ample experimental evidence by van Stralen (1975)
and Akiyama (1969) that as a bubble grows, the dome shape may be approximated as a
truncated sphere with radius R(t), as shown in Figure 3-1. Using cylindrical coordinates,
the local microlayer thickness is denoted by L(r). The radius of the wedge-shaped
interface is denoted by Rb (t), which is typically not equal to R(t). Let
c=Rb(t)/R(t), (3.1)
and the vapor bubble volume Vb (t) can be expressed as
4Ar
Vb(t)-= R3 (t)f(c), (3.2)
3
where f(c) depends on the geometry of the truncated sphere. In the limit c -> 1, the
bubble is a hemisphere and Vb (t) -> (2 r/3)R3 (t). In the limit c -> 0, the bubble
approaches a sphere and Vb (t) (4r / 3)R3(t).
To better focus the effort of the present study on understanding the complex
interaction of the thermal field around the vapor dome, additional simplification is
introduced. The bubble shape is assumed to be hemispherical (c=l) during the growth.
Comparing with the direct numerical simulation technique which solves bubble shape
and fluid velocity field using Navier-Stokes equation, this simplification introduces some
error in the bubble shape and fluid velocity and temperature fields in this study. However,
the hemispherical bubble assumption is generally valid at high Jacob number nucleate
boiling (Mei, et al. 1995a) and at the early stage of low Jacob number bubble growth
(Yaddanapudi and Kim, 2001). A more complete model that incorporates the bubble
shape variation could have been used, as in Mei et al. (1995a); however, the present
model allows for a great simplification in revealing and presenting the existence and the
effects of a thin unsteady liquid thermal boundary layer adjacent to the bubble dome and
the influence of bulk liquid thermal boundary layer on saturated nucleate boiling. The
present simplified model is not quantitatively valid when the shape of the vapor bubble
deviates significantly from a hemisphere.
The energy balance at the liquid-vapor interface for the growing bubble depicted in
Figure 3-1 is described as
dV aT, T,
pVhf d=- -kl Tmi dA+ f k- dAb (3.3)
Sa z=L(r) R'=R (t)
where p, is the vapor density, hfg is the latent heat, k, is the liquid thermal conductivity,
Ti is the temperature of the bulk liquid, Tm, is the temperature of the microlayer liquid, Am
is the area of wedge, Ab is the area of the vapor bubble dome exposed to bulk liquid,
an
is the differentiation along the outward normal at the interface, and R' is the spherical
coordinate in the radial direction attached to the moving bubble. Equation (3.3) simply
states that the energy conducted from the liquid to the bubble is used to vaporize the
surrounding liquid and thus expand the bubble.
3.2.2 Microlayer
The microlayer is assumed to be a wedge centered at r = 0 with local thickness
L(r). Because the hydrodynamics inside the microlayer are not considered, the
microlayer wedge angle 0 cannot be determined as part of the solution. In Cooper and
Lloyd (1969), the angle 0 was related to the viscous diffusion length of the liquid as
Rb (t)tan0 = c, vt in which v, is the kinematic viscosity of the liquid. A small 0
results in
= (3.4)
Rb (t)
Cooper and Lloyd (1969) estimated c, to be within 0.3-1.0 for their experimental
conditions.
A systematic investigation for saturated boiling by Mei et al. (1995b) established
that the temperature profile in the liquid microlayer can be taken as linear for practical
purposes. The following linear liquid temperature profile in the microlayer is thus
adopted in this study
Trrzt)=T+AT{rt -- ,l (3.5)
where ATat (r, t) = T1 (r, z = 0, t) Ts, and T, is the temperature of the solid heater.
3.2.3 Solid Heater
The temperature of the solid heater is governed by the energy equation, which is
coupled with the microlayer and bulk liquid energy equations. Solid heater temperature
variation significantly influences the heat flux into the rapidly growing bubble (Mei et al.
1995a, 1995b). However, in this study, constant wall temperature is assumed so that the
case of Yaddanapudi and Kim (2001) can be directly simulated. Thus,
AT,,t (r,t)= AT=,t = T Ts (3.6)
which can be directly used in Equation (3.5) to determine the microlayer temperature
profile.
3.2.4 On the Bulk Liquid
It was assumed that the vapor bubble is hemispherical in section 2.1. Furthermore,
the velocity and temperature fields are assumed axisymmetric. Unless otherwise
mentioned, spherical coordinates (R', y, (p), as shown in Figure 3-2, are employed for the
bulk liquid.
7/=1
Figure 3-2. Coordinate system for the background bulk liquid.
3.2.4.1 Velocity field
Since there is no strong mean flow over the bubble, the bulk liquid flow induced by
the growth of the bubble is mainly of inviscid nature. Thus the liquid velocity field may
be determined by solving the Laplace equation V2c = 0 for the velocity potential d. In
spherical coordinates, the velocity components are simply given by the expansion of the
hemispherical bubble as
=o yV
R l=1
Figure 3-2. Coordinate system for the background bulk liquid.
3.2.4.1 Velocity field
Since there is no strong mean flow over the bubble, the bulk liquid flow induced by
the growth of the bubble is mainly ofinviscid nature. Thus the liquid velocity field may
be determined by solving the Laplace equation V72 = 0 for the velocity potential 0. In
spherical coordinates, the velocity components are simply given by the expansion of the
hemispherical bubble as
uR, dR R(t) R = 0,u, u=O, (3.7)
dt R' R'
where R = dR(t)
dt
3.2.4.2 Temperature field
By assuming axisymmetry for the temperature fields and using the liquid velocity
from Equation (3.7), the unsteady energy equation for the bulk liquid in spherical
coordinates is
a T, a T (, 12 3 1 a TI
UR I a R s sin (3.8)
The boundary conditions are
a T,
T =0 at = 0, (3.9)
T, = T, at (3.10)
2
T, = Tso, at R' = R(t), (3.11)
T, = T7 (V, t) at R' --- (3.12)
where T is the far field temperature distribution.
To facilitate an accurate computation and obtain a better understanding on the
physics of the problem, the following dimensionless variables are introduced,
t- R' TI Tb
t= ,,R' ,z =T b (3.13)
t' R T -Tb
where t, is a characteristic time chosen to be the bubble departure time, T7 is the initial
solid temperature at the solid-liquid interface, and T, is the bulk liquid temperature far
away from the wall, which equals T,, for saturated boiling.
Using Equation (3.7) and Equation (3.13), Equation (3.8) can be written as
R J 1 R + a 1 k 1o
tR Jr kR2 )R' RR 2 R (3.14)R
(3.14)
a, 1 1 (sin ,
+ --= smw-
RR sinV R'2 (si
In this equation, the first term on the left-hand-side (LHS) is the unsteady term, and the
second term is due to convection in a coordinate system that is attached to the expanding
bubble. The right-hand-side (RHS) terms are due to thermal diffusion.
As shown in Chen et al. (1996) and below, the solution for 0O near R = 1 possesses
Ri
a thin boundary layer when >> 1. Therefore, to obtain the accurate heat transfer
a,
between the bubble and the bulk liquid, high resolution in the thin boundary layer is
essential. Hence, the following grid stretching in the bulk liquid region is applied,
R'=1+(R- -l){1-S, tan 1[(1- 7)tan(l/SR,)]} for 0 77 1
Sr(3.15)
V= S, tan-' tan(l/S)] for 0
where SR, and S, are parameters that determine the grid density distribution in the
R'
physical domain and R= is the far field end of the computational domain along the
R
radial direction.
Typically SR,-0.65 and S,-0.73, and R' ranges from 5 to 25. Figure 3-3 shows a
typical grid distribution used in this study.
Figure 3-3. A typical grid distribution for the bulk liquid thermal field with SR = 0.65,
S = 0.73, andR' =10.
3.2.4.3 Asymptotic analysis of the bulk liquid temperature field during early stages
of growth
To gain a clear understanding on the interaction of the growing bubble with the
background superheated bulk liquid thermal boundary layer, an asymptotic analysis for
non-dimensional temperature 0, is presented, following the work of Chen et al. (1996).
During the early stages, the bubble growth rate is high and expands rapidly so that
RR
A = -- >> 1. (3.16)
a,
Thus, the solution to Equation (3.14) includes an outer approximation in which the
thermal diffusion term on the RHS of Equation (3.14) is negligible and an inner
approximation (boundary layer solution) in which the thermal diffusion balances the
convection. Away from the bubble, the outer solution is governed by
R +oOUt = u 0, (3.17)
R at (R' ) R'
where Oo,"uis the outer solution for 0, in the bulk liquid. The general solution for
Equation (3.17) is
0"'t = F(R(t)(R'3 -1)3), (3.18)
as given in Chen et al. (1996). In the above F is an arbitrary function and it is determined
from the initial condition of 0, or the temperature profile in the background bulk liquid
thermal boundary layer. It is noted that the solution for0,outis described by
R(t)(R~3 -1)13 = const along the characteristic curve.
The initial temperature profile is often written as 00 = f : The solution of
Equation (3.17) is thus expressed as
OUt f Ro co 1+ -R' 1 1) (3.19)
where R0 is the initial bubble radius at t = to << tc. Provided the bubble growth rate is
high, i.e. A >> 1, Equation (3.19) is not only an accurate outer solution for the
temperature field outside a rapidly expanding bubble, but it is also a good approximation
for the far field boundary condition for Equation (3.12).
Near the bubble surface, there exists a large temperature gradient between the
saturation temperature on the bubble surface and the temperature of the surrounding
superheated liquid over a thin region. Therefore, the effect of heat conduction is no
longer negligible in this thin region and must be properly accounted for. For a large value
of A, a boundary layer coordinate Xis introduced,
R'-1
X (3.20)
8 (A)
where 8* (A) << 1 is the dimensionless length scale of the unsteady thermal boundary
layer. Substituting Equation (3.20) into Equation (3.14) results in
R 0a 0 R aX aOO ( + 81 -1 a ol^"-
tR tR t 3X ( (+8 x)2 8* ax
(3.21)
1 (20'n" 28' 0a" 1 1 1 ( 0a ^
+ + --- sin V/
AS*2 KX2 +8*X X A (l+8*X)2 sin/ fa/f K j
Neglecting higher order terms, Equation (3.21) becomes
R 0 _ei 1 20n R X
-R a 1 a2- I+3X 1- R aX a (3.22)
tR ar AS 2 X 2, 3tRX arT 3X
The balance between the convection term and the diffusion term on the RHS of Equation
(3.22) requires
5 = A- = (3.23)
Hence Equation (3.22) becomes
R a I a2oin+ R A lx
t = i+ 3 2 AX (3.24)
Te R cti f th i ) t X
The boundary conditions for the inner (boundary layer) solution are
Tsa-T
0"n = Osat b at X = 0, (3.25)
Tw Tb
0," = 1 at X -- (3.26)
1 R
For << 1, R(t) c t2 and -= Thus, the LHS of Equation (3.24) is small and can
tcR
also be neglected. Equation (3.24) then reduces to
S2 orn anor
/ +3X =0 for r <<1. (3.27)
aX2 IX
The solution for Equation (3.27) is
;" = (1 O,)erf(iX)+Osa,. (3.28)
For << 1, by matching the outer and inner solutions given by Equation (3.19) and
Equation (3.28), the uniformly valid asymptotic solution of the bulk liquid temperature
for the saturated boiling problem considered here is obtained,
S Rocos R 1 ( 3 -1 i -), (3.29)
T -T
where 0s,, = sa b and erfc is the complimentary error function. Equation (3.29) is an
T-T
T. Tb
asymptotic solution for 0, valid for << 1.
The asymptotic solution given by Equation (3.29) for the liquid thermal field
provides an analytical framework to understand: 1) how the temperature field of
background superheated bulk liquid boundary layer influences the temperature 0, near
the vapor bubble through the function f 2) how the bubble growth R(t) and liquid thermal
diffusivity affect the liquid thermal field 0, through the rescaled inner variable X as
defined in Equation (3.20) and Equation (3.23); and 3) how the folding of the temperature
contours near the bubble occurs through the dependence of cos y term in Equation
(3.29). More importantly, from a computational standpoint, it provides: 1) an accurate
measure on the thickness of the rapidly moving thermal boundary layer; and 2) a reliable
guideline for estimating the adequacy of computational resolution in order to obtain an
accurate assessment of heat transfer to the bubble.
3.2.5 Initial Conditions
The computation must start from a very small but nonzero initial time zT, so that
R(,r0) is sufficiently small at the initial stage. To obtain enough temporal resolution for
the initial rapid growth stage and to save computational effort for the later stage, the
following transformation is used,
T=o 2. (3.30)
Thus a constant "time step" Ao can be used in the computation.
The initial temperature profile inside the superheated bulk liquid thermal boundary
layer plays an important role to the solution of 0,, which in turn affects the heat transfer
to the bubble through the dome.
There exist both experimental and theoretical studies that have considered the bulk
liquid temperature profile in the vicinity of a vapor bubble. Hsu (1962) estimated the
temperature profile of the superheated thermal layer adjacent to the heater surface and
found the layer to be quite thin; thus the temperature gradient inside the thermal layer is
almost linear. However, beyond the superheated layer the temperature is held essentially
constant at the bulk temperature due to strong turbulent convection. The experimental
study by Wiebe and Judd (1971) revealed similar results. It was found that the
superheated bulk liquid thermal boundary layer thickness, 6, decreases with increasing
wall heat flux due to enhanced turbulent convection. A high wall heat flux results in
increased bubble generation, and the bulk liquid is stirred more rapidly by growing and
departing vapor bubbles. To estimate the superheated layer thickness, Hsu (1962) used a
thermal diffusion model within the bulk liquid. Han and Griffith (1965) used a similar
model and estimated the thickness to be
8= zat, (3.31)
where tw is the waiting period. The thermal diffusion model often overestimates the
thermal layer thickness, as it neglects the turbulent convection, which is quite strong as
reported by Hsu (1962) and Wiebe and Judd (1971).
Generally, the bulk liquid temperature profile is almost linear inside the
superheated thermal boundary layer, and remains essentially uniform at the bulk
temperature Tb beyond the superheated background thermal boundary layer.
Accordingly, the initial condition for the bulk liquid thermal field used in the numerical
solution is given by
S-- <6
,0 = z (3.32)
[o, z 8
In the asymptotic solution, the discontinuity of 0D, /Dz in the above profile causes
the solution for 0, to be discontinuous. For clarity, the following exponential profile is
employed in representing the asymptotic solution
0 = exp (3.33)
(58
3.2.6 Solution Procedure
An Euler backward scheme is used to solve Equation (3.14). A second order
upwind scheme is used for the convection term and a central difference scheme is used
for the thermal diffusion terms.
After the bulk liquid temperature field is obtained, the solid heater temperature
field is solved, and the bubble radius R(r) is updated using Equation (3.3) and Euler's
explicit scheme. The information for R(r) is a necessary input in Equation (3.14).
Although the solution for R(r) is only first order accurate in time, the O(Ar) accuracy is
not a concern here because a very small Ar has to be used to ensure sufficient resolution
during the early stages. Typically, n = 104 time steps are used.
3.3 Results and Discussions
3.3.1 Asymptotic Structure of Liquid Thermal Field
To gain an analytical understanding of the liquid thermal field near the bubble and
to validate the accuracy of the computational treatment for the thin unsteady thermal
boundary layer, comparison between the computational and the asymptotic solutions for
0, near the bubble surface is first presented. As mentioned previously, the validity of the
outer solution of the asymptotic analysis only requires A >> 1, which is satisfied under
most conditions due to rapid vapor bubble growth. The inner solution is valid for << 1
in addition to A >> 1.
The comparison is presented for bubble growth in saturated liquid with A=14000.
The initial temperature profile follows Equation (3.33) and = 0.5, in which Rc is the
R
bubble radius at t There are 200 and 50 grid intervals along the R' and y -directions,
respectively. The grid stretching factors are S = 0.65 and S, = 0.73 for the
computational case.
Figure 3-4 compares the temperature profiles between the asymptotic and
numerical solutions at z = 0.001, 0.01, 0.1, and 0.3 for V/= 0, 400, and 710. There are
two important points to be noted. First of all, it is seen that the temperature gradient is
indeed very large near the bubble surface because the unsteady thermal boundary layer is
very thin. Secondly, numerical solutions agree very well with the asymptotic solutions at
z = 0.001 and 0.01. The excellent agreement between the numerical and analytical
solutions indicates that the numerical treatment in this study is correct. At z = 0.1 and
0.3, the asymptotic inner solution given by Equation (3.29) is no longer accurate, while
the outer solution remains valid because A >> 1 is the only requirement. At z = 0.1 and
0.3 the numerical solution matches very well with the outer solution. This again
demonstrates the integrity of the present numerical solution over the entire domain due to
sufficient computational resolution near the bubble surface and removal of undesirable
numerical diffusion through the use of second order upwind scheme in the radial
direction for Equation (3.14).
The large temperature gradient near the dome causes high heat transfer from the
superheated liquid to the vapor bubble through the dome. This large gradient results from
the strong convection effect that is caused by the rapid bubble growth (see Equation
(3.14) for the origin in the governing equation and Equation (3.19) for the explicit
dependence on the bubble growth). Thus the bulk liquid in the superheated boundary
layer supplies a sufficient amount of energy to the bubble.
54
09 09 -
07 0 \\ 7 -
6 05- 6 5-
04 0
S=0001 T=0 01
0 3 0 (Numerical) =0 (Numerical)
S- V40 (Numenrcal) 0 (Asymptotic)
=71 (Numerical) 400(Numerical)
S =00 (Asymptotic) y=40 (Asymptotic)
= 40o (Asymptotic) 710 (Numerical)
v=71 (Asymptotic) -- 71 (Asymptotic)
10' 1 201 10 10? '3 101 io' 101
=01 =03
0- 00 (Numerical) \ n3 =00 (Numerical) \ \
y=400 (Numerical) y=400(Numerical) \
v710 (Numerical) v710 (Numerical)
yo00 (Asymptotic) 0FO-- =00(Asymptotic)
y=400 (Asymptotic) y=400 (Asymptotic)
710 (Asymptotic) ----- 710(Asymptotic)
Il I Il
R'-1 R'-I
Figure 3-4. Comparison of the asymptotic and the numerical solutions at r=0.001, 0.01,
0.1 and 0.3 for 100, 400, and 710
To capture the dynamics of the unsteady boundary layer, a sufficient number of
computational grids is required inside this layer. The asymptotic analysis gives an
estimate for the unsteady boundary layer thickness on the order of
6* ~ = 0.025, which agrees with the numerical solution in Figure 3-4. At
V14000
S= 0.01 in Figure 3-4, the discrete numerical results are presented. There are about 23
points inside the layer of thickness 6* = 0.025 This provides sufficient resolution for the
temperature profile in the unsteady thermal boundary layer. In contrast, most
computational studies on the thermal field around the bubble dome reported in the open
literature have insufficient grid resolution adjacent the dome, which leads to an
inaccurate heat transfer assessment.
Figure 3-5 shows the effect of parameter A on the asymptotic solution. When A is
large, the asymptotic and numerical solutions agree very well. The discrepancy between
asymptotic and numerical solutions inside the unsteady thermal boundary layer increases
when A decreases. However, the outer solution remains valid for the far field even when
A becomes small.
09 A=14000
08
10000
00
05 -
04 T=0.01
v=40'
03 I
03 Asymptotic
02 Numerical
01 -
0
103 102 10 100 10'
Figure 3-5. Effect of parameter A on the liquid temperature profile near bubble.
The temperature contours shown in Figure 3-6 are difficult to obtain
experimentally. Only recent progress in holographic thermography permits such
measurements. Ellion (1954) has stated that there exists an unsteady thermal boundary
layer contiguous to the vapor bubble during the bubble growth. Recently, Mayinger
(1996) used a holography technique to capture the folding of the temperature contours
during subcooled nucleate boiling. Although his study considered subcooled nucleate
boiling, the pattern of the temperature distribution near the bubble dome by Mayinger
(1996) is very similar to that shown in the Figure 3-6. It is expected that experimental
evidence of contour folding in saturated nucleate boiling will be reported in the future.
3.3.2 Constant Heater Temperature Bubble Growth in the Experiment of
Yaddanapudi and Kim
In the experiment of Yaddanapudi and Kim (2001), single bubbles growing on a
heater array kept at nominally constant temperature were studied. The liquid used is
FC-72, and the wall superheat is maintained at 22.5 OC, so that Jacob number is 39. The
bubble shape in the early stage appears to be hemispherical. To calculate the heat flux
from the microlayer to the vapor bubble in the present model, the microlayer wedge angle
or constant c, in Equation (3.4) must be determined. Neither 0 or c, has been
measured. However, the authors have reported the amount of wall heat flux from the wall
to the bubble through an equivalent bubble diameter deq assuming that the wall heat flux
is the only source of heat entering the bubble. Since in the present model this heat flux is
assumed to pass through the microlayer, it may be used to evaluate the constant c, via
trial and error. The superheated thermal boundary layer thickness Sof the bulk liquid in
Equation (3.32) is also a required input. The computed growth rate R(t) is matched with
the experimentally measured R(t) in order to determine The simulation is carried out
only for the early stage of bubble growth. This is because after t=6-8x10-4s the base of
the bubble does not expand anymore, and the bubble shape deviates from a hemisphere.
Furthermore, there is the possibility of the microlayer being dried out in the latter growth
stages as a result of maintaining a constant wall temperature, as was observed by Chen et
al. (2003).
S50.01 15 05 50.1
o 05 1 15 2 05 1 15 2
R'-1 R'-1
1 1 Ri
R'-I R'-I
Figure 3-6. The computed isotherms near a growing bubble in saturated liquid at T0.01,
'0.1, 0.3, and 0.9.
Figure 3-7 shows the computed equivalent bubble diameter d,, (t), together with
the experimentally determined equivalent dq (t). In the present model, dq, is calculated
using
p h*dt =- k dA,
P uh i fg dz=L(r)
(3.34)
fz 3
where Vb = deq The heat flux includes only that from the microlayer and this allows
6
c, to be evaluated. For d to match the measured data as shown in Figure 3-7, it
requires c,=3.0.
0 0005 -
000045 - deq(t) present model
A deq(t) measurement
00004
0 00035
00003
0 00025 --
00002
000015 -
0 0001
^ -
5E-05
0 I II II I I I
00 0002 00004 00006 00008
t (s)
Figure 3-7. Comparison of the equivalent bubble diameter dq for the experimental data
of Yaddanapudi and Kim (2001) and that computed for heat transfer through
the microlayer (c, =3.0).
Figure 3-8 compares the computed bubble diameter d(t) = 2R(t) and those
reported by Yaddanapudi and Kim (2001). In Figure 3-8, =30jum is used in addition to
c, =3.0 in matching the predicted bubble growth with measured data. The good
agreement obtained can be partly attributed to the adjustment in the superheated bulk
liquid thermal boundary layer thickness Because the heat transfer to the bubble
(through the microlayer and through the dome) is of two different mechanisms, the good
agreement over the range is an indication of the correct physical representation by the
present model.
Figure 3-9 shows the total heat entering bubble and the respective contribution
from the microlayer and from the unsteady thermal boundary layer. The contribution
from the unsteady thermal boundary layer accounts for about 70% of the total heat
transfer. It was reported by Yaddanapudi and Kim (2001) that approximately 54% of the
total heat is supplied by the microlayer over the entire growth cycle. Since, the simulation
is only carried out for the early stage of bubble growth, it is difficult to compare the
microlayer contribution to heat transfer reported by Yaddanapudi and Kim (2001) with
that predicted by current model. At the end, the bubble expands outside the superheated
boundary layer and protrudes into the saturated bulk liquid. The heat transfer from dome
thus slows down. Hence, the 54% for the entire bubble growth period dose not contradict
a higher percentage of contribution computed from the unsteady thermal boundary layer
during the early stages.
Figure 3-10 shows the computed temperature contours associated with
Yaddanapudi and Kim's (2001) experiment for the estimated Sand c Folding of the
temperature contours is clearly observed in the simulation for saturated boiling.
3.3.3 Effect of Bulk Liquid Thermal Boundary Layer Thickness on Bubble Growth
Since the superheated bulk liquid thermal boundary layer thickness, 6, determines
how much heat is stored in the layer, it is instructive to conduct a parametric study on the
effects bubble growth with varying All parameters are the same as those used in
Yaddanapudi and Kim's (2001) experiment except that 6 is varied. Hence the influence
of the superheated thermal boundary layer thickness Son the bubble growth is elucidated.
d(t) present model
* d(t) measurement
I i I I I I I I I
00002
00006
Figure 3-8. Comparison of bubble diameter, d(t), between that computed using the
present model and the measured data of Yaddanapudi and Kim (2001). Here,
c, =3.0 and 8=30pm.
1 6E-05 -
1 4E-05
1 2E-05
1E-05
8E-06
6E-06
4E-06
2E-06
total heat entering the bubble
heat from micorlayer
heat from bulk liquid thermal boudary layer
I I I
00004
00002
00006
00008
Figure 3-9. Comparison between heat transfer to the bubble through the vapor dome and
that through the microlayer.
00005
00035 -
00003
0 00025
00002
00015
00008
Figure 3-11 shows the effect on the bubble growth rate of varying S(from l1Jm to
100um). The thicker the bulk liquid thermal boundary layer, the faster the bubble grows.
A large Simplies a larger amount of heat is stored in the background bulk liquid
surrounding the bubble. It is also clear that when Japproaches zero, the bubble growth
rate becomes unaffected by the variation of The reason is when Sis small, most of heat
supplied for bubble growth comes from the microlayer and the contribution from the
dome can be neglected.
t5 0.0012ms 5 t=0.12ms o
5 5 15 5 2
RI'-1 R'-1
t=0.36ms t=0.96ms
I 15 2 o5 1 15 2
R'-I R'-
1 5 1_5
Figure 3-10. The computed isotherms in the bulk liquid corresponding to the thermal
conditions reported by Yaddanapudi and Kim (2001).
It is also noted that for =100/m, if the bubble eventually grows to about several
millimeters, the effect of the bulk liquid thermal boundary layer is negligible on R(t) for
most of the growth period except at the very early stages. Physically, this is because the
bubble dome is quickly exposed to the saturated bulk liquid so that it is at thermal
equilibrium with the surroundings. For small bubbles, it will be immersed inside the
thermal boundary layer most of time. Hence the effect of the bulk liquid thermal
boundary layer becomes significant for the bubble growth.
0 0007
8=100 pm
00006 -
00005 =50 m
00004 8=30 pm
00003 =10 pm
6=5 pm
8=1 Pm
00002
00001
0 00002 00004 00006 00008 0001
t(s)
Figure 3-11. Effect of bulk liquid thermal boundary layer thickness Son bubble growth.
The microlayer angle 0 and the superheated bulk liquid thermal liquid boundary
layer thickness Sare the required inputs to compute bubble growth in the present model.
However, neither of these parameters is typically measured or reported in bubble growth
experiments. It is strongly suggested that the bulk liquid thermal boundary layer
thickness 6be measured and reported in future experimental studies. For a single bubble
study, Sin the immediate neighborhood of the nucleation site should be measured.
3.4 Conclusions
In this study, a physical model is presented to predict the early stage bubble growth
in saturated heterogeneous nucleate boiling. The thermal interaction of the temperature
fields around the growing bubble and vapor bubble together with the microlayer heat
transfer is properly considered. The structure of the thin unsteady liquid thermal
boundary layer is revealed by the asymptotic and numerical solutions. The existence of a
thin unsteady thermal boundary layer near the rapidly growing bubble allows for a
significant amount of heat flux from the bulk liquid to the vapor bubble dome, which in
some cases can be larger than the heat transfer from the microlayer. The experimental
observation by Yaddanapudi and Kim (2001) on the insufficiency of heat transfer to the
bubble through the microlayer is elucidated. For thick superheated thermal boundary
layers in the bulk liquid, the heat transfer though the vapor bubble dome can contribute
substantially to the vapor bubble growth.
CHAPTER 4
ANALYSIS ON COMPUTATIONAL INSTABILITY IN SOLVING TWO-FLUID
MODEL
The two-fluid model is widely used in studying gas-liquid flow inside pipelines
because it can qualitatively predict the flow field with a low computational cost.
However, the two-fluid model becomes ill-posed when the slip velocity between the gas
and the liquid exceeds a critical value. Computationally, even before the flow becomes
unstable, computations can be quite unstable to render the numerical result unreliable. In
this study computational stability of various convection schemes for the two-fluid model
is analyzed. A pressure correction algorithm is carefully implemented to minimize its
effect on stability. Von Neumann stability analysis for the wave growth rates by using the
1st order upwind, 2nd order upwind, QUICK (quadratic upstream interpolation for
convection kinematics), and the central difference schemes are conducted. For inviscid
two-fluid model, the central difference scheme is more accurate and more stable than
other schemes. The 2nd order upwind scheme is much more susceptible to instability for
long waves than the 1st order upwind and inaccurate for short waves. The instability
associated with ill-posedness of the two-fluid model is significantly different from the
instability of the discretized two-fluid models. Excellent agreement is obtained between
the computed and predicted wave growth rates, when various convection schemes are
implemented.
The pressure correction algorithm for inviscid two-fluid model is further extended
to the viscous two-fluid model. For a viscous two-fluid model, the diffusive viscous
effect is modeled as a body force resulting from the wall friction. Von Neumann stability
analysis is carried out to assess the performances of different discretization schemes for
the viscous two-fluid model. The central difference scheme performs best among the
schemes tested. Despite its nominal 2nd order accuracy, the 2nd order upwind scheme is
much more inaccurate than the 1st order upwind scheme for solving viscous two-fluid
model. Numerical instability is largely the property of the discretized viscous two-fluid
model but is strongly influenced by VKH instability. Excellent agreement between the
computed results and the predictions from von Neumann stability analysis for different
numerical scheme is obtained. Inlet disturbance growth test shows that the pressure
correction scheme is capable to correctly handle the viscous two-phase flow in a pipe.
4.1 Inviscid Two-Fluid Model
4.1.1 Introduction
Gas-liquid flow inside a pipeline is prevalent in the handling and transportation of
fluids. A reliable flow model is essential to the prediction of the flow field inside the
pipeline. To fully simulate the system, Navier-Stokes equations in three-dimensions are
required. However, it is very expensive to simulate complex two-phase flows in a long
pipe with today's computer capability. To reduce the computational cost and obtain basic
and essential flow properties of industrial interest, such as gas volume fraction, liquid and
gas velocity, pressure, a one-dimensional model is necessary. The two-fluid model is
considered to give a realistic prediction for the gas-liquid flow inside a pipeline.
The two-fluid model (Wallis, 1969; Ishii, 1975), also known as the separated flow
model, consists of two sets of conservation equations for mass, momentum and energy
for the gas phase and the liquid phase. Although it has success in simulating two-phase
flow in a pipeline, the two-fluid model suffers from an ill-posedness problem. When the
slip velocity between liquid and gas exceeds a critical value that depends on gravity and
liquid depth, among other flow properties, the governing equations do not possess real
characteristics (Gidaspow, 1974; Jones and Prosperettii, 1985; Song and Ishii, 2000).
This ill-posedness condition suggests that the results of the two-fluid model under such
condition do not reflect the real flow situation in the pipe. The two-fluid model only gives
meaningful results when the relative velocity between the gas and liquid phase is below
the critical value. However, this critical value coincides with the stability condition of
inviscid Kelvin-Helmholtz instability (IKH) analysis (Issa and Kempf, 2002). Because
the IKH instability results in the flow regime transition from the stratified flow to the slug
flow or annular flow (Barnea and Taitel, 1994a), ill-posedness of two-fluid model has
been interpreted as to trigger the flow regime transition (Brauner and Maron, 1992;
Barnea and Taitel, 1994a).
The computational methods for solving the two-fluid model have been investigated
by many researchers. For computational simplicity, it is further assumed that both liquid
and gas phases are incompressible. This is valid because most stratified flows are at
relatively low speed compared with the speed of sound. To solve the incompressible
two-fluid equations, one approach is to simplify the governing system to only two
equations for liquid volume fraction and liquid velocity and neglect the transient terms in
the gas mass and momentum equations (Chan and Banerj ee, 1981; Barnea and Taitel,
1994b). A more effective method is to use a pressure correction scheme (Patanka 1980).
Issa and Woodburn (1998), and Issa and Kempf (2003) applied the pressure correction
scheme for the two-fluid model and simulated the stratified flow and the slug flow inside
a pipe.
When two-fluid model becomes ill-posed, the solution becomes unstable. A good
discretized model should be capable of capturing the incipience of the instability point.
However, numerical instability may not be the same as the instability caused by the
ill-posedness. Lyczkowski et al. (1978) used von Neumann stability analysis to study a
compressible two-fluid model with their numerical scheme and found that numerical
instability and ill-posedness may not be identical. However, their two-fluid model lacked
the gravitational term and the study focused on one specific discretization scheme and is
thus incomplete. Stewart (1979), Ohkawa and Tomiyama (1995) attempted to analyze the
numerical stability of an incompressible two-fluid model with a simplified model
equation as an alternative. Their study showed that higher order upwind schemes yield a
more unstable numerical solution than the 1st order upwind scheme.
In this study, a pressure correction scheme is employed to solve the two-fluid
model. It is designed to increase the computational stability when the flow is near the
ill-posedness condition. The von Neumann stability analysis is carried out to study the
stability of the discretized two-fluid model with different interpolation schemes for the
convection term. For the wave growth rates using the 1st order upwind, 2nd order upwind,
QUICK, and central difference schemes, the central difference scheme is more accurate
and more stable. Excellent agreement for the wave growth rates is obtained between the
analysis and the actual computation under various configurations.
4.1.2 Governing Equations
The basis of the two-fluid model is a set of one-dimensional conservation equations
for the balance of mass, momentum and energy for each phase. The one-dimensional
conservation equations are obtained by integrating the flow properties over the
cross-sectional area of the flow, as shown in Figure 4-1.
> y Gas velome
Gas phase Gas velocity u h faction a.
Liquid velome
h Y faction a I
Liquid phase Gravity g action
Liquid velocity uz
Pipe cross section
Figure 4-1. Schematic of two-fluid model for pipe flow.
Because the ill-posedness originates from the hydrodynamic instability of the
two-fluid model, only continuity and momentum equations are considered in the inviscid
two-fluid model. Furthermore, no mass and energy transfer occurs between two phases.
Surface tension is also neglected since it only acts on small scales, while the waves
determining the flow structure in pipe flows are usually of long wavelength. The gas
phase is assumed to be incompressible, as the Mach number of the gas phase is usually
very low for the stratified flow. Hence, the mass conservation equations for liquid phase
is
(a,)+ (ua,) = O, (4.1)
at ax
where a, is liquid volume fraction, p, is liquid density, u, is the liquid velocity, t is the
time, x is the axial coordinate.
The liquid layer momentum conservation equation is
S(u,)+a ua= -- -g cos PH, -gg sin (4.2)
at ax I P A x ax
where p, is the pressure at the liquid-gas interface, g is gravitational accelerator, 3 is the
angle of inclination of the pipe axis from the horizontal lane, and H1 is the liquid phase
hydraulic depth. It is defined as
H1- a1 (4.3)
aa, lah a'
where h, is the liquid layer depth. The second term on the right hand side of Equation
(4.2) represents the effect of gravity on the wavy surface of liquid layer.
The gas phase mass conservation equation is
a (a)+a g)=0, (4.4)
where pg, ag, ug are density, volume fraction, and velocity of gas phase. It is noted that
a, + ag = 1. (4.5)
The momentum equation for gas phase is
S(Ugag)+ (ua g c H -g cosH -ggsin (4.6)
at ax Pg x ax
whereHg is the gas phase hydraulic depth. It is defined as
ag a
H =- ,- (4.7)
o g/ah aI
where h, is the gas layer depth,
4.1.3 Theoretical Analysis
4.1.3.1 Characteristic analysis and ill-posedness
It is well known that the initial and boundary conditions need to be imposed
consistently for a given system of differential equations. The condition is well-posed if
the solution depends in a continuous manner on the initial and boundary conditions. That
is, a small perturbation of the boundary conditions should give rise to only a small
variation of the solution at any point of the domain at finite distance from the boundaries
(Hirsch, 1988).
Equations (4.1, 4.2, 4.4 and 4.6) form a system of 1st order PDEs, for which the
characteristic roots, A, of the system can be found. If A's are real, the system is
hyperbolic. Complex roots imply an elliptic system, which causes the two-fluid model
system to become ill-posed because only initial conditions can be specified in the
temporal direction. Any infinitesimal disturbance will cause the waves to grow
exponentially without bound when A's are complex valued.
Let Ube the vector(a1, u,ug p)T Equations (4.1, 4.2, 4.4 and 4.6) can be written
in vector form as
[A] + [B] = [C], (4.8)
at dx
where [A], [B] and [C] are coefficient matrices, given by
1 0 0 0
-1 0 0 0
[A]= u (4.9a)
ui Oa 0 0
Ug 0 a 0
[B] = u1 + gH, cos f 2au, 0 a, (4.9b)
-u +gH cos / 0 2au, a
Pg
0
[C] = (4.9c)
-a,g sin
aOg sin 3
The characteristic roots of the system is determined by solving A from the
following
[A] [B] = 0. (4.10)
where denotes the determinant of the matrix. Substituting Equations (4.9a) and (4.9b)
into Equation (4.10) results in
S- u 0 0
-(2-g) 0 -ag 0
u, ( u, )-gH, cos f a, ( 2u ) 0 -- = 0. (4.11)
-ug ( -ug) gH cos/ 0 ag,(-2ug) g
Pg
After expansion of the above determinant, the characteristic polynomial for A is
obtained:
Pg (A-U)2 + (- 2 g) gcos =0. (4.12)
ag a, a,
The roots are
S ,+ Pgg P1 Pgsin g ( 2
A = a ag (4.13)
AP Pg
a,1 ag
When g = 0, Equation (4.13) can have real roots only if 1 = ug = u,. Otherwise, the
two-fluid model is ill-posed (Gidaspow, 1974). If g 0, the well-posedness with real
roots requires
U2= u )2
A Pg a,
Equation (4.14) gives the critical value AUc for the slip velocity AU between two
phases beyond which the system becomes ill-posed. The two-fluid model stability
criterion from the characteristic analysis is exactly the same as that from the IKH analysis
on two-fluid model by Barnea and Taitel (1994) as shown below.
4.1.3.2 Inviscid Kelvin-Helmholtz (IKH) analysis and linear instability
IKH analysis (Barnea and Taitel, 1994) provides a stability condition for the
linearized two-fluid model as well as useful information on the growth rate of an
infinitesimal disturbance in the two-fluid model.
Splitting the flow variables into the base variables and the small disturbances, such
as a, + a, expressing the disturbances on the form of
a, = exp(I(at kx)), (4.15a)
u, = E exp(I(ac kx)), (4.15b)
Ug = E exp(I(a kx)), (4.15c)
p = exp(I(a kx)), (4.15d)
where "~" denote disturbance value, I = denotes imaginary unit, cis the amplitude
of perturbation, cois the angular frequency of wave and k is the wavenumber.
Substituting them into the differential governing equations (4.1, 4.2, 4.4 and 4.6), and
linearizing the resulting equations, the following system is obtained for the disturbance
amplitudes, E_, eg, Fe ,E
w-ulk -ak 0 0
w -ugk 0 agk 0 E
H k E,
-k gcosf wc-u,k 0 =0. (4.16)
al Pi g
H k
-k H gcosp 0 co-uk ep
a, Pg
For non-trivial solutions to exist, the following dispersion relation between the wave
speed c and the angular frequency Cmust hold
( ag a1 alag
c =-- g (4.17)
k P+ Pg
a, a,
It is note that the negative imaginary part of determines the growth rate of disturbance.
Equation (4.17) is identical to Equation (4.13), only with A being replaced by c. Details
of the derivation for IKH stability condition can be found in Barnea and Taitel (1994).
4.1.4 Analysis on Computational Instability
4.1.4.1 Description of numerical methods
In general, the governing equations (4.1, 4.2, 4.4, and 4.6) are solved iteratively.
The basic procedure is to solve the continuity equation of liquid for the liquid volume
fraction, and the liquid and gas phase momentum equations for the liquid and gas phase
velocities. To obtain a governing equation for the pressure, Equation (4.1) and Equation
(4.4) are first combined to form a total mass conservation,
(ugag)+ (ua,)= 0. (4.18)
ax ax
Substituting the liquid and gas momentum equations into the above leads to
Sa, + ap a 2 (a 2 2\
u +a
9p a 1 (4.19)
+- g cos H, -a,g sin + g cosf3Hg aggsin \.
x ( ax ax
To solve the pressure equation, SIMPLE type of pressure correction scheme (Patanka,
1982; Issa and Kempf, 2002) is used in this study.
A finite volume method is employed to discretize governing equation. A staggered
grid (Figure 4-2) is adopted to obtain compact stencil for pressure (Peric and Ferziger,
1996). On the staggered grids, the fluid properties such as volume fractions, density and
pressure are located at the center of main control volume, and the liquid and gas
velocities are located at the cell face of main control volume. Figure 4-2 shows the
staggered grids arrangement.
SVelocity control
uw Ue volume
Pw Pp PE
Pw Pe
Main control w aP
volume
Figure 4-2. Staggered grid arrangement in two-fluid model.
The Euler backward scheme is employed for the transient term. The discretized
liquid continuity equation becomes
S(a,) (a )P ))+ (a, ) (a,u 2l) = 0, (4.20)
At
where the superscript 0 denotes the values of the last time step. The subscript P refers to
the center of the main control volume, and subscripts e and w refer to the east face and
west face of main control volume, respectively. The liquid velocity on the cell face is
known, and the volume fraction on the cell face can be evaluated using various
interpolation schemes. Among them, central difference (CDS), 1st order upwind (FOU),
2nd order upwind (SOU) and QUICK schemes are commonly used. Equation (4.4) for the
gas phase is similarly discretized.
The liquid momentum equation is integrated on the velocity control volume. Using
similar notations, one obtains
x ((au ,)p (a1u1 )+ (u), (a1U1), (u1 ) (a1u,) =
t ( a, (4 .2 1)
(pI, p )+ ((a, )w (ae ) )H, g cos Ax( ) g sin 8
where P, e, w refer to the center, east face and west face of the velocity control volume,
respectively. The cell face flux is the liquid velocity, which is obtained by using central
difference, and the volume fraction and liquid velocity at the cell face, which are
transported variables, can be interpolated by using different schemes. It is important to
note that the interpolation method used for the Equation (4.21) must be exactly the same
as those for Equation (4.20). For example, if FOU is used in Equation (4.20), the cell face
flux on the east face of velocity control volume in Equation (4.11) is
(a,u, )e (u, )e = (a,u, )pMAX((u, ),O)- (ao,l )EMAX(- (u, )e,0). (4.22)
If CDS is used in Equation (4.20), the cell face flux on the east face in Equation (4.21) is
evaluated as
(aeu, )P + (a1,1 )E
(a,))P + (u) ) (u, )e (4.23)
2
Using similar discretization procedure, the gas phase momentum equation is
integrated:
"- ((agug )p (agu )+ (ug ) (agug (ug (ag )
a ) (4.24)
(p + ((a)w -(a ))Hgcos/- Ax( )Pgsin .
Pg
For convenience, the discretized mass or momentum equations are written in a
general form
Ap~ + A, +E + Aw~, = B, (4.25)
where ( is the variable to be solved, A is the coefficient, B is the general source term.
For the pressure correction scheme, Equation (4.18) is integrated across the main
control volume. The discretized equation is
(aO gUg (ag g) + (aui )j (aou ), = 0. (4.26)
Because Equation (4.18) is obtained by combining Equation (4.1) and Equation
(4.4), the discretization scheme for Equation (4.26) should be exactly the same as those
for Equation (4.20) and the discretized equation of Equation (4.3). For instance, if CDS is
used in Equation (4.20), it must be used in the main control volume for Equation (4.26):
(lg )e (ag)P ( )E (g )w (g )P + (g )w
ge 2 g 2
2 2 (4.27)
(a-(u + E(1 P.
+ )( )e 2 =) \ ) 0.
2 2
The final discretized pressure equation is obtained by substituting these two
momentum equations, Equation (4.21) and Equation (4.24) into Equation (4.26). This
yields
ap p +aeE + awp = b,
(4.28)
a (a,) +(a,)E a (a)p +(a)E, a,- (4.29a)
a ( = -a a, (4.29c)
2 pPT 2 p, (A
(aj, + (a,), a (( + (a,(,,)a
W2 F
2 2 (4.29d)
+ (a l)P + (a ) (aP + (a
2 2
where, p' represents the pressure correction value, u* represents the imbalanced
velocity, and Ap is from the corresponding discretized liquid or gas momentum equation,
Equation (4.15). The flow chart of the pressure correction scheme is shown in Figure 4-3.
Similarly, the pressure correction schemes with FOU, SOU, CDS, and QUICK can
be obtained.
Consistently handling the discretization is critical to the reduction of numerical
diffusion and dispersion. Barnea and Taitel (1994) showed that the viscosity of fluid can
dramatically degrade the stability of two-fluid model through viscous Kelvin-Helmholtz
stability analysis. Although the viscosity in two-fluid model appears as the body force
instead of 2nd order derivative terms in the modified governing equation, it is
hypothesized that the numerical diffusion and dispersion appearing as derivative in the
modified governing equations produce similar impact on the stability of two-fluid model.
Figure 4-3. Flow chart of pressure correction scheme for two-fluid model.
4.1.4.2 Code validation- dam-break flow
The pressure correction scheme is first validated by computing the transient flow
due to dam-break flow (Figure 4-4). The liquid flow is assumed to be over a horizontal
flat surface and the flow is assumed to be one-dimensional. On the left side of the dam is
a body of stationary water in the reservoir with the flat surface of height H. On the right
side of dam is a dry river bottom surface. After the dam breaks suddenly, the water in the
reservoir flows to the downstream due to the gravitational force. If there are no friction
between the fluid and the wall and no viscosity inside the fluid and air pressure is a
constant, an analytical solution for the liquid velocity based on St Venant equation can be
found (Zoppou and Roberts, 2003). The result is shown in Table 2.1.
y
S\Dam
H
H/ Dry river plate
x=- x
Figure 4-4. Schematic for dam-break flow model.
To solve dam-break flow, the pressure at interface, the vapor phase density and
velocity are set to zero. Second order upwind scheme as the cell face interpolation
scheme is implemented in the pressure correction scheme. Figure 4-5 compares water
depth between the present numerical solution and the analytical solution at t=50s. Two
solutions match very well except at the tail end of the liquid, where the numerical
solution is smooth due to a little numerical dissipation. Figure 4-6 compares liquid
velocities between the numerical and analytical solutions at t=50s. Again, these two
solutions match very well except at the leading and tail ends. The discrepancy at the
leading end is due to that the liquid layer is too thin and the numerical result is prone to
error.
Table 4-1. Analytical solution for dam-break flow (Zoppou and Roberts, 2003).
x( position) u (water velocity) h (water depth)
x < -tgH u=0 h=H
-t
x 2t gH u=0 h=0
Although only the dynamics of liquid phase is considered in the dam-break flow, it
is still a solid step for validating the coupling of the pressure and liquid flow (liquid
volume fraction and liquid velocity) in numerical scheme. When both the liquid and the
gas phase present in the flow, the instability in two-fluid model rises due to the
interaction of the liquid and the gas phase, when the slip velocity is large. Numerical
instability of pressure correction scheme emerges and destroys the numerical results
when the two-fluid model near ill-posedness. This numerical instability will be
investigated in the next section and the code will be validated using the theoretical results
of inviscid Kelvin-Helmholtz analysis (Barnea and Taitel, 1994).
10 \ t-osec
9
8
Numerical
7 Analytical
t=50 sec
6
4
3
2
1 -
0 I
0 500 1000 1500 2000
z(m)
Figure 4-5. Water depth at t=50 seconds after dam break.