
Citation 
 Permanent Link:
 http://ufdc.ufl.edu/AA00029735/00001
Material Information
 Title:
 Modeling of bulk evaporation and condensation with internal heating at various gravity levels
 Creator:
 Ding, Zhongtao, 1960
 Publication Date:
 1995
 Language:
 English
 Physical Description:
 xv, 126 leaves : ill. ; 29 cm.
Subjects
 Subjects / Keywords:
 Condensation ( jstor )
Convection ( jstor ) Evaporation ( jstor ) Free convection ( jstor ) Heat transfer ( jstor ) Heating ( jstor ) Liquids ( jstor ) Temperature distribution ( jstor ) Vapor phases ( jstor ) Vapors ( jstor ) Dissertations, Academic  Mechanical Engineering  UF Mechanical Engineering thesis, Ph. D
 Genre:
 bibliography ( marcgt )
theses ( marcgt ) nonfiction ( marcgt )
Notes
 Thesis:
 Thesis (Ph. D.)University of Florida, 1995.
 Bibliography:
 Includes bibliographical references (leaves 119125).
 Additional Physical Form:
 Also available online.
 General Note:
 Typescript.
 General Note:
 Vita.
 Statement of Responsibility:
 by Zhongtao Ding.
Record Information
 Source Institution:
 University of Florida
 Holding Location:
 University of Florida
 Rights Management:
 Copyright Zhongtao Ding. Permission granted to the University of Florida to digitize, archive and distribute this item for nonprofit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
 Resource Identifier:
 34428449 ( OCLC )
022365266 ( ALEPH )

Downloads 
This item has the following downloads:

