A Condensed Review of Nuclear Reactor Thermal-Hydraulic Computer Codes for Two-Phase Flow Analysis by M. Kazimi* and M. Massoud Energy Laboratory Report No. MIT-EL 79-018 February 1980 *Associate Professor, Department of Nuclear Engineering x\ A CONDENSED REVIEW OF NUCLEAR REACTOR THERMAL-HYDRAULIC COMPUTER CODES FOR TWO-PHASE FLOW ANALYSIS by M. Kazimi and M. Massoud Nuclear Engineering Department and Energy Laboratory Massachusetts Institute of Technology Cambridge, Massachusetts 02139 Written: Published: April 1979 February 1980 Sponsored by Boston Edison Company Consumers Power Company Northeast Utilities Service Company Public Service Electric and Gas Company Yankee Atomic Electric Company under MIT Energy Laboratory Electric Utility Program Energy Laboratory Report No. MIT-EL 79-018 -2- ABSTRACT A review is made of the computer codes developed in the U.S. for thermal-hydraulic analysis of nuclear reactors. The intention of this review is to compare these codes on the basis of their numerical method and physical models with particular attention to the two-phase flow and heat transfer characteristics. A chronology of the most documented codes such as COBRA and RELAP is given. The features of the recent codes as RETRAN, TRAC and THERMIT are also reviewed. The range of application as well as limitations of the various codes are discussed. i -3- TABLE OF CONTENTS Page Abstract 2 Table of Contents 3 Abbreviations 5 List of Figures 6 List of Tables 7 1. Introduction 8 2. Classification of Nuclear Reactor ThermalHydraulic Codes 2.1 2.2 Classification According to System Analysis Capability 10 11 2.1.1 Component Codes 11 2.1.2 Loop Codes 17 Classification According to Two-Phase Model 2.2.1 The Homogeneous Equilibrium Model 2.2.1.1 Approximations to the Field 18 18 20 Equations - Component Codes 2.2.1.2 Approximations to the Field 24 Equations - Loop Codes 2.2.2 2.2.3 2.3 Improvements on the Homogeneous Equilibrium Model 33 2.2.2.1 Dynamic Slip Model 34 2.2.2.2 Drift Flux Model 37 The Two-Fluid Model 40 Classification According to Range of Application 47 2.3.1 47 LOCA Codes -4Page 3. Evaluation Model Codes 2.3.1.2 Best Estimate Codes Two-Phase Heat Transfer Models 3.1 4. 2.3.1.1 ........ 48 ........... 49 ........... ........ 55 Heat Transfer Regimes and Correlations Fuel Rod Models ....... 57 ................................. 66 ................................ 70 ............................... 74 Clad Region .................................. 76 4.1 Fuel Region 4.2 Fuel-Clad Gap 4.3 5. Numerical Methods 6. Summary and Conclusions 78 . ........................... 85 .......... ........................... 85 6.1 Summary 6.2 Conclusions ........................... 6.2.1 Component Co des 6.2.2 Loop Codes . 100 ...................... ,ele e e e. e e e e l el 100 . 10 5 References 107 Appendix 1 112 Appendix 2 117 Appendix 2 - References ................................ 127 REPORTS IN REACTOR THERMAL HYDRAULICS RELATED TO THE MIT ENERGY LABORATORY ELECTRIC POWER PROGRAM .................129 .. -5- Abbreviations The following abbreviations are referenced in this report: ATWS: Anticipated Transients Without Scram COBRA: Coolant Boiling in Rod Arrays DBA: Drift Flux Model DSM: Dynamic Slip Model EPRI: Electric Power Research Institute FLECHT: Full Length Emergency Cooling Heat Transfer HEM: Homogeneous Equilibrium Model LOFT: Loss of Flow Transient LOCA: Loss of Coolent Accident MWR: Method of Weighted Residuals NSSS: Nuclear Steam Supply System RIAs: Reactivity Insertion Accidents RETRAN: RELAP4 - TRANsient TRAC: Transient Reactor Analysis Code UHl: Upper Heat Injection WREM: Water Reactor Evaluation Model -6- LIST OF FIGURES No. Page 1 Coolant Centered and Rod Centered 13 2 Lateral Heat Conduction 16 3 Control Volumes for Subchannel Analysis 21 4 Geometry for Mass, Momentum and Energy Equations 26 5 Flow Path Control Volume 28 6 Lumped Parameter Simplification of Multidimensional Flow Regimes 28 7 A Downcomer Representation in RELAP4 31 8 Schematic of RELAP4 Model of a Large PWR 32 9 PWR LOCA Analysis 50 10 BWR LOCA Analysis 51 11 Heat Transfer Regimes Traversed in Blowdown 54 12 Flow and Heat Transfer Regimes in Rod Array with Vertical Flow 55 13 Heat Transfer Regimes at Fixed Pressure 56 14 Effect of Quality on Calculated Boiling Curve Predicted by BEEST Package 61 -7- LIST OF TABLES Page No. 9 1.1 List of Reviewed Thermal-Hydraulic Codes 2.1 Cases Where 1-D HEM is not Acceptable 40 2.2 Two-Phase and Single-Phase Comparison with Respect to the Flow Equations 42 2.3 Two-Phase Flow Model Description 45 2.4 ATWS in BWRs and PWRs 53 3.1 Heat Transfer Regime Selection Logic 67 3.2 Heat Transfer Correlations - Pre-CHF 68 3.3 Heat Transfer Correlations - Post-CHF 69 4.1 Approximations Made to the Equation of Heat Conduction 75 6.1 Component and Loop Codes Comparisons 86 6.2 Models and Methods Used in Some ThermalHydraulic Codes 89 6.3 99 -8- 1. Introduction Numerous computer codes have been written to calculate the thermal-hydraulic characteristics of the reactor core and the primary loop under steady-state and operational transient conditions as well as hypothetical accidents. some of these codes are still to come. New versions of The main purposes of the continuing effort in the development of such computer codes have been improved computational effectiveness and improved ability to predict the response of the core and the primary loop. Therefore, efforts have been continued to incorporate the recent models and methods of analysis in the areas of both hydrodynamics and heat transfer in two-phase flow to the extent that their prediction are reasonably reliable. For example, such a step by step development has been effected in the various versions of COBRA and RELAP Computer Programs. The code users are therefore confronted with the need to develop criteria to choose the most appropriate version to handle a specified case. This is a two pronged decision since it requires not only an evaluation of the models and methods used in each code but also a comparison between the results and experimental data to observe how well these data are predicted. An attempt is made here to address the first step, i.e., comparison of the models and methods. To accomplish this, a study was made on the physical models and numerical methods which have been employed in the WOSUB, RETRAN, TRAC and THERMIT as well as various versions of COBRA and RELAP as listed in Table 1-1. -9- Table 1-1 List of reviewed thermal hydraulic codes. Name of Code Reference Number COBRA-I 1 COBRA- I I 2 COBRA-III 3 COBRA-IIIC 4 COBRA-IIIP 5 COBRA-IV-I 6 COBRA-DF 7 COBRA-TF 8 RELAP2 9 RELAP3 10 RELAP3B-MOD101 11 RELAP4 12 RELAP4-MOD5 13 RELAP4-MOD6 14 RELAP4-MOD7 14 RELAP 4-EM 15 RELAP5 16 WOSUB 17 RETRAN 18 TRAC 19 THERMIT 20 -10- These codes, especially COBRA and RELAP series, are well known thermal hydraulic computer codes and have been extensively used in the nuclear industry. WOSUB and RETRAN introduce a new treatment for the hydrodynamics modeling. TRAC and THERMIT have gone further by applying the most advanced existing treat- ment of the two-phase flow, namely, three-dimensional, twofluid, non-equilibrium model. In the comparison that follows, both the advantages and drawbacks are noted in each code and ultimately it is attempted to assess the capability of each code for handling a specified case. 2. Classification of Nuclear Reactor Thermal-Hydraulic Codes The existing thermal hydraulic codes may be classified under several categories as follows: 1) Capability of the system analysis This contains two different classes of codes, namely, system component codes and loop codes. Basically, the hot channel or the fuel behavior codes are system component codes; however, some of these codes are extended to other situations far removed from subchannel (one channel) geometry. Integration of the down comer, jet pumps (in BWR's), bottom flooding, UHI and the like models into a component codes, makes itta vessel code. As distinct from the loop codes which are devised to analyze the whole primary side including reactor core and the secondary side, a variety of codes ranging from hot channel to vessel codes are called system component codes in this report. -11- 2) Type of two-phase flow modeling This part deals with the mathematical models used in thermal hydraulic codes to calculate the characteristics of the two-phase flow either in the reactor core or in the primary The two pertinent methods in this respect, namely, loop. the homogeneous equilibrium model and the two-fluid model fall in this category. 3) Range of application Since the capability of each code to handle flow and fuel rod calculations depends upon the mathematical models used to represent the physical situations as well as the numerical methods employed, codes can be classified in these respects into steady-state, transient and accident analysis (such as LOCA) codes. Naturally, the more demanding codes in this respect are ATWS and LOCA codes. 4) Type of application Codes may also be classified based upon their types, i.e., Best Estimate (BE) type and Evaluation Model (EM) type. The latter group are basically devised for the purpose of licensing. The type of nuclear reactor for which thermal hydraulic codes are devised category. (such as PWR, BWR and LMFBR) may be another A detailed discussion concerning each mentioned category is presented in the following sections. 2.1 Classification According to System Analysis Capability 2.1.1 Component Codes Core thermal hydraulic assessments necessitates analysis -12- of fluid passing axially along the parallel rod arrays. Such analysis is difficult to conduct due to the degree of freedom associated with parallel rod array and the two-phase flow and heat transfer involved in nuclear reactors. In addition, radial and axial variations of the fuel rod power generation exacerbates this situation. Assumptions have been made to simplify the task of modeling the hydrodynamics and heat transfer characteristics of the rod arrays. Generally, there are three pertinent methods ( 2 1 )used in rod bundle thermal hydraulic analysis of the nuclear reactor core as well as heat exchangers, namely, (a)-subchannel analysis, (b)-porosity and distributed resistance approach and finally (c)-benchmark rod-bundle analysis which uses a boundary fitted coordinate system. The first approach is widely used in the subchannel codes such as COBRA, FLICA, HAMBO and THINC. Whereas the second approach is employed in THERMIT. The subchannel approach will be more elaborated upon here, while a discussion in detail of these three concepts is presented in Ref 21. In the subchannel approach, the rod array is considered to be subdivided into a number of parallel interacting flow subchannels between the rods. The fluid enthalpy and mass velocity is then found by solving the field or conservation equations for the control volume taken around the subchannel. Although a rod-centered system with subchannel boundaries -13I 1 G I - - 1 7 1 ;12 I -O 6 1 I I 8 3 I - -.- 9 I \ 10 }15 14 0-00 ! 7 8 1 9 - axis of sy'rmetry 1 10 'C> IC> 1~ fY61I wS I I I I :/ (a) - ,- IIr_ r - 'r -, 'iIIU',1 Q9V -- I _ - _ I - - ) I QII . -. J -t - ) 0.'1 ('):0,0 - z- , \\I --------- i L _ _ _ _ (b) Fig. 1 [22] (a) Coolant centered subchannel and conventional subchannel numbering. (b) Rod center system with subchannel boundaries. I -14- defined by lines of "zero-shear stress" between rods (Fig. l-b) seems to be a well-defined control volumes, it has become customary to consider a coolant centered subchannel as a control volume (Fig. -a). The number of the above-mentioned control volumes axially is as many as the number of the channel length intervals. Unlike the benchmark rod-bundles approach, the subchannel approach does not take into account the fine structure of both ** velocity and temperature within a subchannel. In other words, there are no radial gradients of flow and enthalpy in the subchannels but only across subchannel boundaries. There- fore, the flow parameters such as velocity, void fraction, and temperature are averaged over the subchannel area. Further- more, the averaged values are assumed to be located at the subchannel centroid. The following example elaborates the latter assumption. (48) The transverse heat conduction in the fluid passing through the subchannels shown in Fig. 2-a becomes T. - T. q"ij= k. 1 [1.a3 13 * This model was first introduced in the Italian subchannel, code (CISE (23). It is especially preferred in modeling the strict annual two-phase flow condition, due to its resemblance to the annual geometry. ** An excellent discussion concerning the fine structure of the flow field within the coolant region is presented in Ref. (24). -15- and T. q jk= kjk - T J Il.b] 1jj ljk where q", k, T and 1 are heat flux, thermal conductivity, averaged temperature and finally centroid-to-centroid distance between the adjacent subchannels respectively. Assuming identical fuel rods, the centroid located averaged subchannel temperature seems to be a valid assumption for subchannel j. However, for subchannels i and k, it is expected that the averaged temperatures are located closer to the gap 1 and gap m respectively. shown in Fig. 2-b. This is also the case for the temperatures The centroid located averaged temperature is a valid assumption for low conductivity coolants and high P ratios, whereas, it is a dramatic assumption for high D * conductivity coolants and tight rod bundles. ___ * A discussion in detail and a suggested method to correct the centroid located averaged values are presented in reference 25. -16- -I Gap 1 I i I I * I I I I I TI , ij ibchannel =ntroid iJ i li_ . . ljk-. , 13 (a) - ( " I a I - (b) Fig. 2 Lateral heat condition. -1 2.1.2 /- Loop Codes Analysis of the whole primary system during transient conditions and hypothetical accidents such as loss of coolant, pump failure and nuclear excursions, necessitates modeling the whole loop components such as pipes, pressurizer (in PWR's), pumps, steam generator, jet pumps (in BWR's), valves and reactor vessel. Also, the effects of the secondary system need to be considered. The thermal hydraulic behavior of the reactor core during the course of a transient is tied to the core nuclear characteristics through the reactor kinetics. HIence, the reactivity feedback should be considered in the process of the primary loop modeling. The RELAP series of computer programs are the well-known transient loop codes which have been extensively used in the nuclear industry. These codes are basically devised to analyze transients and hypothetical accidents in the nuclear reactor loop of LWR's and mainly consist of four major parts as follows: (1) a thermal hydraulic loop part, (2) a thermal hydraulic core part, (3) a heat conduction part, (4) a nuclear part. In these codes, the primary system is divided into volumes and junctions. The fluid volumes serve as control volumes, describe plenums, reactor core, pressurizer, pumps and heat exchangers. Each connection between volumes may be specified as a normal junction, a leak or a fill junction. A fill junc- tion as its name implies, injects water into a well-specified volume. By definition, volumes specify a region of fluid within -18- a given set of fixed boundaries, whereas junctions are the common flow areas of connected volumes.(10) Any fluid volumes may be associated with a heat source or a heat sink, such as fuel rods or the secondary side of a heat exchanger, respectively. While RELAP2 is able to handle only three control volumes with a fixed set of pipes connecting these volumes, representing the whole primary loop, RELAP3B and RELAP4 are capable of handling as many as 75 volumes and 100 junctions or even more, at the expense of more computer core. 2.2 Classification According to Two-Phase Model 2.2.1 Homogeneous Equilibrium Model Flow characteristics in component and loop codes are cal- culated through solving the field or conservation equations written for the well specified control volumes. The basic assumption made in modeling the two-phase flow is representing the two-phase by a pseudo single phase. This method of model- ing is also known as homogeneous equilibrium model (HEM). HEM is extensively used in the thermal hydraulic codes. The The homogeneous assumption implies that the phase velocities are equal and flow in the same direction, also the phase distribution is uniform throughout the control volumes. The equilibrium assumption requires the phase to be at the same pressure and temperature. The one dimensional HEM codes use an approximate set of .*To do this, only the array sizes in the COMMON blocks should be increased. -19- field or conservation equations for the mixture in conjunction with the constitutive relations. The differential form of the conservation equations written for a mixture is as (26) follows: Local mixture continuity equation: [1] *m (PV ) = 0 Pm + 3t Local mixture momentum equations: + 7t PmVm v [m = Om [2] +vm m] (Vm m where the product V V gives an array of nine components. This product can be written as VmVm = (Vi)m (Vk)m (i,k The surface stress tensor, T = 1,2,3,) , is made up of the pressure and the normal and the shear stresses T=P where T m m I - m is the viscous stress tensor and T is a unit tensor. Local mixture energy equation: mt (U + 1/2 V -[V (q m - V ) + [VPm (U [T-V])] + 1/2 + pg-V + m V )V ] = m where qm is heat flux, Qm is the body heating term and Um is the internal energy. These balance equations need to be accompanied by the constitutive equations for Tm' qm, and Qm' the equation of state, and the mixture properties. [3] -20- 2.2.1.1 Approximation to the field equations Component Codes Approximations which are made in solving the conservation equations in the component codes using the homogenous equilibrium model for the two-phase flow will be discussed in this section. Basically, none of the existing subchannel codes use such a generalized three dimensional set of field equations as are given by Equations 1, 2 and 3. Rather simplifying assumptions are made in these equations. For example, in most of the COBRA versions, flow is assumed to have a predominantly axial direction and all the "lateral" flow is lumped into one lateral momentum equation. The reason for such treatment may be justified by considering the noneorthogonal characteristics of subchannel arrangement (Fig. 3) which do not allow treatment of the lateral or transverse momentum equations as rigorously as the axial momentum equation. It is assumed that the interaction between two adjacent subchannels in the transverse direction is through two distinct processes,* namely, diversion cross-flow and turbulent mixing. Axial turbulent mixing between nodes is ignored. The first process, diversion cross-flow is assumed to exist due to local transverse pressure difference in the adjacent subchannels. Such a process transfers mass, momentum *A more general classification is given in Ref. (27) and is referenced in the model making process of WOSUB (17). Also see Ref. (47) for basic notation in subchannel analysis. -21- Control volume for axial mom., equation Exact control volume for transverse mom. equi Fig. 3[21] / Control volume: subchannel ana .. S... -- -22- and energy with the assumption that the cross flow loses its sense of direction when it enters a subchannel(4) Unlike the HEM versions of COBRA, WOSUB which is essentially devised for analysis of ATWS in BWRs does not account for diversion cross flow. The second process, turbulent mixing is assumed to be caused by both pressure and flow fluctuation. In this process, no net mass transfers, only energy and momentum are involved. This is due to the assumption of the equi-mass model.* The magnitude of the turbulent mixing term is determined either by some correlations or by a physical model that includes empirical constant. All the COBRA versions account for a single phase turbulent mixing while the two phase turbulent mixing term was added in the versions following COBRA-II, since COBRA-I does not account for this term. It should be mentioned that forced flow mixing which is caused by some rod spacing methods such as a wire wrap or diverter vanes is taken into account, especially in those codes which are capable of analyzing fluid flow in LMFBRs such as COBRA-IIIC and COBRA-IV-I. Recently, a wire wrap model has been added (2 8 ) to COBRA-III-P which makes it capable of handling LMFBR flow analysis. *The equi-volume model which is based on the change of same volume of flow is used in the MIXER code. For further detail see Ref. (22). -23- The steady-state versions of COBRA, namely, COBRA-I, COBRA-II and COBRA-III do not have any model for forced diversion cross flow. A more complete form of transverse momentum equation is employed in COBRA-IIIC, COBRA-IIIP and COBRA-IV-I which includes the time and space acceleration of the diversion cross flow. As a correction to the homogenous flow assumption, a one dimensional slip flow model which accounts for nonequal phase velocities, is considered in all the COBRA series up to and including COBRA-IV-I. added to these codes. of COBRA-IV-I A subcooled void calculation is also However, COBRA-I and the explicit scheme (to be described) do not have a subcooled void option. In the COBRA codes, the energy equation has been further simplified by assuming the turbulent mixing and convection heat transfer as the unique mechanisms for internal energy exchange. In such treatment, it is assumed that(2 9 ) -- no heat is generated within the fluid, -- changes in kinetic energy is small, -- no work against the gravity field. Neglecting the time change of local pressure, , limits these codes to transients with times that are longer than the sonic propagation time through the channel. (4) Unlike the previous versions, the COBRA-IV-I momentum equations account for the momentum flux term. -24- Further simplifications to the axial momentum equation have been made by neglecting surface tension contribution. This requires equal phase pressures. This basic assumption in addition to the assumed equal phase temperatures are the result of the thermodynamic equilibrium assumption. 2.2.1.2 Approximations to the field equations -Loop Codes Assumptions made to solve the field equations in the loop codes using HEM are discussed here. Except for RELAP5, which is the latest publicly available version of RELAP series, the remaining versions use the HEM for their hydrodynamic modeling. Therefore, a set of conservation equations written for a mixture (Equations 1, 2, & 3) is applicable for theoretical considerations. For the practical purposes, approximations have been made to this generalized set. The RELAP codes, generally have a lumped parameter structure in which the spatial effects are integrated over the control volume for the conservations of mass and energy. For example, the mass balance in its differential form is p = V-(p) at [4] Integrating over the control volume -i v at dr = fjf V- (pV)d-r [5] Now applying the divergence theorem to the right hand side of Equation [5], we get: -25- JJ-T ftV.(PV')dT= fT pV v nds V v V or t n-.(pV)ds + d = 0 [6] V s Using M J as the existing mass in the volume J at time ti and considering the term -ni .(pSV)i = W. which is equal to the inflow from side "i" into volume J (Fig. 4), the mass balance reduces to [7] dt M dt 3 E = i=l W.. 13 Similarly, a simplified form of the energy equation which has been used in RELAP2 and various MOD's of RELAP3 is as follows: n dt= Wij. .h. .+Q i=l ij1 1 j [8 where U. is the internal energy of volume J, hij is the enthalpy of fluid flowing from side i into volume j and finally Qj is the heat input to volume j. The effect of kinetic, potential and frictional energies are neglected in Equation [8]. However, the RELAP4 energy equation accounts for kinetic and potential energy changes. Unlike the mass and energy equations, the momentum equation is written for a shifted control volume as shown in Fig. 4. This method minimizes the extrapolation of boundary conditions. The final form of the momentum equation used in various MOD's of RELAP3 is as follows: -26- A j + 1 Flow channel wall A WJ+ 1 = W L pi+l Momentum cell W = WKL Junction J J ---- X ~I a CJ-1 W lass and energy cell J -l- J-1 AK Fig. 4 Geometry for mass, momentum and energy quations. -27- 1 L c(L 144gc dW. d dt However, P P = i+lj + + AP kjWj Pdz 1 44 jI [9] Pj like the energy equation, an improved form of the momentum equation is implemented in RELAP4 which takes the form: (30) dW. It = I -1-> < L z L. dF - (P +P <-4--> <-5.-> < F L K I1 d(vW) A [10] o > <8 <--9 > -- > dWt dt in equa- It is assumed that the junction inertia term, I tion f> L. L. - 1 7--8- - F 6- KK - K soK < dP ) (P+P <-3- I . ) [10] or the corresponding term d in equation [9] represents the rate of change of momentum everywhere in the selected control volume J. the geometric "inertia" W. (Fig. 5) In equation 10, I is for the flow path and also, = flow rate in junction J, Pk = Pressure in volume K, Pkgj = gravity head contribution for volume K, F = friction terms, v = velocity, A = flow area. The significance of each term in equation [10] is as follows: Term 1 represents the rate of change of momentum, Term 2 and 4 represent the pressure drop between two volumes, Terms 3 and 5 represent gravity, Terms 7 and 8 represent the friction and pressure drop associated with expansion and contraction. Term 6 represents the fanning friction terms. -28- 1k I -7 I I w. Pk L r. . IA I k I IP 1AL WL I i L. 1 I Wk I ! Ko Volume L Volume K Flow path control volume. Fig. 5 L M t t ! K Fig. 626] - . J !- ' Lumped parameter simplification of multidimensional flow- regimes. -29Term 9 represents the momentum flux or spatial acceleration term. Comparing Equations 9, used in RELAP3, and 10, used in RELAP4, it is obvious that the major difference is inclusion of the momentum flux term in Equation 10. The momentum flux term for a homogeneous volume becomes L VLWL-WkV k = d(VW) k _ A _ A Although the inclusion of the momentum flux term has improved the momentum equation, there are some cases in the loop modeling that the code user has to ignore this term. (For example, when the given control volume, J, connects into more than two control volumes, as illustrated in Fig. 6; where the double-ended arrows indicate junction). A similar case might be encountered in the lower plenum modeling when the flow is highly three dimensional. In such cases, it is quite difficult to define the control volume boundaries for the momentum balance at each junction These examples clarify the inability of a lumped parameter approach to model multidimensional regions. This is also the case in calculating the friction factor and heat transfer coefficient which by definition are( 4 8 ) du f = q =-k dT dr dr )dz w) [11] and where w =h(Tw-Tb ) [12] w and u are the viscosity and velocity of the flow and subscripts w and b represent values evaluated at wall and bulk -30- respectively, It is clear from equations 11 and 12 that the derivatives are evaluated at the channel wall. However, since a lumped parameter approach doesn't account for velocity and temperature profile, therefore, the above mentioned derivatives do not make sense. It is this reason which necessitates an input specified friction fractor for the codes using this approach. The junction inertia term is another term in the momentum equation (Equations 9 and 10) which becomes rather ambiguous in modeling the complex geometries. The junction inertia arises from an approximation in the momentum equation to the temporal inertia term as follows:(30) J x2 where x 1 AXx)T A(x) dw dt _ ___ dx dx ldw 2 dx dx ( dw 1I I dt [13] = center of control volume 1 and x 2 = center of control volume 2, and I is the geometric inertia for the flow path defined as: I -(2 dx [14] 1 The geometric inertia for a homogeneous volume, Fig. 4, a b becomes I = 2A + 2A , however for complex geome,trics, K L the inertia term may be determined by using a simplified assumption. The basic assumption which is introduced in this -31- respect (30) is that the inertia of a junction is composed of two independent contributions, one from each connecting volume. For example, Fig. 7 could represent a downcomer region. If we assume that Junction 1 communicates primarily with Junction 3, then with respect to the mentioned basic assumption, the geometric inertia will be: Ij = Ijl + Ij2 + Ij3 [15] or Ij L1 + =2A 1 L2 A22 + L1 13 Where L 1 is the effective length of both Junctions 1 and 3 and the effective length of junction 2 is assumed to be 2L2. \ '- f' L1 I 0 , - - : I O = Junction number --- = Flow Path I I L2 I 1A I Fig. 7(30) - A Downcomer Representation in RELAP4 A schematic of RELAP4 model of a PWR is presented in Fig. 8. It illustrates the complexity of accounting for all volumes and junctions in a LWR plant. -32- Steam Stearn rs sexis Stwnim gensra second 0leam menerator .condary de side Fig. 8[30] Schematic of RELAP4 model of a large PWR. -33- 2.2.2 Improvements on the Homogeneous Equilibrium Model By retaining more field equations, a more realistic approach to analysis of severe transients has become possible in reenct codes. The increased number of field equations enable a code to analyse transients in which the situation is far beyond the capability of the rigid assumption of equal phase velocities and temperatures. In this respect, countercurrent two-phase flow and vapor-liquid phase separation during small break transients and emergency core coolant delivery are notable examples. Since in a non-homogenous flow slip exists between the two phases, there is a relative motion of one phase with respect to the other. This relative motion arises due to density and/or viscosity differences between phases where usually the less dense phase will flow at a higher local velocity than the more dense phase, except for the gravity dominated flow (1 8 ) The general effect of slip is to lower the void fraction below the homogeneous value. V The slip or hold up ratio, s = g should not be VL confused with the slip velocity VsL = VL - Vg, or drift velocity, a concept which is used in the drift flux model. Unlike the one-dimensional HEM, a non-homogeneous, one dimensional flow calculation for a two-phase flow in thermodynamic equilibrium, involves the solution of one -34- equation of state and five differential equations: a mixture energy equation and one continuity and one momentum equations for each phase as it is done in RELAP5. 2.2.2.1 Dynamic Slip Model Simplified assumptions have been made to reduce the number of conservation equations while retaining the improvements over HEM codes. This is done in WOSUB and COBRA-DF by using the concept of diffusion or drift flux model, and in RETRAN by introducing the dynamic slip model. RETRAN computer code is basically developed from the RELAP series of codes. It is a one-dimensional code which solves four field equations written for a fluid volume as follows: Mixture continuity, Momentum, Energy equation, and time dependent behaviour of the velocity difference, obtained by subtracting the momentum equations written for each phase. SD g This additional momentum equation reads: DX (1 DX AgL P9 BgL 1 g sL 1P X a Pk ° 1 cc9P [16] where V,p,a,P represent velocity, density, void fraction and pressure respectively. Also Ag L represents the surface area between vapor and liquid phase per unit volume and BgL represents the friction coefficient between vapor and liquid -35- phases. In deriving equation (16), the following assumptions have been made: 1) The wall friction is nearly equal for the two- phase. 2) The momentum exchange between phases due to mass exchange is small. In addition to inclusion equation (16) in RETRAN, some improvements have been made in the field equations used in RELAP4, as follows: 1) Additional term in the mixture momentum equation with respect to the momentum flux. Mixture momentum flux: ax [A (agPg (V2)g) + al a +a VsL (V2 ) 1] = l 9Plg A Additional Term 2) Additional term in energy equation which accounts for the time rate of change of kinetic energy, at [ U2 [pA(-)]. at 2 -36- 3) Using a flow regime dependent two-phase flow friction multiplier. 2.2.2.2 Drift Flux Model Unlike RETRAN, COBRA-DF which is a vessel code and uses the drift flux model, employs five field equations to determine phase enthalpy, density and velocity*. This code is used exclusively for examination of upper heat injection of water during a LOCA in a PWR. Vapor diffusion or drift flux model is another step toward modeling a non-homogeneous non-equilibrium flow. The basic concept in this model is to consider the mixture of the two-phase as a whole, rather than treating each phase separately. The DFM is more appropriate for the mixture where dynamics of two components are closely coupled, however, it is still adequate where the relatively large axial dimension of the systems gives sufficient interaction time (26) In this model in addition to the three field equations written for the mixture, there is a diffusion * THOR which is developed at BNL (31) uses the DFM and accounts for thermal non-equilibrium of the dispersed phase only. -37- equation written for the dispersed phase which reads: ap O -t [17] -+ gt + V(agpg Vm) = r - V (gpg Vgj) [17 where r is the phase change mass generation, Vgj and Vm are g gm the drift velocity and mixture velocity respectively. WOSUB which is a BWR rod bundle computer code uses the DFM and solves four field equations written for a subchannel control volume, as follows: 1) continuity equation for mixture 2) continuity equation for vapor. at (Pg9ai)i+A where J g qgi (PgJg) i = Ap gii This equation reads: + Pgi qgi [18] = vapor flux = vapor volume flow to subchannel i T. = vapor volume generation in subchannel i per unit volume Equation 18] indicates the fact that the temporal and spatial increase in the mass of vapor in subchannel i is due to vapor generation in the subchannel and vapor addition from the adjacent subchannels. vapor volume generation term, The , appears in Equation [16] due to using the DFM. This term is part of the code constitutive package. It is modeled in WOSUB based on the Bowring's equation which relates -38- Y to the heat flux: i-TAp v where T h fg q [19] is a coefficient depending upon coolant condition, Ph is the heated perimeter and q" is heat flux in the fully developed nucleate boiling region, and A is the flow area. 'The effect of subcooled boiling non-equilibrium condition is considered in the final form of 3) . Mixture axial momentum equation: This is the only momentum equation considered in the code. Therefore it is clear that WOSUB is strictly one dimensional. This may be justified by considering the fact that the intention of creating WOSUB, has been analysing the flow characteristics in encapsuled PWR bundles as well as BWR bundle geometry 1 7) in which, based on a channelwise node, the flow is predominately one dimensional. Nevertheless, the transverse effects are not totally forgotten. In fact a natural turbulence exchange mechanism is considered. Furthermore, vapor diffusion accounts for the tendency of diffusion vapor in the higher velocity regions. -39- These effects are considered in the momentum equation which is given by (1 ~ ap az 9LP (-I)P) az )e+ az)ac [h 7 ): +P aGz t)fr at P + + ( az d )t [20] The last two terms stand for the axial momentum transferred into subchannel i and the turbulent shear stress, respectively. It is also evident that these two terms which connect the subchannel to its neighboring subchannels stem only from flow and pressure fluctuation and not transverse pressure difference as was discussed in Section 2.2.1.1. 4) Mixture energy equation which contains the inflow of enthalpy from adjacent subchannels. Generally, the dynamic slip model, as it is used in RETRAN, has advantages over both slip ratio correlations, as used in most versions of COBRA, and DFMI, as used in WOSUB, as follows (1 8 ) 1) The slip correlations are based on steady- state data whereas the application is for transients. 2) They highly rely on empiricism which may eliminate many mechanistic effects. -40- 3) The slip velocity VsL = V - VL can only assume positive values, hence the possibility of rising liquid and falling vapor cannot be predicted. 2.2.3 The Two-Fluid Model The inability of the simplified methods to treat the multidimensional, non-equilibrium separated and dispersed flows necessitates a better modeling of the two-phase flow. Anticipated reactor transients and postulated accidents like LOCA specially require a more realistic treatment. Those cases in which a one-dimensional HEM is not acceptable are tabulated in Table 2-1. TABLE 2-1(32) CASES WHERE 1-D HEM IS NOT ACCEPTABLE Multidimensional Effects Non-Equilibrium Effects Phase Separation Downcomer region ECC injection Small breaks Break flow entrance Subcooled boiling Steam generator Plena Post-CHF transfer Horizontal pipe flow Steam separators ECC heat transfer Counter current flow Steam generators Low-quality blowdown PWR ECC bypass Reactor core Reflood quench front BWR CCFL -41- The most flexible approach in modeling these cases is through using a two-fluid, full non-equilibrium concept, which is the most sophisticated model employed so far in treating the two-phase flow. The derivation of the field equations in their general tensor form is quite involved. presented in Ref. 33. A detailed derivation is A short-hand representation for the two-fluid model is 2V2T or UVUT which stands for unequal phase velocities and temperatures - whereas lVlT or EVET is used for HEM. The unknowns and equations in this model are summarized in Table 2-2. -42- Table 2-2 Two Phase and Single Phase Comparison with Respect to the Flow Equations Unknowns # of Unkn. Single V 3 Phase P 1 " of Momentum 3 Flow T 1 " of Energy 1 p 1 Equation of State 1 Case Type of Euations # of Egu. Conservation of Mass 6 a(void fra.) V g 1 3 1 6 Liquid Balance Equ. Vapor 1 " " Liquid Mom. " 3 " 3 Two* V1 3 Phase P1 1 Flow Pg 1 P 1 Vapor T1 1 Equation of State T g 1 in Each Phase Vapor " Liquid Energy " 1 " 12 * Table 2-3 gives a more detailed description of various approaches to modeling the two phase flow. 1 2 12 -43- The two-fluid concept is employed in the advanced thermal hydraulic codes such as COBRA-TF, from BNWL, TRAC, KACHINA*, SOLA-FLX, SOLA-DF from LASL and finally THERMIT, which is developed at MIT under EPRI sponsorship. TRAC is the state-of-the art primary loop analysis code. It employs a three-dimensional 2V2T model for the vessel and a one-dimensional drift flux model for the rest of the primary loop. The reactivity feedback is accounted for through coupling the point kinetic equations to the thermal hydraulic model. The same concept of volume and junction defined for RELAP series is used in TRAC as well. A cylindrical coordinate system is used in TRAC for modeling the three-dimensional reactor vessel. This doesn't satisfy the purpose of a common reactor core analysis with its square array pattern governed by the bundle design. THERMIT which is a vessel code, is basically the cartesian version of TRAC. Hence, the same field and constitutive equations used in TRAC is employed in THERMIT as well. A core or a fuel pin analysis in THERMIT essentially is based on treating a whole bundle cross-section as one node where the local details have been smeared throughout the cross-section. Therefore, neither TRAC nor THERMIT account for a turbulent mixing process. Devising a subchannelwise version for THERMIT using a coolant centered control volume * K-FIX and K-TIF are two versions of KACHINA developed at LASL. The first stands for fully Implicit Exchange numerics and the second for Two Incompressible Fields. -44- has been initiated at MIT. This version will include a turbulent mixing model. A summary of the aforementioned two-phase flow models is present in Table II-3. The notations and a description for specifications used in this Table are as follows: a) A partial non-equilibrium model, Tk Tsat, assumes one phase is at saturation, temperature of the other (k) phase computed. b) The notations T,q,r,M and E stand for viscous stress, conduction heat transfer, interphase mass, momentum and energy exchange respectively. c) The notations: Vr, VG - V, VG - J stand for relative velocity, diffusional velocity and the drift term respectively. A glance at this table shows clearly that although the 2V2T model imposes no restriction on the flow condition such as velocity or enthalpy, however it contains the largest number of constitutive equations and it seems that the empiricism which enters in these equations is introduced at a more basic level than the less complicated models such as 1V1T approach. -45- ro 0 IH It W U) m CN ,_m I Q) Hl a r- 0 0 Ero Q l r- a) 0 w~ -46- o 0 1b d or 0 .HQ - I (N 0' -I rd E- 0U) I o -47- 2.3 Classification According to Range of Application 2.3.1 LOCA CODES The major task of thermal-hydraulic LOCA codes is analysis of the severe cases that are encountered by the reactor core or the primary loop during the period of a loss of coolant accident. The four phases of a LOCA, for a PWR double ended cold leg break, in order of occurrence are as follows: 1)*- A blowdown phase which generally lasts for 30 seconds, with 2200 psia initial pressure, and ends when ECCS starts to work. 2)- About sixty seconds after break initiation ECC fills the lower plenum and reaches the bottom of the core (Refill Phase). 3)- The refill phase is followed by a REFLOOD phase which lasts for about 150 seconds, during which the core is fully flooded and quenched by the coolant. 4)- Long-term cooling then follows. It has become customary to call a code a LOCA code even if it is capable of describing only the first phase of the four aforementioned phases. At the same time two codes that are capable of handling the blowdown, may be entirely different with respect to their type. For example, one can be a component code whereas the other a loop code. To avoid any * This step is divided in two periods according to Ref. 31, namely: a) Adiabatic Liquid Depressurization, b) The Blowdown Period. -48- confusion, a classification is necessary with respect to the code's application and type, in addition to the physical models and numerical methods. A usual way to classify the LOCA codes is by categorizing them into two groups as follows. 2.3.1.1 Evaluation Model Codes The first group contains those codes which employ a conservative basis for their physical models. Such conservatism is mandated to satisfy the NRC acceptance criteria. Model. These codes are called the EM-codes for Evaluation They constitute the WREM package which has the capa- bility to analyse the postulated LOCA with ECC injection in accordance with current commission acceptance criteria (35) The codes which constitute the WREM package (36) are the existing computer programs which have been modified to comply with the USNRC criteria. Most of the RELAP series of computer codes are a LOCA Licensing code such as RELAP4-MOD5 and RELAP5-MOD7. Whereas RELAP3B-MOD101 is essentially devised to analyse ATWS and RELAP4-MOD6 and RELAP5 are not based on conservative correlations. RELAP4-EM is the only version of RELAP which is specifically modified to comply with acceptance criteria. The present EM codes comprise an assembly of codes run sequentially. Each member of the sequence is a stand- -49- alone code developed for some special application. With this respect, RELAP4-EM in conjunction with RELAP4-FLOOD (3 and TOODEE2 (3 7 ) constitute the WREM package which perform the PWR LOCA analysis. and MOXY-EM (3 8 analysis. ) 6) Also a combination of RELAP4-EM constitute the WREM package for a BWR LOCA The respected procedure for the above mentioned analysis are presented in Fig. 8 and 9. 2.3.1.2 Best Estimate Codes Most of the LOCA codes lay in this category. basic physical models used in these codes The ely on the best estimate assessment rather than conservative correlations. Unlike the EM-codes, the BE-codes are mostly devised as one large system code consisting of various functions previously performed via the separated stand-alone codes. This guarantees the proper compatibility and continuity between the various calculational phases. As an example, the multi-purpose loop code TRAC can be used for the analysis of the whole phases of a LOCA namely, blowdown, refill and reflood. Unlike the EM-codes which mostly use a homo- geneous equilibrium model in conjunction with a lumped parameter approach for their analysis, the best-estimate codes are much more demanding and the most recent BE-codes employ the state-of-the-art physical models. Therefore they -50- U w Z 0 D 6 I UC~ E-4 ¢0 0 0 Q z I O O Ir E Q 0 0 II 0 E- O D =r Ep w 1Y U] t 'N 4 ,< I U jI 0 E-f E4 .1 r J h E- U O EU 1-4 H pOP h o 0 r= U 0 P4 H U) cn ., Ln H r- N ., 5 I I CNJ C CN 0 12 z :R F:14 w P4 pq 0 ccn i-1 z 0 0 -:: PQ 0 EH4 D E4 p: U ol - | - Er U W HH H W E-1 Z -51- 0 z H HOQ Z 0 H Oz P4 0 0 U H U _D O H A U eqV 14 H Z H U) *Hq U PI i W Ii Eq P 2 0 0 p 0 U H E-4 H d 0 o ml o 0 W >- ,t 5 H P z 3 P4 O r Ln EHH 0< ! P4U C: Cl) W H U] W5: WZ E-4F: H UE H ~D ~H K ] a $4 z z :: o -52- can be used to evaluate the degree of conservatism employed in the licensing (EM) calculations. No codes have specifically been devised to handle the ATWS type of transients. In Table II-4 the causes and consequences of ATWS transients in both PWR's and BWR's are shown. 3. Two-Phase Heat Transfer Model The energy balance of the field equations contains the contribution of the so-called wall heat transfer which accounts for the amount of heat transferred into or out of the control volume through a combination of convection and conduction heat transfer. This requires models for the wall heat transfer. During hypothetical LOCA's, nearly all the two-phase heat transfer regimes are experienced by the coolant in the core of the NSSS, the steam generators, and the pipe of the hydraulic loop (see Fig. 11). This interdependence of the hydrodynamics and the wall heat transfer, as shown in Fig. 12, is accounted for in the thermal hydraulic LOCA codes through using a two-phase heat transfer package. These regimes are elaborated on in a pool boiling curve drawn for a fixed pressure, and shown in Fig. 13. According to Fig. 13, the path ABCDEF is obtained in a temperature -531I a) N a) 5-4 ·H 0 -HO 4)3 4- O) ot -P O to k 0) O )4-0 f H 0) -H - 4P U0) 0 rd a) ) 4 a) a) -H -rq O 0 c UH -*H 0 ,-4I 0) ,-4ol 3 q) -EC 0) to F. 04Id *H 0) U k 4 U) m >(D 4 04-Hr-H IN U) U) Imd O r8 I 3a (a 3 d H f (a 0) 0-H O , O0 k c)-H k I rI 0 OO rO CCf O',-I P4 H 4P 4 kU -P Pm"O 0 E4 3U O O O u ,--I i ,-H 4:3 r4 " 4·-4 r ·d ( d O -HO) r I r, a ,E·,:1 " ):) : rd *H . IN a) i, ,H 4U)-H . i q0 . H) . O O r u - (a 0 -, /L (D ~U) Id 0) OPO g UoH In l.. - .* . O O :O -I I -I - 0 EI ,-I 4 4J4J L R 0 Q) ( Wl k 4 4k k P4 5 3.-) -H UI I 0 . O . O ~ (D O OO OU Oa· mu -) aoU) 4 43 o to CN O E a040H H, 04H1'4 - -4H U; I I O-. k O O,') O U ) 4 E$ rH (0 O w 3u 0> t0 k 04 44 P i- ,-- Itzr En 0)1 (l n ,Ln I 1 a U)tr 0 "- (d E-4 I :>t 44 04J to 43 O ., )O > IH-H -d k4J k) a) O ". Z O Z H U H O-k ¢ -H C UH O #.. H P m4 , a O0 r > 4.U) Q k O si k k P3 O a) 0 w.0)4it U (d 4 u) 0) 5r - U) Id k 3 ) O Q)k a UO ) · . 4- -H k rId a4w 4-H k ) O0IW4Z U 0c4 U4 O aH - 0 r4 > mUU -H 5.-ri d >i ) E O k H)0 O ) i )- --I I > *H )O H4) O Z 0404-U Ok En <. 3 r I C 4u -H Id 3 tm 4--H 4. * 10-PO 4 . c))* > k .,q 0 U 4 HH k )0 ) ) *H 04 nl- >I 0::J > I t 0'~~~ . _ _ U) 0) -i 4 U) 0 U 0l) t7 0- >54 $ I H 04 M U) Q) QH __ 01) *i 0) U U) 54 4 rQ P 04-H 0 u --- H 0 - I 1) -H 404rd .----- 0) O U) a) 0 0U5- O 'a Ena) r. Jm ,-r4 O ro OH U O >) > co 0 CV ( 0 )0 0) -- -C Od O 1 L O U4 r-i I i44 3uz . _ >1 .- 4 k4 4 4J 4J -H . O 0) I -1 (D (0 (00) a) E! 040) 5-4 *H ~ 5-1 C) 5-40 :: 0) a 4. 4.3"O -i,- 0 -H 0 ri t 1f *1 U r 4 En tow 0 m En)-H0() X a: u),B 0-· 0 .- 4 ~-14-4 i-l (d 104 k 0I -H I I t .H to 4 0)00O k - -H -P 4. 04 O 3 0 ri >0 :3 : -iH HO 04-I 4 D 2 -54- I i - i i i i 'DISPERSED FLOW' FILM BOILING '2; SUBCOOLED TRA N5! i 'O11 , 04 BO1LING S W. '* -/0oFORCED ~ '~ Cl) o CL -0 SAT . RAT ED TRAN S!T!ON BO LING 2 / · oo S CONVECTIVE H.T. ro SUPER- Io SlEAlTE 0\ VAPOt cr C, L) z7- SUBCOOLED t3l kn NUCLEATE SAT! BO Li NG 0 CLEATE r.TS . CONVECti TO iLlOUtD LH.T. * *e \' _ SUBCOOLED LIQUID i.'O-,t? /~tx7~~C FORFC II ·A ED .fAS __$ATU ..L. 7 , D-'l .O'-TR Pt' S, A ___ -- _-- _,__I __I i r- CR.ASI,G QUALITTY ............... - I IXTURE EA SUPEH VUPEOR' Iq E - a . APPROXIMATE ROD TRAJECTORY I,, SLO\'¥.DO'w'N TRANSIENT Fig. 11 Heat transfer regimes traversed in blowdown. -55- I( r TreritJ f-r~ I':;I'L1 lle-jo Tron!.', !r ;.. ,; I?.g ir-,.:, - II __. ... . ...... f ;i l I e - P. sE f I to Stein I .Oiivetcti ll to 'icj St:.-rn upt:rhoJ C =· r ··: ·· -· ·- -,r r· L.i qid Deficient ¢oivection end fRadiation to Soturoted Liquid 'I, ?ist, Sprcty or Disper s ed .1 ond I [ . 35 i`. -· - Critical t4sut Flux 4( dryout) 2·r·. Annular or Film I Forced Convection I' Throulh li.lultd Flm ·5` ·· -. r. I-: i ·- · , 'r r·t ·· Slug ,Ch::rr, or Froth I ·- ., I I k23 cc). r. · f 0' I .D I, Nuc!!;te '. nvoct;on I tio.n por Firn 3ih P 9) :·. B3ubbl Froth tr I uot Flux P,, '0 ·: . . * I 9 ,I '." I n oii.) Suhcooled l-Jiling I ·· .'ater \',atr r - L-- n to (:nnve clion to Sinjl..-P P 0 Ji A Floa FI'. ._ Fig. 12 [19] Flow and heat transfer regimes in rod array with vertical upflow. -56- tT VC O a I-,. C', - Ca (A to,0Ew a) as I--F. .Ein Q e C to, n I- ar X XIL X Eo ._1 j LL I t C C lf (SH Z,I./nljE) ,,b 'Xnl .V3H 33Vians -57- controlled surface as the temperature is increased. In general, the same path on cooldown process is not followed in the heat-up process by the boiling mechanics. For example in a heat flux controlled surface, the path ABCC'F will be followed in which point c' indicates the new equilibrium state of the surface at the heat flux value qCHF At steady-state operation conditions, a NSSS fuel rod is a heat flux-controlled surface with a non-uniform axial heat flux distribution. In this case a reduction of the heat flux may be traced on the curve of Fig. 13 by the path EC'D E'BA( 1 8 ). Unlike the steady-state conditions, during a hypothetical LOCA it is not lear which mechanism prevails, since the fuel rods of NSSS's may behave as heat flux-controlled surfaces for some parts of the transient and as temperature controlled surfaces for other parts of the transient. 3.1 Heat Transfer Regimes and Correlations The recent thermal hydraulic LOCA codes have increased their capability of the two-phase heat transfer assessment by inclusion of more distinct heat transfer regimes and using more realistic correlations for calculations of the heat transfer coefficient in each regime. -58- As was discussed in section 3.1, a stand-alone LOCA code is capable of handling only one phase of hypothetical LOCA's, whereas integrated LOCA codes such as RELAP4-MOD7 and TRAC have the capability of calculating both blowdown and refill/reflood phases of a LOCA. This capability is made possible through inclusion of the unique features of bottom flooding (in PWR) and top spray (in BWR) of reflood heat transfer, in the blowdown heat transfer package. Such features are quench front, rewetting and liquid entrainment. Also thermal radiation and dispersed flow film boiling are specially pronounced in reflood heat transfer and are treated explicitly in the reflood heat transfer packages*. The heat transfer package which was used in the early versions of the RELAP series such as RELAP2, is used extensively in the thermal hydraulic codes**. This package is used in various versions of RELAP3 as well as RELAP3BMOD101. Later it was modified by replacing the quality by void fraction to determine the pre-CHF heat transfer regimes and by treating the transition boiling explicitly in which case the heat transfer coefficient is calculated using the MC DOUNOUCH, MILICH and KING correlation. Also * See REFLUX (3 9 ) package which is developed at MIT to analyse the reflood phase of a LOCA. **This package is essentially adopted from the THETA hotchannel code. -59- the Berenson and Groeneveld correlations were added to the formerly existing Dougall-Rohsenow correlation in the film boiling regime. RELAP4-MOD5 and RETRAN use this modified version and a simplified form of this new version was implemented in COBRA-IV-I. RELAP4-EM employs the new version with further modifications to satisfy the acceptance criteria. For example return to nucleate boiling is precluded once CHF happens. Also the GE correlation is added to the CHF correlations as an option to replace the Barnett correlation for BWR analysis. There are however several disadvantages associated with this package( 4 0 ) . 1) There is no CHF scheme to consider CHF during flow reversal or stagnation, which are charactieristic of blowdown in large, cold leg breaks in PWR's. 2) Use of Thom's correlation up to a void fraction equal to 0.8 which corresponds to a quality equal to 0.42 at 2250 Psia, which is above the quality range for which this correlation was verified. 3) Extensive use of correlations whose data base rely on tube or annular geometry, while their application is for rod bundle geometry. -60- 4) Using the correlations which have a steady- state data base for transient conditions. For this reason a heat transfer package called BEEST developed at MIT to overcome these drawbacks. BEEST (4 0 ) stands for BEst ESTimate heat transfer analysis. It is based on best estimate assessment rather than conservative correlations. Several tests of BEEST showed that it is able to construct the complete boiling curve where different heat transfer regimes are smoothly connected (Fig. 14). The heat gransfer selection logic in this package is based on the comparison of the clad surface temperature with the two distinct temperatures on the boiling curve, namely the temperature at the minimum stable film boiling point, TMSFB, and the temperature at the critical heat flux point, TCHF (Fig. 13). This is certainly an unambiguous, efficient and valid criterion for selecting the appropriate heat transfer regime. Once the regime is identified, the second step is to apply a chosen correlation for the heat transfer regime selected. The upflow and downflow heat transfer are treated separately through using the void fraction. The transition boiling in this package is treated in a unique way. This treatment is based upon an interpolation between the Q"MSFB and Q" CHF (which are the heat flux corresponding to the TMSFB and TCHF) with -61- +4. 4- - ILL CO de i1 v iC)' IC0' 1 TwA ~,,~. (°~) boiling curve - calculated 40] Effect lty of on qua Fig. 14 [40] Effect of quality on calculated boiling curve for low flow downflow, predicted by BEEST package. -62- respect to the temperature ratio as follows*: Q" TB = E Q"cF+ (l-E) QM CHF MSFB [21] where E (Twall = TMSFB)/(TCHF TMSFB) In equation (22), E may be interpreted as the fraction of wall area that is wet. BEEST uses the Biasi correlation for the CHF calculations. The Biasi correlation is essentially a dry-out correlation. Therefore it is appropriate for high flows and qualities where the vapor is a dominant factor leading to dry-out. For low flows and qualities the void-CHF correlation developed at MIT is used. The RELAP heat transfer package which was discussed earlier uses the Barnett correlation as well as * This concept was first introduced by W. Kirchner (see Ref. 41), in the form of a Log-Log interpolation: TB CHFTwall wall Log Q"C ( sat T ) where CHF - Log QMSFB Log TCHF - Log TMSFB Kirchner then applied his model in the heat transfer package of TRAC. -63- the modified Barnett and the B&W-2 correlations for the CHF calculations. In the pre-CHF regimes, the Chen correlation is used in the subcooled nucleate boiling, saturated nucleate boiling as well as forced convection vaporization. This correlation has predicted the existing data with reasonable agreement ( 3 8 ) to the other correlations such as: Schrock-Grossman, Bennett et al, Guerrieri - Talty. as compared Dengler-Addams, Sani and finally The Chen correlation is applicable to flow regimes from slug flow through annular flow. While its data base is for low pressures ( 4 2 ) , in most applications it is used at elevated pressures. Also, its dependence on the wall temperature which necessitates an alternative procedures, makes it less desirable. The advantages of the BEEST heat transfer package namely, treating the upflow and downflow separately, using a once through heat transfer regime selection logic, using wall temperature as a heat transfer regime selection tool, using a best estimate assessment and incorporating the new improvements in heat transfer, has made it acceptable to the state-of-the-art LOCA codes. THERMIT uses BEEST with some modifications such as replacement of the void fraction calculated from DFM by that calculated in THERMIT. TRAC heat transfer packages is also very similar to BEEST. In fact it can be considered as an improved version of BEEST -64- with the following additions: 1) Adding two options to the CHF correlation namely, the Bowring and the Zuber Pool boiling correlations. 2) Inclusion of the thermal radiation contribution in the film boiling regime. 3) Using a horizontal film condensation to represent the low flow rates. 4) Inclusion of a vertical film condensation regime. 5) Considering laminar and turbulent flow correlations in steady-state calculation for forced convection to two-phase mixture. A comparison of heat transfer selection logic and correlations used in different thermal-hydraulic codes is present in Tables 1.1, 1.2 and 1.3. The notations and specifications used in these tables are as follows: 1) Thermodynamic quality is represented by h-hf x =h where h represents enthalpy, whereas fg X represents the true quality X = is the mass flow rate. W g Wf+W Wfw ' where W r -65- 2) The notations , Tw, Tf, Tat P,G represent the void fraction, wall, fluid and saturated temperatures, pressure and finally mass flux, respectively. The dimension of the pressure and mass flux are in terms of "Psia" and "lbm/hr-ft 2 respectively. 3) The terms "High" and "Low" flow used in these tables are in accordance with the flooding correlation which read 1 1 J .1 f 2 + mJg = K g [23] where for turbulent flow m is equal to unity and Jf and Jg are dimensionless velocities: * f = Jf 1 2 f f Jg = J pg 1 2 [ gD(pf-P g)] f1 =J [gD [24] (pf-pg)] where D is pipe diameter and K is the flow criteria. For example, in low flow region according to Ref. 40, this criterion is * 1 Jf 2 * 1 - Jg 2 <1.36 for upflow 1[25] Jf 2 + Jg 7 <3.5 f 2 + Jg for downflow -66- 4) The letter x and a in the parentheses in Table 3.2, in front of the Thom and Schrock-Grossman correlations imply that an interpolation is made with respect to the quality and void fraction respectively, i.e. quality or void fraction weighted heat transfer coefficient. 4. Fuel Rod Model Temperature excursions of the fuel rod in case of any transient or accident are a major point of concern in the reactor safety analysis. A high temperature rise following severe transients is a threat to the cladding material whose integrity must be guaranteed in order to prevent any release of radioactive materials. There are four barriers preventing the release of radioactive fission gases to the environment under normal operating conditions namely, the U0 2 fuel, the fuel rod cladding, the reactor primary systems, and finally, the reactor containment building (40). Accordingly fuel melting, threatens the first barrier, and clad rupture violates the second barrier. The ECCS final acceptance criteria requires that failure of these barriers must be avoided under any circumstances. This necessitates a realistic fuel rod modeling, specially for the LOCA codes. In general a fuel rod model consists of an approach to solution of the general three dimensional, time dependent Ln H I O I I (d O Hs II Z 4J I r0 o a) W H 04 0 44 .-I O 0 I ___ __ - m C I4 H14 AI x I Ln acl d., Al X ,, A aI O AC: A\ L: 4I - -_I Al _L_ - I UL) 0 3:o 0H ,- o I. $ I U .0 0 0" rl *H t 01) A 3 0 r0 A 0 .( r--H:: 0- 4-) M rZ4 11 mI r-I H 0 0 P1 I LO X V P4 (N vI t0 *t r* U > s-4 r 4 Cd 0 .-I- -H 0 0 0-H C t A -' N0U4-)P4- --- U) ___ - - *rlI 0) )I ,H a) 5-4 J43 O H a Cd 0 ruxU 4-) > N E' Cd a) 44 v 4-J o o Al A a) P4I- -H I E I I W* O V V vl v H g tC Cd A 0) U) O | r-i S- 4 4-} Cd I Wt A V o v v a) O I) O - o O rl Hl-I 4-) rrd r a 03 MI0) a) Q4 0 4JrLr4 0Z rl __ IIpI i i I d ) r- -> ~4 In t~ I I L - (a Eq v 3 E - - tl)-i I i | 0o U U a I 0 _ I t) ,-I LO o r -.I ":: rZ4 C ::·r 0 H I tr O I H C 4) O a O ErC I r- - 0 0 rZ U X vI ^ Al r-- U > W. 44 x ^I m5 - 3 *H H *.q ] -I P' H -H- 0 I II o ) ____ ___ 0 .~ ^l __ 0 0 I ) o __ I v oE 0o E PL4 X 4d v o V 4 Ln v u I v >1 v o ( - I o J I 1--I rl- > r Q) 0 O Gr O,.QOZ m ; o - z E Al vI 0 Ep 0 U4J v Cd En o Ev A v A a- 4- U] E o -3 o rZ4 H E- I a) U > r- aO s-IO O O rT 4 4 E· I rl I] m/ :d H P 4 o vl vH p re) m o14 P0, P4 M 4 44a0 / ;·PI· u z EH v x. CN C) C ~1 z Ei * 0) I , H o0 o o o o v 04 E-4 V Itr I Z -iL ,¢ H x O H Ix r O > cn W o 2= P H rd ( v v H-4 EH l E4 v 1t 4- Cd Il; _ -68- .-1 44P I k u E: I > 4J a,) kc(d I C O: N N -, *Q rQ m9 ~ ~ - O m I * .- a z = S 4 -H 0 Z-H 1a)H a m U 0a) Cbm 0 0-H H 0 -H O a > I 4J a C 04J 4.J N ro O-H a) a) U >0 o0 0 IJ fd *o< a) 5k 4 N a) Er rd) Ord ) () 4 >Q) -rl 3 r: ar-I rl 4- H > 0 -H C d H a5r: I Cd I rd i 00 am o O0 U) o I10 O Cd X4 1 0 U a E5-4 U U m o () O O a) d O x .- a) Z o O 0 k 5 E3 0O zk a) C 0 0 4 E0 p a 0 O O0 a) a) a) a) rd U) a, :. O O 03 a d a) E zC m E-3 E- a) k I I)a) C 11. u: 0 u4- a G} r 4 4- 0 .,l Q I 0) 5- I 914 o a) 0 il E- 4U ro (-) a). a) O)H I a) 024-) rZ ) 4a 4-) rd rX "O UCU 0 a -H .) .- .,l u) I-----"- LO H I 0)0 0 E 4 P4 P4 ~4 ~-q c4 w 04 w z; z 10 a4 0 N' E W m o C4 m U) 0 o3 U -69uL .. -r i i rd I ro i Ii aO O .. O H I I r4J O d U> Uo a O ----------_--- t - _ -I I 0O rI rd. - ------- E) a) I -_-- -- ._--- - I U U 4) W U) 0) U) 0 OQ) 4-4 U) E-q Cd x Cd -tf + 4U) - O Cm r"- CI: 03 :O0) 0 00 rI 4J I -H U) 0io - -H, I 0 O Ok 4J (PtE- c- uz I O- 0 t rn 0 -H I 3 ..3 o O I *r I o C -H PC e PN e ~ P4 0 I H c, H> - w z-HS-l rd O u0 tN" n' t11 01 0 a H 0s-I m 0:. O U Or O O) 4 OH0, d H 0 -H -H) 01 0 0 - ,- Z7 a) C4 m I 0 k) mO U) -1 0.,,w I.. a 4.) .H0 3H O W) O Cd O MdH g 0 0 H 0 rH ro rd O a) 0 U) . /'1 W I- t1) Iva I ,H o /1 W . I· UU) U Um I Cd k pq H> -H . O H *H - N O r U> i . 0 W1 (aU) .. '' HI oU 0 4-I "U 0) U) 4 H-LO * m z I U) 'O 0 i 3O 03), E o *k ~3 *H U) c4 P4 U) 3 I -l U1 0 I k4 SH d E-q W --- U "G o H 0 O- O, 01) d ·- -- tri. a, q. O ro a 4J m P: > H H -_ -- .H m O C . 00 --- --- ---------- 4) ^--LI----------- : UU) 0 r,4U U) q > I U) l C-;__. I 0 En U001 O 4- Z 0 -H > 4r-,1 Ud 0 U) :d O- r. -W O Q Z , N ,H i -H I H -H PL . Ic z O,0 9>4 O t0 *( -4 CD ___ E- H C) 9 -70- Poisson's equation (heat conduction equation) p c - = [26] ' (K VT) + Q"' where T,t,p and cp represent the temperature, time, In material density and specific heat respectively. this equation K represents the conductivity tensor, K = Kij, i=1,2,3, j=1,2,3 and Q" is the heat source density which represents the amount of heat generated in the material per unit volume per unit time. Generally as the cylindrical shape of fuel rod dictates, a cylindrical coordinate system is chosen to expand the first term in the RHS of Equation 26. 4.1 Fuel Region The expanded form of equation 24 in the fuel region is: T 1 a3T p r 1 aT 3 a r ( k aT 3T) + Q where k in equation 25 is no longer a tensor but a time dependent scaler. This simplification is made possible through the valid assumption of homogeneous, [27] -71- isotropic solid for U0 2 and fuel rod cladding material. The first term in the right hand side of Equation 27 is considered in the more general form in COBRA-IV-1 as follows: r - direction: 1 (ra- k DT By assuming a=2, the cylindrical and a=l, the planar fuel can be treated. The total derivative in the left hand side of Equation 26 is changed to a partial derivative in Equation 27. This simplification is possible as long as a stationary solid is treated. This in turn is a valid assumption since the fuel centerline melting is to be prohibited by design under any circumstances. The azimuthal, or O-direction, conduction is ignored in all the reviewed fuel rod models. This implies an assumption of infinite circumferential heat conduction. The axial conduction, Z-direction, is only considered in COBRA- IV-1 and is ignored in the other codes. Further simplification to Equation 27 is possible by assuming that all physical properties are temperature independent, in addition to the isotropic assumption. is done for example in WOSUB. However, the temperature This -72- dependence of thermal conductivity and heat capacity is considered in TRAC and THERMIT. The latter uses a chebyshev polynomial fitted to the MATPRO(4 3 expressions which represent fits to experimental data for fuel and clad material properties. For example a cubic and a quadratic polynomial is used to fit the temperature dependence of p, cp and k of the fuel, respectively. The Kirchoff's transofrmation is used in COBRA-IV-1 to reduce Equation 27 to a linear partial differential By using this method the temperature dependence equation. of k is taken into account. As for the RELAP series, RELAP2, RELAP3, RELAP3B and RELAP4 use a simplified lumped model for their heat In these codes heat generation conduction calculation. is determined by reactor kinetics routines or by input specified values for power versus time. The fuel rod model used in these codes is patterned after the model used in the HEAT1 code. The final form of the heat conduction equation is presented in Table 4.1, equation 28. In this equation, the average temperature is defined by _ (11). _ TdV pc TVPd /vPC p dV n n =l(PV)n(Cp)n n n E n=l(pV) n (cp) n -73- where V and N represent the fuel volume and fuel pin annulus number respectively. the Also s in Equation 27 is fuel pin surface area. Equation 27 representes the U-tube as well as once- through steam generator modeled in RELAP3B-MOD101. Naturally Equation 27 does not include any heat source density term. Equation 28 reflects the cartesian geometry used Thermal conductivity and in the RETRAN fuel rod model. heat source density temporal and spatial dependency are accounted for. In Table .1, Equation 31, Equation 32 and Equation 33 represent the COBRA IV-1, COBRA-III P, WOSUB, TRAC and THERMIT, one dimensional heat conduction equation. Thermal source density Q"' in COBRA-IIIC & COBRA-III P is calculated as follows: Q = where D, AZ and Q" The total power is D AZ Q" [35] are the fuel rod diameter, the axial interval and the heat flux respectively. Now dividing by the fuel volume gives Q' Z = Q" 2 I TDf 4- Q" 4D -2 Df AZ [36] -74- where Df represents the fuel pellet diameter. The thermal source density in WOSUB is assumed to be spatially uniform but time dependent, whereas in TRAC and THERMIT it depends both on position and time. demonstrates the COBRA-IV-I fuel rod model. Equation 34 This equation is found by using the Kirchoff's transformation. 0 where k k To k(T)dt [37] is the conductivity at reference temperature To . Differentiating Equation 37 with respect to r and t and substituting in Equation 27 we will come up with Equation 34. 4.2 Fuel-Clad Gap The fuel-clad gap heat transfer coefficient is implicitly treated in those models in which clad and the fuel-clad gap, are lumped together. This is done for example in the COBRA-IIIC fuel rod model. However, upon the importance of the fuel-clad gap resistance to the heat flow, it is treated explicitly in RELAP3B-MOD101, WOSUB, TRAC and THERMIT. The gap heat transfer coefficient depends upon the fission gas product in the gap, the radiation heat -75- co tJ 1O (N C14 ,-t O IN (N m td r. rH W4 z o r-l u 4-) ~Lm cc- - ro 0r +"--' ii %, 4-) 1 Cd - k + + 4XX E-4 aIr Er, I Co E- c: 4J td 4 U) -, ' I CH + i H 4 t0,C4- En x -I v Co 0 0S ( N H a) II 0" E- 4-) -1 II 11 E- 0 o 4-) rl ) : Co H, 10 C)4-) U" a C1 rl. 0 0 U4J 0O ! 0d X I-E ,0 0 aO 4 rj1 o C UC 00-M 1) E _a) UH .rda@a) I4-O4-i 4-). o : 4-) O Cd H) C1) f Ii Z 04 4. ·-I 1.0 U 4-) H . r- )(o QI U) -H F m XO d Vsa Pc r/ z~r r'll ct a a C) ·-WI *,1 a) rC a HU1d 4 H ) a)a 0 ro U *H U) a) 90. Hr ZC )Q rdO a I a I -H 4 nl 01- .- Iu o 1* I E H H, HI* Ua) u) I rd V I a I.-1 U UCC= I.,U) HH * . Or I OH 0H 00 rdi 1-1 . E-E- a) a) a) - 4-i C) 0 I I C: -. 4-) r 0 d p 0 a X) r= H 1: 1-Q)Y Ordo O -H a) I-i a) a W d 3w) oo nl 0 II D4 O V dt *H Cd o rd r rbrlo) 0 ro I HO 1C3:> -76- transfer across the gap as well as the fuel-clad contact and fuel-clad pressing( 20 . These are modeled in the GAPCON, MATPRO, FRAP-S and FRAP-T codes. MATPRO model which is implemented in THERMIT. It's In this model all the above mentioned factors are considered. The model used in TRAC ignores the effect of fuel pressing against the clad, whereas it is correlated in THERMIT in terms of the fuel contact pressure against the clad. An effective gap heat transfer coefficient is used in WOSUB. Although this is not as realistic as the models used in TRAC and THERMIT, it still allows nodalization in the clad. COBRA-IIIC and COBRA-IV-1 assume the outer fuel surface and the inner clad surface are in a single node. In this case the conduction equation is written between the fuel and clad exterior surface. The heat transfer coefficient used for this purpose is defined as(4) 1 H 11 H g Yc k c where Yc and kc are the cladding thickness and conductivity and H 4.3 is the fuel-clad gap conductance. Clad Region As mentioned in the previous section 4.2, only WOSUB, TRAC and THERMIT permit clad nodalization. While -77- the clad region is assumed to be heat source free in WOSUB, a metal-water reaction is considered as a heat source in that region in the TRAC fuel rod model as This reaction happens at elevated well as RELAP4. temperatures, below the cladding material melting point, between zirconium and steam and is expressed as : Zr + 2H20 ZrO2 + 2H 2 + Heat [38] Both TRAC and RELAP4 use the parabolic rate low of Baker and Just to represent the rate of this reaction but in a different system of units+ + The mathematical statement of the parabolic rate low reads: dr dt ( a_ b R Jr)exp(. where r, R o, t and T represent the radius at each moment, the initial clad exterior radius, time and temperature respectively. constant values. In this equation a and b are By integrating Equation 39 between the initial and final radii of a time step, the mass of zirconium reacted per unit length during the time step will be found, The amount of heat generated in the clad is then proportional to this reacted mass, and it will be considered as the internal heat source in the clad region. + This exothermic reaction, results in hydrogen gas which poses a threat to the fuel rods in case of accidents by excluding the upper part of the rods to be covered by the coolant. ++ TRAC and THERMIT are the only thermal hydraulic codes using the SI units. -78- 5. Numerical Methods The mathematical models which were discussed in the previous sections are solved numerically in the computer codes, because analytical solutions are impractical. In the subchannel codes, the axial length of the reactor core is divided into several intervals which make each interval the computational control volume. The set of field equations in a finite difference form in conjunction with the constitutive equations are solved for the central volumes. The boundary conditions at the inlet of the core are, uniform or nonuniform pressure and coolant densities and enthalpies. The axial and radial heat flux profile must be specified. The solution is based on reaching a uniform pressure at the core outlet. purpose, a marching technique may be used. For this In a step by step, or marching technique, the calculation starts from the bottom of the core for all subchannels and moves upward. Inlet velocities are first assumed to be known and then solved alternatively through the external iteration loop. At each axial node, the cross-flow is guessed which allows solving the energy equation. A new value for cross-flow is calculated from pressure drop in each subchannel, which in turn is calculated from the momentum equation. -79- This internal iteration on cross-flow is continued until an acceptable pressure balance is reached. By knowing the heat addition into the axial cell and calculated cross-flow for the axial cell, the values of coolant density, velocity, enthalpy, and pressure can be determined at the exit of each computational cell. In turn, these values will be used as the information for the next axial cell. This procedure continues until the top of the core is reached where the criterion of uniform exit pressure is checked. If this criterion is not met, the external iteration loop which covers the whole channel length must be continued, using improved guesses of the flow division among the subchannels at the inlet. COBRA-II and HAMBO. This procedure was employed in The number of external iterations over the core length depends upon the coupling between subchannels. are small), If this is weak (e.g. if the cross-flows a single pass marching solution technique is adequate ( 2 2 ) , otherwise a multipass marching solution is necessary. This concept is used in COBRA-IIIC, in which a pattern of subchannel boundary pressure differentials for all mesh points is guessed simultaneously and then the corresponding pattern of cross-flow is completed using a marching technique up the channel. By updating the pressure differentials during each external -80- iteration loop, the effects of downstream will be propagated upstream. The procedure used in COBRA-IIIP is somewhat different. A new treatment is introduced for the transverse momentum equation which couples the adjacent computational cells. This includes the spatially semi- implicit treatment of the pressure field.+ Using this method guarantees the diagonal dominance of the matrices governing the pressure fields(49) The computed pressure field is then used in the transverse momentum equation to determine the cross-flow distribution. Applying the new concept in COBRA-IIIP has made it capable of increasing the number of computational subchannels markedly, i.e. from 15, in COBRA-IIIC, to 625 in COBRA-IIIP ( 4 4 ). Furthermore, it has increased the compu- tational effectiveness resulting in a shorter running time. ++ By introducing a scalar, value 0, having an arbitrary value between 0 and 1 the pressure field is written: [P] = [Pj] + (1-8) [Pj_l] By introducing this concept into transverse momentum equation, allows the cross-flow distribution to be driven by any combination of the pressure fields that exist at the top and the bottom of each plane of computa'tional cells. Unlike COBRA-IIIC, a double precision is used in computation of pressure field and gradients which is specially pronounced in cross-flow distribution calculation in the vicinity of grids. This in turn will increase the computer running time. -81- While the marching technique determines the flow condition under steady-state situations, in transients the whole procedure will be repeated for each time increment, implicitly in both COBRA-IIIC and COBRA-IIIP. The marching technique is also employed in WOSUB, where for the sake of numerical stability, a backward finite difference form is used in space and The lack of transverse equation and cross-flow time. is compensated by the concept of recirculation loop which is based on the assumption that the net volumetric flow recirculation around closed loops connecting communicating subchannels is zero (17) Simultaneous solutions of the finite difference form of the field equations written for the previously defined computational cell, is used in some subchannel codes another solution technique such as SABRE and COBRA-IV-1. Since the calculated values will be advanced in each time step explicitly, this additional option in COBRA-IV-1 be called the "explicit solution scheme". The following possibilities are available in COBRA-IV-1 solution algorithm: 1) Steady-state and transient calculation using the COBRA-IIIC implicit solution scheme. -82- 2) Implicit steady-state and explicit transient with either a AP or inlet flow boundary condition. 3) Explicit transient calculation based specified initial values with either a AP or inlet flow boundary condition and a zero flow initially. The addition of the explicit numerical scheme with a AP boundary condition makes COBRA-IV-1 capable of handling flow reversal, recirculation and coolant expulsions as well as severe flow blockage. These additional capabilities stem from solving a true boundary value problem rather then dealing with an initial value problem in the marching type solution technique. The additional numerical scheme in conjunction with the boiling curve package and the improved fuel pin-model makes COBRA-IV-1 capable of assessing accidents such as a LOCA, where due to the heat transfer package inability of analysing reflood, the code capability limits to the blowdown phase of a LOCA. It should be realized that due to numerical instabilities and convergence difficulties which mostly result from discontinuities introduced by the physical models, the "explicit scheme" of COBRA-IV-1 uses the strict HEM, i.e., it does not contain the Levy subcooled boiling model or any slip correlation. more the lack of computational effectiveness, as Further- -83- compared to COBRA-IIIC and IIIP, and using the reference pressure concept* which excludes any effect of compressibility and also large heat flux oscillations observed in prediction of depressurization transients ( 4 6 desirable. ) have made COBRA-IV-1 less Overcoming these deficiencies has been presumably the motivation of creating COBRA-DF and TF. As for the loop codes, a fully implicit solution scheme, temporally, is employed in all the RELAP series as well as RETRAN. An automatic time step variation is built in RELAP3B-MOD101. Using this feature, the time-step size increases automatically during slowly varying portions of a transient case of a computer run and vise versa. Both implicit and explicit solution schemes are employed in RETRAN. Unlike the implicit method which is unconditionally stable, the explicit method is conditionally stable in which the so-called "courant criterion" must be respected. This criterion reads: U where U AT AX < I or AT < AX [40] is fluid velocity, AT and AX are time step and The concept of reference pressure which ignores the sound wave propagation effects is employed in all the HEM versions of COBRA as well as WOSUB. This limitation is circumvented in COBRA-DF by using the ACE method, also see Ref. [45] in which this method is applied to the COBRA-IV-1 field equations. -84- mesh spacing respectively. It is clear that for fast transients which involve rapidly changing flow, very small time steps may be needed to resolve the flow evaluation. An explicit numerical scheme may be used for this purpose. Longer time steps are desirable in calculating the mild transients, which necessitates using the implicit numerical scheme. This scheme has not been used in the three dimensional thermal hydraulic codes, because, the fully implicit difference equations are very difficult to solve in more than one space dimension. A marching solution method may be applied to circumvent this difficulty, but as it was discussed earlier, no true boundary value problem can be handled by this method, only initial value problems in which general boundary conditions and local flow reversal cannot be treated. A compromise between the above mentioned techniques has been made in THERMIT by using a "semi-implicit" numerical scheme. As the name implies, both implicit and explicit schemes are employed in such a way that by differencing terms involving sonic propagation implicitly, limitations (U±C)AT I <1 have been eliminated,' whereas the liquid and vapor convection are treated explicitly. -85- Therefore the limitation imposed on the time increment to satisfy courant criteria (Equation 40) still exists. A default value is usually built in the codes which use a temporal explicit scheme to exclude the computational instability. 6. 6.1 Summary and Conclusions Summary A summary of the aforementioned methods and models used in each reviewed code and the range of each code application are presented in Table 6-1 through 6-5. Some of the terms which are used in these tables are further explained as follows: Small breaks (Table 6-1): postulated breaks that are smaller than about 10% of the double-ended break in the discharge flow area. Licensing codes (Table 6-1): notations are in accordance with the notations defined in Table 6-3. Homologous model (Table 6-4): described the centrifugal pumps and specifies relations, connections head, torque, flow rate, and rotational speed. -86- Table 6-1 Component and Loop Codes Comparisons Code Name CHARACTERISTICS APPLICATION TYPE -I O H m X M-i tl 4J X)X 0W W -i Q U) . O 43 1 XI W U Hr d l) W k W 0 U 4 U) Ax i al Forward Lateral Marching Steady State a Axial Forward Lateral Marching Steady State 2 1VlT Pseudo Axial Axial Pseudo Lateral Boundary Condition Steady Steady State 3 1V1T Axial teral Implicit Lateral lVlT COBRA I XX X X X COBRA II X X X X X X 1V1T COBRA III X XX X X X X XX X X X COBRA III-C U COBRA III-P X X X X X X lVlT COBRA IV-I X XX X XX X lVSlT COBRA-DF X X X X XX lVD1T COBRA-TF XX X XXX 2V2TSteady Steady Stat State & Transient SemiSteady Axial Implicit State Forcing Function Transient for cross flow Axial Ilicit Lateral Explicit Steady State & Transient Steady State & Transient State & Transient 4 4 5 6 7 8 -87- Table 6-1 Component and Loop Codes Comparisons (continued) l - TYPE APPLICATION w | l - CHARACTERISTICS - i l - w -IJ cU Code Name U) t7 z rL U ( ci) Q) PQ X ~4 o 0 C "o C -4 -14 0 O E-O) o a) m ¢) U u X O.ci) C WOSUB -P Hrd O q~ 0 rn C) M9 Uaa X X fl J)~ o 5: X U ci (D >O ,1a) 0a) 0 lVD1T (D 4ci) a, Qi U) 03 rJ (d o) z U Axial Lateral Imnplicit Forward Marching ci) U ATWS 17 (1) THERMIT X X X RELAP2 X X X X X X X X X X lV1T i, X X RELAP3 X X I RELAP3B (MOD101) X X X X RELAP4 X X X X X X X X X XX xX X X RELAP4-EM (2) RELAP4-FLOOD X Semi- 3-D x,y,z 2V2T 20 Implicit Lumped Implicit Parameters 9 I I I Lumpe d I Inplicit 10 Lumped i . . Parameters IITplzcit Lumped llct Parameters 11 Lumped Prameters Implicit 12 1V1T |lVlT X XI X I I1V1T | lVlT I 1V1T Lumped I Parameters Lumped lVD1T Lumped Implicit Parameters I 13 LImplicit 14 15 c.t 36 X RELAP4 (MOD6) X X RELAP4 (MOD5) X t ix I I - X x X - fx Txl X X , l : 1VD1T Parameters ) i I l l ii No diversion cross flow In conjunction with TODEE (for PWR) and MOXY (for BWR) . -88- Table 6-1 Component and Loop Codes Comparisons (continued) TYPE APPLICATION CHARACTERISTICS 4J Code Code Name ti~H . -S44~ z C) a) 3 o H :oW:z W oz "a0o za) Hu - a)W J 4 Qq 34 4 u W uM 1 o LH 3 .H - P4 mQ_ ir o ) m_ U) U M RELAP4 (MOD7) X XX XX RELAP 5 X XX RETRAN XX TODEE X MOXY TRAC-P1 1 X 5: u XX X XX X X s X X 4W J d r0 Zl H M rd J M r Q)4 o CD C >1. E k m o > C o ro 0o X lVDlT X X 2VTkTsat k sat X X1VDS1T(1 ) a 3 : X a o Lumped Parameters 1-D l-D Implicit mpli Implicit Implicit Explicit a tH P4 14 16 18 X X XX X 37 X X X X X 38 X X XX X X 2V2T 2 1VT T r, Dynamic slip model from the two-fluid theory 3-D 3-D , z SemiImplicit19 19 -89- Table 6-2 Models and Methods Used in Some Thermal Hydraulic Codes g I CAPABILITY OR MODEL Conservation Equation CODE NAME m CORA ICPmRA I I BNWL BWL 1967 1970 .. * *_. pRA COBRA COB.A COBRA WSUB THEPMIT IIIC IIIP IV-I III BNWL BNWL MIT ] BNWL MIT MIT 1971 1973 1977 1976 1978 1979 I 1 Homogeneous equilibrium model, HEM(lV1T) X X X X 2 One dimensional mass, momentum, energy equation X X X X 3 Drift flux model 4 Separate continuity equation for liquid and vapor phases 5 Separate momentum and energy equation for liquid and vapor phases 6 Turbulent liquid-liquid mixing in the subcooled region, in energy equation 7 Turbulent shear stress in mixture momentum equation X 1 I X I It X X X I, X X, X i Flow Solution-Steady State: I 9 Marching method marching) (forward X i Ii X X i 10 True boundary value method i 11 New treatment of transverse momentum equation i i I X I I I Pseudo boundary value method X x Numerical Scheme 8 I X i IX X I: I x X -90- Table 6-2 Models and Methods Used in Some Thermal Hydraulic Codes (continued) CAPABILITY!~~~~~~~~ OR MODEL CODE NAMEt~~~~~~ CODE NAME CAPABILITY OR MODEL !1 Flow Solution- Transient: COBRA COBRA COBRA OBRA '0BRA COBRA I II III IIIC IIIP IV-I WSUB i 12 Fully implicit 13 Semi-implicit 14 Explicit, arbitrary flow field and boundary conditions (ACE method) THE i i X X X X X X Flow Energy Solution: 15 Spatially explicit 16 Spatially implicit X X X X X X Equation of State 17 Reference pressure 18 Local Pressure 19 Superheated steam properties 20 Steam table that contains the derivative of fluid properties X X X X X X; X X X . X X Transverse Transport Cross Flow Model: 21 Pressure resistance only 22 Transient momentum equation X X X 23 Forced diversion cross flow X X X -91- Table 6-2 Model and Methods Used in Some Thermal Hydraulic Codes (continued) CAPABILITY OR MODEL CODE NAME - COBRA COBRA MRRA lCnOBRA COBRA WDSUB IIIC IIIP IV-I II III I 241 Lateral momentum flux X THERMIT X Two dimensional transverse flow s 2b . X II Turbulent Mixing: i I 261 Single phase turbulent mixing X X X X X 27 X X X X X Two phase turbulent mixing 28 Vapor drift on a volume to volume exchange basis X Accident Analysis 29 Severe flow blockage, coolant expulsion, flow reversal X 30 Recirculation loop X X X X X X X X X X Single Phase Flow 31 Nonuniform channel friction 32 Laminar and turbulent friction correction 33 Hot wall friction correction X X X X X X X X Two Phase Flow One-dimensional slip flow 34 __ _ _ -_ _ ~ ~ - X - I----" X .. .-`--" -92- Table 6-2 Model and Methods Used in Some Thermal Hydraulic Codes (continued) CODE NAME CAPABILITY OR MODEL . _ _ _ . _ _ _ _ _ _ COBRA 'COBRA COBRA COBPA lrnn'DA I II III IIIC IIIP _ " 35 I ^ fflnD . I "' U IV-I X One-dimensional drift flux model (Zuber-Findlay) t 36 Three-dimensional (x,y,z) nonhomogeneous nonthermodynamic equilibrium flow 37 Subcooled voids Model) 38 Vapor generation rate term in subcooled boiling to account for thermodynamic nonequilibrium THERMIT (Levy X X X I X x I X X x ! X X X X X- X Heat Transfer i ! 39 CHF correlation A Boiling curve package q u X X I Heat Transfer Regime Selection Tool: 41 Void fraction and CHF 42 Quality and enthalpy 43 Local clad surface temp. i i II X X ! i i iI X X i x X~~~~~~~~~~~~~~~~~~~~ Heat Conduction-Fluid: I 44 Radial conduction 45 Axial conduction I X X X X X x X X x Fuel Rod Model 46 Specified axial & radial heat flux X x X . . x .. -93- Table 6-2 Models and Methods Used in Some Thermal Hydraulic Codes (continued) CAPABILITY OR MODEL CODE NAME COBRA COBRA COBPA COBPA COBRA COBPA I II III IIIC IIIP IV-I X X 47 One-dimensional heat conduction equ. (r-direc.) 48 Two-dimensional heat conduction equ. (r,z-direc.) 49 Implicit finite difference solution scheme 50 Collocation method 51 Orthogonal collocation technique (MWR) X 52 Temperature dependent thermal conductivity X 53 Transient (time dependent) heat source density 54 Constant fuel-clad gap heat transfer coefficient 55 Thermal radiation and WOSUB THERMIT X X X X X X X X X X X X X X X interfacial contact in the gap heat transfer coefficient 56 Planar or cylindrical fuel X 57 Axial fuel zone X -94- Table 6-2 Models and Methods Used in Some Thermal Hydraulic Codes (continued) CAPABILITY OR MODEL I_ CODE NAME w -- .. RELAP 2 ; (If- 1 n -n V- 7n 'i n Pf-T I n+-i INEL 1968 'n equilibrium model, HEM (1VlT) .Homogeneous 2 1 One-dimensional mass, momentum, energy equation .- I RELAP 3 REIAP 3B RELAP 4 IINEL 1976 INEL 1970 X X X X i RTRAN TRAC-P1 INEL 1977 X X X X X X X X X X X X i 3 Inclusion of K.E. and P.E. in energy equation I i I 4 Consideration of area and density change in momentum equation 5 Dynamic slip model from the two fluid theory to account for nonhomogeneous: flow Three dimensional (r,9,z) flow for vessel 6 X x I II,1x x x One-dimensional flow with drift flux model for the rest of the primary loop 7 LASAL 1978 INEL 1973 Numerical Scheme 8 Fully implicit solution scheme temporally 9 Factor to modify the fullyi implicit scheme X X X X x X X X X iIX 10; Automatic time step variation 11 X Explicit scheme I ' .. . I , _ -95- Table 6-2 Models and Methods Used in Some Thermal Hydraulic Codes (continued) i _ CODE NAME CAPABILITY OR MODEL i Equation of State 12 Local Pressure 13 Steam table that contains the derivative of fluid properties 14 Extension of the steam table above the critical pressure RELAP 2 RELAP 3 RELAP 3B RELAP 4 X X X I EKt'XAN X X T1'AU- I X X X Physical Model Pump Characteristics: X X X X X X X X X X Consideration of inertial effect X X X 21 Consideration of frictional torque X X X 22 Consideration of bearing and windage torque 23 Option for two phase pump 24 Motor torque option, pump stop option, & dimensionless head ratio difference data in two phase pump 15 Homologous pump 16 Only one pump coastdown curve 17 Independent tripping on the independent signals 18 Pump is at a junction 19 Pump is in a volume 20 . X X X X X X X . . X -96- Table 6-2 Models and Methods Used in Some Thermal Hydraulic Codes (continued) I I CODE NAME CAPABILITY OR MODEL i I I RELAP 2 IRELAP 3 RELAP 3B RELAP 4 RETRAN 1 rI n I r I Choked Flow: 1) C . J M- 57 luJuy WI a L.w.-jrllaz: t i cAh V .n .I.I.U , V 111. v V JX v flow model i zb . 27 __ 1I.1- __-.- _ _ -1 - tienry-FtausKe ana extenae Henry-Fauske I _ X X X Sonic choking I Heat Exchanger: I X 28 Non-conduction model (input specified secondary temperature & a constant effective heat transfer X i X X X X i I coefficient) 29 i i Time dependent heat : X X exchanger (Input specified j table of normalized power versus time) I Time dependent secondary temperature X 31i U-tube steam generator (one dimension heat conduction X 30 I equation) 32 X Once through steam generator jSingle Phase Friction !Factor: 33 Laminar friction factor X X X 34 Turbulent friction factor X X X 35: Input specified friction Ir iUUI; ,i -- -- X X X ~i ~ ~~ [.... j~~~~~i -- _ -97- Table 6-2 Model and Methods Used in Some Thermal Hydraulic Codes (continued) i i CODE NAME CAPABILITY OR MODEL - REJAP 2 RAP 3 RELAP 3B REAP 4 Two Phase Frictional Multipliers: - -- / - RETRAN TRAC-P1 - X X 36 Modified Baroczy correlation X 37 New correlation based on modified Baroczy X 38 Beattie correlation using Bennet flow regime map 39 CISE model X 40 Annular flow model X 41 Chisholm model X 42 Homogeneous correlation 43 Armand model X X X X Bubble Rise Model: 44 Linear approximation for the density of bubbles versus height X X X X X X X X X Heat Transfer Package: 45 Correlations for pre and post CHF 46 Ability to construct boiling curve X X X 47 Treatment of transition boiling explicitly X X X 48 Condensation calculation 49 Reflood heat transfer package . X X X X . II -98- Table 6-2 Model and Methods Used in Some Thermal Hydraulic Codes (continued) CODE NAME CAPABILITY OR MODEL -- - - --- - ; Heat Transfer Regime T-- - - RELAP 2 REAP 3 RELAP 3B RELAP 4 RETRAN TRAC-P1 Selection Tool: 50' Quality and CHF 51 X X X Void fraction and CHF X X ~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ,,.,~~~~~~~~~~~~~~~~~~~~~~~ b5z Local caa surface temp. Fuel Rod Model: 53 Lumped approach x X X X X J1- X 54, One-dimensional heat conduction equation 551 Variable gap size during a transient 56 X Thermal radiation and interfacial contact in the gap heat transfer coefficient 57 Fuel-clad gap heat transfer coefficient burn up dependent 58 Thermal conductivity temperature dependent X X X xX X 59! Exothermic metal-water reaction considered as a heat source in the cladding material X X 60 X X X X 61 Explicit numerical scheme i Implicit numerical scheme I X X X Jr -99- l I ko 4) Hd co -100- 6.2 Conclusions A study was made of a number of well-known thermal hydraulic codes. This study attempted to cover both models and methods used in these codes summarizing their basic elements in various tables. The following results are drawn: 6.2.1 Component Code The COBRA series as well as WOSUB and THERMIT fall in this category. As for the COBRA codes, the steady-state versions namely COBRA-I, II, III are certainly obsolete by now and are not suggested for further considerations. COBRA-IIIC and its MIT version COBRA-IIIC/MIT*, which utilize the HEM model of two-phase flow, can be used for both normal operation and transient conditions. The great flexibility of COBRA-IIIC/MIT to simulate core regions, bundles and subchannels at the same time, makes it more desirable to use. An alternative marching solution using implicit numerical scheme, which converges on the cross-flow is used in this code. The fuel rod model in conjunction with Pre-CHF heat transfer correlations make it a fast running code for steady-state as well as mild transients. On the other hand COBRA-IIIP/MIT is capable of handling a larger number of computational subchannels with considerable computational effectiveness, since it deals with the diagonally *The major difference between these two codes is that the latter uses a dynamic data management subroutine which allows the dimensions of the principal arrays as well as the total computer storage requirement to be a function of the problem size. dominant pressure matrix. However it should be realized that these codes, COBRA-IIIC/MIT and IIIP/MIT, are not devised for the purpose of accident analysis such as a LOCA. In fact neither the physical models nor the numerical methods have such capability. For example, lack of a heat transfer package and use of the marching solution technique do not allow any extreme flow as well as reliable fuel rod temperature calculations. A step toward the analysis of severe transients and/or accidents is taken in COBRA-IV-I in which more realistic physical models, with respect to the 2-D fuel rod model and heat transfer package, are implemented. More importantly, a field equation solution technique, explicit solution, is employed in addition to the COBRA-IIIC implicit type solution scheme. It is the simultaneous set of differential equations using explicit solution technique which makes COBRA-IV-I capable of handling severe flow blockages, flow reversal, coolant expulsions and other extreme flow situations. Also, the field equations solution using ACE technique allows specifying a AP boundary condition which relaxes the impractical specification of inlet flow boundary condition in severe transients such as blowdown. -102- Despite these advantages, it should be recalled that several disadvantages are associated with this code which are as follows: 1) It relies extensively on the HEM, 2) Computational ineffectiveness, 3) The time change of local pressure is ignored, 4) Subcooled void is excluded in the explicit scheme. A more realistic two-phase flow model that relaxes the assumption of the HEM is used in the newer versions of COBRA namely, COBRA-DF and TF, still under development. Except for the highlights of the models used in these codes that are presented in the summary tables, there is little additional information available on these codes for the time being. Since analysis of the BWR normal operation and transient conditions is more demanding with respect to the two-phase flow modeling, the recent efforts in the BWR models have been focused on using more realistic assumptions with regard to vapor-liquid momentum exchange, or phase distribution as for example the WOSUB code. -103- Unfortunately, the marching type solution technique, simplifying assumptions such as the reference pressure concept and lack of a complete boiling curve package, in the current WOSUB version, put a limit on its application. For example, no reliable blowdown calculations can be performed with WOSUB in its present form. Furthermore, WOSUB is not supposed to be used for very fast transients. It is now clear that extreme flow situations which are a point of concern in severe transients can only be evaluated by using more physically accurate field equations. The two-fluid concept provides the potential for increased accuracy in modeling the two-phase flow. By implementing this concept in the most recent codes, such as THERMIT, a number of limitations imposed on the flow characteristics are relaxed. Now the motion of two-phase in different directions, having different temperatures, velocity and pressure can be realistically analysed in three dimensions. The best estimate heat transfer package, BEEST, and an improved fuel rod model, specially with respect to the material temperature dependent and fuel-clad gap modeling, which are included in THERMIT, provide a reliable fuel rod temperature calculation and DNBR prediction. A semi-implicit numerical -104- method is used to circumvent the instabilities associated with the explicit numerical scheme. In spite of the above mentioned advantages, it should be remembered (see section 1.3) that the need for mathematical models representing physical phenomena increases with the degree of sophistication of the twophase flow modeling. Furthermore, the difficulty with the general two-fluid approach is that the exchange processes coupling the phases are currently not thoroughly understood [ As a result, despite the possible shortcoming of the HEM, it is not evident that a homogeneous equilibrium model is incapable of predicting adequately some parameters of interest such as vapor flow rate and fuel-clad temperature. Furthermore, in some flow regimes HEM gives surprisingly good results. Therefore, as a final conclusion, for normal operation and mild transients, COBRA-IIIC/MIT and COBRA-IIIP/MIT are still the best available tools. Several shortcomings of these codes such as the fuel-rod model and lack of a heat transfer package and the like may be overcomed by implementing the state-of-the-art models used in the sophisticated codes such as THERMIT. Severe transients and accident analysis are certainly advised to be aalysed by THERMIT. Upon the completion of a subchannelwise version with coolant centered control volume, THERMIT will -105- be the first three dimensional non-homogeneous fully thermodynamic non-equilibrium subchannel code for the nuclear reactor core thermal hydraulic analysis. 6.2.2 Loop Codes The RELAP series as well as RETRAN and TRAC fall in this category. The old versions of RELAP such as RELAP2, and most of the RELAP3 versions are obsolete and need not be further considered. RELAP3B-MOD101 is the only updated version which uses several options for heat exchanger and steam-generator modeling, not even used in TRAC. This version is specially devised for the ATWS analysis. It uses a combination of old models of RELAP3 such as heat transfer package, and new models introduced in RELAP4, such as homologous pump, in addition to several unique features such as variable gap size during a transient and heat exchanger modeling. In light of the detailed information about RELAP4/MOD6 and MOD7 and RELAP5, the available highlights of loop modeling are presented in the related tables. From these tables it is clear that the major step toward the non-homogeneous non-equilibrium modeling of the primary loop is taken in the model making process of RELAP5. Such effort is also done in RETRAN through introducing -106- the unique feature, DSM, "dynamic slip model". RETRAN is specially devised to analyse the reflood phase of a LOCA. The state-of-the-art of the loop codes however is TRAC which is capable of handling all the phases of a hypothetical LOCA. As a result the following codes are suggested for the transient loop calculations. First, RELAP3B-101 which is a one-dimensional, HEM code and it may be used for ATWS transients. Second, RETRAN which uses l-D, DSM and improved physical models, third, TRAC which uses l-D, DFM for the loop calculations, and realistic physical models. -107- REFERENCES 1. D.S. Rowe, "Cross-Flow Mixing Between Parallel Flow Channels During Boiling, Part 1: COBRA: Computer Program for Constant Boiling in Rod Arrays", BNWL371 Pt. 1, Pacific Northwest Laboratory (1967). 2. D.S. Rowe, "COBRA II: A Digital Computer Program for Thermal Hydraulic Subchannel Analysis of Rod Bundle Nuclear Fuel Elements". BNWL-1229, Pacific Northwest Laboratory (1970). 3. D.S. Rowe, "COBRA III: A Digital Computer Program for Steady-State and Transient Thermal Hydraulic Analysis of Rod Bundle Nuclear Fuel Elements", BNWL-B-82, Pacific Northwest Laboratory (1972). 4. D.S. Rowe, "COBRA IIIC: A Digital Computer Program for Steady-State and Transient Thermal Analysis of Rod Cundle Nuclear Fuel Elements", BNWL-1695, Pacific Northwest Laboratory (1973). 5. R.E. Masterson, "Improved Multidimensional Numerical Methods for the Steady-State and Transient Thermal Hydraulic Analysis of Fuel Pin Bundles and Nuclear Reactor Cores", Ph.D. Thesis, MIT (1977). 6. 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, Pacific Northwest Laboratory (1976). 7. C.A. McMonagle, et al., "COBRA Code Development", Fifth Water Reactor Safety Research Information Meeting (1977). 8. M.J. Thurgood, et al., "COBRA-TF Development", Sixth Water Reactor Safety Research Information Meeting (1977). 9. K.V. Moore and W.H. Rettig, "RELAP2 - A Digital Program for Reactor Blowdown and Power Excursion Analysis", 1D0-17262 (1968). -108- 10. K.V. Moore and W.H. Rettig, et al., "RELAP3 A Computer Program for Reactor Blowdown Analysis", IN-1321 (1970). 11. - , "User's Manual for RELAP3B-MOD101, A Reactor System Transient Code", NUREG-0004 (1976). 12. K.V. Moore and W.H. Rettig, "RELAP4 - A Computer For Transient Thermal-Hydraulic Analysis", ANCR-1127 (1973). 13. C.E. Slater, et al., "Small Break Analysis Models in RELAP-4 Computer Code", ANS Trans. Volume 22 (1975). 14. K.R. Katsma, "RELAP4 (MOD6 and MOD7), LWR Reflood Code Development and Verification, Fifth Water Reactor Safety Research Information Meeting (1977). 15. WREM: Water Reactor Evaluation Model, NUREG75/056 (Rev.l), Nuclear Regulatory Commission (1975). 16. V.H. Ransom, et al., "RELAP5 Code Development and Results, Fifth Water Reactor Safety Research Information Meeting" (1977). 17. L. Wolf, et al., "WOSUB, A Subchannel Code for Steady-State and Transient Thermal-Hydraulic Analysis of BWR Fuel Pin Bundles, Volume I", M.I.T. Energy Lab. Electric Power Program (1978). 18. K.V. Moore, et al., "RETRAN-A Program for OneDimensional Transient Thermal-Hydraulic Analysis of Complex Fluid Flow Systems", EPRI NP-408(1977). 19. An Advanced Best Estimate Computer - , "TRAC-Pl: Program for PWR LOCA Analysis", NUREG 1CR-0063 (1978). 20. A Computer W.H. Reed and H.B. Stewart, "THERMIT: Program for Three-Dimensional Thermal-Hydraulic Analysis of Light Water Reactor Cores", Unpublished Draft of an EPRI Report. -109- 21. W.T. Sha, "A Summary of Methods Used in RodBundle Thermal-Hydraulic Analysis", Presented to the Meeting Sponsored by the NATO Advanced Study Institute, Instanbul, Turkey (1978). 22. J. Weisman and R.W. Bowring, "Methods for Detailed Thermal and Hydraulic Analysis of Water-Cooled Reactors", Nuclear Science and Engineering 57, 255-276 (1975). 23. G.P. Gasspari, et al., "Same Consideration on Critical Heat Flux in Rod Clusters in Annular Dispersed Vertical Upward Two Phase Flow", 4th Intern. Conf. Heat Transfer, Paris (1970), also Energio Nucleare, 17, 643 (1977). 24. N.E. Todreas, "Analysis of Reactor Fuel Rod Assembles Isolated Cell, Distributed Parameter Analysis Appraoch", Class Notes at M.I.T. Unpublished. 25. N.E. Todreas, et al., "Single Phase Transport Within Bare Road Arrays at Laminar, Transition and Turbulent Flow Condition", Nuclear Engineering and Design, Volume 30 (1974). 26. S. Fabic, "Review of Existing Codes for Loss-ofCooland Accident Analysis", The Advances in the Nuclaer Science and Technology, Vol. 10, 1977. 27. R.T. Lahey, F.J. Moody,"The Thermal-Hydraulics of a Boiling Water Nuclear Reactor, ANS Monograph Series (1977). 28. F. Seifaee, "An Improved Numerical Method for the Steady-State and Transient Thermal-Hydraulic Analysis of Fast Breeder Reactor Fuel Pin Bundles", M.I.T. Nuclear Engineering Department, NE Thesis (1978). 29. C.W. Stewart, "COBRA-IV: The Model and the Method", BNWL-2214, NRC-4 (1977). 30. C.W. Solbrig and D.J. Barnum,"The RELAP4 Computer Code: Part 1. Application to Nuclear Power - Plant Analysis," Nuclear Safety. Vol. 17, No. 2 (1976). -110- 31. W. Wulff, et al., "Development of a Computer Code for Thermal Hydraulics of Reactors (THOR)", Sixth Quarterly Progress Report, BNL-NUREG-50569 (1978). 32. - , "Multi-Fluid Models for Transient Two-Phase Flow ", NP-618-SR (1978). 33. M. Ishii, "Thermo-Fluid Dynamics Theory of TwoPhase Flow, Eyro Ues (1975). 34. Y.Y. Hsu and R.W. Graham, "Transient Processes in Boiling and Two-Phase Systems", McGraw Hill (1978). 35. 10 CFR Part 50, "Acceptance Criteria for Emergency Core Cooling System for Light Water-Cooled Nuclear Power Plants", Federal Register, Vol. 39, No. 3 (1974). 36. "WREM: Water Reactor Evaluation Model, Revision 1", Division of Technical Review, Nuclear Regulatory Commission, NUREG-751056 (1975). 37. G.N. Lauben, "TOODEE2, A Two-Dimensional Time Dependent Fuel Element Thermal Analysis Program", NRC (1975). 38. D.R. Evans, "MOXY: A Digital Computer Code for Core Heat Transfer Analysis", IN-1392 (1970). 39. W. Kirchner and P. Griffith, "Reflood Heat Transfer in a Light Water Reactor", AICHE SYMPOSIUM SERIES No. 164, Vol. 73. 40. T.A. Bornard, "Blowdown Heat Transfer In a Pressurized Water Reactor", Nuclear Engineering Department, M.I.T., Ph.D. Thesis (1977). 41. W.L. Kirchner, "Reflood Heat Transfer In a Light Water Reactor", U.S. Nuclear Regulatory Commission report NUREG-0106 (1976). 42. J.G. Collier, "Convective Boiling and Condensation", McGraw Hill (1972). 43. "MATPRO-Version 09: A Handbook of Materials Properties for Use in the Analysis of Light Water Reactor Fuel Rod Behavior", Idaho National Engineering Laboratory report TREE-NUREG-1005 (1976). -111- 44. L. Wolf and R. Masterson,"An Efficient Multidimensional Numerical Method for the ThermalHydraulic Analysis of Nuclear Reactor Cores", Nuclear Engineering Department, M.I.T., Tucson, Arizona (1977). 45. L. Wolf, et al., "Review of Iterative Solution Methods for Pressure-Drop Forced BWR Core Analyses", Nuclear Engineering Department, M.I.T. (1977). 46. M. Massoud, "Comparison of Conservative and Best Estimate Heat Transfer Packages with COBRA-IV-I", Nuclear Engineering Department, N.E. Thesis (1978). 47. Z. Rouhani, Chapter 14 of "Two-Phase Flows and Heat Transfer with Application to Nuclear Reactor Design Problsms", edited by J.J. Ginoux, Hemisphere Publishing Corporation. 48. N.E. Todreas, class notes, Department of Nuclear Engineering, M.I.T. -112- APPENDIX 1 In the following, the assumptions which have been made in derivation of the continuity equation in the COBRA codes are discussed. The mass balance for liquid is given by: aat [(l-a)pl] I [(l-a)p (1 1 f) 1 ] = + V -r ( Mass balance for vapor follows at [(Pv r where ) ] + V'[(PvV)]= is the phase change rate. (2) Now using a pseudo single phase concept by assuming: p = (l-a)p (3) + pV = (l-a)p1 V (3) + Pvv and adding terms (1) and (2) we will come up with at )p l + (p)] + V[(1-)pV l + Pv v] = Now using the "averaged values" introduced in (3) in the equation (4), we will come up with the single phase mass balance at( at + V V) = : 0 ,0 Equation (5) in a more involved form becomes: (5) (5) -113- a y(pdxdydz) + at(pddxdz)dx + a (pVdxdz)dy + t DX DY " (6) (pVdxdy)dz = 0 where the element of volume, dn dy dz, represents the control volume used. dz PVxdydz+ .pVxdydz (pVxdydz)dx z dy Fig (1) Now in equation (6) by eliminating dz we'll have: ata(pF- (dxdy) + (Vdxdy) +dxdy) + (pVdydx)= (7) In equation (7) if we assume that the control volume shown in Fig. (1) has only an infinitesimal height, dz, and its cross section normal to the z-axis "finite" instead of "infinitesimal" Dt (i.e., dxdy A) and by substituting a at t((pVA) + 9z pVA = m ) -9a + ~ +t(pA A , we will come up with: + au (pVA) + DX (PVA) = 0 (8) we have: -(m) ay + ) + a(m)= (9) A "subchannel control volume" concept used in COBRA assumes a finite cross section normal to the z-axis which -114- has the same area, A, equal to the physical subchannel size, and an infinitesimal height, dz.t a- .p VA CL Idz t,_j -- , fuel rod t dz w . = 1] I OVAAidetpvA Lde face . 31 [wi ] = lbm/hr-ft Fig (2) In Fig. (3), w*ij represents the 'cross flow" concept. Mass balance written for control volume shown in Fig. (3) becomes: d d pVA + T-(pVA)dz + w* .dz - pVA = (pAdz) dz 13 dt d d dt(pA) + (pVA) + w* = 0 m = pVA (10) Comparing equation (10), used in HEM versions of COBRA, with the pseudo single phase continuity equation, equation (9), we conclude that the cross flow term w*ij represents the two dimensional form of mass flow rate per unit length, i.e., (ax + ay )m + w*.ij (11) This phenomena is only considered in COBRA-IV-I. None Flow Blockage: of the remaining HEM versions of COBRA account for this tThis infinitesimal length becomes the axial interval between the selected axial mesh points in the finite difference solution technique. -115- Therefore in these versions, i.e. COBRA III and feature. COBRA III-C and ..., equation (10) is simplified by assuming aA 0 at 0 + A dt dp p d(m) dz + aP. A Ai am. + t ax ++ w*.. = 0 N i - i where the subscript C w* (12) ij represents the subchannel under consideration. Rewriting again equation (11): d d -- p + -(pV)A dt dz Adt + NYW j=l = 0 (11) where in this equation, the averaged values are emphasized by using a (-) sign. This in fact follows using an integral form of mass equation which reads Jsp( at fTpdT + where T, · n)dA = 0 (12) s and v represent volume, surface and velocity respectively. Define the volume and surface averaged values for density and mass flux: p = <<p>> = 1 - f pdT (13) pV = <pV> = where V = TdT and A = If p(v fsdA . . )dA Also one bracket show area averaged and two bracket shows a volume averaged value. -116- Evaluating equation (12) for the control volume shown in Fig. 2 we have V P+ VAlupper face where w sie faces - VA Ilower face p(v n)dA and (w )bz = O (14) V = A-Az (15) If Az becomes small, in limit we will have equation (11). -117- APPENDIX 2 The details of the correlations introduced in Tables 3-1 to 3-3 are presented here. 1. Pre-CHF Heat Transfer Correlations The heat transfer correlations used in the pre-CHF heat transfer regimes are as follows: Dittus-Bolter(A -l) correlation: 0. ,Gh GDh 0.8 Cpp 0.4 * h = 0.023(-) k k Dh ....Data Base: GDh - Re = Pr - L Dh h greater than 10,000 0.7 to 100 greater than 50 Sieder-Tate ) correlation: GDh 0.8 h = 0.023(-) (Cpp 0.4 * 1. k Dh Data Base: As for Dittus-Boelter correlation. Thom correlation(A-3): P T = T + 0.072 e 1260 w sat Data Base: Vertical upflow of water 1 w . ( )0.14 -118- Round tube: 0.5-in. diameter, 60-in. length Annulus: 0.7-in. ID, 0.9-in. OD, 12-in. length Pressure: 750 to 2000 Psia. Mass flux: 0.77 x 106 to 2.80 x 106 lbm/hr-ft 2 Heat flux: to 0.5 x 106 Btu/hr-ft 2 Schrock-Grossman correlation(A-4): rPV(l-x)Dh1O.8 h = 0.023L C p1 0.4 (- -) 1 [2.5(X 0.75 ) k ]( ) where the inverse of the Martinelli-Lockhart-Nelson Parameter for turbulent flow is 1- x %0.9 Pls 0.5 P( 0.1 tt gs Data Base: Water in round tubes Diameter: 0.1162 to 0.4317 in. Length: 14 to 50 in. Pressure: 42 to 505 Psia Mass flux: 0.175 x 106 to 3.28 x 106 lbm/hr-ft2 Heat flux: 0.06 x 106 to 1.45 x 106 Btu/hr-ft 2 Exit quality: 0.05 to 0.57 Chen* correlation(-A Qw =hNB (Tw- ): 5 Tsat) + hc(Tw- Tf(Z)) where hc F(.023)k0.6G 08 (l-x)0.8c 0.4 -0.4 -0.2 C p 'P h *The values of F and S factors should be found from the corresponding graphs (see Ref. 42 of this report). However, the relations given for them here are the curve fitting which are derived by Butterworth. More details are presented in Ref. 40 of this report. -119- and k0.79 C 0.45p 0.49 0.25 h = 5(0.00122 p 1 C hNcB = S(0.05 0.29h 0.24 0.24 w1heefg Pg where F = L -1 0.736 2.35(Xtt + 0.213) Xt1 < 1.0 tt X > 0.1 where Xtt is the same parameter as introduced in the Schrock-Grossman correlation. The value of S is as follows: [1 + 0.12 (Rp)l.14]S = [1 + 0.42(Rp) 0.78 ]- 0.1 ; Rp < 32.5 ; 32.5 < R p ; Rp < 70 > 70 where RTp is the effective two-phase Reynolds number G(l-x)Dh 1.25 (10-4 Data Base: See Ref. 42 of this report. 2. CHF Correlations The critical heat flux correlations named in Table are as follows: Babcock and Wilcox (B&W-2)(A6) correlation: Q = 1.155 - 0.4 07 (Dh) 8 B A[(0.3 7 0 2 xl 1 0) 6(0.59137x (12.71) (3.054x106G) CHF - 0.15208 hgs G] G) -120- where A = 0.71186 + (2.0729 x 10 ) · (P - ]000) and B = 0.834 + (6.8479 x 10 -4 ) (P - 2000) Data Base: Vertical upflow of water in rod bundles Heated equivalent diameter of subchannels: Heated length: Pressure: Mass flux: 0.20 to 0.50 in. 72 in. 2000 to 2400 Psia 0.75 x 106 to 4.0 x 106 lbm/hr-ft 2 Thermodynamic quality: -0.03 to 0.20 Uniform axial flux distribution. Westinghouse (W-3) (A-7) correlation: QCHF = {[2.02 - 0.430(0.001P)] + [0.172 - 0.000]P exp[18.2x - 0.00413P-x]}(1.16 - 0.87x) [(0.148 - 1.6x + 0.173xlxl) (G/106) + 1.04] [0.266 + 0.836 exp(-3 .15Dh)] [0.826 + 0.0008(ht - h i )]10 Data Base: Diameter: 0.2 to 0.7 in. Length: 10 to 144 in. Pressure: 1000 to 2400 Psia Mass flux: 1 x 106 to 5 x 106 lbm/hr-ft 2 Quality: less than 0.15 -121- Barnett Correlation (A-8) F QCHF 10(6 A(hfg/6 4 9)C +- B(hls-hi) 10 = Z where A = 67. 4 5D0 B = 1.85 08 (G x 10 - 6 ) 0 D1261(Gx10- 1 9 2 {1- 0.74 4 exp [ 0- 0 .5 1 2 Dhy(Gxl 6) 00817 C = 185 *D1.415(GxlO-6)0.212 hy For Annuli the heated and wetted equivalent diameters, Dh and Dhy, are given by Dhy = (D s -D I) and Dh = (D s -D )/D I where D s is the diameter of the shroud and D I is the diameter of the inner rod. Data Base: Vertical upflow of water in annuli geometry Diameter of inner rod: 0.375 to 3.798 in. Diameter of shroud: 0.551 to 4.006 in. Heated length: 24.0 to 108.0 in. Mass flux (9 x 10-6): 0.140 to 6.20 lbm/hr-ft 2 Inlet subcooling: 0 to 412 Btu/lbm Inlet Pressure: 1000 Psia Uniform axial heat flux. )]} -122- Biasi correlation(A-9) The critical heat flux is given as the higher of the two values from the following equations. For the low quality region QCHF 1.883 x 103 f(p) Dn G1/6 hy G/ - 6 x] For the high quality region 3 1.78 x 103S (P) n 0.6 QCHF (l-x) (-x) hy G where n = 0.4 for Dhy > 1 cm n = 0.6 for Dhy < 1 cm f(p) = 0.7249 + 0.099P exp(-0.032P) s(p) = 1.155 + 0.149P exp(0.019P) + 0.889P 889P 10+P2 Data Base: Diameter: 0.3 to 3.75 cm Length: 20 to 600 cm Pressure: 2.7 to 140 bar Mass flux: 10 to 600 g/cm 2-s Quality: 1/(l+Pl/Pg) to 1 VOID-CHF correlation(A- ) -1 QCHF = (1-a)0.9T(24) Data Base: See Ref. A-10. 0.25 0.5 hfpg I[ggcG(Pl-Pg )] -123- 3. Post-CHF Heat Transfer Correlations The transition boiling, film boiling region and thermal radiation heat transfer coefficients introduced in Table 3-3 are as follows: McDonough-Milich and King correlation (A-11) a) As used in RELAP Q = QCHF - h(TwTw, CHF) where h is dependent on Pressure as follows: h 979.2 1180.8 1501.2 P 2000 1200 800 Data Base: Vertical upflow of water in round tubes Diameter: 0.152 in. Length: 12.5 in. Mass flux: 0.2 x 106 to 1.4 x 106 lbm/ft 2-hr Wall temperature: less than 1030°F 800, 1200 and 2000 Psia Pressure: b) As used in RETRAN For pressure greater than 1200 Psia, P-1200 h = hl 1 2 0 0 - (hl1 2 0 0 - hi 2 0 0 0 ) ( 800) and for P < 1200 Psia, h = l1200 + hl800 - hl200) ( 1 400- P -125- Bjornard-Griffith correlation (A12); QT.B. h Tw -T 1 where QT.B. = QCHF + (1 -E)QMSFB where £ = [(Tw TMSFB )/(TCHF - TMSFB ) - ] Data Base (see Ref. A-12). Dougall-Rohsenow correlation (A-13) h = DG (p (l-x)+x)3] 0.8 [P rg0. 0.4kg 0.023[(DG) () Pg P rg D The physical properties are evaluated at saturation conditions. If n < 0.0, the term [ (l1-x)+xJ is taken equal to 1.0 which causes the correlation reduces to the Dittus Boelter correlation. Data Base: See Ref. A-13 Groeneveld 5.7 correlation(A-14): h = 0.052 k GD [ Dh 9 p 0.688 1.06 -1.06 x + i (l-x)] P Y Pi rl and the modified Groenveld 5.7 correlations is given by k h = 0.052 where G'D -x 0.688 [ I Dh P gadfm 1.06 -1.06 P rw Y -125-4 Y = 1 - 0.1O(p 1lPj9 g1 - 1) ( 1-x) Data Base Diameter: 0.06 - 0.25 in. Pressure: 500 - 1400 Psia Mass flux: 6 x 106 Quality: 0.1 - 0.9 Heat flux: 1.4 x 105 - 7 x 105 Btu/hr-ft2 x 106 lbm/hr-ft 2 Radiation heat transfer coefficient 4 -T h = oF(T 1 4 )/(Tw-T ) 1 where 1 £ 1 + 1 a 1 c = emissivity of the wall, and a = absorptivity of the coolant jl -126- NOMENCLATURE FOR APPENDIX 2 C : specific heat Dh: heated equivalent diameter Dhy: wetted equivalent diameter h: heat transfer coefficient hfg: heat of vaporization k: thermal conductivity (evaluated at bulk temp. of coolant) L: length P: pressure QW surface heat flux Tsat: saturation temperature Tw: wall temperature x: quality Greek pi: density of the liquid phase Pg: density of the vapor phase Pl ,: density of saturated liquid p9 : density of saturated vapor gs a: surface tension P: viscosity viscosity of the liquid phase g19: viscosity of the vapor phase -127- Appendix 2 - References A-1: F.W. Dittus and L.M.K. Boelter, "Heat Transfer in Automobile Radiators of the Tubular Type," University of California Publications in Engineering, 2, 443-461, 1930. A-2: E.N. Seider and G.E. Tate, Ind. Eng. Chem., 28 (1936) 412a-35. A-3: J.R.S. Thom, W.M. Walker, T.A. Fallon and G.F.S. Reising, "Boiling in Subcooled Water During Flow Up Heated Tubes or Annuli," Proc. Instn. Mech. Engrs., 180 Pt 3c, 226-246, 1966. A-4: V.E. Schrock and L.M. Grossman, "Forced Convection Vaporization Project," TID 14632, 1959. A-5: J.C. Chen, "A correlation for boiling heat transfer to saturated fluids in convective flow." Paper presented to 6th National Heat Transfer Conference, Boston, 11-14 Aug. 1963. ASME Preprint 63-HT-34. A-6: J.S. Gellerstedt et al, "Two Phase Flow and Heat Transfer in Rod Bundles," ASME Conference, Los Angeles, California, Nov. 18, 1969. A-7: L.S. Tong, Boiling Crisis and CriticalIHeat Flux, TID-25887 (1972). A-8: P.G. Barnett, A Comparison of the Accuracy of Some Correlations for Burnout in Annuli and Rod Bundles, AEEW-R558 (1968). A-9: L. Biasi et al, "Studies on Burnout," Part 3, Energia Nuclear, 14, S (1967). A-10: P. Griffith, T.C. Avedisian and J.P. Walkush, "Countercurrent Flow Critical Heat Flux," Annual H.T. Conference, San Francisco, California, August 1975. A-ll: J.B. McDonough, W. Milich, E.C. King, Partial Film Boiling with Water at 2000 Psig in a Round Vertical Tube, MSA Research Corp., Technical Report 62 (1958) (NP-6976). A-12: Bjornard, T.A., "Blowdown Heat Transfer in A Pressurized Water Reactor," Ph.D. Thesis, M.I.T., August 1977. -128- A-13: R.S. Dougall and W.M. Rohsenow, Film Boiling on the Inside of Vertical Tubes with Upward Flow of the Fluid at Low Qualities, MIT-TR-9079-26 (1963). A-14: D.C. Groenveld, "Post-Dryout Heat Transfer at Reactor Operating Conditions," AECL-4513 (1973). -129- REPORTS IN REACTOR THERMAL HYDRAULICS RELATED TO THE MIT ENERGY LABORATORY ELECTRIC POWER PROGRAM A. Topical Reports A.1 A.2 A.3 A.4 A.1 (For availability check Energy Laboratory Headquarters, Room E19-439, MIT, Cambridge, Massachusetts, 02139) General Applications PWR Applications BWR Applications LMFBR Applications M. Massoud, "A Condensed Review of Nuclear Reactor Thermal-Hydraulic Computer Codes for Two-Phase Flow Analysis", MIT Energy Laboratory Report MIT-EL-79-018, April 1979. J.E. Kelly and M.S. Kazimi, "Development and Testing of the Three Dimensional, Two-Fluid Code THERMIT for LWR Core and Subchannel Applications", MIT Energy Laboratory Report, MIT-EL-79-046, December 1979. A.2 P. Moreno, C. Chiu, R. Bowring, E. Khan, J. Liu, N. Todreas, "Methods for Steady-State Thermal/Hydraulic Analysis of PWR Cores", Report MIT-EL-76-006, Rev. 1, July 1977 (Orig. 3/77). J.E. Kelly, J. Loomis, L. Wolf, "LWR Core Thermal-Hydraulic Analysis--Assessment and Comparison of the Range of Applicability of the Codes COBRA IIIC/MIT and COBRA IV-1", Report MIT-EL-78-026, September 1978. J. Liu, N. Todreas, "Transient Thermal Analysis of PWR's by a Single Pass Procedure Using a Simplified Nodal Layout", Report MIT-EL-77-008, Final, February 1979, (Draft, June 1977). J. Liu, N. Todreas, "The Comparison of Available Data on PWR Assembly Thermal Behavior with Analytic Predictions", Report MIT-EL-77-009, Final, February 1979, (Draft, June 1977). -130- A.3 L. Guillebaud, A. Levin, W. Boyd, A. Faya, L. Wolf, "WOSUB A Subchannel Code for Steady-State and Transient ThermalHydraulic Analysis of Boiling Water Reactor Fuel Bundles", Vol. II, User's Manual, MIT-EL-78-024, July 1977. L. Wolf, A. Faya, A. Levin, W. Boyd, L. Guillebaud, "WOSUB A Subchannel Code for Steady-State and Transient ThermalHydraulic Analysis of Boiling Water Reactor Fuel Pin Bundles", Vol. III, Assessment and Comparison, MIT-EL-78-025, October 1977. L. Wolf, A. Faya, A. Levin, L. Guillebaud, "WOSUB - A Subchannel Code for Steady-State Reactor Fuel Pin Bundles", Vol. I, Model Description, MIT-EL-78-023, September 1978. A. Faya, L. Wolf and N. Todreas, "Development of a Method for BWR Subchannel Analysis", MIT-EL 79-027, November 1979. A. Faya, L. Wolf and N. Todreas, "CANAL User's Manual", MIT-EL 79-028, November 1979. A.4 W.D. Hinkle, "Water Tests for Determining Post-Voiding Behavior in the LMFBR", MIT Energy Laboratory Report MIT-EL-76-005, June 1976. W.D. Hinkle, ed., "LMFBR Safety & Sodium Boiling - A State of the Art Report", Draft DOE Report, June 1978. M.R. Granziera, P. Griffith, W.D. Hinkle, M.S. Kazimi, A. Levin, M. Manahan, A. Schor, N. Todreas, G. Wilson, "Development of Computer Code for Multi-dimensional Analysis of Sodium Voiding in the LMFBR", Preliminary Draft Report, July 1979. -131- B. Papers B.1 B.2 B.3 B.4 General Applications PWR Applications BWR Applications LMFBR Applications B.1 J.E. Kelly and M.S. Kazimi, "Development of the Two-Fluid Multi-Dimensional Code THERMIT for LWR Analysis", accepted for presentation 19th National Heat Transfer Conference, Orlando, Florida, August 1980. B.2 P. Moreno, J. Liu, E. Khan and N. Todreas, "Steady-State Thermal Analysis of PWR's by a Simplified Method," American Nuclear Society Transactions, Vol. 26, 1977, p. 465. P. Moreno,J. Liu, E. Khan, N. Todreas, "Steady-State Thermal Analysis of PWR's by a Single Pass Procedure Using a Simplified Nodal Layout," Nuclear Engineering and Design, Vol. 47, 1978, pp. 35-48. C. Chiu, P. Moreno, R. Bowring, N. Todreas, "Enthalpy Transfer Between PWR Fuel Assemblies in Analysis by the Lumped Subchannel Model," Nuclear Engineering and Design, Vol. 53, 1979, 165-186. B. 3 L. Wolf and A. Faya, "A BWR Subchannel Code with Drift Flux and Vapor Diffusion Transport," American Nuclear Society Transactions, Vol. 28, 1978, p. 553. B.4 W.D. Hinkle, (MIT), P.M. Tschamper (GE), M.H. Fontana, (ORNL), R.E. Henry, (ANL), and A. Padilla, (HEDL), for U.S. Department of Energy, "LMFBR Safety & Sodium Boiling," paper presented at the ENS/ANS International Topical Meeting on Nuclear Reactor Safety, October 16-19, 1978, Brussels, Belgium. M.I. Autruffe, G.J. Wilson, B. Stewart and M.S. Kazimi, "A Proposed Momentum Exchange Coefficient for Two-Phase Modeling of Sodium Boiling", Proc. Int. Meeting Fast Reactor Safety Technology, Vol. 4, 2512-2521, Seatle, Washington, August, 1979. M.R. Granziera and M.S. Kazimi, "NATOF-2D: A Two Dimensional Two-Fluid Model for Sodium Flow Transient Analysis", Trans. ANS, 33, 515, November 1979.