DEVELOPMENT OF A Tl'10-~FLUID, TWO-PHASE MODEL FOR LIGHT WATER REACTOR SUBCHANNEL ANALYSIS by JOHN EDWARD KELLY B.S.N.E., University of Michigan, (1976} S.H., Massachusetts Institute of Technology, (1978) SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY at the MASSACHUSETTS INSTITUTE OF TECHNOLOGY AUGUST 1980 @ Massachusetts Institute of Technology 1980 Signatur.e of Author - ~~~nat~-~e ..... _9_<1~8LIIIIC.,..te_d _____ _ Depart~t of Nuclear Engi~ng, August 8, 1980 Signature redacted Certified by ,,,- ""' -------o~-------0--#---T-h_e_s_i_s_S--up-e-.r-v_i_s_o_r Signature redacted Accepted by ___ C_h_a_i_rm_a_n-,-b--D-e-~_..~-rt--~-e-rn-;-C-;nnn-----g-t-t-.e-e_o_n_G_r_a_d_u_a_t_e_S_t_u_d~e-n_t_s IRcHtV((J rt-""'"'A "' "· •Mu.) o~~T~~~1J&gi51;;u~ NOV 3 19SO 1.IBRARIES - 2- DEVELOPMENT OF A TWO-FLUID, TWO-PHASE MODEL FOR LIGHT WATER REACTOR SUBCHANNEL ANALYSIS by John Edward Kelly Submitted to the Department of Nuclear Engineering on August 8, 1980 in partial fulfillment of the requirements for the Degree of Doctor of Philosophy in Nuclear Engineering ABSTRACT The-problem of developing and assessing the two-fluid model computer code THERMIT for light water reactor (LWR) subchannel analysis has been addressed. The developmental effort required a reformulation of the coolant-to-fuel rod coupling so that THERMIT is now capable of traditional coolant-centered subchannel analysis. A model that accounts for mass, momentum and energy transport between mesh cells due to turbulent mixing for two-phase conditions has also been introduced. This model is the first such attempt in a two-fluid context. The liquid-vapor interfacial exchange terms in the two-fluid model have been modified for improved accuracy. A systematic evaluation of the exchange models has been performed. The mass and momentum exchange rates between the vapor and the liquid for pre-CHF conditions were evaluated by comparison to void fraction data in over 30 one-dimensional steady-state experiments reported in the open literature. The liquid- vapor energy exchange rate for post-CHF conditions was assessed using 15 steady-state, one-dimensional wall temperature measurements also found in the open literature. The two-phase mixing model has been evaluated using G.E. and Ispra BWR and PWR rod-bundle measurements. Comparisons with these measurements have shown the appropriateness of this model. The assessment of the wall-to-coolant heat transfer model used steady-state, one-dimensional as well as transient, three-dimensional measurements. The comparisons resulted in a modification of the Biasi CHF correlation to imptove 'its predictions. THERMIT has been shown to accurately predict the thermal- hydraulic behavior of rod-bundles. Thus, it represents the first two- fluid computer code with this proven capability. Thesis Supervisor: Mujid S. Kazimi Title: Associate Professor of Nuclear Engineering Thesis Reader: John E. Meyer.. Title: Professor of Nuclear Engineering - 3- ACKNOWLEDGEMENTS The author would like to express his sincere gratitude to all those who assisted him in the preparation of this thesis. In particular, special recognition is given to the following individuals. To Professor Mujid Kazimi who, as my thesis advisor has generously given his time and advice. His valuable comments and suggestions have greatly contributed to this thesis. To Professor John Meyer who, as my thesis reader, has been very helpful during the course of this work. To Professor Lothar Wolf, whose encouragement and advice have been greatly appreciated. To Don Dube and Andrei Schor, who have shared with me many instructive discussions concerning THERMIT. To Ping Kao, Samir Mohammed, Jim Loomis and Mahmoud Massoud who have assisted me in the computer related activities. To Gail Jacobson and Riva Esformes who have typed this thesis. To my wife, Sue, and all of my family whose constant encouragement and support have led to the successful completion of this work. Finally, the financial support of Boston Edison Company, Consumers Power Company,.Northeast Utilities Service Company, Public Service Electric and Gas.Company, and Yankee Atomic Electric Company under the sponsorhsip of the M.I.T. Energy Laboratory is greatly appreciated. -4- TABLE OF CONTENTS Section ABSTRACT ACKNOWLEDGEMENT S TABLE OF CONTENTS LIST OF FIGURES LIST OF TABLES NOMENCLATURE CHAPTER 1 - INTRODUCTION 1.1 Background 1.2 Research Objective 1.3 Development Approach CHAPTER 2 - REVIEW OF PREVIOUS WORK 2.1 Introduction 2.2 Rod Bundle Analysis Techniques 2.2.1 Subchannel Analysis 2.2.2 Distributed Resistance Models 2.3 Two-Phase Flow Models 2.4 Description of THERMIT 2.4.1 Background 2.4.2 General Characteristics 2.4.3 Two-Fluid Model Conservation Equations 2.4.4 Finite Difference Equations 2.4.5 Constitutive Equations CHAPTER 3 - DEVELOPMENT OF THERMIT SUBCHANNEL ANALYSIS CAPABILITY 3.1 Introduction Page 2 3 4 7 12 14 17 17 19 21 23 23 29 29 32 33 35 35 35 38 40 49 52 52 -5- Section 3.2 Geometrical Modeling Capability 3.3 Two-Phase Turbulent Mixing Model 3.3.1 Background 3.3.2 Model Formulation 3.3.2.1 Background 3.3.2.2 Analytical Formulation and Discussion 3.3.2.3 Numerical Scheme CHAPTER 4 - DEVELOPMENT AND ASSESSMENT OF THE LIQUID- VAPOR INTERFACIAL EXCHANGE MODELS 4.1 Introduction 4.2 Assessment Strategy 4.3 Interfacial Mass Exchange 4.3.1 Background 4.3.2 Subcooled Vapor Generation Model 4.3.3 Droplet Vaporization Model 4.4 Interfacial Energy Exchange 4.5 Interfacial Momentum Exchange CHAPTER 5 - ASSESSMENT OF THE TWO-PHASE MIXING MODEL 5.1 Introduction 5.2 G.E. 9 Rod Bundle Tests 5.2.1 Test Description 5.2.2 Single-Phase Comparisons 5.2.3 Uniformly-Heated Cases 5.2.4 Evaluation of Mixing Parameters 5.2.5 Non-Uniformly Heated Cases Page 53 56 56 60 60 62 69 71 71 72 75 75 80 94 105 114 138 138 145 145 150 152 156 161 -6- Section 5.3 Ispra BWR Tests 5.3.1 Test Description 3.3.2 Results 5.4 Ispra 16 Rod PWR Cases 5.4.1 Test Description 5.4.2 Results 5.5 Conclusions CHAPTER 6 - THE HEAT TRANSFER MODEL 6.1 6.2 6.3 CHAPTER 7.1 Introduction Modifications Assessment 6.3.1 Steady-State Results 6.3.2 Transient Results 7 - CONCLUSIONS AND RECOMMENDATIONS Conclusions 7. 2 Recommendations REFERENCES APPENDIX A DERIVATION APPENDIX B TWO-PHASE M APPENDIX C HEAT TRANSFI OF THERMIT GOVERNING EQUATIONS IXING MODEL ASSESSMENT RESULTS ER CORRELATIONS 163 163 167 175 175 177 186 188 188 194 199 200 217 232 232 239 242 247 267 280 -7- LIST OF FIGURES Figure 2.1 2.2 2.3 3.1 4.1 4.2 4.3 4.4 4.5 214-9-3 Fraction versus Enthalpy for Maurer 214-3-5 Fraction versus Enthalpy for Marchaterre 168 Fraction versus Enthalpy for Marchaterre 184 Temperature Comparisons for Bennett 5332 Temperature Comparisons for Bennett 5253 Temperature Comparisons for Bennett Case 4.6 Void Case 4.7 Void Case 4.8 Void Case 4.9 Wall Case 4.10 Wall Case 4.11 Wall Case 5442 4.12 Comparison of Wall Temperature Predictions Using Various r Models for Bennett Case 5332 4.13 Temperature Distributions for Subcooled Boiling and Droplet Vaporization Page 31 42 48 57 77 79 82 83 Coolant-Centered and Rod-Centered Layouts Typical Fluid Mesh Cell Showing Locations of Variables and Subscripting Conventions Typical Rod Arrangement in Transverse Plane Illustration of Fuel Rod Modeling Boiling Regimes in Two-Phase Flow in a Vertical Ttbe with Heat Addition Illustration of Vapor Generation Rate, F, versus Equilibrium Quality, X Illustration of Vapor Bubble Nucleation and Growth Typical Temperature Distributions in Subcooled Boiling Void Fraction versus Enthalpy for Maurer 101 102 103 104 107 90 91 92 93 -8- Figure Page 4.14 Predicted Liquid and Vapor Temperatures for l11 Maurer Case 214-3-5 4.15 Predicted Liquid and Vapor Temperatrues for 112 Bennett Case 5336 4.16 Typical Flow Patterns in Two-Phase Flow 115 4.17 Void Fraction versus Enthalpy for Maurer 129 Case 214-3-4 4.18 Void Fraction verius Enthalpy for Christensen 130 Case 12 4.19 Void Fraction versus Enthalpy for Christensen 131 Case 9 4.20 Void Fraction versus Enthalpy for Marchaterre 132 Case 185 4.21 Void Fraction versus Enthalpy for Christensen 134 Case 12 4.22 Comparison of Predicted Slip Ratios for 135 Christensen Case 12 4.23 Vapor Superficial Velocity versus Void Fraction 137 for Christensen Data 5.1 Cross Sectional View of G.E. 9 Rod Bundle Used 146 in Mass Velocity and Enthalpy Measurements 5.2 Radial Peaking Factors for Non-Uniformly Heated 148 Cases 5.3 Comparison of Measured and Predicted Mass 151 Velocities for G.E. Isothermal Tests 5.4 Comparison of Measured and Predicted Exit Quality 153 in Corner Subchannel for G.E. Uniformly Heated Cases 5.5 Comparison of Measured and Predicted Exit Quality 154 in Edge Subchannel for G.E. Uniformly Heated Cases 5.6 Comparison of Measured and Predicted Exit Quality 155 in Center Subchannel for G.E. Uniformly Heated Cases -9- Figure 5.7 Comparison of Measured and Predicted Mass Velocities for G.E. Uniformly Heated Tests 5.8 Comparison of G.E. 2E Cases with 6M Varied from I to 10 5.9 Comparison of THERMIT Predictions for Case 2E2 with Variations in KM 5.10 Comparison of Measured arid Predicted Quality for G.E. Non-Uniformly Heated Tests 5.11 Comparison of Measured and Predicted Mass Velocities for G.E. Non-Uniformly Heated Tests 5.12 Cross Sectional View of Ispra BWR Test Section 5.13 Comparison of Measured and Predicted Exit Quality in Subchannel 2 - Ispra BWR Tests 5.14 Comparison of Measured and Subchannel I versus Bundle Ispra BWR Tests 5.15 Comparison of Measured and Subchannel 5 versus Bundle Ispra BWR Tests 5.16 Comparison of Measured and Subchannel 4 versus Bundle Ispra BWR Tests 5.17 Comparison of Measured and Velocities for Subchannels BWR Tests Predicted Quality for Average Quality - Predicted Quality for Average Quality - Predicted Quality for Average Quality - Predicted Mass 1 and 2 - Ispra 5.18 Comparison of Measured and Predicted Mass Velocities for Subchannel 4 and 5 - Ispra BWR Tests 5.19 Cross Sectional View of Ispra PWR Test Section 5.20 Comparison of Measured and Predicted Exit Quality for Subchannel 1 versus Bundle Average Quality - Ispra PWR Tests 5.21 Comparison of Measured and Predicted Exit Quality for Subchannel 2 versus Bundle Average Quality - Ispra PWR Tests Page 157 159 160 162 164 165 169 170 171 172 173 173 176 179 180 -10-- Figure Page 5.22 Comparison of Measured and Predicted Exit 181 Quality for Subchannel 3 versus Bundle Average Quality - Ispra PWR Tests 5.23 Comparison of Measured and Predicted Exit 182 Quality for Subchannel 4 versus Bundle Average Quality - Ispra PWR Tests 5.24 Comparison of Measured and Predicted Exit 183 Quality for Subchannel 5 versus Bundle Average Quality - Ispra PWR Tests 5.25 Comparison of Measured and Predicted Mass 184 Velocities for Ispra PWR Tests 5.26 Comparison of Measured and Predicted Mass 185 Velocities for Ispra PWR Tests 6.1 Typical Boiling Curve 189 6.2 BEEST Heat Transfer Logic 193 6.3 Modified Heat Transfer Logic 198 6.4 Typical Wall Temperature Distribution for 201 Bennett Case 5273 6.5 Comparison of Measured and Predicted Pre-CHF 203 Wall Temperatures for Bennett Case 5332 6.6 Comparison of Measured and Predicted Pre-CHF 204 Wall Temperatures for Bennett Case 5276 6.7 Comparison of Measured and Predicted Pre-CHF 205 Wall Temperatures for Bennett Case 5253 6.8 Comparison of Measured and Predicted Pre-CHF 206 Wall Temperatures for Bennett Case 5394 6.9 Comparison of Measured and Predicted Pre-CHF 207 Wall Temperatures for Bennett Case 5451 6.10 Illustration of Mass Velocity and Quality 213 Dependence of Biasi CHF Correlation 6.11 Cross Sectional View of G.E. 9 Rod Bundle Used 219 in Transient CHF Tests 6.12 Comparison of MCHFR Predictions versus Time for 221 Case 9-151 -11- Figure Page 6.13 Comparison of MCRFR Predictions versus 222 Time for Case 9-170 6.14 Comparison of MCEFR Predictions versus 223 Time for Case 9-175 6.15 Comparison of MCHFR Predictions versus 224 Time for Case 9-179 6.16 Comparison of Biasi Correlation Predictions 226 for Case 9-151 with Modified Grid Coefficients 6.17 Comparison of Biasi Correlation Predictions with 228 a Tighter Convergence Criterion (Case 9-179) 6.18 Comparison of Measured and Predicted Maximum 229 Wall Temperatures for Case 9-179 6.19 Comparison of Measured and Predicted Maximum 231 Wall Temperatures for Case 9-151 A.1 Illustration of Control Volume Containing 252 Liquid, Vapor and Solid -12- LIST OF TABLES Tableage 2.1 Features of Some Thermal-Hydraulic Computer Codes 24 2.2 Thermal-Hydraulic Code Classification Criteria 25 2.3 Summary of Two-Phase Flow Models 27 2.4 Summary of Transport Processes 50 3.1 Implicit Heat Transfer Algorithm 55 4.1 Summary of Assessment P.rogram 74 4.2 Test Conditions for One-Dimensional Steady-State 89 Data 4.3 Test Conditions used to Develop Saha Correlation 97 for Post-CHF 4.4 Bennett Test Conditons for CHF in Tubes 100 4.5 Summary of Liquid-Vapor Intezfaciai Forces 121 4.6 Comparison of Viscous Force Coefficients 124 4.7 Comparison of Inertial Force Coefficients 126 5.1 Test Conditions for Rod-Bundle Experiments 143 6.1 Summary of Heat Transfer Correlations 192 6.2 CHF Comparisons for Bennett Test Cases 211 6.3 CHF Comparisons for Bennett Test Cases with 216 Corrected Biasi Correlation 6.4 Summary of Transient CHF Cases 218 7.1 Summary of Modifications and Improvements Made 233 in THERMIT A.1 Summary of Terms Used in Conservation Equations 250 B.1 Test Conditions for- Rod-Bundle Experiments 268 B.2 Comparison of Measured and Predicted Exit Mass 269 Velocities for Isothermal Tests in 9 Rod G.E. Tests Table B.3 Comparison of Measured and Predicted Exit Quality Distributions for Uniformly Heated 9-Rod G.E. Cases B.4 Comparison of Mass Velocity Heated Cases B.5 Comparison of Mass Velocity BWR Cases Measured and Predicted Quality and Distributions for G.E. Non-Uniformly Measured and Predicted Quality and Distributions for Ispra 16-Rod c0.1 Heat Transfer Correlations Page 270 274 276 281 NOMENCLATURE A Area Cd Drag Coefficient C Specific Heat D Diameter e Internal Energy F Gravitational Force F Vapor-Liquid Interfacial Momentum Exchange Rate Ft Turbulent Momentum Exchange Rate F Wall Frictional Force g Gravitational Constant G Mass Flux H Heat Transfer Coefficient i Enthalpy J v Superficial Vapor Velocity k Thermal Conductivity K Mixing Model Parameter M L Length P Pressure Pr Prandtl Number Qi Interfacial Heat Transfer Rate Qt Turbulent Heat Transfer Rate QW Wall Heat Transfer Rate i Power q"I Heat Flux R b Bubble Radius -15- Nomenclature (continued) Re Reynolds Number S Slip Ratio (V /V S ij Gap Spacing Between Coolant Channels t Time T Temperature Td Bubble Departure Temperature V Velocity VR Relative Velocity (Vv - V) W Volume W' Turbulent Mixing Rate WI" Turbulent Mass Flux t W t Mass Exchange Due to Turbulence X Quality a Void Fraction a Surface Tension Parameter Defined in Eq. (3.33) r Vapor Generation Rate p Density 11 Viscosity 0 Mixing Model Parameter E/ Turbulent Velocity Generalized Mixing Rate Term T Turbulent Shear Force 6 Droplet Diameter Nomenclature (continued) Subscripts i,j,k f g s v w Nodal Locations Saturated Liquid Saturated Vapor Liquid Saturation Vapor Wall Superscripts x,y,z Spatial Directions Anticipated Transient Without Scram Boiling Water Reactor Critical Heat Flux Critical Heat Flux Ratio Critical Power Ratio Departure from Nucleate Boiling Loss of Coolant Accident Light Water Reactor Minimum Stable Film Boiling Pressurized Water Reactor Acronyms_ ATWS BWR CHF CHFR CPR DNB LOCA LWR MSFB PWR - 17- 1.0 INTRODUCTION 1.1 Background Light water reactor safety research is ultimately aimed at ensuring that the public will not be adversely affected if any of a variety of anticipated or postulated reactor accidents should occur. This require- ment is met by specifying operational design limits that are based on conservative assumptions for the behavior of the reactor. Reactor safety research is primarily concerned with validating the appropriateness of these limits as well as assessing the margins- present in these limits. In order to study the normal and abnormal transient behavior of nuclear reactors, many complex-phenomena and systems need to be analyzed. One of the major areas which must be investigated.:is the thermal-hydraulics of the reactor system. Included here are the reactor core heat removal system, the secondary heat removal system (if present) as well as any auxillary systems which are related to removal of heat from the reactor. Since most of the radioactive inventory is contained within the reactor core, the preservation of the core integrity is essential. Moreove;,.the most likely radioactivity release mechanisms result from a thermally induced failure of the fuel rod cladding. Thus, the thermal-hydraulic behavior of the core is generally the most important consideration of reactor safety analysis. In order to meet the objective of accurately predicting the thermal- hydraulic field in the reactor core a number of analytical tools have been developed. These range from simple one equation models, used to predict a particular phenomenon, to large computer codes which attempt -18- to analyze the entire reactor system. Typically, the most widely used and generally the most useful tools are the thermal-hydraulic computer codes. Simply stated, these codes attempt to numerically solve the mass, momentum and energy conservation equations for a particular geometrical configuration and for the conditions of interest. Since the conservation equations must be supplemented by empirical correlations needed to describe specific phenomena, the thermal-hydraulic computer codes are engineering analysis tools which combine basic physics with empirical models. In the past few years, the need for improved analysis of nuclear reactor safety has lead to the rapid development of advanced methods for multidimensional thermal-hydraulic analysis. These methods have become progressively more complex in order to account for the many physical phenomena which are anticipated during both steady-state and transient conditions. In particular, the modeling of two-phase flow, which is required for both BWR and PWR systems, is especially complex. In two- phase flows, both thermal and mechanical non-equilibrium between the two phases can exist. These non-equilibrium effects take the form of sub- cooled boiling, vapor superheating and relative motion of the two phases. In order to have realistic calculations, these physical phenomena must be accounted for in the numerical method. The numerical methods must also be capable of analyzing the many flow patterns which occur in postulated transients. For example, in a loss of coolant accident (LOCA) or a severe anticipated transient without scram (ATWS), flow reversal or counter-current flow may occur in the reactor core. Elaborate solution techniques have been developed specifi- cally to be able to describe fluid fields with no restriction on speed or direction. Although it is improbable that a particular computer code be appli- cable for all transients, it is necessary that a code be able to analyze all anticipated flow conditions in problems for which the code is applied. The only practical way to realize the needed flexibility is to combine realistic physical models with unrestricitive numerical solution techniques. Hence, the trend in current thermal-hydraulic safety research is to pursue the development of such codes. 1.2 Research Objective As discussed in the previous section, the thermal-hydraulic computer codes play a key role in LWR safety analysis. However, due to the limita- tions of present day computers, precise details of the thermal-hydraulic behavior can only be determined for a relatively small region of the core. The response of the entire reactor can be determined if large control volumes are used. However, within these volumes information about the temperature and flow distribution is lost. If these distribu- tions are important for assessing the safety of the reactor, then detailed modeling is required. By using smaller control volumes, for example sub- channels, sufficient temperature and flow information can be determined, but only for limited regions of the core. For instance, the largest region which might be analyzed on a subchannel basis would be one BWR 8 x 8 assembly. Nevertheless, if the limiting region of the core can be identified, then this type of detailed analysis is sufficient to evaluate the safety of the reactor. A number of power and flow transients do require detailed subchannel modeling, particularly in the hottest part of the core. However, -20- previous computer codes, which have used subchannel modeling, have either lacked a realistic two-phase flow model (e.g. COBRA IV [1]) or lacked an unrestricted solution technique (e.g. COBRA-IIIC [2]). Consequently, the applicability of the previous codes is somewhat limited. In view of the shortcomings of the previous codes, a new code which does not suffer from these deficiencies has been developed. Using the computer code THERMIT [31 as a framework, the present developmental effort has expanded the capabilities of THERMIT such that the new version of THERMIT can successfully analyze subchannel geometry [4]. THERMIT has been selected for this project due to its two-phase flow model and solution technique. The two-fluid, two-phase model which is used in THERMIT realistically allows for thermal and mechanical non- equilibrium between the phases. This feature permits description of the complex phenomena encountered during transients. The solution technique is a modification of the ICE method [5,6], and is capable of predicting the flow conditions with minimum restrictions. The primary application of this new version is transient analysis of LWR rod bundles on a subchannel basis. Although depressurization tran- sients (i.e. LOCA) have not been excluded as possible applications of this code, these transients are not the primary type of transient under consid- eration. Rather, anticipated or near-operational power and flow transients are the main focus of the present development considerations. By concen- trating on non-depressurization transients, the code can be validated for several practical conditions. Furthermore, the proper analysis of a LOCA generally requires modeling the entire reactor system and THERMIT has been designed for core analysis only. Consequently, the applications for this new tool are limited to cases which can be analyzed by modeling only the core and supplying appropriate core flow boundary conditions. Neverthe- less, these cases represent a large number of problems with practical interest. With the ability to analyze subchannels, THERMIT is the first two- fluid model code with this capability. Due to the advanced treatment of the two-phase flow and reliable solution method, this code represents a significant addition to the area of rod-bundle thermal-hydraulic analysis. Other multi-fluid codes that may be used for subchannel applications are still under development at Argonne National Laboratory (COMMIX-2 (7-1) and Battle Pacific Northwest Laboratory (COBRA-TF [8]). 1.3 Development Approach The development of this new version of THERMIT has been accomplished using the following strategy: (1) Modify the code structure and numerical method as necessary, (2) Verify, extend, and assess the constitutive models, (3) Assess numerical properties of the code, and (4) Implement improved models as necessary. This strategy is actually iterative in nature. That is, as the need for improved models is found, code modifications and assessment are subsequent- ly required. Hence, the above steps overlap one another. This development can be divided into two main steps. The first step has involved modifying the original version of THERMIT so that subchannel geometry could be analyzed. This modification has affected both the geo- metrical modeling capability as well as the physical modeling. The geo- metrical modeling changes were required so that the traditional coolant- centered subchannels might be analyzed with-THERMIT. The changes in the -21- -22- physical modeling wer necessary to account for turbulence effects in single phase and two-phase flows for rod-bundle analysis. After reviewing previous work in Chapter 2, the significant work related to this modifica- tion effort is discussed in Chapter 3. After implementing the capability for subchannel analysis, the second step has been the validation and assessment of the code. A strategy has been adopted which allows for independent assessment of the various con- stitutive models using open literature experimental measurements. Measure- ments typical of expected subchannel conditions have been compared with the code's predictions in this effort. These comparisons are useful for both validating the predictive capability of the code as well as identi- fying areas which require improvement. The net result of this assessment effort is that the code can be used with confidence for subchannel appli- cations. The results of this assessment are discussed in Chapters 4, 5 and 6. A listing of the actual computer code will not be given here due to its length. Rather, the interested reader is referred to reference 9 which contains detailed information on the usage of this new version of THERMIT. Sample problems as well as input instructions are given in this reference. -23- 2.0 REVIEW OF PREVIOUS WORK 2.1 Introduction Nuclear reactor thermal-hydraulic safety research encompasses both experimental and analytical investigations. The experimental research attempts to measure and identify the important variables in both single- phase and two-phase flows. The analytical research attempts to develop methods which numerically solve the equations describing the heat transfer and fluid dynamics in a reactor. Elabocate numerical methods have evolved which rely heavily on the use of digital computers. Conceivably, if all of the significant physical phenomena are considered in the computer code, then accurate predictions of the flow conditions can be obtained. These-methods can also analyze conditions which could not be directly measured. The only practical limitation of these methods is the problem size which a computer can accomodate in terms of both storage and execution times. Since this thesis has been concerned with the development of the thermal-hydraulic computer code, THERMIT, it is instructive to review the general characteristics of nuclear reactor thermal-hydraulic codes. The key features of a few of these codes are presented in Table 2.1. As discussed by Massoud [13j, it is possible to classify the codes according to the criteria summarized in Table 2.2. The first major division is re- lated to the code's capability to handle one component or the entire by- draulic loop. Loop codes analyze a number of components simultaneously and, TABLE 2.1 Features of Some Thermal-Hydraulic Computer Codes Computer Type of Method of Two-Phase Solution Technique Code Analysis Analysis Flow Model COBRA IIIC [2] Component Subchannel Homogeneous Equilibrium Marching Method COBRA IV [1] Component Subchannel Homogeneous Equilibrium Marching Method or I.C.E. Method WOSUB [10] Component Subchannel Drift Flux Marching Method COMMIX-2 [7] Component Distributed Two-Fluid I.C.E. Method Resistance THERMIT [3] Component Distributed Two-Fluid I.C.E. Method Resistance TRAC [12] Loop Distributed Two-Fluid or Drift Flux I.C.E. Method Resistance to -25- TABLE 2.2 Thermal-Hydraulic Code Classification Criteria 1. System Analysis Capability A. Loop Codes B. Component Codes i. Subchannel Analysis ii. Distributed Resistance Analysis iii. Distributed Parameter Analysis 2. Two-Phase Flow Model A. Homogeneous Equilibrium Model B. Drift-Flux Model C. Two-Fluid Model -26- consequently, analysis is not as detailed as found in the individual com- ponent codes. However, the component codes must have appropriate boundary conditions supplied from external calculations. This prevents accurate modeling of the coupling between the component and the rest of the system. Naturally, loop codes do not suffer from this problem. Component codes, specifically those intended for rod-bundle analysis, can be further classified according to their analysis method. This topic has been reviewed by Sha [14]. Three types of methods can be identified; subchannel analysis, distributed resistance analysis and dis- tributed parameter analysis. Each of these analytical techniques has certain advantages and disadvantages relative to the other methods. Sub- channel analysis techniques permit fairly detailed analysis of the flow, but are limited by inherent assumptions concerning the flow. Distributed resistance methods can analyze either large or small regions but require the accurate determination of the flow resistances. The distributed parameter analysis method gives the most detail of the flow structure, but is limited to small regions. All gf the core component codes use one of these three analysis techniques. The second major division is the type of two-phase flow model. The important types of models which have been incorporated into thermal- hydraulic codes include the homogeneous equilibrium model (HEM), the drift-flux model, and the two-fluid model. Essentially, the type of two-phase model refers to the number of conservation equations which are used to describe the two-phase flow, as summarized in Table 2.3. As the number of conservation equations increases, the number of constitutive models also increases. However, with more equations, accurate results are more likely to be predicted for severe conditions. The more general TABLE 2.3 Summary of Two-Phase Flow Models [151 Two-Phase Conservation Equations Constitutive Laws Imposed Restrictions On: Flow Model Mass Energy Momentum F Q r Q Fi Phasic Phasic Phatic W W Temperatures Velocities Distributions Homogeneous 1 1 1 1 1 0 0 0 Equal Equal Uniform Drift Flux 2 1 1 1 1 1 0 0 T -T Slip None Relahion Relation 4 Equation 1 2 1 1 2 0 1 0 None Slip Uniform Relation Models 1 1 2 2 1 0 0 1 T -T2 None Uniform Relation Drift Flux 2 2 1 1 211 1 0 None Slip None Relation 5 Equation 2 1 2 2 1 1 0 1 T -TR None None Relation Models 1 2 2 2 2 0 111 None None Uniform NJ TABLE 2. 3 (continued) Two-Phase Conservation Equations Constitutive Laws Imposed Restrictions On: Flow Model Mass Energy Momentum F 1 P Q F Phasic Phasic Phasic Temperatures Velocities Distributions Two-Fluid 2 2 2 2 2 1 1 1 None None None Three-Fluid 3 3 3 3 3 2 2 2 None None None (Liquid, Vapor and Liquid Drops) Go -29- equations also allow better physical modeling which is essential for the description of two-phase flow. Since the present work is concerned with the application of the two- fluid model code THERMIT to detailed rod-bundle analysis, it is useful to discuss in detail both the type of analysis techniques and the two-phase flow models. The analysis techniques are discussed in Section 2.2 while the two-phase flow models are discussed in Section 2.3. Following this discussion, the specifics of the THERMIT computer code will be given in Section 2.4. 2.2 Rod-Bundle Analysis Techniques As discussed in the previous section, three types of techniques are available for rod-bundle analysis. These include subchannel analysis, distributed resistance analysis and distributed parameter analysis. The distributed parameter methods are limited to very small regions and will not be discussed here. The other two methods, however, are very useful for analyzing the entire rod-bundle and are discussed in detail. 2.2.2 Subchannel Analysis Of all the methods developed for analyzing the thermal-hydraulic behavior in complex rod-bundle geometry, the subchannel method has been found to be particularly well-suited. Weisman and Bowring [16] and Rouhani (17] have reviewed this type of analysis and present the following view of traditional subchannel analysis. In this method, the rod-bundle cross section is subdivided into a number of parallel interacting flow subchannels. Conventionally these subchannels are defined by lines joining the fuel rod center (see Figure -30- 2.la). This choice is somewhat arbitrary and other choices are possible, such as the lines of zero shear stress (see figure 2.lb) [12]. This latter type of subchannel is referred to as a rod centered subchannel while the former type is called a coolant centered subchannel. Once the radial plane has been defined, each subchannel is divided axially into a number of intervals (nodes) which are typically between 8 and 30 cm long. For each node, which can be thought of as a control volume, a set of mass, energy, axial momentum and transverse momentum conservation equations are written and solved with an iterative technique. The main assumptions of this method are: (1) The detailed velocity and temperature destributions within a subchannel are ignored; (2) The transverse nomentum equation is simplified due to the assumption of predominantly axial flow. The first assumption reflects the fact that only spatially averaged parameters are contained in the conservation equations. Consequently, the distributions within the control volume can not be calculated. The second assumption means that, due to the predominance of the axial flow, the transverse momentum exchange can be crudely represented without introducing significant errors. Hence, the transverse momentum equation is usuallr much simpler than the axial momentum equation. A number of- computer codes have been developed which use the sub- channel analysis method. Among these are included COBRA IIIC [2], COBRA-IV [1] and WOSUB (10]. These codes treat most of the important phenomena in the same way and in each code a marching type solution method is utilized. (COBRA-IV also contains a modified I.C.E. method [19] for transient analysis.) The marching method begins the Figure 2.la: Coolant-Centered Subchannel Layout - I - Figure 2.1b: Rod-Centered Subchannel Layout .IMF- - nowalm- ----- Immmwmr -32- calculation of the flow parameters at the core inlet and moves upward, in a stepwise manner, simultaneously solving the conservation equations for all subchannels, at each axial level. Typically, more than one sweep through the core will be required to obtain a converged solution. Therefore, the marching method is basically an iterative technique. For steady-state, single-phase conditions, the subchannel codes can generally predict the correct flow distributions in rod-bundles [16]. However, for two-phase conditions or severe transients, the use of the subchannel codes may not be strictly valid. For example, comparisons of COBRA-IIIC with steady-state two-phase flow measurements have indicated that the correct flow and enthalpy distributions could not be calculated [10]. Also, if a strong perturbation causes large lateral flow, then the basic assumption in these codes is violated. Furthermore, if reverse flow should occur, then the marchiAg type solution method will fail un- less appropriately modified. Consequently, although useful for many rod-bundle problems, the subchannel codes are limited in their applica- tions. 2.2.2 Distributed Resistance Models In order to eliminate the assumptions and restrictions of the sub- channel methods, distributed resistance models have been developed. These models also referred to as porous body models, use orthogonal coordinates and geometrically similar control volumes. The name for this method is due to the fact that frictional resistances are distributed throughout each of the control volumes. Quasi-continuum governing equations are written for the conservation of mass, energy or momentum and no simpli- fications of the transverse momentum equations are made. Consequently, -33- no restrictions are placed on the flow conditions. However, as in the case of subchannel analysis, the details of the flow structure within a control volume cannot be determined. By employing Cartesian coordinates, the geometrical noding of the rod-bundle can be the same as for subchannel analysis (i.e., for square bundles). Hence, the information obtained by the distributed resistance method is at least as detailed as that found in subchannel methods. Of course, the governing equations in the distributed resistance methods are more general than those in subchannel analysis methods. However, the key to successful use of this method is the correct formulation of expressions for the transport processes in the control volume (i.e., heat transfer, friction, etc.). These processes can be described for most conditions, but completely general formulations are noc possible. However, these processes can usually be defined for many cases of practical interest. A number of computer codes use the distributed resistance method. Among these are COMMIX-1 [11], TRAC [121, and THERMIT [3]. Each of these codes use some form of the I.C.E. solution technique [6]. This technique coupled to the full three-dimensional representation of the distributed resistance method allows for the calculation of flow reversal, recirculating flow and even counter-current flow in two-phase conditions. With the ability to analyze these conditions, the above codes represent powerful tools for steady-state and transient thermal-hydraulic analysis of rod-bundles. 2.3 Two-Phase Flow Models Aside from the choice of modeling technique for the flow, the other -34- important feature to be defined is the two-phase flow model. A wide variety of possibilities exist for describing the two-phase flow. These range from describing the two-phase flow as a pseudo single-phase flow to a multi-component flow (e.g., liquid, vapor and droplets). The various possibilities are summarized in Table 2.3. Generally, as the two-phase flow model becomes more complex (i.e., more equations), more constitutive equations are required to represent the various interactions between the phases. The homogeneous equilibrium model (HEM) is the simplest of these models. It assumes that the vapor and liquid are in thermal equilibrium and that there is no relative velocity between the two phases. These assumptions are clearly limiting, but may be adequate for certain flow conditions. Extensions of this model to include relative velocity (slip) and thermal non-equilibrium (subcooled boiling) effects are possible using empirical models. The drift flux models, either the four or five equation models add some complexity to the two-phase flow description. By treating the vapor and liquid phases as separate streams still in thermal equilibrium, these methods allow for accurate velocity predictions. In the two-fluid model, separate conservation equations are written for the vapor and liquid phases. This model allows a very general des- cription of the two-phase flow. However, it also introduces a large number of constitutive equations. The most important relations are those which represent the transfer of mass, r, transfer of energy, Qi, and transfer of momentum, F, ,across liquid-vapor interfaces. The advan- tages of using this model is that physically based mechanistic models can be formulated for these terms which should be valid over a wide range -35- of conditions. An extension of the two-fluid model is the three-fluid model in which the three-fluid fields are the vapor, liquid and droplet fields. COBRA-TF [8 ] is an example of such a model. This formulation, while introducing more constitutive models, seems to contain the necessary capability to analyze complex flow situations such as the reflood stage of a LOCA. However, for all but reflood analysis, the two-fluid model is probably general enough to describe the important non-equilibrium effects. Con- sequently, THERMIT which uses the two-fluid model, is expected to provide a good description of most two-phase conditions. The complexity of the two-phase flow model is seen to depend on the assumptions concerning the non-equilibrium phenomena. Mixture models, that is either the homogeneous equilibrium or the drift flux models, contain one or more restrictions on either the thermal or mechanical non-equilibrium in the flow. Only when both phases are represented with separate conservation equations can all the non-equilibrium effects be modeled. 2.4 Description of THERMIT 2.4.1 Background It is instructive to review the key characterisitcs of THERMIT prior to the description of the modifications involved in the present work. The characteristics include the conservation equations, finite difference equations and constitutive models. 2.4.2 General Characteristics The thermal-hydraulic computer code, THERMLT, originally developed -36- at MIT under EPRI sponsorship, solves the three-dimensional, two-fluid equations describing the two-plir;e flow and heat transfer dynamics. This two-fluid model uses separate partial differential equations expressing conservation of mass, momentum, and energy for each individual fluid phase. By using this two-fluid model, thermal and mechanical non- equilibrium between the phases can exist, only requiring that mathematical expressions for the exchange of mass, momentum and energy be available. Such a formalism allows very general and physically reasonable modeling 'of relative motion of the phases and of thermal non-equilibrium. A second important feature of the THERMIT fluid dynamics is the three-dimensional representation of flow is x-y-z geometry. Previous codes (e.g., COBRA-IIIC) have used a subchannel model which assumes predominantly axial flow. The rectangular coordinate system in THERMIT is well-suited for either core-wide or subchannel analyses. THERMIT also offers the choice of either pressure of velocity boundary conditions at the top ana botuom of the core. This feature permits realistic modeling of the core boundary conditions and is important for reactor transient analysis. A third important feature of THERMIT is the heat transfer modeling. A radial heat conduction model (with gap conductance between the fuel pellet and cladding) is used with a continuous general boiiing curve describing heat transfer to the fluid both below and above the critical heat flux. The boiling curve is based on recommendations by Bjornard [20] and consists of five basic regimes: convection to liquid, nucleate boiling, transition region, stable film boiling, and convection to vapor. The heat flux is modeled as a heat transfer coefficient times a wall- fluid temperature difference in all regimes except in the transition -37- region, where the heat flux is computed for each phase. The final important feature of THERMIT is the numerical method used to solve the fluid dynamics equation. A semi-implicit technique is used which is a modified version of the I.C.E. method [5,6]. As such, the method has a stability restriction in the form of a maximum allowable time step: At < (Az/Vk min (2.1) where Az is the mesh and Vk is the larger of the phase velocities. How- ever, the method is not restricted by the direction or speed of the flow. Furthermore, convergence can be obtained at each time step if the time step is sufficiently small.. Consequently, this numerical method is ideally suited for severe transient analysis. Although coarse mesh sizes had originally been envisioned when using THERMIT, there is no intrinsic reason to prohibit the code's application to small mesh size problems; up to a point. From a numerical point of view, the solution method does not explicitly restrict the size of the mesh. However, due to stability considerations, a linear mesh size smaller than 0.2 mm may lead to numerical problems [21]. Since this limit is at least 30 times smaller than subchannel size, no instabilities would be expected for subchannel applications. Overall.then, it can be stated that THERMIT is a very powerful analytical tool. This code contains an advanced two-phase model and a fairly unrestricted solution.technique. Also since the code is theor- etically not restricted to large mesh applications, THERMIT would seem to be well-suited for subchannel applications. However, as will be dis- cussed in Chapter 3, the original version of THERMIT had certain geo- -38- metrical and physical characteristics which prevented accurate sub- channel analysis. Hence, the code needed to be modified to permit this type of application. 2.4.3 Two-Fluid Model Conservation Equations The governing equations of the two-fluid model in THERMIT, which are the mass, energy, and momentum conservation equations for each phase, can be derived from local, instantaneous conservation equations. The general procedure is to average the equations over time and then average them over an arbitrary volume. The result is a set of time and space averaged conservation equations which contain a number of integral terms. Examples of this type of derivation can be found in references 22 and 23. The THERMIT conservation equations are derived in Appendix A. This derivation begins by applying the appropriate time and space averaging operators to the local, instantaneous balance equations. The assumptions required to obtain the THERMIT equations are given and,by suitable re- arrangement, the appropriate two-fluid model equations are obtained. The major simplifying assumptions are: (1) that viscous stress and energy dissipation can be neglected, and (2) that the liquid and vapor pressures are assumed to be equal within any control volume. The assumption concerning the viscous and energy dissipation terms is appropriate due to the relatively small value of these terms. The assumption of uniform pressure is also appropriate provided the size of the volume is not too large. Following the derivation in Appendix A, one obtains the following set of equations: Conservation of Vapor Mass T (pv) + -(pV) = 1-tv Conservation of Liquid Mass S[(1-a) p9 ]atk Conservation - (p vev) + = Q +Q-io Conservat ion + V[(1-a)pv] = -F- W of Vapor Energy V(ap e V) + P 0-(V + P v v v v at Qtv of Liquid Energy [(l-a)pk] + V-[(1-a)p e V ] + P V.(i-a)o] aQQ - Pa = Q - Q - Qtk Conservation of Vapor Momentum V + + + + ap --- X+czpY *V ?V +ctVP =- F - F v at v v v v wv iv 4. 4+ + cpvg - Ftv Conservation of Liquid Momentum (+ p+ (1-a) VP = - -F + (1-a g + w -9 + (1-a)P gz -39- (2. 2a) (2.2b) (2.2c) (2.2d) (2.2e) (2.2f) -40- The notation for these equations can be found in the Nomenclature section. A few important characteristics of these equations should be discussed. First, it is seen that all the important transport mechanisms are included. In particular, the terms describing the turbulent transport effects are included in these equations. In the original version of THERMIT, these terms had been neglected. However, for subchannel applications, as well as large fluid plena applications, it is imperative that these terms be included. The turbulent transport terms are discussed in detail in section 3.3. A second feature of the equations concerns the representation of the interfacial heat and momentum exchange terms. As written, these terms include the effects of mass transfer between the phases. That is, the interfacial heat transfer term, Q., includes heat conduction between the phases as well as the heat transfer due to mass transfer (e.g., evap- oration). Similarily, the interfacial momentum exchange term, Fi, includes the momentum exchange due to mass transfer. In the original version of THERMIT both of these mass transfer effects had been neglected. The absence of the momentum exchange contribution is probably appropriate due to its relatively small value for most problems. However, the energy exchange contribution is comparAble with the other terms and, hence, has now been included in the present formulation. Further details on these models may be found in Chapter 4. 2.4.4 Finite Difference Equations The finite difference equations which approximate the above conser- vation equations, without the newly added terms, have been presented in -41- reference 3. The procedure for obtaining the difference equations is to approximate the temporal and spatial derivatives by difference operators. Since a semi-implicit differencing method is used, the temporal deriva- tives are replaced by a forward difference operator. The other terms are treated either implicitly or explicitly depending on the term. New time variables are represented with the superscript n+l while old time variables are superscripted with an n. Source terms are treated as implicitly.as possible, but do contain some variables evaluated at the old time. Consequently, source terms are superscripted with n+1/2 to indicate their semi-implicit formulation. The spatial discretization of the equations requires a three-dimen- sional grid to be overlayed on the geometrical configuration under con- sideration. Rice this grid has been defined the locations of the variables are determined. The convention for associating the variables with a particular mesh cell is illustrated in Figure 2.2. All unknowns except the velocities are associated with cell centers. The velocity components normal to that face are defined. On cell faces, subscripts for the cell- centered quantities are i, j, and k, while on cell faces half incre- ments are used (e.g., 1+1/2, j, k). In order to simplify the following equations, only the half integral subscripts will be retained (e.g., Pi+1 /2 refers to P i+1/2,j,k)' With this background information, the finite difference equations will now be given. In the mass and energy equations control volume flux balances are used to approximate the divergence terms. With this approxi- mation the equations are: tV4 k-+ v x Vi. l i,j,k fv / / / / / / / Vi+ Cell Centered Quantities: P, Tv' T' eV . ' PV' P ' (, Figure 2.2: Typical Fluid Mesh Cell Showing Locations of Variables and Subscripting Conventions / x i+ -43- Vapor Mass (ap V)n+ l apV)n 1n x l v v -+ LA/ PQv(v n+1 1+1/2At - [A(ap)"(V) n+i-1/2 + [A(apv) n n+1 j+1/2 v v j-1/2 + v k+1/2 - [A(ap )nz(Vz n+1 v v k-1/2j Sn+l/2 n+1/2 tv Liquid Mass ( (1-cx) p dn+1 n-(l~cPz) At +1 -[A( (,a) P.) n yn+1 -/2+ - [A((1--)p )n)n+1 1/2 + - [A((1-a)pz) n 1n+1] _ _]n+1/2 k-1/2 [A((1-a)p(V)n+ i+/ [A((1-a) p nY) n+1j+1/2 [A((1-a) nVz) n n + 1/ n+1/2 wt 9 Vapor Energy (czp e )n+1 v v n S(czv e +1 [ipn + +/ n1+1/2At - [Pn + (pv -1/2]v[Acn ( +1 i-1/2 + [pn + (pve)?+!2[An( n+ljl,2 - [ + J-1/2 (2.3c) continued (2.3a) (2.3b) -44- +Pfl + (pev)n I[Acf(VZ)fl+l]I v v +l/2 v k+1/2 (2.3c) concluded - [Pn + (pe1) - 2 ][Aa?(Vz) n+11 ] + Pn a+1 - an n+1/2 n+i/2 _ n+1/2 At WvQtv Liquid Energy ((i-a)pPe2 )n+1_ -((l-a)pke )n (2.3d) At + Pn + (P e ) +1/2][A(l-a)n(Vt) n+1, 1 - [Pn + (p e ) -n/2][A(la)n xn+l - + [Pn + (P e ) +11 2 ][A(l-)n l n+V) +l / - [Pn + (pzet)?-j/2I[A(l-c)n n+1j-/]2 j-l/2 + [Pn + (P e )k+n 2][A(1-)n zn+lk++(pe)k1/22P, k+l/2 - Pn + (P e ) -i2][A(ct)n (Vn+ 1 k-l/2 n (+1 - an n+1/2 n+1/2 n+1/2 - P At )%9Jz N t2. For the momentum equations, the equation for a particular direction is differenced between the centers of the two appropriate cells. Conse- quently, tbe mesh used in the momentum equations is different than that used in the mass and energy equations. Since there are a total of six momentum equations all having the same form, only the z-direction vapor equation will be given. The other equations are found by permutation. -45- The vapor z-direction momentum equation is given by (zyn+l zynz[(V,) l - (Vv) ( Pj] 1 2 x ( sv) ap)n v v k+1/2+ a) n xV)Ax v(k+1/2 At + v k+1/2 v k+1/2 Axk+1/2 + (VZV))A zV l v k+1/2 Ay k+1/2 + (V )+/2 z k+l/2 Sk+ k + zk+1/2 Az n+1 /2 n al k~ -Pk~l (Fz ) n+l/2 (F z )f+l/2 k+1 2 [AZkk+1/2 l= - (w k+l/2 - iv k+1/2 -Qap g) n (F n (2.4) v k+1/2 tv k+1/2 .A few important features of these equations need to be highlighted. First of all, it is seen that values are required for the unknowns at locations other than those defined by the noding convention. For example in the mass equation, the quantity (pv)1+ 1 /2 is required. For all such terms in the mass and energy equations the donor cell logic is used. Mathematically this can be expressed as Ci+1 i+1/2 <0 Ci+/ = (2.5) C i+12 - C if Vi+1/2 >e0 where C is the quantity of interest (i.e., a, P, Pv' p etc.). In the momentum equations, no such general rule exists for specifying variables at locations other than the noding convention. Instead, each required term is specified separately. The expression for (cPv)k+1/2 is given by (apv k+1/2 " Xk+1/2 pvk+1/2 (2.6) -46- where Ck+l AZk+ + (AZk 'k+1/2 =AZk+ + AZk(2.7) and (pv k+l AZk+l + v AZk (2.8) v k+1/2 - AZk+1 + AZk Every velocity except (V Z)k+1/2 needs to be difined since they are not at the- aprior defined location. These velocities are defined by: k+1/2 j-1/2 + v)j+1/2 v+(v j-1/2, k-l + (V j+/2, k+ (2.9) (1?) =!1I(VX) +(x) (X vk+1/2 4 v 1-1/2 + (V 1+1/2 v -/2, k+ + (V')5+ 1/2 k+ 1 (2.10) Finally for the convective operators, which use the donor cell logic, the following expressions are used: (V)i+1 , k+11/ - () (Vkx ( A V z A X -1 / 2 v k +/ 2 < 0x v_ . (V)k+1/2 - (-Vv)i 1 , k+1/2 if (it) 0 AX-1/2 vk+1/2 (2.lla) Vv ( +1V k+1/2 - k+1/2 A V z AY J+1/2y v_ AY k+1/2 (VZ)k+1/2 j- k+l1/2 AY j-1/2 (VZ)k+3/2 -(Vv)k+l/2(v z AZ k+1ZAv AZ k+kl2(zvz . (VZ)k+1/29V;)k-1/2 AZ k where the mesh (AX) i+1/2 i f (Vbk) ,2 <0v k+1/2 ' (2.lb) if (VY) 1 2 0 if (V z) 1 ,2 <0v k+1/2 (2.lIc) if (VZ) / 0 v k+lI/ 2 .0" spacing are given by (AXi+1 1+ A) 2 (AYJ+1 +AY.) - j+12 (2.12a) (2. 12b) The second important feature of the equations is the definition of the transverse flow areas (i.e., AX and AY). As the rows of rods are transversed, the cross-sectional area normal to the x (or y) direction changes with x (see Figure 2.3). Since the momentum equation control volumes do not coincide with the mass and energy control volumes' it is necessary to carefully define these areas. To be consistent, the cell averaged transverse flow areas must be used [3]. This requirement is the origin of the concept of a distributed resistance approach in which -48- Momentum Equation Control Volume - 0 L (D 0..00 U .I- lljlllllll IU00_________ Mass or Energy Equation Control Volume x Figure 2.3: Typical Rod Arrangement in Transverse Plane 0 0 0 0 it I -49- the structure and associated resistances are averaged over the control volume. By using volume averaged flow areas, the transverse velocity and flow areas are consistent so that continuity of mass and energy can be achieved. 2.4.5 Constitutive Equations The two-fluid formulation of the conservation equations introduces terms which represent the transfer of mass, energy or momentum in a given control volume. These transport processes occur at one of the four types of interfaces found in a control volume. These interfaces include: (1) Wall-Liquid Interfaces within cell volume (2) Wall-Vapor Interfaces within cell voldme (3) Liquid-Vapor Interfaces within cell volume (4) Inter-Cell Interfaces at cell boundary Table 2.4 summarizes the transport mechanisms which occur at each interface. The wall friction and wall heat transfer terms are common to all thermal- hydraulic codes. However, f9 r the two-fluid model, the total friction or heat transfer must be apportioned into liquid and vapor components. A unique feature of the two-fluid equations is that the transport of mass, energy and momentum; across liquid-vapor interfaces must be modeled explicitly. These interfacial exchange terms while presenting complex interactions, do allow for general modeling of phasic non-equil- ibrium. Across the interchannel interfaces turbulent eddy transport leads to the transfer of mass, energyand momentum. The terms which represent these transport mechanism are referred to as the turbulent mixing terms. These terms account for the coolant-coolant interactions which occur due . I -50- TABLE 2.4 Summary of Transport Processes Wall to Coolant F Fwy QuV Liquid to Vapor F1 Qi - Wall Frictional Force on the Liquid - Wall Frictional Force on the Vapor Wall Heat Transfer to the Liquid - Wall Heat Transfer to the Vapor - Interfacial Mass Transfer Rate - Interfacial Momentum Exchange Rate - Interfacial Heat Exchange Rate Inter-Cell w tv WtR, Qtv Qt F tv Ft P. - Turbulent - Turbulent - Turbulent - Turbulent - Turbulent - Turbulent Vapor Mass Exchange Rate Liquid Mass Exchange Rate Vapor Energy Exchange Rate Liquid Energy Exchange Rate Vapor Momentui Exchange Rate Liquid Momentum Exchange Rate -51- to property gradients from one channel to the next. All of the above terms need to be specified by correlations in order to specify the variable of the two-fluid model. The specific correlations are discussed along with their assessment in Chapter 4, 5, and 6. In general, mixture model correlations are apportioned to determine the wall- coolant transport processes while mechanistic models are used to describe the interfacial and turbulent transport processes. -52- 3.0 DEVELOPMENT OF THERMIT SUBCHANNEL ANALYSIS CAPABILITY 3.1 Introduction While containing many capabilities for thermal-hydraulic analysis, THERMIT, as originally written, was not acceptable for analyzing sub- channels for two reasons. The first was the limitation of being able to model only one fuel rod per unit cell (i.e., per coolant channel). This restriction prevented the analysis of traditional coolant centered sub- channels in which up to four fuel rods per channel need to be modeled. Since many experimental measurements are based on this coolant centered geometry, for validation purposes, the code had to be able to analyze coolant centered subchannels. The second deficiency in THERMIT was the lack of description of turbulent mixing or for that matter any coolant-coolant interactions at channel boundaries. For large control volumes, the omission of tur- bulent mixing is probably justified but for subchannel analysis these effects are very significant. Hence, in its original form, the use of THERMIT for subchannel analysis seemed questionable at best. In view of the above mentioned deficiencies, a developmental effort has been undertaken to improve the capabilities of THERMIT. This effort has required a reformulation of the capability to describe the inter- actions of the coolant and the fuel rods. Additionally, a turbulent mixing model has been added to THERMIT. This model has been formulated to be applicable for both single and two-phase coolant conditions. With these modifications, THERMIT now contains consistent thermal-hydraulic models capable of traditional coolant-centered subchannel analysis. A detailed description of the modifications is given in this chapter. -53- 3.2 Geometrical Modeling Capabity The original version of THERMIT permitted the modeling of only one fuel rod per coolant channel which is adequate for coarse mesh (core) analysis, where only the behavior of an average fuel rod can be determined. However, if subchannel size control volumes are used, then only rod cen- tered subchannels can be exactly analyzed (see Figure 2.1). If coolant centered subchannels are to be analyzed then the fuel rod sections within a channel would have had to be lumped together. Clearly, this lumping causes the loss of all information about the actual clad temperature distribution. Since one of the reasons for performing subchannel analyses is to determine the clad temperature distribution, the restriction of one rod per channel is not compatible with subchannel analysis. Furthermore, a second aspect of the geometry is related to the valida- tion of the code. This validation relies on comparing the code pre- dictions to experimental measurements. Many measurements in rod arrays have been made based on a coolant centered subchannel. Consequently, realistic calculations and comparisons are feasible only if the same geometry is used. In view of these considerations, THERMIT has been modified so that coolant centered subchannels could be modeled. The first step in this reformulation has been the modification and expansion of the coupling between the coolant and fuel. This coupling occuts through the heat flux which can be written a . n+l " l n n+l n3 q =1H (T -Tf ) (3.1) wf The heat flux couples the temperature calculations in the fuel to the thermal-hydraulic calculations in the coolant. The significant feature -54- of this coupling is that it is implicit in time. In order to have this implicit coupling, the special algorithm described in Appendix E of reference 3 is needed. This algorithm, outlined in Table 3.1, solves the fuel rod conduction equation in a two-step procedure. By doing this, the wall temperature is found iteratively, thus preserving the implicit coupling between the wall and the fluid temperatures. This coupling has been maintained and expanded so that up to four fuel rods can contribute to the power input of a subchannel. The power released by the ith rod to the jth subchannel is given by q n+1 = H.n(Twn+1 -T n+1)P .AZ (3.2) j w. f hi 2. j Where P . is the heated perimeter of the ith rod which faces the adjacenthi subchannel. The power input to the jth subchannel is then given by 4 n+l 4H n+ln+l E q.nH = H(T n+-T f +) (AZ(3.3) = i=l $ wi f PhiA With this modification, a given subchannel can be coupled to as much as four fuel rods, consistent with coolant centered subchannel analysis. However, a complication is introduced with this formulation. Namely, four clad temperatures are required for each fuel rod. (Alternatively an average clad temperature for each rod could be defined, but would not be consistent with the implicit coolant-to-fuel coupling.). This requirement, while increasing the complexity of the heat transfer calculations, allows for detailed fuel rod modeling. Since accurate fuel rod temperatures are of interest, the increased complexity caused by adding this capability is certainly welcome. This capability has resulted in a more general thermal modeling of the fuel rods. In particular, a given fuel rod does not necessarily -55- TABLE 3.1 Implicit Heat Transfer Algorithm 1. Calculate Hn using previous time step wall and fluid conditions. 2. Set up fuel rod conduction equation using the boundary condition q" = Hn n+1 _Tn+1 W f n+1 n at this stage the assumption T = Tn is made. f f 3. Forward Elimination of the rod conduction problem yields both an initial guess for new wall temperature, n+1,O0 n+1 n+1T and3T /13T . w w f 4. Solve the fluid dynamics equations using q" -. n n+1(0) _Tn+1) + Hn n+1/3Tn+1 n+1Tn w f w f f f n+1 n+15. Once T is found, T is calculated usingf w Tn+1 -=Tn+1,(0) n+1 n+1 - TnW w w f f f 6. Complete the backward substitution step of the rod conduction equation. -56- have to be modeled as a single rod. That is, the rod may be divided into four quarters with each quarter being analyzed separately. The only restriction of this method is that each rod section must be adjacent to only one coolant channel. Hence, as illustrated in Figure 3..1 a rod may be divided into four sections, two sections or one section. In each section, a complete heat transfer calculation is performed so that the temperature distribution throughout the section is calculated. Thus, the clad temperature for each rod section are determined as required by the expanded fuel-to-coolant coupling. A disadvantage of this approach is that for any given rod modeled as four sections, there will be four centerline temperatures calculated which are not necessarily equal. This is not always accurate due to azimuthal heat conduction effects which are neglected here. For all cases of practical importance that have been run, negligible differences in the centerline temperature were calculated. Furthermore, uncertain- ties in the physical properties and voiding of the fuel near its center outweigh this numerical approximation. Another minor disadvantage is that the computational time will be increased, but this increase should not be excessive. Therefore, on the whole, the fuel pin expanded model- ing together with the coolant-centered subchannel capability provide THERMIT with the geometrical flexibility required for subchannel analysis. 3.3 Two-Phase Turbulent Mixing Model 3.3.1 Background One of the most important phenomenon that must be accounted for in subchannel analyses is the exchange or mixing of coolant between adjacent subchannels. This mixing is important as it leads to the transfer of -57- 112 34 9 110 11 12 C2 5 6 78 13 14 15 16 314 Each Rod Modeled In 4 Sections Each Rod Modeled In 2 Sections Figure 3.1: Illustration of Fuel Rod Modeling I -58- mass, energy and momentum between adjacent subchannels. As discussed by Rogers and Todreas [24], who have reviewed this subject for the case of single-phase flows, this mixing can be either forced or natural. Forced mixing is caused by mechanical protuberances, such as grids or wire wrap, which either randomly break up the flow or actually divert the flow in a preferred direction. Natural mixing, which occurs in the absence of such protuberances, consists of turbulent mixing and diversion cross-flow. The main distinction between these two types of natural mixing is whether the mixing occurs with or without pressure gradients. Turbuleht mixing results from the natural eddy transport between subchannels, while diversion cross-flow is caused by radial pressure gradients. Even though these various types of mixing have been identified for single-phase flows, they will also exist in two-phase flows. However, an additional mixing phenomenon has been postulated to occur in two- phase flows. This phenomenon, known as vapor diffusion, has been pos- tulated in order to describe the experimental measurements which could not be explained with single-phase concepts of turbulent mixing [25]. Specifically, the void fraction profile in a rod bundle geometry, referred to as the fully developed distribution, is such that the more open areas have the larger vapor (void) concentrations. In other words, the vapor tends to diffuse to unobstructed regions. This observation cannot be predicted using turbulent mixing alone. Hence, vapor diffusion must also be included as a separate type of natural mixing in two-phase flows. Another important difference between single and two-phase flows is the mechanism for eddy transport. In single-phase flow the conventional approach is to assume equal mass exchange between two cells such that no net mass transfer occurs due to turbulent mixing. However, in two- -59- phase systems, the equal mass exchange model must be replaced by an equal volume exchange model in order to explain energy mixing [263. in this latter model a volume of vapor in one channel exchanges with an equal volume of liquid in the adjacent channel. In this manner a net transfer of mass occurs. However, this mass transfer is reasonable as evidenced by experimental findings [26]. Furthermore, if one considers energy exchange in saturated boiling, then with an equal mass exchange model no energy can be transferred. On the other hand, energy is clearly transferred in the equal volume model. It should also be noted that in single-phase flows the two exchange models are essentially equivalent. Consequently, in view of its superior physical interpretatien of two- phase flow, the equal volume exchange model is preferred. If all these types of mixing are present, then they must be accounted for in the analytical model. In the original version of THERMIT, only the diversion cross-flow type of mixing has been explicitly considered. The omission of the forced mixing is justified since the use of wire-wrap of flow diverters is not common in light water reactor rod bundles. However, the absence of turbulent mixing would only be justified if the scale of the mesh size is much larger than the scale of the turbulent eddy transport in the flow. Practically, this means that the dimensions of the computational cell must be larger than the Prandtl mixing length (taken here as the scale of the turbulence). This condition is met when subassembly size cells are used. However, for subchannel applications the mixing length is on the order of the subchannel size, therefore, the effects of turbulence must be included. -60- 3.3.2 Model Formulation 3.3.2.1 Background In view of the need to include the effects of turbulent mixing, a model suitable for use in THERMIT has been developed. Before discussing the machematical formulation of this model, its important physical features will be discussed. The first important feature is the two-phase flow regime dependence of the model [27]. By considering the mixing rate as a function of the flow quality, the general characteristics of the model are qualitatively described as follows. . For single-phase liquid or single-phase vapor, the model naturally defaults to a single-phase tur- bulent mixing model appropriate for the phase present. For two-phase conditions it is found that the mixing rate is enhanced above the single- phase rate. The peak ixing rate depends on the flow regime. Near the slug-annular transition point the mixing is found to be greatest and at this point the model reaches its maximum. The approach to this maximum from the two single-phase limits is approximated with simple functions. These functions, while not exactly corresponding to measured data, do represent the observed trends in the two-phase mixing rate, as will be discussed in Chapter 5. Aside from the flow regime dependence of the model, a second impor- tant feature is that,for tUo-phase conditions the effects of both tur- bulent mixing and vapor diffusion are included. Although these two phenomena are different, they do share similar characteristics. Both phenomena result in the transfer of mass, energy and momentum between adjacent subchannels. Furthermore, both phenomena are postulated to occur in the absence of pressure gradients. In view of these similari- -61- ties the two phenomena are combined into a net two-phase mixing model. The mathematical formulation of the two-phase mixing model has been adapted from previous work for mixture models (28,29] in order to fit into the framework of THERMIT's two-fluid model. The previous work has dealt with simpler two-phase flow models such as homogeneous or drift- flux. The extension of the semi-mechanistic models which attempt to represent the important physical processes to the two-fluid model is not unambiguous. Many possibilities exist for proportioning a given model when going from a homogeneous model to a two-fluid model. However, the one constraint which must be observed is that when the liquid and vapor portions are added together the sum is the homogeneous model result. The adaptation of these previous models is, therefore, to some extent arbitrary and the choice given here is justified by the validation with experiments discussed in Chapter 5. Another consideration of the mathematical formulation is that the equations must not introduce any numerical instabilities. Both the numerical and physical description of a particular phenomenon are impor- tant to the overall structure of THERMIT. Since THERMIT is not fully implicit, a number of temporally explicit terms are included in the equations. These explicit terms require stability restrictions such that numerical errors do not grow with time. The inclusion of a new phen- omenon, such as two-phase mixing, may introduce limits which did not previously exist. Hence, the numerical representation of any phenomenon must be done with some care. The two-phase mixing model can be discussed in terms of its physical as well as its numerical nature. In order to simplify this description, the physical nature of the model will be given first. The salient -62- numerical aspects of the model will be given next. 3.3.2.2 Analytical Formulation and Discussion The terms which represent the effects of turbulent mixing have been identified in Section 2.4. These terms, which represent the exchange of mass, energy and momentum due to turbulent interactions, can be written as follows: S.. W =V wi=" =- ' (").. (3.4) tv tv Ai tv ij S.. Wz VWw" A = E 1 (W"z) (3.5) tz tz A tz ij S.. Q =7 - q" = ZU(q" ). (3.6)tv cv A tv ij S.. Q =V -q"t = E I (q" ).. (3.7) tk tz A tzi 1 S.. F =V- T = E (T ) (3.8) tv tv A. tv ij S F =V' T = E (T )(3.9) tt tv A. tzi 1 where S is the gap between rods. The divergence operator has also been approximated by its control volume equivalent. In each of these equations, there are still the terms which represent the flux of the specific property involved. In the mass equations, W"tv and W't are the vapor and liquid mass fluxes due to mixing. The energy equation terms q" and q" are the vapor and liquid heat fluxes due to tv tz mixing. Finally, in the momentum equations the terms T tv and rttZ are the shear stresses due to mixing. It should be noted that only transverse mixing is considered which is appropriate due to the larger interaction area in the transverse direction. The set of equations which describe these terms are similar in form and, hence, only one of these terms will -63- be discussed in detail. If one considers single-phase flow, then the turbulent shear stress term, also referred to as the Reynolds stress, is usually written as T = Ep(3.10) where C is the eddy diffusivity. This term can be approximated as follows: ,E (G. -G.) ( 3)(3.11) where is the effective mixing length. This form is convenient to use provided the term E/ can be determined. This term, which has the dimensions of a velocity and is sometimes referred to as the turbulent velocity, can be related to measured mixing through the following equation E/k = W'/Pj (3.12) The mixing rate, W', has been measured by a number of authors [27,30], and is usually expressed as a function of the Reynolds number. Once E/ is determined, the Reynolds stress is easily calculated. The extension of the above equations to two-phase flow requires the addition of two physical effects. First, the vapor diffusion phenomena must be added to the model. Following the work of Lahey [25], the two- phase Reynolds stress is written as TTP = E/2(G - G. - (G. - G )F) (3.13) tP u i j i FDu where the subscript, FD denotes the fully-developed distribution. This form of the shear stress term accounts for both turbulent mixing and -64- vapor diffusion. The justification for this form is based on the experi- mental work of Gonzalez-Santalo and Griffith [26]. These authors have shown that the net two-phase mixing is proportional to the nonequilibrium void fraction gradient. Specifically, they have been able to correlate their vapor mixing rate data using the following equation W' = K p [(a. - a.) - (a. - a.) ] (3.14) v i j . j3FD where K is an empirically determined constant. By analogy, this form can be applied to the other mixing terms so that the shear stress can be written as in Equation 3.13 [25]. The second physical effect which must be included is the dependence of the mixing rate on the flow regime. This effect is included by writing the two-phase turbulent velocity, (c/Z)TP, as a function of the single- phase value, (e/9)5 . Mathematically, following the work of Beus [27], this can be expressed as (e/)TP = (6/SP e (3.15) where 6 is a "two-phase multiplier" which depends on the quality. As indicated above, the mixing rate (and hence E/) reaches a maximum at the slug-annular transition point. The criterion for this transition point is the Wallis model [31] which is in terms of the superficial velocities: J = 0.4 + 0.6 (3.16) where * 1/2 - -1/2.(3,17) v v v sD -65- and j 1/2 (p _p )g D]-1/2jz ~ z JvP (9 jv = XG/pv j = (1-X)G/pZ (3.18) (3.19) (3.20) By rearranging Equation 3.16, an expression for the quality at the transi- tion point can be obtained: XM = [0.4(p9 (P2 - Pv D) 1 1 2/G+ 0.6]/[(p%/pyyl/ 2 +0.6] (3.21) Again following Beus, the function for 0 is assumed to increase linearly between X = 0 and X = XM. For qualities greater than XM, 6 is assumed to decrease hyperbollically. At X = XM,0e = 0M i.e., 8 reaches a maximum. These conditions are expressed as follows: 0 = 1 + (M- 1) t/X1 ( Xo/XM0 = 1 + (0M - 1) if X < if X> xm x - 0.57 Re0.0417 (3.24) as correlated by Beus. If a value for 6M can be prescribed, then the function describing 6 is complete. and (3.22) (3.23) -66- The only remaining unknown in model for (s/2)TP is the value for (E/Z)P. As given in Equation 3.12, (E/t) S can be expressed as a function fo the mixing rate. The correlation for the single-phase mixing rate used in this model is that of Rogers and Rosehart [30]. By making appropriate substitutions, their correlation can be expressed as a single- phase turbulent velocity which is given by (E/) = 1 A Re~0'1[1 + (2)1.5 DiGi (3.25) SP 2 Di DFS P where .. -1.46 X = 0.0058 I -(3.26) D'.S Both the Reynolds number and the density are based on the two-phase mixture and DFS is the fuel rod diameter. With this correlation, the description of the two-phase turbulent velocity is complete. It should be noted that this velocity is assumed ,to be the same for both vapor and liquid phases. Clearly, there may be some differences in this velocity for each phase, but within the scope of the overall model these would have small, if any, impact. - Returning to the mixing terms, by analogy these terms can be written as (w") = IR. [ap a. - [&K)a - QMQ v FD1 (3.27) tv ij v ijiF (W't)i. = E/ [((-a)9P) - ((l-cP)t -(1-0)P) - ((l-a)pz)IF] (3.28(3.28) -67- (q").. = E/Z [(ape). - (ap e ). - [(ape). tv i3 V v i vv jV v - (apve v j FD] (3.29) (q" = E/Z [((l-a)p e) - ((l-c)pe2,) - [((l-a)pZe ). - ((1-a)p e ) IFD (3.30) (T v).. = 61/, [G .i - G .j - [G .i - G .j] F (3.31) tv 1 vi 3 vijFD F (Ttz. =ij E/Z [G U - G Zj - [Gi - G Zi]FD] (3.32) These may be written more concisely as = [(P) - (PW) - [P) - (P)]FD](3.33) where $ is the generalized mixing flux term and = 1 in the mass equations - e in the energy equations - V in the momentum equations p = apv in the vapor equations = (l-a)p in the liquid equations. The only terms ,yet to be discussed are the fully-developed dis- tributions. The assumption used by Lahey [25,32], is also used here as a basis for developing expressions for these terms. This assumption states that the fully-developed void fraction distribution is propor- tioned to the fully-developed mass velocity distribution. Mathematically this is expressed as follows: -68- (cc. - a) =K(G. - G.) (3.34) 1 jEFD ' jEFD Physically, this equation reflects the observed trend that the void fraction is higher in the channels with higher velocity. Recently, Drew and Lahey [32] have analytically derived this expression. While the fully-developed mass velocity distribution is also not known, it is assumed that this distribution is proportional to the calculated mass velocity distribution. By using these assumptions and extending the assumption concerning the fully-developed void fraction distribution to the other fully- developed distributions, the following set of equations are obtained. (G. - G.) ((a - (a ) KM (G +YG )(apvi - (ap)) (3.35) (((l-a)p ) - ((c-t)ZP)0)FD = (G. - G.) = -KM (G) G) + (p ))(3.36) M(C. + GC ) (c, i tP. ((ctp e ). - (ctp e ).) vv i vvjFD (G. - G.) = KM (G.+ C.) ((apvev + (apvev)) (3.37) 1 3 (((l-a)p- e ) - ((l-a)pke ) )FD = (G. - G.) =-KML + ((ap e) + (ape).) (3.38) - G ( +C)J j -69- (G. - G.) ((G).-(GK)) = M- (Gv. + G . (3.39) v i vjFD i(G.+G.) vi vj 1 J ((Gz)O - )(G0)FD = (G. - G.) KM (Gi + G) PiV ) + (ap ZV) )(3.40) With these equations the description of the two-phase mixing model is complete except for the specification of eM and KM. These two parameters are treated as constants and can be estimated from experimental measurements. For example, measured two-phase mixing rates show that 0M can vary from 2 to 10 for typical BWR conditions [27]. The value for KM should be approximately equal to unity [291. Using typical values, comparisons have been made with subchannel data and recommendations for these parameters have been made as discussed in Chapter 5. 3.3.2.3 Numerical Scheme As a final note about the two-phase mixing model, the numerical aspects of this model will be discussed. Due to the structure of the solution method in THERMIT, implicit coupling between adjacent channels is only possible through the pressure. This restriction prevents the implicit treatment of either the void fraction or internal energy for adjacent channel terms in the two-phase mixing model. Consequently, a fully implicit formulation of the two-phase mixing model is not possible with the current solution technique. In view of this restriction, two different numerical formulations of the two-phase mixing terms have been investigated. The first attempt -70- has been a fully explicit formulation which may be written as $n4= (/Z)n[(p)_ (p)? GP) -n (3.41) Si - j) FD (.1 However, when this formulation was used, the code developed numerical instabilities. Apparently, the inclusion of the mixing effects intro- duces an additional stability limit. Hence, this formulation was un- acceptable. The above explicit formulation has been modified to be at least partially implicit. This modification consists of treating the (4p). term in a fully implicit manner. Rewriting the generalized mixing term using this more implicit formulation one obtains: n+1/2= (1) n n+l nn $ (/t [$p. - ($p). - [(#p). - ($p).] ] (3.42)iJi jW) FD It should be noted that the fully-developed distribution term is still fully explicit. This formulation is more implicit in time, but is not exactly conservative. However, this should not be a problem since long time step sizes cannot be used due to the convective stability limitation (equation 2.1). Hence, this formulation represents the most acceptable combination of implicit and explicit terms and is the current choice for the two-phase mixing model. -71- 4.0 DEVELOPMENT AND ASSESSMENT OF THE LIQUID-VAPOR INTERFACIAL EXCHANGE MODELS 4.1 Introduction While the two-fluid model and modified ICE solution technique contained in THERMIT are very flexible and well-suited for analyzing transient two-phase flow, the accuracy of THERMIT is strongly dependent on the choice of constitutive models. As discussed in Section 2.4.5, these models are required to represent the various transport processes which occur in the two-phase flow. Careful definition of these models is essential for the accurate prediction of the complicated non-equili- brium effects which occur in two-phase flows. The most important non-equilibrium phenomena are subcooled boiling, vapor superheating and relative motion of the two phases (i.e., vapor slip). Both subcooled boiling and vapor slip are important for operational conditions, while vapor superheating becomes important only after the critical heat flux (CHF) has been exceeded. In addition to modeling these non-equilibrium phenomena, accurate constitutive models are required to predict the correct two-phase flow distribution for subchannel applications. These models account for the transport of mass, energy, and momentum due to turbulent mixing and vapor diffusion. Finally, for accurate wall temperature prediction, the wall-to- coolant heat transfer model must also be carefully developed. This model needs to account for the heat transfer mechanisms in both the pre-CHF -72- and post-CHF regimes. The program for developing and validating the models in THERMIT has been undertaken with two goals: (1) to define and develop necessary models and, (2) to validate the predictive capabilities. The emphasis of this effort has been on evaluating the code for subchannel applications. In particular, the liquid-vapor interfacial exchange models, the two-phase mixing model and the heat transfer model have all been carefully reviewed, modified as necessary and evaluated. The overall predictive capabilities of THERMIT have been judged primarily based on comparisons with experimental measurements. 4.2 Assessment Strategy In order to meet the goals of this program, an orderly progression of tests and ccmparisons has been performed. These include comparisons with both one-dimensional and three-dimensional experimental data. The order in which these comparisons have been made is structured so that iidividual models could be assessed and validated in a logical manner. Ideally, this procedure consists of selecting a set of eperimental data which can be used to evaluate a specific model independent of the other models. Then, once a model has been judged appropriate, it can be used with some confidence in the effort to validate the other models. In this manner the data base for the code is built up systematically. -73- Of course, it is not always possible to evaluate each model separately. For these cases, engineering judgement is required to interpret the results of the comparisons with experimental data. Nevertheless, the present evaluation strategy represents a viable method for validating THERMIT. Another consideration in the effort is that experimental measurements suitable for model evaluation need to be relatively simple in order to successfully validate individual models. Consequently, the measurements used in this study are straightforward and well-documented. Simple steady-state, one-dimensional measurements as well as steady-state and transient three-dimensional data have been used. The experimental measurements and models which have been assessed are summarized in Table 4.1. It is seen that, typically, more than one model is evaluated using a given set of data. This is possible, though, due to the logical order in which the comparisons have been made. In this systematic procedure, the initial evaluations have been performed using steady-state, one-dimensional void fraction data. The void fraction measurements of Maurer [33], Christensen [34], and Marchaterre [35] have been used in this study. These measurements are utilized to evaluate the interfacial mass exchange model, F, and the interfacial momentum exchange rate, F.. Steady-state, three-dimensional measurements have also been used in this study. These measurements include both mass velocity and quality distributions in subchannel geometry. Both the nine rod G.E. bundle data [36] and the sixteen rod Ispra bundle data [37,38] have been used. These measurements are useful for validating the two-phase mixing model. -74- TABLE 4.1 Summary of Assessment Program Constitutive Model Measurement F Q Turbulent Mixing Model Steady-State 1 - D Void Fraction X X Steady-State 1-D Wall Temperature X X Steady-State 3-D Mass Velocity and Quality X Transient 3-D LCS X 9 -75- The final comparisons have been made with steady-state, one-dimension- al wall temperature data [391 and transient three-dimensional CHF data [40]. These measurements permit both the transient capabilities as well as the heat transfer model to be assessed. The results of the comparisons will be discussed here in terms of the individual models. That is, each model, with the exception of the two-phase mixing model already described in Section 3.3, will first be described in detail and then the results of the assessment effort will be given. The interfacial exchange models are discussed in this chapter, the two phase-mixing results are presented in Chapter 5 and the wall-to- coolant heat transfer model assessment is discussed in Chapter 6. 4.3 Interfacial Mass Exchange 4.3.1 Background In the two-fluid model, the exchange of mass across liquid-vapor interfaces must be explicitly modeled. In reactor applications, this exchange usually takes the form of vapor generation so that the mass exchange model is also referred to as the vapor generation model. Physi- cally, this exchange of mass is strongly dependent on the flow conditions. The mass exchange actually occurs on mi:roscopic scales with vapor being produced at interfaces which are in constant motion. Hence, to describe this phenomenon on the microscopic scale would be virtually imp~ssible. Fortunately, the mass exchange needed in THERMIT is the net exchange which occurs in a given control volume. Consequently, only the integral of all the microscopic effects is considered in the formulation. However, even the integral or macroscopic vapor generation rate is difficult to define due to the various types of vaporization which can -76- occur. For BWR conditions at least three different vaporization regimes can be identified as illustrated in Fig. 4.1. The first is termed subcooled boiling due to the fact that vapor is generated even though the bulk liquid is subcooled. For this vaporization mechanism, vapor bubbles are formed at nucleation centers on the heater surface. These bubbles grow and detach when the bulk liquid temperature is above a certain value (referred to as the bubble departure temperature). Of course, since the liquid is subcooled, condensation of the vapor in the bulk fluid may also occur. Consequently, for subcooled boiling condi- tions, both the vapor generation on the heater walls as well as the vapor condensation need to be modeled. The second type of vapor generation is referred to as saturated flow boiling. As its name implies, this type of vaporization occurs when the bulk liquid is at saturated conditions. For these conditions, a liquid film is assumed to coat the heater surface. Heat is transferred directly through the liquid film so that vaporization occurs at the liquid-vapor interface. Since both phases are assumed to be at saturation, for steady-state conditions, all of the wall heat flux produces vapor (i.e., neither phase temperature is increased). Hence, if the wall heat flux is known, then the determination of the vapor generation rate follows from an energy balance. The third type of vapor generation is that which occurs when a superheated vapor transfers heat to liquid droplets thus evaporating the droplets. This form of vaporization primarily occurs after CHF has been exceeded. For these conditions, the liquid can no longer wet the heater surface and, therefore, the entire wall heat flux is transferred directly to the vapor. Due to the relatively low conductivity of the vapor, not -77- Ar Oct Oo 69 '0 so a 00 Of a c-Jo 0 I 00 ,o a0 TUBE INLET Figure 4.1: J DROPLET VAPORIZATION SATURATED FLOW BOILING SUBCOOLED BOILING Boiling Regimes in Two-Phase Flow in a Vertical Tube with Heat Addition I -78- all of the transferred heat will vaporize the remaining liquid (which takes the form of entrained liquid droplets). Rather, a portion of the heat flux will superheat the vapor with the remainder vaporizing the liquid droplets. This type of vaporization leads to substantial vapor superheat [41] and is clearly a non-equilibrium process. These three types of vaporization represent the primary vapor generation regimes for steady-state and non-depressurization transient conditions. Both subcooled and saturated boiling occur for steady-state BWR conditions, while only subcooled boiling occurs for most steady-state PWR conditions. Certain transients in either reactor type may result in all three regimes. In order to decide which regime dominates for a given set of conditions, the range of application of each must be carefully defined. As sketched in Fig. 4.2, the vapor generation rate may be considered to be a function of the equilibrium quality. It is seen that one clear dividing point is the CHF quality (XCHF). For pre-CHF conditions either subcooled or saturated boiling will occur. However, for post-CHF conditions only droplet vaporization will occur. Hence, the droplet vaporization mechanism will be postulated to occur only after CHF has been exceeded. The division between subcooled and saturated boiling is not very clear due to the gradual transition from one to the other. Hence, it is advantageous to describe both types of boiling in a single, continuous model so that the gradual transition from subcooled to saturated boiling is well represented. Furthermore, this model would be used for all pre-CHF conditions. Hence, the choice of vapor generation regime depends directly on the heat transfer regime. -79- Subcooled Boiling Saturated Boiling Droplet Vaporization 0 EQUILIBRIUM Figure 4.2: XCHF QUALITY (X ) e 1.0 Illustration of Vapor Generation Rate, r , versus Equilibrium Quality, X . ( r = a / i )e e %r fg elr' r, 0. -I------II -80- This approach of using the heat transfer regime to determine the appropriate vapor generation rate has been incorporated into THERMIT. The use of this simple selection scheme eliminates the need to have a more elaborate flow regime map. Two models are used to represent the vaporization phenomena. The first, referred to as the subcooled boiling model, is employed for all pre-CHF conditions. For the post-CHF droplet vaporization regime, the Saha [42] vapor generation model is used. In the original version of THERMIT this droplet vaporization had not been considered so that the addition of the Saha model represents a significant improvement to the code's capability. Furthermore, the subcooled boiling model was not described in the original report on THERMIT. Detailed descriptions of the physical bases as well as the results of the assessment effort are presented in the following two sections. 4.3.2 Subcooled Vapor Generation Model The subcooled vapor generation model in THERMIT accounts for both subcooled and saturated boiling in the pre-CHF regime. Since it is relatively easy to formulate a model to describe saturated boiling, the main difficulty in formulating this model is in representing vapor generation for subcooled conditions. The difficulty is that in subcooled boiling vaporization occurs at discrete sites along the heater wall for highly non-equilibrium conditions. Furthermore, the vaporization rate is found to be strongly dependent on the bulk fluid conditions as well as the wall heat flux. On a microscopic scale, the subcooled vapor generation can be directly related to the vapor bubble rate of growth. Vapor bubbles are -81- formed at nucleation sites on the heated surface as illustrated in Fig. 4.3. These nucleation sites are only activated when the wall temperature is greater than the saturation temperature. However, once the wall temperature exceeds the saturation temperature, bubbles can begin to form. The temperature distribution in the liquid permits slightly superheated liquid to exist near the wall (see Fig. 4.4). This superheated liquid is easily vaporized provided a vapor bubble site exists. As the bubble liquid temperature increases, so does the region of superheated liquid which in- turn allows the bubbles to grow. The bubbles will remain attached to the wall until the bulk liqsid temperature reaches the bubble departure temperature, Td. Once Td has been exceeded the bubbles detach and flow into the main flow stream. Hence, the bubble departure temperature (which of course is less than the saturation temperature) represents the bulk liquid temperature at which vaporization may begin. The value for Td is found to be strongly dependent on the heat flux. For the same bulk liquid temperature, as the heat flux is increased the region of superheated liquid-near the wall increases and Td will decrease. Hence, as the heat flux is increased, boiling will begin at higher subcoolings because the amount of superheated liquid is greater (see Fig. 4.4). The evaluation of the bubble departure temperature is seen, then, to be very important. This parameter has been correlated by many authors [43, 441. The correlation of Ahmad [43] has been selected for use in THERMIT. In this correlation, Td is related to the heat flux through a heat transfer coefficient. The expression for this relationship is given by -82- FLOW LIQUID VAPOR BUBBLES WALL Figure 4.3: Illustration of Vapor Bubble Nucleation and Growth -83- NOMINAL CASE T w Superheated Liquid Region 4J T sat T 0 Distance from Wall Tw HIGHER HEAT FLUX CASE Superheated Liquid Region T sat w T -- 134 bulk TD a f W 0 r Distance from Wall Typical Temperature Distributions in Subcooled BoilingFigure 4.4: Td =T -"/HA (4.1) The beat transfer coefficient HA has been correlated using a large number of experimental measurements and is given by k 1/2 1 / 3 1'i 1/3 i 1/3 (4.2) H A= - 2.44 Re Pr1 A D Li fifJ It is seen that for given flow conditions (i.e., H A and Ts constant) if the heat flux is increased, then Td decreases as expected. Hence, the proper trend of the microscopic picture is obtained in this correlation using bulk flow parameters. Even though Td is well defined by correlation, the problem of obtaining the vaporization rate based on bulk flow properties still remains. Again following Ahmad, the following physical picture can be constructed. For bulk liquid temperature below Td, bubbles do not detach and the net vaporization rate is zero. At the other limit, that is T z Ts, all of the wall heat flux leads directly to vapor generation so that the equilibrium vaporization rate, re, may be written as e ifg where is the power transferred to the coolant and i is the heat of vaporization. Ahmad then assumes that the vaporization rate increases linearly from Td to Ts. With these assumptions the vapor generation rate may be written as -85- 0 if Tk< T T -Tdk = d if T < T < T (4.4) T - T e d z s s d re ifT > T e 2-s It is seen that this model correctly defaults to the saturated boiling model once the liquid becomes saturated. Although the assumption concerning che linear increase in F may not be strictly valid for all cases, it is appropriate for most cases of interest. Hence, this model is able to realistically describe the vapor generation rate. However, this model is still not complete since vapor condensation has not been included. If the bulk liquid is subcooled, some of the vapor bubbles which detach from the wall may be condensed. Hence, the loss of vapor due to condensation must be accounted for. The model used to represent the condensation is relatively simple, but appropriate. The condensation rate, rc, is modeled as a conductibn term divided by the heat of vaporization. This can be written as rC= AiH1(Ti - Tv)/itg (4.5) The term A represents the interfacial area and if one assumes spherical bubbles of radius Rb then Ai may be written as A = (3a/R4.6) -86- Hence, if the average bubble radius can be calculated then the inter- facial area is obtained. The value for Rb is given by a modified form of the Ahmad correlation [43]: 9czR[11/3 <0.1(4.7) 9 C-( a >_ 0.1 and R o= 0.45 Z1-0V [1+1.34((1 - c)V) ] (4.8) The interfacial heat transfer coefficient, Hi, also needs to be defined. Based on the effective conductivity of the two phases, H is given by k~T T (49) O.OlR 0 kz + 0.015RD k v v ( With these definitions, the condensation model is complete. It is seen that if one assumes the vapor to be saturated, then the driving force for condensation is the amount of liquid subcooling. It should also be noted that for subcooled conditions (Tt < T '), c is negative as expected. Both the vaporization and condensation terms can be combined to obtain the net vaporization rate: 0 ifT < Td T -Td T -%Tdr + A H (T - TVT)if if T T From the previous discussion it is seen that this model is mechanistically based on the physical phenomena involved. Nevertheless the model needs to be validated by comparisons with experimental measure- ments. The main characteristics of the model include the boiling incipient point and the vapor generation rate for subcooled conditions. Both of these characteristics have been assessed using steady-state, one-dimensional void fraction measurements. The boiling incipient point, which corresponds to the bubble departure temperature, can be clearly identified in the measurements which makes the assessment of this characteristic rather straightforward. The vapor generation rate in subcooled conditions can be directly related to the void fraction if the assumption of no slip is used. This assumption will be appropriate for void fractions at low quality and high pressure. For these cases the expression for the void fraction is I -X P Iv -(4.11) x -88- The quality, in turn can be related to the vapor generation rate via the vapor mass equation (simplified for one-dimensional, steady-state conditions): ax r (4.12) 3z G Hence, the vapor generation rate can be assessed with one-dimensional, steady-state void fraction measurements. In the assessment effort, over 30 void fraction comparison cases have been made. The data of Maurer [33], Marchaterre [35], and Christensen [34], have -been used in this study. These data cover a wide range of flow conditions as seen in Table 4.2. For assessing the vaporization rate, only comparisons at low qualities have been used. Excellent agreement has been found in these comparisons for both the boiling incipient point and the subcooled void fraction. Typical comparison cases, covering a range of pressures, are illustrated in Fig. 4.5 - 4.8. The point where boiling begins is seen to be well-predicted in each case. This good agreement indicates the appropriateness of using Ahmad's correlation for the bubble departure temperature. The void fraction for subcooled conditions is also well predicted by THERMIT. This result strongly supports the use of the mechanistic subcooled vapor generation model. Other comparison cases also exhibit this good agreement, even though a wide range of conditions have been considered. In view of the above comparisons and the inherent physical attributes of the subcooled vapor generation model, it can be stated that the model satisfactorily predicts subcooled boiling. Extending this model to -89- TABLE 4.2 Test Conditions for One-Dimensional Steady-State Data Test Pressure Hydraulic Mass Flux Heat Flux Inlet Range Diameter Range Range Subcooling 2 2 Range (MPa) (mm) (kg/m s) (kW/m2) (kJ/kg) Maurer 8.3-11.0 4.1 540-1220 280-1900 150-350 Christensen 2.7-6.9 17.8 630-950 190-500 9-70 Marchaterre 1.8-4.2 11.3 600-1490 45-250 9-63 o Data -- THERMIT 0 0 0 0 1400 ENTHALPY (kJ/kg) Figure 4.5: Void Fraction versus Enthalpy for Maurer Case 214-9-3 -90- 1.0 TEST CONDITIONS p = 11.03MPa 2 G = 555kg/m-s q" = 0.306MW/m2 ai =sub 183kJ/kg D = 4.1mm e 0.75 1- 0.50 0.25 0. (13I 0/ Aofo i 1200 1300 1500 1600 F I -91- 0.50- 0.25- 1000 1200 1400 1600 ENTHALPY (kJ/kg) Figure 4.6: Void Fraction versus Enthalpy - Maurer Case 214-3-5 0.75 z0 H 0 TEST CONDITIONS 0 Data THERMIT p = 8.3MPa 2 G = 1223kg/m-s q = 1.89MW/m 2 Ssub 336kJ/kg D = 4.lmm e -- 0( 0. 1800 1.01- i -92- TEST CONDITIONS p = 4.2MPa 2 G = 842kg/m-s q = 0.24MW/m2 isub = 58.kJ/kg D = 11.2mm e 00 01 1100 o Data --- THERMIT 1150 ENTHALPY (kJ/kg) Figure 4.7: Void Fraction versus Enthalpy for Marchaterre Case 168 1. 0.75I-- 0.5 0 Pr4 0 0.25 0 . 1050 1200 u L---. I -93- 1.0 TEST CONDITIONS 0 Data P = 1l8MPa THERMIT 2 G = 862kg/m -s o.75 - q" = 0.16MW/m2 Sisub = 30.kJ/kg SD = 11.2mm 0 0.50 .. 0. 0150~ 0.25 001 850 900 950 1000 Enthalpy (kJ/kg) Figure 4.8: Void Fraction versus Enthalpy for Marchaterre Case 184 -94- three-dimensional cases also seems to be appropriate due to its mechan- istic nature. That is, the actual vaporization mechanism has been modeled on a local basis independent of surrounding control volumes. It should be noted that this model will approach the correct saturated boiling limit as the liquid becomes saturated. Therefore, the subcooled vapor generation model should be applicable for all pre--CHF conditions except for depressurization transients (in which flashing becomes the significant vaporization mechanism). 4.3.3 Droplet Vaporization Model For post-CHF conditions, the predominant form of vapor generation is evaporation of entrained liquid droplets. The reason for this is that after CHF has been exceeded the wall temperature will rapidly increase and in a short period of time the minimum stable film boiling temperature, Tfb, will be exceeded. Once this temperature has been attained, the liquid can no longer receive heat directly from the wall. Instead, the vapor is in contact with the wall and only by vapor-to-liquid heat transfer can the liquid be heated and evaporated. Hence, the rate of vapor generation is directly dependent on the rate of heat transfer from the vapor to the liquid. However, due to the low conductivity of the vapor, the vapor-to-liquid heat transfer is not very efficient. Consequently, the vapor becomes superheated by a significant amount (e.g., 150 *K.441]). The key to predicting the correct vaporization rate is to carefully model the heat transfer between phases. Once this heat transfer rate is determined, the vaporization rate is found by simply dividing by the heat of vaporization: -95- P = AiH (Tv - T ) / i fg(4.13) where A and H are the appropriate interfacial area and effective heat transfer coefficient. The interfacial area is dependent on the droplet diameter while the heat transfer coefficient will depend on the flow conditions, droplet diameter and vapor conductivity. The temperature difference T - T Zmay be written as T - T if it is assumed that the vV 5 liquid is saturated. Hence the vapor must be superheated in order for vaporization to occur. However, if r is zero the vapor will superheat since it receives heat from the wall without losing any of it to the liquid. Consequently, there is a direct coupling between the amount of vapor superheat and the vaporization rate. The difficulty in determining r is in defining relations for A. and H . As discussed by Saha [42] each of these parameters may be written as a function of the flow variables, but ultimately a correlation is required to complete the function. The interfacial area per unit volume may be written as A =6(1-a)(4.14)i 6 where 6 is the droplet diameter. This diameter is strongly dependent on the flow conditions and is, therefore, usually empirically correlated. The interfacial heat transfer coefficient is correlated as a Nusselt number based on the droplet diameter: W96- k (P(Vvz)6)0 5 5 0.33 H, - 2 + 0.459(J Pr (4.15) Again 6 needs to be determined from correlation. Hence, both A and H must be correlated as functions of the flow conditions. In view of this difficulty, Saha has combined the two parameters, A and Hi, into a single parameter K which is then correlated as a function of the flow conditions. This approach eliminates the need to use two correlations which may be difficult to determine separately. A wide range of conditions have been used in developing this correlation as illustrated in Table 4.3. The final form of this vaporization rate correlation is given by 2r 2-1/2 -2 pVyD k (T -T) P'-=6300 1 - -P-vj v v v s (4.16) 1 Pcr_ L a - D 2 1 ifg The droplet diameter has been assumed to be proportional to the hydraulic diameter, D. The interfacial area per unit volume is seen to be inversely proportional to D with the heat transfer coefficient being proportional to k /D. As the vapor velocity increases, the droplets become smaller, v increasing the interfacial area and increasing P. Hence, this model apparently contains sufficient physical characteristics to predict the vaporization of liquid droplets. Obviously, the important quantity which this model is intended to predict is the rate of vapor generation for post-CHF conditions (or whenever vaporization of liquid droplets is significant). Unfortunately, this rate cannot be directly measured. Consequently, the assessment of the Saha model has required indirect methods. This assessment relied on TABLE 4.3 Test Conditions used to Develop SaIha Correlation For Post-CHF F Pressure 1.5 MPa and 6.9 MPa Diameter 14.9 mm and 12.6 mm Mass Velocity 393-2600 kg/M2S Heat Flux 0.045 to 0.127 MW/M2 Equilibrium Quality 0.18 to 1.50 -98- the tight coupling between the amount of vapor superheat and the vaporiza- tion rate. For post-CHF conditions, the entire wall heat flux is received by the vapor. Part of this raises the vapor temperature with the remainder evaporating the liquid. Where the wall heat flux is known, an indirect assessment of the vaporization rate can be made if the fraction of the heat flux which raises the vapor temperature can be determined. This fraction is easily calculated if the vapor temperature is known. Unfortunately, the vapor temperature is not easily measured [41]. However, the vapor temperature can be inferred from wall temperature measurements using the known heat flux and an appropriate heat transfer coefficient. This method is straightforward provided the heat transfer coefficient is judiciously chosen. Since the heat transfer mechanism is primarily forced convection to the vapor, Saha recommends the use of a single-phase vapor forced convection heat transfer correlation for the heat transfer coefficient. Many such correlations are available that are reasonably accurate fpr high pressure conditions. Hence, the procedure for assessing the post-CHF vapor generation model is to make comparisons between measured wall temperatures and the predictions of THERMIT for post-CHF conditions. The wall temperature comparisons give a direct indication of the vapor temperature predictions which, in turn, relate to the vapor generation rate. Though this proce- dure is indirect, it is the only viable method for assessing the post-CRF vapor generation rate, which is not measured directly. For this study, the steady-state, one-dimensional wall temperature measurements of Bennett [391 have been used. Test sections having lengths -99- of either 5.5m or 3.3m have been used. Both pre-CHF and post-CHF temperatures were measured. A wide range of conditions, especially for the mass velocity, have been used in these tests as summarized in Table 4.4. It is seen that these conditions are close to those found in an operational BWR. A number of representative cases have been simulated with THERMIT. Comparison of measured and predicted post-CHF wall temperatures are in overall good agreement even though a range of conditions have been considered. Example comparisons are illustrated in Figures 4.9 - 4.11 for a range of mass velocities. It should be noted that the predicted CHF location has been adjusted to correspond with the measured value in order to accurately assess the post-CHF regime (assessment of the CHF predictive capabilities are discussed in Chapter 6). In each case, the trend in the data is correctly predicted and good agreement is found except near the test section exit. These discrepancies are probably due to axial conduction effects which are not modeled in THERMIT. As a means of comparison, two other vapor generation models which represent limiting values are compared with the Saha model and measurements in Fig. 4.12. The first model, termed the equilibrium model, assumes that all of the wall heat flux leads to vapor production and may be written as I' = (4.17) fg -100- TABLE 4.4 Bennett Test Conditions for CHF in Tubes Outlet Pressure 6.9 MPa Diameter 12.6 mm 2 Mass Velocity 664 to 5180 kg/m s 2 Heat Flux 0.56 to 1.77 MW/M Inlet Subcooling 72 to 146 kJ/kg Tube Length 3.66 and 5.56 m -101- 1100 o Data - THERMIT 1000 900 0 800 TEST CONDITIONS P = 6.9MPa 700 2 G = 665kg/m-s q" = 0.65MW/m2 600 -i = 137.kJ/kg D = 12.6mm e 500 3.5 4.0 4.5 5.0 5.5 AXIAL HEIGHT (m) Figure 4.9: Wall Temperature Comparisons for Bennett Case 5332 (Length = 5.56m) -102- 1100 o Data - THERMIT 1000 900 -- 0 O 0 800-E-4 rI TEST CONDITIONS 700 P = 6.9MPa 2 G = 1356kg/m s q" = 0.90MW/m2 600 a isub = 102 kJ/kg D = 12.6mm 500 3.5 4.0 4.5 5.0 5.5 AXIAL HEIGHT (m) Figure 4.10: Wall Temperature Comparisons for Bennett Case 5253 (Length = 5.56m) -103- 1100 o Data - THERMIT 1000 ~o 0 0 900 800 TEST CONDITIONS P = 6.9MPa 2 700 G = 4815kg/rn-s q"1 = 2.07MW/m2 ai sub = 86.kJ/kg D = 12.6mm 600 - 0-- 500 2.0 3.0 4.0 AXIAL HEIGHT (i) Figure 4.11: Wall Temperature Comparisons for Bennett Case 5442 (Length = 3.66m) -104- 1100 o Data -- THERMIT (Saha) 1000 - Frozen Quality - - Thermal Equilibrium 900 oo .0 800 TEST CONDITIONS 700 P = 6.9MPa 2 G = 665kg/m - s q" = 0.65MW/m 2 600 Ai=sub 133.kJ/kg D = 12.6mm e 500 - . - -- 3.5 4.0 4.5 5.0 5.5 AXIAL HEIGHT (m) Figure 4.12: Comparison of Wall Temperature Predictions Using Various P Models for Bennett Case 5332 -105- When this model is used, the vapor does not superheat until all the liquid is evaporated and consequently the wall temperatures are under- predicted. The second model, called the frozen quality model, assumes that F is zero after the CHF point which prevents the qrality from changing. When this model is used, no evaporation is allowed so that all the wall heat is transferred to the vapor. Consequently, the vapor superheat as well as the wall temperatures are significantly overpredicted. The comparisons with these two limiting models as well as with the data demonstrate the appropriateness of the Saha model for post-CHF vapor production. This model has been mechanistically developed and correlated with data which cover a wide range of conditions. Consequently, based on the physically motivated formulation and good agreement with data it is appropriate to use this model for post-CHF conditions. 4.4 Interfacial Energy Exchange A second interfacial phenomenon which must be modeled is the inter- facial energy exchange rate. This energy transfer is directly related to phasic temperatures and, hence, controls the thermal non-equilibrium. Since the ability to predict thermal non-equilibrium is a key feature of THERMIT, appropriate modeling of the interfacial energy exchange is essential. The interfacial energy exchange rate represents the rate of energy transfer from one phase to the other. This transfer can be due to either conduction, which is a function of the temperature distribution, or mass transfer. The physical picture for this transfer can be explained with -106- the aid of Fig. 4.13. If the interface is assumed to be of a discrete, but infinitesimal, size and at saturated conditions then the energy transfer can be modeled. For example, in Fig. 4.13a, which illustrates the thermal field for subcooled boiling conditions, the liquid adjacent to the vapor bubble interface is superheated while the vapor is saturated (or slightly subcooled). Defining the energy transfer as positive when the vapor receives the energy, the energy transfer rate may be written as: Qi= Hi (T - Ts ) + if = Hiv (T- Tv ) + ri (4.18) where H. is the liquid-to-interface heat transfer coefficient and Hi is the vapor-to-interface heat transfer coefficient. This equation shows that the rate of energy transfer from the liquid to the interface is the same as the- energy transport rate from the interface into the vapor. In view of the equivalence of energy transfer rates, one may use either form. A second example of the interfacial energy exchange is illustrated in Fig. 4.13b. In this case, which represents the liquid droplet vaporization regime (i.e., post-CHF), the vapor is superheated while the liquid is saturated. Even with this very different temperature profile, Eq. 4.18 still describes the interfacial energy exchange. Hence, Eq. 4.18 describes the general form of the interfacial energy exchange. Notwithstanding the generality of Eq. 4.18, a practical problem remains in choosing the appropriate form of the interfacial energy exchange (i.e., either liquid-to-interface or interface-to-vapor). The choice of either form is dictated by the assumed temperature distribution -107- LIQUID T i--, T- INTERFACE T S VAPOR T Figure 4.13a: Temperature Distribution Near Vapor Bubble for Subcooled Boiling Figure 4.13b: Temperature Distribution Near Liquid Droplet for Droplet Vaporization LIQUID INTERFACE VAPOR f g T I a -108- for a given set of flow conditions. In THERMIT, the flow conditions have been classified as either pre-CHF or post-CHF, corresponding to the mass exchange rate models. Again the pre-CHF regime includes both subcooled and saturated boiling, while in the post-CHF regime only vaporization of entrained liquid is considered. The choice of interfacial energy exchange model will be determined by the assumed temperature distributions for these types of boiling conditions. In subcooled and saturated boiling conditions, a slightly superheated liquid exists at the interface and transfers energy to the vapor. However, the bulk liquid temperature is subcooled so that it would be difficult to use the liquid-to-interface energy transfer mechanism to represent the interfacial energy exchange without doing a detailed analysis of the temperature distribution in the liquid. Hence, the liquid-to-interface energy transfer will not be considered here. On the other hand, the interface-to-vapor energy transfer can be appropriately modeled by considering the vapor to be at saturated conditions for both of these types of boiling conditions. In order to maintain the vapor at saturated conditions when the bulk liquid is subcooled, a relatively high rate of heat transfer across the interface must be assumed. This high rate can be interpreted as a large value for the vapor-to-interface heat transfer coefficient Hiv. Hence, if Hiv is chosen sufficiently large, the vapor will be maintained at saturation. Consequently, for subcooled and saturated boiling conditions the interfacial energy exchange is modeled as an interface-to-vapor energy transfer mechanism. This exchange rate can be written as -109- Q = H (T - T) +P9i (4.19) where Hiv is set to a very large value (1011W/m3) in order to force the vapor to be saturated. It should be noted that since the bulk liquid temperature is not used in this equation, the liquid is unconstrained and may, therefore, be subcooled. Hence, the use of Eq. 4.19 for the interfacial energy exchange rate permits appropriate modeling of both the bulk liquid and vapor temperatures in subcooled and saturated boiling. For post-CHF conditions, where droplet vaporization is the form of mass exchange, the superheated vapor is assumed to transfer heat by conduction to the interface while receiving energy due to the vaporiza- tion of the liquid. In this case, modeling of the vapor-to-interface energy transfer is difficult unless the detailed vapor temperature distribution is known. However, the liquid-to-interface energy exchange -can be-adequatelyi-odeled-since-the -iiquid-is-assumed-to-be-at--or-near--- saturation. Therefore, by simply choosing a value for H which is sufficiently large, the liquid will be forced to saturated conditions. Consequently, for the droplet vaporization regime, the interfacial energy exchange is modeled as a liquid-to-interface energy transfer mechanism. This exchange rate may be written as Q= F if - Hi (Ts - Tz) (4.20) where Hit is set to a large value (101k7/m3) in order to force the liquid to saturation. The bulk vapor temperature is not constrained by this -1 10- equation which allows che vapor to superheat. Hence, this model allows for the appropriate liquid and vapor temperatures to be predicted for the droplet vaporization regime. In spite of the mechanistic nature of these two interfacial energy exchange models, assessment of these models is still required. However, validation of either model is not possible since the interfacial energy exchange cannot be directly related to a measureable quantity. Therefore, these models can only be assessed qualitatively by inference which means, when used, the models should produce the expected results. For example, in subcooled conditions the bulk liquid temperature should be subcooled while the vapor should be saturated. Alternatively, for droplet vaporiza- tion, the vapor should be superheated with the liquid saturated. If these results are predicted, then the interfacial energy exchange rate is at least qualitatively correct. These models have been used in all of the mass exchange rate valida- tion studies and have yielded the expected results in all cases. A typical temperature profile is illustrated in Fig. 4.14 for one of the void fraction comparison cases. It is seen that the vapor temperature follows the saturation temperature which is decreasing due to the pressure drop. The liquid temperature is initially subcooled, but eventually reaches saturation near the end of the test section. Hence, for subcooled and saturated boiling conditions the interfacial energy exchange rate given by Eq. 4.19 seems to be an appropriate choice. For post-CHF conditions similar results are obtained. The temperature distributions for one of the Bennett cases is illustrated in Fig. 4.15. At the inlet, the liquid is subcooled but quickly becomes saturated and -111- 600 VAPOR TEMPERATURE - LIQUID TEMPERATURE 550 500 450 - 0. 0.25 0.50 0.75 1.0 RELATIVE AXIAL HEIGHT Figure 4.14: Predicted Liquid and Vapor Temperatures for Maurer Case 214-3-5 / / - - VAPOR TEMPERATURE -- LIQUID TEMPERATURE 700 650 > 600 550 500 4.0 5.0 AXIAL HEIGHT (i) Figure 4.15: Predictcd Liquid and Vapor Temperatures for Bennett Case 5336 -112- / / / / // / / / / / / / / / / / / / 0. 1.0 I 2.0 3.0 j Ii I I -113- remains so along the entire heated length. The vapor remains at satura- tion before CHF, but quickly superheats after CHF has been attained. These predictions are the expected result so that the interfacial energy exchange model given by Eq. 4.20 seems to be an appropriate choice for post-CHF conditions. The above assessment has been based on the expected value for the liquid and vapor temperatures. As such, only the conduction term in the energy exchange model could be assessed. The mass exchange term (either Fi. or fi is also important, but not for steady-state conditions. That is, in steady-state the same temperature distribution is obtained whether or not the mass exchange term is included. However, for transient conditions it is essential that the mass exchange term be included. The reason for this can be understood by considering a case in which CHF occurs and the liquid is still subcooled (e.g., DNB type CHF). If the mass exchange is not included, then Q before CHF is Q,= Hi (T - Tv) (4.21) and after CHF Q- Hi (Tz - TS) (4.22) The pre-CHF expression is large and positive, i.e., T > Tv, while the post-CHF expression is large and negative. This sudden change in the value for Q represents a severe discontinuity and prevents convergence of the code. -114- If the mass exchange term is included, then Q will still change value but the discontinuity is not as severe. In this case, the code can converge despite the discontinuity. Hence, proper modeling of the interfacial energy exchange rate is essential for both steady-state and transient conditions. 4.5 Interfacial Momentum Exchange The third type of interfacial exchange phenomena which must be modeled in THERMIT is the interfacial momentum exchange. This exchange represents the transfer of momentum from one phase to the other and cor-rols therelative velocity of the two phases. As in the case of the other interfacial exchange phenomena, the interfacial momentum exchange is strongly dependent on the flow conditions, since the structure of the two-phase flow changes with the flow conditions. As illustrated in Fig. 4.16, if the vapor concentration is low, then a bubbly flow is expected in which vapor bubbles move through a continuous liquid medium. As the vapor concentration increases the bubbles agglo- morate in the center of the flow channel. At higher concentrations, an annular flow is found in which liquid coats the wall with the vapor forming a continuous central core. Of course, the possibility of liquid droplets in a continuous vapor phase also exists after CHF. In each case, as the flow pattern changes, the interfacial area changes. Since the momentum exchange is directly proportional to the interfacial area, the flow conditions are seen to have a strong influence on the interfacial momentum exchange. In attempting to model the interfacial momentum exchange, it is necessary to consider the various forces which can act between the two Figure 4.16: Typical Flow Patterns in Two-Phase Flow (3 0 0 Q 00 o 0 0 o 00 oQ 0 a '3 0 00 00 0~ 0 0 0 0 a 0 C UQJQooS a* * C)) .-0 I * Increasing Void Fraction ANNULAR 0 6 0 Ut O 0 0 00oo to 0 0 I, 0 C 0 0 t FLOW 0 0 a a 0 3 0 0 a a 80 0 0 0 I, 0 0 0 0 a 0 BUBBLY i -116- phases. At least five different forces can be postulated to exist. These may be divided into steady flow and transient flow forces. The steady flow forces include viscous, inertial and buoyancy forces while the transient flow forces include the Basset and virtual mass force [31]. Each of these forces will be significant for certain conditions and, hence, it is important to understand the characteristics of each force. The viscous force, which arises due to the viscous shear stress and is only significant at low relative velocities, is approximately described by Stokes law. This law, originally derived for the force on a sphere moving in a viscous fluid, has since been modified to account for the motion of droplets of one fluid in a continuous second fluid. The flow of vapor bubbles in a liquid and the flow of liquid drops in a vapor are examples of such motion. The force on a single solid sphere given by Stokes law can be written as F =3wTrpDsV (4.23) p cs r where Ds is the sphere diameter, pc is the viscosity of the continous phase and Vr is the relative velocity. As discussed by Soo [45], modifi- cations of this equation are required for systems in which the droplet (or sphere) is deformable (such as vapor bubbles in liquid). An example of such a modification is given by Levich [46]: F = 6rc D V (4.24) y c dsr where the subscript c refers to the continuous phase and the subscript d refers to the dispersed phase. This expression is similar to other -117- expressions [45] and is valid for many practical droplet or bubble flow situations. The force given by Eq. 4.24 represents the force on a single droplet. In order to convert this to a force per unit volume, Eq. 4.24 must be divided by the volume of a droplet and multiplied by the void fraction. Performing this operation yields 36pa V F = c r (4.25) y 2 d This expression represents the interfacial force due to viscous effects within a given control volume. The second type of force is that due to inertial effects. This force also referred to as the drag force, represents the momentum loss due to the motion of two continuous fluid streams relative to one another. Hence, this force tends to dominate in annular flow regimes. Following Wallis [31], the shear stress between the phases may be written as Ti = C 2PV (4.26) Sd v r Since the diameter of the vapor core is given by D = DvW (4.27) c the interfacial force per unit volume is - 118- 2C PV 2 F. = d yr (4.28) D where Cd is the interfacial drag coefficient. Values for Cd, appropriate for annular flow, have been formulated with Wallis recommending the following value [31], Cd = .005(1+ 75(1-ct)) (4.29) Using this coefficient the interfacial drag force can be evaluated. As indicated above, this force will be significant for annular flow and when the relative velocity is greater than zero. The third type of force is that due to buoyancy effects. This force arises due.to, the difference in densities of the two phases. In a gravitational field (e.g., vertical flow), this density differences causes a force between the two phases. This force may be written as: Fg = a( - a)O(p- Pv)g (4.30) The buoyancy force will only be significant for low velocity flows or when the other forces are small relative to this force. The fourth force is that due to virtual mass effects. This force arises from the apparent increase in mass of an accelerated particle. When a particle is accelerated relative to the surrounding fluid a potential flow field possessing kinetic energy will be established. The particle effectively accelerates this surrounding fluid, termed virtual -119- mass. Since the virtual mass is accelerated, it represents an additional force on the particle. Hence, as a particle is accelerated its mass appears to increase which leads to an increase in the interfacial force. The virtual mass interfacial force has been discussed by Wallis [31] and Cheng et al. [47] and may be written as [48]. F 1 [f+21jjc d(V dV C)(4.31) vm 2d. 1- ad c dt where the subscript d refers to the dispersed phase and c refers to the continuous phase. It is seen that this force depends on the rate of change of the relative velocity and teniei to decrease the lag between the phase velocities. This force will only become significant when one phase is accelerated more rapidly than the 'ther. The final force,potentially important in rapidly accelerating ftows, is the Basset force. This force arises from the fact that as a particle is accelerated, a viscous flow field is established around the particle. This flow field introduces boundary layer development which tends to increase the drag on the particle. Unfortunately, this force is difficult to calculate, since it depends on the previous flow history. For laminar flow Basset [49] has derived the following analytical expression: -t 9ct P4 t3(V -V.) 9d c c c d dt (4.32)F /= V Basset Dd JT td 3 tt -120- where again the subscript d refers to the dispersed phase while the subscript c refers to the continuous phase. This force represents an instantaneous flow resistance with the previous flow history contained in the time integral. Wallis [31] shows that if laminar flow and constant acceleration are assumed, then the ratio of the Basset force to the steady drag force is given by FBasset /F Drag = Dd/ /t (4.33) Ec At very small times, this ratio will be large and, hence, the Basset force will be significant. Consequently, for rapidly accelerating flows and short times the Basset force represents a significant interfacial force. As indicated by the above descriptions of the important interfacial forces, each force models a specific interaction and is significant only for certain conditions. These characteristics are summarized in Tabla 4.5. Although not explicitly indicated, all forces except the inertial force have been formulated based on dispersed or bubble flow conditions. It is also worth noting that, in the THERMIT interfacial momentum exchange model, only the viscous and inertial forces, which will dominate the interfacial force term for steady flow or near-steady reactor conditions have been included. For rapidly accelerating flows which are anticipated in blowdown transients both the virtual mass and Basset forces may need to be included. However, since this type of transient is not the primary application of the current research, the exclusion of these transient forces is probably justified for cases of practical interest. The buoyancy force is expected to become significant when the relative -121- TABLE 4.5 Summary of Liquid-Vapor Interfacial Forces Force Equation No. Applicable Range InT ein Viscous 4.25 Low relative velocity, Yes Bubbly or Droplet Flow Inertial 4.28 Annular Flow Yes Buoyancy 4.30 Low flow conditions, NO Bubbly Flow Virtual Mass 4.31 Rapidly accelerating NO flow Basset 4.32 Rapidly accelerating NO flow, Short times -122- velocity is small. Such cases exist for low flow conditions or when the wall friction is very small which are not typical of reactor conditions. Hence, for typical reactor operating conditions this force may be neglected. In view of the simplifications to the interfacial force, it is important to remember that arbitrary application of THERMIT to different conditions is not warranted unless the appropriate interfacial forces are included. However, the forces have been clearly identified and, therefore, extention of the interfacial momentum exchange rate can be easily accomplished. Turning now to the actual interfacial momentum exchange model in THERMIT, it is important to reiterate that only the viscous and inertial forces have been included in this model. Furthermore, in order to avoid the use of a flow regime map, this model has been formulated to be continuous for all flow regimes. Flow regime maps based on the void fraction have been used in other two-fluid model codes (for example TRAC [12]). In these maps, the void fraction determines the flow regime which then defines the appropriate interfacial momentum exchange rate to be used. However, this type of formulation is probably not warranted for the present applications. Hence, a continuous interfacial momentum exchange model has been developed. In this model, the coefficients of both the viscous and inertial forces have been approximated by simple functions of the void fraction. This approximation produces the desired numerical result while yielding appropriate values for the coefficients. Taking both forces together, the interfacial momentum exchange model in THERMIT is -123- F = - V + a v r r(4-34) aD k r aD 2 where c = max(0.1,ca) and D = hydraulic diameter V = V - V r v Z The reason for the restriction on a is to prevent a singularity when a = 0. From the previous discussion, it should be obvious that the first term in this expression represents the viscous force while the second term represents the inertial drag force. Comparing the viscous term with Eq. 4.25, one finds that the following approximation has been made: 36a ~ l -4V D 2aD (4.35) where Dv is the vapor bubble diameter appropriate for bubbly flow. Since this force is only significant in bubbly flow regime, the approximation here is only appropriate for low void fractions. This fact is illustrated in Table 4.6 where the two coefficients are compared for a range of void fractions assuming representative values for the diameters. Only for void -124- TABLE 4.6 Comparison of Viscous Force Coefficients 36ct D2 *1* 1~ 4.0 5.0 5.8 6.4 6.8 7.3 x10 5 x10 5 x 10 x 10 x 10 0 m Assumptions D = 0.01 m D =2Qcx 4Tr 1/3 v 3 with N = 10 bubbles/m3 0.05 0.10 0.15 0.20 0.25 0.30 2 [1.- otD 8.1 x10 5 58.l xlO 3.2 x 10 1.6 x 10 9.0 x 104 5.0 x 104 -125- fractions of approximately 0.15 and less is the approximated coefficient comparable with the Levich model coefficient. However, this range corresponds to the conditions for which the viscous force is important. Furthermore, since the Levich model is only typical of the viscous force, the THERMIT viscous force term seems to be appropriate for its intended use. The inertial force term in the interfacial momentum exchange model can be compared with Eq. 4.28. In order to equate the two expressions, the following approximation must be made: 2C 2 1 a(4.36)d 2at These two coefficients are compared in Table 4.7 over a range of void fractions. It is seen that at low void fractions the THERMIT model predicts a higher coefficient which is necessary to have continuity between the viscous and inertial regimes. However, at higher void fractions the two are approximately the same. Since annular flow would be expected for a>0.6, the approximated inertial drag coefficient in THERMIT seems to be appropriate. Hence, the formulation of the interfacial momentum exchange model seems to be satisfactory in spite of the approximations which have been made. However, as in the case of the other models, validation and assessment of the interfacial momentum exchange rate is the key to successful use of THERMIT. The assessment of the interfacial momentum exchange model has employed the same one-dimensional void fraction measurements used to -126- TABLE 4.7 Comparison of Inertial Force Coefficients cc0.01(1 + 75(1 -c)) U 2a 0.4 0.29 0.75 0.5 0.27 0.50 0.6 0.24 0.33 0.7 0.20 0.21 0.8 0.14 0.13 0.9 0.08 0.06 -127- assess the interfacial mass transfer rate. While the verification of the mass exchange model was concerned with the low quality void fractions, assessment of the momentum exchange rate has relied on the high quality data. The reason for this is that only for thermal equilibrium conditions (i.e., non-subcooled conditions), can the void fraction measurements be used to independently assess the momentum exchange rate. This fact can be illustrated by considering the definition of the void fraction: X CLX=V(4.37) K + (1 - K) ____ P ZVz For a given pressure, the void fraction is seen to depend on the flow quality and the slip ratio, S (S = Vv/Vz). The flow quality has been shown to depend on the vapor generation rate by Eq. 4.12, while the slip ratio depends on the interfacial force. If the flow quality is not known, then the slip and, hence, the interfacial force cannot be determined from the void fraction alone. Fortunately, for thermal equilibrium conditions the flow quality can be calculated since the vapor generation rate simply becomes r = qw/if (4.38) with the wall heat transfer term, c, already known. Hence, the flow quality can be determined from an energy balance so that the momentum exchange rate can be assessed with void fraction measurements. As indicated in Section 4.3.2., a large number of void fraction comparison cases have been made. For assessing the interfacial momentum -128- exchange rate, only the higher quality data have been used. Generally, the THERMIT, predictions agree rather well with the measured void fraction values over the range of flow conditions considered here. Typical comparison cases, covering a wide range of pressures, are illustrated in Figures 4.17 - 4.20. It is seen that in the higher quality regimes, the measured void fraction values are satisfactorily predicted in each case. Hence, considering the range of flow conditions which have been analyzed, the interfacial momentum exchange model can be expected to be appropriate for most cases of practical interest. The only minor deviation between the measured and predicted values occurs for some of the lower pressure cases. In these cases, THERMIT tends to underpredict the void fraction. This indicates that the slip ratio is too high or, in other words, the interfacial momentum exchange rate is too low. In order to assess these deviations an alternative interfacial momentum exchange model has been added [4]. This model, referred to as the LASL model [50] due to its usage in many of the LASL codes, is similar to the THERMIT model in that only viscous and inertial forces are considered, but the coefficients are different. This model is given as: 3 12p;+ PIV vj F[=- - + (V - V) (4.39) i 8 2 where P= a p +(l-cz)Pp -12 9- o Data - THERMIT 00 TEST CONDITIONS P = 8.3 MPa 2 G = 799.kg/m -s o if = 1.09 MW/rm2 ai sub= 148 kJ/kg D = 4.1mm0 I -~ 4J 1600 ENTHALPY (kJ/kg) 1800 2000 Figure 4.17: Void Fraction versus Enthalpy for Maurer Case 214-3-4 1.0 0.75k z 0 0.50 0 0.25 1- 0 . 1200 1400 I I I m -130- 1.0 0 Dat'a TEST CONDITIONS p = 4.1 MPa - THERMIT 0.75 G = 929kg/m2- s q" = 0.35 MW/m 2 4 isub = 34 kJ/kg . - 0.0 0 4j 0.25 1050 1100 1150 1200 ENTHALY (kJ/kg) Figure 4.18: Void Fraction versus Enthalpy for Christensen Case 12 1.0- 1 Data TEST CONDITIONS - THERMIT P = 2,8 MPa 2 0.75 G = 641 kg/m s q" = 0 .21 MW/rm2 6Ai = 11 kJ/kg Z sub D =17.8 mm HO 0.50 - 0 0 0 H 0 o 0) 0.25 *rI 950 1000 1050 1100 ENTHALPY (kJ/kg) Figure 4.19: void Fraction versus Enthalpy for Christensen Case 9 -131-- -132- 1.0 0 Data TEST CONDITIONS - THERMIT P = 1.8 MPa 2 0.75 G = 863kg/m-s q" = 0.16MW/m 2 a isub = 30kJ/kg 0 D e =11.2mm 0.50 0 0 000.25 0. 850 900 950 1000 ENTHALPY (kJ/kg) Figure 4.20: Void Fraction versus Enthalpy for Marchaterre Case 185 -133- I (i = (ap/pc + (1-COP/P)P a a1/3 2/3 n r = (6/a )1/3 n 3 N 2 10 /M3 =c aif a < 0.5 1-c aif a > 0.5 The most significant difference between these two models is the density used in the inertial force term. In THERMIT, the vapor density is used while in the LASL model the mixture density is used. The mixture density can be as much as 10 times greater than the vapor density, so that the LASL momentum exchange rate can be significantly larger. The LASL model has been used in THERMIT and a number of the void fraction cases have beexn repeated. In all cases, the LASL model predicts higher void fractions than THERMIT and typically overpredicts the measured data. One such comparison case is illustrated in Fig. 4.21. The vcid fraction at low qualities is approximately the same for both models. However, at high qualities the LASL model predicts significantly higher void fractions. The slip ratio predictions for this case using both of these models is illustrated in Figure 4.22. It is seen that with the LASL Fi model the slip ratio quickly attains a value of "'1.24 and remains constant along the remainder of the tube. When the THERMIT F model is used, the slip ratio increases rapidly at first, but then.levels off near the end of the channel. -134- 1.0 TEST CONDITIONS o Data p = 4.1 MPa -- THERMIT 2 G = 929kg/m-s THERMIT with 2 0.75 q = 0.35MW/m LASL F. Model 3. Aisub = 34kJ/kg D = 17.8mm e H .50 r 0 . - U) / / 0.2 1050 1100 1150 1200 INTHALPY (kJ/kg) Figure 4.21: Void Fraction versus Enthalpy foi Christensen Case 12 -135- - THERMIT - - THERMIT with LASL F. Model 3. 2.0 -. - Homogeneous Model 0 H 1.5 100 0. 0.25 0.50 0.75 1.0 RELATIVE AXIAL HEIGHT Figure 4.22: Comparison of Predicted Slip Ratios for Christensen Case 12 -136- Of course, if a homogeneous model is used, a slip ratio of 1.0 is an imposed restriction. Since a lower slip ratio means a higher void frac- tion, the homogeneous model void fraction predictions will be the highest for a given quality. However, since the measurements lie between the predictions of the two models, the homogeneous model void fraction pre- dictions will be too high. Hence, an advantage of the two-fluid model is the ability to predict the velocity profiles for each phase. These profiles are needed for accurate void fraction predictions. In order to evaluate these two models in a more consistent manner, the measured and predicted values for Christensen's tests have been plotted on a superficial vapor velocity (Jv) versus a graph. As seen in Figure 4.23, the data tend to fall within a band between the two models. The THEERMIT model forms the lower void fraction edge while the LASL model forms the upper edge. The large amount of scatter makes it difficult to definitively state which model is better. Clearly, the THERMIT model should have an increased momentum transfer rate at high Jv. In conclusion, it can be asserted that the:momentum exchange rate in THERMIT appears to be appropriate for most steady and near-steady flow reactor conditions. Rapid transients or very low flow cases may not be adequately analyzed, since the appropriate forces have been neglected. Also the use of the current model for analyzing droplet flow may not be strictly valid. Nevertheless for a large number of cases of practical interest, the current interfacial momentum exchange rate seems to be quite satisfactory. -137- 0.4 VOID FRACTION Figure 4.23 : Vapor Superficial Velocity versus Void Fraction for Christensen Data 2.0 1.5 1.0 02PI * Data /, THERMIT F. Model 1 0 S**. a. 0.5 I 0. 0. 0.2 0.6 F L q -138- 5.0 ASSESSMENT OF THE TWO-PHASE MIXING MODEL 5.1 Introduction As discussed in Chapter 3, the original version of THERMIT has been modified so that the effects of turbulent mixing between cells are now accounted for. This modification consisted of adding a two-phase mixing model that has been formulated on physical bases using the recommendations of previous work [28, 29, 30]. Specific details of the model formulation have been presented in Section 3.3. An integral part k. of the model developmental effort has been the assessment of the model. This assessment has been made by comparing experimental measurements in rod-bundles for typical BWR and PWR conditions to the predictions of - THERMIT. In this chapter, the results of these comparisons are presented and discussed. As a prelude to the discussion of the results of the assessment, it is instructive to review the important phenomena being represented. The effects of turbulence, that is the transport of mass, energy and momentum due to turbulent eddy diffusion, are, on a LWR assembly scale, very localized. Even though the effect of the transport is on a scale of the order of subchannel sizes, turbulence can play a significant role for subchannel applications. To analytically describe the motion of the eddies would conceivably be the best way to model the effects of turbu- a. lence. However, this type of model is beyond the scope of the current work and a simpler, but practical engineering approach has been adapted. In this approach, the integral effect of the localized eddy transport mechanisms are embodied in the two-phase mixing model. In this context, -139- the integral effect means that all turbulence effects are lumped into a net mixing rate. Hence, a macroscopic modeling approach is used to represent the very complicated small-scale turbulent motion. In the integral approach two important phenomena are represented. The first, termed turbulent mixing, results from the motion of turbulent eddies in the flow. As the turbulence level increases, so does mixing between adjacent channels, resulting in a more uniform flow. Turbulent mixing tends to eliminate flow and enthalpy gradients by promoting flow exchange between the channels. It should also be noted that turbulent mixing is important for both single-phase and two-phase flow. The second phenomenon, termed vapor diffusion, represents the observed tendency of the vapor to migrate to the unobstructed (i.e. more open) regions of the rod bundle. Since this transport mechanism occurs in the absence of pressure gradients, the vapor is said to "diffuse" to the unobstructed regions. This apparent diffusion of the vapor leads to the transport of energy and momentum in proportion to the rate of vapor transport. Both turbulent mixing and vapor diffusion occur in the absence of pressure gradients and lead to mass, energy and momentum transport. In view of this, the analytical description of both phenomena can be included in a unified two-phase mixing model. As presented in Section 3.3, the generalized mixing terms can be written as $ = /1 ( - (00)J- ((0)i - (0$) )FD1 (3.33) where 0 = 1 in mass equations - e in energy equations v-140- = V in momentum equations p = ap. in vapor equations = (1-a)p in liquid equations In the mass equations this term represents an additional mass exchange rate between adjacent channels. In the energy equations this term represents the apparent heat flux at channel boundaries. Finally, in the momenum equation this term is an apparent shear stress acting along the channel boundary. The term E/, which represents the turbulent velocity, has been assumed to be the same for both the liquid and vapor. This assumption means that the mixing rate for each phase is determined from the turbulence level of the total flow and not from the turbulence level of each phase. With this assumption, if the liquid and vapor equations are added together, then the liquid and vapor mixing terms combine to yield the mixture model formulation of Lahey [25]. This result represents the appropriate limiting value and it is important that the THERMIT two-phase mixing model reduce to this limit. It should also be noted that for single-phase conditions, the model reduces to the correct single-phase limit. As discussed in Section 3.3, there are two parameters in this model which need to be specified to complete this model. The first.is the parameter 6M which is the value of the peak-to-single-phase mixing rate. This parameter must be specified in order to define the turbulent velocity, E/Z. The second is the parameter KM which is the proportionality constant between the fully-developed void fraction profile and the mass velocity -141- distribution. This parameter is important for determining the fully- developed distributions of the mass, energy and momentum. Appropriate values for these two parameters need to be specified to complete the formulation of the two-phase mixing model. Faya [29],in an investigation of mixing parameters for a two-phase mixture model (i.e. not a two-fluid model), recommends the following values for these two parameters: 0 m 5 (5.la) KM = 1.4 (5.1b) The selection of these values is based on numerical comparisons of the CANAL computer code [29] with experimental measurements. The magnitude of these parameters seems to be appropriate for BWR conditions [27, 32]. Since these two parameters have the same physical meaning in the THERMIT two-phase mixing model, the values given in Equation 5.1 are used for the reference case. It should be noted that while the value of KM may in fact be fairly constant, 6M is expected to depend on the flow conditions. As the flow rate increases, the value for 0M has been experimentally observed to decrease [27]. H nce', the assumption of using a constant value for 6M may not be appropriate if the flow rate rapidly changes. However, as will be discussed in Section 5.2, the sensitivity of the predicted results to variations in 0M is found to be small. Therefore, the assumption of constant 60 , which simplifies the two-phase mixing model, does not appear to adversely affect the predicted results. -142- One of the difficulties in developing a two-phase mixing model is that the mixing terms cannot be directly compared to experimental measurements on a local basis. The problem is that these terms have not (and probably cannot) be measured. Rather, isokinetic measurements of the exit flow and enthalpy distributions in rod bundles have traditionally been used to infer the appropriateness of the above terms. These exit distributions reflect the integration of the mixing effects along the entire test section length. Appropriate modeling of the transport mechanisms due to mixing is required in order to calculate these exit distributions. Consequently, the two-phase mixing model can be assessed with measurements of this type. The validation and assessment of the two-phase mixing model has been performed using rod-bundle measurements. In particular, three sets of experimental measurements have been examined in this study: the GE 9 rod bundle tests [16], the Ispra 16 rod BWR tests [37] and the Ispra 16 rod PWR tests [38]. In each of these tests, a large number of exit mass velocity and enthalpy distributions have been isokinetically measured. Test conditions for these tests are given in Table 5.1. The first two tests have been run at a pressure of approximately 6.9 MPa in geometries typical of a BWR, while the third test has been run at 16.0 MPa with the geometry being typical of a PWR. Hence, these tests cover the expected conditions for subchannel applications. In the isokinetic measuring technique, flow samples are extracted from the various subchannels (defined on a coolant centered basis) at the test section exit. The sample flow rate is adjusted to match the pressure distribution that exists when the sampling device is not present. In this -143- TABLE 5.1 Test Conditions for Rod-Bundle Experiments G.E. 9-Rod Ispra 16-Rod BWR Ispra 16-Rod PWR P (MPa) 6.9 7.0 16.0 G (kg/m s) 650 to 2200 1000 to 2000 2500 to 3500 q" (MW/M2 ) 0.71 to 2.1 0.12 to 0.77 0.07 to 0.11 Aisub (kJ/kg)j 67 to 525 30 to 180 250 to 400 xout % 3 to 22 2 to 31 -20 to 20 D (mm) 12.1 13.3 10.7 e Length (m) 1.83 3.66 3.66 Spacer Type Pin Grid Grid Radial Power Uniform and Uniform Uniform Distribution Non-Uniform Bom. WMM -144- way, the sample flow rate should equal the actual flow rate for the subchannel. The enthalpy and flow rate of the sampled subchannel are then measured by calorimetry to complete the experimental procedure. All characteristic subchannel types in the rod-bundle should theoretically be simultaneously sampled in order to insure that mass and energy are conserved. In practice, however, thorough sampling is not always done. In certain tests not all subchannel types are sampled. In other tests, all subchannel types are sampled, but not simultaneously. Consequently, significant mass and energy conservation errors are typi- cally found in measurements of this type. These errors present a problem when comparing the measurements to the code predictions since the predictions always conserve mass. Consequently, perfect agreement between the code predictions and measurements cannot and should not be expected. Nevertheless the assessment of the two-phase mixing model has been accomplished through comparisons of the measured and predicted exit mass velocity and quality (or enthalpy) distributions for selected cases from the above test sets. Since these comparisons have been made subsequent to the one-dimensional void fraction comparisons discussed in Chapter 4, it has been assumed that both the subcooled boiling model and the interfacial momentum exchange model are appropriate for these cases. If this assumption is made, then the bundle exit distributions can be used to assess the two-phase mixing model directly. The general procedure for analyzing these cases can be described as follows. First the geometrical details of the test section, as provided by the experimental reports, are used to define coolant-centered sub- channels in the test section. The reason for using a coolant-centered -145- approach is that the isokinetic measurement technique samples subchannels defined on this bais. Next, the measured inlet mass flow rate, assumed to be radially uniform at the test section inlet, and the measured inlet temperature are used to define the inlet boundary condition. The measured outlet pressure, also assumed to be radially uniform, is used as the exit boundary condition. Finally the measured power level is specified to determine the wall-to-coolant heat flux boundary condition. With these boundary conditons, the THERMIT calculations are performed until a steady-state solution is obtained. The steady-state results can then be compared with the experimental values and are discussed in this chapter. 5.2 GE 9 Rod Bundle Test 5.2.1 Test Description The mass velocity and quality distributions at the outlet of a 9-rod electrically heated test section have been measured. Both the two-phase flow conditons and geometry of the test section are similar to those found in BWR rod bundles. The pressure for these cases is 6.9 MPa. The average exit quality ranges from 3% to 22% and mass velocities ranging from 650 to 2200 kg/m2-s have been used. These conditions are in the range-of operating BWR conditions. The rod diameter and rod-to-rod pitch also closely resemble those in a BWR rod bundle. A detailed cross sectional view of the 9 rod bundle is illustrated in Figure 5.1. While the measurements have been taken primarily for two-phase flow conditions, some measurements have been made for single-phase conditions. These single-phase cases (cases lB to 1E in GE notation) have been run using isothermal conditions. For the two-phase cases both radially -146- 1 2 2' 3 2' 3 1' 2 Center Subchan 2 3 3 2 nel Edge 2 Subchannal 2 Corner Subchannel Geoemtrical Details Rod Diameter 14.478 Rod-Rod Gap 4.267 Rod-Wall Gap 3.429 Radius of Corner 10.2 Heated Length 1829. mm mm mm mm mm Figure 5.1 Cross Section View of G.E. 9 Rod Bundle Used in Mass Velocity and Enthalpy Measurements . mma I -147- uniform (cases 2B to 2E) and radially non-uniform (cases 3B to 3E) power distributions have been utilized. For each of the two-phase cases the axial power distribution is uniform. The radial non-uniform power distribution is given in Figure 5.2. Before discussing the comparisons, it is instructive to review the important characteristics of these measurments. As seen in Figure 5.1, there are three distinct types of subchannels; namely the corner, edge, and center subchannels, Measurements have been made in each of these subchannel. types and the results will be discussed in terms of the subchannel types. For the isothermal cases only mass velocity measurements have been made. These data are useful for assessing the single-phase characteris- tics of the two-phase mixing model. Since the two-phase mixing rate is a function of the single-phase value, it is important that an appropriate single-phase mixing rate be used. These measurements allow this aspect of the model to be assessed. For the uniformly heated cases both mass velocity and quality measure- ments have been made. The most significant phenomenon observed in these measurments is that the quality in the corner subchannel is much lower than the bundle average. This behavior occurs in spite of the fact that the power-to-flow ratio is highest for the corner channel. This observa- tion indicates that vapor is transported preferentially away from the corner subchannel to the more open (central) subchannels. Other than this paculiar behavior in the corner channel, the quality measurements in the other subchannels closely follow the bundle average behavior. The center subchnnel always has the highest quality and is slightly above the -148- Hot Corner Hot Edge 1 '2 9 1.305 gm -1.158 -1.001 3 ( 7 - 1.156 1.047 -C- 0.82 I 6 -- 1.005-- 0.815 ~~~ 0.67 to 6' 8 9 -- 4 Cold Edge 5 Cold Corner Figure 5.2 Radial Peaking Factors for Non-Uniformly Heated Case Hot Center ---- W=moro -149- bundle average. The edge subchannel is usually at or slightly below the average quality. In the non-uniformly heated cases, mass velocity and quality measurements have been made in five different subchannels (see Figure 5.2). Actually, for a non-uniform power distribution, there are ten distinct subchannel types. Sampling has been done in only five of these ten subchannel types because of experimental difficulties [51]. Consequently mass and energy balances could not be evaluated. The highest quality is found in the hot center subchannel while the cold side sub- channel shows the lowest quality. It is interesting to note that, although the hot corner subchannel has a higher than average quality, its quality is less than that in the hot center subchannel. The rod-bundle used in this experiment did not contain grid spacers. Rather, spacer pins were used to prevent rod motion. Frictional losses due to these pins have been reported [36]. These losses are important in determining the flow distribution and have been included in the THERMIT hydraulic modeling. A thorough error analysis has not been performed for these measure- ments. Consequently it is difficult to judge the quality of the data. However, there are a few relevent points that should be considered. The first is that sampling of the subchannels was not done simultaneously. Consequently, accurate mass balances were not always obtained. In fact, continuity errors as large as 5% have been reported. A second point is that the repeatability of the measurements was not reported. Consequently, it is difficult to quantify the error of a single measurement. Finally, for the non-uniform cases, only five subchannels were sampled of the ten - 150- characteristic types. Hence, mass and energy errors could not even be evaluated for these cases. Furthermore, only a limited number of non-uniform power cases (4 to be exact) were made making it difficult to identify any trends in the data. Estimates of errors in individual measurements have been made. Errors in individual mass velocity measurements are estimated to be 3%, while errors in the quality measurements are estimated to be 0.02 in quality [36]. It should be noted, however, that, due to the difference in flow areas of the various subchannels, a 3% error in the center subchannel velocity has a much larger effect on the total continuity error than a 3% error in the corner subchannel. Hence, the continuity error is very sensitive to errors in the center subchannel mass velocity measurement. This point must be considered when evaluating the results of the experimental comparisons. 5.2.2 Single-Phase Comparisons The comparisons of the single-phase mass velocity measurements with the predictions of THERMIT have been found to be in overall good agree- ment. Tabulated mass velocity comparisons are presented in Appendix B with Figure 5.3 graphically illustrating the good agreement between data and predictions. It is seen that the mass velocity in the center and edge subchannels is well-predicted by THERMIT over the entire range of average mass velocities. The mass velocity in the corner subchannel is satisfactorily predicted although some minor deviations are observed at low and higher average mass velocities. These deviations are not too significant considering the estimated error in the measurements. Hence, the two-phase mixing model reduces to the single-phase limit as expected. r-151- o Data -- THERMIT CENTER 0 EDGE CORNER 000 I 1000 1500 2000 2 Gav (kg/rn .s) 2500 Figure 5.3: Comparison of Measured and Predicted Mass Velocities for G.E. Isothermal Tests 1.5 -'- 1.0 I- 03 CD N CD 0.5 - 0.0 500 3000 II I -152- Furthermore, the single-phase mixing rate seems to be appropriate for these cases. 5.2.3 Uniformly-Heated Cases For the uniformly heated cases, comparisons have been made for a wide spectrum of mass velocities and qualities. Test conditions along with the measured and predicted exit mass velocity and quality distribu- tions are tabulated in Appendix B. It should also be noted that for a few cases THERMIT predictions have been made without using the two-phase mixing model. These predictions demonstrate the importance of the two- phase mixing. On the whole, the agreement between the measured and predicted values is rather good. The quality distributions as a function of the bundle average quality are illustrated in Figure 5.4 to 5.6. For the corner subchannel it is seen that the quality is significantly less than the bundle average over the entire quality range. The predicted values follow the trend in the data very well, although a slight underprediction of the data is seen at the highest quality. Without the two-phase mixing model, the quality is greatly overpredicted and these results are in very poor agreement with the data. For the edge subchannel the quality is also satisfactorily predicted by THERMIT. The majority of the data are slightly below the bundle average and this trend is predicted by THERMIT. However, in the highest quality case the predictions tend to be less than the measured value. For the center subchannel, it is seen that the measured quality is consistently greater than the average. THERMIT successfully predicts this tren and the results are in very good agreement with the measured values. -153- 0.5 OAA 0.4 03 0.2 2 G = 725 kg/m 's 0010 Data2 G =1450 kg/m s Data 0.1 ---. . THERMIT 00..0. THERMIT-No Mixing X, ave 0.00 0. .1 0.2 0.3 0.4 0.5 BUNDLE AVERAGE QUALITY Figure 5.4: Comparison of Measured and Predicted Exit Quality in Corner Subchannel for G.E. Uniformly Heated Cases -154- 0.5 0.4 / 0.3 0.2 -G = 725 kg/m -as rd Data2 -G = 1 4 5 0 kg/m -s Data - - -THERMIT 0.1 - -THERMIT-No Mixing 2 ave 000 0. 0.1 0.2 0.3 0.4 0.5 BUNDLE AVERAGE QUALITY Figure 5.5: Comparison of Measured and Predicted Exit Quality in Edge Subchannel for G.E. Uniformly Heated Cases -155- 0.5 0.4 E-11 0.3 - U 4 - 2 S 0G = 25 kg/m -s 2 Data/7G l45Okg/rn s - - -THERMIT 0.1 - -"THERMIT-No Mixing 3 ave 0.0 0. 0.1 0.2 0.3 0.4 0.5 BUNDLE AVERAGE QUALITY Figure 5.6: Comparison of Measured and Predicted Exit Quality in Center Subchanael for G.E. Uniformly Heated Cases -156- Without the two-phase mixing model, the predictions at high qualities (> 15%) tend to deviate significantly from the data. Overall, it is seen that the quality distribution is well-predicted by THERMIT. The mass velocity comparisons are illustrated in Figure 5.7. The predicted mass velocity values for the center and edge subchannels are seen to be in fairly good agreement with the data. The predicted corner mass velocity values are found to be in satisfactory agreement although a few points are underpredicted by as much as 10%. However, as seen in Figure 5.7 if the two-phase mixing effects are not included, then substantial underprediction of the data is found. Compared with these predictions, the results for the corner subchannel are seen to be quite good. Furthermore, in view of the uncertainties in the measurements, the predicted results seem to be satisfactory. Hence, due to the good agree- ment in both the quality and mass velocity distributions, the current formulation of the two-phase mixing model seems to be appropriate. 5.2.4 Evaluation of Mixing Parameters As discussed in Section 5.1, the two mixing parameters, 0M and KM' have been set to the recommended values of 5.0 and 1.4, respectively. The choice of these values has been based on the work of Faya [29]. While these values are probably also appropriate for THERMIT, different values ..may, in fact, be more optimal. Therefore, the sensitivity of the predic- tions to variations in 6M and KM for the 2E series cases have been studied. Values of 0g ranging from 1 to 10 along with variations of Km from 1 to 2 have been used in this study. These variations cover the range of possible values for these parameters and illustrate the sensitivity of the code predictions- -157- 800. 1 400) 4(00. 800. 1200. 1600. Measured Mass Velocity (kg/m2es) Figure 5.7: Comparison of Measured and Predicted Mass Velocities for G.E. Uniformly Heated Tests 1600. 1200. 0* 9 4-) IL) U) 4J) C) ari $4 124 -- Corner+% A +10% O Edge / / 6 Center /'G/ * Corner No Mixing -10% 0/ /0 --. A -- 3 / / / / -158- Physically, 6M controls the level of mixing for two-phase conditions. As 6M is increased, the mixing rate increases leading to a more uniform flow and quality distribution. The effect of varying 6M for the 2E cases, which represent typical cases, is illustrated in Figure 5.8 (KM = 1.4). While the variation of 6M affects all subchannels, the corner subchannel results are seen to be the most sensitive to the value of 6M' As 6M is increased from 1 to 10, the enthalpy (quality) in the corner channel decreases while the mass velocity increases. The changes in the enthalpy are relatively small. However, the mass velocity is seen to show large variations. The large variations in the corner subchannel mass do not significantly affect the mass continuity so that the mass velocities in the edge and center subchannels do not change by a large amount. The difference between 0m = 5 and 0m = 10 cases is seen to be relatively small. For these particular cases it may appear that a large value of 0M may be warranted. However, for the majority of cases (both here and in the following sections) a value of 0m = 5 seems to be satisfactory. The parameter KM has also been varied from 1.0 to 2.0 (with 0M = 5). Physically, as Km increases ti-e fully-developed void fraction profile becomes steeper and the amount of vapor diffusion will be greater. Conse- quently, more vapor will be transported from the corner to the center subchannel. This effect is illustrated in Figure 5.9 for a typical case (2E2). The results show that as KM is increased the quality in the corner subchannel decreases and the mass velocity increases. The center sub- channel shows the opposite trend while the quality in the edge subchannel is virtually unchanged. For this case, a high value of KM would tend to bring the corner subchannel results into better agreement, but will make 0 Corner o Edge t Center Data\ 0 -10- 14 N m o = 5 1M Da 'Ot itaN = m 10- G\ =4 ) m = 1~ 1400 Data M = 10 GM = O m 1600 ENTHALPY (kJ/kg) Figure 5.8: Comparison of G.E. 2E Cases with 0M Varied from 1 to 10 -159- 17009. - 2E_ 2E2 2E3 1500.e -CD > H- 1300 ..- 11(J 0 . 1200 Ir% f% -160- 0 Corner o Edge 6 Center 1500- Data U K2= 2 1300 -- KM=1 M 110i 1300 1400 1500- ENTHALPY (kJ/kg) Figure 5.9: Comparison of THERMIT Predictions for Case 2E2 with Variations in KM -1S1- the center subchannel results worse. Hence, it is not possible to improve the agreement by varying KM alone. It should also be noted that the changes in the mass velocity distri- butions are greater for 6M variations than for KM variations. Of course, 6M is varied over a larger range. The changes in the quality distributions are found to be nearly the same for either 6M or KM variations. The results of this sensitivity study indicate the importance of 8M and KM on the mass velocity and quality distributions. The masa&velocity distribution is most affected by eM. However, the 6M = 5 and 6M = 10 cases give similar results which supports the use of the 0m = 5 assumption. The quality distribution is virtually unaffected by variations in either eM or KM. Since neither the quality or mass velocity is overly sensitive to variations in KM, the use of the KM = 1.4 assumption is appropriate. 5.2.5 Non-Uniformly Heated Cases Four tests using the non-uniform radial power distribution (illustated in Figure 5.2) have been made. In these tests exit mass velocity and quality measurements have been taken in five of the ten characteristic subchannel types. Due to the limited amount of data and because mass and energy balances could not be evaluated, it is difficult to assess the two- phase mixing model with this data. In spite of the shortcoming of the measurements, comparisons between the measured and predicted distributions have been made. The results of the comparisons are listed in Appendix B. Figure 5.10 shows that the quality distribution can be satisfactorily predicted in each case. Nearly all of the predictions lie within 0.025 (in quality) of the measured value. Only in case 3E1 is the hot center subchannel -162- o Hot Corner 0 Hot Edge A Hot Center * Cold Corner 5 Cold Edge Hot Corner o.6 - (No Mixing) +0.025 / 0.5 / / -0.025 0.4 / 0 0.3 / / 0.N 0/ 0.1 .A 0.20 0. 0.1 0.2 0.3 0.4 0.5 0.6 Measured Quality Figure 5.10: Comparison of Measured and Predicted Quality for G.E. Non-Uniformly Heated Tests -163- underpredicted by more than 0.05. On the whole, and in view of the limitations of these measurements, the predicted quality distributions are found to be in good agreement with the data. Comparisons between the measured and predicted mass velocities are illustrated in Figure 5.11. It is seen that a majority of the data is predicted within 10%. In particular, the hot corner, hot edge and cold edge velocities are well-predicted. However, the agreement for the hot center and cold corner is only fair. As discussed above, the mass flow error for these cases could not be evaluated so that is is difficult to assess the accuracy of the measurements. Hence, although the agreement in the mass velocities is not as good as expected, the comparisons indicate that, overall, THERMIT is satisfactorily predicting the exit mass velocity and quality distributions for the non-uniformly heated cases. 5.3 Ispra BWR -tests 5.3.1 Test Description Mass velocity and quality distribution have been measured in a 16-rod electrically heated test section at Ispra [37].. As in the case of the G.E. tests, both the two-phase flow conditions as well as the geometrical characteristics of the test section closely matched those found in a BWR. The nominal pressure was 6.9 MPa, exit qualities ranged ... 2 from 2% to 31% and mass velocities of 1000, 1500, and 2000 kg/msec have been used. The geometrical details of the test section are illustrated in Figure 5.12 and it is seen that the rod diameter and rod-to-rod pitch are similar to those found in a BWR bundle. Approximately 225 test cases have been reported in reference 37. For -164- o Hot Corner U Hot Edge A Hot Center * Cold Corner +10% 3 Cold Edge 20.00 / / / -10% 1600 / / / >10 r44 /0/ 1200 202 /0 / r) P4) 049 800 7 Measured Mass Velocity (kg/m2S) Figure 5.11: Comparison of measured and Predicted Mass Velocities for G.E. Non-Uniformly Heated Tests -165- Geometrical Details Rod Diameter 15 mm Rod-Rod Clearance 4.5 mm Rod-Wall Clearance 3.37 mm Corner Radius 104.1 mm Heated Length 3660 mm Figure 5.12 Cross Sectional View of Ispra BWR Test Section 2 1 3 1 2 1' 5 ' 6 5 :1 31 6 4 : 6 ' 3 1' 5 ' 6 : 5 1 2' 1 3 ; -166- each case, simultaneous isokinetic sampling at the test section exit in four of the six distinct subchannel types have been performed. The sampled subchannels along with the numbering scheme are indicated in Figure 5.12. A uniform power distribution has been used in all cases. Three important differences between the Ispra test section and G.E. test section can be identified. First, the Ispra test section was 3.66m long, which is the nominal BWR bundle length, while the G.E. bundle was only 1.83m long. Since the two test sections were operated at approxi- mately the same mass velocity, inlet temperatures and outlet qualities, the heat flux in the Ispra bundieis much less than that used in the G.E. bundle. The second difference is the number of rods in each bundle. The Ispra bundle has 16 rods while the G E. bundle has 9 rods. With 16 rods there are six distinct subchannel types rather than three types in a 9 rod bundle. The third difference is that the Ispra bundle used grid spacers while the G.E. bundles used spacer pins. The loss coefficients associated with the grids are much larger than those for the pins and are very important for proper modeling of the hydraulic character- istics of the bundle. Unfortunately, these loss coefficients were not measured and only estimated values have been reported. Since additional information isnot available for the grid loss coefficients, the estimated values have been used in the subsequent analyses. The general trends in the quality measurements show many of the same characteristics found in the G.E. tests. For example, the quality in the corner subchannel is below the average even though the mass velocity in this subchannel is the lowest. The quality in the other subchannels closely follows the bundle average behavior, again as seen in the G.E. -167- tests. In the mass velocity measurements one main difference between the G.E. and Ispra tests can be seen. The difference is that, while in both tests the corner mass velocity is below the average, in the Ispra tests it is much lower than average. This result is attributed to the use of grid spacers which have a higher loss coefficient in the corner subchannel. Besides this difference, the mass velocities in the other subchannels show the same qualitative behavior in both test sections. The estimated errors in the individual mass velocity and quality measurements are on the order to 3% [37]. However, since only four of the six characteristic subchannels were sampled, a mass balance could not be calculated. Attempts to estimate the continuity error indicate that errors as high as 8% may be possible [37]. This error is rather large, but is compensated by the fact that numerous measurements were made for the same conditions. Consequently, the data show significant spread, bit this would be expected for measurements of this type. With these measurements, evaluation of the two-phase mixing model is easier since the trends in the data are much clearer. 5.3.2 Results For the Ispra BWR tests a total of eleven representative cases have been simulated with THERMIT. These cases cover a range of qualities at each of the three mass velocity values (i.e. 1000, 1500, and 2000 kg/m2 -sec). It should be noted that for all of these cases, 0M is set to 5.0 while KM is set to 1.4. Also as in the case of the G.E. tests, uniform inlet velocity and uniform outlet pressure boundary conditions have been used. -168- Comparisons between the measured and predicted quality and mass velocity distributions have been made. In Appendix B, the test conditions, measured and predicted values as well as the predictions using the original version of THERMIT (i.e. no mixing) are listed. Through these comparisons, the two-phase mixing has been assessed. On the whole, the quality distribution is well-predicted for each subchannel. Figures 5.13 to 5.16 show the comparisons between the measured and predicted quality values for the various subchannels. In each figure, the shaded areas represents the spread in the actual measured values. For the corner subchannel, the predicted values are in very good agreement with the measured values (see Figure 5.13). If the mixing model is not included, then the quality is significantly overpredicted. Hence, as in the G.E. cases, the inclusion of the effects of mixing are essential for accurate predictions of the corner subchannel quality. It should also be noted that COBRA IIIC cannot predict the correct quality behavior of the corner subchannel [53]. In the edge subchannel, good agreement between the measured and predicted qualities is also found (see Figure 5.14). For both of the center subchannels, THERMIT tends to predict qualitites which are on the high side of the data at high bundle average quality (see Figure 5.15 and 5.16). However, on the whole, the quality distribution is very well predicted by THERMIT for these cases. The mass velocity comparisons are illustrated in Figure 5.17 and 5.18. Figure 5.17 shows the comparisons for the edge and corner subchannel results. It is seen that the THERMIT predictions for the corner subchannel are consistently high. This result may be due to one of two reasons. -169- 0.4 G = 1000 kg/m2DS G = 1500 kg/m2 9S 0.3- G = 2000 kg/m2 e X 2 .2 THERMIT 0.1 0.0 - 0. 0.1 0.2 0.3 0.4 BUNDLE AVERAGE QUALITY Figure 5.13: Comparison of Measured and Predicted Exit Quality of Subchannel 2 - Ispra BWR Tests -170- 0.4 G = 1000 kg/m2. S G = 1500 kg/m2-s - G = 2000 kg/m2s2/ 0.3/ x/ $ 0.2 ./ /THERMIT 0.1 0.0- 0. 0.1 0.2 0.3 0.4 BUNDLE AVERAGE QUALITY Figure 5.14: Comparison of Measured and Predicted Quality for Subchannel 1 versus Bundle Average Quality Ispra BWR Tests -171- 0.4 OA' 2 frYG =lOO00kg/m-s // G= 1SO0 kg/rn-s 2 ~ G=2000 kg/m -s 0.3G THERMIT x 5 0.2 0.1 -- 0.0 / 0. 0.1 0.2 0.3 0.4 BUNDLE AVERAGE QUALITY Figure 5.15: Compaison of Measured and Predicted Quality for Subchannel 5 versus Bundle Average Quality Ispra BWR Tests -172- 0.4 +v/ G = 1000 gms G = 1500 kg/m2es 2 N/ 0.3 - G=2000 kg/m2-s x4 THERMIT 0.2 0.1 0.0 0. 0.1 0.2 0.3 0.4 BUNDLE AVERAGE QUALITY Figure 5.16: Comparison of Measured and Predicted Quality for Subchannel 4 versus Bundle Average Quality Ispra BWR Tests- -173- 1.2 1.0 0.8 0.6 0 Figure 5.17: 1.2 1.0 0.8 0.6 0. 0.1 0.2 0.3 BUNDLE AVERAGE QUALITY Comparison of Measured and Predicted Mass Velocities for Subchannels 1 and 2 - Ispra BWR Tests 0.1 0.2 BUNDLE AVERAGE QUALITY 0.3 Figure 5.18: Comparison of Measured and Predicted Mass Velocities for Subchannels 4 and 5 - Ispra BWR Tests a, CD N *rI CD Subchannel 1 S h THERMIT Subchannel 2 U) CD N *rI CD Subchannel 4 THERMIT Subchannel 5 I .a -174- First, the mixing rate may be too high which, as seen in Section 5.2.4, leads to higher mass velocity values. However, these cases have also been run without including the effects of mixing and the mass velocity is still overpredicted. The second reason may be that the grid coeffi- cients are not correct. As discussed by Lahey [25], proper modeling of the grids is essential for predicting the correct flow distribution. In fact experimental measurements have shown that the trends in flow distribution can be significantly different in a rod-bundle with grids compared to one without [52]. Since THERMIT models-the grids as simple pressure losses, it may be anticipated that the deviations in the corner subchannel mass velocity are due to improper grid modeling. The mass velocity in the edge subchannel is seen to be slightly underpredicted. However, the agreement is satisfactory. The deviations in the predictions may again be due to the grid modeling. The mass velocity comparisons for the two center subchannels are illustrated in Figure 5.18. It is seen that the predictions in both subchannels. are nearly the same. The agreement for subehannel 5 is very good in all cases. The predictions for subchannel 4 are also found to be quite satisfactory. In summary, it has been shown that the predicted quality distribution is in overall good agreement with the dati for these cases. Furthermore, the lower than average behavior of the corner subchannel is well predicted by THERMIT. This trend cannot be predicted by COBRA IIIC [531. The mass velocity distribution is found to be satisfactorily predicted in view of the uncertainties in the grid spacer modeling. The corner subchannel shows the largest deviations, but this subchannel will also be the most sensitive I - .- . .1 . - -175- to the grid coefficient. 5.4 Ispra 16 Rod PWR Cases 5.4.1 Test Description Until recently, there have been no consistent experiments to determine the steady-state mass velocity and quality distributions in typical PWR geometry and for typical PWR operating conditions. The usual problem with most previous work is that only a few of the charac- teristic subchannels were sampled [54]. However, recent experiments performed at Ispra [38] represent the first consistent effort to determine the mass velocity and quality distributions for PWR conditions. These experiments have been conducted in a 16 rod electrically heated test section. Rod diameters and rod-to-rod spacing closely resemble those found in a PWR rod-bundle. Details of this test section are given in Fig. 5.19. The operating conditions for these tests also simulated PWR conditions. The nominal pressure was 16.0 MPa, the mass velocity ranged from 2500 to 3000kg/m2.sec, and the exit quality ranged from -20% to +20%. The high end of the quality range served to simulate two-phase flow conditions which might occur in transients. As in the case of the Ispra BWR tests, grid spacers have been used to maintain the rod spacing. Unfortunately, details of the spacer design are not available. Consequently, the grid loss coefficient could only be estimated. However, as indicated in the last section, the mass velocity distribution is strongly dependent on the grid modeling. Therefore, some error in the predictions will most likely result from the lack of detailed -176- Geometrical Details Rod Diameter 10.75 mm Rod-Rod Clearance 3.55 mm Rod-Wall Clearance 2.93 mm Corner Secant 8.3 mm Heated Length 3660 mm Figure 5.19 Cross Sectional View of Ispra PWR Test Section S3 1 ' 56' 6 - 6 1 4 1 6 '3 5 1 6 1 5 : 1 2'1 Q 31 1 ;2 -177- grid modeling, The experimental measurements indicate that the subchannel exit mass velocity and quality distributions closely follow the bundle average values. Unlike the BWR cases, the quality in the corner subchannel is not significantly below the average. This result shows that the pressure plays an important role in determining the flow and enthalpy distributions. The higher mixing rate in the BWR bundle seems to be a direct result of 33the larger specific volume of the vapor (for BWR pv = 36.1 kg/rn3 ; for PWR P =85.4 kg/mn). The most likely reason the typical BWR corner subchannel behavior is not found in the PWR cases is that the increased pressure alters the flow regime. The slug-annular transition quality varies from approximately 0.1 at 7.0 MPa to 0.21 at 16.0 MiPa. Physically, this result indicates that a bubbly flow can be maintained up to higher qualities at the higher pressure. In terms of the vapor diffusion rate, if annular flow does not occur then it is more difficult to transport the vapor from the corner subchannel to the center subchannel. Hence, the vapor diffusion rate is lower at the high pressure. Error estimates for these measurements are comparable to those found for the Ispra BWR tests. That is, individual flow and enthalpy errors are on the order to 3% [38]. However, in the PWR tests five of the six characteristic subchannels have been sampled, which allows for more accurate assessment of the continuity and energy errors. 5.4.2 Results Comparisons between the measured and predicted mass velocities and qualities have been made for the PWR cases. As in the BWR cases, the -178- values for the mixing parameters bM and I are 5.0 and 1.4 respectively. Using these values, a few representative cases have been analyzed. however, due to the proprietary nature of these measurements, only graphical illustrations of the comparisons will be given. Overall, the quality predictions of THERMIT agree quite well with the measurements. Figures 5.20 to 5.24 illustrate these comparisons. As discussed above, it is seen that the quality in each subchannel is very near the bundle average value, even for the corner subchannel. THERMIT is able to predict this trend rather well. Good agreement between the measured and predicted values is also found for each of the subchannels. COBRA IIIC/MIT also predicts the correct trend in the data. This indicates that the effects of mixing are not as important for these cases. The mass velocity measurements are found to be satisfactorily predicted by THERMIT. Figures 5.25 and 5.26 show the comparison between the measured and predicted values. Except for the corner subchannel, the predicted mass velocities are always within 10% of the data. The corner subchannel shows a significant underprediction for subcooled conditions. However, since the qualities are in such good agreement, these deviations are probably due to the grid modeling. On the whole, and in view of the approximate grid modeling, the mass velocity distribution is satisfactorily predicted by THERMIT. The good overall agreement shows that the mechanistic formulation of the two-phase mixing model is also appropriate for PWR conditions. This is an important result since it verifies the assumption that mechanistic-models can be extended beyond their data base range. -179- THERMIT COBRA IIIC/MIT +10 XI -10 -10 Subchannel 1 -- -/ N'K +10 x ave Figure 5.20: Comparison of Measured and Predicted Exit Quality for Subchannel 1 versus Bundle Average Quality - Ispra PWR Tests (COBRA IIIC/MIT Results form Reference 38) I I -180- -THERMITSubchannel. 2 --..- COBRA-IIIC/MIT +20 x 2 +10 - -10 +10 X ave -10- Figure 5.21: Comparison of Measured and Predicted Exit Quality for Subchannel 2 versus Bundle Average Quality - Ispra PWR Tests (COBRA IIIC/MIT Results from Reference 38) -181- THERMIT - COBRA IIIC/MIT +201 x3 +101 - - -19mz-Ih %L 7: _ _ -10 / / / N / / / / -10 / N / N / '7 / / Subchannel 3 .. / M/ +10 x ave Figure 5.22: Comparison of Measured and Predicted Exit Quality for Subchannel 3 versus Bundle Average Quality - Ispra PWR Tests (COBRA IIIC/MIT Results from Reference 38) I I I -182- - THERMIT - -- COBRA IIIC/MIT +2 0| x4 +14~ N -10 -1C - I +10 x ave Figure 5.23: Comparison of Measured and Predicted Exit Quality for Subchannel 4 versus Bundle Average Quality - Ispra PWR Tests (COBRA IIIC/MIT Results form Reference 38) Subehannel 4 -/ a -183- THERMIT - - COBRA IIIC/MIT x5 +10 /j -10 -10 I --/ Subchannel 5/ -- N +10 x ave -- Figure 5.24: Comparison of Measured and Predicted Exit Quality for Subchannel 5 versus Bundle Average Quality - Ispra PWR Tests (COBRA IIIC/MIT Results from Reference 38) I -184- -10 0 +10 +20 Bundle Average Quality Subchannel 2 -7 THEPM IT I I -10 0 Bundle Average Quality +10 +20 Figure 5.25: Comparison of Measured and Predicted Mass Velocities for Ispra PWR Tests 1.10 1.05 1 CD N Ci 1.00 1 0.95 L 0.90 L -2( Subchannel 1 THERMIT . . - I. - .- .' 0 1.05 1.00 0.95(5N .rI(5 0.90 L- 0.85 -20 -. M -185- Subchannel 3 1.0 THERMIT 0.95 0.901 -20 -10 0 +10 +20 Bundle Average Quality 1.05 - Subchannel 4 THERMIT a1.00 - N rq 0 .9 5 - 0.901 -20 -10 0 +10 +20 Bundle Average Quality Subchannel 5 THERMIT 1.05- 1.00 0.95' -20 -10 0 +10 +20 Bundle Average Quality Figure 5.26: Comparison of Measured and Predicted Mass Velocities for Ispra PWR Tests -186- Consequently, in addition to assessing the two-phase mixing model for PWR conditions, this study has shown that the two-phase mixing model may be used for a wide range of conditions. 5.5 Conclusions The above discussion of the two-phase mixing model assessment has shown that the current formulation of this model is appropriate for a wide range of conditions. Both BWR and PWR subchannel quality distribu- tions have been successfully predicted. This good agreement includes predicting the lower than average quality in the corner subchannel for BWR conditions as well as predicting the more or less uniform behavior for PWR conditions. It should be emphasized that the same mixing parameters have been used for both types of conditions. Thus, the formula- tion of the two-phase mixing model seems to be well-founded. A review of the predicted mass velocity distributions indicates that although, on the whole, the predictions are good, the mass velocity in the corner subchannel for the Ispra tests is usually not well predicted. In the BWR tests the corner mass velocity is overpredicted, while in the PWR tests, it is underpredicted. These deviations cannot be explained in terms of the mixing model so it has been assumed that improper grid modeling may be the cause of the errors. However, it should be reiterated that all mass velocities are usually predicted to within 10% of the measurements which seems to be satisfactory considering the inherent experimental errors. The variations in the mixing parameters 6M and KM have illustrated the dependence of the flow and enthalpy distributions to these parameters. Increases in either 6M or KM tend to increase the quality and decrease the -187- mass velocity in the center subchannels, while decreasing the quality and increasing the mass velocity in the corner subchannels. The pre- dicted distributions are not overly sensitive to these parameters. Therefore, in view of the good agreement over the wide range of conditions studied here, it may be concluded that the two-phase mixing model is appropriate for both BWR and PWR conditions. -188- 6.0 THE HEAT TRANSFER MODEL 6.1 Introduction Heat transfer from the fuel rod to the coolant must be carefully modeled. This heat transfer not only represents the energy input to the coolant, but also serves as a boundary condition for the fuel rod temperature calculations. In the energy conservation equation, the source term, Q%, is that due to the rod heat transfer. On the other hand, the solution of the time-dependent heat conduction equation within the fuel rod involves the rod heat transfer as a boundary condition for the temperature profile in the fuel rod. A fundamental difficulty arises when seeking to model heat transfer phenomena over a wide range of conditions. Typically, heat transfer modeling involves application of empirical correlations to determine the heat transfer coefficient, H, which relates the heat flux q" to the temperature difference between the wall and fluid: q"1= H(T - T ) (6.1) A typical boiling curve, seen in Fig. 6.1, illustrates that the relation- ship between the heat flux and the temperature difference is very complicated so that the heat transfer coefficient is not a simple function. Nevertheless, for a limited range of application, heat transfer coeffi- cients can be found and correlated as functions of the flow conditions. Within the limits for which a correlation has been developed, the accuracy of the correlation is generally found to be satisfactory. -189- Single Phase Transition Film Vapor [Boiling Boiling q'' CHF -Nucleate Boiling qmsfb Single Phase Liquid TCHF WALL TEMPERATURE Figure 6.1: Typical Boiling Curve -190- Therefore, for a given set of flow conditions accurate heat transfer analysis can be performed provided an appropriate heat transfer coeffi- cient is chosen. However, for the wide range of flow conditions which may be anticipated for LWR transients, no one correlation could possibly be accurate over the entire range of conditions. Nevertheless, THERMIT should be equipped to analyze the various situations. Consequently, this requirement means that the heat transfer model will actually contain a number of correlations which collectively should cover the range of anticipated conditions. However, since these correlations will be used in a computer code, it is necessary that a logical scheme be developed to select the appropriate coefficient for a given set of conditioni. Hence, the heat transfer model, also referred to as the heat transfer package, consists of a number of heat transfer correlations and a logic system which dictates the choice for the application. A heat transfer model, proposed by Bjornard [20], has been originally incorporated into THERMIT. This model, referred to as the BEEST (best estimate) heat transfer package had been developed for PWR blowdown heat transfer analysis. As such, correlations covering high and low pressures, pre-CHF and post-CEF conditions as well as a wide range of anticipated mass velocities have been included. The basic idea behind this model is that a complete boiling curve can be constructed for any location of interest. Then, based on the local flow conditions an appropriate heat transfer coefficient can be selected. In order to have this ability to construct a complete boiling curve, correlations for all of the anticipated heat transfer regimes must be -191- included. The BEEST model considers five major heat transfer regimes which include forced convection to liquid, nucleate boiling, transition boiling, film boiling and forced convection to the vapor. Within each of these regimes, except transition boiling, further division can be made depending on the flow conditions so that a total of ten regimes are actually modeled (see Table 6.1 for summary of heat transfer correlations). The other important feature of the BEEST model is its heat transfer regime selection system. This system uses local instantaneous values for the flow as input to the selection scheme. As illustrated in Fig. 6.2, there are four main checks used in the selection system. First, if the quality is greater than 0.99, then single-phase vapor heat transfer is assumed. Second, if the wall temperature is less than the saturation temperature, no boiling can occur so that single-phase liquid heat transfer is assumed. For conditions between the first two checks, that is two-phase flow conditions, the two remaining checks compare the wall temperature to the minimum stable film boiling temperature, Tmsfb, and the critical heat flux temperature, TCHF. If the wall temperature is greater than T , then film boiling is assumed. If the wall temperature msfb' is less than TCHF, then nucleate boiling is assumed. For wall temperatures between Tmsfb and TCHF transition boiling is assumed. Thus, the selection system is simple and computationally efficient. However, as stated before, the BEEST package had originally been developed for PWR blowdown heat transfer analysis. As such certain simplifications were made in developing the model. These simplifications may not be appropriate for non-blowdown transients which are the primary focus of the present work. Hence, the heat transfer model has -192 TABLE 6.1 Summary of Heat Transfer CorrelatIons Regime Correlation 1. Forced convection to single-phase liquid Sieder 2. Natural convection to single-phase liquid McAdams 3. Subcooled boiling Chen 4. Nucleate boiling Chen 5. Transition boiling Interpolation between qCF and isfb 6. High P, high G film boiling Groeneveld 5.7* 7. Low P, high G film boiling Modified Dittus- Boelter* 8. Low G film boiling Modified Bromley plus either McAdams vapor or high flow film boiling* 9. Forced convection to single-phase vapor Sieder-Tate 10. Natural convection to single-phase vapor McAdams Regime Checkpoints Correlation 1. Critical Heat Flux Biasi, W-3**, CISE-4** 2. Minimum Stable Film Boiling Temperature Henry * Correlation deleted in this research ** Correlation added in this research .I -193-- Figure 6.2 CALCULATE FLOW FIELD I X > 0.99 NO T < Ts w sat __JNO BEEST Heat Transfer Logic YES YES > Tfb iYES NO T 0.99 NO Modified Heat Transfer Logic YES r1 T < T -FR w - sat NO T >T 1 YES w msfb NO T < T CEF NO TRANSITION BOILING YES SINGLE-PHASE VAPOR SINGLE-PHASE LIQUID SINGLE-PHASE VAPOR NUCLEATE BOILING -199- 6.3 Assessment The overall model as described above has been assessed by making comparisons with both steady-state and transi'nt measurements. The purpose of this effort has been to identify the reliability of the predictions of this model. Two different sets of experimental measurements have been used in this study. The first set is the Bennett [39] test measurements while the second is the G.E. transient CHF measurements in a 9-rod bundle [39]. The Bennett tests, already discussed in connection with the verification of the droplet vaporization model, are steady-state measurements in which CF occurs at locations upstream of the exit. Since wall temperature measurements were made along the entire heated length, both pre-CHF and post-CHF measurements are available. Consequently, this data is well- suited for evaluating the entire heat transfer model for steady-state conditions. The second set of measurements used in this study are the G.E. transient CHF test data. In these tests, simulated BWR fuel rods are subjected to a variety of flow transients with the time and location of the critical heat flux being measured. This data allows for the CHF predictive capability to be evaluated for transient conditions. Together, these two sets of experimental measurements permit the assessment of the main features of the heat transfer model. The steady- state results are discussed in the next Section while the transient results are discussed in Section 6.3.2. -200- 6.3.1 Steady-State Results The Bennett tests have been used to assess the THERMIT heat transfer model for steady-state conditions. In these tests, wall temperature measurements have been made along vertically heated tubes having lengths of either 3.6m or 5.5m. A wide range of conditions have been covered in these tests as summarized in Table 4.4. Both pre-CHF and post-CHF temperatures have been measured with typical data illustrated in Fig. 6.4. There are two distinguishing regimes in this figure which can be identified. The first is the pre-CHF regime which extends from the inlet to the CHF location. In this regime, the heat transfer mechanism is predominantly the nucleate boiling type so that measurements in this regime can be used to assess the nucleate boiling heat transfer correlation in THERMIT. The second regime is the post-CHF regime. As discussed in Section 4.3.3, the wall temperatures in this regime are utilized to assess the droplet vaporization model which, strictly speaking, prevents independent validation of the heat transfer correlations. However, the appropriateness of the correlations in this regime can still be assessed. Between the two regimes lies the CHF location. Due to the obvious location of this point, it is rather easy to assess the CHF correlations using this data. In view of the nature of the data, the assessment effort has been divided into three categories: (1) pre-CHF regime assessment (2) CHF location assessment (3) post-CHF regime assessment. -201- 1150 1050 - +-- II 950- 850- 750- 650- CHF Location 550 - 0 1.0 2.0 3.0 4.0 5.0 AXIAL HEIGHT (m) Figure 6.4: Typical Wall Temperature Distribution for Bennett Case 5273 -202- Comparisons between the measurements and THERMIT predictions have been made for each of these categories with the following results being found. For the pre-CHF regime, the comparisons have indicated that THERMIT consistently overpredicts the wall temperatures when using the Chen correlation [59]. This behavior is illustrated in Figures 6.5 to 6.9. It is seen that, except for the lower heat flux case, the wall temperatures are typically overpredicted by 10K. This overprediction indicates that the Chen heat transfer coefficient is too low for these cases. However, it should be noted that there is considerable scatter in the data and no error analysis is reported for these measurements. Some data points are obviously in error (for example those below the saturation temperature) and the fluctuations from one point to the next are physically too large. Two features of this experiment may account for these observations. First, the primary intent of these experiments was to measure post-CHF temperatures, so that there was little emphasis on accurately measuring the pre-CHF temperatures. Second, the thermocouples were placed on the outside of the tube. Consequently, the inside wall temperatures were not directly measured, but were calculated. These two features of the experiment may account for the relatively poor data for the pre-CHF conditions. In order to assess the Chen correlation In greater detail, the above cases have been reanalyzed using two different heat transfer correlations. For this study the Thom correlation [60] and the Jens-Lottes correlation [61] have been used. Both of these correlations are nucleate boiling correlations, while the Chen correlation is a combination of a forced convection vaporization and nucleate boiling correlation. The o Data -- Chen Thom Jens-Lottes - --- T sat -o 0o - * 0 0 0 0 0 a 0 560. 0. 1.0 2.0 3.0 4.0 Figure 6.5: AXIAL. HEIGHT (m) Comparison of Measured and Predicted Fce-CHF Wall Temperatures for Bennett Case 5332 -203- TEST CONDITIONS P = 6.9MPa 1 2q" = 0.65MW/m G = 665kg/m- Ai s = 137kJ/kg sub De = 12.6mm 580 . r- 570. 5.0 18 - - mm dlN -- IN- dmv II I I -204-- a Data Chen TEST CONDITIONS Thom p = 6.9MPa Jens-Lottes 2 tG = 1017kg/m-s Tsat G=2 --l q" = 0.79MW/m Aisub = 121.kJ/kq D = 12.6mm e 0 U - . - U '- . .a 40 0 a 0 0 I* 1.0 2.0 3.0 4.0 5.0 AXIAL HEIGHT (m) Figure 6.6: Comparison of Measured and Predicted Pre-CHF Wall Temperatures for Bennett Case 5276 580. 570. 560. 0o 0 a 0. I II -205- 0 Data Chen ----- Thom 2 G = 1356kg/m- Jens-Lottesq" = 0.9MW/m2 Tsat 4isub = 102kJ/kg D = 12.6mm e 00 00 00 0 0 1.0 2.0 3.0 AXIAL HEIGHT (M) 0 4.0 5.0 Figure 6.7: Comparison of Measured and Predicted Pre-CHF Wall Temperatures for Bennett Case 5253 TEST CONDITIONS P = 6.9MPa s 580. --- 570. Is 0 560. r- 0. I II I -206- O0 Data Chen Thom Jens-Lottes T ais = 72.kJ D = 12.6m e 00 -. 0 . 0 Oa 0- 0 0 a 0 5.00. 1.0 2.0 3.0 4.0 AXIAL HEIGHT (m) Figure 6.8: Comparison of Measured and Predicted Pre-CHF Wal1 Temperatures for Bennett Case 5394 590. 580. k-. TEST CONDITIONS P = 6,9MPa 2G = 5181kg/m-s q" 1.75MW/m2 /kg rn 134 0 570. -- 560. - ova aft-b ll" I I I i -207- o Data - Chen TEST CONDTIONS -- Thomp = 6.9M Jens-Lottes - sa-TG = 1343kg q" = 1-16M A isub 130kJ/ D = 12.G6 e 0 0 0 0 00 o 0 - 0 a S I - - . I - a - aa 1.0 2.0 3.0 4.0 AXIAL HEIGHT (i) Figure 6.9: Comparison of Measured and Predicted Pre-CHF Wall Temperatures for Bennett Case 5451 (Length = 3.66m) Pa 2 LW/m 'kg mfl 580 1- 570 560 I- 0. 5.0 I I I aI II -208- predictions using these other two correlations are also illustrated in Figures 6.5 to 6.9. It is seen that the predictions with either the Thom or Jens-Lottes correlation are in much better agreement with the data. This good agreement indicates that the primary heat transfer regime is nucleate boiling and that the nucleate boiling part of the Chen correlation is not appropriate for the high heat flux cases. As discussed by Lahey and Moody [25], the form of the nucleate boiling correlation may be given as q"1= K(T - T ) (6.3) w S where K is usually a pressure dependent constant. For the three correla- tions here, the values for m are m = 4 Jens-Lottes m = 2 Thom (6.4) m = 0.99 Chen The large value of m in the Jens-Lottes correlation explains why this correlation always predicts the lowest wall temperatures. On the other hand, the low exponent in the Chen correlation explains why this correlation predicts the greatest wall temperatures. However, at the lowest heat flux it is the Chen correlation which is in best agreement with the measurements. This fact tends to indicate that the nucleate boiling part of the Chen correlation is appropriate only for the lower heat flux cases. However, the lower heat flux values are typical of average BWR conditions which are the conditions of interest for this -209- research. One final observation can be made for these comparisons. While the Thom and Jens-Lottes correlations show linear temperature profiles, the Chen correlation temperature distribution shows an initial increase to a maximum followed by a decrease. This behavior illustrates the influence of the forced convection vaporization. As the quality is increased, the liquid film on the wall decreases and the more efficient forced convection vaporization heat transfer mechanism begins to dominate. The heat transfer coefficient for this regime is higher than the nucleate boiling coefficient so that the wall temperatures decrease. The predicted trend is seen very clearly in Fi3. 6.8. It should be noted that, even with the large scatter in the data, the measured values also tend to decrease at the higher qualities (i.e., at higher axial heights). Only the Chen correlation is able to predict this trend. In summary, it has been shown that for a majority of the cases, the Chen correlation overpredicts the measured wall temperatures. Both the Thom and Jens-Lottes correlations can satisfactorily predict the wall temperature distribution for these cases. However, only the Chen correla- tion can predict the decrease in the wall temperatures due to the forced convection vaporization heat transfer regime. Also, due to the poor quality of these measurements, it is difficult to assess any correlation with much certainty. Finally it should be noted that the Chen correlation is recommended by both Lahey and Moody [25] and Collier [62]. Thus, even with the overpredictions at high heat fluxes, THERMIT should satisfactorily predict the wall temperatures for typical BWR conditions using the Chen correlation. -210- The evaluation of the CHF correlations has focused on assessing the Biasi and CISE-4 correlations. (The W-3 correlation could not be assessed with the Bennett measurements because the quality was outside the range of applicability of the W-3 correlation.) Predictions using both the Biasi and CISE-4 correlations have been compared with the observed CHF location for representative cases. In these comparisons, the location of CHF has been used as a means of comparing the measurements with the predictions. A range of mass velocities and heat fluxes have been examined in this study. Table 6.2 2lists the results of these comparisons. For the G = 665 kg/ms cases, it is seen that the Biasi correlation overpredicts the CHF location while the CISE-4 correlation underpredicts it. At the highest heat flux (case 5336), the Biasi correlation is seen to be in excellent agreement with the data. The CISE-4 correlation shows no apparent dependence on the heat flux and typically underpredicts the CHF location in terms of the boiling length by 15%. The Biasi correlation shows a maximum overpredic- tion of 7% for these cases. For the G 1000 kg/m2.s cases, the trend of the two correlatibns is 'the same as above. The Biasi overpredicts the CHF location by as much as 25% while the CISE-4 correlation underpredicts the CHF location by 8%. It is seen that in this case the CISE-4 correlation shows better agreement than the Biasi correlation. Also, there is no improvement in the predictions of the Biasi correlations for the higher heat fluxes. For the G = 1356 kg/m2.s cases, the CISE-4 correlation is seen to be in good agreement with the data, having a maximum error of 6%. The Biasi correlation.shows poor agreement since CHF is not predicted f'r these cases. -211- TABLE 6.2 CHF Comparisons for Bennett Test Cases Case No. G20q'2 Z Measured ZBiasi ZCISE-4 kg/m -s Mw/mm () 5325 5332 5336 5266 5276 5273 5245 5247 5253 5391 5394 5396 665 665 665 1004 1017 1017 1356 1356 1356 5181 5181 5181 0.572 0.649 0.822 0.704 0.795 0.922 0.804 0.845 0.903 1.66 1.76 1.84 5.16 4.39 3.56 5.00 4.39 3.56 4.95 4.70 4.39 5.00 4.55 4.17 5.44 4.69 3.56 No CHF 5.31 4.40 No CHF No CHF No CHF 4.43 3.91 3.57 4.41 3.81 2.92 4.62 4.06 3.42 4.77 4.60 4.12 3.59 3.16 2.91 I a I - 1 --1 __ __L __ __ _ __J I -212- 2 Finally, at G = 5181 kg/ms, it is seen that the Biasi correlation now underpredicts the CHF location. The maximum error is approximately 15%. The CISE-4 correlation significantly underpredicts the data for these cases (> 30% error) and, hence, the agreement here is poor. By collectively examining these comparisons the following trends can be seen. First, the CISE-4 correlation, while always underpredicting the CHF location, is in reasonably good agreement with the data except at high mass velocities. The Biasi correlation is found to overpredict the CHF location for all but the highest velocity cases. The agreement between the Biasi correlation and the data is seen to be best at the low and very high mass velocity cases, but is poor to fair at intermediate velocities. Hence, neither correlation shows good agreement over the entire range of conditions considered here. The poor agreement found with Biasi correlation at intermediate velocities can be understood by considering the behavior of this correlation as a function of the mass velocity and quality. Figure 6.10 illustrates the predicted critical heat flux versus quality for a number of mass velocities. It is seen that there are two distinct regimes; the low quality regime (X< 0.4) and the high quality regime. The measured data are also shown and it is seen that for the low quality regime the correlation seems to be appropriate. However, at the high qualities the measured points lie beneath their respective curves. This indicates that the Biasi is not predicting a low enough critical heat flux at the higher qualities. The reason for the overprediction of the critical heat flux is probably due to the fact that all curves go to zero at a quality of 1.0. -213- 5.0 C G = 665kg/m2 s G = 665 kl-. 2 /52 0 G = 1000kg/m -s 24.0 -t G = 1356kg/n2a SG =5181kg/m2-s N kg G=1356 3.0 m.s G = 1000 kg m m -s 2.0 a G= 5200 kg2 100 _m .s r . H1 G0=0500 k 0. 0.2 0.4 0.6 0.8 1.0 QUALITY Figure 6.10: Illustration of Mass Velocity and Quality Dependence of Biasi CHF Correlation -214- This behavior, while not incorrect, does not account for entrainment effects. That is, if the entrainment is high enough, then CHF will occur at qualities of less than 1.0. Furthermore since the entrainment rate depends on the mass velocity, with a higher the mass velocity CHF will occur earlier. This would explain shift of the higher mass velocity data to lower qualities for the same heat flux. These observations suggest a possible modification of the Biasi correlation. For high qualities the correlation is proportional to (1- X): q"CHF cc(1-X) (6.5) To account for entrainment either the quality could be effectively increased, due to the lower liquid flow rate on the heater surface, or the functional form of the quality dependence could be modified. The latter option is probably the better approach. For example, one possibility is to write q"CHF 01(A(G)-X) (6.6) Where A(G) is a parameter which is less than 1.0 and depends on the mass velocity. For the Bennett data, A(G) has been found to be A(G) = 4.8G-.2 5 (6.7) where G is in kg/m2as -215- for 2 2 650 kg/m.s < G < 1500 kg/m s Using this modification the previously discussed cases have been reanalyzed and the results are presented in Table 6.3. It is seen that the predictions of the modified Biasi correlation are in much better agreement with the data over the entire range of mass velocities. It should be noted that this modification only affects the high quality CF regime. Of course, the functional form of A(G) is not expected to be valid outside the above range. Nevertheless, this modification is justified on a physical basis and with additional data a more appropriate form could be obtained. The results of the post-CHF regime assessment have already been discussed in Section 4.3.3. As shown in that section, if the vapor temperature can be accuratelycalculated,_then-the--wall--temperatures-can be satisfactorily predicted if a single-phase vapor heat transfer correlation is used. Typical comparisons have been shown in Figures 4.9 to 4.11. The agreement between the measured and predicted temperatures is quite good over the entire range of mass Velocities. The only minor negative trend is that the code slightly underpredicts many of the wall temperature measurements. This result may indicate the need to evaluate the assumption concerning the use of a single-phase vapor correlation. However, an equilibrium film boiling correlation will generally, for the same conditions, yield a higher heat transfer coefficient and, thus, lead to even lower wall temperatures. It seems unlikely, then, that a -216- TABLE 6.3 CHF Comparisons for Bennett Test Cases with Corrected Biasi Correlation Case No. G2 2 Z Measured zBiasi (corrected) kg/m -s MW/m (m) m) _ _ _ _ I _ _ I _ _ _ __ I 5325 5332 5336 5266 5276 5273 5245 5247 5253 5391 5394 5396 665 665 665 1004 1017 1017 1356 1356 1356 5181 5181 5181 0.572 0.649 0.822 0.704 0.795 0.922 0.804 0.845 0.903 1.66 1.76 1.84 5.16 4.39 3.56 5.00 4.39 3.56 4.95 4.70 4.39 5.00 4.55 4.17 5.16 4.44 3.37 5.07 4.39 3.61 4.97 4.78 4.29 4.43 3.91 3.57 -217- film boiling correlation should be used. Thus, a possible remedy is to use a slightly modified single-phase vapor correlation. Nevertheless, the degree of agreement between the measurements and the present code predictions is quite satisfactory. 6.3.2 Transient Results Transient CHF measurements taken in a 9-rod bundle have been compared with THERMIT predictions in order to assess the transient CUF predictive capability of THERMIT. In the experimental tests, the rod bundle was initially brought to a steady flow condition with exit qualities on the order of 25%. The transient was then initiated by instantaneously halving the inlet flow rate using a fast closing value. The flow rate was maintained at this lower value until CHF occurred. Both the time-to-CHF and the approximate CHF location have been measured. By comparing these measurements with the predictions of THERMIT, the CHF predictive capability could be assessed for this severe flow transient. Representative cases from those given in Reference 38, have been simulated with THERMIT. In each case a steady-state solution is first obtained after which the transient is initiated. To generate the steady-state solution the initial conditions given in Table 6.4 have been used. The geometrical features of the test section are illustrated in Fig. 6.11 and it should be noted that the axial and radial power distributions are uniform in each case. The cases listed in Table 6.4 have been analyzed with THERMIT using both the Biasi and CISE-4 correlations. Hence both correlations could be evaluated for these cases. (Once again the W-3 correlation is not applicable for these cases). -218- TABLE 6.4 Summarytof Conditions for Transient CHF Cases Case No. 9-151 9-170 9-175 9-179 G(kg/m2s) 1366 854 844 1383 Power (kW) 1093 834 812 841 Aisub (kJ/kg) 70 96 107 97 Xout (%) 23.0 27.2 26.1 14.4 0 Low!GHigh 0.49 0.46 0.47 0.3 Side Subchannel /2I 2 3 1- 2' 2i 1 1 2 ii1 Center Subchannel Corner Subchannel Geometrical Details Rod Diameter 14.3 mm Rod-Rod Gap 4.42 mm Rod-Wall Gap 3.51 mm Radius of Corner 10..2 mm Heated Length 1829 mm Figure 6.11 Cross Sectional View of G.E. 9 Rod Bundle Used in Transient CUF Tests -219- 3 3 2 2 '%W-AMA I -220- The results of these comparisons show that both CF correlations predict premature CHF in all cases. This can be clearly seen in Figures 6.12 to 6.15 where the minimum critical heat flux ratio (MCHFR) versus time is illustrated. In each case, the time-to-CHF is substantially underpredicted by the CISE-4 correlation. The Biasi correlation predicts the time-to-CHF fairly well in all cases except case 9-179. The location of the CHF point is consistently measured to occur on the corner rod facing the corner channel. This result seems surprising in view of the meadsured mass velocity and quality profiles for the G.E. 9-rod bundle discussed in Section 5.2. However, the test section used here had grid spacers which tend to distort the "clean" bundle quality and mass velocity distributions. In fact, a completely opposite trend in the measured mass velocities is observed [52]. In the "clean" 9-rod bundle the mass velocity is highest in the center, with the lowest mass velocity occurring in the corner subchannel. In the bundle with grids the mass velocity is lowest in the center and highest in the corner. This redistribution is attributed to the fact that the grid loss coefficient is highest in the center subchannel and lowest in the corner subchannel. Since the CHF is inversely proportional to the mass velocity, the corner subchannel with the highest mass velocity would be expected to reach CHF first. However, in THERMIT proper modeling of this grid effect is currently not included. Consequently, even though CHF is first predicted to occur on the corner rod, it does so on the side facing the center subchannel. In order to verify that this prediction is indeed due to the grid modeling, case 9-151 has been reanalyzed with non-uniform grid -221- - Biasi 3.0 - - CISE-4 2H Ei 1 .0 0. 1.0 2.0 TIME (seconds) Figure 6.12: Comparison of MCHFR Predictions versus Time for Case 9-151 -222- - Biasi - - CISE-4 3.0 E i2.0 tCHF 1.0 N O'. 0.1 0.2 TIME (seconds) Figure 6.13 : Comparison of MCHFR Predictions versus Time for Case 9-170 MMMNMMM-dMWl -223- -- Biasi 3.0 -- CISE-4 0 2.0 - H 1.0 0. 1.0 2.0 TIME (seconds) Figure 6.14 : Comparison of MCHFR Predictions versus Time for Case 9-175 -224- -- Biasi 3.0 - - CISE-4 0 H E < 2.0 z Ht HH 1.0 t'CH 0. 0.1 0.2 TIME (seconds) Figure 6.15: Comparison of MCHFR Predictions versus Time for Case 9-179 -225- coefficients (i.e., lowest in the corner and highest in the center). As seen in Fig. 6.16, the time-to-CHF prediction is essentially the same, but the first occurrence of CHF is now in the appropriate location. This result clearly indicates the large sensitivity to the proper modeling of the grid spacers. However, the grid characteristics are not always known in advance which limits the modeling capability. The problem with underpredicting the measured time-to-CHF has also been observed by Van Haltern [63]. In his work, which involved analyzing flow decay transients in a 16-rod bundle, a number of well-known CHF correlations were compared to the data. In every case, substantial underprediction of the time-to-CHF has been found for all correlations except the one formulated with the actual data (i.e., a test section specific correlation). Leung [64] has also analyzed a number of flow transients using various CHF correlations. For the range of experiments which were analyzed, the Biasi correlation was found to satisfactorily predict the data. However, all of these studies were for tubes (i.e., one- dimensional), so that the three-dimensional effects found in the G.E. tests were not present. Thus, these results indicate that the Biasi correlation is appropriate for transient CHF predictions provided the correct flow distribution is calculated. Consequently, Leung has concluded that the use of steady-state CHF correlations for transient application is appropriate. In order to study the sensitivity of the predicted results to the two-phase flow modeling case 9-179 has been reanalyzed with a smaller value of the convergence criterion. The result of this variation is Base Case - -- Modified Grid Coefficients I' CHF I 1.0 TIME (seconds) Figure 6.16: Comparison of Biasi Correlation Predictions for Case 9-151 with Modified Grid Coefficients -226- 3.0 1- 0) H x- 8w z 2.0 1.0 0. 2.0 -227- illustrated in Fig. 6.17. Using the Biasi correlation, it is seen that with a tighter criterion, slightly different MCHFR's are predicted, but the time-to-CHF is virtually unchanged. However, an interesting result is obtained when the tighter criterion is used. With the tighter criterion, THERMIT predicts a rewetting of the rods. That is, after the initial CHF occurrence, the MCHFR increases to a value larger than 1.0 and then oscillates until 2.0 seconds. After 2.0 seconds the rewetting ends and the clad temperature increases rapidly. The reason for this resetting, which may reflect the actual situation, is THERMIT's ability to calculate the local pressure field. After CHF occurs, the wall friction, which previously had been solely apportioned to the liquid, is now apportioned to vapor. This introduces a local perturbation in the pressure distribution which tends to decrease the void fraction. The decrease in void fraction decreases the quality which leads to a higher critical heat flux. Eventually, rhdn, the MCHFR becomes greater than 1.0. However, once this occurs the pressure distribution changes again and CHF is again predicted to occur. This oscillatory behavior continues for a short time until the quality becomes so high that the perturbations in the pressure field can no longer suppress the CHF condition. Apparently a tightly converged pressure distribution is required to observe this phenomenon. The net effect of this rewetting process is to maintain the wall temperature near the steady-state value. Since CHF is experimentally observed by noting the rapid increase in the wall temperature, the code predicted CHF based on this criterion is in good agreement with the data. Figure 6.18 illustrates the wall temperature profiles for case 9-179. -228- -- Base Case - Tighter Convergence Criterion 1.0 \\ /f t CHF 2.0 TIME (seconds) Figure 6.17: Comparison of Biasi Correlation Predictions with a Tighter Convergence Criterion (Case 9-179) 3. 0 H H E- 2. I3 1.0- 0. If'( % / p -229- 600 1 590 [ 630 620 610 0. 1.0 THERMIT - THERMIT Tighter Convergende Data -e -- - -- - .~~1 / / / / / / / / 2.0 TIME (seconds) Figure 6.12: Comparison of Measured and Predicted Maximum Wall Temperatures for Case 9-179 / / 9,0M- 580 570 I I -230- It is seen that, when the tight criteria is used, the wall temperaure does not increase rapidly at the first occurrence of CHF. The predicted time-to-CHF based on the wall temperature is also seen to be in good agreement with the data. This observation can also be seen for the case 9-151 results illustrated in Fig. 6.19. Again, based on the temperature profiles the predicted and measured CHF values are found to be in good agreement. In summary, the key features of the transient CHF capability of THERMIT have been identified. The CISE-4 correlation consistently under-predicts the time-to-CHF, while the Biasi correlation is found to be in good agreement with the data provided accurate flow conditions are calculated. The flow conditions are found to be very sensitive to the use of grid spacers. Clearly, more work is required to develop an appropriate grid model for THERMIT. Finally, it is seen that if the tight convergence is used, then a rewetting phenomena is predicted. This rewetting delays the apparent time-to-CHF which results in better comparisons between the measured and predicted data. -231- THERMIT THERMIT - Tighter Convergence Data -A 000 .. in 2.0 TIME (seconds) Figure 6.19: Comparison of Measured and Predicted Maximum Wall Temperatures for Case 9-151 6301 620 1-- 610 L. rzi E- 600 [ 580 5701 1 0 1.0 a I 590 1 k -232- 7.0 Conclusions and Recommendations 7.1 Conclusions This research has addressed the problem of developing a two-fluid model suitable for two-phase flow analysis of LWR rod-bundles on a sub- channel basis. The approach used in this work has been to modify arid improve the two-fluid model computer code THERMIT and then assess the predictive capabilitites of the code for subchannel applications. The final result of this research is that the first two-fluid model code capable of subchannel analysis has now been developed. The work relating to the modification of THERMIT has been directed towards improving the code's predictive capabilities with the primary emphasis being placed on LWR subchannel analysis. To meet this objective several modifications and improvements of THERMIT have been made with the major ones being listed in Table 7.1. The two improvements specifically required for subchannel analysts are the expansion of the geometrical modeling capability and the addition of a two-phase mixing model. The geometrical improvements were necessary to allow both coolant-centered subchannel analysis and detailed fuel rod modeling. The addition of the two-phase mixing model was necessary to account for the inter-channel exchange processes arising from turbulent mixing and vapor diffusion. The other modifications have improved the overall predictive capabilities of THERMIT. The addition of droplet vaporization model and CHF correla- tions has eliminated modeling deficiencies which existed in the original version of THERMIT. The modification of the interfacial energy exchange model and the wall heat transfer logic has replaced previously simplified -233- TABLE 7.1 Summary of Modifications and Improvements Made in THERMIT 1. Expanded Geometrical Modeling to Allow - Coolant-Centered Subchannel Analysis - Detailed Fuel Rod Analysis 2. Added Two-Phase Mixing Model 3. Added Droplet Vaporization Model 4. Improved Liquid-Vapor Interfacial Energy Exchange Model 5. Modified Wall Heat Transfer Model by - Adding Additional CHF Predictive Capability - Modifying Post-CHF Heat Transfer Selection Logic S h -234- modeling with more realistic modeling. Hence, the modifications made in THERMIT have expanded the analytical capability to allow subchannel analysis as well as substantially improving the two-phase flow modeling. The assessment effort, performed in conjunction with the developmental work, has evaluated the important models in THERMIT. Since the analysis of LWR rod-bundles on a subchannel basis is the primary application of this research, the models have been assessed for conditions which are consistent with those found in LWR rod-bundles. The method used in this assessment has been to compare experimental measurements with the code predictions. The experiments have been chosen to resemble either BWR or PWR flow conditions with the geometry being either rod-bundle or tube. (A subchannel is geometrically similar to a tube.) The order in which these tests have been performed has been chosen so that separate assess- ment of the individual models could be made. Although certain assumptions were necessary to perform the evaluation, this assessment strategy is an appropriate method for evaluating the models. The models which have been assessed include: the liquid-vapor inter- facial exchange models, the two-phase mixing model and the wall heat transfer model. From this assessment effort the following conclusions can be made about these models. The interfacial exchange models control the transfer of mass, energy and momentum between the liquid and vapor and determine the degree of thermal and mechanical non-equilibrium in the two-phase flow. The assessment of these models has been concerned with verifying that the proper non-equilibrium conditions are predicted by THERMIT. Each model has been assessed separately and the overall conclusion of this -235- evaluation work is that these models are appropriate for LWR rod-bundle analysis. The interfacial mass exchange model has been shown to predict the proper rate of mass exchange between the liquid and vapor for both pre-CHF and post-CHF conditions. (Depressurization transients have been excluded from this assessment effort.) In the pre-CHF regime, comparisons with void fraction measurements have illustrated the appropriateness of this model for subcooled as well as saturated boiling conditions. For post-CHF conditions, THERMIT accurately predicts the amount of vapor superheat and the wall temperatures which implies that the droplet vapori- zation rate is properly modeled. These results show that the interfacial mass exchange model in THERMIT can accurately analyze the various types of vaporization mechanisms anticipated for LWR rod-bundle steady-state and non-depressurization transients. Using the current formulation of the interfacial energy exchange model, the proper liquid and vapor temperature distributions are predicted for thermal non-equilibrium conditions in two-phase flow. For subcooled boiling, the bulk liquid temperature is predicted to be subcooled while the liquid is saturated. For post-CHF conditions, the vapor is predicted to be superheated while the liquid is saturated. These results illustrate that the code can predict the appropriate temperature distributions for thermal non-equilibrium conditions. Also, for transient applications, the current formulation of the interfacial energy exchange model is necessary to prevent numerical instabilties. Hence, the interfacial energy exchange rate is now properly modeled for both steady-state and transient conditions. -236- Even though the interfacial momentum exchange model in THERMIT treats the interfacial forces in a simplified manner, the model is still able to satisfactorily predict the momentum exchange between the liquid and vapor. Comparisons of the code predictions with void fraction measurements have indicated that the appropriate relative velocity is predicted. This result can be used to infer that the interfacial momentum exchange rate is proper. It should be noted that this conclusion is based on steady-state saturated boiling measurements. For droplet flows, very low flows, or rapidly accelerating flows, the model has not been assessed and should be used with some caution for these flow condi- tions. However, for steady-state and near-operational transient rod bundle analysis, the interfacial momentum exchange model should be appropriate. The important consequence of the above results is that both thermal and mechanical non-equilibrium between the phases are accurately predicted by THERMIT. The non-equilibrium phenomena are accurately represented by the interfacial exchange models. Furthermore, these models have been formulated in a mechanistic manner and are expected to be applicable over a wide range of flow conditions. Hence, non-equilibrium between the phases is accurately predicted in a realistic manner which is an important advantage of the two-fluid formulation used in THERMIT. In addition to evaluating the interfacial exchange terms, this research has also addressed the problem of developing an appropriate model to account for the effects of turbulence. A two-phase mixing model has been incorporated into THERMIT for this purpose and has been assessed using a number of rod-bundle experimental measurements. This work has -237- shown that with the two-phase mixing model, THERMIT can accurately predict the flow and enthalpy distributions in both BWR and PWR rod- bundles. One important result is that TbERMIT can correctly predict the measured trend in the corner subchrniael quality. This trend is that the quality is much lower than the bundle average for BWR conditions, while being near the average for PWR conditions. THERMIT accurately predicts this behavior using the same mixing model parameters in each case. Hence, the two-phase mixing model is valid over the range of pressures which are typical of BWR and PWR rod-bundle conditions. Another important result of this assessment is that the effects of grid spacers must be carefully modeled in order to predict the correct flow distribution. The grid spacers can significantly alter the flow distribution in a non-uniform manner. Proper modeling of the grid is needed to predict the appropriate trends in the mass velocity measurements. However, the quality distribution is rather insensitive to the grid spacer modeling. The final model to be assessed is the wall heat transfer model. Three aspects of the model have been investigated: pre-CHF correlations, post-CHF correlations and steady-state and transient CHF predictive capability. Overall, the model is able to satisfactorily predict the experimental data to which it has been compared. However, certain areas of the model may need to be improved. For pre-CHF conditions the Cheti correlation is found to underpredict the heat transfer coefficient except at low heat fluxes. Although this result leads to conservative wall temperature predictions which are -238- probably satisfactory for many applications, the use of an alternative heat transfer correlation may be needed. Either the Thom or Jens-Lottes correlation can appropriately predict the heat transfer coefficient for all cases which have been studied. However, neither of these two correla- tions can calculate the heat transfer coefficient for forced convection vaporization while the Chen correlation is able to calculate this mode of heat transfer. This type of heat transfer is anticipated to be important for BWR conditions. It should also be noted that for typical BWR heat fluxes, the Chen correlation should also calculate the appropriate heat transfer coefficient. Hence, for BWR rod-bundle analysis, the Chen correlation should predict satisfactory results. For post-CIF conditions, accurate wall temperature predictions are more dependent on the vapor temperature calculation than on the heat transfer correlation. Consequently, if the droplet vaporization model is used, accurate wall temperatures can be predicted using a single-phase vapor heat transfer coefficient. This type of modeling is currently included in THERMIT and it has been shown that the appropriate post-CHF temperature distributions are predicted. In evaluating the CHF predictive capability of THERMIT, both the CISE-4 dorrelation and Biasi correlation have been assessed for steady- state arid transient conditions. The CISE-4 correlation has been found to underpredict the critical heat flux (or more appropriately the critical power) in all cases. The Biasi correlation while underpredicting the critical heat flux for transient conditions, usually overestimates the critical heat flux for the steady-state tests which were studied. The agreement is poorest for high qualities and low heat flux cases and is -239- probably due to the failure of the correlation to properly account for en- trainment. A simple modification of this correlation for high qualities has significantly improved the predictions. On the whole, both correla- tions will probably be satisfactory for BWR steady-state or transient cases. For the transient CF analysis, THERMIT is able to predict a rewetting of the wall after an initial occurrence of CHF. This prediction is a direct result of THERMIT's ability to calculate the local pressure distribution. At the initial CHF location the pressure distribution is perturbed causing the flow to accelerate locally, thus lowering the quality and rewetting the wall. Consequently, the predicted time of CHF based on the wall temperature excursion is increased such that the predic- tions are in good agreement with the data. In summary, it can be concluded that THERMIT can now successfully analyze LWR rod bundles on a subchannel basis. The geometrical and physical modeling capability needed for this type of analysis has been added to the code. Assessment of the important models for conditions typical of LWR rod-bundles has shown that appropriate results are predicted by the code. Hence, the main objective of this research has been accomplished, since THERMIT is the first two-fluid model code which has been developed and tested for LWR subchannel applications. Specific details on the programming and usage of THERMIT may be found in Ref. 9. 7.2 Recommendations In the course of this research, a number of possibilities for future study have become apparent. First, even though the interfacial momentum exchange model has been found to be appropriate for pre-CHF conditions, further evaluation of this model for post-CHF applications is probably -240- warranted. Also, if the range of application of the code is to be expanded then additional interfacial forces may need to be included. In fact, it may be necessary to introduce a flow regime dependent interfacial momentum exchange model although for the present applications this is apparently not needed. Second, the wall heat transfer model should be examined in greater de- tail and perhaps modified if necessary. In particular, the nucleate boiling regime may need to be expanded by adding additional, and more accurate, correlations especially for high heat flux cases. The CHF predictive capability also needs to be evaluated for PWR conditions where the W-3 correlation has been implicitly assumed to be accurate. This assumption needs to be tested for both steady-state and transient conditions. However, it must be remembered that heat transfer modeling is limited by the restrictions placed on individual correlations. Thus, since THERMIT can analyze a wide range of two-phase flow conditions, it is necessary that compatible heat transfer correlations be included in the code. Third, the importance of accurate grid spacer modeling for rod- bundle analysis dictates that improved modeling of the grid effects be included in THERMIT. The grids can significantly alter the flow distribution which becomes important for critical heat flux predictions in rod-bundles. Unfortunately, much of the grid spacer information is proprietary which makes it difficult to develop a general model. Fourth, it has been implicitly assumed throughout this work that the wall friction model is appropriate. This assumption, while probably valid for steady-state conditions, needs to be investigated for transient -241- applications. Fifth, although assessed for steady-state conditicns, most constitutive models have been assumed to be valid for transient applica- tions. This assumption has been verified for some conditions [64] but cannot be assessed for all transient conditions due to the lack of experimental measurements. However, as more transient data becomes available it is important that further assessment of the transient capabilities be performed. Sixth, with the ability to accurately model the flow distribution in a rod-bundle, it is possible to compare CHF predictions using one-dimen- sional modeling to those obtained using three-dimensional modeling. This evaluation would indicate the appropriateness of bundle average CHF calculations. Finally, the computational efficiency of THERMIT should be improved. This improvement is especially needed for steady-state analysis of rod- bundles when small mesh sizes are used. There are probably a number of techniques for improving the efficiency and these should be investigated. -242- REFERENCES 1. C. L. Wheeler et al., "COBRA-IV-I: An Interim Version of COBRA for Thermal-Hydraulic Analysis of Rod Bundle Nuclear Fuel Elements and Cores," BNWL-1962, (1976). 2. D. S. Rowe, "COBRA IIIC: A Digital Computer Program for Steady-State and Transient Thermal Hydraulic Analysis of Rod Bundle Nuclear Fuel Elements," BNWL-1695, (1973). 3. W. H. Reed and H. B. Stewart, "THERMIT; A Computer Program for Three-Dimensional Thermal-Hydraulic Analysis of Light Water Reactor Cores," M.I.T. Report prepared for EPRI, (1978). 4. J. E. Kelly and M. S. Kazimi, "Development and Testing of the Three- Dimensional, Two-Fluid Code THERMIT for LWR Core and Subchannel Applications, "M.I.T. Energy Laboratory, Report No. MIT-EL 79-046, (1979). 5. D. R. Liles and W. H. Reed, "A Semi-Implicit Method for Two-Phase Fluid Dynamics," Journal of Computational Physics, 26, p. 390, (1978). 6. F. H. Harlow and A. A. Amsdem, "A Numerical Fluid Dynamics Calcula- tion Method for All Flow Speeds, " Journal of Computational Physics, 8, p. 197, (1971). 7. V. L. Shah, et al., "Some Numerical Results with the COMMIX-2 Computer Code," NUREG-CR-0741, (1979). 8. M. J. Thurgood and J. M. Kelly, "COBRA TF: Model Description," NRC Advanced Code Review Group Meeting, Silver Spring, Maryland, (1980). 9. J. E. Kelly and M. S. Kazimi, "THERMIT-S: A Two-Fluid Model Computer Code for LWR Rod Bundle Analysis," M.I.T. Energy Laboratory Report, (to be published). 10. L. Wolf, et al., "WOSUB, A Subchannel Code for Steady-State and Transient Thermal-Hydraulic Analysis of BWR Fuel Pin Bundles," MIT-EL-78-025, (1977). 11. W. T. Sha., et al., "COMMIX-1, A Three-Dimensional Transient Single Phase component Computer Program 'for Thermal-Hydraulic Analysis," NUREG-0415, (1978). 12. J. F. Jackson, et al., "TRAC-PIA: An Advanced Best-Estimate Computer .Computer Program for PWR LOCA Analysis," NUREG/CR-0665, (1979). 13. M. Massoud, "A Condensed Review of Nuclear Reactor Thermal-Hydraulic Computer Codes for Two-Phase Flow Analysis," M.I.T. Energy Laboratory, Report No. MIT-EL 79-018, (1980). -243- 14. W. T. Sha, "A Summary of Methods Used in Rod-Bundle Thermal-Hydraulic Computer Codes for Two-Phase Flow Analysis," Turbulent Forced Convec- tion in Channels and Rod Bundles Conference, Instanbul, (1978). 15. J. A. Boure, "Mathematical Modeling and the Two-Phase Constitutive Equations," European Two-Phase Flow Group Meeting, Haifa, Israel, (1975). 16. J. Weisman and R. W. Bowring, "Methods of Detailed Thermal and Hydraulic Analysis of Water-Cooled Reactors," Nuclear Science and Engineering, 57, p. 255, (1975). 17. J. J. Ginoux, Two-Phase Flows and Heat Transfer with Application to Nuclear Reactor Design Problems, McGraw-Hill, New York, (1978). 18. G. P. Gaspari, A. Hassid and G. Vanoli, "Some Considerations on Critical Heat Flux in Rod Clusters in Annular Dispersed Vertical Upward Two-Phase Flow," Energia Nucleare, _7, p. 643, (1970). 19. C. W. Stewart, et al, "COBRA IV: The Model and the Method," BNWL-2214, (1977). 20. T. A. Bjornard, "Blowdown Heat Transfer in a Pressurized Water Reactor," Ph.D. Thesis, Department of Nuclear Engineering, M.I.T., (1977). 21. H. B. Stewart, "Stability of Two-Phase Flow Calculations Using Two-Fluid Models," Journal of Computational Physics, 33, No. 2, p. 259, (1979). 22. J. G. Delhaye and J. L. Achard, "On the Averaging Operators Introduced in Two-Phase Flow Modeling," Specialists' Meeting on Transient Two-Phase Flow, Toronto, (1976). 23. J. F. McFadden, et al, "An Evaluation of State-of-the-Art Two- Velocity Two-Phase Flow Models -and Their Applicability to Nuclear Reactor Transient Analysis," EPRI NP-143, (1976). 24. J. T. Rogers, and N. E. Todreas, "Coolant Interchannel Mixing in Reactor Fuel Rod Bundles Single-Phase Coolants," in Heat Transfer in Rod Bundles, ASME, (1968). 25. R. T. Lahey and F. J. Moody, The Thermal Hydraulics of a BWR, ANS Monograph, (1975). 26. J. M. Gonzalez-Santalo and P. Griffith, "Two-Phase Flow Mixing in Rod Bundle Subchannels," ASME, 72-WA/NE-19, (1972). 27. D. S. Rowe and C. W. Angle, "Crossflow Mixing Between Parallel Flow Channels During Boiling," BNWL-371, (1967). 28. S. G. Beus, "Two-Phase Turbulent Mixing Model for Flow in Rod Bundles," WAPD-T-2438, (1970). -244- 29. A. J. Faya, "Development of a Method for BWR Subehannel Analysis," Ph.D. Thesis, M.I.T. Department of Nuclear Engineering, (1979). 30. J. T. Rogers and R. G. Rosehart, "Mixing by Turbulent Interchange in Fuel Bundles, Correlations and Inferences," ASME, 72-HT-53, (1972). 31. G. B. Wallis, One Dimensional Two-Phase Flow, McGraw-Hill Book Company, New York, (1969). 32. Drew and R. T. Lahey Jr., "An Analytic Derivation of a Subchannel Void-Drift Model," Transactions of the ANS, 33, (1979). 33. G. W. Maurer, "A Method of Predicting Steady-State Boiling Vapor Fraction in Reactor Coolant Channels," WAPD-BT-19, (1960). 34. H. Christensen, "Power-to-Void Transfer Functions," ANL-6385, (1961). 35. J. F. Marchaterre, et al, "Natural and Forced-Circulation Boiling Studies," ANL-5735, (1960). 36. R. T. Lahey Jr., et al, "Two-Phase Flow and Heat Transfer in Multirod Geometries: Subchannel and Pressure Drop Measurements in a Nine-Rod Bundle for Diabatic and Adiabatic Conditions," GEAP-13049, (1970). 6 37. H. Herkenrath and W. Huf schmidt, "Experimental Investigation of the Enthalpy and Mass Flow Distribution Between Subchannels in a BWR Cluster Geometry (PELCO-S)," EUR-6585-EN, (1979). 38. H. Herkenrath, W. Hufschmidt and L. Wolf, "Experimental Investigation of the Enthalpy and Mass Flow Distributions in Subchannels of a 16-Rod PWR Bundle (EUROP)," European Two-Phase Flow Group Meeting, Glasgow, (1980). 39. A. W. Bennett, et al, "Heat Transfer to Steam-Water Mixtures Flowing in Uniformly Heated Tubes in Which the Critical Heat Flux Has Been Exceeded," AERE-R-5373, (1967). 40. B. S. Shiralkar, et al, "Transient Critical Heat Flux - Experimental Results," GEAP-13295, (1972). 41. S. Nijhawan, et al, "Measurement of Vapor Superheat in Post-Critical Heat Flux Boiling," 18th National Heat Transfer Conference, San Diego, (1979). 42. P. Saha, "A Non-Equilibrium Heat Transfer Model for Dispersed Droplet Post-Dryout Regime," International Journal of Heat and Mass Transfer, 23, p. 481, (1980). -245- 43. S. Y. Ahmad, "Axial Distribution of Bulk Temperature and Void Fraction in a Heated Channel with Inlet Subcooling," Journal of Heat Transfer, 92, p. 595, (1970). 44. P. Saha and N. Zuber, "Point of Net Vapor Generation and Vapor Void Fraction in Subcooled Boiling," Proceeding of the 5th International Heat Transfer Conference, (1974). 45. S. L. Soo, Fluid Dynamics of Multiphase Systems, Blaisdell Publishing Co., Waltham, MA, (1967). 46. V. G. Levich, Physicochemical Hydrodynamics, Prentice-Hall, Inc., Englewood Cliffs, NJ, (1962). 47. L. Y. Cheng, D. A. Drew, and R. T. Lahey, Jr., "Virtual Mass Effects in Two-Phase Flow," NUREG/CR-0020, (1978): 48. M. Ishii and N. Zuber, "Drag Coefficient and Relative Velocity in Bubbly, Droplet or Particulate Flows," AIChE Journal, 25, No. 5, p. 843, (1979). 49. A. B. Basset, Hydrodynamics, Dover Publications Inc., New York, (1961). 50. W. C. Rivard and M. D. Torrey, "Numerical Calculation of Flashing from Long Pipes Using a Two-Fluid MoeIel," LA-6104-MS, (1975). 51. R. T. Lahey, Jr., "Out-of-Pile Subchannel Measurements in a Nine-Rod Bundle for Water at 1000 PSIA," Progress in Heat and Mass Transfer, 6, p. 345, (1972). 52. R. T. Lahey, Jr., et al, "Deficient Cooling, 7th Quarterly Progress Report," GEAP-10221-7, (1971). 53. H. Herkenrath, W. Hufschmidt, and F. Lucchini, "Experimental Subchannel Investigation in a 16 Rod Test Section by Means of the Isokinetic Sampling Technique," 2nd Multi-Phase Flow and Heat Transfer Symposium-Workshop, Miami Beach, (1979). 54. F. S. Castellana and J. E. Casterline, "Subchannel Flow and Enthalpy Distributions at the Exit of a Typical Nuclear Fuel Core Geometry," Nuclear Science and Engineering, 22, p. 3, (1972). 55. L. Biasi, et al, "Studies on Burnout, Part 3," Energia Nucleare, 14, No. 9, p. 530, (1967). 56. L. S. Tong, Boiling Crisis and Critical Heat Flux, TID-25887, AEC Critical Review Series, (1972). 57. G. P. Gaspari, A. Hassid and F. Lucchini, "A Rod-Centered Subchannel Analysis with Turbulent (Enthalpy) Mixing For Critical Heat Flux Prediction in Rod Clusters Cooled by Boiling Water," 5th International Heat Transfer Conference, Tokyo, (1979). -246 58. T. A. Bjornard and P. Griffith, "PWR Blowdown Heat Transfer," Thermal and Hydraulic Aspects of Nuclear Reactor Safety - Vol. 1, Light Water Reactors ASME Symposium, (1977). 59. J. C. Chen, "A Correlation for Boiling Heat Transfer to Saturated Fluids in Convective Flow," ASME, 63-HT-34, (1963). 60. J. R. S. Thom, et al, "Boiling in Sub-cooled Water During Flow Up Heated Tubes and Annuli," Proceedings of the Institution of Mechanical Engineers, 180, part 3C, p. 226, (1965). 61. W. H. Jens and P. A. Lottes, "Analysis of Heat Transfer, Burnout, Pressure Drop and Density Data for High Pressure Water," ANL-4627, (1951). 62. J. G. Collier, Convective Boiling and Condensation, McGraw-Hill Book Co., New York, (1972). 63. M. L. Van Haltern, "A Sensitivity Study of Thermal-Hydraulic Correla- tions in the Computer Code MEKIN," S.M. Thesis, Department of Nuclear Engineering, M.I.T., (1980). 64. J. C. Leung, "Transient Critical Heat Flux and Blowdown Heat Transfer Studies, " Ph.D. Thesis, Department of Chemical Engineering, Northwestern University, (1980). -247- Appendix A: Derivation of the THERMIT Conservation Equations A.1 Introduction The procedure for deriving the governing equations in THERMIT from the local, instantaneous balance equations is presented in this section. This derivation consists of a straightforward application of time and space averaging operators along with clear identification of the resulting terms. Examples of the derivation approach for a two- fluid model may be found in references 22 and 23. Besides identifying the assumptions used in the THERMIT conserva- tion equations, two important terms have now been retained in the equations which were not included in the original THERMIT equations. The first is the term describing the effects of turbulent mixing. This term originates from the averaging procedure and must be included for subchannel or plena applications. The second term is that describing the liquid-to-vapor transport of either energy or momentum due to mass exchange across liquid-vapor interfaces. This term arises from the application of the volume averaging operator and is important for transient applications. The inclusion of these terms enhances the physical modeling capability of THERMIT. Before presenting the actual derivation, the necessary mathematical theories and notation will be discussed. Since both time and volume averaging are required in this derivation, it is convenient to intro- duce notation to represent these terms. The time average of an arbi- trary variable B is defined by B T Bk dt (A.l)k T k -248- where T is a characteristic time interval for time averaging. In a similar way, the volume average of the variable B is defined by - (A.15) the first term can be written as { (Pk k)dW ak y JPkkTVi-n dA (A.16) at at W JAkk The Gauss theorem Eq. (A.6) is now applied to the second term of Eq. (A.14 ) to obtain WjV'(PkTk k + k Jk VLxk Pkrk Vk Jk>) (A.17) 1 - + 4+ -+ W (pk k k k)-n dAA The area over which the surface integrals are applied includes all internal surfaces within volume W. Hence, these integrals can be divided into two components; one for the surfaces between phases and one for the fluid-solid surfaces. The general surface integral can be written as -254- B dA B dA + B d A (A.18) A -A A S where A is the area between phases and As is the fluid-solid surface area. This property of the surface integrals will be used to simplify the conservation of energy and momentum equations. Combining Eqs. (A.14, A.16 & A.17) and using Eq. (A.15) one obtains E(k k) + k kT kk + = 0). Substituting these expressions into Eq. (A.17) and simplifying yields the following expression -255- 9t(cakY) + 7.(ctky KVk> +) - akYP k >k - < A(P I< V Vk) J >-.ndAW1 +k - k +k> Wk + ->())- An dA = Ic( k k k k V>+k Vk k k - k - fA (A.21) k k k k k k - ck k k - k This expression may be simplified as follows. Since THERMIT impli- citly accounts for volume porosity in the formulation of its finite difference equations, the porosity factor , Y, can be set equal to unity. Further simplifications can be made by assuming that the fluctuating density terms are zero: Pk' = 0 (A.22) Physically this assumption means that the phase density is uniform in the control volume and is valid provided the volume is not too large. By applying this assumption and lumping together all non-integral fluctuating terms into M'kEq. (A.21) may be written as -256- at (k >kk> + Ve(cCk( + qk> 1 J -kkik + +; +(A.23) =- (