Full Text 
MODELING OF BULK EVAPORATION AND CONDENSATION WITH INTERNAL
HEATING AT VARIOUS GRAVITY LEVELS
By
ZHONGTAO DING
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 1995
UNIVERSITY OF FLORIDA LIBRARIES
Copyright 1995
by
Zhongtao Ding
To my wife Zhaohui Cheng
ACKNOWLEDGMENTS
It has been a wonderful experience to pursue further education and academic development at the University of Florida, where I have known so many nice professors and friends who have helped me help in many ways. From them I learned that the way to achieve goals through hard work, determination, and effort.
First of all, I appreciate the insightful guidance, constant support and encouragement provided by Professor Samim Anghaie. Without his enthusiasm and commitment to this challenging project, my exploration in this work could not have been accomplished. His help in finding solutions to my problems both personal and professional will never be forgotten.
I sincerely thank Professor James F. Klausner, Chairman of the supervisory committee, for all his spiritual encouragement and technical assistance during the course of this research. His helpful suggestions and comments are greatly appreciated.
I would like to express my special thanks to Professor Wei Shyy with greatest sincerity. I owe much of my knowledge on computational modeling of phase change problems to his direct and sincere help. I have learned a lot from the discussions with him and his students, from his class, his books and his many papers.
I also want to extend my appreciation to Professors C. K. Hsieh and R. Mei for their serving on the supervisory committee. Their suggestions have always been constructive and helpful.
I would like to express appreciation for all spiritual encouragement and financial support from Professor Nils J. Diaz.
iv
Finally, I want to thank my parents and sister who have supported me unconditionally and expect nothing from me in return. I wish to thank my wife, Zhaohui Cheng. I can not imagine that I could have finished my research and dissertation without her understanding, support and patience throughout the program.
TABLE OF CONTENTS
pa e
ACKNOW LEDGM ENTS ........................................ ......... ................ iv
LIST OF TABLES ........................................ viii
LIST O F FIG U R E S..................................................... .............. ix
N O M EN CLA TU R E ....................................... ................. ...... ... ........... xi
A B S T R A C T ........................................ ................... ................... ................. xiv
CHAPTERS
1 IN TROD U CTION ........................................... 1
2 LITERATU RE SURVEY ........................................ ............ ...............9
2.1 SolidLiquid Phase Change Problems ................................. 9
2.2 Numerical Modeling of SolidLiquid Phase Change Problems ................... 11
2.3 LiquidVapor Phase Change Problems ........................... .................. 12
2.4 Natural Convection in the Single Phase ....................................... 13
2.5 Phase Change with Convection ...................... ..................... 14
3 1D CONDUCTION BULK EVAPORATION AND CONDENSATION........ 15
3.1 The Internal Energy Formulation ...................................... ............... 19
3.2 A nalysis and A ssum ptions ........................................ ................ 22
3.3 The Vapor Phase Fraction Update ........................................ 28
3.4 Properties of the Control Volume Containing Two Phases ........................ 32
3.5. Results and D iscussion ................................... ................ 33
3.5.1 A ccuracy A spect ........................... ..................... ....... ........... 33
3.5.2 Comparison of H Formulation and E Formulation ............................ 38
3.5.3 Comparison of EBased and TBased Update Methods ...................38
3.5.4 Bulk Evaporation and Condensation ............... .......... . 40
3.6 Concluding Remarks ................ 41
vi
4 BULK EVAPORATION/CONDENSATION WITH INTERNAL HEATING
AT ZERO GRAVITY ................................... ....... . ..... .......... 45
4.1 The Governing Equations ...... ................ ................. ................. 47
4.2 Conductivity at Control Volume Surfaces ..................................... ... 49
4.3 Results and Discussion ........................................ ....... ................. 51
5 NATURAL CONVECTION WITH INTERNAL HEATING IN
SINGLEPHASE AND TWOPHASE SYSTEMS .......................... 59
5.1 External Natural Convection in a SquareCavity Benchmark Problem ....... 60 5.2 Natural Convection in a SquareCavity with Internal Heating .................... 67
5.3 Natural Convection in LiquidGas Systems without Internal Heating ......... 74
5.3.1 Gas Convection Induced by Liquid Convection ................................ 75
5.3.2 B oth Phases H eated ...................................................................... 76
5.4 Natural Convection in LiquidGas Systems with Internal Heating .............. 80
5.4.1 Case 1: Rectangular Cavity ........... . ............ ................ 80
5.4.2 Case 2: Cylindrical Container ... ......... ................ ....... ........ 81
6 BULK EVAPORATION/CONDENSATION WITH INTERNAL HEATING
AT MICROGRAVITY AND NORMAL GRAVITY LEVELS ..................... 84
6.1 Brief Preview ................ .................. 84
6.2 The Governing Equations ....... .. ............................... 85
6.3 N um erical Treatm ent ....... ......... ............................... 88
6.4 Results of M icroGravity (10  g) ........................ ...................... 90
6.5 Results of Normal Gravity (1 g ) ................................ .................... 106
7 CONCLUSIONS AND RECOMMENDATIONS ....................................... 115
REFEREN CES .............................. ............................................... ....... 119
BIOGRAPHICAL SKETCH ........................................................................... 126
vii
LIST OF TABLES
Table page
2.1 Brief summary of major characteristics of phase change problems .................. 10
5.1 Comparison of present solutions with benchmark solutions ............................... 63
viii
LIST OF FIGURES
Figure page
1.1 Schematic diagram of a typical nuclear thermionic fuel element ...................... 2
1.2 Schematic diagram of the multiphase nuclear fuel system .............................. 4
3.1 Schematic diagram of the bulk evaporation and condensation ..................... 17
3.2 The control volume containing two phases ........................ ................ 33
3.3 Solution features for bulk evaporation ...... ................ .................. 35
3.4 Comparison of two formulations for evaporation ........................................ 36
3.5 Comparison of two update methods for evaporation ..................................... 37
3.6 Bulk evaporation and condensation ............................................ 42
4.1 Schematic diagram of the twophase cell operated at 0g .............................. 46
4.2 The control volume contacting other material ...................... ................ 50
4.3 Parameter for heat generation rate ............. ........... ................ 52
4.4 The evolution of liquidvapor interface (0g) .......................................... 56
4.5 The evolution of temperature field (0g) ....................... ................... 57
4.6 Temperature surf of the system at steady state (0g) ............ ............ 58
4.7 Variation of saturation temperature (0g) ..................................................... 58
5.1 Domain and boundary conditions for the benchmark problem ...................... 61
5.2 Flow patterns and isotherms of buoyancy driven natural convection ............ 64
5.3 Nusselt number distribution on side walls ......... ................ ............. 65
5.4 Temperature profile at horizontal centerline ................... ...................... 66
ix
5.5 Domain and boundary conditions for internal heating induced convection ...... 70 5.6 Flow pattern and isotherms of internal heating induced convection ( case 1 ) ...71 5.7 Flow pattern and isotherms of internal heating induced convection ( case 2) ...72 5.8 Nusselt number distribution on the side walls .......................................... 73
5.9 Natural convection of twophase without internal heating ( case 1) .............. 78
5.10 Natural convection of twophase without internal heating ( case 2 ) ............... 79
5.11 Natural convection of twophase with internal heating ( case 1 ) ................... 82
5.12 Natural convection of twophase with internal heating ( case 2) .................. 83
6.1 Schematic diagram of mass reflux ..................................... 96
6.2 The evolution of liquidvapor interface ( 103 g ) .................................... 98
6.3 The evolution of isotherm ( 10 g ) ...... ............................................. 99
6.4 Temperature surf at x = 2.0 ( 10 g ) ..................... ... ...................... 100
6.5 Variation of saturation temperature ( 10 g ) ......................................... 100
6.6 Solutions with different grid sizes and schemes ........................................ 101
6.7 The evolution of liquidvapor interface with coarse grid .............................. 101
6.8 Flow patterns in vapor and liquid phases r = 2.0 ( 103 g ) ........................ 102
6.9 Velocity U profile at different axial positions (10 g ) ............................... 103
6.10 The evolution of liquidvapor interface ( 1g ) .......... ............ .. ............. 109
6.11 The evolution of isotherm s ( g ) ......................................... ............... 110
6.12 Tem perature surf at = 2.0 ( 1g ) ............................................................. 111
6.13 Variation of saturation temperature ( 1g ) ................................. 111
6.14 Flow patterns in vapor and liquid phases r = 2.0 ( 1g) ............................. 112
6.15 Velocity U profile at different axial positions ( 1g) ............... ................. 113
6.16 Comparison of interface and temperature distribution at different g levels ..... 114
X
NOMENCLATURE A, C, D constants E internal energy [J] G internal heat generation rate [J/kg] H enthalpy [J] Q heat quantity [J] T temperature [K] U dimensionless velocity V volume [m3] or dimensionless velocity W work [J] X, R dimensionless coordinates c, specific heat at constant pressure [J/kg/K] cv specific heat at constant volume [J/kg/K] e specific internal energy [J/kg] h specific enthalpy [J/kg] f vapor phase fraction g gravity acceleration [m/s2
1 length scale [m] p pressure [bar] t time [s]
xi
u, v velocity [m/s] x, r cylindrical coordinates [m]
Nondimensional Characteristic Numbers Gr Grashof number Pe Peclet number Pr Prandtl number Re Reynolds number Ra Rayleigh number St Stefan number
Greek symbols
a thermal diffusivity [m2/s], or thermal diffusivity ratio e phase change temperature interval K thermal conductivity [W/K/m], or its ratio O order of magnitude
0 dimensionless temperature p density [kg/m'], or its ratio T dimensionless time Yr stream function
xii
Subscripts
f interface position, mixed phase i interface I liquid phase r characteristic scale s, sal saturation point v, vap vapor phase
Superscripts
k iteration level n time level arrow bar vector bar  nondimensional quantity
xiii
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 OF BULK EVAPORATION AND CONDENSATION WITH INTERNAL HEATING AT VARIOUS GRAVITY LEVELS
By
Zhongtao Ding
December, 1995
Chairman: James F. Klausner CoChairman: Samim Anghaie Major Department: Mechanical Engineering
A comprehensive numerical model is developed to simulate the dynamic processes of bulk evaporation and condensation in an encapsulated container, associated with internal heat generation under different gravitational strengths. A complete set of equations governing the conservation of mass, momentum and energy is solved. The thermal performance of a multiphase nuclear fuel element at zero gravity, microgravity, and normal gravity is investigated. In the simulation, the numerical solutions have revealed much useful information, including the evolution of the bulk liquidvapor phase change process, the evolution of the liquidvapor interface, the formation and development of the liquid film covering the side wall surface, the temperature distribution, the convection flow field, and the strong dependence of the thermal performance of such multiphase fuel cells on the gravity conditions.
xiv
In regard to the development of computational techniques, the internal energy formulation for the constant volume phase change process is proposed. The heat transfer mechanism for this type of phase change problem, including the convection effect induced by the density change, is analyzed and discussed. The computational performance of the internal energy formulation is evaluated. The performance of two update methods for the vapor phase fraction, the Ebased and Tbased methods, is investigated. Simulation of a onedimensional conductiondriven bulk evaporation and condensation problem is performed and the evolution of the phase change process is obtained.
Some natural convection problems with or without internal heating in the singlephase and twophase systems are studied. The solutions agree quite well with benchmark solutions and experimental results. It is revealed that interactions between the convection cells in liquid and vapor phases can substantially complicate the flow fields and the computational modeling.
XV
CHAPTER I
INTRODUCTION
The thermionic energy converter is a nonmechanical gaseous electronic device for converting heat energy directly into electric potential through thermionic electron emission. In a thermionic diode, electrons are emitted from a hot electrode, which is the emitter, and collected by a colder electrode, which is the collector, at a higher potential energy (lower electrical potential). Part of the heat removed from the emitter due to evaporating electrons is rejected to the collector by condensing electrons. The remaining heat is converted into electric power in the load as electrons return to the emitter potential. The physics governing thermionic energy conversion is well understood (Hatsopoulos and Gyftopoulos, 1973), and this energy conversion technique is attracting wide interest in its application toward a variety of power conversion systems (Lee et al., 1993, Young et al., 1993). In a review paper on thermionics, Rasor (1991) outlined the history, application options, and ideal performance of thermionic energy converters and described the basic plasma types associated with various modes of converter operation.
Various types of thermionicconverter systems, such as solar thermionic generators, radioisotope thermionic generators, chemical thermionic generators and nuclear reactor thermionic generators, have been or are currently being developed. The nuclear reactor thermionic generators are simply thermionic converters supplied with heat produced by nuclear fission. The fission heat is conducted from the reactor core to the emitter, and the rejected heat is conducted from the collector and carried away via a coolant. A typical schematic configuration of the thermionic fuel element (TFE) fueled with solid nuclear fuel U02 is shown in Figure 1.1. The nuclear fuel resides in the center
1
2
of the TFE. The space between the emitter and collector is filled with an easily ionizable rarefied medium, typically cesium vapor, in order to enhance the electrical performance. A liquid metal coolant flows in the annulus channel surrounding the TFE. Recently a variety of investigations involving such a TFE have been carried out that consider the thermal and electrical performance, system simulation and design, and application as a space power system, Lewis et al. (1991), McVey and Rasor (1992), Young et al. (1993).
insulator gap
void
coolant
channel
fuel
\cladding
emitter collector
Figure 1.1 Schematic diagram of a typical nuclear thermionic fuel element
3
The thermal performance of a thermionic fuel element is tightly coupled with its electrical performance. The electrical output current density is directly related to the emitter temperature, that is, the current density increases significantly with emitter temperature, which can be illustrated by a simple example of an ideal diode thermionic converter as follows
J= ATE2 exp[(V+ )/ kTE] (1.1) where J is the current density, A and k are constants, V is the output, c is the collector work function, and TE is the emitter temperature.
Consequently a nuclear fueled TFE, which can provide high temperature, is well suited for this application. The maximum operating temperature of solid nuclear fuel UO2 could be as high as 1800 K. In order to further enhance the TFE output, a higher operating temperature is required. Such higher temperatures could be realized if liquid nuclear fuel is employed. If nuclear fuels such as UF4 or UF4/UO02 in a liquid or vapor state were used, they could sustain a temperature of 2600 K or higher. The large increases in operating temperature will raise the power output significantly. The use of a multiphase nuclear fuel reactor in conjunction with a thermionic converter provides the potential for a compact, higher power output electrical generator that is well suited for space applications.
Certain problems related to the thermal performance of TFE arise with regard to the application of multiphase nuclear fuel reactors. These include the phase change (PC) process between the liquid and vapor, the behavior of evaporation and condensation associated with internal heating and external cooling, the location and motion of the liquidvapor interface, the temperature distribution, the influences of gravitational fields, and so on. With such a system, particular attention must be paid to limiting the vapor pressure that is safe for the containment vessel because evaporation in a closed container
4
will raise the vapor pressure very rapidly. A thorough investigation of these issues will be required to fully apply this concept.
insulated
vapor liquid film
conden ation cooled
evaporation liquid Qg,, pool
insulated Figure 1.2 Schematic diagram of the multiphase nuclear fuel system
A cylindrical container filled with multiphase nuclear fuel, such as UF4, is shown in Figure 1.2. The fuel experiences liquidvapor phase change processes: evaporation and
5
condensation. The internal heat provided by the nuclear fission reaction evaporates the liquid, while the side wall cooling tends to condense the vapor. Globally, the evaporation of liquid and condensation of vapor take place simultaneously. The evaporation and condensation processes could be driven by either pure heat conduction or by both heat conduction and convection, depending mainly on the gravity field. Under a zerogravity condition, buoyancydriven convection is nonexistent. However, under microgravity or onegravity fields, buoyancyinduced convection is present. The location of the liquidvapor interface may vary depending on the mode of heat transfer, the internal heat generation intensity and spatial distribution, and the gravity field. Such an interface could be defined as a free and moving interface. To solve the phase change problem analytically or numerically, the free and moving interface is not known in advance, so that its position has to be determined as a part of the solution. Thus the solution to this heat transfer problem is quite complex.
With such background, specific issues are to be investigated. For the phase change processes between liquid and vapor, such as bulk evaporation and condensation involved in the multiphase system, no similar investigation has been found. Limited investigations have been conducted. Investigations of this problem in different parametric regimes have been performed by Shyy (1994b). To the best of the author's knowledge, no prior researcher has considered the general bulk condensation and evaporation processes in an enclosed container, under various gravity fields, with internal heat generation. Liquidvapor phase change problems are also characterized by a large step change in the density across the phasic boundary. The singular jump of the density at the interface, other than the jump of temperature gradient which exists unavoidably for all phase change problems, brings substantial challenges to computational fluid dynamics and heat transfer.
Since the investigations for the solidliquid phase change problems, such as melting and solidification, are numerous, some successful experiences and lessons can be drawn from those efforts and applied to the study of the bulk evaporation and condensation
6
problems. Appropriate modeling techniques, such as formulation and numerical treatment, must be developed. For the solidliquid phase change problems, the thermodynamic and transport properties between two phases vary moderately during phase change processes. These properties, however, vary significantly or dramatically during the phase change processes between liquid and vapor. This causes some difficulties. Scaling and order of magnitude analysis need to be conducted in order to make simplifications that can lead to a concise and efficient solution procedure for the governing equations. Some special treatment and closure relationships, such as those for vapor phase fraction, the vapor density, the saturation pressure and temperature, have to be applied and implemented into the computation.
Like all phase change processes, evaporation and condensation are also inevitably subject to the effect of convection induced by either buoyancy or surface tension variations. The existence of the convection can substantially affect the flow field, the temperature field, and subsequently the performance of the phase change process. The effect of the convection on the evaporation and condensation processes is to be investigated. For the conventional solidliquid phase change problems, heat is usually provided externally from the boundary to the system. For the nuclear system, heat is generated internally via a nuclear fission reaction. The problems of phase change induced by internal heating are different from those induced by external heating. The internal heating, like external heating, results in temperature and density gradients, and induces the buoyancydriven flow. Such natural convection induced by internal heating is defined as internal natural convection, to be distinguished with the classical natural convection caused by external heating/cooling. This internal natural convection has drawn certain attention for some investigations (Acharya and Goldstein, 1985; Lee and Goldstein, 1988). The behavior of the internal natural convection and its influence on the phase change processes, in particular on the bulk evaporation and condensation processes, are to be investigated.
The objective of this dissertation is to develop a numerical method to model the dynamic processes of bulk evaporation and condensation, with internal heat generation and natural convection, in an encapsulated container under different gravitational levels. Based on this model, the following topics will be studied in depth: 1) the bulk evaporation and condensation processes with different modes of heat transfer, that is, conduction and convection, 2) the movement of the liquidvapor interface, 3) the evolution of the temperature field, and 4) the influence of various gravitation fields on the performance of the phase change process. The modeling and results can be applied to validate the concept of a multiphase nuclear system, including the multiphase TFE, and to predict the performance and assist in devising strategies for its thermal management and control for space applications.
The current research status concerning these problems mentioned above is reviewed in Chapter 2. The modeling of a onedimensional, externally heated or cooled and conductiondriven bulk evaporation/condensation problem is described in Chapter 3. This phase change problem is well suited to develop the mathematical formulation, the closure relationships, and the updating of the vapor phase fraction. The results of the model are used to investigate the basic behavior of the evaporation/condensation processes, and to examine the computational efficiency and accuracy.
The numerical approach developed in Chapter 3 is then expanded to investigate the bulk evaporation/condensation processes with internal heat generation at zero gravity level, which is presented in Chapter 4. The evolution of the phase change processes, the temperature field, and the motion of the liquidvapor interface topology are demonstrated. Chapter 5 describes the modeling of the combined internal and external natural convection in the singlephase system and twophase system without phase change. To fully understand the mechanism and behavior of the internal natural convection and to establish the confidence of the computational accuracy, the external natural convection is also computed and comparisons with the benchmark solutions are made. In Chapter 6, the
8
bulk evaporation/condensation processes with internal heat generation and natural convection under different gravitational fields are modeled. A complete set of equations governing the conservation of mass, momentum and energy is solved. The effects of the buoyancyinduced convection via internal heating and external cooling, and different gravitational fields, on the phase change processes are studied. Finally, the conclusions summarizing the achievements of this work and recommendations for further research are provided.
CHAPTER 2
LITERATURE SURVEY
Previous investigations of phase change problems are numerous. Those studies, which can be traced back as early as Clausius, Clapeyron, Stefan, Nusselt, and others, have been conducted both experimentally, analytically and are currently being attempted using numerical simulations and approximation methods. In their annual review of heat transfer papers, Eckert et al. (1990, 1991 and 1992), suggested that most of these current works have focused on melting and solidification, pool and flow boiling, film and droplet evaporation, filmwise and dropwise condensation, with applications in thermal energy storage, solidification of alloys and metals, casting processes, crystal growth, heat pipes, heat pumps, boilers, heat exchangers, and nuclear or chemical reactors. The major characteristics of the phase change problems and the current research status are briefly summarized in Table 2.1.
2.1 SolidLiquid Phase Change Problems
Early analytic work on phase change between solid and liquid include studies by Lame and Clapeyron (1831) and Stefan (1891). The fundamental features of this type of problem and some basic analytical solution methods were discussed by Ozisik (1980), Crank (1984) and Yao et al. (1989). Crank (1984) presented a broad but reasonably detailed account of the mathematical solution, both analytical and numerical, of these problems. Due to the existence of the interface, difficulties were encountered in solving
9
10
these problems to the extent that only a handful of very simple cases can be theoretically analyzed. Other efforts on the analytical and numerical treatment for the interface include the source and sink method that has been used and developed by Hsieh et al. (1993) and Li (1993). In this method, the phase change conditions at the interface are introduced into the heat conduction equation as an impulse internal heat source or sink. Then the equations for both regions can be solved simultaneously.
Table 2.1 Brief summary of major characteristics of phase change problems
Moving Boundary  Phase Change Problems
PC Between Solid  Liquid Liquid  Vapor flow or film or bulk evaporation Process Melting and pool drop and
solidification boiling conden condensation sation
Heating externally from externally externally or /Cooling the system boundary from boundary internally heating Thermo constant or
physical change moderately change significantly properties
Conduction conduction conduction conduction and conduction conduction
and/or only and convection only and convection convection convection Solution analytical numerical experimental Approach numerical empirical numerical numerical
empirical analytical
Research numerous investigations numerous very few found, Status investigations being conducted
11
At the University of Florida, fruitful investigations on various phase change problems have been carried out by Shyy and Chen (1991), Shyy and Rao(1992), Shyy et al. (1992, 1993a, b, 1996), Shyy (1994a, b), Hsieh et al. (1992, 1993), Li (1993), and Hung et al. (1995).
2.2 Numerical Modeling of SolidLiquid Phase Change Problems
Since the 1970s, greater attention and interest have been paid to purely numerical analysis and simulations, which is reflected in the large number of publications covering these topics. Shyy (1994a) and Shyy et al. (1996) summarized a substantial portion of their research activities, and presented the various computational elements important for the prediction of complex fluid flows and interfacial transport.
The numerical attempts can be divided into two groups, based on the choice of dependent variables to be used. In the first group, temperature is the sole dependent variable, and energy equations are written separately in the solid and liquid regions. At the interface, the temperature is continuous, while its gradient is not continuous due to the release or absorption of latent heat. The boundary conditions at the interface are used to couple and connect the temperature in the two regions. The temperature gradient jump is difficult to capture and requires extraordinary effort. In the second group, enthalpy is used as a dependent variable along with temperature. This formulation is termed the enthalpy method, and has gained popularity. In the enthalpy model, the interface is eliminated from consideration in the calculations, and the problem is made equivalent to one of nonlinear heat conduction without phase change. The shortcoming of this method is that it is difficult to determine the interface position due to its oscillations during the numerical solution process. Shamsundar and Sparrow (1975), Shyy (1994a) and Shyy et al. (1996) summarized these features and demonstrated the usefulness of the enthalpy formulation.
12
Many investigators have employed the enthalpy formulation, including Shyy et al. (1992a, 1993b, 1996), Shyy (1994a, b), Voller and Prakash (1987), Voller and Swaminathan (1991), Yao and Prusa (1989).
2.3 LiquidVapor Phase Change Problems
In recent years, interest in the other forms of phase change processes, such as evaporation and condensation, has grown. Carey (1992) systematically summarized the state of the art in boiling and/or condensation phenomena. His book focuses on some fundamental elements and the basic physical mechanism of condensation and vaporization processes and provides a broad review of the research activities in these areas. Among the investigations on condensation and evaporation, those on film condensation are well developed, particularly on laminar film condensation, due to its simplicity. The development of the current theory and the prediction of heat transfer rates associated with film condensation have its origin in the analyses presented by Nusselt (1916). Subsequently, many modifications of this analysis have appeared, each one further relaxing the restrictive assumptions of the Nusselt analysis. Considerations have been given to gravityflow condensation, pure forced convection condensation in the presence of a body force, such as in Sparrow et al. (1967). The influence of entrained noncondensable gas in the vapor and the associated interfacial resistance, that is, a temperature jump at the interface between the condensate and the vapor, on these phase change processes have also been investigated by Sparrow and Eckert (1961), Sparrow and Lin (1964), Sparrow et al. (1967), and Rohsenow (1973). Recently some studies on general film condensation transients and the interfacial condensation model have been conducted by Flik and Tien (1989) and Gerner and Tien (1989). With the development of computational fluid dynamics and numerical heat transfer, more and more attempts have been made to model and simulate condensation phenomena, for example, Jones and Renz
13
(1974), Cossmann et al. (1982), and Gaultier et al. (1993). However, most of these activities, either analytical or numerical, have focused on flow condensation problems with open plates, channels or open containers where the vapor pressure, saturation temperature and latent heat are maintained constant. These problems, where forced convection plays a dominant role, are usually parabolictype problems, and consequently are much simpler than those in the present study in which the saturation temperature, pressure and latent heat continuously change. A few computational simulations of evaporation in containers have been carried out recently. Jang el al. (1990) modeled a heat pipe startup from the frozen state by using onedimensional equations for different phases. Harley and Faghri (1994) performed numerical analysis of thermosyphons including the falling condensate film. Shyy (1994b) developed a computational method for predicting the twophase transient flow and heat transfer characteristics within a constant pressure reservoir of a capillarypumped loop and investigated the phase change of liquid and vapor, under various conditions.
2.4 Natural Convection in the Single Phase
Considerable attention has been given to the study of natural convection of a single phase fluid in enclosures. The temperature gradients across the system that are supplied externally induce density gradients and consequently the buoyancy driven flows. A great deal of research, both numerical and experimental, on various types of convective heat transfer problems has been conducted, including that of Ostrach (1982, 1983, 1988), Gebhart et al. (1988), de Vahl Davis (1983), Hortmann et al. (1990), Shyy et al. (1992), Shyy (1994a), Bejan (1984), Carpenter et al. (1989), Ramaswamy and Jue (1992), Bergman and Ramadhyani (1981). Ostrach (1982, 1983, 1988), Gebhart et aL (1988) and Bejan (1984) reviewed and summarized those activities and achievements. Among the large number of the natural convection problems, the buoyancydriven flow in a square
14
cavity with vertical side walls heated differentially, so called a benchmark problem for natural convection, is a suitable vehicle for testing and validating computer codes. Accurate benchmark solutions have been published by de Vahl Davis (1983) and by Hortmann et atl. (1990) and are available for comparison with the present computations.
Natural convection induced by external heating is of importance in many applications. Natural convection in an enclosure in which the internal heat generation is present is also of prime importance in certain technological applications. Examples are postaccident heat removal in nuclear reactors and some other geophysical or astrophysical problems. In contrast to the extensive investigations on the conventional external heating convection problems, natural convection induced simultaneously by external heating/cooling and internal energy sources has received limited attention, such as Archarya and Goldstein (1985), Lee and Goldstein (1988), and Fusegi et al. (1992). The presence of internal heat generation provides an additional dynamic in overall convective flow systems and offers a complicated yet challenging aspect of the convection under consideration.
2.5 Phase Change with Convection
With no exception, convection induced by either buoyancy or surface tension variations will also be involved in evaporation and condensation like all phase change problems. The involvement of convection can substantially affect the phase change process. Numerous investigations on the solidliquid phase change problems with convection can be found in Ostrach (1983), Crank (1984), Brice (1986), Brown (1988), Gau and Viskanta (1986), Nadarajah and Narayanan (1990), Shyy and Chen (1990,1991), Shyy and Rao (1992, 1994), Shyy et al. (1992, 1996) and Shyy (1994a, b).
CHAPTER 3
1D CONDUCTION BULK EVAPORATION AND CONDENSATION
In recent years, phase change processes such as melting and solidification, boiling and condensation have attracted increasing interest in their development and applications. Among these activities, numerical modeling and simulation of the phase change processes have been the focus of many investigators.
Numerical modeling of heat transfer problems involving phase change has been discussed by a large number of investigators, such as Crank (1984), Yao and Prusa (1989), Voller and Prakash (1987), Voller and Swaminathan (1991), Samarskii et al. (1993), Prakash et al. (1987), Shamsunder and Sparrow (1975), Hsieh and Choi (1992), Hsieh and Akbari (1993), Li (1993), Hung et al. (1995), Shyy and Chen (1990, 1991), Shyy and Rao (1992, 1994), Shyy et al. (1993a, b, 1994, 1996), Shyy (1994a, b). Shyy (1994a) and Shyy et al. (1996) presented a comprehensive discussion of the various computational elements important for the prediction of complex fluid flows and interfacial transport. Most of the noted investigators have approached the solidliquid phase change problems, such as those involved in melting and solidification processes. There has also been an ongoing interest in the investigation of liquidvapor phase change problems, such as evaporation and condensation. More attempts have been made to numerically model and simulate condensation phenomena. In particular, a number of investigators have developed numerical models for flow condensation problems, for example, Jones and
15
16
Renz (1974), Cossman et al. (1982), Gautier et al. (1993), Tournier and ElGenk (1992). Most of these numerical modeling efforts have focused on parabolic flow condensation problems where forced convection plays a dominant role. In this type of condensation problem, the vapor pressure, saturation temperature and latent heat have been treated as constant properties. In the present study, the saturation temperature, and pressure are modeled as dependent variables. In a more recent work, Jang et al. (1990) have modeled a heat pipe startup from the frozen state by using a system of simplified one dimensional equations. Shyy (1994b) has successfully developed a computational method for predicting the twophase transient flow and heat transfer characteristics within a constant pressure reservoir of a capillarypumpedloop and investigated the phase change of liquid and vapor, under various conditions. The present approach, to be discussed later in this chapter, has expanded the mathematical modeling of the phase change problems by including the pressure as a dependent variable and using the internal energy formulation.
The numerical models for the phase change problem can be divided into two groups, based on the choice of dependent variables. In the first group, temperature is the sole dependent variable, and energy equations are written separately for the solid and liquid phases. In the second group, enthalpy is used as the main dependent variable along with the temperature. This formulation which is widely used by many investigators is known as the enthalpy method. In the enthalpy model, the phasic interface is eliminated from the formulation and the problem is made equivalent to a nonlinear heat conduction problem without phase change. Shamsunder and Sparrow (1975) and Shyy (1994a) summarized these features and demonstrated the usefulness of the enthalpy formulation.
17
The enthalpy formulation has been employed in many investigations, for example, Shamsunder and Sparrow (1975), Voller and Prakash (1987), Voller and Swaminathan (1991), Shyy (1994a, b), Shyy and Rao (1994), Shyy et al. (1996), Yao and Prusa (1989).
heated heat input
insulated
vapor / vapor insulated
liquid\ liquids
insulated cooled heat output
(a) Evaporation (b) Condensation
I 1 I
VT g Vp
Figure 3.1 Schematic diagram of the bulk evaporation and condensation
The purpose of the work in this chapter is to develop a consistent numerical method for the computation of the bulk evaporation and condensation processes involved in a specific class of liquidvapor phase change problems under constant volume and no
18
buoyancy conditions. To develop a generic numerical model and obtain the basic features of this phase change process, a simplified, one dimensional, conduction driven problem is particularly designed.
This problem can be depicted with the case presented in Fig.3.1. A rigid, cylindrical container is partially filled with vapor and liquid. Initially, the vapor and liquid are kept saturated at a certain pressure and uniform temperature. Well selected boundary conditions on the walls are imposed to activate the pure bulk evaporation or condensation, and also to avoid the buoyancy effect and bubble generation. While keeping the side wall insulated, only heating the top wall or cooling the bottom wall induces pure evaporation or condensation, respectively. The volumetric expansion or contraction due to the density change between two phases induces convection in the vapor phase. The scaling analysis conducted in this chapter indicates that this convection effect is quite week and negligible. It can be reasonably assumed that the phase change processes are driven entirely by the conduction and the release or absorption of latent heat at the phase change interface. Similar situations are also possible for the multidimensional phasechange problems under zero gravity condition.
No expansion work is involved in this constant volume phase change process. In the absence of mechanical work, the first law of thermodynamics indicates that the heat exchange between the system and the ambient equals the change of the internal energy of the system. The released or absorbed latent heat associated with the phasechange is equal to the change of internal energy. Based on the physical conditions of this system, a compact internal energy formulation of the energy equation is proposed. This formulation,
19
similar to the enthalpy formulation, is a single region formulation, where one set of governing equations can describe both phases.
3.1 The Internal Energy Formulation
The first law of thermodynamics, for any thermodynamic system, is
6Q=dE +6W=dH(pdV +Vdp)+6W (3.1) where E is the internal energy, H is enthalpy, Q is the heat input, W is the work output, p and V are the pressure and volume of the system, respectively. For the phase change process involved in the rigid and enclosed container with constant volume, there is no mechanical work. Therefore the first law for this closed and rigid system reduces to SQ = dE= dH  Vdp (3.2) which indicates that the heat exchange between the system and the ambient causes the internal energy, the enthalpy and pressure of the system to change. Therefore the released or absorbed latent heat associated with the phase change induces the changes of internal energy, enthalpy and pressure. The quantity of the heat exchange is equal to the internal energy change, but not equal to the enthalpy change. The conventional enthalpy formulation of the energy equation, commonly used for the phase change between solid and liquid, should include a transient pressure term for this phase change between liquid and vapor. The weak form of the energy equation based on enthalpy formulation is Dh Dp
p = +V(KVT) (3.3) DI Dt
20
This formulation including the transient pressure term introduces additional complexity and does not seem to be the best choice for this problem. A more attractive alternative, the internal energy based formulation for the energy equation which does not include the pressure term, is proposed and developed in this work. Similar to the enthalpy formulation used for melting and solidification problems, the internal energy formulation does not require to explicit tracking of the interface position.
The weak form of the energy equation with internal energy as the primary variable can be expressed as follows:
De
p = V (V T) (3.4) Dt
where the internal energy e may be expressed, in the form of a summation of sensible heat and latent heat, for the liquid, vapor and mixed phases as: e =c ,T f=0 e=c., T, + fAe O
21
PD (c, = V.(fVT)y f =0
Dt
D Df
Pf (c,,'T) = V (~,VT)pfAef O < <1 (3.6)
Dt
PIDt(cgT)= V(VT) =1
where the latent heat now appears as a source term. The form of this internal energy based equation is identical to that of the enthalpy based equation, except for the pressure term appearing in the latter one. Just like the enthalpy formulation, this internal energy formulation is a single region formulation, where one set of governing equations can describe both phases. A major advantage of the internal energy formulation is that it is not necessary to split the domain into separate subdomains consisting of different phases. A fixed grid can be employed to facilitate the computations.
Both enthalpy formulation and internal energy formulation are implemented and examined. The vapor phase fraction, f, may be defined in two different ways, which can lead to different update methods.
Some additional closure relations are required for the multiphase condition in the container. Property variations are based on a simplified searching routine. The different phases are differentiated based on the local temperature or internal energy at each nodal point, within the frame work of the temperature based or internal energy based update method. Using a temperature based update, the liquid or the vapor phase is identified by the local temperature, wherever the temperature at any point is lower than the saturated liquid temperature or higher than the saturated vapor temperature, that point is located in the liquid or vapor region, respectively. Similarly, if the internal energy based update is
22
used, the liquid or the vapor phase is identified by the local internal energy. Either case allows the vapor phase fraction in a computational cell to be calculated.
Using the ClausiusClapeyron equation, the saturation temperature as a function of the system pressure is calculated:
C
Tsat  (3.7) DInpsat
where C and D are constants.
The saturation pressure, or the vapor pressure, can be obtained from the approximate equation of state for an ideal gas, and the mass conservation in which the bulk vapor density and temperature are used:
p,, = p, = RT ,, mrotal ,= m+ mi = Constant (3.8) Vop
where R is the gas constant for the vapor phase. Eqs.(3.7) and (3.8), have been proposed and used as the basis for the model developed by the authors, Ding and Anghaie (1994a, b, 1995a, b, c), Anghaie and Ding (1995).
3.2 Analysis and Assumptions
To simplify and generalize the formulation for the bulk evaporation and condensation processes, these governing equations and closure relationships should be normalized. Scaling and analysis of order of magnitude are conducted to develop simplifying assumptions that lead to a concise and efficient solution procedure for the governing equations. For the different phases of any material, the thermodynamic and
23
transport properties are quite different, and they may be in the same order or may have a difference of several orders of magnitude. For the same phase, some properties are independent or slightly dependent on temperature or pressure, but others may be strongly dependent. The following relations are found about the orders and variations of these properties for the liquid and vapor phases far from the thermodynamic critical point, within a moderate phasechange range.
c, = constant, c,, _ constant, c,
~ 0(1),
C,
C, a constant, , _ constant, K
K = O(10)  O(102), (3.9)
p, a constant, p, > constant, P= (10) (1 0),
a,t constant, av constant, a= a, / a, O (10') In Eq.(3.9) ' O ' represents the order of magnitude.
Therefore all the properties, except the vapor density and vapor thermal diffusivity, can be assumed constant. The liquid specific heat and vapor specific heat are assumed equal. In the following analysis and computation, the values of the ratios such as S ' 10,p = = , = 10' are used, where the subscript 0 denotes a
K, Pv.o avo
referential state. Other values of these ratios can also be used. Computer experiments reveal that for the cases with higher property ratios, that is, p = 103, K = 102, the computation becomes more difficult to converge, especially when convection is also considered. It is due to the larger discontinuity encountered at the liquidvapor interface.
24
As discussed in some literature, for example, Crank (1984), Ozisik (1980), Ostrach (1983), due to the density difference between two different phases, the volumetric expansion effect induces convection in the phase of lower density. Let U, denote the speed for the motion of the interface and U, be that for the motion of the vapor phase resulting from the density change. The mass conservation at the interface should be satisfied as follows
pU, = p,(U,  U,) (3.10) The energybalance equation at the interface is
K, , +p,e,U, = (p,e,  p,e)U, (3.11) where e, and el denote the internal energies for the saturated vapor and liquid phases, respectively.
Eliminating U, from Eq.(3.11) via Eq.(3.10) yields
C ; = p,AeU, (3.12) which can also be rearranged as
Ui, = (, , )= St(  ) (3.13) p,Ae & & L dX Lp IX where L is the thermal diffusion length scale which is the length of the container, St=c,To/Ae, 0 = T/T X =x/L. p =p/p. There are two possible velocity scales at the interface, U = St(a/ll) and U2 = St(all)/p. Their ratio is, UI/U2 = K/cK, >1, thus U,>U2. It is known that the lower velocity scale U2 controls the movement of the interface.
25
Eq.(3.13) can also be expressed as
U, = U, U (3.14) dX dX
The velocity, U2 = S(a/L)/p, includes the density change effect which actually reveals the constraints of energy balance and mass conservation at the interface, and the volumetric expansion or contraction during evaporation or condensation, respectively. It is lower than the pure thermal diffusion velocity in liquid, alL and that in vapor aJL. Thus U2 is to be chosen as the characteristic velocity scale where the density change effect, as well as the St number, appears explicitly. The dimensionless velocity of the interface movement is
U = K U= (3.15) U, dx 6W
For the case of evaporation, for example, heat is supplied to the system from the top wall which is the boundary for the vapor phase. If the appropriate amount of heat flux is provided to ensure O(1), it is found that U ~ O(1), thus Ui  U, = St(a, / L) / p. Numerical results which will be shown in the next section to support this approximation.
The characteristic time scale is
L L2p
t,  (3.16) U, a St
Obviously this time scale is larger than those of thermal diffusion in vapor and L2 L2 L2
liquid phases, that is,  p >  > , which implies that the motion of the interface a, a, a,
26
controls the rate of this heat transfer process. From the above analysis, it can be concluded that the motion of the interface is affected by the density difference between liquid and vapor phases, thus its velocity is very low. During the evaporation process within a constant system volume, a density change increases the vapor pressure and the saturation temperature, which in turn suppresses the evaporation.
From Eqs.(3.10) and (3.15), the velocity of vapor convection induced by the 1pa
density change is U, = (1 p)Ui  St, and since p >>1, U, (c/L)St , which is p L
still very low. Therefore the Reynolds number and Peclet number for the convection are Re= U, 1, I/v=(1/L)(St/Pr), Pe= (IJL) St. , where 1, is the vapor motion length scale. Since the velocity of the vapor convection is very low, and the convection just occurs near the interface, the vapor convection length scale will be very small compared to the thermal diffusion length scale, say, 1,L < 0(1). For most vapors, Pr ~ 0(1), and if moderate St is considered, that is, St ~ 0(1), then Re < 0(1), and Pe < 0(1). For the problems in which Re and Pe are less than the order of unity, the convection effect is insignificant. Therefore heat conduction dominates the heat transfer process and the convection effect induced by the density change is negligible.
On the other hand, the convection effect caused by the density change during phase change can be much less significant than the convection induced by the buoyancy effect. The characteristic velocity of the natural convection is U, = (Gr)'12 v/L where Gr is the Grashof number, the corresponding Reynolds number Re= U,L/v = (Gr)"2. If moderate Grashof number is considered, for example, Gr = 104  106, then Re = 10 ~ I10, which is much stronger than the density change effect.
27
From this analysis, the convection effect caused by the density change can be neglected reasonably and the bulk evaporation and condensation processes can be assumed entirely driven by heat conduction in the absence of gravity. As a result, the governing equations (3.6) can be reduced to the the form of pure heat conduction equations:
p, (CtT) = V. (K'VT) f=o
P, (c.,f)= V.(KfVT) pfAe 0
p, (c T) = V  (rVT) f = 1
Eq.(3.17) can also be written in terms of dimensionless parameters as:
St o9
 aV2 f= pdr
St l f
= aV o< O <1 (3.18) pdr pcr
St 19
 V2 f=1 pOr
Here St is the Stefan number defined as cT,o/Ae, 0 = T/T,o where T,o is the initial saturation temperature. T,.o and the length of the container are chosen as the characteristic temperature scale and length scale, respectively. The characteristic time scale is shown in Eq.(3.16). The relationships presented in Eqs. (3.7) and (3.8) are nondimensionalized as
A
Os  AIn(s ) Ps = v(3.19) where p, and p, are dimensionless saturation pressure and vapor density, respectively.
28
3.3 The Vapor Phase Fraction Update
To solve the physical model described by Eq.(3.18), a closure relationship between the vapor phase fraction f and the temperature 0 or the internal energy e should be formulated. For a pure substance undergoing evaporation or condensation, the total internal energy is a discontinuous function of the temperature:
e
e =SO 0
(3.20)
e
e S,9+l > esat Ae Sat
where O,,, is the saturation temperature. However, from a computational viewpoint, discontinuities are difficult to track and it is often necessary to smear the phase change over a small temperature range to attain numerical stability. Thus we can have
e=S,O 0
e = S,  9, + (9 O,) / (0,  O,) 01 < 0 < ov (3.21)
J=S,0+1 Ov <
where 0=0,r, Ov=q.at+, and e is the phase change interval. As a result, the discontinuity is replaced by a small window, 2c, over which the phase change occurs. For a pure substance, this window is a purely numerical artifact and must be chosen as small as possible in order to accurately model the physical system. Substitution of Eq.(3.21) into the normalized form of Eq.(3.5) yields
f= 1 0 1 (3.22) 0,o0 2E
29
which is used to iteratively update the phase fraction from the computed temperature field. This technique has been successfully used by Shyy (1994a, b), Shyy et al. (1994a) and was referred to as the Tbased update method. This technique has proven to be very effective for achieving convergence. An alternative formulation to the Tbased update method is to express the phase fraction as a function of the total enthalpy rather than the temperature. This update is referred to as the Hbased update method. More information about this method can be found in Shyy (1994a) and Shyy et al. (1996).
A third formulation expressing the vapor phase fraction as a function of the total internal energy, rather than the temperature or enthalpy, is developed in this work. According to the classical thermodynamic definition, the relationship between the vapor phase fraction and the internal energy is
f =  e (3.23) This is used for the iterative update of the vapor phase fraction. This update procedure is referred to as the Ebased update method and is analogous to the isothermal case discussed by Voller and Swaminathan (1991). In particular, it is noted that =0O is allowable in this formulation, thereby allowing for modeling phase change of a pure substance with higher accuracy. This is the most important advantage of the Ebased method compared to the Tbased method where a significant effort must be made to determine the phase change temperature interval e. For the phase change problems involving complex interface geometry, the normal temperature gradient varies with position along the interface. The determination of the phase change intervals for such problems is a tedious task. The Ebased method does not need to use a phase change
30
window and can be conveniently used for many phase change problems including the ones with complicated geometry. Another advantage of the Ebased method is the fact that if the latent heat is large relative to the sensible internal energy, the vapor phase fraction varies slowly with respect to the total internal energy. This enhances the stability of the update procedure. For the liquidvapor phase change occurring far away from its critical point, the latent heat is relatively large and the application of the Ebased update method even for simple geometry is recommended. In the present work, the Tbased and Ebased methods are examined in detail from the viewpoint of accuracy and computational efficiency.
A straightforward implementation of the Tbased method is to discretize Eq. (3.19) for each computational cell:
appk = ka,,b k +bk  AV [fpk' _fp'] (3.24) pAr
where the subscript P represents the nodal point of the control volume to be solved, and nb represents the neiboring nodal points. This equation is solved according to the Tbased update procedure given in Eq.(3.22). As reported by Shyy and Rao (1994) and Prakash et al. (1987), for Stefan number less than unity (large values of latent heat), this form of the iterative update is prone to numerical instability in the form of oscillations. For moderate values of Stefan numbers, the use of an underrelaxation approach for the calculation of the source term facilitates the solution process. The main cause of the oscillatory, nonconvergent behavior is the extreme sensitivity of the vapor fraction to small changes in temperature which leads to a large negative feedback effect through the source term. Substituting Eq.(3.22) into the discrete form of the conservation law and moving the
31
resulting coefficients of , into the ap term in Eq.(3.24), yields a set of stable equations similar to those formulated by Shyy and Rao (1994). These equations are: AV
aBpk = kaB, +bk  p[f" _ k,k] < O or pk> 9
AV 1 AV 0< (
[ap + ]Op = kab + bk+  t +fpn1] < ,Ov (3.25)
pAr 26 P +b + c
where the superscript n denotes the current time step and k denotes the current iterative value. Using this form of the update procedure, the temperature and the vapor phase fraction fields can be computed simultaneously and are therefore coupled at every step of the iterative process. This numerical procedure converges quickly with no loss of accuracy, granted that a proper value of E must be used.
The implementation of the Ebased method is relatively straightforward and does not need any preconditioning. The internal energy is updated first, with the current values of saturation temperature and latent heat, Fp = SOpk +f 0 < fk < 1 (3.26) Then the vapor phase fraction is updated with the current internal energy
fk = k ,
3.4 Properties of the Control Volume Containing Two Phases
The control volume which contains two phases is shown in Fig.3.2. Assume the vapor phase fraction of this control volume is f The vapor phase occupies a volume of
32
fAxAy, while the liquid occupies (1]) AxAy. For this one dimensional case, the interface is normal to the x direction. The effective conductivity of this control volume can be obtained through a heat flux analysis. Assuming the control volume is at a steady state, that is, the interface does not move and the latent heat is not involved, the heat flux along the x axial direction is conservative AT AT A T
q, =   =K = K, (3.28) S'A fAx ( f)Ax where Ic is the effective conductivity, and the temperature is continuous
AT = T  T, = (T,  T,) + (T.,  T) = AT, + AT (3.29) Hence
IKC f 1f (3.30) K, KI
which is in a harmonic average form.
According to the thermodynamic definition of the specific volume per unit mass for a mixture system, it is easy and straightforward to determine the effective density of this control volume as
1 1 1
P (3.31) v1 vf +(1f)v, f 1 f Pv P1
where v is the specific volume per unit mass. It is similar to Eq.(3.30).
33
fAX (1Axf) Tt T
T vapor liquid k,
kv Ay
inter ace
Ax
Figure 3.2 The control volume containing two phases
3.5 Results and Discussion
3.5.1 Accuracy Aspect
A test case of pure bulk evaporation is selected to examine the accuracy, and computational performance of each formulation with different update methods. A fixed volume container with one fourth of volume is filled with saturated liquid is chosen. The rest of the container's volume is filled with the saturated vapor. The initial temperature of the saturated fluid is 1.0. When the top wall is heated and kept at a temperature of 1.5,
34
while the bottom wall is insulated. Such an orientation of heat flow is required to avoid the superheat in liquid, and the bubble generation and boiling, while the bulk evaporation is allowed. A mesh of 8 1x6 nonuniform and fixed grids is adopted for the liquid and vapor within the container. The time step size Ar = 104 is chosen. Second order central differences are used for all spatial derivative terms and fully implicit differencing is employed for time marching.
Figs.3.3(a) and (b) show the evolution of the evaporation process. At the time instant of 0.03, about 18% of liquid is vaporized, and the saturation temperature is about 1.38. A little wavy or nonuniform variation can be seen due to the nature of the discretization and the nonuniform release or absorption of latent heat. A significant part of the latent heat content in a computational cell is absorbed when it starts to undergo evaporation. Thus the interface movement is slowed down when phase change in this cell starts and subsequently speeds up. The grids are clustered near the liquidvapor interface, and the finer grids are placed near the interface, covering the range of interface movement. Different minimum grid sizes, such as Ax,,,, = 1/120, 1/200, 1/400, are selected to examine the solution accuracy and the independence of solution on the grid size. The results with different minimum grid sizes are basically identical, while the finer grid size yields a smoother solution. In the following computations, the grid size with AX,,, = 1/200 is used.
35
0.18
0
0.16 8 0.14 S0.12 S0.1 . 0.08 0.06 0.04
0.02
0
0 0.005 0.01 0.015 0.02 0.025 0.03 Time (t)
1.4
1.35
1.3
0 1.25 S1.2 S1.15
S1.1
1.05
0 0.005 0.01 0.015 0.02 0.025 0.03 Time (7)
Figure 3.3 Solution features for bulk evaporation
(Grids 81x8; solid line: AXi=1/120; dashed line: AX,=1/200;
dashed and solid line: AX,=1/400)
36
0.16
0.14
0.12
S 0.1
0.08 0.06
0.02
0 0.005 0.01 0.015 0.02 0.025 0.03 Time( )
(a) evolution of vaporization
1.4
135
1.3
1.25
S1.15 " 1.1
1.05
0 0.005 0.01 0.015 0.02 0.025 0.03 Time ( t )
(b) variation of saturation temperature
1,100
1,000 8 900
700
0 500
400
S300
200
100
0 0.005 0.01 0.015 0.02 0.025 0.03 Time( )
(c) computation efficiency performance
Figure 3.4 Comparison of two formulations for evaporation
(The solid line: the internal energy formulation; the dashed line: the enthalpy formulation)
37
0.1
0.16
0.14 0.12
0.1
g 0.08
0.06
0.02
0
0 0.005 0,01 0.015 0.02 0.025 0.03 Time ( T )
(a) evolution of vaporization
1.4
1 1.35
1.3 1.25
1.15
1.05
0.005 0.01 0.015 0.02 0.025 0.03 Time ( s )
(b) variation of saturation temperature
1,400 1,200
' 1,000
600
Z 200
0 0.005 0.01 0.015 0.02 0.025 0.03 Time ( T )
(c) computation efficiency performance
Figure 3.5 Comparison of two update methods for evaporation
(The solid line: the Ebased method; the dashed line: the Tbased method)
38
3.5.2 Comparison of H Formulation and E Formulation
The same test case is computed to examine the performance of different formulations, the enthalpy formulation and internal energy formulation. Figs.3.4(a) and (b) show the evaporated liquid fraction and the variation of the saturation temperature during the evaporation process. Fig.3.4(c) shows the efficiency aspect of the computations, that is, the number of iterations required at each time step for the residual to go below 10 4. The enthalpy formulation requires a few more iterations to converge. It is found that the internal energy formulation and enthalpy formulation yield nearly identical results, with comparable computational efficiency. It indicates that both of these formulations are satisfactory for this problem and have similar performance efficiency. The only advantage in the internal energy formulation is the more concise form without an explicit pressure term. In all further investigations, only the internal energy formulation is adopted.
3.5.3 Comparison of EBased and TBased Update Methods
The Tbased update method and the Ebased update method are examined for the same testing case. Figures 3.5(a) and 3.5(b) show that the Tbased method and Ebased method yield nearly identical results. It is seen that convergence for the Ebased method requires about half the iterations which is needed for the Tbased method. Thus the Ebased method is computationally more efficient than the Tbased method, as shown in Fig.3.5(c). The performance efficiency of these two update methods for the evaporation
39
and condensation problems are different from those for melting and solidification problems, where the Tbased is better than the Hbased method as reported by Shyy and Rao (1994).
The reason for the improved efficiency of the Ebased method is that the phase change temperature window is used for the Tbased method, which is not required for the Ebased method. When the Tbased method is used and a purely artificial phase change window has to be adopted, it is not easy to use the temperature to identify the phase interface since the temperature is always continuous at the interface. For melting and solidification problems, the phase change temperature is usually close to a constant. The saturation temperature, however, changes continuously for evaporation and condensation processes in a constant volume system. This difference causes distinctive computational efficiency performances. For each iteration at every time step, the saturated vapor and liquid temperature as well as the saturation temperature should be updated. Therefore, the Tbased method converges more slowly. When the Ebased method is used, the internal energy has a sharp jump at the interface, and it is easy to use internal energy to identify the phases. As has been mentioned before, the relatively large amount of latent heat will result in a more stable update procedure and faster convergence rate. Furthermore, as discussed in the section of update methods, if the Tbased method is used, there will be some difficulty to determine the phase change temperature intervals for the complicated phase change problems where the interface could be of any arbitrary shape. Use of the Ebased method can completely circumvent the need for a phase change window. Therefore, the Ebased method is superior to the Tbased method, in terms of computational efficiency and
40
formulation simplicity; also it is capable of handling more complicated phase change geometry. Based on this finding,the Ebased update method is used to study the remaining cases.
3.5.4 Bulk Evaporation and Condensation
In the following work, the internal energy formulation and the Ebased update method are used to investigate the bulk evaporation and condensation processes and also to analyze the interface motion. With other conditions remaining the same, such as the initial temperature, the grid distribution and the time step size, different boundary conditions are imposed at the top and bottom walls. First, the container's top wall is heated with a constant nondimensional heat flux q, = ~'x0= 1.0, and the bottom wall BX
is insulated. After the time instant of 0.02, the top wall is insulated while the bottom wall is cooled with a constant heat flux q, = ic!x=.= 10. Correspondingly, the system undergoes an evaporationcondensation process. Fig.3.6(a) shows the variation of the saturation temperature during the transient. The motion of the interface, which is respectively downward and upward during evaporation and condensation, is exhibited in Fig.3.6(b). The axial temperature distributions at different time instants during the evaporation process are plotted in Fig.3.6(c), and those of condensation process are plotted in Fig.3.6(d).
41
The velocity of the interface motion can be calculated from results presented in Fig.3.6(b). During the evaporation and condensation processes, the interface velocity is always constant because only the constant heat flux at either the top wall or bottom wall is provided. During the evaporation process which takes place with a constant heat flux at the top wall, the interface velocity is a constant about 2.5, which is in the order of 1.0. While during the condensation process with a constant heat flux at the bottom wall, the interface velocity still remains constant at the level of about 0.8. This indicates clearly that the nondimensional interface velocity is consistently of the order of 1.0. In other words, it can be concluded that U, ~ 0(1), thus Ui ~ Ur = S,(a, / ,) / p. This analysis verifies the phenomenological analysis presented in previous part of this paper. This also indicates that the velocity scale Ur = S, (a, / l,) /p has been chosen appropriately so that the velocity of the interface motion can be calculated directly.
3.6 Concluding Remarks
A numerical method is developed to model the bulk evaporation and condensation processes in an encapsulated container. An internal energy formulation for the constant volume phase change process is proposed and successfully implemented. Some assumptions regarding the convection heat transfer rate at the vaporliquid interface are made through scaling and order of magnitude analysis. The velocity and time scales of the interface movement are obtained through scaling analysis. Depending on the configuration and boundary conditions, the convection effect induced by the density change is analyzed
42
and found negligible. The dominating mechanism for the heat transfer at the interface is conduction.
1.5
S1.45 1.4
S1.35
a 1.3 3 1.2
o 1.15
1.05
1
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07
Time (r )
(a) Variation of saturation temperature during evaporation and condensation
0.82
0.81 0.8
0.79
S0.78
0.77 0.76
0.75 I I
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Time (c)
(b) Motion of interface during evaporation and condensation
Figure 3.6 Bulk evaporation and condensation
43
2.1
0=0.01
1.9   0
1.7 C=0.02
1.2
1.1
0 0.10.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Axial position X
(c) Temperature distribution during evaporation
1.8
1.7  0 ,=.03
 1.6 T=0,O4
1.5
1.4 . .0
1.3 t06 05
1 =0.06 T=0.05
S1.2
1.1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Axial position X
(d) Temperature distribution during condensation
Figure 3.6 ( continued) Bulk evaporation and condensation
44
The computational performances of the wellestablished enthalpy formulation and the proposed internal energy formulation are evaluated and compared. Both formulations yield identical results with similar computational efficiencies. The internal energy formulation has the advantage of having a more compact form. The performances of two update methods for the vapor phase fraction, the Ebased and Tbased methods, are investigated. The Ebased method is superior to the Tbased method in terms of computational efficiency and formulation simplicity for evaporation and condensation problems. The Ebased method is also capable of handling more complicated phase change problems, because it does not require the artificial phase change window. The internal energy formulation and the Ebased method are successfully implemented and used to compute bulk evaporation and condensation processes under different conditions. The evolution of the bulk evaporation and condensation processes is simulated. The variations of the temperature distribution, the liquid content and the saturation temperature with time for a variety of boundary conditions are obtained.
Further effort will be highly worthwhile and strongly recommended to extend the current model for liquidvapor phase change processes to a unified threephase model. All phase change processes between these three phases, that is, solidliquid, liquidvapor, solidvapor, or solidliquidvapor, could be simulated.
CHAPTER 4
BULK EVAPORATION/CONDENSATION WITH INTERNAL HEATING AT ZERO GRAVITY
Based on the modeling procedure developed in the previous chapter, the bulk evaporation/condensation processes with internal heat generation at zero gravity level are studied. The system undergoing the phase change process and at a steady state is shown as Figure 4.1. Since the system is designed to operate under a 0g environment, there are no buoyancy effect and no natural convection. Initially, the vapor and liquid are kept saturated at a certain pressure. The top wall and bottom wall are insulated. The side wall is cooled with a constant heat flux. The internally generated heat evaporates the liquid, while the cold side wall removes the heat from the liquidphase and condenses the vaporphase. For this system with such an orientation and operated under 0g, there is no bubble or droplet generation on the wall surface. The temperature along the liquidvapor interface does not change significantly, because the curvature of the interface is moderate. The liquidvapor interface can be reasonably assumed to be isothermal, consequently the thermal capillary effect can be neglected. Hence the phasechange process is entirely driven by the internal heat generation and heat conduction. The rate of evaporation or condensation is determined by the local temperature distribution. The time dependent and nonuniform internal heat generation results in a transient temperature distribution and causes the motion and deformation of the interface.
45
46
insulated
;vapor condenIation va or cooled
evapqration
liquid Qgen,l
liquid insulated
(a) schematic diagram of the phasechange (b) at the final steady state
process
Figure 4.1 Schematic diagram of the twophase cell operated at Og
47
4.1 The Governing Equations
The internal energy formulation of the energy equation for this phasechange process with the constraint of constant volume is used. The weak form of the energy equation for conduction is
p = V (VT) +Q, (4.1) The source term Qg,, is the internal heat generation rate which varies with time and space according to the following relation:
Q , = pB(t)  sin(L) (4.2) where B is a time dependent parameter, p is the local phase density, L is the length of the container.
The internal energy e for different phases is also given by Eq.(3.5). Substituting Eq.(3.5) into the governing Eq.(4.1) yields
p(c ,,'T) = V(KjVT)+ Qg,,
p (c,P)= V(K,VT)pAe + Q. (4.3)
p (c,T)= V ( (VT) + Q,,
48
where the latent heat now appears as a source term. The vapor phase fraction is defined as
f = e e (4.4) e,  e,
This is the Ebased update method discussed in chapter 3 and used for the iterative update of the vapor fraction. The Tbased method does not work well for this twodimensional bulk evaporation and condensation problem, while the Ebased always works quite well. Some additional closure relations, the same as Eqs.(3.78) and (3.19), are required for the multiphase condition in the container.
Some assumptions and simplifications have to be made based on the order and scaling analysis as discussed in chapter 3. All the properties, except the vapor density and vapor thermal diffusivity, are assumed constant. To simplify the computation, the liquid specific heat and vapor specific heat are assumed equal. In the following computation, these values of the ratios such as =   10, Pj   = 100 are used, where the K, P,.o
subscript 0 denotes a referential state.
The governing equation (4.4) is normalized as:
P 1 V2 +g f=O
dr
p = V2 e pf O
P =St Or
p,,=V 6+Q, 0
where St is the Stefan number defined as c,T,oAe, The nondimensional heat generation rate is
49
..,, = AG(r)sin(irX) (4.6) where G(t) = B(t)(tr/c,T,). The initial saturation temperature Ts,o and the length of the container are chosen as the characteristic temperature and length scales, respectively. The characteristic time scale is tr = L2/a which is the thermal diffusion time scale in the vapor phase.
4.2 Conductivity at Control Volume Surfaces
The finite control volume difference method is employed to discretize the differential equations. For different materials, that is, the solid wall and the fluids in the container, their thermal and transport properties such as densities, heat conductivities, are quite different. Due to the nature of this discretization, the interface between two different media is always located on the control volume surface. This arrangement requires that the effective properties at the interface, namely the control volume surface, be determined. Patankar (1980) discussed this problem and presented the determination of harmonic mean effective conductivity at the interface for the rectangular coordinate system.
Here in this work, the general determination of the conductivity at the surfaces of any control volume whether it contacts to an interface or not, for a cylindrical coordinate system, is developed through the heat flux analysis. As shown in Figure 4.2, P is the nodal point of the control volume, where the value off at P is to be calculated. E, W, N and S are nodal points surrounding the control volume, and e, w, n and s represent the surfaces enclosing this control volume.
50
\IN
/   e
r
Figure 4.2 The control volume contacting other material
Considering the two dimensional, axisymmetrical coordinates, the linear heat flux along radial direction is
q, = 2r  = 2r (T  T) (4.7) c+ ln(rN / rp) and also
qr = 2 (T.  T) 2 p (T  T) (4.8) q, =2tKN =27r ) (4.8) In(r,/r.) In(r./rp)
Using the relation
TTN = T T + T T (4.9) the heat conductivity at the north surface, ,, is
51
KP = IrN (4.10) cp In(rN / r,) + Nln(r, / r,)
The determination of heat conductivity at the east surface is simpler and more straightforward (see Patankar, 1980), =(1f, f (4.11) KP, KE
wheref, is a ratio defined
f, (&), (4.12)
(&)V+ + (&),
in terms of the distances shown in Figure 4.2. The heat conductivity at the south and west surfaces could be obtained similarly. This approach used to determine the transport properties at the control volume surfaces provides consistent relations for properties that can be applied for any interfacial boundary, and results in higher accuracy with low computation cost.
4.3. Results and Discussion
The second order central difference scheme for the spatial differential terms and the fully implicit marching scheme for the temporal term are adopted to solve Eq. (4.5). A uniform mesh of 101x41 grids covers the liquid phase and vapor phase within the container, and another non uniform mesh with 101x5 grids covers the side wall region. The outer surface of the side wall is cooled by a coolant flow at constant flux which is negative unity.
52
Because of symmetry, only the right half of the container is considered. Initially, the lower half of the container is filled with the liquid phase and the upper half is occupied by the vapor phase. The aspect ratio of the container is 2.5. The initial saturated vapor and liquid temperature is 1. At r = 0, the side wall is cooled at a constant heat flux, and the internal heat generation begins. The time dependent parameter for the heat generation rate expressed in Eq.(4.7) is given by Eq.(4.13) and shown in Figure 4.3.
G(r) = 12r 0< r< 1.5
G(r) = 12x 1.5 = 18 r > 1.5
G (r)
Figure 4.3 Parameter for heat generation rate
The evolution of the evaporation/condensation process is shown in Figure 4.4. Using the phase identifier which is the local vapor phase fraction, the distribution of liquid and vapor phases and the position of the interface are determined. The topology of the liquidvapor interface and the temperature distribution at several time instants are
53
displayed in Figure 4.4 and Figure 4.5, respectively. The internal heat generation evaporates the liquid, while the side wall cooling tends to condense the vapor. Globally, the evaporation of liquid and condensation of vapor take place simultaneously, as shown schematically in Figure 4.1. The local phasechange process, whether evaporation or condensation, is dominated by the local internal heating or by the wall cooling effect, since the internal heat generation is a spatial function. The evaporation is dominant near the central axial position where the heat generation intensity is higher and evaporates more liquid. The condensation caused by the side wall cooling is dominant near the axial end and side wall since the heat generation density is lower there. The simultaneous evaporation and condensation cause motion and deformation of the interface. Because of the nature of time dependent heat generation, the local evaporation and condensation, the motion and deformation of interface continue.
As the process evolves, a certain amount of liquid is evaporated into vapor due to the internal heat generation which increases at the early time period, as Eq.(4.13) indicates. Meanwhile, some amount of vapor near the side wall is condensed into liquid due to the cooling of the side wall. At a very earlier time period, a samll portion of liquid film is formed along the side wall owing to the condensation, and the liquidvapor interface becomes concave because of the evaporation. With increasing time, the interface keeps moving and deforming. The liquid film grows longer and thicker, and the position of the interface at the centerline of the container descends lower and lower. When r = 1.0, the liquid film grows up to the top end and covers the entire side wall. The liquid phase which is initially filled in the lower part of the container now forms a long film covering
54
and residing along the entire side wall. The liquid film becomes thicker and thicker at the top end where the internal heat generation is zero. At r =1.5, the interface at the low central side descends quite close to the bottom end. After r =1.5, the internal heat generation rate is a constant as indicated by Eq. (4.13) and the cooling heat flux at the side wall is still kept constant, the system reaches a steady state and the interface no longer moves. The topology of the interface strongly depends on the temperature distribution and the internal heat generation rate distribution which now is independent of time but spatially dependent. The shape of the interface approximately follows the heat generation distribution in the form of chopped sine function.
During the transient process, the temperature distribution varies with time. The temperature contours for the entire system gradually become denser and denser, since the internal heat generation is increasing. After steady state is reached, the axial temperature distribution follows approximately the form of chopped sine function, due to the nature of the heat generation distribution as shown in Eq. (4.2). There is a temperature gradient jump at the interface due to the release or absorption of latent heat associated with phase change. The discontinuity of the temperature gradient can be observed more clearly in a three dimensional temperature distribution displayed in Figure 4.6. There is an obvious fold on the temperature surf where the liquidvapor interface is exactly located. The axial temperature profile on the side wall is also in an approximate form of chopped sine function, but much flatter than that in the vapor region. The vapor is well superheated and its temperature gradient is greater than that of the liquid, because the liquid is a much better conductor than the vapor. This isothermalization of the liquid near the side wall is
55
expected to satisfy the design feature of the multiphase fuel cell to be possibly applied in space. The formation of the liquid film and its full coverage on the side wall can prevent the side wall from being overheated, which is the expectation for the application of such multiphase fuel cell.
The transient characteristic of this system can also be seen from the variation of the saturation temperature with time shown in Figure 4.7. The saturation temperature decreases a little bit in the very early time period when the internal heating is very weak and the side wall cooling and condensation are dominant. As the internal heating rate increases, the evaporation becomes dominant and the saturation temperature rises. After r = 1.5, the system gradually reaches steady state and the evaporation and condensation are balanced and the saturation temperature maintains constant.
The operation in a zero gravity environment is ideal and unrealistic. However, the study conducted in this chapter provides a basis for further exploriation of the fundamental performance and features of such a multiphase and internally heated system. The present results will be used as reference for comparison with the performance and features of this system operated under microgravity conditions, which will be presented in chapter 6.
56
vapor
liquid
= 0.2 0.4 0.6 0.8
r= 1.0 1.2 1.4 1.6
Figure 4.4 The evolution of liquidvapor interface (Og)
57
0.2 0.4 0.6 0.8
= 1.0 1.2 1.4 1.6
Figure 4.5 The evolution of temperature field (0g)
58
Figure 4.6 Temperature surf of the system at steady state (Og)
1.6
1.5 1.4
1.3
1.2
.2 1.1
0.9
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Time (x)
Figure 4.7 Variation of saturation temperature (0g)
CHAPTER 5
NATURAL CONVECTION WITH INTERNAL HEATING IN SINGLEPHASE AND TWOPHASE SYSTEMS
Like all phase change processes in a finite gravitational field, evaporation and condensation are inevitably subject to the effect of convection induced by buoyancy. The existence of convection can substantially affect the temperature field, and subsequently the performance of the phase change process. Internal heating, like external heating, results in temperature and density gradients, and results in a buoyancy effect in nonzero gravity. Such convection induced by internal heating is defined as internal natural convection which wull be investigated. Natural convection is involved in both liquid and gas phases. The coexistence of natural convection in both phases causes interaction of the convection cells between two phases.
In this chapter, some natural convection problems with internal heating in the singlephase and twophase systems are investigated. A benchmark problem, natural convection in a square cavity of singlephase without internal heating, is computed. Comparisons of the present solutions with accurate published benchmark solutions are made to verify the computational accuracy. The natural convection in a square cavity of singlephase with internal heating is also computed. Qualitative comparisons with the published experimental results are also made. Then natural convection in a liquidgas system without phasechange, with or without internal heating, is investigated. The interactions between these convection cells in both phases are revealed.
59
60
5.1. External Natural Convection in a SquareCavity Benchmark Problem
Considerable attention has been given to the study of natural convection in enclosures. The temperature gradients supplied externally across the system induce the density gradients and consequently the buoyancydriven flow. A great deal of research, both numerical and experimental, on this classical type of convective heat transfer process has been conducted, Ostrach (1982, 1983, 1988), Gebhart et al. (1988), de Vahl Davis (1983), Hortmann et al. (1990), Shyy et al. (1992b), Shyy (1994a), Bejan (1984), Carpenter et al. (1989), Ramaswamy et al. (1992), Bergman et al. (1981), etc.. Ostrach (1982, 1983, 1988), Gebhart et al. (1988) and Bejan (1984) reviewed and summarized those activities and achievements. Among the large number of natural convection problems, buoyancydriven flow in a square cavity with vertical side walls heated differentially, a socalled benchmark problem for natural convection, is a suitable vehicle for testing and validating computer codes. Accurate benchmark solutions have been published by de Vahl Davis (1983) and Hortmann et al. (1990) and are available for comparison with the present computations.
The problem being considered is that of twodimensional flow of a Boussinesq fluid of Prandtl number 0.71 in an upright square cavity. The horizontal walls are insulated, and the vertical sides are at temperature Th and T,. Figure 5.1 shows the domain and boundary conditions of the problem in schematic form.
61
= 0,u=v=0
T=T, T=T
u=v=O g u=v=0
S=O,u=v=O
Figure 5.1 Domain and boundary conditions for the benchmark problem The governing equations in dimensionless form are
+ = 0 (5.1) U+V = + Pr VU (5.2)
U+V  = +Pr V2V + RaPr (5.3)
d di'di'
62
U +V = V29 (5.4) OIX IY
where the external Rayleigh number is defined as RaE T) (5.5) au
and f is the thermal expansion coefficient which is used in the Boussinesq approximation p= P0[l f(T To)] (5.6) In obtaining the above equations, the following dimensionless variables are used
X = x/L, Y = y/L, U = u/(a/L), V= v/(a/L), 9 = (TTd)/(TT) (5.7) where (a/L) is the thermal diffusion velocity scale. Other velocity scales, such as buoyancy velocity scale, etc., can also be used to yield different forms of the governing equations. Solutions have been obtained using a secondorder centraldifference scheme for the momentum and energy equations. Two grid sizes, 41x41 and 81x81, are adopted for cases ofRaE = 103 and RaE = 106 , respectively.
Two external Rayleigh numbers, 103 and 106, are chosen for the calculations. Figure 5.2 shows the steady state computed results. An important feature of the solution is the skewsymmetry of the velocity and the temperature field with respect to the center of the cavity. Fluid is heated at the left wall, decreases in density and rises up, whereas fluid adjacent to the cold right wall loses heat and increases in density and falls down. This establishes a clockwise recirculating flow inside the box. With increasing Rayleigh number, it can be observed that the convection strength increases and the isotherms are nearly horizontal in the core flow. The heat transfer rate is observed to be enhanced with increasing Ra as can be deduced from the isotherm spacing at the left and right walls.
63
Increasing the Ra increases the vorticity generation in the flow field which in turn
increases the convection strength and results in an enhanced overall heat transfer rate
through the domain. The enhancement of the heat transfer rate with the increasing Ra can
be observed more clearly in Figure 5.3 which shows the Nusselt number distribution on
the hot and cold walls. The Nusselt number on the wall is defined as follows
hL c0
Nu = I = _ (5.8) K OI
which actually is the nondimensional temperature gradient or heat flux at the wall.
The present results are compared with the benchmark solution of de Vahl Davis
(1983). This comparison is shown in Table 5.1 for RaE = 103, 106. As may be observed
from this table, the present solutions agree quite well with the benchmark solutions.
Table 5.1 Comparison of present results with benchmark solutions of de Vahl Davis (1983)
Ra=10 (41x4Grid) Ra=106 81x81Grid)
Davis present discre Davis present discre(1983) pancy (1983) pancy
Y1,,0 1.174 1.167 0.5% 16.32 16.08 1.5%
V/ma    16.75 16.57 1%
XY 0.151, 0.543 0.151, 0.540
Um, 3.634 3.629 0.1% 64.63 64.3 0.5%
Y 0.813 0.813 0.855 0.852
V,,, 3.679 3.674 0.1% 219.36 219.1 0.1%
X 0.178 0.178 0.0379 0.0375
Nuo., 1.505 1.515 0.6% 17.92 17.38 2.7% Y 0.092 0.088 0.0378 0.0375
Nuo,,,, 0.692 0.693 0.1% 0.989 1.023 3.3%
Y 1 I 1 1
64
(a) Ra = 10, Pr =0.71 (41x41 Grid) Av=0.1 AO= 0.05
(b) Ra= 10, Pr = 0.71 (81x81 Grid) AY= 1 AO= 0.05
Figure 5.2 Flow patterns (left) and isotherms (right) of buoyancy driven natural convection
65
0.9  hot wall
0.8 cold wall
0.7
0.6
0.5
0.4 0.3
0.2 0.1
0
0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6
Nusselt number distribution on hot and cold walls
(External natural convection with Ra=103, Pr0.71)
0.9  hot wall
0.8
0.7  cold
0.6
0.5
0.4
0.3
0.2
0.1
0 I I
0 2 4 6 8 10 12 14 16 18 Nu (Ra=106)
Figure 5.3 Nusselt number distribution on side walls
66
An examination of the accuracy and dependence of the solution on the grid size is performed. The temperature profile at the horizontal centerline of the Ra = 106 case is chosen to be tested with different grid sizes, for example, AX = 1/20, 1/40, 1/80 and 1/120. Figure 5.4 shows the numerical solutions with AX = 1/40, 1/80 and 1/120 are almost identical and the solution with Ad = 1/20 is not quite accurate. With AX < 1/40, the solutions are spatially accurate and independent of grid size for the benchmark problem with Ra = 106.
0.8
0.6
0.4
0.2
0I
0 .1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x
Figure 5.4 Temperature profile at horizontal centerline
(Ra = 106 ; Solutions of different grid sizes; solid lines : 121x121 and 81x81, dots: 41x41, dashed line  : 21x21)
67
5.2 Natural Convection in a SquareCavity with Internal Heating
Natural convection induced by external heating is of importance in many applications. Natural convection in an enclosure in which the internal heat generation is present is also of prime importance in certain technological applications. Examples are postaccident heat removal in nuclear reactors and some other geophysical or astrophysical problems. Natural convection induced simultaneously by external heating and internal energy sources, however, in contrast to the extensive investigations on the conventional external heating convection problems, has received very limited attention, for example, Archarya and Goldstein (1985), Lee and Goldstein (1988), Fusegi et al. (1992). The presence of internal heat generation provides an additional dynamic in overall convective flow systems and offers a complicated yet challenging aspect of the convection under consideration. The introduction of internal heat generation brings forth an additional nondimensional parameter Raj, which was conveniently termed as internal Rayleigh number in Archarya and Goldstein (1985), Lee and Goldstein (1988), Fusegi et al. (1992). This is in addition to the conventional Rayleigh number, RaE, which is based on the externally imposed temperature difference.
The work described in this section is to explore the effect and significance of the internal heat generation induced convection and qualitatively verify the results with previous investigations such as Archarya and Goldstein (1985), Lee and Goldstein (1988), Fusegi et al. (1992). In addition, further details of flow and heat transfer characteristics
68
will be displayed. The physical situations, the domain and boundary conditions, to be considered are shown in Figure 5.5. Two boundary conditions on temperature are considered. One is that all the walls are kept at temperature Tw. Another one is that the top and bottom walls are insulated while the side walls are maintained at Th and T The continuity and momentum equations are the same as Eqs.(5.1)(5.3). The only difference is the appearance of an additional term containing Rai in the energy equation as shown here,
U +V = V+ Ra, (5.9) a aY RaE
where RaE, the external Rayleigh number, is still the same as defined in Eq.(5.5) and the internal Rayleigh number, Raj, is defined as
Ra, = gOGL (5.10) avic
where G is the volumetric internal heat generation rate.
The results are shown in Figures 5.65.8. For case I where all walls are maintained isothermal at T,, the natural convection is purely induced by the internal heat generation. The flow field and temperature field are all symmetric about the vertical centerline due to the uniform temperature along all walls and the uniform internal heat generation. In the presence of internal heating, a pair of counterrotation convective rolls is formed. The flows move upward from the central lower region, divides and move downward separately and symmetrically along both side walls. The convection induced by the internal heating shifts the isotherms toward the top wall so that the position of the maximum temperature is higher than that of pure conduction. When the internal heating is not strong, that is, Raj
69
< RaE, the convection effect is not significant. The isotherms are closely symmetric about the center of the cavity. As Raj increases, the position of T_,, moves higher. When Ra is increased to 106, there are two extreme values of temperature in the fluid and two extreme values of temperature gradient on the top wall. The flow pattern and temperature distribution agree with the experimental results of Lee and Goldstein (1988) qualitatively quite well. As the Ra, increases, the flow velocity increases, and the strength of the convection is increased. The heat transfer is increased significantly as the internal heating increases, which is displayed in Figure 5.8.
In case 2 the characteristics of natural convection, as shown in Figure 5.7, are quite different since both external and internal heating exist. In the presence of internal heating (Raj > 0), the aforementioned trends in flow and heat transfer persist if Rai < RaE. However, when the magnitude of Raj considerably exceeds that of RaE( e.g., RaE = 10, Raj = 10S) the flow field takes on a different pattern (Figure 5.7 (a)). It now consists of two counterrotating convective rolls spanning the entire length normal to the adiabatic boundaries. The maximum dimensionless temperature occurs between the hot and cold walls which explains why the fluid rises in the central region of the cavity and sinks along both the hot and cold walls. With fixed Raj, as the RaE is increased to 104, the flow field varies significantly. Three eddies are present, as shown in Figure 5.7 (b). The largest eddy is the counterclockwise flow moving downward along the cold wall and upward along the lower portion of the hot surface. The other two eddies are formed at the topright comer and bottomright comer of the box with the flow directed downward along the hot wall. The eddy at the topright corner is faster than the one at bottomright comer although
70
they are both sluggish relative to the largest eddy. The flow pattern and temperature distribution agree qualitatively very well with the numerical results of Archarya and Goldstein (1985) under the same conditions.
T= T,,u= v= 0
T= Tg T=T u=v=0 u=v=0
T= T,,u=v=O
case 1
(a) all walls at Tw
 =0,u = v = 0
T=T g T=Th u=v=0 u=v=0
=0,u = V = 0
case 2
(b) top and bottom walls insulated, left wall at T, and right wall at Th
Figure 5.5 Domain and boundary conditions for internal heating induced convection
71
(a) Rai = 104 (left: A V= 0.04; right: A9= 0.01)
(b) Ra; = 10 (left: A= 0.4; right: AO= 0.07)
(c) Rai = 106 (left: A = 1; right: AO= 0.5) Figure 5.6 Flow pattern (left) and isotherms (right) of internal heating induced convection ( case 1: all walls are at T~, RaE = 104, Pr = 0.71)
72
(a) RaE = 10'
AV= 0.4 AO= 1
(b) RaE = 104
Ay/= 0.5 AO= 0.1
Figure 5.7 Flow pattern (left) and isotherms (right) of internal heating induced convection
( case 2: left wall at 0 =0, right wall at h =1, Rai = 105, Pr = 0.71)
73
1 l
0.9  Ra l
0.8
0.7  Ra,= 10'
0.6
0.5
0.4 = 10
0.3
0.2
0.1
0 I I
0 5 10 15 20 25 30 35 40
Nu on side wall (RaE=104)
0.9 cold wall hot wall 0.8 RaE= 10 R 10
0.7
0.6
0.5
0.4
0.4 cold wall wall 0.3 RaE= 0 E=10
0.2
0.1
0
80 60 40 20 0 20 40 60 80
Nu on hot and cold walls (Ra 1=105)
Figure 5.8 Nusselt number distribution on the side walls
( Top: case 1; Bottom: case 2)
74
5.3 Natural Convection in LiquidGas Systems Without Internal Heating
In liquidgas systems, natural convection takes place in both phases. The natural convection in one phase may affect the convection in another phase. This causes the interaction between the convection cells in liquid phase and gas phase, and consequently affects the thermal performance of the system. The interaction between these convection cells also makes the computation more difficult to converge. In this section and the next section, a twophase system without phase change, with or without internal heating, is considered in order to study its convection and heat transfer characteristics. These investigations have not been found in the literature. New findings are obtained through numerical modeling experiments.
A rectangular cavity, with an aspect ratio of 2, filled with liquid and gas is shown in Figure 5.9 (a). The liquid and gas occupy the bottom and top halves of the cavity, respectively. Both the top wall and the bottom wall are adiabatic. The ratio of the liquid density to the gas density is assumed to be 10, the conductivity ratio of liquid to gas is assumed to be unity. The Prandtl numbers for both phases are assumed to be 1. The Grashof number for the gas phase, Gr, which is based on gas properties, is 10'.
The governing equations for different phases in nondimensional form are
, + = 0 (5.11) OX 6Y
c(p,UU) d(AUV) c]B 1
+  + V.(AVU) (5.12) f Ol Y cX Gr,
75
9(_,UV) 1w(AV) c I
OX + Y + VK v.( VV) + A.0o (5.13)
(U + (V) ,V2 (5.14) dX aY Prf, [The forms of Eqs (5.11)  (5.14) are different from those of (5.1) to (5.4), and for now the natural convection velocity scale in the gas phase is given by u= AV g, (5.15)
and is used instead of the heat diffusion velocity scale, where vg is viscousity for gas phase. This formulation has the advantage that all the variables can be maintained order 1, which may benefit the computation. The second order central difference scheme and a mesh with 81x41 grids are used.
5.3.1 Gas Convection Induced by Liquid Convection
The upper side walls, contacting the gas phase, are insulated. The lower side walls covering the liquid phase are maintained isothermal with T = 2 and T, = 1. Only the liquid phase is heated externally. Figures 5.9 (b) through (d) show the convection flow patterns in the liquid and gas phases, and the temperature distribution in the entire system, respectively.
In the liquid phase, a single clockwise rotating convection cell named as cell 1 is formed, moving upward on the left side and downward on the right side. Since no heat is provided from the side walls to the gas phase, no active convection which is directly
76
induced by buoyancy differences is expected to appear in the gas phase. Two counterrotating convection cells, however, are found in the gas phase, as shown in Figure 5.9(c). A counterclockwise convection cell, named cell 2, appearing on the right side of the gas phase is obviously induced by the convection cell 1 in the liquid phase, through the shear force interaction at the phase interface. Consequently, cell 2 drives another relatively weaker clockwise cell on the left side of the gas phase, which is cell 3. The natural convection cells 2 and 3 in the gas phase are induced, through the shear interaction with convection cell 1 in the liquid phase. Therefore convection in the gas phase is purely passive. The intensities of passive convection cells, cell 2 and cell 3, are much less than that of cell 1, therefore the flow pattern in the liquid phase is influenced slightly. Due to the weak convection effect, the heat transfer through the gas phase is very weak which can be seen from the isotherms. The heat transfer in the liquid phase is strongner.
5.3.2 Both Phases Heated
In this case, the entire side walls covering both liquid and gas phases are maintained isothermal, that is, T = 2 and T, = 1. Other conditions are the same. The convection flow patterns in the liquid and gas phases, and the temperature distribution in the entire system are shown in Figures 5.10(a), (b) and (c), respectively.
The temperature difference between the two side walls induces natural convection in the gas phase, as well as in the liquid phase. Three major convection cells appear. A
77
single clockwise cell in the liquid phase is also named cell 1, and the other two counterrotating cells in the gas phase are called cells 2 and 3, as treated in the previous case.
The counterclockwise rotating cell 2 in the gas phase is also induced by the convection cell 1 in the liquid phase. Cell 3 is caused by the temperature difference between side walls and the shear effect of cell 2. Thus cell 3 is generated both actively and passively. Due to the duel effect, the strength of this cell 3 is greater than that purely passive cell 3 of the previous case. On the other hand, cell 3 also affects cell 2, thus the convection of cell 2 is enhanced a little bit. The interactions of these convection cells in both phases are clear in this case. The intensities of passive convection cells 2 and 3, are much less than that of cell 1, therefore the flow pattern in the liquid phase is not changed significantly. Due to its weak convection strength, the heat transfer in the gas phase basicaly remains conduction characteristics. The temperature distribution in the liquid phase remains significant characteristics of natural convection because of its stronger convection strength.
Another similar case, where the Grashof number for the gas phase is increased to 105, is examined but not shown. The main features of the flow patterns are similar. Three convection cells are formed, two in gas phase and one in liquid phase. The temperature distribution is changed. Because of the higher Gry, the convection and heat transfer are enhanced. The interactions between these convection cells are also increased significantly. The enhanced interactions between these convection cells cause severe difficulty for the computational stability and convergence. An unsteady solution procedure, which is not required for the previous cases, has to be employed to obtain convergent results.
insulated
gas
cell 1
insulated
(a) Schematic of the cavity (b) Convection in liquid
cell 3
cell 2
c) Convection in gas (d) Isotherms
Figure 5.9 Natural convection of a twophase without internal heating
(case 1: upper side wall insulated, lower side wall isothermal, Gr = 103)
79
cell 3
cell I
(a) (b) (c) Convection in liquid: Convection in gas Isotherms
Cell 1: + Cell 2: ; cell 3: +
Figure 5.10 Natural convection of twophase without internal heating
( case 2: both side walls isothermal, Grg,=103 )
80
5.4 Natural Convection in LiquidGas Systems with Internal Heating
The flow and heat transfer characteristics of a twophase system without phase change but with internal heating is investigated in this section. Two cases, one rectangular cavity, one cylindrical container, are considered.
5.4.1 Case 1: Rectangular Cavity
The rectangular cavity with aspect ratio of 2 is filled with liquid and gas. The volume ratio of liquid to gas is 1. The energy equation, where internal heat generation is involved, is
1 Gr
V*(),VO) Pr r 2 (5.16) PrgGr g ' ' (GgE)12
Other equations retain the same forms as Eqs. (5.11)  (5.13). In Eq.(5.16) the constant GrgI is the internal Grashof number. The volumetric heat generation rate is proportional to its mass density and is uniformly distributed in each phase.
The top wall, bottom wall and the left side wall are insulated. The right side wall is kept isothermal with T,= 1. The density ratio of liquid to gas is 10, and the conductivity ratio of liquid to gas is unity. The external Grashof number based on the gas properties is 10'. The internal Grashof number is 316.
The results are shown in Figure 5.11. Three convection cells are formed, one in the liquid region and two in the gas region. Similar to the previous cases, one passive convection cell on the right side in the gas phase is induced by the convection cell in the
81
liquid phase. Another clockwise convection cell in the gas phase is induced by the internal heating and the shear of the counterclockwise cell. These convection cells interact and affect each other. The intensity of the convection in the liquid phase is much greater than that in the gas phase due to the higher mass density and higher heat generation rate in the liquid phase. The convection heat transfer in the liquid region is much stronger than that in the gas region where heat diffusion is dominant.
5.4.2 Case 2: Cylindrical Container
The results for the case of the cylindrical container are displayed in Figure 5.12. The container with aspect ratio of 2.5 is filled with half liquid and half gas. The forms of the governing equations are similar to those of Eqs. (5.11)  (5.13) and (5.16), except that the cylindrical coordinate system is used here. The top wall and bottom wall are insulated, and the side wall is isothermal with T, = 1. Only the right half of the container is considered owing to symmetric boundary conditions, and the axial centerline of the container is adiabatic. The external gas propertybased Grashof number is 10s, and the internal Grashof number which is also based on gas properties is 1.58x106.
The main flow patterns and heat transfer characteristics are similar to those of the previous case. Three convection cells are generated, one clockwise rotating cell in liquid and two counterrotating cells in gas. Since the centerline of the container is adiabatic, heat transfer and convection in the region close to the centerline are quite weak. The natural convection and heat transfer are relatively strong in the liquid region near the side
82
wall. In the liquid phase, the convective heat transfer is dominant. The heat diffusion is still quite significant in the gas phase due to the weak convection.
0007
0.
(a) (b) (c)
Convection in liquid: Convection in gas Isotherms
Cell 1: + Cell 2: ; cell 3: +
Figure 5.11 Natural convection of twophase with internal heating
(case 1: cavity, left wall adiabatic, right wall isothermal, GrgE =10' , Grg, =316)
83
1.56
(a) (b) (c)
Convection in liquid: Convection in gas: Isotherms
Cell 1: + /,,=2.0x10O4 V,,,=8.4x0l4 Ay=5x10 Cell 2: ; cell 3: +
Figure 5.12 Natural convection of twophase with internal heating
(case 2: cylindrical container, side wall isothermal, GrgE =105 , Grg, =1.58x106)
CHAPTER 6
BULK EVAPORATION/CONDENSATION WITH INTERNAL HEATING
AT MICROGRAVITY AND NORMAL GRAVITY LEVELS
6.1 Brief Preview
So far, only conduction dominated effects under limited conditions have been taken into account for the bulk evaporation and condensation processes. Like all phase change processes, evaporation and condensation are also inevitably subject to the effect of convection induced by either buoyancy or surface tension variations. The existence of the convection can substantially affect the flow field, the temperature field, and subsequently the performance of the phase change process. Numerous investigations on phase change with convection have been conducted, which can be found in Ostrach (1983), Crank (1984), Brice (1986), Brown (1988), Gau and Viskanta (1986), Nadarajah and Narayanan (1990), Shyy et al. (1991, 1992a, b, 1994a, 1995, 1996) and Shyy (1994a, b, 1995), etc.
In this Chapter, the modeling of the bulk evaporation and condensation processes, associated with internal heat generation and under various gravity levels, is presented. The system to be modeled is shown in Figure 1.2 similar to that of Figure 4.1. Initially the liquid and vapor are assumed to be saturated at the liquidvapor interface, while it is also assumed that the liquid is subcooled and the vapor is superheated in their own regions, respectively. The top wall and bottom wall are insulated. The side wall is cooled externally by a constant outward heat flux which removes heat from the system and tends to 84
85
condense the vapor phase. The internal heat generation provides heat to evaporate the liquid phase. Since the system is to be operated under different gravity conditions, there is convection, both in the liquid and vapor phases, induced by the buoyancy force due to the density variations. The temperature along the liquidvapor interface does not change significantly, because the curvature of the interface is moderate. The liquidvapor interface can be reasonably assumed as isothermal, thus the thermal capillary effect can be neglected. The phase change processes are controlled by conduction and convection. The performances of the phase change can be quite different at various gravity condition, as can be expected. The flow field, the temperature distribution and the motion of the interface are strongly dependent on the glevels, as will be demonstrated later.
6.2 The Governing Equations
The equations governing the thermal energy transport and phase change processes, with convection and heat generation, are similar to Eq. (4.4), except for the additional convection terms. The equations describing natural convection in both phases are similar to those of Eqs (5.1) (5.3), while the vapor density in a certain initial state is chosen as a reference basis for both phases.
The governing equations in cylindrical coordinate system for different phases in dimensional form are
+ V . p ) = 0 (6.1)

Full Text 
xml version 1.0 encoding UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID E9296YWM7_FHLQB6 INGEST_TIME 20150325T20:34:02Z PACKAGE AA00029735_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES
PAGE 1
MODELING OF BULK EVAPORATION AND CONDENSATION WITH INTERNAL HEATING AT VARIOUS GRAVITY LEVELS By ZHONGTAO DING 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 1995 UNIVERSITY OF FLORIDA LIBRARIES
PAGE 2
Copyright 1995 by Zhongtao Ding
PAGE 3
To my wife Zhaohui Cheng
PAGE 4
ACKNOWLEDGMENTS It has been a wonderful experience to pursue further education and academic development at the University of Florida, where I have known so many nice professors and friends who have helped me help in many ways. From them I learned that the way to achieve goals through hard work, determination, and effort. First of all, I appreciate the insightful guidance, constant support and encouragement provided by Professor Samim Anghaie. Without his enthusiasm and commitment to this challenging project, my exploration in this work could not have been accomplished. His help in finding solutions to my problems both personal and professional will never be forgotten. I sincerely thank Professor James F. Klausner, Chairman of the supervisory committee, for all his spiritual encouragement and technical assistance during the course of this research. His helpful suggestions and comments are greatly appreciated. I would like to express my special thanks to Professor Wei Shyy with greatest sincerity. I owe much of my knowledge on computational modeling of phase change problems to his direct and sincere help. I have learned a lot from the discussions with him and his students, from his class, his books and his many papers. I also want to extend my appreciation to Professors C. K. Hsieh and R. Mei for their serving on the supervisory committee. Their suggestions have always been constructive and helpful. I would like to express appreciation for all spiritual encouragement and financial support from Professor Nils J. Diaz. iv
PAGE 5
Finally, I want to thank my parents and sister who have supported me unconditionally and expect nothing from me in return. I wish to thank my wife, Zhaohui Cheng. I can not imagine that I could have finished my research and dissertation without her understanding, support and patience throughout the program.
PAGE 6
TABLE OF CONTENTS page ACKNOWLEDGMENTS j v LIST OF TABLES viii LIST OF FIGURES ix NOMENCLATURE xj ABSTRACT xiv CHAPTERS 1 INTRODUCTION j 2 LITERATURE SURVEY 9 2.1 SolidLiquid Phase Change Problems 9 2.2 Numerical Modeling of SolidLiquid Phase Change Problems 11 2.3 LiquidVapor Phase Change Problems 12 2.4 Natural Convection in the Single Phase 13 2.5 Phase Change with Convection 14 3 1 D CONDUCTION BULK EVAPORATION AND CONDENSATION 15 3.1 The Internal Energy Formulation 19 3.2 Analysis and Assumptions 22 3.3 The Vapor Phase Fraction Update 28 3.4 Properties of the Control Volume Containing Two Phases 32 3.5. Results and Discussion 33 3.5.1 Accuracy Aspect 33 3.5.2 Comparison of H Formulation and E Formulation 38 3.5.3 Comparison of EBased and TBased Update Methods 38 3.5.4 Bulk Evaporation and Condensation 40 3.6 Concluding Remarks 41
PAGE 7
4 BULK EVAPORATION/CONDENSATION WITH INTERNAL HEATING AT ZERO GRAVITY 45 4.1 The Governing Equations 47 4.2 Conductivity at Control Volume Surfaces 49 4.3 Results and Discussion 51 5 NATURAL CONVECTION WITH INTERNAL HEATING IN SINGLEPHASE AND TWOPHASE SYSTEMS 59 5.1 External Natural Convection in a SquareCavity Benchmark Problem 60 5.2 Natural Convection in a SquareCavity with Internal Heating 67 5.3 Natural Convection in LiquidGas Systems without Internal Heating 74 5.3.1 Gas Convection Induced by Liquid Convection 75 5.3.2 Both Phases Heated 76 5.4 Natural Convection in LiquidGas Systems with Internal Heating 80 5.4.1 Case 1: Rectangular Cavity 80 5.4.2 Case 2: Cylindrical Container 81 6 BULK EVAPORATION/CONDENSATION WITH INTERNAL HEATING AT MICROGRAVITY AND NORMAL GRAVITY LEVELS 84 6.1 Brief Preview 84 6.2 The Governing Equations 85 6.3 Numerical Treatment 88 6.4 Results of MicroGravity (10" 3 g ) 90 6.5 Results of Normal Gravity (1g ) 106 7 CONCLUSIONS AND RECOMMENDATIONS 115 REFERENCES 119 BIOGRAPHICAL SKETCH 126 vii
PAGE 8
LIST OF TABLES Table page 2.1 Brief summary of major characteristics of phase change problems 10 5 . 1 Comparison of present solutions with benchmark solutions 63 viii
PAGE 9
LIST OF FIGURES Figure page 1 . 1 Schematic diagram of a typical nuclear thermionic fuel element 2 12 Schematic diagram of the multiphase nuclear fuel system 4 3 . 1 Schematic diagram of the bulk evaporation and condensation 17 3.2 The control volume containing two phases 33 3.3 Solution features for bulk evaporation 35 3.4 Comparison of two formulations for evaporation 36 3.5 Comparison of two update methods for evaporation 37 3.6 Bulk evaporation and condensation 42 4.1 Schematic diagram of the twophase cell operated at 0g 46 4.2 The control volume contacting other material 50 4.3 Parameter for heat generation rate 52 4.4 The evolution of liquidvapor interface (0g) 56 4.5 The evolution of temperature field (0g) 57 4.6 Temperature surf of the system at steady state (0g) 58 4.7 Variation of saturation temperature (0g) 58 5 1 Domain and boundary conditions for the benchmark problem 61 5.2 Flow patterns and isotherms of buoyancy driven natural convection 64 5.3 Nusselt number distribution on side walls 65 5.4 Temperature profile at horizontal centerline 66
PAGE 10
5 . 5 Domain and boundary conditions for internal heating induced convection 70 5 6 Flow pattern and isotherms of internal heating induced convection ( case 1 ) ...71 5.7 Flow pattern and isotherms of internal heating induced convection ( case 2) ... 72 5.8 Nusselt number distribution on the side walls 73 5.9 Natural convection of twophase without internal heating ( case 1) 78 5.10 Natural convection of twophase without internal heating ( case 2) 79 5.11 Natural convection of twophase with internal heating ( case 1 ) 82 5.12 Natural convection of twophase with internal heating ( case 2) 83 6.1 Schematic diagram of mass reflux 96 6.2 The evolution of liquidvapor interface ( 10" 3 g ) 98 6.3 The evolution of isotherm ( 10" 3 g ) 99 6.4 Temperature surf at x = 2.0 ( 10" 3 g) 100 6.5 Variation of saturation temperature ( 10" 3 g) 100 6.6 Solutions with different grid sizes and schemes 101 6.7 The evolution of liquidvapor interface with coarse grid 101 6.8 Flow patterns in vapor and liquid phases x = 2.0 ( 10" 3 g ) 102 6.9 Velocity U profile at different axial positions (10' 3 g) 103 6.10 The evolution of liquidvapor interface ( 1g) 109 6.11 The evolution of isotherms ( 1g) HO 6.12 Temperature surf at x = 2.0 ( 1g) HI 6.13 Variation of saturation temperature ( 1g ) Ill 6. 14 Flow patterns in vapor and liquid phases x = 2.0 ( 1g ) 112 6.15 Velocity U profile at different axial positions (1g) 113 6.16 Comparison of interface and temperature distribution at different g levels 1 14
PAGE 11
NOMENCLATURE A.C.D constants E internal energy [J] G internal heat generation rate [J/kg] H enthalpy [J] Q heat quantity [J] T temperature [K] U dimensionless velocity V volume [m ] or dimensionless velocity W work [J] X,R dimensionless coordinates c p specific heat at constant pressure [J/kg/K] c v specific heat at constant volume [J/kg/K] e specific internal energy [J/kg] h specific enthalpy [J/kg] f vapor phase fraction 8 gravity acceleration [m/s 2 ] 1 length scale [m] P pressure [bar] t time [s] XI
PAGE 12
u,v velocity [m/s] x, r cylindrical coordinates [m] Nondimensional Characteristic Numbers Gr Grashof number Pe Peclet number Pr Prandtl number Re Reynolds number Ra Rayleigh number Si Stefan number Greek symbols a thermal diffusivity [m 2 /s], or thermal diffusivity ratio Â£ phase change temperature interval k thermal conductivity [W/K/m], or its ratio O order of magnitude 9 dimensionless temperature p density [kg/m 3 ], or its ratio r dimensionless time w stream function xii
PAGE 13
Subscripts / interface position, mixed phase i interface / liquid phase r characteristic scale s, sat saturation point v, vap vapor phase Superscripts k iteration level n time level arrow bar vector bar nondimensional quantity xiii
PAGE 14
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 OF BULK EVAPORATION AND CONDENSATION WITH INTERNAL HEATING AT VARIOUS GRAVITY LEVELS By Zhongtao Ding December, 1995 Chairman: James F. Klausner CoChairman: Samim Anghaie Major Department: Mechanical Engineering A comprehensive numerical model is developed to simulate the dynamic processes of bulk evaporation and condensation in an encapsulated container, associated with internal heat generation under different gravitational strengths. A complete set of equations governing the conservation of mass, momentum and energy is solved. The thermal performance of a multiphase nuclear fuel element at zero gravity, microgravity, and normal gravity is investigated. In the simulation, the numerical solutions have revealed much useful information, including the evolution of the bulk liquidvapor phase change process, the evolution of the liquidvapor interface, the formation and development of the liquid film covering the side wall surface, the temperature distribution, the convection flow field, and the strong dependence of the thermal performance of such multiphase fuel cells on the gravity conditions.
PAGE 15
In regard to the development of computational techniques, the internal energy formulation for the constant volume phase change process is proposed. The heat transfer mechanism for this type of phase change problem, including the convection effect induced by the density change, is analyzed and discussed. The computational performance of the internal energy formulation is evaluated. The performance of two update methods for the vapor phase fraction, the Ebased and Tbased methods, is investigated. Simulation of a onedimensional conductiondriven bulk evaporation and condensation problem is performed and the evolution of the phase change process is obtained. Some natural convection problems with or without internal heating in the singlephase and twophase systems are studied. The solutions agree quite well with benchmark solutions and experimental results. It is revealed that interactions between the convection cells in liquid and vapor phases can substantially complicate the flow fields and the computational modeling.
PAGE 16
CHAPTER 1 INTRODUCTION The thermionic energy converter is a nonmechanical gaseous electronic device for converting heat energy directly into electric potential through thermionic electron emission. In a thermionic diode, electrons are emitted from a hot electrode, which is the emitter, and collected by a colder electrode, which is the collector, at a higher potential energy (lower electrical potential). Part of the heat removed from the emitter due to evaporating electrons is rejected to the collector by condensing electrons. The remaining heat is converted into electric power in the load as electrons return to the emitter potential. The physics governing thermionic energy conversion is well understood (Hatsopoulos and Gyftopoulos, 1973), and this energy conversion technique is attracting wide interest in its application toward a variety of power conversion systems (Lee et al, 1993, Young et al., 1993). In a review paper on thermionics, Rasor (1991) outlined the history, application options, and ideal performance of thermionic energy converters and described the basic plasma types associated with various modes of converter operation. Various types of thermionicconverter systems, such as solar thermionic generators, radioisotope thermionic generators, chemical thermionic generators and nuclear reactor thermionic generators, have been or are currently being developed. The nuclear reactor thermionic generators are simply thermionic converters supplied with heat produced by nuclear fission. The fission heat is conducted from the reactor core to the emitter, and the rejected heat is conducted from the collector and carried away via a coolant. A typical schematic configuration of the thermionic fuel element (TFE) fueled with solid nuclear fuel U0 2 is shown in Figure 1.1. The nuclear fuel resides in the center
PAGE 17
of the TFE. The space between the emitter and collector is filled with an easily ionizable rarefied medium, typically cesium vapor, in order to enhance the electrical performance. A liquid metal coolant flows in the annulus channel surrounding the TFE. Recently a variety of investigations involving such a TFE have been carried out that consider the thermal and electrical performance, system simulation and design, and application as a space power system, Lewis el al. (1991), McVey and Rasor (1992), Young et al. (1993). insulator gap coolant channel cladding collector Figure 1 1 Schematic diagram of a typical nuclear thermionic fuel element
PAGE 18
The thermal performance of a thermionic fuel element is tightly coupled with its electrical performance. The electrical output current density is directly related to the emitter temperature, that is, the current density increases significantly with emitter temperature, which can be illustrated by a simple example of an ideal diode thermionic converter as follows J = AT E 2 exp[(V + c )/kT E ] (1.1) where J is the current density, A and k are constants, V is the output, c is the collector work function, and T E is the emitter temperature. Consequently a nuclear fueled TFE, which can provide high temperature, is well suited for this application. The maximum operating temperature of solid nuclear fuel U0 2 could be as high as 1800 K In order to further enhance the TFE output, a higher operating temperature is required. Such higher temperatures could be realized if liquid nuclear fuel is employed. If nuclear fuels such as UF 4 or UF4/UO2 in a liquid or vapor state were used, they could sustain a temperature of 2600 K or higher. The large increases in operating temperature will raise the power output significantly. The use of a multiphase nuclear fuel reactor in conjunction with a thermionic converter provides the potential for a compact, higher power output electrical generator that is well suited for space applications. Certain problems related to the thermal performance of TFE arise with regard to the application of multiphase nuclear fuel reactors These include the phase change (PC) process between the liquid and vapor, the behavior of evaporation and condensation associated with internal heating and external cooling, the location and motion of the liquidvapor interface, the temperature distribution, the influences of gravitational fields, and so on. With such a system, particular attention must be paid to limiting the vapor pressure that is safe for the containment vessel because evaporation in a closed container
PAGE 19
will raise the vapor pressure very rapidly. A thorough investigation of these issues will be required to fully apply this concept. insulated liquid film 'vapor condensation evaporation liquid pool Zgen.l insulated cooled Figure 1 ,2 Schematic diagram of the multiphase nuclear fuel system A cylindrical container filled with multiphase nuclear fuel, such as UF 4 , is shown in Figure 1.2. The fuel experiences liquidvapor phase change processes: evaporation and
PAGE 20
condensation. The internal heat provided by the nuclear fission reaction evaporates the liquid, while the side wall cooling tends to condense the vapor. Globally, the evaporation of liquid and condensation of vapor take place simultaneously. The evaporation and condensation processes could be driven by either pure heat conduction or by both heat conduction and convection, depending mainly on the gravity field. Under a zerogravity condition, buoyancydriven convection is nonexistent. However, under microgravity or onegravity fields, buoyancyinduced convection is present. The location of the liquidvapor interface may vary depending on the mode of heat transfer, the internal heat generation intensity and spatial distribution, and the gravity field. Such an interface could be defined as a free and moving interface. To solve the phase change problem analytically or numerically, the free and moving interface is not known in advance, so that its position has to be determined as a part of the solution. Thus the solution to this heat transfer problem is quite complex. With such background, specific issues are to be investigated. For the phase change processes between liquid and vapor, such as bulk evaporation and condensation involved in the multiphase system, no similar investigation has been found. Limited investigations have been conducted. Investigations of this problem in different parametric regimes have been performed by Shyy (1994b). To the best of the author's knowledge, no prior researcher has considered the general bulk condensation and evaporation processes in an enclosed container, under various gravity fields, with internal heat generation. Liquidvapor phase change problems are also characterized by a large step change in the density across the phasic boundary. The singular jump of the density at the interface, other than the jump of temperature gradient which exists unavoidably for all phase change problems, brings substantial challenges to computational fluid dynamics and heat transfer. Since the investigations for the solidliquid phase change problems, such as melting and solidification, are numerous, some successful experiences and lessons can be drawn from those efforts and applied to the study of the bulk evaporation and condensation
PAGE 21
problems. Appropriate modeling techniques, such as formulation and numerical treatment, must be developed. For the solidliquid phase change problems, the thermodynamic and transport properties between two phases vary moderately during phase change processes. These properties, however, vary significantly or dramatically during the phase change processes between liquid and vapor. This causes some difficulties. Scaling and order of magnitude analysis need to be conducted in order to make simplifications that can lead to a concise and efficient solution procedure for the governing equations. Some special treatment and closure relationships, such as those for vapor phase fraction, the vapor density, the saturation pressure and temperature, have to be applied and implemented into the computation. Like all phase change processes, evaporation and condensation are also inevitably subject to the effect of convection induced by either buoyancy or surface tension variations. The existence of the convection can substantially affect the flow field, the temperature field, and subsequently the performance of the phase change process. The effect of the convection on the evaporation and condensation processes is to be investigated. For the conventional solidliquid phase change problems, heat is usually provided externally from the boundary to the system. For the nuclear system, heat is generated internally via a nuclear fission reaction. The problems of phase change induced by internal heating are different from those induced by external heating. The internal heating, like external heating, results in temperature and density gradients, and induces the buoyancydriven flow. Such natural convection induced by internal heating is defined as internal natural convection, to be distinguished with the classical natural convection caused by external heating/cooling. This internal natural convection has drawn certain attention for some investigations (Acharya and Goldstein, 1985; Lee and Goldstein, 1988) The behavior of the internal natural convection and its influence on the phase change processes, in particular on the bulk evaporation and condensation processes, are to be investigated.
PAGE 22
The objective of this dissertation is to develop a numerical method to model the dynamic processes of bulk evaporation and condensation, with internal heat generation and natural convection, in an encapsulated container under different gravitational levels. Based on this model, the following topics will be studied in depth: 1) the bulk evaporation and condensation processes with different modes of heat transfer, that is, conduction and convection, 2) the movement of the liquidvapor interface, 3) the evolution of the temperature field, and 4) the influence of various gravitation fields on the performance of the phase change process. The modeling and results can be applied to validate the concept of a multiphase nuclear system, including the multiphase TFE, and to predict the performance and assist in devising strategies for its thermal management and control for space applications. The current research status concerning these problems mentioned above is reviewed in Chapter 2. The modeling of a onedimensional, externally heated or cooled and conductiondriven bulk evaporation/condensation problem is described in Chapter 3. This phase change problem is well suited to develop the mathematical formulation, the closure relationships, and the updating of the vapor phase fraction. The results of the model are used to investigate the basic behavior of the evaporation/condensation processes, and to examine the computational efficiency and accuracy. The numerical approach developed in Chapter 3 is then expanded to investigate the bulk evaporation/condensation processes with internal heat generation at zero gravity level, which is presented in Chapter 4. The evolution of the phase change processes, the temperature field, and the motion of the liquidvapor interface topology are demonstrated. Chapter 5 describes the modeling of the combined internal and external natural convection in the singlephase system and twophase system without phase change. To fully understand the mechanism and behavior of the internal natural convection and to establish the confidence of the computational accuracy, the external natural convection is also computed and comparisons with the benchmark solutions are made. In Chapter 6, the
PAGE 23
bulk evaporation/condensation processes with internal heat generation and natural convection under different gravitational fields are modeled. A complete set of equations governing the conservation of mass, momentum and energy is solved. The effects of the buoyancyinduced convection via internal heating and external cooling, and different gravitational fields, on the phase change processes are studied. Finally, the conclusions summarizing the achievements of this work and recommendations for further research are provided.
PAGE 24
CHAPTER 2 LITERATURE SURVEY Previous investigations of phase change problems are numerous. Those studies, which can be traced back as early as Clausius, Clapeyron, Stefan, Nusselt, and others, have been conducted both experimentally, analytically and are currently being attempted using numerical simulations and approximation methods. In their annual review of heat transfer papers, Eckert et al. (1990, 1991 and 1992), suggested that most of these current works have focused on melting and solidification, pool and flow boiling, film and droplet evaporation, filmwise and dropwise condensation, with applications in thermal energy storage, solidification of alloys and metals, casting processes, crystal growth, heat pipes, heat pumps, boilers, heat exchangers, and nuclear or chemical reactors. The major characteristics of the phase change problems and the current research status are briefly summarized in Table 2.1. 2. 1 SolidLiquid Phase Change Problems Early analytic work on phase change between solid and liquid include studies by Lame and Clapeyron (1831) and Stefan (1891). The fundamental features of this type of problem and some basic analytical solution methods were discussed by Ozisik (1980), Crank (1984) and Yao et al. (1989). Crank (1984) presented a broad but reasonably detailed account of the mathematical solution, both analytical and numerical, of these problems. Due to the existence of the interface, difficulties were encountered in solving
PAGE 25
these problems to the extent that only a handful of very simple cases can be theoretically analyzed. Other efforts on the analytical and numerical treatment for the interface include the source and sink method that has been used and developed by Hsieh et al. (1993) and Li (1993). In this method, the phase change conditions at the interface are introduced into the heat conduction equation as an impulse internal heat source or sink. Then the equations for both regions can be solved simultaneously. Table 2.1 Brief summary of major characteristics of phase change problems Moving Boundary Phase Change Problems PC Between Solid Liquid Liquid Vapor flow or film or bulk evaporation Process Melting and pool drop and solidification boiling condensation condensation Heating externally from externally externally or /Cooling the system boundary from boundary internally heating Thermoconstant or physical change moderately change significantly properties Conduction conduction conduction conduction and conduction conduction and/or only and convection only and convection convection convection Solution analytical numerical experimental Approach numerical empirical empirical numerical analytical numerical Research numerous investigations numerous very few found, Status investigations being conducted
PAGE 26
At the University of Florida, fruitful investigations on various phase change problems have been carried out by Shyy and Chen (1991), Shyy and Rao(1992), Shyy et al. (1992, 1993a, b, 1996), Shyy (1994a, b), Hsieh et al. (1992, 1993), Li (1993), and Hunger al. (1995). 2.2 Numerical Modeling of SolidLiquid Phase Change Problems Since the 1970s, greater attention and interest have been paid to purely numerical analysis and simulations, which is reflected in the large number of publications covering these topics. Shyy (1994a) and Shyy et al. (1996) summarized a substantial portion of their research activities, and presented the various computational elements important for the prediction of complex fluid flows and interfacial transport. The numerical attempts can be divided into two groups, based on the choice of dependent variables to be used In the first group, temperature is the sole dependent variable, and energy equations are written separately in the solid and liquid regions. At the interface, the temperature is continuous, while its gradient is not continuous due to the release or absorption of latent heat The boundary conditions at the interface are used to couple and connect the temperature in the two regions. The temperature gradient jump is difficult to capture and requires extraordinary effort. In the second group, enthalpy is used as a dependent variable along with temperature. This formulation is termed the enthalpy method, and has gained popularity. In the enthalpy model, the interface is eliminated from consideration in the calculations, and the problem is made equivalent to one of nonlinear heat conduction without phase change. The shortcoming of this method is that it is difficult to determine the interface position due to its oscillations during the numerical solution process. Shamsundar and Sparrow (1975), Shyy (1994a) and Shyy et al. (1996) summarized these features and demonstrated the usefulness of the enthalpy formulation.
PAGE 27
Many investigators have employed the enthalpy formulation, including Shyy et al. ( 1 992a, 1993b, 1996), Shyy (1994a, b), Voller and Prakash (1987), Voller and Swaminathan (1991), Yao and Prusa (1989). 2.3 LiquidVapor Phase Change Problems In recent years, interest in the other forms of phase change processes, such as evaporation and condensation, has grown. Carey (1992) systematically summarized the state of the art in boiling and/or condensation phenomena. His book focuses on some fundamental elements and the basic physical mechanism of condensation and vaporization processes and provides a broad review of the research activities in these areas. Among the investigations on condensation and evaporation, those on film condensation are well developed, particularly on laminar film condensation, due to its simplicity. The development of the current theory and the prediction of heat transfer rates associated with film condensation have its origin in the analyses presented by Nusselt (1916). Subsequently, many modifications of this analysis have appeared, each one further relaxing the restrictive assumptions of the Nusselt analysis. Considerations have been given to gravityflow condensation, pure forced convection condensation in the presence of a body force, such as in Sparrow et al. (1967). The influence of entrained noncondensable gas in the vapor and the associated interfacial resistance, that is, a temperature jump at the interface between the condensate and the vapor, on these phase change processes have also been investigated by Sparrow and Eckert (1961), Sparrow and Lin (1964), Sparrow et al. (1967), and Rohsenow (1973). Recently some studies on general film condensation transients and the interfacial condensation model have been conducted by Flik and Tien (1989) and Gerner and Tien (1989). With the development of computational fluid dynamics and numerical heat transfer, more and more attempts have been made to model and simulate condensation phenomena, for example, Jones and Renz
PAGE 28
13 (1974), Cossmann et al (1982), and Gaultier et al. (1993). However, most of these activities, either analytical or numerical, have focused on flow condensation problems with open plates, channels or open containers where the vapor pressure, saturation temperature and latent heat are maintained constant. These problems, where forced convection plays a dominant role, are usually parabolictype problems, and consequently are much simpler than those in the present study in which the saturation temperature, pressure and latent heat continuously change. A few computational simulations of evaporation in containers have been carried out recently. Jang el al. (1990) modeled a heat pipe startup from the frozen state by using onedimensional equations for different phases Harley and Faghri (1994) performed numerical analysis of thermosyphons including the falling condensate film. Shyy (1994b) developed a computational method for predicting the twophase transient flow and heat transfer characteristics within a constant pressure reservoir of a capillarypumped loop and investigated the phase change of liquid and vapor, under various conditions. 2.4 Natural Convection in the Single Phase Considerable attention has been given to the study of natural convection of a single phase fluid in enclosures. The temperature gradients across the system that are supplied externally induce density gradients and consequently the buoyancy driven flows. A great deal of research, both numerical and experimental, on various types of convective heat transfer problems has been conducted, including that of Ostrach (1982, 1983, 1988), Gebhart el al. (1988), de Vahl Davis (1983), Hortmann et al. (1990), Shyy et al. (1992), Shyy (1994a), Bejan (1984), Carpenter et al. (1989), Ramaswamy and Jue (1992), Bergman and Ramadhyani (1981). Ostrach (1982, 1983, 1988), Gebhart et al. (1988) and Bejan (1984) reviewed and summarized those activities and achievements. Among the large number of the natural convection problems, the buoyancydriven flow in a square
PAGE 29
cavity with vertical side walls heated differentially, so called a benchmark problem for natural convection, is a suitable vehicle for testing and validating computer codes. Accurate benchmark solutions have been published by de Vahl Davis (1983) and by Hortmann et al. (1990) and are available for comparison with the present computations. Natural convection induced by external heating is of importance in many applications. Natural convection in an enclosure in which the internal heat generation is present is also of prime importance in certain technological applications. Examples are postaccident heat removal in nuclear reactors and some other geophysical or astrophysical problems. In contrast to the extensive investigations on the conventional external heating convection problems, natural convection induced simultaneously by external heating/cooling and internal energy sources has received limited attention, such as Archarya and Goldstein (1985), Lee and Goldstein (1988), and Fusegi et al. (1992). The presence of internal heat generation provides an additional dynamic in overall convective flow systems and offers a complicated yet challenging aspect of the convection under consideration, 2.5 Phase Change with Convection With no exception, convection induced by either buoyancy or surface tension variations will also be involved in evaporation and condensation like all phase change problems. The involvement of convection can substantially affect the phase change process. Numerous investigations on the solidliquid phase change problems with convection can be found in Ostrach (1983), Crank (1984), Brice (1986), Brown (1988), Gau and Viskanta (1986), Nadarajah and Narayanan (1990), Shyy and Chen (1990,1991), Shyy and Rao (1992, 1994), Shyy et al. (1992, 1996) and Shyy (1994a, b)
PAGE 30
CHAPTER 3 1D CONDUCTION BULK EVAPORATION AND CONDENSATION In recent years, phase change processes such as melting and solidification, boiling and condensation have attracted increasing interest in their development and applications. Among these activities, numerical modeling and simulation of the phase change processes have been the focus of many investigators. Numerical modeling of heat transfer problems involving phase change has been discussed by a large number of investigators, such as Crank (1984), Yao and Prusa (1989), Voller and Prakash (1987), Voller and Swaminathan (1991), Samarskii et al. (1993), Prakash et al. (1987), Shamsunder and Sparrow (1975), Hsieh and Choi (1992), Hsieh and Akbari (1993), Li (1993), Hung et al. (1995), Shyy and Chen (1990, 1991), Shyy and Rao (1992, 1994), Shyy et al. (1993a, b, 1994, 1996), Shyy (1994a, b). Shyy (1994a) and Shyy et al. (1996) presented a comprehensive discussion of the various computational elements important for the prediction of complex fluid flows and interfacial transport. Most of the noted investigators have approached the solidliquid phase change problems, such as those involved in melting and solidification processes. There has also been an ongoing interest in the investigation of liquidvapor phase change problems, such as evaporation and condensation. More attempts have been made to numerically model and simulate condensation phenomena. In particular, a number of investigators have developed numerical models for flow condensation problems, for example, Jones and
PAGE 31
Renz (1974), Cossman et al. (1982), Gautier et al. (1993), Tournier and ElGenk (1992). Most of these numerical modeling efforts have focused on parabolic flow condensation problems where forced convection plays a dominant role. In this type of condensation problem, the vapor pressure, saturation temperature and latent heat have been treated as constant properties. In the present study, the saturation temperature, and pressure are modeled as dependent variables In a more recent work, Jang et al. (1990) have modeled a heat pipe startup from the frozen state by using a system of simplified one dimensional equations. Shyy (1994b) has successfully developed a computational method for predicting the twophase transient flow and heat transfer characteristics within a constant pressure reservoir of a capillarypumpedloop and investigated the phase change of liquid and vapor, under various conditions. The present approach, to be discussed later in this chapter, has expanded the mathematical modeling of the phase change problems by including the pressure as a dependent variable and using the internal energy formulation. The numerical models for the phase change problem can be divided into two groups, based on the choice of dependent variables. In the first group, temperature is the sole dependent variable, and energy equations are written separately for the solid and liquid phases. In the second group, enthalpy is used as the main dependent variable along with the temperature. This formulation which is widely used by many investigators is known as the enthalpy method. In the enthalpy model, the phasic interface is eliminated from the formulation and the problem is made equivalent to a nonlinear heat conduction problem without phase change. Shamsunder and Sparrow (1975) and Shyy (1994a) summarized these features and demonstrated the usefulness of the enthalpy formulation
PAGE 32
The enthalpy formulation has been employed in many investigations, for example, Shamsunder and Sparrow (1975), Voller and Prakash (1987), Voller and Swaminathan (1991), Shyy (1994a, b), Shyy and Rao (1994), Shyy et al. (1996), Yao and Prusa (1989). heated heat input liquidinsulated vapor W444 insulated insulated (a) Evaporation I 1 vr g liquid vapor 4^+4 insulated cooled heat output (b) Condensation i Vp Figure 3.1 Schematic diagram of the bulk evaporation and condensation The purpose of the work in this chapter is to develop a consistent numerical method for the computation of the bulk evaporation and condensation processes involved in a specific class of liquidvapor phase change problems under constant volume and no
PAGE 33
buoyancy conditions. To develop a generic numerical model and obtain the basic features of this phase change process, a simplified, one dimensional, conduction driven problem is particularly designed. This problem can be depicted with the case presented in Fig.3.1. A rigid, cylindrical container is partially filled with vapor and liquid. Initially, the vapor and liquid are kept saturated at a certain pressure and uniform temperature. Well selected boundary conditions on the walls are imposed to activate the pure bulk evaporation or condensation, and also to avoid the buoyancy effect and bubble generation. While keeping the side wall insulated, only heating the top wall or cooling the bottom wall induces pure evaporation or condensation, respectively. The volumetric expansion or contraction due to the density change between two phases induces convection in the vapor phase. The scaling analysis conducted in this chapter indicates that this convection effect is quite week and negligible. It can be reasonably assumed that the phase change processes are driven entirely by the conduction and the release or absorption of latent heat at the phase change interface. Similar situations are also possible for the multidimensional phasechange problems under zero gravity condition. No expansion work is involved in this constant volume phase change process. In the absence of mechanical work, the first law of thermodynamics indicates that the heat exchange between the system and the ambient equals the change of the internal energy of the system. The released or absorbed latent heat associated with the phasechange is equal to the change of internal energy. Based on the physical conditions of this system, a compact internal energy formulation of the energy equation is proposed. This formulation,
PAGE 34
similar to the enthalpy formulation, is a single region formulation, where one set of governing equations can describe both phases. 3 . 1 The Internal Energy Formulation The first law of thermodynamics, for any thermodynamic system, is 8Q = dE + 8W=dH(pdV + Vdp) + 8W (3.1) where E is the internal energy, H is enthalpy, Q is the heat input, Wis the work output, p and V are the pressure and volume of the system, respectively. For the phase change process involved in the rigid and enclosed container with constant volume, there is no mechanical work. Therefore the first law for this closed and rigid system reduces to SQ = dE = dHVdp (3.2) which indicates that the heat exchange between the system and the ambient causes the internal energy, the enthalpy and pressure of the system to change. Therefore the released or absorbed latent heat associated with the phase change induces the changes of internal energy, enthalpy and pressure. The quantity of the heat exchange is equal to the internal energy change, but not equal to the enthalpy change. The conventional enthalpy formulation of the energy equation, commonly used for the phase change between solid and liquid, should include a transient pressure term for this phase change between liquid and vapor. The weak form of the energy equation based on enthalpy formulation is P Â— = Â— + V(kV7") (3 3) H Dt Dt K ' y '
PAGE 35
20 This formulation including the transient pressure term introduces additional complexity and does not seem to be the best choice for this problem. A more attractive alternative, the internal energy based formulation for the energy equation which does not include the pressure term, is proposed and developed in this work. Similar to the enthalpy formulation used for melting and solidification problems, the internal energy formulation does not require to explicit tracking of the interface position. The weak form of the energy equation with internal energy as the primary variable can be expressed as follows: De PÂ— = V(kVT) (3.4) where the internal energy e may be expressed, in the form of a summation of sensible heat and latent heat, for the liquid, vapor and mixed phases as: (3.5) where T sat is the saturation temperature, Ae is the internal energy difference between saturated vapor and saturated liquid (latent heat), and the constant volume specific heat c,'s are functions of temperature provided by Reynolds (1979). The vapor phase fraction/ is zero in the region filled with the liquid phase and unity in the region occupied by the vapor phase. The vapor phase fraction lies between zero and unity when the control volume is undergoing phase change. Substituting Eq.(3.5) into the governing Eq.(3.4) yields ec vl T /=o e = c v.,?;Â„+/Ae 0
PAGE 36
P,^.,T) = y(K,VT) / = o Pfjfa&mV.WTiptjJjL 0
PAGE 37
used, the liquid or the vapor phase is identified by the local internal energy. Either case allows the vapor phase fraction in a computational cell to be calculated. Using the ClausiusClapeyron equation, the saturation temperature as a function of the system pressure is calculated: C T sat=Â— (3.7) D ~ ln Psat where C and D are constants. The saturation pressure, or the vapor pressure, can be obtained from the approximate equation of state for an ideal gas, and the mass conservation in which the bulk vapor density and temperature are used: Psat=P,a, = ~RT bvap , m total = m v + mf = Constant (3.8) vap where R is the gas constant for the vapor phase. Eqs.(3.7) and (3.8), have been proposed and used as the basis for the model developed by the authors, Ding and Anghaie (1994a, b, 1995a, b, c), Anghaie and Ding (1995). 3.2 Analysis and Assumptions To simplify and generalize the formulation for the bulk evaporation and condensation processes, these governing equations and closure relationships should be normalized. Scaling and analysis of order of magnitude are conducted to develop simplifying assumptions that lead to a concise and efficient solution procedure for the governing equations. For the different phases of any material, the thermodynamic and
PAGE 38
transport properties are quite different, and they may be in the same order or may have a difference of several orders of magnitude. For the same phase, some properties are independent or slightly dependent on temperature or pressure, but others may be strongly dependent. The following relations are found about the orders and variations of these properties for the liquid and vapor phases far from the thermodynamic critical point, within a moderate phasechange range. c , = constant, c = constant, c v . C~ 0(1)i k, = constant, k % constant, k, *: = ^~O(10)O(10 2 ), (3.9) p, = constant, pÂ„ * constant, p, p = Â— ~O(10 2 )O(10 3 ), a, = constant, a, * constant, a= ail a v ~0 (10" 1 ) In Eq (3.9) ' O ' represents the order of magnitude. Therefore all the properties, except the vapor density and vapor thermal diffusivity, can be assumed constant. The liquid specific heat and vapor specific heat are assumed equal. In the following analysis and computation, the values of the ratios such as K= = ' Â° Pa ~Â— = 1 2 "Co = Â— = 1 (T ' are used, where the subscript denotes a S AYo "v.o referential state Other values of these ratios can also be used. Computer experiments reveal that for the cases with higher property ratios, that is, p = 10 3 , k 10 2 , the computation becomes more difficult to converge, especially when convection is also considered. It is due to the larger discontinuity encountered at the liquidvapor interface.
PAGE 39
As discussed in some literature, for example, Crank (1984), Ozisik (1980), Ostrach (1983), due to the density difference between two different phases, the volumetric expansion effect induces convection in the phase of lower density. Let U, denote the speed for the motion of the interface and <7 V be that for the motion of the vapor phase resulting from the density change The mass conservation at the interface should be satisfied as follows P t U,=pMU v ) (3.10) The energybalance equation at the interface is K ~^K ,^ + P^JJ, = {p^~p,e,)U l (3.11) where e v and e, denote the internal energies for the saturated vapor and liquid phases, respectively. Eliminating U, from Eq.(3. 1 1) via Eq.(3.10) yields ST. ar K >Â£K *& = P> 6eU > (312) which can also be rearranged as U, & < Â«&*l?~ S *L7KLt*> (313) where L is the thermal diffusion length scale which is the length of the container, St=cJ(/M 6 = T/To, X =x/L p =p/pÂ„. There are two possible velocity scales at the interface, U, = St(afl) and U 2 = St(aSlJ/p . Their ratio is, U,/U 2 = k/k, > 1 , thus U,>U 2 . It is known that the lower velocity scale U 2 controls the movement of the interface.
PAGE 40
25 Eq.(3. 13) can also be expressed as ., n ae, 00, U, = U, Â— {/, Â— n mi 'ax 2 ax ( ' The velocity, U 2 = St(aJL)/p, includes the density change effect which actually reveals the constraints of energy balance and mass conservation at the interface, and the volumetric expansion or contraction during evaporation or condensation, respectively. It is lower than the pure thermal diffusion velocity in liquid, a/L and that in vapor aJL. Thus U 2 is to be chosen as the characteristic velocity scale where the density change effect, as well as the St number, appears explicitly. The dimensionless velocity of the interface movement is u L= ae, ae v u~ K ax ax U,=^= kÂ— (3.15) For the case of evaporation, for example, heat is supplied to the system from the top wall which is the boundary for the vapor phase. If the appropriate amount of heat flux is provided to ensure Â— * ~ 0(1), it is found that U,~0(1), thus u , ~ U, = St(cc v I L)l p. Numerical results which will be shown in the next section to support this approximation. The characteristic time scale is L L 2 p Obviously this time scale is larger than those of thermal diffusion in vapor and Z, 2 L 2 L 2 liquid phases, that is, Â— p > Â— > Â— , which implies that the motion of the interface
PAGE 41
controls the rate of this heat transfer process. From the above analysis, it can be concluded that the motion of the interface is affected by the density difference between liquid and vapor phases, thus its velocity is very low. During the evaporation process within a constant system volume, a density change increases the vapor pressure and the saturation temperature, which in turn suppresses the evaporation. From Eqs.(3.10) and (3.15), the velocity of vapor convection induced by the density change is U, = (1 p)U t Â£~fSt , and since p Â» 1, U, (aJL)St , which is p L still very low Therefore the Reynolds number and Peclet number for the convection are Re{/Â„ /Â„ lv v =(lJL)(StlPr), Pe= (IJL) St. , where /Â„ is the vapor motion length scale. Since the velocity of the vapor convection is very low, and the convection just occurs near the interface, the vapor convection length scale will be very small compared to the thermal diffusion length scale, say, //Z, < 0(1). For most vapors, Pr ~ 0(1), and if moderate St is considered, that is, Sf ~ 0(1), then Re < 0(1), and Pe < 0(1). For the problems in which Re and Pe are less than the order of unity, the convection effect is insignificant. Therefore heat conduction dominates the heat transfer process and the convection effect induced by the density change is negligible. On the other hand, the convection effect caused by the density change during phase change can be much less significant than the convection induced by the buoyancy effect. The characteristic velocity of the natural convection is V, = (Gr)" 2 v/L where Gr is the Grashof number, the corresponding Reynolds number Re= U r L/v = (Gr)" 2 . If moderate Grashof number is considered, for example, Gr = W ~ 10 6 , then Re = 10 2 ~ IO 3 , which is much stronger than the density change effect.
PAGE 42
27 From this analysis, the convection effect caused by the density change can be neglected reasonably and the bulk evaporation and condensation processes can be assumed entirely driven by heat conduction in the absence of gravity. As a result, the governing equations (3.6) can be reduced to the the form of pure heat conduction equations: St Pf^(c VJ T) = V(K f Vr)~p f Ae^ 0 (319 > where p s and p v are dimensionless saturation pressure and vapor density, respectively
PAGE 43
3.3 The Vapor Phase Fraction Update To solve the physical model described by Eq.(3.18), a closure relationship between the vapor phase fraction / and the temperature 9 or the internal energy e should be formulated. For a pure substance undergoing evaporation or condensation, the total internal energy is a discontinuous function of the temperature: e= Â— = S,e e6 ^ where 6U is the saturation temperature. However, from a computational viewpoint, discontinuities are difficult to track and it is often necessary to smear the phase change over a small temperature range to attain numerical stability. Thus we can have e = S, 9
PAGE 44
which is used to iteratively update the phase fraction from the computed temperature field. This technique has been successfully used by Shyy (1994a, b), Shyy et al. (1994a) and was referred to as the Tbased update method. This technique has proven to be very effective for achieving convergence. An alternative formulation to the Tbased update method is to express the phase fraction as a function of the total enthalpy rather than the temperature. This update is referred to as the Hbased update method. More information about this method can be found in Shyy (1994a) and Shyy et al. (1996). A third formulation expressing the vapor phase fraction as a function of the total internal energy, rather than the temperature or enthalpy, is developed in this work. According to the classical thermodynamic definition, the relationship between the vapor phase fraction and the internal energy is f = ee, (3 23) This is used for the iterative update of the vapor phase fraction. This update procedure is referred to as the Ebased update method and is analogous to the isothermal case discussed by Voller and Swaminathan (1991). In particular, it is noted that Â£=0 is allowable in this formulation, thereby allowing for modeling phase change of a pure substance with higher accuracy. This is the most important advantage of the Ebased method compared to the Tbased method where a significant effort must be made to determine the phase change temperature interval e. For the phase change problems involving complex interface geometry, the normal temperature gradient varies with position along the interface. The determination of the phase change intervals for such problems is a tedious task The Ebased method does not need to use a phase change
PAGE 45
30 window and can be conveniently used for many phase change problems including the ones with complicated geometry. Another advantage of the Ebased method is the fact that if the latent heat is large relative to the sensible internal energy, the vapor phase fraction varies slowly with respect to the total internal energy. This enhances the stability of the update procedure. For the liquidvapor phase change occurring far away from its critical point, the latent heat is relatively large and the application of the Ebased update method even for simple geometry is recommended. In the present work, the Tbased and Ebased methods are examined in detail from the viewpoint of accuracy and computational efficiency. A straightforward implementation of the Tbased method is to discretize Eq. (3.19) for each computational cell: Â«A' = ZajJ + 6* ~CA W "A"] (3.24) where the subscript P represents the nodal point of the control volume to be solved, and nb represents the neiboring nodal points This equation is solved according to the Tbased update procedure given in Eq.(3.22). As reported by Shyy and Rao (1994) and Prakash et al. (1987), for Stefan number less than unity (large values of latent heat), this form of the iterative update is prone to numerical instability in the form of oscillations. For moderate values of Stefan numbers, the use of an underrelaxation approach for the calculation of the source term facilitates the solution process. The main cause of the oscillatory, nonconvergent behavior is the extreme sensitivity of the vapor fraction to small changes in temperature which leads to a large negative feedback effect through the source term. Substituting Eq.(3.22) into the discrete form of the conservation law and moving the
PAGE 46
resulting coefficients of 6 P into the a P term in Eq.(3.24), yields a set of stable equations similar to those formulated by Shyy and Rao (1994). These equations are: AV a F e; = xajj +b k Â—[//' //' ] **' < e, or e P k ' > a [a " +J pti ]0 r k = ^"^ + *' + ^7 [ ^ + / '" ] e > < *" < * (325) where the superscript Â« denotes the current time step and k denotes the current iterative value. Using this form of the update procedure, the temperature and the vapor phase fraction fields can be computed simultaneously and are therefore coupled at every step of the iterative process. This numerical procedure converges quickly with no loss of accuracy, granted that a proper value of e must be used. The implementation of the Ebased method is relatively straightforward and does not need any preconditioning. The internal energy is updated first, with the current values of saturation temperature and latent heat, e p k =S,0 ;+/'' 0
PAGE 47
fAxAy, while the liquid occupies (I J) AxAy. For this one dimensional case, the interface is normal to the x direction. The effective conductivity of this control volume can be obtained through a heat flux analysis. Assuming the control volume is at a steady state, that is, the interface does not move and the latent heat is not involved, the heat flux along the x axial direction is conservative Ar AT AT, tx = ~ K r = *. ~ = K, n 781 1 Ax "/Ax '(l/)Ax (S M) where K f is the effective conductivity, and the temperature is continuous AT = T w T. = (7; TÂ„)+ (TÂ„ T.) = A?; + AT, (3.29) Hence 1 (3.30) which is in a harmonic average form. According to the thermodynamic definition of the specific volume per unit mass for a mixture system, it is easy and straightforward to determine the effective density of this control volume as 1 1 P/ = ^ = vj + o/)v, = T^Z (3 ' 31) P, Pi where v is the specific volume per unit mass. It is similar to Eq.(3.30).
PAGE 48
Figure 3.2 The control volume containing two phases 3.5 Results and Discussion 3.5.1 Accuracy Aspect A test case of pure bulk evaporation is selected to examine the accuracy, and computational performance of each formulation with different update methods. A fixed volume container with one fourth of volume is filled with saturated liquid is chosen. The rest of the container's volume is filled with the saturated vapor. The initial temperature of the saturated fluid is 1.0. When the top wall is heated and kept at a temperature of 1.5,
PAGE 49
while the bottom wall is insulated. Such an orientation of heat flow is required to avoid the superheat in liquid, and the bubble generation and boiling, while the bulk evaporation is allowed. A mesh of 81x6 nonuniform and fixed grids is adopted for the liquid and vapor within the container. The time step size Az = 10" 4 is chosen. Second order central differences are used for all spatial derivative terms and fully implicit differencing is employed for time marching. Figs. 3. 3(a) and (b) show the evolution of the evaporation process. At the time instant of 0.03, about 18% of liquid is vaporized, and the saturation temperature is about 1.38. A little wavy or nonuniform variation can be seen due to the nature of the discretization and the nonuniform release or absorption of latent heat. A significant part of the latent heat content in a computational cell is absorbed when it starts to undergo evaporation. Thus the interface movement is slowed down when phase change in this cell starts and subsequently speeds up. The grids are clustered near the liquidvapor interface, and the finer grids are placed near the interface, covering the range of interface movement. Different minimum grid sizes, such as Ax mi Â„ = 1/120, 1/200, 1/400, are selected to examine the solution accuracy and the independence of solution on the grid size. The results with different minimum grid sizes are basically identical, while the finer grid size yields a smoother solution. In the following computations, the grid size with AX mi Â„ = 1/200 is used.
PAGE 50
0.005 0.01 0.015 0.02 0.025 0.03 Time (t) 0.005 0.01 0.015 0.02 0.025 0.03 Time (x) Figure 3.3 Solution features for bulk evaporation (Grids 81x8; solid line: AX min =l/120; dashed line: AX mi Â„=l/200; dashed and solid line: AX min =l/400)
PAGE 51
0.01 0.015 0.02 Time ( I ) (a) evolution of vaporization a> 1.4 1.35 1.3 Â§ 1.25 8 1.2 = jj 1.15 a C/3 1.05 C 0.005 0.01 0.015 Time ( i 0.02 0.025 0.0 (b) variation of saturation temperature 0.005 0.01 0.015 0.02 0.025 Time ( t ) (c) computation efficiency performance Figure 3.4 Comparison of two formulations for evaporation ( The solid line: the internal energy formulation, the dashed line: the enthalpy formulation )
PAGE 52
0.005 0.01 0.015 0.02 0.025 0.03 Time ( t ) (a) evolution of vaporization 0.005 0.01 0.015 0.02 0.025 0.03 Time( i ) (b) variation of saturation temperature 0.005 0.01 0.015 0.02 0.025 0.03 Time ( t ) (c) computation efficiency performance Figure 3 5 Comparison of two update methods for evaporation (The solid line: the Ebased method; the dashed line: the Tbased method)
PAGE 53
3.5.2 Comparison of H Formulation and E Formulation The same test case is computed to examine the performance of different formulations, the enthalpy formulation and internal energy formulation. Figs.3.4(a) and (b) show the evaporated liquid fraction and the variation of the saturation temperature during the evaporation process. Fig. 3. 4(c) shows the efficiency aspect of the computations, that is, the number of iterations required at each time step for the residual to go below 10" 4 . The enthalpy formulation requires a few more iterations to converge. It is found that the internal energy formulation and enthalpy formulation yield nearly identical results, with comparable computational efficiency. It indicates that both of these formulations are satisfactory for this problem and have similar performance efficiency. The only advantage in the internal energy formulation is the more concise form without an explicit pressure term. In all further investigations, only the internal energy formulation is adopted. 3.5.3 Comparison of EBased and TBased Up d ate Methods The Tbased update method and the Ebased update method are examined for the same testing case. Figures 3.5(a) and 3.5(b) show that the Tbased method and Ebased method yield nearly identical results. It is seen that convergence for the Ebased method requires about half the iterations which is needed for the Tbased method. Thus the Ebased method is computationally more efficient than the Tbased method, as shown in Fig.3.5(c). The performance efficiency of these two update methods for the evaporation
PAGE 54
and condensation problems are different from those for melting and solidification problems, where the Tbased is better than the Hbased method as reported by Shyy and Rao (1994). The reason for the improved efficiency of the Ebased method is that the phase change temperature window is used for the Tbased method, which is not required for the Ebased method. When the Tbased method is used and a purely artificial phase change window has to be adopted, it is not easy to use the temperature to identify the phase interface since the temperature is always continuous at the interface. For melting and solidification problems, the phase change temperature is usually close to a constant. The saturation temperature, however, changes continuously for evaporation and condensation processes in a constant volume system. This difference causes distinctive computational efficiency performances. For each iteration at every time step, the saturated vapor and liquid temperature as well as the saturation temperature should be updated. Therefore, the Tbased method converges more slowly. When the Ebased method is used, the internal energy has a sharp jump at the interface, and it is easy to use internal energy to identify the phases. As has been mentioned before, the relatively large amount of latent heat will result in a more stable update procedure and faster convergence rate. Furthermore, as discussed in the section of update methods, if the Tbased method is used, there will be some difficulty to determine the phase change temperature intervals for the complicated phase change problems where the interface could be of any arbitrary shape. Use of the Ebased method can completely circumvent the need for a phase change window Therefore, the Ebased method is superior to the Tbased method, in terms of computational efficiency and
PAGE 55
formulation simplicity; also it is capable of handling more complicated phase change geometry. Based on this finding.the Ebased update method is used to study the remaining 3.5.4 Bulk Evaporation and Condensation In the following work, the internal energy formulation and the Ebased update method are used to investigate the bulk evaporation and condensation processes and also to analyze the interface motion. With other conditions remaining the same, such as the initial temperature, the grid distribution and the time step size, different boundary conditions are imposed at the top and bottom walls. First, the container's top wall is heated with a constant nondimensional heat flux q M = Â— \ Xs0 = 1.0, and the bottom wall is insulated. After the time instant of 0.02, the top wall is insulated while the bottom wall is cooled with a constant heat flux q bw = k~\ Xs]0 = 10. Correspondingly, the system undergoes an evaporationcondensation process. Fig.3.6(a) shows the variation of the saturation temperature during the transient. The motion of the interface, which is respectively downward and upward during evaporation and condensation, is exhibited in Fig.3.6(b) The axial temperature distributions at different time instants during the evaporation process are plotted in Fig.3.6(c), and those of condensation process are plotted inFig.3.6(d).
PAGE 56
The velocity of the interface motion can be calculated from results presented in Fig. 3. 6(b). During the evaporation and condensation processes, the interface velocity is always constant because only the constant heat flux at either the top wall or bottom wall is provided. During the evaporation process which takes place with a constant heat flux at the top wall, the interface velocity is a constant about 2.5, which is in the order of 1.0. While during the condensation process with a constant heat flux at the bottom wall, the interface velocity still remains constant at the level of about 0.8. This indicates clearly that the nondimensional interface velocity is consistently of the order of 1.0. In other words, it can be concluded that U, . ~ 0(1), thus UÂ» ~ Ur = S, (aÂ„ / /, ) / p . This analysis verifies the phenomenological analysis presented in previous part of this paper. This also indicates that the velocity scale Ur = 5,(a v //,)/ p has been chosen appropriately so that the velocity of the interface motion can be calculated directly. 3.6 Concluding Remarks A numerical method is developed to model the bulk evaporation and condensation processes in an encapsulated container. An internal energy formulation for the constant volume phase change process is proposed and successfully implemented. Some assumptions regarding the convection heat transfer rate at the vaporliquid interface are made through scaling and order of magnitude analysis. The velocity and time scales of the interface movement are obtained through scaling analysis Depending on the configuration and boundary conditions, the convection effect induced by the density change is analyzed
PAGE 57
and found negligible. The dominating mechanism for the heat transfer at the interface is conduction. 0.01 002 0.03 0.04 0.05 0.06 07 Time(x) (a) Variation of saturation temperature during evaporation and condensation 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Time (x ) (b) Motion of interface during evaporation and condensation Figure 3.6 Bulk evaporation and condensation
PAGE 58
0.10.20.30.40.50.60.70.80.9 1 Axial position X (c) Temperature distribution during evaporation 0.10.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Axial position X (d) Temperature distribution during condensation Figure 3.6 ( continued ) Bulk evaporation and condensation
PAGE 59
The computational performances of the wellestablished enthalpy formulation and the proposed internal energy formulation are evaluated and compared. Both formulations yield identical results with similar computational efficiencies. The internal energy formulation has the advantage of having a more compact form. The performances of two update methods for the vapor phase fraction, the Ebased and Tbased methods, are investigated. The Ebased method is superior to the Tbased method in terms of computational efficiency and formulation simplicity for evaporation and condensation problems. The Ebased method is also capable of handling more complicated phase change problems, because it does not require the artificial phase change window. The internal energy formulation and the Ebased method are successfully implemented and used to compute bulk evaporation and condensation processes under different conditions. The evolution of the bulk evaporation and condensation processes is simulated. The variations of the temperature distribution, the liquid content and the saturation temperature with time for a variety of boundary conditions are obtained. Further effort will be highly worthwhile and strongly recommended to extend the current model for liquidvapor phase change processes to a unified threephase model. All phase change processes between these three phases, that is, solidliquid, liquidvapor, solidvapor, or solidliquidvapor, could be simulated.
PAGE 60
CHAPTER 4 BULK EVAPORATION/CONDENSATION WITH INTERNAL HEATING AT ZERO GRAVITY Based on the modeling procedure developed in the previous chapter, the bulk evaporation/condensation processes with internal heat generation at zero gravity level are studied. The system undergoing the phase change process and at a steady state is shown as Figure 4. 1 . Since the system is designed to operate under a 0g environment, there are no buoyancy effect and no natural convection. Initially, the vapor and liquid are kept saturated at a certain pressure. The top wall and bottom wall are insulated. The side wall is cooled with a constant heat flux. The internally generated heat evaporates the liquid, while the cold side wall removes the heat from the liquidphase and condenses the vaporphase. For this system with such an orientation and operated under 0g, there is no bubble or droplet generation on the wall surface. The temperature along the liquidvapor interface does not change significantly, because the curvature of the interface is moderate. The liquidvapor interface can be reasonably assumed to be isothermal, consequently the thermal capillary effect can be neglected. Hence the phasechange process is entirely driven by the internal heat generation and heat conduction. The rate of evaporation or condensation is determined by the local temperature distribution The time dependent and nonuniform internal heat generation results in a transient temperature distribution and causes the motion and deformation of the interface.
PAGE 61
insulated evapqration liquid ; O gm ,i cooled liquid insulated (a) schematic diagram of the phasechange (b) at the final steady state process Figure 4. 1 Schematic diagram of the twophase cell operated at 0g
PAGE 62
4. 1 The Governing Equations The internal energy formulation of the energy equation for this phasechange process with the constraint of constant volume is used. The weak form of the energy equation for conduction is />Â§ = v(Â«vr)+e^, (4.i) The source term Q ge Â„ is the internal heat generation rate which varies with time and space according to the following relation: fix Q sm = pB(t)sMÂ—) ( 4.2) where B is a time dependent parameter, p is the local phase density, L is the length of the container. The internal energy e for different phases is also given by Eq.(3.5). Substituting Eq (3.5) into the governing Eq.(4.1) yields P/fM = VOt,VT)p / As. + ^ (4.3) Pvfe,r) = V(^ v vr) + (2 s Â«Â„, v
PAGE 63
where the latent heat now appears as a source term. The vapor phase fraction is defined , ee. This is the Ebased update method discussed in chapter 3 and used for the iterative update of the vapor fraction. The Tbased method does not work well for this twodimensional bulk evaporation and condensation problem, while the Ebased always works quite well. Some additional closure relations, the same as Eqs.(3.78) and (3.19), are required for the multiphase condition in the container. Some assumptions and simplifications have to be made based on the order and scaling analysis as discussed in chapter 3. All the properties, except the vapor density and vapor thermal diflrusivity, are assumed constant. To simplify the computation, the liquid specific heat and vapor specific heat are assumed equal. In the following computation, these values of the ratios such as k = %= 10, p = & = 1 00 are used where the *v A..0 subscript denotes a referential state. The governing equation (4.4) is normalized as: 09 (4.5) where St is the Stefan number defined as cJJAe. The nondimensional heat generation rate is ri i \z gen / = o 08 _ , 1 cf Pf St&z 0
PAGE 64
49 Â£U, = P,G(z)sm{KX) (4.6) where G(t) = B(t)(tr/c v T sQ j. The initial saturation temperature T sM and the length of the container are chosen as the characteristic temperature and length scales, respectively. The characteristic time scale is tr = L 2 /a, which is the thermal diffusion time scale in the vapor phase. 4.2 Conductivity at Control Volume Surfaces The finite control volume difference method is employed to discretize the differential equations. For different materials, that is, the solid wall and the fluids in the container, their thermal and transport properties such as densities, heat conductivities, are quite different. Due to the nature of this discretization, the interface between two different media is always located on the control volume surface. This arrangement requires that the effective properties at the interface, namely the control volume surface, be determined. Patankar (1980) discussed this problem and presented the determination of harmonic mean effective conductivity at the interface for the rectangular coordinate system. Here in this work, the general determination of the conductivity at the surfaces of any control volume whether it contacts to an interface or not, for a cylindrical coordinate system, is developed through the heat flux analysis. As shown in Figure 4.2, P is the nodal point of the control volume, where the value of/ at P is to be calculated. E, W, N and S are nodal points surrounding the control volume, and e, w, n and s represent the surfaces enclosing this control volume.
PAGE 65
5 I \\ \\\\NN\\\\\\ \ \\\\\\ \ \ / w ,p / : / / 17 / / /!"// / / / / / / / , / f // '///. t/ s (5x) e _ Figure 4.2 The control volume contacting other material Considering the two dimensional, axisymmetrical coordinates, the linear heat flux along radial direction is and also q, = IntKÂ— = Ijik ( Tp T ^ dr "\a(r N lr P ) lnfe/zJ ln(rÂ„//>) Using the relation (4.7) (4.8) (4.9) the heat conductivity at the north surface, tcÂ„, is
PAGE 66
K K K p ln(r N lr F ) K n = aÂ— t Â— ^ Â— ^ (4 101 K p \n{r N lr n ) + K N \n(rJr p ) The determination of heat conductivity at the east surface is simpler and more straightforward (see Patankar, 1980), ,c e = ( ^A + Ay (411) fC p K E where / e is a ratio defined /,S (*)1*W (412) in terms of the distances shown in Figure 4.2. The heat conductivity at the south and west surfaces could be obtained similarly. This approach used to determine the transport properties at the control volume surfaces provides consistent relations for properties that can be applied for any interfacial boundary, and results in higher accuracy with low computation cost. 4.3. Results and Discussion The second order central difference scheme for the spatial differential terms and the fully implicit marching scheme for the temporal term are adopted to solve Eq. (4.5). A uniform mesh of 101x41 grids covers the liquid phase and vapor phase within the container, and another non uniform mesh with 101x5 grids covers the side wall region. The outer surface of the side wall is cooled by a coolant flow at constant flux which is negative unity.
PAGE 67
52 Because of symmetry, only the right half of the container is considered. Initially, the lower half of the container is filled with the liquid phase and the upper half is occupied by the vapor phase. The aspect ratio of the container is 2.5. The initial saturated vapor and liquid temperature is 1 . At l = 0, the side wall is cooled at a constant heat flux, and the internal heat generation begins. The time dependent parameter for the heat generation rate expressed in Eq.(4.7) is given by Eq.(4. 13) and shown in Figure 4.3. G(r) = 12r G(z)= 12x1.5= [I 01.5 (4.13) G(T) ' / jT r Figure 4.3 Parameter for heat generation rate The evolution of the evaporation/condensation process is shown in Figure 4.4. Using the phase identifier which is the local vapor phase fraction, the distribution of liquid and vapor phases and the position of the interface are determined. The topology of the liquidvapor interface and the temperature distribution at several time instants are
PAGE 68
displayed in Figure 4.4 and Figure 4.5, respectively. The internal heat generation evaporates the liquid, while the side wall cooling tends to condense the vapor. Globally, the evaporation of liquid and condensation of vapor take place simultaneously, as shown schematically in Figure 4.1. The local phasechange process, whether evaporation or condensation, is dominated by the local internal heating or by the wall cooling effect, since the internal heat generation is a spatial function. The evaporation is dominant near the central axial position where the heat generation intensity is higher and evaporates more liquid. The condensation caused by the side wall cooling is dominant near the axial end and side wall since the heat generation density is lower there. The simultaneous evaporation and condensation cause motion and deformation of the interface. Because of the nature of time dependent heat generation, the local evaporation and condensation, the motion and deformation of interface continue. As the process evolves, a certain amount of liquid is evaporated into vapor due to the internal heat generation which increases at the early time period, as Eq.(4.13) indicates. Meanwhile, some amount of vapor near the side wall is condensed into liquid due to the cooling of the side wall. At a very earlier time period, a samll portion of liquid film is formed along the side wall owing to the condensation, and the liquidvapor interface becomes concave because of the evaporation. With increasing time, the interface keeps moving and deforming. The liquid film grows longer and thicker, and the position of the interface at the centerline of the container descends lower and lower. When z = 1.0, the liquid film grows up to the top end and covers the entire side wall. The liquid phase which is initially filled in the lower part of the container now forms a long film covering
PAGE 69
54 and residing along the entire side wall. The liquid film becomes thicker and thicker at the top end where the internal heat generation is zero. At r =1.5, the interface at the low central side descends quite close to the bottom end. After r "1.5, the internal heat generation rate is a constant as indicated by Eq. (4. 13) and the cooling heat flux at the side wall is still kept constant, the system reaches a steady state and the interface no longer moves. The topology of the interface strongly depends on the temperature distribution and the internal heat generation rate distribution which now is independent of time but spatially dependent. The shape of the interface approximately follows the heat generation distribution in the form of chopped sine function. During the transient process, the temperature distribution varies with time The temperature contours for the entire system gradually become denser and denser, since the internal heat generation is increasing. After steady state is reached, the axial temperature distribution follows approximately the form of chopped sine function, due to the nature of the heat generation distribution as shown in Eq. (4.2). There is a temperature gradient jump at the interface due to the release or absorption of latent heat associated with phase change. The discontinuity of the temperature gradient can be observed more clearly in a three dimensional temperature distribution displayed in Figure 4.6. There is an obvious fold on the temperature surf where the liquidvapor interface is exactly located. The axial temperature profile on the side wall is also in an approximate form of chopped sine function, but much flatter than that in the vapor region. The vapor is well superheated and its temperature gradient is greater than that of the liquid, because the liquid is a much better conductor than the vapor. This isothermalization of the liquid near the side wall is
PAGE 70
55 expected to satisfy the design feature of the multiphase fuel cell to be possibly applied in space. The formation of the liquid film and its full coverage on the side wall can prevent the side wall from being overheated, which is the expectation for the application of such multiphase fuel cell. The transient characteristic of this system can also be seen from the variation of the saturation temperature with time shown in Figure 4.7. The saturation temperature decreases a little bit in the very early time period when the internal heating is very weak and the side wall cooling and condensation are dominant. As the internal heating rate increases, the evaporation becomes dominant and the saturation temperature rises. After r = 1.5, the system gradually reaches steady state and the evaporation and condensation are balanced and the saturation temperature maintains constant. The operation in a zero gravity environment is ideal and unrealistic However, the study conducted in this chapter provides a basis for further exploriation of the fundamental performance and features of such a multiphase and internally heated system. The present results will be used as reference for comparison with the performance and features of this system operated under microgravity conditions, which will be presented in chapter 6.
PAGE 71
vapor liquid / 7=0.2 0.4 0.6 0.8 7=1.0 1.2 1.4 1.6 Figure 4.4 The evolution of liquidvapor interface (0g)
PAGE 72
51 r=0.2 0.4 0.6 0.8 r=1.0 1.2 1.4 1.6 Figure 4.5 The evolution of temperature field (0g)
PAGE 73
Figure 4.6 Temperature surf of the system at steady state (0g) 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Time ( x ) 1.8 Figure 4.7 Variation of saturation temperature (0g)
PAGE 74
CHAPTER 5 NATURAL CONVECTION WITH INTERNAL HEATING IN SINGLEPHASE AND TWOPHASE SYSTEMS Like all phase change processes in a finite gravitational field, evaporation and condensation are inevitably subject to the effect of convection induced by buoyancy. The existence of convection can substantially affect the temperature field, and subsequently the performance of the phase change process. Internal heating, like external heating, results in temperature and density gradients, and results in a buoyancy effect in nonzero gravity. Such convection induced by internal heating is defined as internal natural convection which wull be investigated. Natural convection is involved in both liquid and gas phases. The coexistence of natural convection in both phases causes interaction of the convection cells between two phases. In this chapter, some natural convection problems with internal heating in the singlephase and twophase systems are investigated. A benchmark problem, natural convection in a square cavity of singlephase without internal heating, is computed. Comparisons of the present solutions with accurate published benchmark solutions are made to verify the computational accuracy. The natural convection in a square cavity of singlephase with internal heating is also computed. Qualitative comparisons with the published experimental results are also made. Then natural convection in a liquidgas system without phasechange, with or without internal heating, is investigated. The interactions between these convection cells in both phases are revealed. 59
PAGE 75
60 5.1. Extern al Natural Convection in a SquareCavity Benchmark Problem Considerable attention has been given to the study of natural convection in enclosures. The temperature gradients supplied externally across the system induce the density gradients and consequently the buoyancydriven flow. A great deal of research, both numerical and experimental, on this classical type of convective heat transfer process has been conducted, Ostrach (1982, 1983, 1988), Gebhart et al. (1988), de Vahl Davis (1983), Hortmann et al. (1990), Shyy et al. (1992b), Shyy (1994a), Bejan (1984), Carpenter et al. (1989), Ramaswamy el al. (1992), Bergman et al. (1981), etc.. Ostrach (1982, 1983, 1988), Gebhart et al. (1988) and Bejan (1984) reviewed and summarized those activities and achievements. Among the large number of natural convection problems, buoyancydriven flow in a square cavity with vertical side walls heated differentially, a socalled benchmark problem for natural convection, is a suitable vehicle for testing and validating computer codes. Accurate benchmark solutions have been published by de Vahl Davis (1983) and Hortmann et al. (1990) and are available for comparison with the present computations. The problem being considered is that of twodimensional flow of a Boussinesq fluid of Prandtl number 71 in an upright square cavity. The horizontal walls are insulated, and the vertical sides are at temperature T h and T c . Figure 5. 1 shows the domain and boundary conditions of the problem in schematic form.
PAGE 76
61 Â® Â„ = 0,Z/= V= <$/ T=T h T=T c u = v = g u = v = " ar = 0,!< = V = * Figure 5. 1 Domain and boundary conditions for the benchmark problem The governing equations in dimensionless form are + Â— = ax at UÂ— + V~r = + PrV 2 U ax ay ax ax ar ar + PrV 2 V + RaÂ„Prd (5.1) (5.2) (5.3)
PAGE 77
ax at K ' where the external Rayleigh number is defined as (5.5) ^tfazag and /? is the thermal expansion coefficient which is used in the Boussinesq approximation p = Po [\P(TT )\ (5.6) In obtaining the above equations, the following dimensionless variables are used X = x/L, Y = y/L, U = u/(ct/L), V = v/(o/L), = (TTJ/(T h T c ) (5.7) where (a/L) is the thermal diffusion velocity scale. Other velocity scales, such as buoyancy velocity scale, etc., can also be used to yield different forms of the governing equations. Solutions have been obtained using a secondorder centraldifference scheme for the momentum and energy equations. Two grid sizes, 41x41 and 81x81, are adopted for cases ofRa E = 10 3 and Ra E = 10 6 , respectively. Two external Rayleigh numbers, 10 3 and 10 6 , are chosen for the calculations. Figure 5.2 shows the steady state computed results. An important feature of the solution is the skewsymmetry of the velocity and the temperature field with respect to the center of the cavity. Fluid is heated at the left wall, decreases in density and rises up, whereas fluid adjacent to the cold right wall loses heat and increases in density and falls down. This establishes a clockwise recirculating flow inside the box. With increasing Rayleigh number, it can be observed that the convection strength increases and the isotherms are nearly horizontal in the core flow. The heat transfer rate is observed to be enhanced with increasing Ra as can be deduced from the isotherm spacing at the left and right walls.
PAGE 78
Increasing the Ra increases the vorticity generation in the flow field which in turn increases the convection strength and results in an enhanced overall heat transfer rate through the domain. The enhancement of the heat transfer rate with the increasing Ra can be observed more clearly in Figure 5.3 which shows the Nusselt number distribution on the hot and cold walls. The Nusselt number on the wall is defined as follows hi 00 , Nu(5.8) k 3X '"""' which actually is the nondimensional temperature gradient or heat flux at the wall. The present results are compared with the benchmark solution of de Vahl Davis (1983). This comparison is shown in Table 5.1 for Ra E = 10 3 , 10 6 . As may be observed from this table, the present solutions agree quite well with the benchmark solutions. Table 5 1 Comparison of present results with benchmark solutions ofde Vahl Davis (1983) Ra=10 3 (41 x41Grid) Ra=10 6 (8 1x81 Grid) Davis (1983) present discrepancy Davis (1983) present discrepancy W>Â«id 1.174 1.167 0.5% 16.32 16.08 1.5% Vmrnax XJ ~ ~ 16.75 0.151, 0.543 16.57 0.151, 0.540 1% Umax, Y 3.634 0.813 3.629 0.813 0.1% 64.63 0.855 64.3 0.852 0.5% Vmm X 3.679 0.178 3.674 0.178 0.1% 219.36 0.0379 219.1 0.0375 0.1% Nu , ma Y 1.505 0.092 1.515 0.088 0.6% 17.92 0.0378 17.38 0.0375 2.7% Y 0.692 1 0.693 1 0.1% 0.989 1 1.023 1 3.3%
PAGE 79
64 (a)Jto= 10 3 , Pr = 0.71 (41x41 Grid) Ay/=QA 4(9=0.05 (b) Ra = 10 6 , Pr = 0.71 (81x81 Grid) Ay/=\ zl0=O.O5 Figure 5.2 Flow patterns (left) and isotherms (right) of buoyancy driven natural convection
PAGE 80
0.6 0.7 0.8 0.9 1.1 1.2 1.3 1.4 1.5 1.6 Nusselt number distribution on hot and cold walls (External natural convection with Ra=103, Pr=0.71) 6 8 10 12 Nu(Ra=10 6 ) 14 16 18 Figure 5.3 Nusselt number distribution on side walls
PAGE 81
66 An examination of the accuracy and dependence of the solution on the grid size is performed. The temperature profile at the horizontal centerline of the Ra = 10 6 case is chosen to be tested with different grid sizes, for example, AX = 1/20, 1/40, 1/80 and 1/120. Figure 5.4 shows the numerical solutions with AX = 1/40, 1/80 and 1/120 are almost identical and the solution with AX = 1/20 is not quite accurate. With AX < 1/40, the solutions are spatially accurate and independent of grid size for the benchmark problem with Ra = 10 6 . 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 5.4 Temperature profile at horizontal centerline (Ra = 10 6 ; Solutions of different grid sizes; solid lines : 121x121 and 81x81, dots: 41x41, dashed line Â— : 21x21)
PAGE 82
5.2 Natural Convection in a SquareCavity with Internal Heating Natural convection induced by external heating is of importance in many applications. Natural convection in an enclosure in which the internal heat generation is present is also of prime importance in certain technological applications. Examples are postaccident heat removal in nuclear reactors and some other geophysical or astrophysical problems. Natural convection induced simultaneously by external heating and internal energy sources, however, in contrast to the extensive investigations on the conventional external heating convection problems, has received very limited attention, for example, Archarya and Goldstein (1985), Lee and Goldstein (1988), Fusegi et al. (1992). The presence of internal heat generation provides an additional dynamic in overall convective flow systems and offers a complicated yet challenging aspect of the convection under consideration. The introduction of internal heat generation brings forth an additional nondimensional parameter Ra h which was conveniently termed as internal Rayleigh number in Archarya and Goldstein (1985), Lee and Goldstein (1988), Fusegi et al. (1992). This is in addition to the conventional Rayleigh number, Ra E , which is based on the externally imposed temperature difference. The work described in this section is to explore the effect and significance of the internal heat generation induced convection and qualitatively verify the results with previous investigations such as Archarya and Goldstein (1985), Lee and Goldstein (1988), Fusegi et al. (1992). In addition, further details of flow and heat transfer characteristics
PAGE 83
68 will be displayed. The physical situations, the domain and boundary conditions, to be considered are shown in Figure 5.5. Two boundary conditions on temperature are considered. One is that all the walls are kept at temperature T w . Another one is that the top and bottom walls are insulated while the side walls are maintained at T h and T c . The continuity and momentum equations are the same as Eqs (5.1)(5.3). The only difference is the appearance of an additional term containing Ra, in the energy equation as shown here, tf* +r *.V*+A (59) dX. m Ra E y ' where Ra E , the external Rayleigh number, is still the same as defined in Eq.(5.5) and the internal Rayleigh number, Ra, , is defined as auK where G is the volumetric internal heat generation rate. The results are shown in Figures 5.65.8. For case 1 where all walls are maintained isothermal at 7* w , the natural convection is purely induced by the internal heat generation. The flow field and temperature field are all symmetric about the vertical centerline due to the uniform temperature along all walls and the uniform internal heat generation. In the presence of internal heating, a pair of counterrotation convective rolls is formed. The flows move upward from the central lower region, divides and move downward separately and symmetrically along both side walls. The convection induced by the internal heating shifts the isotherms toward the top wall so that the position of the maximum temperature is higher than that of pure conduction. When the internal heating is not strong, that is, Ra,
PAGE 84
< Ra E , the convection effect is not significant. The isotherms are closely symmetric about the center of the cavity. As Ra, increases, the position of T max moves higher When Ra, is increased to 10 6 , there are two extreme values of temperature in the fluid and two extreme values of temperature gradient on the top wall. The flow pattern and temperature distribution agree with the experimental results of Lee and Goldstein (1988) qualitatively quite well. As the Ra, increases, the flow velocity increases, and the strength of the convection is increased. The heat transfer is increased significantly as the internal heating increases, which is displayed in Figure 5.8. In case 2 the characteristics of natural convection, as shown in Figure 5.7, are quite different since both external and internal heating exist. In the presence of internal heating (Ra, > 0), the aforementioned trends in flow and heat transfer persist if Ra, < Ra E . However, when the magnitude of Ra, considerably exceeds that of Ra E ( e.g., Ra E 10 3 , Ra, 10 5 ) the flow field takes on a different pattern (Figure 5.7 (a)). It now consists of two counterrotating convective rolls spanning the entire length normal to the adiabatic boundaries. The maximum dimensionless temperature occurs between the hot and cold walls which explains why the fluid rises in the central region of the cavity and sinks along both the hot and cold walls. With fixed Ra,, as the Ra E is increased to 10 4 , the flow field varies significantly. Three eddies are present, as shown in Figure 5.7 (b). The largest eddy is the counterclockwise flow moving downward along the cold wall and upward along the lower portion of the hot surface. The other two eddies are formed at the topright corner and bottomright corner of the box with the flow directed downward along the hot wall. The eddy at the topright corner is faster than the one at bottomright corner although
PAGE 85
70 they are both sluggish relative to the largest eddy. The flow pattern and temperature distribution agree qualitatively very well with the numerical results of Archarya and Goldstein (1985) under the same conditions. T = T w ,u = v = T=T w u = v = % T=T w u = v = 1 Â•= Â£,K = V = case 1 (a) all walls at T w = 0,!/= v= T= T c u = v = g T=T h u = v = t W Â— = 0,ii = v = V case 2 (b) top and bottom walls insulated, left wall at T c and right wall at T h Figure 5.5 Domain and boundary conditions for internal heating induced convection
PAGE 86
71 (a) Ra, = 10 4 ( left: Ay/ = 0.04; right: A6= 0.01) (b)Ra, = 10 5 ( left: Ay/ = 0.4; right: A0= 0.07) (c)Â«a ; = 10 6 ( left: 4^=1; right: ^fl=0.5) Figure 5.6 Flow pattern (left) and isotherms (right) of internal heating induced convection ( case 1 : all walls are at T w , Ra E = 10 4 , Pr = 0.71)
PAGE 87
72 Ai//=OA (a)Ra E = 10 3 A9=\ Figure 5.7 Flow pattern (left) and isotherms (right) of internal heating induced convection ( case 2: left wall at C =0, right wall at 6 h 1, Ra, = 10 5 , Pr = 0.71)
PAGE 88
73 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 A ~~ftj,^J0 6 _ \Ra,= 10 5 \ \ _Ra,= \0 4 / / / / y s s' iÂ— T 1 1 1 1 1 10 15 20 25 30 35 Nu on side wall (Ra E =10 4 ) 40 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 cold wall Ra E =K 80 60 40 20 20 40 60 80 Nu on hot and cold walls (Ra j=10 5 ) Figure 5.8 Nusselt number distribution on the side walls ( Top: case 1; Bottom: case 2)
PAGE 89
74 5.3 Natural Convection in LiquidGas Systems Without Internal Heating In liquidgas systems, natural convection takes place in both phases. The natural convection in one phase may affect the convection in another phase. This causes the interaction between the convection cells in liquid phase and gas phase, and consequently affects the thermal performance of the system. The interaction between these convection cells also makes the computation more difficult to converge. In this section and the next section, a twophase system without phase change, with or without internal heating, is considered in order to study its convection and heat transfer characteristics. These investigations have not been found in the literature. New findings are obtained through numerical modeling experiments. A rectangular cavity, with an aspect ratio of 2, filled with liquid and gas is shown in Figure 5.9 (a). The liquid and gas occupy the bottom and top halves of the cavity, respectively Both the top wall and the bottom wall are adiabatic. The ratio of the liquid density to the gas density is assumed to be 10, the conductivity ratio of liquid to gas is assumed to be unity. The Prandtl numbers for both phases are assumed to be 1. The Grashof number for the gas phase, Gr e which is based on gas properties, is 10 3 . The governing equations for different phases in nondimensional form are SpJU dpâ€¢ Â„ ^T + ^F = ( 511 ) d(pXJU) d(p,UV) dp 1 sx ft ex JGr
PAGE 90
75 d(p,lJV) d(p,W) dp 1 _ ____ =_ Â„ i + J ^r i = ^+^v.( ft vF) + /; iA / (5.13) ax fff dl pr s ^ J+ ^^=?^^;^ (514) The forms of Eqs (5.11) (5.14) are different from those of (5.1) to (5.4), and for now the natural convection velocity scale in the gas phase is given by ", = yV^ (5.15) and is used instead of the heat diffusion velocity scale, where v g is viscousity for gas phase. This formulation has the advantage that all the variables can be maintained order 1, which may benefit the computation. The second order central difference scheme and a mesh with 81x41 grids are used. 5.3 1 Gas Convection Induced bv Liquid Convection The upper side walls, contacting the gas phase, are insulated. The lower side walls covering the liquid phase are maintained isothermal with T, = 2 and T r = 1 . Only the liquid phase is heated externally. Figures 5.9 (b) through (d) show the convection flow patterns in the liquid and gas phases, and the temperature distribution in the entire system, respectively. In the liquid phase, a single clockwise rotating convection cell named as cell 1 is formed, moving upward on the left side and downward on the right side. Since no heat is provided from the side walls to the gas phase, no active convection which is directly
PAGE 91
76 induced by buoyancy differences is expected to appear in the gas phase. Two counterrotating convection cells, however, are found in the gas phase, as shown in Figure 5.9(c). A counterclockwise convection cell, named cell 2, appearing on the right side of the gas phase is obviously induced by the convection cell 1 in the liquid phase, through the shear force interaction at the phase interface. Consequently, cell 2 drives another relatively weaker clockwise cell on the left side of the gas phase, which is cell 3. The natural convection cells 2 and 3 in the gas phase are induced, through the shear interaction with convection cell 1 in the liquid phase Therefore convection in the gas phase is purely passive. The intensities of passive convection cells, cell 2 and cell 3, are much less than that of cell 1, therefore the flow pattern in the liquid phase is influenced slightly. Due to the weak convection effect, the heat transfer through the gas phase is very weak which can be seen from the isotherms. The heat transfer in the liquid phase is strongner. 5.3.2 Both Phases Heated In this case, the entire side walls covering both liquid and gas phases are maintained isothermal, that is, T, = 2 and T r = 1 . Other conditions are the same. The convection flow patterns in the liquid and gas phases, and the temperature distribution in the entire system are shown in Figures 5.10(a), (b) and (c), respectively. The temperature difference between the two side walls induces natural convection in the gas phase, as well as in the liquid phase. Three major convection cells appear. A
PAGE 92
77 single clockwise cell in the liquid phase is also named cell 1, and the other two counterrotating cells in the gas phase are called cells 2 and 3, as treated in the previous case. The counterclockwise rotating cell 2 in the gas phase is also induced by the convection cell 1 in the liquid phase. Cell 3 is caused by the temperature difference between side walls and the shear effect of cell 2. Thus cell 3 is generated both actively and passively. Due to the duel effect, the strength of this cell 3 is greater than that purely passive cell 3 of the previous case. On the other hand, cell 3 also affects cell 2, thus the convection of cell 2 is enhanced a little bit. The interactions of these convection cells in both phases are clear in this case. The intensities of passive convection cells 2 and 3, are much less than that of cell 1, therefore the flow pattern in the liquid phase is not changed significantly. Due to its weak convection strength, the heat transfer in the gas phase basicaly remains conduction characteristics. The temperature distribution in the liquid phase remains significant characteristics of natural convection because of its stronger convection strength. Another similar case, where the Grashof number for the gas phase is increased to 10 , is examined but not shown. The main features of the flow patterns are similar. Three convection cells are formed, two in gas phase and one in liquid phase. The temperature distribution is changed. Because of the higher Gr v , the convection and heat transfer are enhanced. The interactions between these convection cells are also increased significantly. The enhanced interactions between these convection cells cause severe difficulty for the computational stability and convergence. An unsteady solution procedure, which is not required for the previous cases, has to be employed to obtain convergent results.
PAGE 93
78 insulated insulated (a) Schematic of the cavity (b) Convection in liquid c) Convection in gas (d) Isotherms Figure 5.9 Natural convection of a twophase without internal heating (case 1: upper side wall insulated, lower side wall isothermal, G/g =10 3 )
PAGE 94
7 a (a) (b) Convection in liquid: Convection in gas (c) Isotherms Cell 1 : + Cell 2:; cell 3: + Figure 5.10 Natural convection of twophase without internal heating ( case 2: both side walls isothermal, Gr g =\0 3 )
PAGE 95
5.4 Natural Convection in LiquidGas Systems with Internal Heating The flow and heat transfer characteristics of a twophase system without phase change but with internal heating is investigated in this section. Two cases, one rectangular cavity, one cylindrical container, are considered. 5.4.1 Case 1: Rectangular Cavity The rectangular cavity with aspect ratio of 2 is filled with liquid and gas. The volume ratio of liquid to gas is 1. The energy equation, where internal heat generation is involved, is Other equations retain the same forms as Eqs. (5.11) ~ (5.13). In Eq.(5.16) the constant Gr gl is the internal Grashof number The volumetric heat generation rate is proportional to its mass density and is uniformly distributed in each phase. The top wall, bottom wall and the left side wall are insulated. The right side wall is kept isothermal with T r = 1. The density ratio of liquid to gas is 10, and the conductivity ratio of liquid to gas is unity. The external Grashof number based on the gas properties is 10 3 . The internal Grashof number is 3 16 The results are shown in Figure 5.11. Three convection cells are formed, one in the liquid region and two in the gas region. Similar to the previous cases, one passive convection cell on the right side in the gas phase is induced by the convection cell in the
PAGE 96
81 liquid phase. Another clockwise convection cell in the gas phase is induced by the internal heating and the shear of the counterclockwise cell. These convection cells interact and affect each other. The intensity of the convection in the liquid phase is much greater than that in the gas phase due to the higher mass density and higher heat generation rate in the liquid phase. The convection heat transfer in the liquid region is much stronger than that in the gas region where heat diffusion is dominant. 5.4.2 Case 2: Cylindrical Container The results for the case of the cylindrical container are displayed in Figure 5.12. The container with aspect ratio of 2.5 is filled with half liquid and half gas. The forms of the governing equations are similar to those of Eqs. (5. 1 1) ~ (5. 13) and (5.16), except that the cylindrical coordinate system is used here The top wall and bottom wall are insulated, and the side wall is isothermal with T r = 1. Only the right half of the container is considered owing to symmetric boundary conditions, and the axial centerline of the container is adiabatic. The external gas propertybased Grashof number is 10 5 , and the internal Grashof number which is also based on gas properties is 1.58xl0 6 . The main flow patterns and heat transfer characteristics are similar to those of the previous case. Three convection cells are generated, one clockwise rotating cell in liquid and two counterrotating cells in gas. Since the centerline of the container is adiabatic, heat transfer and convection in the region close to the centerline are quite weak. The natural convection and heat transfer are relatively strong in the liquid region near the side
PAGE 97
82 wall. In the liquid phase, the convective heat transfer is dominant. The heat diffusion is still quite significant in the gas phase due to the weak convection. (a) (b) (c) Convection in liquid: Convection in gas Isotherms Cell 1: + Cell 2:, cell 3: + Figure 5.11 Natural convection of twophase with internal heating (case 1: cavity, left wall adiabatic, right wall isothermal, Gr gE =10 3 , Gr gI =316 )
PAGE 98
( a ) (b) (c) Convection in liquid: Convection in gas: Isotherms Cell 1 : + (( / ma ,= 2.0xl04 V Â„ m = 8.4xl(r 4 Ay/=SxW 5 Cell 2:; cell 3: + Figure 5. 12 Natural convection of twophase with internal heating (case 2: cylindrical container, side wall isothermal, Gr gE =10 5 , Gr gI =1.58xl0 6 )
PAGE 99
CHAPTER 6 BULK EVAPORATION/CONDENSATION WITH INTERNAL HEATING AT MICROGRAVITY AND NORMAL GRAVITY LEVELS 6.1 Brief Preview So far, only conduction dominated effects under limited conditions have been taken into account for the bulk evaporation and condensation processes. Like all phase change processes, evaporation and condensation are also inevitably subject to the effect of convection induced by either buoyancy or surface tension variations. The existence of the convection can substantially affect the flow field, the temperature field, and subsequently the performance of the phase change process. Numerous investigations on phase change with convection have been conducted, which can be found in Ostrach (1983), Crank (1984), Brice (1986), Brown (1988), Gau and Viskanta (1986), Nadarajah and Narayanan (1990), Shyy etal. (1991, 1992a, b, 1994a, 1995, 1996) and Shyy (1994a, b, 1995), etc. In this Chapter, the modeling of the bulk evaporation and condensation processes, associated with internal heat generation and under various gravity levels, is presented. The system to be modeled is shown in Figure 1.2 similar to that of Figure 4.1. Initially the liquid and vapor are assumed to be saturated at the liquidvapor interface, while it is also assumed that the liquid is subcooled and the vapor is superheated in their own regions, respectively._The top wall and bottom wall are insulated. The side wall is cooled externally by a constant outward heat flux which removes heat from the system and tends to 84
PAGE 100
85 condense the vapor phase. The internal heat generation provides heat to evaporate the liquid phase. Since the system is to be operated under different gravity conditions, there is convection, both in the liquid and vapor phases, induced by the buoyancy force due to the density variations. The temperature along the liquidvapor interface does not change significantly, because the curvature of the interface is moderate. The liquidvapor interface can be reasonably assumed as isothermal, thus the thermal capillary effect can be neglected. The phase change processes are controlled by conduction and convection. The performances of the phase change can be quite different at various gravity condition, as can be expected The flow field, the temperature distribution and the motion of the interface are strongly dependent on the glevels, as will be demonstrated later. 6.2 The Governing Equations The equations governing the thermal energy transport and phase change processes, with convection and heat generation, are similar to Eq. (4.4), except for the additional convection terms. The equations describing natural convection in both phases are similar to those of Eqs (5.1) (5.3), while the vapor density in a certain initial state is chosen as a reference basis for both phases. The governing equations in cylindrical coordinate system for different phases in dimensional form are a L + V.(p,F) = (6i)
PAGE 101
86 dp.u dp ^+V.( A KÂ«) = ^ + V.( ft V M )A g (6.2) ^ + V.( / l F , v) = J + V.(/i i Vv) (6.3) J ^ + V Â»( P yc,T)= K ,V 2 T + p,G(t)An(Kj)he^^heV Â»( P yf) (6.4) where the subscript ;=/ and v, respectively denotes the liquid and vapor phases, pÂ», is the reference vapor density, G(t) is the volumetric internal heat generation rate. The thermal expansion coefficient, /?, is introduced into the umomentum equation to account for buoyancy, according to the Boussinesq approximation AZkBACTÂ©] (6.5) Substituting Eq.(6.5) into the last term in Eq.(6.2) and upon expansion yields Pig = gP,AP,
PAGE 102
B7 dpâ€¢ dp ^T+VÂ«(Afr) = ^+PrÂ„V.(SVF) (6.9) 85,6 , _ 1 d(pf) 1 ^ + V.(p^^) = ^V 2 e + ^G(T)sin(^)^V.(^r/) (6.10) where V is nondimensional velocity vector, and g is the gravity factor representing the fraction of earth gravity the system will be subjected to. The Stefan number is also defined as St = c v T / Ae where Ae is the latent heat. The Raleigh number, based on oneg and vapor properties, is defined Ra v ^MH (611) The nondimensional internal heat generation rate is G(r)= p2(6.12) where tr=L 2 /a v , is the thermal diffusion time scale. The Prandtl number for vapor phase Pr v =^ (6.13) The liquid Raleigh number and Prandtl number are also implicitly included in the governing equations. In obtaining the above equations, the following dimensionless variables are used e = , pBjT Â°' ?,.Â»' Pv.oK.o/^) 2
PAGE 103
88 x= t â€¢= i' T =rwk*)> (6H) U = Â— , V = , V=Ut+Vt (a v0 /L) Ko/i) where u r =a v ,{/L is the thermal diffusion velocity scale in the vapor phase, e^and e r axe the axial and radial unit vectors, respectively. Other velocity scales, such as buoyancy velocity scale, etc., can also be used to yield different forms of the governing equations. It is chosen that the initial density ratio of liquid to vapor, p t /fK, , is 100, and the initial conductivity ratio, kj/kVi0 , is 10, the vapor Raleigh number RaÂ„ is 1.2xl0 6 , vapor Prandtl number Pr, is 0.53, Stefan number St is 0.67. 6.3 Numerical Treatment Fully implicit time marching is used for the temporal difference. The hybrid scheme is used for the spatial difference, where the first order upwind scheme for the convection terms and the second order central difference scheme for diffusion terms are used accordingly. The computation is unstable and difficult to converge if only the central difference scheme is used for the convection terms, during the computing of the transient process. One reason is that there are interactions between the natural convection cells in the liquid phase and vapor phase, and the interaction is particularly strong near the liquidvapor interface, as will be shown later. These interactions make the computation to be unstable and easily divergent. Another reason is that there is a computational grid limitation of cell
PAGE 104
89 Reynolds number, Re M , or cell Peclet number, Pe^. When the Rayleigh number or Grashof number is considerably large, a very thin grid has to be adopted, which is time consuming and expensive. For example, the natural convection velocity scale is u, = (GrJ v/L, Gr v is the Grashof number based on vapor properties. The cell Reynolds number is defined as Re M = u,Ax/v = (Grj" 2 AX. It is required that Re M be not greater than 2 for computational stability. To satisfy this, if Gr, = 10 6 , AX must be smaller than 1/500. Such a dense mesh is not realistic for the computation capacity of ordinary workstations. By adopting underrelaxation and more innerloop iterations for convection equations, etc., higher Re^ , i.e, 10 25, can be tolerated for some simple problems of natural convection in a single phase system, such as those of cases in chapter 5. Associated with the phase change and the interactions of convection cells between two phases, the computational stability becomes more difficult to be guaranteed. Very dense meshes, for example, 41x101, 81x201, were tested for the cases of g = 10"', Rav = 10 6 , where the Re^ is only 1020, the central difference scheme failed. A multischeme procedure is designed for the computation The first order upwind scheme has to be used for the convection terms to guarantee the computational stability and convergence during computing the transient process. When the system is close to the steady state, it is switched to the second order central difference scheme to improve the accuracy of results. The price paid is the lower accuracy for the transient process. In the following computations, the mesh with 46x101 uniform grids is used. The time step size is 10" About 200 ~ 700 iterations are needed for convergence at each time step. With such mesh and 2000 time steps, the typical CPU time in the Workstation of SunSparc 20 is
PAGE 105
90 about 10 days. Coarser mesh, 51x26, is also used and the steady state's results with different grid sizes and different schemes are compared. 6.4 Results of MicroGravity HO" 3 g ) The evaporation and condensation processes under microgravity and normal gravity conditions, that is g = 10 3 and 1, are modeled. The dependence of the phase change performance on the gravity level is shown clearly. In the following results, the evolution of the liquidvapor interface topology, the temperature distribution, the flow pattern of natural convection in liquid and vapor phases, are demonstrated. Initially, the lower half of the container is filled with saturated liquid, and the upper part is occupied by vapor. The aspect ratio is 2.5. Because of the axial symmetry, only the right half of the container with the liquid and vapor is modeled. The internal heat generation provides heat both in the liquid and vapor phases according to the following formulae G(t)=\2t 01.5 (6 15) At t = 0, the side wall also starts to be cooled with an outward constant heat flux. The concurrent internal heating and external cooling lead to the simultaneous bulk evaporation and condensation. As time increases, more and more heat is generated until X = 1.5 when the heat generation keeps constant. The results of the evaporation/condensation process is shown in Figures 6.2 through 6.9. Using the phase identifier which is the local vapor phase fraction, the
PAGE 106
91 distribution of the liquid and vapor phases and the position of the interface are determined. The topology of the liquidvapor interface and the temperature distribution at various time instants are displayed in Figures 6.2 and 6.3, respectively. The flow fields in the vapor and liquid phases are presented in Figures 6.8 and 6.9. The internal heat generation evaporates the liquid, while the side wall cooling tends to condense the vapor. Globally, the evaporation of liquid and condensation of vapor take place simultaneously. The local phasechange process, whether evaporation or condensation, is dominated by either the local internal heating or the wall cooling effect, since the internal heat generation is spatially dependent. The evaporation is dominant near the central axial position where the heat generation intensity is higher and evaporates more liquid. The condensation caused by side wall cooling is dominant near the side wall. The simultaneous evaporation and condensation cause the motion and deformation of the interface Because the heat generation is time dependent, the local evaporation, condensation, and motion and deformation of interface continue to evolve in time. As the process evolves, a certain fraction of liquid is evaporated into vapor due to the internal heat generation which increases at the early time period, as indicated by Eq.(6.15). Meanwhile, some amount of vapor near the side wall is condensed into liquid due to the cooling of the side wall. At an earlier time, a short piece of liquid film is formed along the side wall owing to condensation, and the liquidvapor interface becomes concave because of evaporation. As time increases, the interface continues moving and deforming The liquid film grows longer and thicker, and the position of the interface at the centerline of the container descends lower and lower. When r = 1.0, the liquid film
PAGE 107
92 grows to the top end and covers the entire side wall At x= 1.5, the interface at the low central side descends to the lowest point. The liquid phase, which initially fills the lower part of the container, now forms a long film covering and residing along the entire side wall. After r = 1.5, the internal heat generation rate is a constant, as indicated by Eq. (6.15), and the cooling heat flux at the side wall is kept constant. The system reaches a steady state and the interface no longer moves. The topology of the interface strongly depends on the temperature distribution and the internal heat generation rate distribution which becomes time independent and remains space dependent. The shape of the interface does not follow the heat generation distribution which is a chopped sine function, as it does in the case of zero gravity. The evolution of the interface topology is different from that of the zero gravity case where no buoyancy convection is considered. Under the microgravity conditions, and due to the weak gravity effect, the buoyancy induced convection is not strong but it does change the overall characteristics of the phasechange process. During the transient process, the temperature distribution varies with time. The temperature contours for the entire system gradually become denser and denser, since internal heat generation is increased. When the process approaches steady state, at r2.0, the axial temperature distribution does not follow exactly the form of a chopped sine function, as it does in the 0g case. The maximum temperature position is shifted upward due to the buoyancy convection. The heat conduction is still quite significant since the natural convection is relatively weak in this microgravity environment. There is a temperature gradient jump at the interface due to the release or absorption of latent heat
PAGE 108
93 associated with phase change and the conductivity difference between vapor and liquid. The discontinuity of the temperature gradient can be observed more clearly in the three dimensional temperature distribution displayed in Figure 6.4 There is an obvious fold on the temperature surf where the liquidvapor interface is precisely located. The temperature gradient in the liquid region is greater than that in the vapor region, since the thermal power density is proportional to the mass density and more heat is generated in the liquid phase. The axial temperature profile on the side wall is much flatter than that in the vapor region. Liquid is a much better conductor than vapor, and thus the liquid film covering the side wall is more isothermalized. This isothermalization effect is expected to satisfy the design constraint of a multiphase fuel cell for potential application in space. The formation of the liquid film and its full coverage on the side wall can prevent the side wall from being overheated, which is also the expectation for the application of such multiphase fuel cell. The temperature distribution at each time instant is different from that in the case of zero gravity, since the convection induced by gravity, though very weak in this case, affects the temperature field. The convection patterns are discussed later. The transient characteristic of this system can also be seen from the variation of the saturation temperature with the time shown in Figure 6.5. The saturation temperature decreases a little bit in the early time period when the internal heating is very weak, the side wall cooling effect is dominant. As the internal heating rate increases, the evaporation becomes dominant and the saturation temperature rises. After x= 1.5, as the system tends to reach steady state when the evaporation and condensation rates are balanced and the saturation temperature remains constant.
PAGE 109
94 Different grid sizes are used to examine the spatial accuracy of the results, as shown in Figures 6.5 through 6.7. As shown in Figure 6.7, the results of variation of saturation temperature with two grid sizes, 101x41 and 51x21, are quite close. The deviation is less than 1.0%. The results with coarser grid size are less smooth. As has been mentioned earlier, a multischeme procedure has to be adopted to guarantee the convergence and gain accuracy. The first order upwind scheme for convection terms is used during the transient period, and it is switched to the second order central difference scheme when the steady state is nearly reached. At T1.95 when the process has reached steady state, the second order central scheme is invoked and used until r = 2.0. The axial temperature profile in the position of half radius at r = 2.0 is chosen to examine the deviations between the solutions with different schemes and grid sizes, as shown in Figure 6.6. The solutions with the same fine grid size, 101x41, but different schemes, are almost identical. It indicates that the solution with finer grid size is quite accurate. Coarser grid size, 51x21, yields less smooth and less accurate results, and the deviation is about 1%. Figure 6.7 shows the evolution of the liquidvapor interface with a coarse mesh, 51x21. Compared to that in Figure 6.2, the resultant topology of the liquidvapor interface with a coarse mesh is less smooth. Figure 6.8 displays the flow patterns of the natural convection at r = 2.0 in the vapor and liquid phases, respectively. In order to be shown more clearly, not all data, only 25x20 of 101x41, are used to plot the velocity vector. Vapor and liquid are two phases with largely different properties, and their flow and heat transfer characteristics are
PAGE 110
sufficiently different that they need to be studied separately. The values of liquid and vapor velocities are quite different and are not even on the same order of magnitude. The vapor phase basically moves upward due to evaporation and natural convection. The mass transfer across the liquidvapor interface due to condensation and evaporation can be observed. In the region near the horizontal portion of the interface, vapor moves upward faster because of a higher evaporation rate. Its axial velocity decreases with increasing axial distance. In the region near the almost vertical portion of the interface, vapor moves toward the interface where condensation is dominant. The global body force, gp,,o in Eq.(6.6) or gRg,Pi; p ffi ( ) in Eq.(6.8), Pv*0 always drives the mass moving downward. The liquid phase of larger mass density is subjected to a greater body force than the vapor phase. The vapor phase will naturally occupy the space well above the liquid phase in the case of no bulk condensation. In this case when the condensation occurs along the side wall due to its cooling, some amount of liquid is accumulated at the side wall. The continuing condensation and the accumulation of liquid lead to the formation of a thin and long liquid film along the side wall. The liquid in the film is drawn by the body force downward and back to the liquid pool. A reflux cycle about the phase transformation, natural convection, and mass transfer is formed, which is schematically demonstrated in Figure 6.1. Since the liquid film is thin, the temperature difference across it is not large enough to induce natural convection. Thus liquid basically moves downward in the film. The temperature difference across the liquid pool, however, is large enough to induce natural convection. A clockwise rotating convection cell is formed in the liquid pool.
PAGE 111
96 vapor natural convection evaporation condensation liquid film liquid pool reflux Figure 6.1 Schematic diagram of mass reflux The axial velocity profiles at different axial positions, as shown in Figure 6.9, can more clearly elucidate the flow field. The position of X = 0.1 is completely located in the liquid pool, velocity is distributed smoothly. Liquid moves upward in the central region and downward near the side wall, forming a typical natural convection cell due to the buoyancy effect, which is also shown in Figure 6.8. At positions of X = 2.0 or higher, where both vapor and liquid films are present, the velocity profiles are no longer smooth. There is a sudden change in velocity gradient at the liquidvapor interface due to the property changes. The shear force balance at the interface is expressed as follows: 3U 'at fi, _d(J_ (6.16) Since liquid viscosity is larger than vapor viscosity, the velocity gradient at the vapor side is greater than that at the liquid side.
PAGE 112
97 At these axial positions, vapor velocity is basically positive and liquid velocity is negative, except for the region near the liquidvapor interface. A countercurrent twophase flow is formed. Vapor velocity decreases at higher axial positions, indicating that the strength of natural convection decreases with increasing X. The liquid velocity at different axial positions does not change as much as vapor velocity. The convection flow patterns and the interactions between the motion of liquid and vapor can be observed. At X = 0.2, liquid moves upward near the interface and downward near the side wall. On the one hand, the vapor velocity is quite high and drives the liquid upward through shear interaction at the interface. On the other hand, the liquid film is quite thick, and the temperature difference across it is large enough to induce natural convection in the liquid phase At X= 0.3, both liquid and vapor velocities are close to zero at the interface. At the positions of X = 4 or higher, vapor convection becomes weaker, and the liquid film becomes thinner. The temperature difference across the liquid film is not large enough to induce natural convection The downward motion of liquid due to the larger body force becomes more dominant. Near the interface, vapor moves downward due to its interaction with liquid film. The temperature gradient in the wider vapor region is large enough to cause natural convection. This becomes more significant at higher axial positions. Compared to the results of the zero gravity case, the process of the liquid film formation is slower and less liquid is accumulated in the film for the microgravity case. It is due to the gravitational force and the reflux effect. With the action of the gravity force, the liquid can not move upward as freely as in the 0g case.
PAGE 113
98 vapor liquid J O' X = 0.2 0:4 0.6 r = 1.0 1 2 1.4 2.0 Figure 6.2 The evolution of liquidvapor interface ( 10" 3 g )
PAGE 114
99 ^ t = 0.2 0.4 0.6 0.8 r = 1.0 1.2 1.4 2.0 Figure 6.3 The evolution of isotherm (10" 3 g)
PAGE 115
Figure 6.4 Temperature surf at r = 2.0 ( 10" 3 g ) 1 2 8. Figure 6.5 Variation of saturation temperature ( 10" 3 g ) ( solid line: 101x41 grid; dashed line: 51x21 grid )
PAGE 116
B u 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Axial position ( X ) Figure 6.6 Solutions with different grid sizes and schemes (Axial temperature profile in half radius at T 2.0. solid line: 1st order upwind scheme with 101x41 grid; dots: 2nd order central differences scheme with 101x41 grid; dashed line: 1st order upwind scheme with 51x21 grid) 0.2 0.6 1.0 2.0 Figure 6.7 The evolution of liquidvapor interface with coarse grid ( 51x21 grid )
PAGE 117
(a) vapor phase ( U max = 5.8 ) (b) liquid phase ( U max = 0.85 ) MM Hi. tUtttH, 'Ultttii, tWtttn unttti \UUttt, ftlltllft. Figure 6.8 Flow patterns in vapor and liquid phases at x = 2.0 ( 10" 3 g )
PAGE 118
103 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Figure 6.9 Velocity U profiles at different axial positions ( x = 2.0, 10" 3 g)
PAGE 119
104 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 R Figure 6.9 Velocity U profiles at different axial positions ( 10" 3 g , continued )
PAGE 120
105 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 R 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.9 0.8 0.7 r~ 0.6 o 0.5 X 0.4 Â— ' 01 S 0.2 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Figure 6.9 Velocity U profiles at different axial positions (10' 3 g , continued )
PAGE 121
106 6.5 Results of Normal Gravity ( 1g ) All the boundary and initial conditions are maintained the same as those for the case of 10" 3 g. The gravity factor g is 1 in Eq (6.8). The results are shown in Figure 6. 10 through Figure 6.15. The topology of the liquidvapor interface and the temperature distribution at various time instants are displayed in Figures 6.10 and 6.11, respectively. The flow fields in the vapor and liquid phases are presented in Figures 6.14 and 6.15. The global evolution of the phase change process is similar to the 10" 3 g case. As time lapses, a piece of thin and long liquid film is formed and more and more liquid is accumulated along the side wall due to the cooling effect. After r = 1.5, the system gradually tends reach a steady state. Significant differences appear in this normal gravity case, compared to the 10~ 3 g case. The liquid film grows more slowly than it does with 10~ 3 g. Even at about r = 2.0, the liquid film does not move up to the top end and cover the entire surface of the side wall. The liquid film is also much thinner. Since the gravity effect is greater in the present case, the liquid film more readily falls down. More liquid is drawn from the film back to the liquid pool and less liquid is accumulated in the film. Consequently the position of the horizontal surface of the liquid pool is higher than that in the 10" 3 g case. Furthermore, the temperature distribution is quite different from those at 0g and 10" g. The axial temperature gradient across the liquid phase is much larger than that across the vapor phase, owing to the much intensive internal heat generation rate in the
PAGE 122
107 liquid phase. The radial temperature gradient is much smaller than that of 10" 3 g, indicating that more heat is removed outward owing to stronger convection. The isotherms are nearly horizontal in the region far from side wall. The natural convection is very strong and dominant, which is quite similar to the results of benchmark problems shown in Figure 5.2. In particular, the temperature gradient across the liquid film is larger than that across the vapor phase, since the liquid film is thinner. Such a thin liquid film with a large temperature gradient is a very good conductor for heat transfer. The internally generated heat in the vapor phase, which is less intensive than in the liquid phase, is barely accumulated and easier to be removed. Therefore the axial temperature gradient across the vapor region is smaller. Figure 6.14 shows the flow fields in the vapor and liquid phases. Compared to Figure 6.8, the flow patterns are different. The convection intensity is much larger in 1g than that in 10' 3 g, due to the larger gravity effect. The natural convection cells, in which some amount of liquid moves upward and downward, are still formed in the liquid pool. The reflux of some amount of liquid back to the liquid pool can be seen more clearly. Axial velocity profiles at several axial positions are plotted in Figure 6.15. They are a little similar to those of 10" 3 g case. The velocity distribution in the liquid pool at a position of X = 0.3 is quite smooth. Upward motion and downward motion of liquid form a typical natural convection cell. At higher axial positions where both vapor and liquid films are present, there is a sudden change of velocity gradient at the liquidvapor interface due to the property change. The buoyancy effect and the interaction between the motion of vapor and liquid can also be observed as they are in the 10" 3 g case. The velocity values are
PAGE 123
108 much greater than those of 10" 3 g case owing to the much larger body force and stronger convection effects. The differences and the dependence of the phase change performance on the different gravity levels can be observed more clearly in Figure 6.16 where the topology of the liquidvapor interface, the temperature distribution at the final steady state under three cases, 0g, 10" 3 g and 1g, are plotted together. The gravity influences on the thermal performance of the bulk liquidvapor phase change process in the multiphase nuclear fuel cell are displayed explicitly. Under all three gravity conditions, the liquid film is formed and covers the entire side wall, and the liquid film covering the side wall is more isothermalized at the wall surface, which can prevent the side wall from being overheated. At 0g, the thermal performance is controlled only by the heat transfer, such as internal heat generation and side wall cooling. At nonzero gravity conditions, its thermal performance is controlled by both heat transfer and gravity effects. As the gravity increases, the liquid film is thinner, the temperature gradient is larger across the liquid film and smaller across the vapor phase. The model developed in this work can be used to provide very useful, and at least qualitative, information for the design and application of multiphase nuclear fuel elements to be operated under different gravity conditions, particularly when the experiments at microgravity environment are expensive and not easy to be conducted. The current model can also be applied to some twophase flow problems with or without phase change, for example, film condensation, the liquidvapor or liquidgas flow problems in heat exchangers and material processing equipment.
PAGE 124
109 % = 0.2 0.4 0.6 0.8 r = 1.0 1.2 1.4 2.0 Figure 6.10 The evolution of liquidvapor interface ( 1g )
PAGE 125
110 3 02 0.4 0.6 0.8 1.0 1.2 1.4 2.0 Figure 6 1 1 The evolution of isotherm ( 1 g )
PAGE 126
Ill ' *=^ as Figure 6. 12 Temperature surf at % = 2.0 ( 1 g ) Figure 6.13 Variation of saturation temperature ( 1 g )
PAGE 127
112 (a) vapor phase ( U mca = 210 ) (b) liquid phase (  = 65) I t f \Uttrr?/fs*. lUlllll///////^ WWlfllf!///^ IlltilllllllUs munm\\\\, immmnu, 'tttfttfffff///A mmmttWt. f 1: *i >" H Figure 6. 14 Flow patterns in vapor and liquid phases at x = 2.0 ( 1 g )
PAGE 128
113 200 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Figure 6. 1 5 Velocity U profiles at different axial positions ( x = 2.0, 1 g )
PAGE 129
114 (a) The topology of the liquidvapor interface at z= 2.0 0g lO" 3 g (b) The temperature distribution at x= 2.0 1g Figure 6. 16 Comparison of interface and temperature distribution at different g levels
PAGE 130
CHAPTER 7 CONCLUSIONS AND RECOMMENDATIONS Based on theoretical modeling and computational results of the present work, the following conclusions can be drawn: 1. A comprehensive model is successfully developed to simulate the dynamic process of bulk liquidvapor phase change associated with internal heat generation and natural convection under various gravity conditions. 2. The internal energy formulation for the constant volume phase change process is proposed and successfully implemented. The heat transfer mechanism for this type of phase change problems, including the convection effect induced by the density change, and velocity and time scales of the interface movement, is analyzed and discussed. The convection effect induced by the density change is analyzed and found negligible, and the dominating mechanism for the heat transfer is conduction under certain heating/cooling conditions. 3. The computational performance of the internal energy formulation and enthalpy formulation are evaluated. Both formulations yield identical results, while the internal energy formulation has the advantage of having more compact form. The performance of two update methods for the vapor phase fraction, the Ebased and Tbased methods, are investigated The Ebased method is superior to the Tbased method in terms of 115
PAGE 131
116 computational efficiency and formulation simplicity for evaporation and condensation problems. The Tbased method is intractable for twodimensional bulk evaporation and condensation problems, while the Ebased method always works quite well. 4. A one dimensional conductiondriven bulk evaporation and condensation problem is simulated. The evolution of the bulk liquidvapor phase change, the variations of the temperature distribution and the saturation temperature with time for a variety of boundary conditions are investigated. 5. The investigation of bulk evaporation and condensation processes associated with internal heating under zero gravity, which is also conduction driven, is performed. It displays the motion and deformation of the liquidvapor interface, the variation of saturation temperature, the dependence of temperature distribution and interface topology on the internal heat generation rate. 6. Some natural convection problems associated with external or internal heating in a singlephase system are modeled. The solutions agree quite well with benchmark results and experimental results. The natural convection caused by external or internal heating, in the twophase system without phase change, is studied. The convection flow patterns in both phases are obtained. It is revealed that there exist interactions between the convection cells in liquid and vapor phases. Such interactions complicate the flow fields and the modeling procedure. 7. Based on the above investigation, a complete model is established to simulate the bulk liquidvapor phase change process associated with internal heat generation, gravity effect and natural convection. A complete set of equations governing the
PAGE 132
117 conservation of mass, momentum and energy for both phases is solved. The thermal performance of a multiphase nuclear fuel element under microgravity and normal gravity conditions is investigated. The numerical simulation elucidates the evolution of the bulk liquidvapor phase change process. This includes the evolution of the liquidvapor interface, the formation and development of the liquid film covering the side wall surface, the temperature distribution and the convection flow field. 8. Numerical modeling results show the strong dependence of the thermal performance of such multiphase nuclear fuel cell on the gravity conditions. The differences of the phase change performance of the multiphase fuel element at various gravity levels are quite significant. Under all three gravity conditions, the liquid film is formed and covers the side wall. The liquid film covering the side wall is more isothermalized at the wall surface, which can prevent the side wall from being overheated. At 0g, the performance is controlled only by the heat transfer, such as internal heat generation and side wall cooling. At nonzero gravity conditions, its performance is controlled by gravity effect, as well as by heat transfer. As gravity increases, the liquid film is thinner, the temperature gradient is larger across the liquid film and smaller across the vapor phase. The followings are recommended for further research: 1. The model developed in this work can be used to predict and provide very useful, and at least qualitative, information for the potential design and application of multiphase nuclear fuel elements to be operated under different gravity conditions,
PAGE 133
118 particularly when the experiments at microgravity environment are expensive and not easy to be conducted. 2. Further effort will be very worthwhile and strongly recommended to extend the current model for liquidvapor phase change problems to a unified threephase model. All phase change processes between these three phases, such as solidliquid, liquidvapor, solidvapor, or solidliquidvapor, could be simulated. 3 . The current model can also be applied to some twophase flow problems with or without phase change, for example, the film condensation, the liquidvapor or liquidgas flow problems in heat exchangers and material processing equipment. 4. Most of the thermophysical properties in the present work have been assumed constant. In the future, variable properties can be implemented to get more accurate results with more computational expense. 5. Relative experiments, using the materials with accurate properties, should be conducted to verify and improve the model.
PAGE 134
REFERENCES Acharya, S. & Goldstein, R. J. (1985), Natural convection in an externally heated vertical or inclined square box containing internal energy sources, J. Heat Transfer . 107 . 855866 Anghaie, S. (1992), Thermophysical properties of UF4 at high temperatures (Internal Report), Innovative Nuclear Space Power and Propulsion Institute, Univ. of Florida Anghaie, S. & Ding, Z. (1995), Thermal performance of a high temperature phase change thermionic fuel cell , paper presented at 1995 National heat transfer Conf., Portland, Oregon Bejan, A (1984), Convection Heat Transfer . lohn Wiley & Sons, New York Bergman, T. L. & Ramadhyani, S. (1981), Combined buoyancyand thermocapillarydriven convection in open square cavities, Numer. Heat Transfer . 9, 44 1 45 1 Brice, J. C (1986), Crystal Growth Processes . Blackie, London Brown, R. A. (1988), Theory of transport processes in single crystal growth from the melt, A.I.Ch.E. J . 34, 881911 Buckle, V. & Peric, M. (1992), Numerical simulation of buoyant and thermocapillary convection in a square cavity, Numerical Heat Transfer. Part A . 21 . 121141 . Carey, V. P. (1992), LiquidVapor PhaseChange Phenomena . Hemisphere Pub. Corp., New York Carpenter, B. M. & Homsy, G. M. (1989), Combined buoyantthermocapillary flow in a cavity . J. Fluid Mech .. 207 . 121132 119
PAGE 135
120 yCossmann, R., Odenthal, H.P. & Renz, U. (1982), Heat and mass transfer during partial condensation in a turbulent pipe flow, U. Grigul (Eds), Proc. 7th Int. Conf. of Heat Transfer , Hemisphere Publishing Corp., Washington, DC. Crank, J. (1984), Free and Moving Boundary Problems. Clarendon Press, Oxford de Vahl Davis, G. (1983), Natural convection of air in a square cavity: a bench mark numerical solution, Intl. J. for Numerical Meth. in Fluids . 3, 249264 Ding, Z. & Anghaie, S. (1994a), Numerical investigation of the twophase equilibrium state in a cylindrical nuclear fuel cell under zero gravity condition , paper presented at the 6th AJAA/ASME Joint Thermophysics and Heat Transfer Conf, Colorado Springs, CO Ding, Z & Anghaie, S. (1994b), Modeling of R12 bulk evaporation in an encapsulated container , paper presented at the ASME Intl. Mechanical Engineering Congress & Exhibition, Chicago, IL Ding, Z. & Anghaie, S , Numerical modeling of bulk evaporation and condensation with constant volume, Intl. J. Numerical Methods in Engineering (in press) Ding, Z. & Anghaie, S. (1995b), Modeling of bulk evaporation , paper presented at the 30th AIAA Thermophysics Conf, San Diego, CA Ding, Z. & Anghaie, S., Analysis and formulation of bulk evaporation and condensation problems, Intl. J. Heat & Mass Transfer . (Submitted) Eckert, E R. G, Goldstein, R. I, et al. (1990), Heat transferA review of 1989 literature, Intl. J. Heat and Mass Transfer . 33, No. 12, 23492437 Eckert, E. R. G, Goldstein, R. J., et al. (1991), Heat transferA review of 1990 literature, Intl. J. Heat and Mass Transfer . 34. No. 12. 29313010 Eckert, E. R. G, Goldstein, R. J , et al. (1992), Heat transferA review of 1991 literature, Intl. J. Heat and Mass Transfer . 35, No. 12 Flik, M. I. & Tien, C. L.(1989), An Approximate analysis for general film condensation transients, J. Heat Transfer . 111. 511517.
PAGE 136
121 Fusegi, T., Hyun, J. M. & Kuwahara, K. (1992), Natural convection in a differentially heated square cavity with internal heat generation, Numer Heat Transfer . 2]_, 215229 Gau, C. & Viskanta, R (1986), Melting and solidification of a pure metal on a vertical wall, J. Heat Transfer . 108 . 174181 Gaultier, M, Lezaun, M. & Vadillo, F.(1993), A problem of heat and mass transfer: proof of the existence condition by a finite difference method, Intl. J. for Numerical Methods in Fluids . 16, 87104 Gebhart, B., Jaluria, Y., Mahajan, R. L. & Sammakia, B.(1988), BuoyancyInduced Flows and Transport . Hemisphere, Washington D. C. Gerner, F. M. & Tien, C. L.(1989), Axisymmetric Interfacial Condensation Model, J. Heat Transfer . 111 . 503510 Hatsopoulos, G. N. & Gyftopoulos, E. P. (1973) Thermionic Energy Conversion Vol.1: Process and Devices . MIT Press, Cambridge, MA Harley, C & Faghri, A. (1994), Complete transient twophase closed thermosyphons including the falling condensate film, J. Heat Transfer . 116 . 418426 Hortmann, M., Peric, M. & Scheuerer, G. (1990), Finite volume multigrid prediction of laminar natural convection: benchmark solutions, Intl. J. Numerical Meth. Fluids . 11. 189207 Hsieh, C. K., Akbari, M. & Li, H. (1993), Solution of inverse Stefan problems by a sourceandsink method, Intl. J Numerical Methods Heat and fluid Flow . 2, No.5, 391406 Hsieh, C. K. & Choi, C.Y. (1992), Solution of oneand twophase melting and solidification problems imposed with constant or timevariant temperature and flux boundary conditions, J. Heat Transfer . 114 . 524528 Hung, C. I., Shyy, W. & Ouyang, H. (1995), Transient natural convection and conjugate heat transfer in a crystal growth device, Int. J. Heat Mass Transferm. 38, 701712
PAGE 137
122 Jang, J. H., Faghri, A., Chang, W. S. & Mahefkey, E. T.(1990), Mathematical modeling and analysis of heat pipe startup from the frozen state, J. Heat Transfer . 112 , 586594 Jones, W. P. & Renz, U. (1974), Condensation from a turbulent stream onto a vertical surface, Int. J. Heat Mass Transfer , 17, 10191028 Kimura, S. & Bejan, A. (1983), The "heatline" visulization of convective heat transfer, J. Heat Transfer , 105 . 916919 Lame, G. & Clapeyron, (1831) Ann. Chem. Phvs .. 47, 250256 Lee, H., Lewis, B. R & Klein, A. C. (1993) System modeling for the advanced thermionic initiative single cell thermionic space nuclear reactor, Proc. 10th Symposium on Space Nuclear Power and Propulsion . Conf930103, M. S. ElGenk & M. D. Hoover, eds., American Institute of Physics, New York. Lee, J. & Goldstein, R. J. (1988), An experimental study on natural convection heat transfer in an inclined square enclosure containing internal heat sources, J Heat Transfer . 110 . 345349 Lewis, B. R., Powlowski, R. A., Greek, K. J. & Klein, A. C. (1991) Advanced thermionic reactor system design code, Proc. 8th Symposium on Space Nuclear Power Systems . Conf9101 16, M. S. EIGenk & M. D. Hoover (Eds), American Institute of Physics, New York Li, H. (1993), SourceandSink method of Solution of Twoand ThreeDimensional Stefan Problems, Ph.D. Dissertation. Mechanical Engineering Dept., Univ. of Florida McVey, J. B. & Rasor, N. S. (1992) The TECMDL thermionic converter computer model, Proc. 27th Intersocietv Energy Conversion Engineering Conference . Society of Automotive Engineers, San Diego, CA: 2125 Nadarajah, A. & Narayanan, R. (1990), Comparison between morphological and RayleighMarangoni instabilities, Meinkoen, D. & Haken, H. (Eds), Dissipative
PAGE 138
123 Structures in Transport Processes and Combustion . 215228, SpringerVerlag, New York Nusselt, W. Z. (1916), Die oberflachenkondensation des wasserdampfes, Z. Ver. Distch Ing., 60, 541569 Ostrach, S. (1982), Lowgravity fluid flows, Ann. Rev. Fluid Mech .. 14. 313345 Ostrach, S. (1983), Fluid mechanics in crystal growthThe 1982 freeman scholar lecture, J. Fluids Engineering, 105 . 520 Ostrach, S. (1988), Natural convection in enclosures, J. Heat Transfer . 110. 1 1751 190 Ozisik, M. N. (1980), Heat Conduction . John Wiley, New York Patankar, S. V. (1980), Numerical Heat Transfer and Fluid Flow . Hemisphere Publication Corp., New York Poirier, O. & Salcudean, M (1988), On numerical methods used in mathematical modeling of phase change in liquid metals, J. Heat Transfer . 110 . 562570 Prakash, C, Samonds, M. & Singhal, A. K. (1987), A fixed grid numerical modeling methodology for phase change problems involving a moving heat source, Int. J. Heat Mass Transfer. 30, 269094 Ramaswamy, B. & Jue, T C. (1992), Analysis of thermocapillary and buoyacyaffected cavity flow using FEM, Numer. Heat Transfer . 22, 379399 Rasor, N. S. (1991) Thermionic energy conversion plasmas, invited review paper, IEEE Transactions on Plasma Science . 19, No. 6 Reynolds, W. C. (1979), Thermodynamic Properties in SI . Dept. of Mechanical Engineering, Stanford University Rohsenow, W. M.(1973), "Film Condensation" In Handbook of Heat Transfer (Edited by W. M. Rohsenow & J. P. Hartnett), Set. 12A. McGrawHill, New York Samarskii, A. A , Vabishchevich, P. N, Iliev, O. P. & Churbanov, A. G. (1993), Numerical simulation of convection/diffusion phase change problemsA Review, Int. J. Heat Mass Transfer . 36, 40954106
PAGE 139
124 Shamsundar, N. & Sparrow, E. M. (1975), Analysis of multidimensional conduction phase change via the enthalpy model, J. Heat Transfer , 97, 333340 Shyy, W. (1994a), Computational Modeling for Fluid Flow and Interfacial Transport . Elsevier, Amsterdam, The Netherlands Shyy, W. (1994b), Modeling of transient twophase heat transfer for spacecraft thermal management, Microgravitv Sci. Tech .. 7, 219227 Shyy, W. & Chen, M. (1990), Steadystate natural convection with phase change Int J Heat Mass Transfer . 33, 25452563 Shyy, W. & Chen, M. (1991), Interaction of thermocapillary and natural convection flow during solidification: normal and reduced gravity conditions, J. Crystal Growth . 108 . 247261 Shyy, WÂ„ Pang, Y., Hunter, G. B., Wei, D. Y. & Chen, M.H. (1992a), Modeling of turbulent transport and solidification during continuous ingot casting, Int. J. Heat Mass Transfer . 35. 12291245 Shyy, W. & Rao, M. (1992b), Convection treatment for high Rayleigh number, laminar, natural convection calculation, Numerical Heat Transfer. Part B . 22, 367374 Shyy, W , Udaykumar, H. S. & Liang, S.J. (1993a), An interface tracking method applied to morphological evolution during phase change, Int. J. Heat Mass Transfer . 36. 18331844 Shyy, W., Gringrich, W. K., Krotiuk, W. J & Fredley, J. E. (1993b), Transient twophase heat transfer in a capillarypumpedloop reservoir for spacecraft thermal management, Proc. 10th Svmp.m on Space Nuclear Power and Propulsion . Conf930103, American Institute of Physics Shyy, W & Rao, M. M. (1993c), Simulation of transient natural convection around an enclosed vertical channel, J. Heat Transfer . 115 , 946954 Shyy, W. & Rao, M. M. (1994), Enthalpy based formulations for phasechange problems with application to gjitters, Microgravitv Sci. Technology , 7, 4149
PAGE 140
125 Shyy, W., Liang, S.J. & Wei, D. Y. (1994), Effect of dynamic perturbation and contact condition on edgedefined fiber growth characteristics, Int. J. Heat Mass Transfer. 37, 977987 Shyy, W., Udaykumar, H. S., Rao, M. M. & Smith, R. W. (1996), Computational Fluid Dynamics with Movina Boundaries . Talyor & Francis, Washington, DC. Sparrow, E. M. & Eckter, E. R. G. (1961), Effects of superheated vapor and noncondensable gases on laminar film condensation, A.I.Ch.E. J. . 7, No. 3 Sparrow, E. M. & Lin, S. H. (1964), Condensation heat transfer in the presence of noncondensable gas, J. Heat Transfer . 86, 430436 Sparrow, E. M., Minkowycz, W. J. & Saddy, M. (1967), Forced convection condensation in the presence of noncondensables and interfacial resistance, Intl . J. Heat Mass Transfer . 10. 18291845. Stefan, J. (1891), Ann. Phvs. Chemie . 42, 269286. Tournier, J.M. & ElGenk, M. S. (1992), "HPTAM" Heatpipe transient analysis model: an analysis of water heat pipes, Proc. 9th Svmp. Space Nuclear Power Systems Voller, V. R. & Prakash, C. (1987), A fixed grid numerical modelling methodology for convectiondiffusion mushy region phasechange problems, Int. J. Heat Mass Transfer . 30. 17091719 Voller, V R. & Swaminathan, C. R. (1991), General sourcebased method for solidification phase change, Numerical Heat Transfer. Part B . 19, 175189 Watanabe, Y. & Anghaie, S. (1993) Thermophysical properties of gas phase uranium tetrafluoride, paper presented at the AIAA 28th Thermophvsics Conf Yao, L. S. & Prusa, J. (1989), Melting and freezing, in Advances in Heat Transfer . Hartnett (Eds.), Academic Press, Inc., San Diego, J9, 2385 Yong, T. J., Thayer, K. L. & Ramalingan, M. L. (1993) Performance simulation of an advanced cylindrical thermionic fuel element with a graphite reservoir, paper presented at the AIAA 28th Thermophvsics Conf Orlando, FL
PAGE 141
BIOGRAPHICAL SKETCH Zhongtao Ding was born on October 29, 1960, as the second child of two to Di Hua and Chengliang Ding in Hangzhou, China He was enrolled in Northwestern Polytechnic University in Xian, China where he received his B.S degree in aerospace propulsion engineering in 1982. He entered the Nanjing Aeronautical Institute ( NAI, currently Nanjing University of Aeronautics and Astronautics ) in Nanjing, China. He earned his MS. degree from the aerospace propulsion engineering department in 1985. He then worked in NAI as an assistant lecturer from 1985 to 1987 and a lecturer from 1987 to 1991, teaching courses and conducting research in the areas of thermodynamics, heat transfer, combustion and fluid flow. He married Zhaohui Cheng in 1989. He attended the University of Florida to pursue Ph.D. in 1991. He joined the Ph.D. program in the Mechanical Engineering Department and worked as a research assistant at Innovative Nuclear Space Power and Propulsion Institute (INSPI). He is the member of American Society of Mechanical Engineers (ASME) and American Institute of Aeronautics and Astronautics (AIAA). He earned his Ph.D. in December 1995 from the University of Florida.
PAGE 142
I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. ner, Chair I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Samim Anghaie, CoChair Professor of Nuclear Engineering Sciences I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. A&L C.K. Hsieh Professor of Mechanical Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy n Wei Shyy Professor of Aerospace Engineering, Mechanics, and Engineering Science
PAGE 143
I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Renwei Mei Associate Professor of Aerospace Engineering, Mechanics, and Engineering Science This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy December, 1995 /* Winfred M. Phillips Dean, College of Engineering Karen A. Holbrook Dean, Graduate School
PAGE 144
LD 1780 1995 oniveRS'IT OF FLOWO* s n Â»gS554 8732

