Citation
Modeling of bulk evaporation and condensation with internal heating at various gravity levels

Material Information

Title:
Modeling of bulk evaporation and condensation with internal heating at various gravity levels
Creator:
Ding, Zhongtao, 1960-
Publication Date:
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 )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1995.
Bibliography:
Includes bibliographical references (leaves 119-125).
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 non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
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 Solid-Liquid Phase Change Problems ................................. 9
2.2 Numerical Modeling of Solid-Liquid Phase Change Problems ................... 11
2.3 Liquid-Vapor 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 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 E-Based and T-Based 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
SINGLE-PHASE AND TWO-PHASE SYSTEMS .......................... 59

5.1 External Natural Convection in a Square-Cavity Benchmark Problem ....... 60 5.2 Natural Convection in a Square-Cavity with Internal Heating .................... 67
5.3 Natural Convection in Liquid-Gas 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 Liquid-Gas 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 MICRO-GRAVITY 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 icro-Gravity (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 multi-phase 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 two-phase cell operated at 0-g .............................. 46

4.2 The control volume contacting other material ...................... ................ 50

4.3 Parameter for heat generation rate ............. ........... ................ 52

4.4 The evolution of liquid-vapor interface (0-g) .......................................... 56

4.5 The evolution of temperature field (0-g) ....................... ................... 57

4.6 Temperature surf of the system at steady state (0-g) ............ ............ 58

4.7 Variation of saturation temperature (0-g) ..................................................... 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 two-phase without internal heating ( case 1) .............. 78

5.10 Natural convection of two-phase without internal heating ( case 2 ) ............... 79

5.11 Natural convection of two-phase with internal heating ( case 1 ) ................... 82

5.12 Natural convection of two-phase with internal heating ( case 2) .................. 83

6.1 Schematic diagram of mass reflux ..................................... 96

6.2 The evolution of liquid-vapor interface ( 10-3 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 liquid-vapor interface with coarse grid .............................. 101

6.8 Flow patterns in vapor and liquid phases r = 2.0 ( 10-3 -g ) ........................ 102

6.9 Velocity U profile at different axial positions (10- -g ) ............................... 103

6.10 The evolution of liquid-vapor interface ( 1-g ) .......... ............ .. ............. 109

6.11 The evolution of isotherm s ( -g ) ......................................... ............... 110

6.12 Tem perature surf at = 2.0 ( 1-g ) ............................................................. 111

6.13 Variation of saturation temperature ( 1-g ) ................................. 111

6.14 Flow patterns in vapor and liquid phases r = 2.0 ( 1-g) ............................. 112

6.15 Velocity U profile at different axial positions ( 1-g) ............... ................. 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 Co-Chairman: 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 multi-phase nuclear fuel element at zero gravity, micro-gravity, and normal gravity is investigated. In the simulation, the numerical solutions have revealed much useful information, including the evolution of the bulk liquid-vapor phase change process, the evolution of the liquid-vapor 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 multi-phase 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 E-based and T-based methods, is investigated. Simulation of a one-dimensional conduction-driven 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 two-phase 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 non-mechanical 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 thermionic-converter 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 multi-phase 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 multi-phase 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 liquid-vapor 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 multi-phase nuclear fuel system


A cylindrical container filled with multi-phase nuclear fuel, such as UF4, is shown in Figure 1.2. The fuel experiences liquid-vapor 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 zero-gravity condition, buoyancy-driven convection is nonexistent. However, under micro-gravity or one-gravity fields, buoyancy-induced 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 multi-phase 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 solid-liquid 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 solid-liquid 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 solid-liquid 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 buoyancy-driven 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 liquid-vapor 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 multi-phase nuclear system, including the multi-phase 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 one-dimensional, externally heated or cooled and conduction-driven 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 liquid-vapor interface topology are demonstrated. Chapter 5 describes the modeling of the combined internal and external natural convection in the single-phase system and two-phase 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 bench-mark 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 buoyancy-induced 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 Solid-Liquid 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 Solid-Liquid 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 Liquid-Vapor 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 gravity-flow 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 parabolic-type 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 start-up from the frozen state by using one-dimensional 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 two-phase transient flow and heat transfer characteristics within a constant pressure reservoir of a capillary-pumped 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 buoyancy-driven 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 post-accident 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 solid-liquid 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
1-D 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 solid-liquid phase change problems, such as those involved in melting and solidification processes. There has also been an ongoing interest in the investigation of liquid-vapor 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 El-Genk (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 start-up from the frozen state by using a system of simplified one dimensional equations. Shyy (1994b) has successfully developed a computational method for predicting the two-phase transient flow and heat transfer characteristics within a constant pressure reservoir of a capillary-pumped-loop 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 liquid-vapor 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 multi-dimensional phase-change 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 phase-change 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 Clausius-Clapeyron equation, the saturation temperature as a function of the system pressure is calculated:

C
Tsat - (3.7) D-Inpsat

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 phase-change 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 liquid-vapor 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 energy-balance 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 1-pa
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 - A-In(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 T-based update method. This technique has proven to be very effective for achieving convergence. An alternative formulation to the T-based 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 H-based 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 E-based 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 E-based method compared to the T-based 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 E-based 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 E-based 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 liquid-vapor phase change occurring far away from its critical point, the latent heat is relatively large and the application of the E-based update method even for simple geometry is recommended. In the present work, the T-based and E-based methods are examined in detail from the viewpoint of accuracy and computational efficiency.

A straight-forward implementation of the T-based 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 T-based 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 under-relaxation 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 +fpn-1] < ,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 E-based method is relatively straightforward and does not need any pre-conditioning. 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 1-f (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 +(1-f)v, f 1- f Pv P1

where v is the specific volume per unit mass. It is similar to Eq.(3.30).





33





fAX (1-Axf) 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 non-uniform 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 non-uniform variation can be seen due to the nature of the discretization and the non-uniform 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 liquid-vapor 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 E-based method; the dashed line: the T-based 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 E-Based and T-Based Update Methods


The T-based update method and the E-based update method are examined for the same testing case. Figures 3.5(a) and 3.5(b) show that the T-based method and E-based method yield nearly identical results. It is seen that convergence for the E-based method requires about half the iterations which is needed for the T-based method. Thus the Ebased method is computationally more efficient than the T-based 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 T-based is better than the H-based method as reported by Shyy and Rao (1994).

The reason for the improved efficiency of the E-based method is that the phase change temperature window is used for the T-based method, which is not required for the E-based method. When the T-based 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 T-based method converges more slowly. When the E-based 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 T-based 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 E-based method can completely circumvent the need for a phase change window. Therefore, the Ebased method is superior to the T-based 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 E-based 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 E-based 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 evaporation-condensation 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 vapor-liquid 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 well-established 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 E-based and T-based methods, are investigated. The E-based method is superior to the T-based method in terms of computational efficiency and formulation simplicity for evaporation and condensation problems. The E-based 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 E-based 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 liquid-vapor phase change processes to a unified three-phase model. All phase change processes between these three phases, that is, solid-liquid, liquid-vapor, solid-vapor, or solid-liquid-vapor, 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 0-g 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 liquid-phase and condenses the vaporphase. For this system with such an orientation and operated under 0-g, there is no bubble or droplet generation on the wall surface. The temperature along the liquid-vapor interface does not change significantly, because the curvature of the interface is moderate. The liquid-vapor interface can be reasonably assumed to be isothermal, consequently the thermal capillary effect can be neglected. Hence the phase-change 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 non-uniform 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 phase-change (b) at the final steady state
process







Figure 4.1 Schematic diagram of the two-phase cell operated at O-g





47





4.1 The Governing Equations



The internal energy formulation of the energy equation for this phase-change 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 E-based update method discussed in chapter 3 and used for the iterative update of the vapor fraction. The T-based method does not work well for this twodimensional bulk evaporation and condensation problem, while the E-based always works quite well. Some additional closure relations, the same as Eqs.(3.7-8) 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

T-TN = 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), =(1-f, 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 liquid-vapor 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 phase-change 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 liquid-vapor 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 liquid-vapor 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 multi-phase and internally heated system. The present results will be used as reference for comparison with the performance and features of this system operated under micro-gravity 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 liquid-vapor interface (O-g)





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 (0-g)





58



























Figure 4.6 Temperature surf of the system at steady state (O-g)




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 (0-g)













CHAPTER 5
NATURAL CONVECTION WITH INTERNAL HEATING IN SINGLE-PHASE AND TWO-PHASE 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 non-zero 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 co-existence 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 single-phase and two-phase systems are investigated. A bench-mark problem, natural convection in a square cavity of single-phase without internal heating, is computed. Comparisons of the present solutions with accurate published bench-mark solutions are made to verify the computational accuracy. The natural convection in a square cavity of single-phase with internal heating is also computed. Qualitative comparisons with the published experimental results are also made. Then natural convection in a liquid-gas system without phase-change, 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 Square-Cavity 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 buoyancy-driven 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, buoyancy-driven flow in a square cavity with vertical side walls heated differentially, a so-called 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 two-dimensional 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 = (T-Td)/(T-T) (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 second-order central-difference 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 skew-symmetry 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, Pr-0.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 Square-Cavity 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 post-accident 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.6-5.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 counter-rotation 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 counter-rotating 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 top-right comer and bottom-right comer of the box with the flow directed downward along the hot wall. The eddy at the top-right corner is faster than the one at bottom-right 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 Liquid-Gas Systems Without Internal Heating


In liquid-gas 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 two-phase 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 counter-clockwise 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 counter-clockwise 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 two-phase 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 two-phase without internal heating
( case 2: both side walls isothermal, Grg,=103 )





80




5.4 Natural Convection in Liquid-Gas Systems with Internal Heating


The flow and heat transfer characteristics of a two-phase 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 counter-clockwise 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 property-based 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 counter-rotating 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 two-phase 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 two-phase with internal heating
(case 2: cylindrical container, side wall isothermal, GrgE =105 , Grg, =1.58x106)













CHAPTER 6
BULK EVAPORATION/CONDENSATION WITH INTERNAL HEATING
AT MICRO-GRAVITY 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 liquid-vapor 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 liquid-vapor interface does not change significantly, because the curvature of the interface is moderate. The liquid-vapor 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 g-levels, 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 UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID E9296YWM7_FHLQB6 INGEST_TIME 2015-03-25T20: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 Solid-Liquid Phase Change Problems 9 2.2 Numerical Modeling of Solid-Liquid 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 E-Based and T-Based 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 SINGLE-PHASE AND TWO-PHASE SYSTEMS 59 5.1 External Natural Convection in a Square-Cavity Benchmark Problem 60 5.2 Natural Convection in a Square-Cavity with Internal Heating 67 5.3 Natural Convection in Liquid-Gas 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 Liquid-Gas 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 MICRO-GRAVITY 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 Micro-Gravity (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 multi-phase 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 two-phase cell operated at 0-g 46 4.2 The control volume contacting other material 50 4.3 Parameter for heat generation rate 52 4.4 The evolution of liquid-vapor interface (0-g) 56 4.5 The evolution of temperature field (0-g) 57 4.6 Temperature surf of the system at steady state (0-g) 58 4.7 Variation of saturation temperature (0-g) 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 two-phase without internal heating ( case 1) 78 5.10 Natural convection of two-phase without internal heating ( case 2) 79 5.11 Natural convection of two-phase with internal heating ( case 1 ) 82 5.12 Natural convection of two-phase with internal heating ( case 2) 83 6.1 Schematic diagram of mass reflux 96 6.2 The evolution of liquid-vapor 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 liquid-vapor 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 liquid-vapor interface ( 1-g) 109 6.11 The evolution of isotherms ( 1-g) HO 6.12 Temperature surf at x = 2.0 ( 1-g) HI 6.13 Variation of saturation temperature ( 1-g ) Ill 6. 14 Flow patterns in vapor and liquid phases x = 2.0 ( 1-g ) 112 6.15 Velocity U profile at different axial positions (1-g) 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 Co-Chairman: 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 multi-phase nuclear fuel element at zero gravity, micro-gravity, and normal gravity is investigated. In the simulation, the numerical solutions have revealed much useful information, including the evolution of the bulk liquid-vapor phase change process, the evolution of the liquid-vapor 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 multi-phase 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 E-based and T-based methods, is investigated. Simulation of a one-dimensional conduction-driven 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 two-phase 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 non-mechanical 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 thermionic-converter 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 multi-phase 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 multi-phase 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 liquid-vapor 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 multi-phase nuclear fuel system A cylindrical container filled with multi-phase nuclear fuel, such as UF 4 , is shown in Figure 1.2. The fuel experiences liquid-vapor 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 zero-gravity condition, buoyancy-driven convection is nonexistent. However, under micro-gravity or one-gravity fields, buoyancy-induced 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 multi-phase 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 solid-liquid 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 solid-liquid 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 solid-liquid 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 buoyancy-driven 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 liquid-vapor 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 multi-phase nuclear system, including the multi-phase 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 one-dimensional, externally heated or cooled and conduction-driven 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 liquid-vapor interface topology are demonstrated. Chapter 5 describes the modeling of the combined internal and external natural convection in the single-phase system and two-phase 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 bench-mark 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 buoyancy-induced 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 Solid-Liquid 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 Solid-Liquid 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 gravity-flow 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 parabolic-type 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 start-up from the frozen state by using one-dimensional 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 two-phase transient flow and heat transfer characteristics within a constant pressure reservoir of a capillary-pumped 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 buoyancy-driven 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 post-accident 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 solid-liquid 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 1-D 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 solid-liquid phase change problems, such as those involved in melting and solidification processes. There has also been an ongoing interest in the investigation of liquid-vapor 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 El-Genk (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 start-up from the frozen state by using a system of simplified one dimensional equations. Shyy (1994b) has successfully developed a computational method for predicting the two-phase transient flow and heat transfer characteristics within a constant pressure reservoir of a capillary-pumped-loop 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 W-444 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 liquid-vapor 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 multi-dimensional phase-change 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 phase-change 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 = 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 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 e-c vl T /=o e = c v.,?;„+/Ae 0
PAGE 36

P,^.,T) = y(K,VT) / = o Pfjfa&mV.WTi-ptjJjL 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 Clausius-Clapeyron 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 phase-change 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 liquid-vapor 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,=pM-U v ) (3.10) The energy-balance 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 *L7K--L-t*> (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 T-based update method. This technique has proven to be very effective for achieving convergence. An alternative formulation to the T-based 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 H-based 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-e, (3 23) This is used for the iterative update of the vapor phase fraction. This update procedure is referred to as the E-based 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 E-based method compared to the T-based 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 E-based 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 E-based 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 liquid-vapor phase change occurring far away from its critical point, the latent heat is relatively large and the application of the E-based update method even for simple geometry is recommended. In the present work, the T-based and E-based methods are examined in detail from the viewpoint of accuracy and computational efficiency. A straight-forward implementation of the T-based 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 T-based 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 under-relaxation 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 E-based method is relatively straightforward and does not need any pre-conditioning. 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 non-uniform 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 non-uniform variation can be seen due to the nature of the discretization and the non-uniform 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 liquid-vapor 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 E-based method; the dashed line: the T-based 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 E-Based and T-Based Up d ate Methods The T-based update method and the E-based update method are examined for the same testing case. Figures 3.5(a) and 3.5(b) show that the T-based method and E-based method yield nearly identical results. It is seen that convergence for the E-based method requires about half the iterations which is needed for the T-based method. Thus the Ebased method is computationally more efficient than the T-based 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 T-based is better than the H-based method as reported by Shyy and Rao (1994). The reason for the improved efficiency of the E-based method is that the phase change temperature window is used for the T-based method, which is not required for the E-based method. When the T-based 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 T-based method converges more slowly. When the E-based 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 T-based 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 E-based method can completely circumvent the need for a phase change window Therefore, the Ebased method is superior to the T-based 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 E-based 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 E-based 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 evaporation-condensation 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 vapor-liquid 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 well-established 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 E-based and T-based methods, are investigated. The E-based method is superior to the T-based method in terms of computational efficiency and formulation simplicity for evaporation and condensation problems. The E-based 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 E-based 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 liquid-vapor phase change processes to a unified three-phase model. All phase change processes between these three phases, that is, solid-liquid, liquid-vapor, solid-vapor, or solid-liquid-vapor, 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 0-g 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 liquid-phase and condenses the vaporphase. For this system with such an orientation and operated under 0-g, there is no bubble or droplet generation on the wall surface. The temperature along the liquid-vapor interface does not change significantly, because the curvature of the interface is moderate. The liquid-vapor interface can be reasonably assumed to be isothermal, consequently the thermal capillary effect can be neglected. Hence the phase-change 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 non-uniform 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 phase-change (b) at the final steady state process Figure 4. 1 Schematic diagram of the two-phase cell operated at 0-g

PAGE 62

4. 1 The Governing Equations The internal energy formulation of the energy equation for this phase-change 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) Pv|fe,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 , e-e. This is the E-based update method discussed in chapter 3 and used for the iterative update of the vapor fraction. The T-based method does not work well for this twodimensional bulk evaporation and condensation problem, while the E-based always works quite well. Some additional closure relations, the same as Eqs.(3.7-8) 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/z-J 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 liquid-vapor 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 phase-change 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 liquid-vapor 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 liquid-vapor 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 multi-phase and internally heated system. The present results will be used as reference for comparison with the performance and features of this system operated under micro-gravity 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 liquid-vapor interface (0-g)

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 (0-g)

PAGE 73

Figure 4.6 Temperature surf of the system at steady state (0-g) 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 (0-g)

PAGE 74

CHAPTER 5 NATURAL CONVECTION WITH INTERNAL HEATING IN SINGLE-PHASE AND TWO-PHASE 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 non-zero 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 co-existence 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 single-phase and two-phase systems are investigated. A bench-mark problem, natural convection in a square cavity of single-phase without internal heating, is computed. Comparisons of the present solutions with accurate published bench-mark solutions are made to verify the computational accuracy. The natural convection in a square cavity of single-phase with internal heating is also computed. Qualitative comparisons with the published experimental results are also made. Then natural convection in a liquid-gas system without phase-change, 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 Square-Cavity 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 buoyancy-driven 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, buoyancy-driven flow in a square cavity with vertical side walls heated differentially, a so-called 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 two-dimensional 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(T-T )\ (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), = (T-TJ/(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 second-order central-difference 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 skew-symmetry 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 Square-Cavity 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 post-accident 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.6-5.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 counter-rotation 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 counter-rotating 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 top-right corner and bottom-right corner of the box with the flow directed downward along the hot wall. The eddy at the top-right corner is faster than the one at bottom-right 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 Liquid-Gas Systems Without Internal Heating In liquid-gas 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 two-phase 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 counter-clockwise 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 counter-clockwise 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 two-phase 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 two-phase without internal heating ( case 2: both side walls isothermal, Gr g =\0 3 )

PAGE 95

5.4 Natural Convection in Liquid-Gas Systems with Internal Heating The flow and heat transfer characteristics of a two-phase 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 counter-clockwise 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 property-based 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 counter-rotating 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 two-phase 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 two-phase 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 MICRO-GRAVITY 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 liquid-vapor 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 liquid-vapor interface does not change significantly, because the curvature of the interface is moderate. The liquid-vapor 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 g-levels, 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 u-momentum equation to account for buoyancy, according to the Boussinesq approximation A-ZkB-ACT-©] (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 one-g 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 = -, p--BjT °' ?,.»' 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 liquid-vapor 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 under-relaxation and more inner-loop 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 10-20, the central difference scheme failed. A multi-scheme 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 Micro-Gravity HO" 3 -g ) The evaporation and condensation processes under micro-gravity 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 liquid-vapor 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 phase-change 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 liquid-vapor 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 micro-gravity conditions, and due to the weak gravity effect, the buoyancy induced convection is not strong but it does change the overall characteristics of the phase-change 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 0-g 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 micro-gravity 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 liquid-vapor 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 multi-scheme 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 liquid-vapor interface with a coarse mesh, 51x21. Compared to that in Figure 6.2, the resultant topology of the liquid-vapor 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 liquid-vapor 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 liquid-vapor 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 liquid-vapor interface. A counter-current 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 micro-gravity 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 0-g 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 ( 1-g ) 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 liquid-vapor 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 0-g 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 1-g 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 liquid-vapor 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 liquid-vapor interface, the temperature distribution at the final steady state under three cases, 0-g, 10" 3 -g and 1-g, are plotted together. The gravity influences on the thermal performance of the bulk liquid-vapor 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 0-g, the thermal performance is controlled only by the heat transfer, such as internal heat generation and side wall cooling. At non-zero 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 multi-phase nuclear fuel elements to be operated under different gravity conditions, particularly when the experiments at micro-gravity environment are expensive and not easy to be conducted. The current model can also be applied to some two-phase flow problems with or without phase change, for example, film condensation, the liquid-vapor or liquid-gas 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 liquid-vapor 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 liquid-vapor interface at z= 2.0 0-g lO" 3 -g (b) The temperature distribution at x= 2.0 1-g 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 liquid-vapor 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 E-based and T-based methods, are investigated The E-based method is superior to the T-based method in terms of 115

PAGE 131

116 computational efficiency and formulation simplicity for evaporation and condensation problems. The T-based method is intractable for two-dimensional bulk evaporation and condensation problems, while the E-based method always works quite well. 4. A one dimensional conduction-driven bulk evaporation and condensation problem is simulated. The evolution of the bulk liquid-vapor 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 liquid-vapor 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 single-phase 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 two-phase 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 liquid-vapor 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 multi-phase nuclear fuel element under micro-gravity and normal gravity conditions is investigated. The numerical simulation elucidates the evolution of the bulk liquid-vapor phase change process. This includes the evolution of the liquid-vapor 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 multi-phase nuclear fuel cell on the gravity conditions. The differences of the phase change performance of the multi-phase 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 0-g, the performance is controlled only by the heat transfer, such as internal heat generation and side wall cooling. At non-zero 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 multi-phase nuclear fuel elements to be operated under different gravity conditions,

PAGE 133

118 particularly when the experiments at micro-gravity 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 liquid-vapor phase change problems to a unified three-phase model. All phase change processes between these three phases, such as solid-liquid, liquid-vapor, solid-vapor, or solid-liquid-vapor, could be simulated. 3 . The current model can also be applied to some two-phase flow problems with or without phase change, for example, the film condensation, the liquid-vapor or liquid-gas 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 . 855-866 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, 881-911 Buckle, V. & Peric, M. (1992), Numerical simulation of buoyant and thermocapillary convection in a square cavity, Numerical Heat Transfer. Part A . 21 . 121-141 . Carey, V. P. (1992), LiquidVapor Phase-Change Phenomena . Hemisphere Pub. Corp., New York Carpenter, B. M. & Homsy, G. M. (1989), Combined buoyant-thermocapillary flow in a cavity . J. Fluid Mech .. 207 . 121-132 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, 249-264 Ding, Z. & Anghaie, S. (1994a), Numerical investigation of the two-phase 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 R-12 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 transfer-A review of 1989 literature, Intl. J. Heat and Mass Transfer . 33, No. 12, 2349-2437 Eckert, E. R. G, Goldstein, R. J., et al. (1991), Heat transfer-A review of 1990 literature, Intl. J. Heat and Mass Transfer . 34. No. 12. 2931-3010 Eckert, E. R. G, Goldstein, R. J , et al. (1992), Heat transfer-A 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. 511-517.

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 . 174-181 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, 87-104 Gebhart, B., Jaluria, Y., Mahajan, R. L. & Sammakia, B.(1988), Buoyancy-Induced Flows and Transport . Hemisphere, Washington D. C. Gerner, F. M. & Tien, C. L.(1989), Axisymmetric Interfacial Condensation Model, J. Heat Transfer . 111 . 503-510 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 two-phase closed thermosyphons including the falling condensate film, J. Heat Transfer . 116 . 418-426 Hortmann, M., Peric, M. & Scheuerer, G. (1990), Finite volume multigrid prediction of laminar natural convection: bench-mark solutions, Intl. J. Numerical Meth. Fluids . 11. 189-207 Hsieh, C. K., Akbari, M. & Li, H. (1993), Solution of inverse Stefan problems by a source-and-sink method, Intl. J Numerical Methods Heat and fluid Flow . 2, No.5, 391-406 Hsieh, C. K. & Choi, C.-Y. (1992), Solution of oneand two-phase melting and solidification problems imposed with constant or time-variant temperature and flux boundary conditions, J. Heat Transfer . 114 . 524-528 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, 701-712

PAGE 137

122 Jang, J. H., Faghri, A., Chang, W. S. & Mahefkey, E. T.(1990), Mathematical modeling and analysis of heat pipe start-up from the frozen state, J. Heat Transfer . 112 , 586-594 Jones, W. P. & Renz, U. (1974), Condensation from a turbulent stream onto a vertical surface, Int. J. Heat Mass Transfer , 17, 1019-1028 Kimura, S. & Bejan, A. (1983), The "heatline" visulization of convective heat transfer, J. Heat Transfer , 105 . 916-919 Lame, G. & Clapeyron, (1831) Ann. Chem. Phvs .. 47, 250-256 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 . Conf-930103, 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 . 345-349 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 . Conf-9101 16, M. S. EI-Genk & M. D. Hoover (Eds), American Institute of Physics, New York Li, H. (1993), Source-and-Sink method of Solution of Twoand Three-Dimensional 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: 21-25 Nadarajah, A. & Narayanan, R. (1990), Comparison between morphological and Rayleigh-Marangoni instabilities, Meinkoen, D. & Haken, H. (Eds), Dissipative

PAGE 138

123 Structures in Transport Processes and Combustion . 215-228, SpringerVerlag, New York Nusselt, W. Z. (1916), Die oberflachenkondensation des wasserdampfes, Z. Ver. Distch Ing., 60, 541-569 Ostrach, S. (1982), Low-gravity fluid flows, Ann. Rev. Fluid Mech .. 14. 313-345 Ostrach, S. (1983), Fluid mechanics in crystal growth-The 1982 freeman scholar lecture, J. Fluids Engineering, 105 . 5-20 Ostrach, S. (1988), Natural convection in enclosures, J. Heat Transfer . 110. 1 175-1 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 . 562-570 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, 2690-94 Ramaswamy, B. & Jue, T C. (1992), Analysis of thermocapillary and buoyacy-affected cavity flow using FEM, Numer. Heat Transfer . 22, 379-399 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. McGraw-Hill, New York Samarskii, A. A , Vabishchevich, P. N, Iliev, O. P. & Churbanov, A. G. (1993), Numerical simulation of convection/diffusion phase change problems-A Review, Int. J. Heat Mass Transfer . 36, 4095-4106

PAGE 139

124 Shamsundar, N. & Sparrow, E. M. (1975), Analysis of multidimensional conduction phase change via the enthalpy model, J. Heat Transfer , 97, 333-340 Shyy, W. (1994a), Computational Modeling for Fluid Flow and Interfacial Transport . Elsevier, Amsterdam, The Netherlands Shyy, W. (1994b), Modeling of transient two-phase heat transfer for spacecraft thermal management, Microgravitv Sci. Tech .. 7, 219-227 Shyy, W. & Chen, M. (1990), Steady-state natural convection with phase change Int J Heat Mass Transfer . 33, 2545-2563 Shyy, W. & Chen, M. (1991), Interaction of thermocapillary and natural convection flow during solidification: normal and reduced gravity conditions, J. Crystal Growth . 108 . 247-261 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. 1229-1245 Shyy, W. & Rao, M. (1992b), Convection treatment for high Rayleigh number, laminar, natural convection calculation, Numerical Heat Transfer. Part B . 22, 367-374 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. 1833-1844 Shyy, W., Gringrich, W. K., Krotiuk, W. J & Fredley, J. E. (1993b), Transient twophase heat transfer in a capillary-pumped-loop 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 , 946-954 Shyy, W. & Rao, M. M. (1994), Enthalpy based formulations for phase-change problems with application to g-jitters, Microgravitv Sci. Technology , 7, 41-49

PAGE 140

125 Shyy, W., Liang, S.-J. & Wei, D. Y. (1994), Effect of dynamic perturbation and contact condition on edge-defined fiber growth characteristics, Int. J. Heat Mass Transfer. 37, 977-987 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, 430-436 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. 1829-1845. Stefan, J. (1891), Ann. Phvs. Chemie . 42, 269-286. Tournier, J.-M. & El-Genk, M. S. (1992), "HPTAM" Heat-pipe 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 convection-diffusion mushy region phase-change problems, Int. J. Heat Mass Transfer . 30. 1709-1719 Voller, V R. & Swaminathan, C. R. (1991), General source-based method for solidification phase change, Numerical Heat Transfer. Part B . 19, 175-189 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, 23-85 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, Co-Chair 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