(V - V )-)-n)dA + -|

T '(V - V)W k k i k k W k k i k - Jk')-ndA - M'k This equation is the general form of the time and volume averaged conservation equation. The mass, momentum and energy equations will now be obtained by substituting the appropriate values for Tk' k and $k into Eq. (A.23). A. 4 Conservation of Mass The conservation of mass equation is obtained by substituting T= 1, $k 0 k= 0 into Eq. (A.23). Performing this substitution, one obtains at k k) + S7(tkPk> 4k) A J (Vi )-ndA - M'k (A.24) where M'k represents the mass exchange due to turbulent fluctuations. The surface integral can be identified as the interfacial mass exchange rate, r'k. Therefore, Eq. (A.24) can be written as (ak k>)+ & kk> k- k (A.25) This equation is the same as that used in THERMIT with the time and space average notation dropped and using -257- a =ct v 1' =1 ' v ' = -1 (A.26) M' = W v tv M'=W St% With these expressions the vapor and liquid mass conservation equations are given by: Vapor Mass Equation (ap ) + V-(ctp V ) = - W (A.27)3t v v V tv Liquid Mass Equation ((1- a)p ) + V-((l- a) pt V) = - - Wt (A.28) A.5 Conservation of Momentum The conservation of momentum equation, which is actually a vector =4. ' - nEq A2) equation, is found by using k=k' Ok - g and 1k=k k in Eq. (A.23). Since viscous stress terms are assumed to be small and are not included in the present THERMIT formulation the viscous stress term Tk is set to zero. -258- --t- (a k)) + V- (CA k) + VC k< Pk > - a k gk k vk) ckkkkkak YE -N++ V~c +p +> - + - [(

(V V )on

ndA(A.29) WA k k k k W JAkk k(V V n Pk'n] dA - M'k The area integrals can be simplified by dividing the integrations into the liquid-vapor interfacial and fluid-structure conponents. On fluid-structure surfaces, the velocities are assumed to be zero. Also, on the liquid-vapor interfaces, the fluctuating velocity Vk' is assumed to be zero. Using these assumptions the area integrals can be rearranged as follows 3- f 1 -H + + -- + (

(V -)-n -

n] dA W A k k k k +- (<> V 'V- V )on- P n~dA W A k k *i -- k k (.0 =- [

(V- V )-ndA-- n dA WJA k k i- k W A. 'k' - -A < PI > n dA The first and second integrals of the right hand side representithe momentum exchange across the liquid-vapor interface and are replaced with the term .ik The third integral represents the fluid-structure inter- actions, or wall friction, and is replaced by F .wk By applying Eq. (A.6) with B-1, the fourth term can be written as -259- A > ndA = < P -k < > k (A.31) Combining these expressions, the momentum equations are given by (a< >V> + Vs(0k

) + ak SKPk> (A.32) -Fwk - F ik +cak g-Mk Again by dropping the time and space average notation identifying K. ' as the momentum exchange due to turbulence, F the momentum K tkI equation may be written in THERMIT notation as tP 'k) + V(czkPkVk Vk) + ctk VPkat k k + k k k k + k k (A.33) -30.+ . . --wk i - Fk akpg- Ftk As discussed in Reference 3, Eq. (A.33) represents the conservative form of the momentum equation, while the non-conservative form is used in THERMIT. To obtain the non-conservative form, the first two terms of Eq. (A.33) are differentiated by parts and then simplified using the mass Eq. (A.25). Performing this operation one obtains k + + + + akPk -5;-+ akk VkV Vk + ak k = wk - Fik + akk (A.34) 4. 4 - Ftk -(r -W)V tkk tw kc -260- It is now assumed that the momentum exchange due to mass exchange can be included in the interfacial momentum exchange term, F , and that the pressure, Pk' is the same for each phase (i.e., P V =PjP). Then, the vapor and liquid momentum equations may be written as Vapor Momentum Equation 3V 4.++ ap --- + a VVV +aVP=-F -iv +cap g -F (A.35) v t v V V wv iv v tv Liquid Momentum Equation (U-a)p + (1-a)p * -VV + (l-cX)VP (A.36) 4. ++ + =-F9 - Fi + (l-c)pzg - F. , A.6 Conservation of Energy The conservation of energy equation is obtained by using Tk= ek' + 4. . =4 $k= 0 and k =qk+ Pk k -Tk Vk Once again the viscous terms are neglected so that r = 0. Substituting these terms into Eq. (A.23) yields at (a k) + Ve(cak(Pk> +)) 1 (---4+ + + -- + 'd - W - Vk)-ndA (A.37) 1 A +1++ - + + e 'V

can be neglected due to its small value compared to the other terms. Next the area integrals are divided into liquid-vapor and fluid-solid components. As in the case of the momentum equations, it is assumed that the velocity on fluid-solid interfaces is zero. The integrals can then be rearranged and identified as follows: 1 + + + -A+ k -kk k kkik (A.38)fA lJ -- + + 1--r + + -- Oc

V - n WdA=-k

V -ndA =

3-- (A.39) A k k W. k jk kat WJ -qk-n dA % k (A.40) Combining these expressions the energy conservation equation may be written as -26 2- (akp+ < e> +

Va kVk) + czk V. (<_k 'k'> (A.46) This term represents the mass exchange due to both temporal and spatial fluctuations and, hence, it is appropriate to associate this term with the turbulent mass exchange rate. In THERMIT this term is approximated as Wtk =kk i kkkj)) (A. 47 ) where E/1 is the turbulent velocity (see Section 3.3 for more details). In the momentum equation, the first assumption is that the viscous shear forces can be neglected. This assumption is valid for reactor conditions because the viscous force is small compared to the other forces. The second assumption is that the pressure is the same for each phase -264- (i.e., P=P =Pv). This assumption is appropriate provided the control volumes are not too large. Assumptions concerning the integrals and fluctuating terms have also been made. The wall friction term and interfacial friction term are: F kW J n dA (A.48) s + 1 -- + + + -- Fk k k[Pk>V>(V-V)n - Pk ndA .WAi (A.49) + (I -Wtk) The final term which is identified is the F . In THERMIT, only the Z-direction component is included. As in thetv case of the mass equation, the temporal and spatial fluctuating terms give rise to the turbulence effects. This term can be written as Ftk k k kV( Vk)+ ak V* ) (A.50) In THERMIT, this term is approximated by Ftk = -(ickk i - (ckk V )) (A.51) In the energy equation it is again assumed that viscous effects can be neglected and that the pressure is the same for each phase with a control volume. It is also assumed that heat conduction between channels and the work dissipation term can be neglected due to their relatively small values. -265- Assumptions concerning the various integrals and fluctuating terms have also been made. The wall heat transfer term has been associated with the fluid-solid area integral of the heat conduction term; Jk A qk' n dA (A.52) The energy exchange between phases due to mass transfer and conduction is associated with the interfacial energy exchange rate; Q= ((V- - k k Vk)-n dA (A.53) The final term to be identified is the energy exchange due to turbulence. This term is associated with the fluctuating terms and may be written as Qtk kk&kVk + c k V () (A.54) where fluctuations in both the density and internal energy have been neglected. In THERMIT this term is approximated as Qtk kj k ek i k k ekj(A.55) By using all of the above assumptions, the THERMIT conservation equations have been obtained from the local, instantaneous balance equations. This discussion has attempted to identify the major simplifying assumptions to obtain these equations. The form of the -266- equations used in THERMIT can now be understood in terms of their origin and restrictions caused by neglecting certain phenomena. -267- Appendix B Two-Phase Mixinz Model Assessment Results The tabulated results of the two-phase mixing model assessment effort are presented in this Appendix. The test conditions for the various experiments used in this effort are listed in Table B.1. The G.E. isothermal test comparisons are presented in Table B.2. For the heated G.E. tests, Table B.3 contains uniformly heated test comparisons while the non-uniformly heated test comparisons are given in Table B.4. Finally, the Ispra BWR test comparisons are presented in Table B.5. For each of these tables the measured and predicted exit mass velocity and quality (except for the isothermal tests) distributions are given. The bundle average mass velocity and exit quality are also listed. For certain cases, the code predictions without the two-phase mixing model are given as a means of comparison. Further details of these comparisons have been given in Chapter 5. -268-- TABLE B.1 Test Conditions for Rod-Bundle Experiments G.E. 9-Rod Ispra 16-Rod BWR Ispra 16-Rod PWR P (MFa) G (kg/m2s) q" (MW/M2 Lib (kJ/kg) sub X % out D (mm) e Length (m) Spacer Type Radial Power D .stribution 6.9 650 to 2200 0.71 to 2.1 67 to 525 3 to 22 12.1 1.83 Pin Uniform and Non-Uniform 7.0 1000 to 2000 0.12 to 0.77 30 to 180 2 to 31 13.3 3.66 Grid Uniform 16.0 2500 to 3500 0.07 to 0.11 250 to 400 -20 to 20 10.7 3.66 Grid Uniform -269- TABLE B.2 comparison of Measured and Predicted Exit Mass Velocities For Isothermal Tests in 9 Rod G.E. Tests IC Case G Mass Flow G1 G2 G3 Number (kg/m2s) Error (%) (kg/m2 2s)-g/r *s) (kg/m2* lB Data 650 -1.6 422 627 713 THERMIT 454 605 755 THERMIT (No Mixing) 440 584 736 IC Data 1343 +0.7 951 1274 1560 THERMIT 1003 1261 1525 THERMIT (No Mixing) 926 1259 1607 ID Data 2048 +0.46 1485 1954 2292 THERMIT 1535 1923 2323 THERMIT (No Mixing) 1399 1906 2435 1 E Data 2672 +1.06 2197 2591 2970 THERMIT 2010 2518 3035 THERMIT (No Mixing) 1829 2491 3182 TABLE B.3 Comparison of Measured and Predicted Exit Quality -and Mass Velocity Distributions for Uniformly Heated 9-Rod G.E. Cases Case G Xout Mass Quality X1 I x2 Ix3 G1 G2 G3 Number (kg/m2*s) Error (%) Error (kg/m2-s) (kg/m -s) (kg/m2S) 2B2 Data 719 0.029 -1.5 -0.01 0.003 0.014 0.03 505 707 732 THERMIT 0.037 0.025 0.028 550 690 809 2B4 Data 726 0.176 -0.0 +0.015 0.133 0.18 0.22 711 701 760 THERMIT 0.109 0.156 0.21 617 701 800 2C1 Data 1438 0.042 -0.05 -0.006 0.029 0.018 0.059 1309 1446 1461 THERMIT 0.037 0.039 0.046 1161 1391 1571 I TABLE B.3 (continued) F r ___________ r v -1 G2(kg/rn2aS) out Mass Error (%) I 4- 4- 1 1449 0.075 Quality Error Xl x2 X3 91 2 (kg//m2s) G2 (kg/m*s) G32 (kg/m2s) - 4I 4- - + * +0.05 +0.0091 0.063 0.056. 0.075 0.069 0.10 0,084 1313 1221 1394 1411 1552 1572 2D1 IData 732 0.110 +0.74 -0.002 0.083 0.105 0.117 576 760 754 THERMIT 0.091 0.101 0.121 602 713 798 2D3 Data 732 0.318 -0.03 +0.019 0.26 0.33 0.36 665 722 764 THERMIT 0.19 0.28 0.38 603 713 802 2El Data 1465 0.035 2.8 0.0 0.004 0.025 0.05 1288 1495 1576 THERMIT 0.039 0.03 0.039 1250 1417 1567 THERMIT (No Mixing) 0.20 0.02 0.032 695 1485 1655 Case Number 2C2 Data THERMIT TABLE B.3 (continued) Case ( G 2 s out Mass Quality X 1 X2 X3 G1 G2 G3 Number (kg/m-S) Error (%) Error (kg/m s) (kg/m *s) (kg/M -s) 2E2 Data 1465 0.106 3.24 -0.007!0.04910.t97 0.105 1418 1462 1600 THERMIT ;0.07610.096 0.11 1259 1431 1559 THERMIT (No Mixing) 0.30 0.081 0.106 783 1477 1646 2E3 Data 1438 0.215 2.6 -0.007:0.16 0.18 0.25 1309 1466 1527 THERMIT '0.13810.19 0.25 1219 1392 1554 THERMIT (No Mixing) 10.48 0.20 0.19 855 1370 1687 2GI Data 1451 0.038 -4.9 0.003 0.031! 0.044 0.042 1196 1313 1549 THERMIT 0.049' 0.033 0.041 1142 1408 1570 THERMIT (No Mixing) 0.29 0.017 0.032 640 1493 1632 N) N) TABLE B.3 (continued) Case G out Mass Quality X X2 X3 G G2 G3 Number (kg/m 2s) Error (%) Error (kg/m2.s) (kg/m2. s) (kg/2. 2G2 Data 1465 0.09 2.5 -0.008 0.02 0.068 0.110 1356 1507 1533 THERMIT 0.075 0.084 0.10 1241 1432 1563 THERMIT (No Mixing) 0.37 0.06 0.089 721 1515 1620 2G3 I Data 1451 0.16 4.1 -0.009 0.074 0.127 0.1761 1173 1535 1573 THERMIT 0.11 0.146 0.185 1273 1411 1549 THERMIT (No Mixing) . - 0.46 0.12 0.165 799 1482 1604 N) TABLE B.4 Comparison of Measured and Predicted Quality and Mass Velocity Distributions for G.E. Non-Uniformly Heated Cases Case G X I X2 X3 4 5 G G2 G3 G4 G5 Numberkk_(k 3B2 Data 726 0.032 0.08 0.042 0.108 -.0431 .009 543 688 753 685 434 THERMIT 0.064 0.057 0.058 -.015 -.017 538 673 802 696 564 THERMIT (No Mixing) 0.192 0.059 0.06 -.03 -.02 409 802 840 653 518 3Dl Data 739 0.084 0.123 -- -- -.037 .0241 437 -- -- 852 454 THERMIT 0.113 0.124 0.135 +.009 -.005 547 654 756 791 682 THERMIT (No Mixing) 0.49 0.132 0.127 -.0091 .012 384 638 786 797 615 1%) -Is TABLE B.4 (continued) Case G X X X X X X G G G G0 1 2 3 4 5 1 2 C3 4S Number kAkAkgL . &k.k 2 2 2 2 2 2 m -s M .s m .s m -s Im.s 1M-S 3E1 Data 1465 0.035 0.106 -- 0.163 -. 036 .002 1077 -- 1156 1255 1940 THERMIT 0.0770.078 0.076 -.017 -.02 986 1177 1381 1625 1815 THERMIT (No Mixing) 0.44 0.095 0.08 -.016 -.038 555 1091 1363 1467 1960 3E2 Data 1438 0.10 0.16 0.167 0.227 .034 .075 1085 1024 1207 2000 1275 THERMIT 0.099 0.122 0.143 .05 .029! 1137 1293 1441 1584 1469 THERMIT (No Mixing) 0.55 0.183 0.16 .01 .036 648 1102 1368 2015 1307 ~1 -J Lri TABLE B.5 Comparison of Measured and Predicted Quality and Mass Velocity Distributions for Ispra 16-Rod BWR Cases Case G X X X2 X X5 G G2 G4 G5 Number (kg/m2s) (kg/m*s) (kg/m2.s) (kg/m2.s) (kg/mS) 130.3 Data 998 .085 .077 .04 .08 .086 995 724 1007 978 THERMIT .073 .054 .094 .092 962 882 1038 1026 THERMIT (No Mixing) .075 .11 .086 .085 911 663 1110 1111 131.2 Data 1000 .148 .128 .066,.144 .164 979 667 979 983 THERMIT .13 .0991.178 .168 944 888 1023 1008 THERMIT (No Mixing) .13 .22 1.146 .144 899 573 1095 1095 Ch3 I' TABLE B.5 (continued) Case2 X0 1 2 4 5 G G2 G4 G5 Number (kg/m2. s) (kg/M2. s) (kg/m *s)g (kg/m -s) (kg/m2S) 107.3 Data' ~1017 -.15 .~142 .0671.172 .152 987 686 1047 1012 THERMIT .134 .103'.18 .173 972 917 1055 1041 THERMIT (No Mixing) .139 .23 .152 .152 936 595 1135 1134 99.3 Data 1000 .219 .176 .083j.189 .22 1027 648 1037 973 THERMIT .171 .127'.24 .23 960 891 1055 i 1043 THERMIT (No Mixing .21 .51 .18 .19 886 537 1152 1136 109.6 Data 999 .284 .26 .1771.30 .30 947 629 1012 978 THERMIT .23 .158 .33 .32 971 883 1057 1 1048 THERMIT (No Mixing), .27 .59 .27 .27 910 575 1126 1118 ____ ___ ___ __'No_____ ___ig)'__ 126_ 1118 TABLE B.5 (continued) Case G X X X X X . G G2 G 50 1 2 4 5 1 25 Number (kg/m2 *s) (kg/m2 -*s) (kg/m2 -s) (kg/m2 -s) (kg/m2 *s) 125.4 Data 1524 .04 .045 021 .037 .05 1448 1048 1592 1489 THERMIT .037 .029 .045 .044 1464 1324 1604 1583 THERMIT (No Mixing) .033 .096 .043 .041 1424 841 1672 1690 124.4 Data 1520 .084 .086 .051 .08 .092 1408 1029 1710 1621 THERMIT .075 .058 .0941.092 1472 1365 1575 1559 THERMIT (No Mixing) .074 .199 .0841.086 1411 803 1686 1675 141.8 Data 1528 .144 .162 .117 .133 .166 1416 934 1572 1430 THERMIT .127 .101 .16 .158 1484 1405 1596 1575 THERMIT (No Mixing) .13 .45 .142 .144 1421 762 1693 1692 co I, TABLE B.5 (continued) Case G X X X X X G G G G 0 1 2 4 51 2 4 52 . 12 2 2 2Number (kg/m -s) (kg/rn -s) (kg/m *s) (kg/m *s) (kg/m *s) 118.2 Data 1984 .029 .028 .015 .039 .037 1901 1334 1975 1970 THERMIT .027 .021 .031 .0311 1894 1707 2086 2057 TTHERMIT (No Mixing) .027 .031 .03 .032 1832 1331 2195 2187 113.4 Data 1994 .083 .074 .038J.084 .0741 1869 1296 2064 2000 THERMIT .074 .058 .092 .089! 1939 1804 2069 2049 THERMIT (No Mixing) .072 .22 .083 .0841 1851 1015 2196 2187 115.4 Data 1976 .122 .139 .075 .106 .106 1845 1239 2152 2039 THERMIT .101 .078 .124 .121 1938 1821 2051 2033 THERMIT (No Mixing)' .096 .29 .113 .118j 1869 1016 2195 2162 Ia -280- Appendix C Heat Transfer Correlations The heat transfer model in THERMIT has been discussed in Chapter 6. In this Appendix the heat transfer correlations are presented. A summary of the correlations used in the model is given in Table C.1. The logic for using these correlations is presented in Figure 6.3. The correlations which follows and their associated logic system are used in the current version of THERMIT. Sieder-Tate (vapor or liquid) hT 0.023 k Re'8 Pr3 (p/Y )4(C.1)ST Re P (Fluid properties at bulk fluid temperature, except P1w at T ) McAdams (vapor or liquid) h = 0.13 k [p2 g (T-T) Pr / P 3 (C.2) (Fluid properties should be at a fluid film temperature; T is either T or T ) Chen Chen = hfc (TW - T ) + hnbT w - Tsat) (C.3) k .8 .4 h 0.023-i"kRe Pr'o F (0.4fc D k k -281- TABLE C.1 Heat Transfer Correlations Regime Correlation 1. forced convection to Sieder-Tate single-phase liquid 2. natural convection to McAdams single-phase liquid 3. subcooled boiling Chen 4. nucleate boiling Chen 5. transition interpolation between qCHFand qMSFB 6. forced convection to Sieder-Tate single-phase vapor 7. natural convection to McAdams single-phase vapor Regime Checkpoints Correlation 1. Critical Heat Flux Biasi, W-3, CISE-4 2. Minimum Stable Film Boiling Henry Temperature -282- [kcj.5 r m h = 0.00122 S nb 1 Pr-.29 .25Pri 2.35 (X~ + .213).736tt (C.6) 0.1 (plp v ) 5 v /P '=x/(1-x) [1 + .12Re5 I-4TP [1 + .42 Re* 78 0.1 Re < 32.5 32.5 < ReTP < 70 Re TP >70 -4 1.25 D / pRe = 10 F (1-ct) pk DY, lJk Biasi (Critical Heat Flux) Use the second expression below for G < 300 kg/m2-sec; for higher G, use the larger of the two values: = 2.764 x 10 (100D)-n C1/6 [1. 468 F (P.bar) C1/6-xl (C.10) qBiasi = 15.048 x 10 (100 D)-n C-6 H (Pbar) [1-x] F(Pbar) = 0.7249 + 0.099 P bar exp (-0.032 P bar) (C.5) < 0.1 X- xtt (0.7) (0.8) (0.9) (C.11) (C.12) a _ 75 c p.(T w-T sat .24 ( Pw -P) h p -283- H(Pbar) = -1.159 + 0.149 Pbarxp(-.019 P bar + 9 P (10 + P2 -1 (C.13)bar bar Note: P =bar 1- p n = .4 D < .01 m .6 D > .01 m (C.14) W-3 (Critical Heat Flux) The W-3 correlation for a uniform heat flux is a function of the pressure, hydraulic diameter, mass flux, quality and inlet subcooling. This correlation is given as U CHF {(2.022 - 0.0004302 P) + (0.1722 - 0.0000984 P) 106 x exp [(18.177 - 0.004129 P) X]} x [(0.1484 - 1.596 X + 0.1729X JXI) (G/106) + 1.037] x (1.157 - 0.869 X) x [0.2664 + 0.8357 exp(-3.151 DE)] x [0.8258 + 0.000794 (h - hIN)] Btu/hr ft2 (C.15) For a non-uniform heat flux the critical heat flux is given as NU U /F (C.16) CHF CHF -284- where and C 9 - exp(-cZ)) 0 f [q(z) exp(-c(Z - z))]dz 0.15 (1 - X ) 4.31 C 6 C -1 (G/1O6 )7 (C.17) (C.18) k is the channel location and Xc and q(Z) are the quality and heat flux at the location Z. It should be noted that English units are used in this correlation. CISE-4 (Critical Heat Flux) The CISE-4 correlation, which is used to calculate the critical power, depends on the pressure, hydraulic diameter, mass flux and boiling length. The correlation is given as (0.19) QC aLb GAh b + Lb where a (1+1.481 x 10 (1 - P/PC) G) a = (1 - P/Pc) (G/1000)1/3 if C < G* if G > G* (C.20) (C.21) -285- b - 0.199 (P/P - 1)a4 G Dh14 Ch G*= 3375 (1-P/Ps C3 Lb-bboiling length (m) Q-critical power (w) Pc - critical Pressure (Pa) CHF-Void q CHF .1178 (1-a) hv g ~g v) Pv21025 Minimum Stable Film Boiling Temperature TMSFB HN + HN a T) [(pkc ) /I(Pkep)w .5- (P) T - 581.5 + .01876 (P - 010 630.39 + .004321 (P - PO 5 Po 68.95 x 105 Pa (P)- 0 1273 - 26.37-10-5 P Note: (pkcP)w above refers to properties clad surface material properties. P < 4.8260105 Pa (C.27) P > 4.826-105 Pa of the wall itself, i.e. (C.22) (C.23) (C.24) (C.25) (C.26) P < P0 -0 P > P0