THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING APPLICATION OF ADIABATIC APPROXIMATION TO THE SPACE-DEPENDENT KINETICS OF BOILING WATER REACTORS IIideaki Nishihara A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan Department of Nuclear Engineering 1965 May, 1965 IP-705

Doctoral Committee: Professor George L. West, Jr., Chairman Assistant Professor Ziya A. Akcasu Professor Frederick G. Hammitt Professor William Kerr Associate Professor M. Rasin Tek

ACKNOWLEDGEMENTS I wish to express my sincere gratitude to the members of my doctoral committee for their guidance in the course of the thesis research. In particular, I wish to acknowledge stimulating discussions with Professors G. L. West, Jr., Chairman of the doctoral committee, G. L. Gyorey and Z. Akcasu. The use of the IBM-7090 computer at the University of Michigan Computing Center is hereby acknowledged. I am grateful for the financial support to the present investigation given by the School of Engineering in the form of scholarships and by the Michigan Memorial-Phoenix Project in the form of a fellowship. I am also indebted to the United States Government, whose Fulbright scholarship made it possible for me to initiate studies in this country. The assistance of the Industry Program of the College of Engineering in printing the manuscript is appreciated. Finally I wish to thank my wife, Kyoko, for her encouragement and assistance during this work. ii

TABLE OF CONTENTS Page ACKNOWLEDGEMENTS................................................... ii LIST OF TABLES...................................................... v LIST OF FIGURES...................................... vi NOMENCLATURE...................................... viii CHAPTER I INTRODUCTION................................... 1 A. M otivation............................................ 1 B. Historical Background.......................... 1 C. Adiabatic Approximation............................... 3 D. Present Research........................... 8 II FEEDBACK DYNAMICS OF MULTIREGION BOILING WATER REACTORS... 11 A. General Remarks....................................... 11 B. Basic Equations of Two-Phase Flow..................... 12 C. Integration of Basic Equations....................... 13 D. Slip Ratio............................................ 21 E. Inlet Water Velocity.................................. 22 F. Boiling Boundary Dynamics............................. 28 G. Heat Transfer from Fuel to Coolant.................... 29 H. Pressure Dynamics.................................... 33 III FORMULATION OF MATHEMATICAL MODEL BY METHOD OF ADIABATIC APPROXIMATION............................................ 35 A. General Remarks....................................... 35 B. Step-I -- Steady State Diffusion Calculations......... 36 C. Step-II -- Correlation Process....................... 41 D. Step-III -- Computer Program of Mathematical Model.... 47 IV EXAMINATION OF ADEQUATE NUMBER OF REGIONS................. 57 A. General Remarks....................................... 57 B. Nuclear Characteristics............................. 58 C. Feedback Characteristics........................ 66 D. Conclusion.......................................... 69 iii

TABLE OF CONTENTS (CONT'D) Page V SPACE-TIME BEHAVIOR OF BOILING WATER REACTOR............ 70 A. General Remarks...................................... 70 B. Stepwise Change in Reactivity....................... 71 C. Sinusoidal Change in Reactivity...................... 75 D. Sinusoidal Change in Vertical Acceleration........... 78 E. Discussion of Results............................... 84 VI CONCLUDING REMARKS................................... 89 A. Summary and Conclusion..........................~.... 89 B. Future Work................................0*..... 92 APPENDICES..........................,............. 94 A. Derivation of Reactor Kinetic Equations from Transport Theory.................................... 95 B. Vessel Pressure Dynamics............................ 101 C. Sample Boiling Water Reactor......................... 108 D. Simplified Diffusion Code (SDC)...................... 124 E. Main Program Listing.......................,,.. 134 REFERENCES........................*..,.,..**....*~* * *. *.... 0* 155 iv

LIST OF TABLES Table Page I Regression Coefficients for Shape Function and Reactivity... 45 II Properties of the Sample Reactor........................... 70 III Comparison of Some Analytical Models for Boiling Water Reactor Kinetics Studies.................................... 88 IV Composition of the Sample Reactor Core...................... 109 V Cross Sections for Wigner-Wilkins Spectra................... 109 VI Effective Fast Group Cross Sections......................... 115 VII Cross Sections for Fast Fission Calculation................. 117 VIII Resonance Integral and Resonance Escape Probability......... 121 v

LIST OF FIGURES Figure Page 1 Division of the Core Along the Axis, and the Trapezoidal Approximation of the Steam Void Fraction Distribution........ 16 2 Pressure Difference in the Flow Loop......................... 25 3 Mean of Void Distribution Used for Steady State Calculations. Data Spread Is the Standard Deviation........................ 37 4 Mean of Flux Distributions Obtained by the Diffusion Code. Data Spread Is the Standard Deviation........................ 40 5 Initial Value Extrapolations................................ 51 6 Flow Diagram of the Program.............................. 54 7 Continuously Distributed and Sectionally Averaged Reactor Constants............................ 59 8 Flux Distribution Calculated by FOG.......................... 60 9 Flux Distribution of Multiregion Reactors................... 61 10 Eigenvalue vsO Number of Regions............................. 63 11 Flux Distribution Calculated by SDC......................... 64 12 Computer Time vso Number of Regionso Computation Terminated After 300 Iterations..................................... 65 13 Behavior of Multiregion Reactor with Space-Independent Nuclear Characteristics....................... 67 14 Space- and Time-Dependent Behavior of Void and Flux Distribution. At t=O, 90 Cents of Reactivity Is Added in Steps...... 72 15 Behavior of a Natural Circulation Boiling Water Reactor to a Step Increase in Reactivity of 90 Cents...................... 73 16 Frequency Response of Flux Level to External Reactivity. pext=10 sin 2,f t (in Cents).................................. 76 Pext=10 sin 2tft (in Cents).76 17 Frequency Response of Void Fraction to External Reactivity. Pext = 10 sin 2ntft (in Cents)................ 77 18 Responses of the Reactor to a Sinusoidal Change in Reactivity ext = 10 sin 2jt (in Cents)................................. 79 vi

LIST OF FIGURES (CONT'D) Figure Page 19 Responses of the Reactor to a Sinusoidal Change in Vertical Acceleration. g = go(l+0.5 sin tt)......................... 83 20 Comparison of Different Models in a Sinusoidally Changing Acceleration Field. g = g (1+0.5 sin t)...................85 21 Ratio of Amplitude Variations Obtained by One-Region Model to Space-Dependent Model for an External Perturbation of g = go(l+0.5 sin 2 ft)...................................... 86 22 Reactor Vessel Pressure Dynamics............................ 102 23 Slow Group Constantss........................... 112 24 Fast Group Constants........................................ 116 25 Four Factors and koo................................... 122 26 Nodal Diffusion Approximation.......................1...... 125 vii

NOMENCLATURE A Area a Constant B Regression coefficient B2 Transverse buckling b Thickness of fuel slab Ci(t) Space-integrated delayed neutron precursor density of i-th group ci(r,t) Delayed neutron precursor density of i-th group c Heat capacity D Diffusion coefficient De Hydraulic diameter Dh Constant E Neutron energy e Internal energy Fi Integrated friction force of i-th region F(t) Time-dependent quantity, defined in Appendix A i(u'- u) Scattering frequency f Thermal utilization f(x,t) Steam void fraction f. Steam void fraction at i-th node f Isothermal friction factor iso fnon-iso Non-isothermal friction factor f Frequency G Mass flow rate viii

Gi Integrated gravity force of i-th region g0 Acceleration due to gravity gdi Delayed neutron spectrum of i-th group gf Prompt fission neutron spectrum gt Averaged fission neutron spectrum H Heat capacity of coolant per unit core height h Specific enthalpy J Conversion factor J_ Neutron current in positive direction J Neutron current in negative direction K Constant k MMultiplication factor kg Thermal conductivity of fuel L2 Diffusion area M Mass N Total number of regions Ni Atomic density of i-th nuclei N(t) Time function of space-time dependent neutron flux density PO Collision probability p Pressure p Resonance escape probability (in Appendix C) Q Integrated power Q' Effective integrated power Q* Integrated power given to boiling portion of incipient boiling region q Power released in coolant ix

R Frictional resistance Re Reynolds number RH Hydraulic radius RI Resonance integral r Radius r Position vector S(r,_,u;t) Neutron source S Ratio of thickness of moderator slab to that of fuel slab S/M Surface-mass ratio of fuel rod s Laplace transform variable T Coolant temperature t Time u Lethargy V Volume v Velocity X Steam quality x Axial coordinate in reactor core a Volume fraction of holes a28 Capture to fission ratio of U238 92 ~.i Delayed neutron fraction of i-th group Y Slip ratio A Axial length of one core-region 6nb Fraction of nonboiling length in incipient boiling region E Fast fission factor x

TQ ~ Regeneration factor hFuel temperature A(t) Effective neutron lifetime x. Decay constant for i-th group delayed neutron precursors o0 Mean free path Dynamic viscosity o0 Average of cosine of scattering angle v Average number of neutrons produced per fission t|> Average lethargy gain per collision p Density p(t) Reactivity vG Macroscopic cross section a Microscopic cross section T Neutron age Th Heat transfer time constant.(r, 2,u;t) Angular dependent shape function of neutron flux density p(r, 2,u;t) Neutron flux density 0(r,u;t) Shape function of neutron flux density *(r, u) Adjoint flux 00 o(r,u ) Angular dependent adjoint flux _Q ~ Solid angle Sub script a Absorption b Boiling C Core c. Capture xi

dc Downcomer eff Effective ext External f Fission fast Fast removal in Core inlet inel Inelastic int Internal lsp Liquid single phase M Moderator nb Nonboiling out Core outlet r Riser res Resonance s Steam sat Saturated tot Total tp Two-phase tr Transport U Fuel v Vaporization w Water 0 Steady state or at some reference point I Subcooled water II Saturated water III Saturated steam xii

Superscripts H Hydrogen 1/v Unit 1/v-absorber 1 Fast neutron group 2 Slow neutron group 25 U-235 28 U-238 ( ) Average ( ) First time derivative xiii

CHAPTER I INTRODUCTION A. Motivation In solving reactor kinetics problems, it is generally assumed that the variations of space-dependent quantities, such as the neutron flux, can be described for a representative point in the reactor. This assumption of space and time separability is justified when the reactor size is small. In recent years, large boiling water reactors, such as the Dresden reactor,* are being built. It is reported that a tendency exists in large boiling water reactors for spatial oscillation to develop independently of xenon poisoning.**(1) A series of experiments with the Dresden reactor, however, does not indicate the spatial oscillation.(2) Questions immediately arise: Is the space-dependent behavior of a large boiling water reactor significant? If so, how large is this effect? How can we interpret this effect? These questions motivated the present investigation. B. Historical Background The art of analytical studies of the power reactor kinetics has evolved from the point reactor models to the space and time nonseparable models. A brief historical review on the evolution of the boiling water reactor models will serve as a background for the present investigation. Only a few representative models will be reviewed here. * The height of the Dresden reactor core is 103.6 inches (cold), and the circumscribed diameter of the core is 129 inches. 3 ** Superior numerals correspond to the numbered references at the end of the thesis. -1

-2An early model by Beckjord(4) is a linear, point reactor model. By comparing with the experiment, p:ceper values are chosen for the constants that appear in the model. This model was used successfully in predicting the stability of the Experimental Boiling Water Reactor at high power levels. In order to include more accurate description of feedback effects, the feedback kinetics can be made space-dependent, while keeping the point-reactor approximation in the nuclear part. This provides a good analytical model when the shape of the neutron flux distribution is practically constant. Akcasu included space-dependent feedback in his model of the boiling water reactor. (6) A nonlinear model was introduced by Fleck who derived the feedback equations based on two(7) phase flow balance equations.(7) His model, however, assumes a constant void distribution form. The Kjeller model(8) uses space-dependent, nonlinear equations to describe the feedback dynamics, but the neutron flux shape is still kept constant. The space-dependent feedback characteristics are only reflected in evaluating the reactivity to be fed to the neutron dynamics The ul.timai e st-age of the e.volution is the space and time nonseparable models. The analytical approaches to the nonseparable problems can be classified in two major categories.' One is the "modal" analysis approach, where the time-dependent meft.ron flx and othed an r pertinent ensity functions are expressed in terms of a. series of eigenfunctions with time-dependent coefficients. Each eigenfn'lction is associated with a given spatial mode, starting with the fundamerntal mode and going on to the progressively higher modes.

-3The other approach is the "nodal" treatment of the problem. Here the reactor is divided into sections and a set of kinetic equations is given to each section. With proper inter-section coupling, the system will be solved for a given perturbation. The recent works of Crowther and Wolf, (1011) and of Nahavandi and Von Hollen(l2) are of this category. They each programmed for a digital computer a set of time- and space-dependent diffusion equations. Since the convergence of diffusion calculation is generally slow for large reactors, the calculation takes a large amount of time, even when the feedback loop is open.* Therefore, a more practical approach to the space-dependent boiling water reactor kinetics studies is desired. C. Adiabatic Approximation 1. Basic Concept The adiabatic approximation** assumes that the flux distribution. the reactivity and other nuclear properties of a reactor at an instant of time are approximated by the steady state values determined by the distributions of the reactor constants at that time. The reactor is made critical by fictitiously changing the average number of neutrons born per fission. The time-dependent behavior of the reactor is thus approximated by a succession of steady states. * In one case the computer time was 5 minutes on the IBM-7094 for the real time of 0.17 second.(12) ** The term "adiabatic approximation" has no relationship with thermodynamics. Rather, it is used in quantum mechanics as a method for solving the Shroedinger equation when the Hamiltonian depends on time. If the Hamiltonian is assumed to change infinitely slowly, and if "the system is initially in an eigenstate of H(0), it will, at time t, have passed into the eigenstate of H(t), that derives from it continuously. This is known as the adiabatic theorem."(13)

-4In 1949 Hurwitz suggested the possibility of this approximation. In his derivation of the reactor kinetic equations, he showed that the time- and space-dependent fission rate can be written as a product of two functions.(4) One is a function that depends only on time; the other depends on time and space, normalized in space. When the properties of the reactor are changing slowly compared with the prompt neutron generation time, Hurwitz showed that the latter of the two functions can be approximated by the neutron flux distribution due to the prompt neutrons only. More recently, Henry(l5) derived the reactor kinetic equations from the transport theory by using the static adjoint flux, and laid a more rigorous basis to the approximation made by Hurwitz. in Appendix A, Henry's derivation of the reactor kinetic equations is reviewed. 2. Mathematical Foundation The space- and time-dependent neutron flux density cp(r,_,u;t) can be written as a product of two functions ( (_r, _,, t) = (r. l;t) N(t) (1-1) where 0 is a space- and time-dependent function and N is a function of time only. If the reduced fluxes are defined as #(ri;t)= (r L,t) d (1-2) < ) =J fl ^ rL au)c 2 d^(1-3)'r.

-5where 4r is the time-independent adjoint flux,(16,17) and if the fluxes are nomalized such that At Ji ro L' A ) di(, C d 3c = 0 (1-4) u r where v(u) is the speed of neutrons with lethargy u, then the timedependent transport equation for cp reduces to ~_~_I,> =?(~- ~..!~_/Xl,:(t^+ 7C(.Cdt A)t = A(t)and the precursor densities, Ci, satisfy d t A (t) The details of the derivation of Equations (1-5) and (1-6) are found in Appendix A. It should be noted that Equations (1-5) and (1-6) have the same forms as the usual reactor kinetic equations.(18) N(t), the timedependent part of the neutron flux, is called the time function and D(r,_,u;t), or 0(r,u;t), is called the shape function. The shape function 0(r,u;t) and the accompanying quantities p(t), D(t) and A(t) can be obtained by the adiabatic approximation. When the reactor properties are changing slowly, it can be assumed that the neutron flux shape and other nuclear quantities, i.e., p(t), P(t) and A(t), which are caused by the variations in the core properties, can be given by the steady state values based upon the reactor conditions at

-6that time. This is possible since the shape function is determined mostly by the prompt neutron distribution which adjusts itself to the changing reactor conditions in a few prompt neutron lifetimes. In the adiabatic approximation, the space-dependent effect of delayed neutron precursors is neglected. It was pointed out by Kaplan and Margolis that the approximation becomes poor for a very large reactor, since the effect of delayed neutrons becomes important in large reactors (in terms of the neutron migration length). (9) Yet for reactors of moderate size (20 migration lengths or less), this method proves to be very useful. 3. Method of Application In evaluating the approximate values of 0(t), p(t), T(t) and A(t), it is possible to perform the steady state diffusion calculation at every time interval. As shown later in Chapter IV, computer time is quite long even in the steady state diffusion calculation, and it is preferred to avoid such computations. In the present investiga(20,21) tion, therefore, a conventional method adopted by Henry and Curlee will be followed. In this conventional method, it is not necessary to perform the diffusion calculation during the reactor transient simulation on the computer. Instead, the shape function and the accompanying functions are given as the interpolations of some pre-calculated values. This is accomplished in the following steps: Step-I Choose an arbitrary set of reactor conditions that the reactor may assume during the transient. They are, for example, temperature distributions or void distributions of a specific reactor to be analysed.

-7For each of the assumed space-dependent reactor conditions, a steady state flux distribution is calculated. For an arbitrary reactor condition, the reactor is not at the steady state, but the average number of neutrons produced per fission, v, may be adjusted so that the reactor becomes critical. The fictitious v is related to the reactivity. Other pertinent quantities, namely 3 and A are also calculated. Step-II In the above, 0, p, P and A are calculated only for a limited number of reactor conditions. In order to obtain 0, p, and A for any reactor condition that may occur during the transient, express these quantities as functions of the reactor conditions. These functions give the values of 0, p, P and A for any reactor condition as interpolations of the values already calculated in the previous step. In this process, the reactor is divided into several regions, and a power series fit is made of the desired quantity. The independent variables are the reactor conditions of these regions. Step-III Solve the kinetic equations for the time function N(t), using a group of time-dependent equations in which the reactor conditions at an instant of time yield 0, p, 5 and A. p, B and A in turn yield N(t). The space- and time-dependent reactor power, N(t)> 0(r,t), allows computations of the reactor conditions through the feedback equations. The space-dependent reactor condition changes with the spacedependent reactor power, and the new reactor condition gives 0, p,, and A that will be used in the next time-step.

-8D. Present Research The objectives of the present research are (1) to develop an analytical model for the space-dependent boiling water reactor kinetics by the method of adiabatic approximation, (2) to solve space-time problems of a boiling water reactor by this model, and (3) to investigate the importance of the space-dependence from these results. In the present investigation only the axial space-dependence is considered, since the void transit phenomenon along the core axis may have a significant bearing upon the space-dependent behavior of boiling water reactors. Chapter II develops an analytical model that describes the hydraulic and the thermal behavior of boiling water reactors. The reactor core is divided axially into a number of sections. Balance equations of mass, energy and momentum for the two-phase flow are integrated in each of these regions. From the balance equations of mass and energy, equations for the void fraction and the water velocity in each region are derived. These are new expressions, derived as an extension of Fleck's model(7) of a one-region reactor to multiregion reactors. The integration of the momentum equation around the core loop yields the inlet water velocity. Levy's friction coefficient for the two-phase (22) (23) flow() and Haywood's correlation for the slip ratio3) are used in the numerical calculations performed in Chapter V. Brown's formulation for the vessel pressure dynamics(24) is adopted in the present model. The heat transfer from the fuel to the coolant is described by a model with a single time constant. The dynamics of the boiling boundary is derived also as an extension of Fleck's model.

-9Chapter III is devoted to the formulation of a space- and time-dependent analytical model of boiling water reactors by the method of adiabatic approximation. The three steps discussed in Section C of this chapter are followed here for the axially space-dependent boiling water reactor system. Numerical calculations are performed for a specific reactor, but the present technique is applicable to any boiling water reactor. In Step-I, over one hundred void distributions were assumed for a sample reactor, and the flux shape and the other nuclear properties discussed in Section C are calculated, assuming that the reactor is at the steady state. In order to relate the flux shape and the reactivity to the void distribution (Step-II), a stepwise regression computer program(25) was found to be a powerful tool. In Step-III, the void to flux shape and the void to reactivity correlations obtained in Step-II and the feedback equations developed in Chapter II are combined. A computer program was written in MAD (Michigan Algorithm Decoder)(26) to carry out the numerical work. The program was intended for the solution of the sample reactor behavior, but can be easily modified for any boiling water reactor. Chapter IV discusses the problem of the adequate division of the core. The more regions the core has, the more exactly the model can describe the system but the more technical difficulties are encountered. The problem was investigated numerically for the sample reactor. For this particular reactor, a division of the core into five regions gave a satisfactory response compared with the models with more number of regions.

Chapter V examines the space-time behavior of the saiple reactor by using the computer code developed in Chapter III. External perturbations are introduced to the samp2e reactor which is otherwise 7naintained at the steady state. The in-vestigated external disturbances xtere (1) a large step increase in reactivity, (2) sinusoidally varying reactivity perturbations and (3) sinusoida;lly varying v-erical accelerations. Little space-dependence in the flux va-riation was observed in the latter two calculations. Chapter VI gives conclusions to the present i nvtigation.

CHAPTER II FEEDBACK DYNAMICS OF MULTIREGION BOILING WATER REACTORS A. General Remarks The reactivity of a reactor is affected by various internal causes. To mention only a few of such causes, they are the content of void in the moderator liquid, the temperature of the moderator, or the content of neutron poisons in the reactor core. These are in turn dependent on the reactor power, and the reactor power changes with the reactivity change. The cause of the internal variation of reactivity is referred to as the feedback mechanism. In boiling water reactors, the main contribution to the nuclear properties is from the steam void. In the present investigation, therefore, only the effect of the steam void is considered as the cause of the changes in the nuclear properties. Furthermore, a small contribution due to the subcooled boiling is neglected, and only the bulk boiling steam void is considered. A set of time-dependent differential and algebraic equations that describes the feedback dynamics of the multiregion reactor core is presented in this chapter. Of these equations, void fraction equations, water velocity equations and water temperature equations are originally derived for the multiregion core as an extension of Fleck's (7) equations for the one-region core.(7) The slip ratio correlation, the loss factors of two-phase flow and the equation of the vessel pressure dynamics are borrowed from other authors. -11

-12B. Basic Equations of Two-Phase Flow A set of partial differential equations can be written for the two-phase flow.in a channel. They are balance equations of mass, energy and momentum. These basic relations were investigated by many authors, (7,27,28,29) If it is assumed that (1) the flow is one-dimensional in the core, (2) the kinetic energy of the flow is negligible for the energy balance, and (3) the flow is continuous, the balance equations of mass, energy and momentum become, respectively a, P d-o-a1 P+ { (I- f)v + -v5} =0~ (2-1) at t a iPt -^eW(]-f^ [^ -)V} +99x ttha~c -fv wteh,;f9 5}= (2-2) ax -^ R IPw (- f )+ PS (2-3) where e = internal energy f = void fraction h = specific enthalpy p = pressure q = heat power per unit volume

-13t = time v = fluid velocity R = friction resistance p = density The subscripts w and s indicate water and steam respectively. The above equations are given for the two-phase mixture, and the relation between the two phases is not obtainable from these equations. It may be had by writing two momentum balance equations, one for the liquid phase and the other for the vapor phase. However, an accurate knowledge of the shear force between the two phases is essential. Beckjord and Harker give formulas for the shear stresses between the two phases, but they were checked only with static experimental results.(30) Another way of getting the relation between the two phases is by the slip ratio defined as vS y _ Vs (2-4) The values for 7 are measured by many authors and found in the (22,23,31,32) literature. (3 ) The slip ratio correlations are discussed in Section D. C. Integration of Basic Equations If the inlet water velocity, position of the boiling boundary and the power distribution are given, the void fraction and the water velocity in the boiling regions are determined. In this section, equations for the void fraction and the water velocity are derived.

-14Integrating Equations (2-1), (2-2) and (2-3) from xi to Xi+l, t e' —^-f+Psf.t- T I- tw,ip V (2 5) d-t I W S es AC Wtl[ Wi i l A'' ^c <1~ P~ e~ f I —+PP - -f )v (2-6) -f) AW+ Vs dx = It li 1+1 lStlj w( x~^1~~Ax - (P - )t F; + Gi (2-7) where r Xii Qi - $ A- dx Fi = -AcJ R d x Xi X-L'1,

-15The above equations are applicable even if the region contains a non-boiling portion. The void fraction, f, simply vanishes in such a portion. In the above the friction term Fi is not given in an explicit form. Should a reverse flow occur, care must be taken in interpreting this term. The momentum, which is a vector, has a positive direction in the general steady state flow direction. 1. Void Fraction If a trapezoidal void distribution is assumed as shown in Figure 1, the void fraction at x becomes, in a totally boiling region, ~~~~F (x = - rr+ +.....- X < X < X< - Aii ~- 1 + X;'.1AT - x- <x- (2-8) and in the incipient boiling region, 0X. X Xb \( + Ix) — = - (2-9) -- X* — --- - -- -- x X z, where xb is the position of the boiling boundary, and tb -Xi( l -Xb (2-10) Then Equation (2-5) becomes, in a totally boiling region,

-160; B o 4-P -4- H 4 - H H O4- HO H 0 + + 0 I.H.*H. r- *-:.r'l'*H +4~~~~I *-b0 ~~~r-I P4'C1 t x H 74 W1 AW1~~~~~ 0 0.,-I *rl 0) S 0 *H ).1 p L XH Q C^ I o s ^43o I H I o - H d Q I V- +o a, oD c 0 -hD V 0 dX *H 0 00_ H 0 E- O * H 1 0' ~ X $ F4 5 ~) V4- -.'H CH -H Ie- -

-17t ii t (l + V (l )} + t I( + it -(t i) (I2I-11) Similarly, in the incipient boiling region, ~f = 2 |~VWi fI-f VW)vwi +i-I {Ab+ + t,,,1,,V (2-12) Equation (2-6) can be written as follows by using Equation (2-8) for f(x) and Equation (2-5) for the term f px(l-f) dx xi Pw If the region is totally boiling +h Ps 2 2 ) L (PSff A + v-0_5 t(_ +, - -W_ _1) (,)vt wP Rss LiP,~tv/t-itw t (2- 3PS 2 Ad _- K^- - (2-13) As L- 1

-18From Equations (2-11) and (2-13) fIt, - -l+ Acf i 1V (2-14) Ac Ps ~v It j,,, where Q' = Q;-A<A-^h ^ [(l-f;)|] + Pw[-P} PS f, ~+ Ps;l- w, fP} (2-15) Similarly, in the incipient boiling region, 2 | b Q i t (/t - Z ) | (2-16) * fi1 &6 IAc Ps h, ~l WitL 2I where Q' is the effective power given to the boiling portion of this region, and is approximately,'^ ^^, (2-17) It is seen that Equation (2-16) is the previously obtained equation (7,33) for the one-region reactor. 2. Water Velocity In Equations (2-14) and (2-16), water velocity at Xi+l is not specified. Expressions for vw will be derived here both for "i+l

-19the totally boiling region and for the incipient boiling region. If Equation (2-5) is multiplied by ew, the result is I i - - - x I fr t i A -; sj - -t- e-.-' (i - *;' l -i -e {j | i,.,;, - - < i a- as e. -. X,. - --;! ( {-t ^) v ^' (v <.- t-....-+ | (2-18),. i' "-.~'V In Equation (2-18), it is assumed that the densities and enthalpies of the two phases do not change with position in the core, but their representative values may change with time. In the boiling regions this assumption was proved to be good for reactors at high pressures. (2,33) For the sake of simplicity, the same assumption will be employed in the non-boiling regions. It will also be assumed that the subcooled boiling is negligible, and that all the heat power released in the boiling region transforms water into steam. By subtracting Equation (2-18) from Equation (2-6), and multiplying by ew, the following equation yields.'j... _ - i t ^., a (, 91 x - - * - +. - X);>', l_, i - *..- 1- - (-i t -. - I. I --''' -Ky - t )V,, - -,),V 3 i 7;>^ - - * L^

-20Similarly, by multiplying Equation.(2-5) by e S.fXL+ti{ ^'r 6( e/s tl x,-f} x- e( +f —f) *'x^e I -' x?'-dP J 1 jii~w it 1 wl ~- f-a, V vs - fV, (2-20) From Equations (2-6) and (2-20), 4 f ficfl+ (i-{)(er-e5) cId + peJ.~~~~~~~~ xi^ fs + 1(-)d ( es ew +(h)efs) d ) P f~w -es kV vvt.I - (i f~~)vwjj+(f.. vs~i~ v5~ (2~22)(+ f VPf-t (2-21) Subtracting Equation (2-19) from Equation (2-21)), (Pw^-P~r)i kit(es-e,,f d/ 4-,^ (.e)~ dx f - PJ (e,.-e,)P~dx+(^es-p J ( (-es-( w xLe PS~~~~~~~X (,-tl, W Xr + ),,, t Pe )d ew.-f,)^] (f, -f ) j (2-22) X Ps, O f)~~I S,

-21Rearranging v tf-+l( i ( h Ac ( )-,, (T-' vi (2-23) If the boiling starts in the interested region, fi vanishes and, ~ t,-) h, PShvA, PW Vj (2-24) i'W it Q' Q A P\ \ D. Slip Ratio Slip ratio is defined as the ratio of the steam-phase velocity to the water-phase velocity. It depends on the conditions of the two-phase flow, namely the system pressure, steam quality, flow rate of the two-phase mixture, and geometry of the flow channel.(31) These conditions determine the flow pattern of the two-phase flow, and it is the flow pattern that is directly related to the slip phe(34) nomenon.( Accordingly, the slip ratio differs in vertical and (23) horizontal channels. Several authors proposed correlations of the slip ratio to (22,23,31,32) the flow conditions. Any one of the appropriate models may be used in the present study. Since one of the objectives of the present study is to investigate the effects of vertical accelerations, the model by Haywood, (23) et al., will be utilized in the calculations of Chapter V. In fact, only this model gives distinction between

-22vertical and horizontal channel flows. As an approximation, the horizontal case may be considered as the zero-gravity case. The slip ratio in the horizontal channel is close to unity. This is in agreement (35) with the assumption made in a previous study. ) Haywood, et al., did not explain the cause of the difference between the two cases, but mentioned that it may be because of the buoyancy effect. For varying gravity, the slip ratio may be obtained approximately by interpolating or extrapolating these two cases. Haywood gives the slip ratio in the following form: I t KX X S 0.1: =< J (.2-25) 1 O.I K0 X >. I where K is a function of pressure through a linear combination of the specific volumes of steam and water at saturation, and X is the steam quality. At a pressure of 1000 psig,'7.5 f-or hor;on'ital chdaynel [ /13. f(-or vertiCcal chanrne (-26) E. Inlet Water Velocity The inlet water velocity is governed by the momentum balance around the flow loop. For natural convection, the loop consists of the reactor core, the riser and the downcomer. In Equation (2-7), the left hand side is the time rate of change of momentum in the i-th region,

-23which is represented by the gain and loss terms on the right hand side. In order to evaluate the circulation velocity, the momentum balance of the entire flow loop is considered. Let mi = momentum in the i-th region per unit area mi = (-4)V ^ Vs d (2-27) then Equation (2-7) becomes - =- it - ) Vs2 - (iI- f )- V - v v -f } (-M Pk) (-f (2-28) I Ii _I AC If the momentum of the entire core is represented by Mcore, or -, Mcore = mL (2-29) i = then, adding Equation (2-28) for all i,,T VL COYe ), v sit "c"e = l(' - 1 V < - I - N; WNH Ni ~- (PN,-913 ~ A - P E (tl-f;)h (2-30) -I AC IIt should be noted that the pressure difference term (PN+l - P) in Equation (2-30) is evaluated inside the channel of the

-24reactor core (see Figure 2). If the same term is evaluated from outside the core, i.e., from the momentum balance in the riser and the downcomer, -it can be written as P - PN+l = - (core outlet expansion loss) - (riser outlet expansion loss) - (riser friction loss) - (core inlet contraction loss) - (downcomer friction loss) + (density difference buoyancy gain) - (inertia pressure drop in the downcomer) - (inertia pressure drop in the riser) (2-31) The last two terms in Equation (2-31) can be interpreted as the rates of change of momentum in the downcomer and in the riser, repectively. Thus if the momenta are denoted by mdc and mriser' Equation (2-30) reduces to the momentum balance equation for the entire flow loop. Defining Mtct = M 4P-. " +iser md (2-32) cor! "Tiser d.c ( 2 Equation (2-30) takes the following form: Ji^-~ o f a } - (frictioal loss) - (rictio loss) - v (frictional losz)- (frriser i=l - (frictional loss)ac - (contraction iand. expansion losses) + (buoyancy due to density difference) (2-33) The friction loss terms in Equation (2-33) will now be examined.

-25fNh RISER N+I TRAPEZOIDAL NN^X VDVOID XN N DISTRIBUTION DOWN- I COMER \ i CORE Xb X2 —- 2'p Figure 2. Pressure Difference in the Flow Loop.

-261. Non-Boiling Frictional Loss The frictional pressure loss is well established for isothermal one-phase flow. For the Reynolds number Re ranging between 0.5 x 105 and 5 x 105, and for the smooth channel wall, a correlation(36) - 0.2;50 =: x 0046 ( Re ) (2-34) will be used in the pressure drop equation. d.x.);5 2 ^ x ~e (2-35) ISO 2?~ De where G = mass flow rate De = hydraulic diameter of the channel. The pressure drop for the non-isothermal flow may also be given by Equation (2-34), with a modification due to the temperature difference in the flow profile. For the non-isothermal flow, _| /wv, but k tem pe-ra tute i non-so,/, woll temperature where = dynamic viscosity of water w and n ranges between 0.14 and 0.4. Small value of n is more widely used.2)

-272. Two-Phase Frictional Loss Levy derived an analytical correlation for the two-phase flow (22) frictional factor based on the momentum exchange model. ) The twophase pressure drop per unit length, (, is given by the following. (GP) -M(j) Id (2-37) ( X tp ( X ATr d (tp ) where f is the steam void fraction, and - is the pressure dx ~sp gradient of the liquid phase, if it were assumed to flow alone. /~~] d' x:(Ii)2. (2-38) where f sp is the friction drop factor for the liquid phase alone. Equation (2-34) can be substituted for fs if the Reynolds number for ~ sp the liquid phase alone is used. 4X ^ ^<xd9. 6(R'q::= 4 9 ^ o.b C+(hvO ) e4(2-39) Thus, for example, the friction loss term for the i-th region of the core can be taken as A 0^0 v (R J;' L )2 (2-40) Ac where R is the hydraulic radius. H

-28In the above it is assumed that the direction of the flow is in the positive or upward direction in the core. It may occur, during the transient of the reactor, that the flow is reversed in some regions of the sectioned flow loop. In this case the resistance force is in the positive direction and a negative loss occurs. If such a reversed flow is possible, the computer program should be written to handle the situation properly. F. Boiling Boundary Dynamics The boiling boundary is defined as a position in the coolant channel at which the bulk boiling starts. The bulk boiling starts when the water temperature reaches the saturation temperature. In a nonboiling region, water temperature is described by Equation (2-41) H~ i aT(x,,'t) A,1(Xt) D i T(X\ t) + T y ( at = Ac xt) (2-41) H at < W J where T(x,t) = temperature of the coolant at position x and time t Hw = heat capacity of the coolant per unit length of the core. If a trapezoidal temperature distribution is assumed in the non-boiling region, T= Tt + (T,,-T, ) x, < x < xi, (2-42) T Ti 1

-29then -3 =T - o 1+ -X — atl (2-43) at i at a, at Substituting Equation (2-43) in Equation (2-41) and integrating over Ai 3TI+1 2 ( w ( T.1 9 T at A H2 vwl t ( 1 -) (2-44) Equation (2-44) determines the temperature at each node in the nonboiling region. In the incipient boiling region, T T t X- X ( T t-Ti) -nb (2-45) By using Equation (2-45) in Equation (2-41) and integrating over Anb - Equation (2-46) is obtained, which together with a suitable initial condition determines the behavior of the boiling boundary in the reactor core. d /,2 /s D Tst D dt nb =- ~ - 1 - 2 I t+t + 2 VW1 (2-46) G. Heat Transfer from Fuel to Coolant The thermal power generated in the fuel rods is transferred to the liquid coolant if it is subcooled, or used to convert water into steam if it is saturated. Heat is transferred by conduction through

-30the interior of the fuel element to its outer surface. In simplest terms it can be considered as being transferred to the bulk of the coolant across a relatively stagnant thin liquid film. Since the time constant of heat transfer in the ceramic fuel rod is rather large and the boiling heat transfer is rapid, the heat transfer time constant through the film is often negligible.(32) It must also be pointed out that the boiling region may probably extend to the subcooled region because of the likely existence of subcooled boiling. The heat transfer coefficient in this region is generally about as good as in the bulk boiling region.(3) In evaluating the large heat transfer time constant of the fuel rod it is inevitable that some uncertainties will enter, since the conductivity of uranium dioxide depends largely on the temperature and is an experimental value which differs from one author to another. Also some assumptions such as the uniform distribution of power in the fuel rod and the neglect of heat transfer resistance across the interface between fuel pins and cladding, which is assumed to be uniform, are embodied in the analytical representation of the heat transfer model; otherwise the model becomes too involved. Thus, in this study, a model that has a single time constant for the fuel heat transfer is used. This value can be varied within the physically reasonable range if necessary. In the bulk boiling region, therefore, it will be assumed that the heat generated in the fuel rods is brought into the water with a certain heat transfer time constant which is a characteristic

-31of the fuel element. In the subcooled region it will be assumed that the fuel surface temperature is at the saturation temperature, since there is usually subcooled boiling and the temperature drop across the film is much less than across the fuel rod. The amount of heat transferred to the water raises the temperature until the bulk water temperature reaches the saturation point. Bulk boiling starts at this position. The transfer function of heat transfer inside a cylinder is well established in a simple analytical form when some assumptions are made. 37) With the assumptions that 1) heat conduction in the axial direction is negligible, 2) heat capacity of the cladding is negligible, 3) power density distribution is uniform over the cross section of the rod, 4) physical properties of the rod are constant, and 5) heat transfer resistance between fuel pins and cladding is negligible, ransfer function of the power given to the coolant due to the generated in the fuel rod becomes (2-47)!,,', h - l I + -C1 s where Aq(s) - the Laplace transform of increment of the power given to the coolant An(s) = the Laplace transform of increment of the power generated in the fuel rod

-32Dh = a constant dependent on the properties of the rod, and ~ Dh = 1 h=l Th GU k = thermal conductivity p = density c = heat capacity, and rUOh = roots of Equation (2-48) Jo ( ) ~ (2-48) where rU = fuel rod radius. When only the first time constant is taken, the transfer function becomes,,, c _(s) I /\ n () 1 _ + (2-49) n (s) St +T S or, in the time region, jt)-(tI (2-50) rjt~~

-33H. Pressure Dynamics Although the effect of the local pressure gradient on the overall system behavior may be negligible in a boiling reactor with high operating pressure, the time rate of change of the operating pressure itself may have considerable effects on the system. This effect was (27) investigated in the balance equation form by Kanai, et al. In the present model, the pressure effect is included in the modified power term of Qi of Equations (2-23) and (2-24). This can be rewritten, if it is assumed that the densities and the enthalpies are functions of pressure only,, p scpl - dP (2-51) Lw i. Q.*F;~ 1 -[j dt D and Q, becomes dependent on the time rate of change of pressure. The pressure variation can be derived by the balance equations, and by using the first law of thermodynamics for the entire reactor system. Among a number of these models, (4,2433) Brown's (24) model will be used in the present analysis. It must be remembered that all of these models use rather crude assumptions and generally the effects of the plant load are not considered. The resulting equation for the system pressure is

-34-m - - 4^-~ —--- n ____ ( i -J. h ____ d t b ~ )(s b - Y. 9) I- N - t V Ms DV_ (2-52) The details of the derivation are given in Appendix B.

CHAPTER III FORMULATION OF MATHEMATICAL MODEL BY METHOD OF ADIABATIC APPROXIMATION A. General Remarks In Chapter I, the theoretical foundation of the adiabatic approximation was discussed and the steps to be taken in the conventional adiabatic method were indicated. These steps are followed here for the axially space-dependent boiling water reactor. The numerical results contained in this chapter are for a specific sample reactor whose characteristics are discussed in Appendix C. The same technique is readily applicable to other boiling water reactors. For the sample boiling water reactor a five-region model is assumed. The adequacy of the five core regions is examined in Chapter IV. In Step-I a set of arbitrary void fraction distribution is chosen and for each void distribution the flux shape, 0, the reactivity, p, the effective prompt neutron lifetime, A, and the effective delayed neutron fraction, P, are calculated by a one-dimensional code. The latter two quantities are found insensitive to the void distribution. Therefore only the shape function and the reactivity are expressed as functions of the void distribution in Step-II. The reactor kinetic equations, the feedback equations and the correlations obtained in Step-II are solved simultaneously in Step-III. A computer program that does the numerical solution of these equations is written for Step-III, Actual solutions to various problems are carried out in Chapter V. -35

-36B. Step-I -- Steady State Diffusion Calculations 1. Choice of Void Distributions An arbitrary set of void distributions is chosen according to the following simple rules: 1. The void fraction takes only discrete values of 0, 10, 20,... 60 per cent, and is assumed constant in each region. 2. There is only one peak in the void distribution. 3. No void is assigned to the first or the bottom core region. 4. The void in the second region may not exceed 30 per cent. 5. The difference between the void fractions of two adjacent regions may not exceed 20 per cent. 6. The total of the void fractions in the four regions is not less than 50 per cent. One hundred and seventeen cases satisfy the above conditions. Figure 3 gives the average values of these void fractions in the four core regions. The standard deviations of the assumed void fractions in the four regions are also shown in Figure 3. It is desired that the void fraction distribution during the transient will lie in the range covered by the chosen set of the void distributions, because in the transient calculations of Step-III the flux shape and the reactivity are evaluated as interpolations of the values calculated for the chosen set of void distributions at this step. 2. Diffusion Calculations Diffusion calculations were performed for the sample boiling water reactor by the FOG code.(38) In the two-group approximation, the diffusion equations become:

-37200i 5 160 - 4 c.) 120- 0: oo m m w0. 0: I 4 / VOID FRACTION Figure 3. Mean of Void Distributions'Used for Steady State Calculations. Data Spread Is the Standard Deviation.

-38(3-1) t ot( ( ) D2 v2 +1jj42 -= P2I An(3-2) where tot D Bi+tast;tj D2 (' )2Z + ~2~ ~ ~ 2 B = transverse buckling X = eigenvalue and the superscripts 1 and 2 indicate fast and thermal groups, respectively. The transverse bucklings are obtained from the radius of the reactor core adjusted for reflector savings. The core constants for the sample reactor are given in Appendix C as functions of the void fraction. Essential output data produced by the FOG code are the group neutron fluxes and the eigenvalue. In the present method, the condition

-39of Equation (1-4) requires that the flux be normalized according to Equation (3-3). i ) Lr J -du d4r = c s (3-3) (u) a r In the sample boiling water reactor the adjoint flux used for the normalization was the one for the void fraction distribution of 0, 10, 25, 30 per cent in the second, third, fourth and fifth region, respectively. Figure 4 shows the average distribution of the normalized flux shape functions based on the chosen void distributions. As expected 3 and A were found very insensitive to the variation of the void distribution. Thus the values for the reference distribution A = 0.0000650 (3-4) P = 1.211 P (3-5) are used throughout this study. In evaluating these values, the thermal neutron speed was estimated to be 4.0 x 105 cm/sec. The latter figure is based on the cross section of the unit 1/v absorber* in the interested range of the void content. al/ for the hardened neutron spectrum is about 55 per cent of the cross section at 2.2 x 109 cm/sec neutron speed. * A unit 1/v absorber is a conceptual material whose cross section obeys the 1/v-law, and is one barn for neutrons with a speed of 2200 m/sec.

-40200 160 120 I — LJ IJ 0r 80 - 0 40 0 - 0 0.1 0.2 0.3 0.4 FLUX SHAPE Figure 4. Mean of Flux Distributions Obtained by the Diffusion Code. Data Spread Is the Standard Deviation.

-41C. Step-II -- Correlation Process In Step-I a set of void distributions was chosen and for each void distribution the flux shape and the reactivity were computed by the diffusion approximation. These values were obtained only for an arbitrarily chosen set of void distributions. It is necessary to obtain the flux shape and the reactivity for any void distribution that may take place during the reactor transient. They can be obtained as interpolations of the values already calculated for a finite number of cases. One way of carrying out the interpolation is by the regres(25) sion analysis technique.(2 In this section the regression program is described and the correlations of the shape functions.i(i = 1,2,..., 5) and the reactivity to the void fractions fi(i = 2,3,4,5) are described. 1. Stepwise Regression Program The correlations of the flux shape and the eigenvalue to the void distribution were established by the REGRESSION code. (5) The program provides the best estimate of the dependent variable as a linear combination of functions of independent variables in much the same manner as the least mean squares fit. The functions of the independent variables for this subroutine must be programmed by each user, In the present study, these independent variables are the void fractions in four regions (since no void is assumed in the first region), and the dependent variables are the flux shapes in five regions and the eigenvalue.

-42The regression program deals with a single dependent variable Y and a set of independent variables X1, X2,..., XK. The objective of this program is to produce a relation of the form Y = B + BlXl + 3 + + BKXK (3-6) The coefficients Bi are so determined that the expected value obtained by Equation (3-6) fits the actual value of Y in the sense of the least squares method. If it is expected that the dependent variables is a function of, say, X1, X2, X1 X2... etc., then the independent variables can be changed in the following manner. - X zS -x * Z's rather than X's are now considered as the independent variables. Therefore, whether or not the conversion option is used, the assumption of Equation (3-6) is the general form of the stepwise regression program. In forming Equation (3-6), only the significant terms are retained in the equation. Coefficient of insignificant terms are put to zero. Following steps are taken in the program. First, it looks

-43for the Xi which alone can best predict Y. In other words, it chooses an independent variable which has the largest correlation coefficient with Y. Suppose that the best candidate for Equation (3-6) is X3, then the first predicting equation is Y = B + B3Y3 o 33 where Bo and B3 are found by the least squares method. In succeeding steps, the program sorts X's into two sets Xi,1 and X 2 The X are the terms which are already included i^l i, i,l in the predicting equation, whereas Xi2 are not. For each one of 1,2 Xi 1, the program calculates an importance factor, which is a measure of the relative contribution of the variable to the predicted equation. If the variable with the smallest importance factor is less important than the user's prescription, this variable is eliminated from the equation. When all the importance factors exceed the minimum value required for retention of the set Xi.l' the program examines the set X. 2' The potential importance factor, the importance factor that a variable in Xi 2 would produce if it were included in Xi 1, is calculated for every member of the set Xi,2. The variable with the largest factor is considered the next best candidate for the inclusion in the predicting equation. A test is made whether the probability of inserting this variables incorrectly is less than the user's prescription. If the candidate has passed the test, it is included in the predicting equation. If so, the program goes back to test the significance of Xi. It should be noted that a variable once entered into the equation may be excluded from the equation, if a better combination of variables is found at a later step and the old variable loses its significance.

-442. Resulting Equations It is assumed that the mean flux levels of five regions and the reactivity are functions of the following form: /Y _ SBo > I) 3^,(3-7) L }2(i)+1 ( - ~ ~tt ^ -2(l A-)t,- / S (3-7) where fk is the average void fraction of the k-th region (k = 2,3,4,5). In order to obtain better correlations, the range of the void distribution is divided into two groups. Group I includes the cases with f5 = 0.4, 0.5, 0.6. Void distributions with f5 = 0.4 are included in both groups to make the correlations continuous. Correlations based on Group I are used when fT < 0.4 and those based on Group II are used when fV' 0.4. Table I shows the coefficients obtained by the regression code. The standard deviations of the predicted values from the actual values were 1.87 x 10-4 for the reactivity, and 1.90 x 10-3 for the fifth region flux. The latter was the largest of the flux deviations. For the magnitude of the nomalized flux distribution, see Figure 4. The maximum absolute deviation in the reactivity was 4.68x10, and in the flux it was 0.0077, which occurred in the fifth region when the average void fractions were 20, 40, 30 and 10 per crnt in the second, third, fourth and i'ifth regions, respectively. The l -er deviation corresponded to 5.09 per cent of the magnitude of the normalized flux.

-45C\ C r- C0 " -. C r'- O.' L \ n n' n C\O -.o _:\.t- \kv.O 1- H V-~ o L. H -- \0D C H n\ L C o..0 o. \o \. r-l o' - O H- co: nCO)'1 O kO 3 l O Oc O _- - riO - - I O r LC 1 OO - \o co rn \ m no o \ t00- o0 j- on _j_-p CM n M - \-j-\O H CM q CM \OO t H CM H0 O C 0cr O HC OC0 CMO - -0 HO O I! I I I I COCO!l c C:OO CC- X r ir- \ o C_:-N O CO c O ~om H - C C O\ t — C\ O II- CU _OO C- CO o N 0 C0 c n 0 -ct'\LR \ n 0O nM CO H. 10 Co 0.- m-n 0 rO- \0 c0 0 Co 0 0 H -zf >3 \ n )n m O M-1- CO O r L( O n H CN O O O -z n f 0 H a rJ 0 0 0 0 O O 0w E —I i io I Io PW; ~ -~ M tHdl — \ S \ M0'T0 — 000 O0 O Co \0 O - r- H O C0 O, CT eq C\J Qm o H N \ L( \ maLn 0 CX: 0c n o oo O O\ 0 -o0o, — c\ \ no0 co M (n \0 cO M \19 n _t co C\ N N C\J 3 o o1oo oCcXc), o o oz -o-,,L,. o C oMo0o'-J.. _.~ C~ O *-* O * * *. _ C* * O * O * * * C O O*,j 0 0 1 Q 0 0 C - 0 LfC C 0 M0 CM C CM 0O |.. - CMOOOOOO 0'0CM — 000000.- - 0 iH... * *.. EH ro1 1 0 V *, \ 0D o o o-d _C M00 O \O n 0 \0 H q > rj- oa 00 0 H N C\J L C0O Lr\ 0- 0 \ C) C)00 0,q -G t — o a1 o - - C. _- O \O -tH O 4I HT l. -1 Hc \~ - OOO r r -l - \ O CU oL-\ 0'- J -I c - - t -z o OtC O N,- - o O oC O \DL OO0I O O 0 C Ho o c-r-C 0.-."OC ~ H M. O0 E j 0 00 0 \0 0 0 H 0 0 0 \ 000 0 0 0 H - I II Eq r-C cx X O \00 CO 0 c\\ _J LA nO n H I'O LA (OOLn O O t O- O C l\ c\ m \j (n 0 a I\ c - CV t CM r -Ht- - O CV CO _:- H 0c 0 I HO O O o 0 o 0 oO rool MO O O O O H I. lo 1..1.. 1. I 1 a 1 Q) H _ r Hr -|OH r Jo 00 -CM c.- L\O L-00co o0 _O \00 - Q CG\ o i H C r* C HO.H( =1 q q m pq pq pq pq pq pq pq pq pq pq pq pq q pq pq 0 cH o 0 ^- ^ 0,~

-46Cm cO- Oo t O - oo C — co Co -\1 - L r O COC H M r O O O O flr O O _HHO O GI oo L\f0 n o o Ow - H- O 0 \O O Co n c.- -- O r- -- t —n - rl \.O. O n o L nO C n o t- r - O Co'\ r- i r-O\ o Ic- O AS IO O c O c O C0 X- H _C J O0 \-' \ o — O - L- -- -o 0O o D C\j r-IMHr —IOOOO OOLC\OOOOOOO_ Ic m co I - On O 10 rI C O \O O * ^*o __- ~O h -b-!Pi o —r —i.o o r - N cO oco M L0R_- CH-\ \ tM L- kOO r-1 \\ r H OODn rh c C\jO O c) O 01 r- 0 00 O0 r- 00O"0 CM0 C00J0r-I 0m0 r) _- mO ILr. — t n h-0 - \ O0 u'2O1 0\H CO_- O- C-)- L\ r-I O C\-Ic) O OOO HOO x OOO O O O O \00OO! i ci 0\000-J-coo CM o h- co c o0 0 H o,-!0 on 1 - 0 co O — oO 0 0- cO \C\0 C(i M'i c1 O Ln o oo ct no cO o r-o _o o- cO\'t ~- O.- O n kOh-O Ln c n -n c- m -o i od 0 o Z r CMHO 0 0 0 0 HOOQ c 0 C 0 0 0 0\ H H I Lno i m'o\ ooHo CLR-O\ HOC) —-t ( D H- OH r E-e q~~LRLC\H nH- Q\L - HO.coGD\ — CO \ CTh cCO O\ \'h "- - > O ~ — n C0 -" -CM -- C' OU H ON-O cOC Od H OO C n LC CM CM Ho H-CM 0 C n:\-L 0\-O L0 lo0 O O C)o H O 0 0 O C O 0 0tOCM O O Oo LnO I H 0 0 o —oj r t0oI I I I i I OO Go -c LNcO O cV OJ OX-N- Ln" H —- "c ON r- 0t C N-O 0 0CMG CO --- O 0 co L n Lr pr pq 01 Ao o A A A -A coo A Co Hp CO LCNL O C "DG H!-k' O.-.O KO 0n HL Do O O- Lr\ H- O L O " O C0 O O O CO OO O O O O O O O O- O-OOO O O O CH 0 a H CH A O r ckl 0) l W) F 0 ) r1 I 0 r0 r0 Ir\1O rk o k k () o - V

-47D. Step-III -- Computer Program of Mathematical Model The computer program is written in the MAD language and carries out all the arithmetic operations on the equations described in Chapter I, II and in Step-II above. Additional techniques are needed in the numerical calculations to complete the program. The program is explained in this section. 1. Principles The program consists of two portions. First, the steady state of the system is calculated. The reactor condition is assumed to be at the steady state before t = 0. At t > 0, external perturbations of reactivity, inlet water temperature, or gravity are given to the system. The time-dependent computations are the second part. Runge-Kutta's method(39) for solving the differential equations is used in the second part of this program, since the Runge-Kutta subroutine for solving differential equations (RKDEQ.) is available at the Computing Center, University of Michigan, As the step-size of the independent variable can be varied in this subroutine, a considerable time saving is possible when the magnitudes of time rates of change of variables vary in one computer run, e.g., transient simulation for a step change in reactivity. 2. Steady State Calculations Since it is not possible to solve the steady state in a single series of MAD statements, some iteration processes are incorporated until the time derivatives vanish or become small enough to be considered as the steady state. For a boiling region, Equation (2-l14) reduces to

-48(0= — - - ( yr V -YV, ) (3-8) at the steady state. The velocity Equation (2-23) is Assuming that Qi and 7i are constant, and given v i and fi, quations (3-8). and (3-9) are easily solved for fi+l an&. vi+ ZThun, if a boiling boundary, power distribution and the inlet water velocity are assumed, the water velocity distribution and the void disw-ribution are easily calculated. Then for the newly obtained quantities, the shape function and the slip ratios are re-evaluated. The above calculations are repeated using these new values of the sha-pe function and the: slip ratios. When the inlet water velocity is governed by natural circulation, the former mu st be selected so that the derivativ:e of the water momentumn vanishes for the steady state. The inlet water velocity is increased or decreased until this cundition is satisfied. When a constant pump head is added to the derivative of the momrntum, the head must be adjusted to obtain the desired mass flow rate. Other initial values are either easily obtained or given. Then the program checks if all the time derivatives are negligibly small, and proceeds O, the solution of the differential equations.

_493. Time-Dependent Calculations For the five-region reactor model the time-dependent differential equations and algebraic equations listed below are solved simultaneously. The differential equations are: 1) five equations for the void fractions in the core 2) five equations for the void fractions in the riser 3) five equations for the heat transfer 4) an equation for the system pressure 5) an equation for the total momentum 6) five equations for the water temperatures 7) an equation for the boiling boundary 8) an equation for the neutron flux level 9) six equations for the delayed neutron precursor densities. Important algebraic equations are: 1) an equation for the inlet water velocity 2) equations for the slip ratios 3) correlations for the flux shape 4) a correlation for the reactivity. The void fraction is always positive, and the subcooled temperature is always below the saturation temperature when the local pressure variation is negligible. As the boiling boundary shifts from one region to another, the derivatives of the void fraction and the subcooled temperature may no longer stay zero at the inter-regional boundary. When the boiling boundary approaches the core inlet crossing

-50an inter-regional boundary, the void fraction at that point starts growing from zero. The derivative of this void fraction may fluctuate to produce an unstable situation. This is easily understood from Equation (2-16). To prevent such situations, an amply large initial value must be given to the void fraction at the inter-regional boundary. A similar phenomenon is observed when the boiling boundary moves upward crossing an inter-regional boundary. In this case the temperature at this point starts decreasing. It can be observed from Equation (2-46) that such a shifting of the boiling boundary may create an unstable situation. In order to avoid these instabilities, an appropriate initial condition is given to the void fraction or to the temperature. Depending on the position of the boiling boundary, three b.anchings are made in the program. By defining snb - (3-10) in the incipient boiling region, the situation is divided into three cases. In the following, cnb is a small positive number. Reference is made to Figure 5. Case A: 6nb > 1 - nb nb nb When the boiling boundary lies in this range, either the boiling boundary has recently moved into this region, or the boiling boundary is about to move up out of the region. In either case, f. is given by extrapolating fi+l as Ab (3-11)' all A +Ab

-51meoI5'' — 0 - *H O H H~~ HH |U fC~j r~ F4 v I? II *H n ^i * VU Is ~ I-r4 E z^H H H I O. m < ~rl~~~~ ^s - g ~ ^ S^^ X~~~~~ M oC5c

-52Case B: nb < 6n 1-- nb nb ~ nb Since enb is a small number, there is still a possibility that the boiling boundary has just moved down into the present region from the region just above it. In such a case, an initial value must be given to the void fraction by Equation (3-11). Or it may be that the boiling boundary has just moved up from the i-lst region. In this case, the derivative of the void fraction, fil, vanishes but the void fraction itself has a certain non-zero value obtained at the preceeding time point. This must be made zero. Similar initial value settings are employed for the subcooled water temperature, if it is found that the boiling boundary entered in this region at the present time step. When the difference of the power distribution in the i-th and the i+lst regions is large, there may be a noticeable difference between the extrapolated void fraction and the "true" void fraction as the boiling boundary moves from one region to another. A small difference in void fraction may cause a noticeable discontinuity in the reactivity, since the voids are the main contributor of the internal reactivity. To prevent this situation, a transient region is set up, in which the void fraction fi+l is given as the weighted average of the extrapolated and the "true" values. Case C: 0 < nb cnb Case C is symmetric to Case A, in that the void fraction and the water temperature are reversed. Here, -the water temperature is extrapolated in contrast to the extrapolated void fraction in Case A,

-53in order to guarantee that the water temperature Twi never exceeds the saturation temperature. T,- Tsot - (Tst -T, ) (3-12) 4. P rogr a In Appendix E, the computer program is listed with references to the symbols used. The flow diagram of this program is shown in Figure 6 in a simplified form. Only the main features of the program are illustrated in this figure, A number of options are made available in the program to facilitate various uses of the program. An option card containing a vector of 16 integers, N(1) through N(16), is read before the other data cards. These options are explained below. Definitions of variables used in the program are given at the end of Appendix E. N(1) = 0 constant shape function 1 variable shape function N(2) = 0 one-group delayed neutron 1 six-group delayed neutron N(3) = 0 Pext 0 next =0 1 next = P0 sin 2~t t 0 O 2 Pext { at t 20 OPO t ~

-54CALCULATE DI START 1 <READ DATA — ~ CALCULATE SLIP RATIOS _, VW AND F VW AND F AND ONVERGE... I..L FLUX SHAPE | YES CHANGE INLET WATER THS SPEED,RISER NATURAL CIR HEIGHT, OR CULATI PUMP HEAD CALCULATE OTHER INITIAL VALUES CALCULATE _ W STEP OF RKDO TIME SOLVE <RKDEQ. (0O=~' DIFFERENTIAL...... [EQUATIONS NO INTEGRATED < EQ.O SYE VALUES ARE STORED PRINT NO IME>TL M -| AND/OR Vs/ L~ ~RESULTS PROCEED TO THE NEXT PROBLEMES Figure 6. Flow Diagram of the Program.

-55N(4) = O Pint is through a constant void-reactivity coeff icient. 1 Pint is evaluated by Equation (3-7) N(5) = 0 slip ratio is constant 1 slip ratio is variable and a function of g 2 slip ratio is variable but independent of g N(6) = O inlet water speed is constant 1 inlet water speed is controlled by natural circulation. Both HEAD and RISER are constant 2 same as above, but only HEAD is constant 3 same as 1, but only RISER is constant N(7) = 0 inlet water temperature is constant 1 inlet water temperature is varied stepwise N(8) not used at the present time N(9) = 0 boiling boundary is constant 1 boiling boundary is variable N(10) = 0 P is constant 1 P is variable N(ll) = 0 riser is treated as one 1 riser is treated in regions N(12) = 0 g is constant 1 g is varied sinusoidally N(13) = 0 constant step in RKDEQ. subroutine 1 step is varied in RKDEQ. subroutine

-56N(14) = 0 calculate the transient 1 calculate the steady state only N(15) = 0 do not print iterations in the steady state calculations 1 print the iteration in the steady state calculations N(16) =-2 plot deviations from the steady state; also print Y -1 plot and print Y 0 print Y and Z 1 print Y.

CHAPTER IV EXAMINATION OF ADEQUATE NUMBER OF REGIONS A. General Remarks In Chapter III, the sample reactor was divided into five regions without giving logical justification to this number. Here it will be examined how good an approximation the five-region model gives, and how it is compared with other core-divisions. Naturally, if the number of the space nodes is large, the approximation improves. However, if the number of nodes becomes too large, the computer time easily exceeds the economic limit, hence a compromise is required. The problem is investigated in two different parts: one is the problem concerning the nuclear characteristics and the other is that of the feedback dynamics. The problem of the adequate space nodes in the nuclear part is studied for a reactor at the steady state, since the adiabatic approximation treats the dynamic behavior of the flux shape and the reactivity as a succession of steady states. By using steady state diffusion codes, it was observed how the shape function and the eigenvalue changed as the number of regions was varied. The problem in the feedback part is dynamic in nature. A flat shape function and a constant void-reactivity coefficient were used in a computer program similar to the main program discussed in Chapter III. For different number of regions the response of the reactor was studied as the reactivity was changed in steps. -57

-58B. Nuclear Characteristics As discussed in Chapter III, the flux shape and the reactivity were solved by the two-energy group diffusion calculation. In each region of the reactor the nuclear properties were considered constant, since the available diffusion code, FOG,28) is provided with constant input parameters for each region. In examining the adequacy of the fiveregion division only one energy group was used. By using all the forty available regions and almost all the mesh points (200 out of the available 239), a reference calculation was made for a reactor with continuously distributed parameters. Figure 7 shows the linearly varying parameters used in this calculation. The choice of the parameters was rather arbitrary, but the general characteristics of the sample reactor were kept. The results of multiregion reactor calculations were then compared with that of the reference calculation. In multiregion reactors, the core constants were averaged over each region, Figure 8 is the results of these calculations. The plots are made for the average fluxes in each region. For the given reactor it is seen from this figure that five or more nodal divisions will give a good approximation to the uniformly distributed reference reactor. For a large number of divisions these average values coincide with the continuous case. Figure 9 shows the flux distribution of the two-region and the five-region approximations plotted for each mesh point. The flux distribution of the five-region reactor is almost identical to the continuously distributed case (Figure 8). The deviations of the plots from the curve in Figure 8 may therefore be attributed to the flux averaging process in

-59CORE CORE BOTTOM TOP -I I CM 0.060 0.055 0.053 0.050 0.045 - -= 0.045 CM 0.3.D 10.3 0.2 Distributed Parameter 0.1 Region Averaged Parameter (4-Region Reactor ) 0 I I I 0 50 100 150 200 CORE HEIGHT,CM Figure 7. Continuously Distributed and Sectionally Averaged Reactor Constants.

-60+3.1..p -p -S-p p < -. o o o o) ) o o Cd C)Cd C) Cd C) C). H.H.H H H H H tb hD b h hO bD bO N I8 I - =ln 3a03 O / CD ( W 0W ~ (.INn AVa8elIJ8V) xnf- NOJll.3N

-610 gw O 0 0 W~~~~~ 0 0- 0 o < w w c r z 0 z C.) O. O ~3 da'3 a o ~ ~~~~~O ~~~~~~~~~~~~~ O wo ^^ >^ u- > ILINn A8Vl.188iV'' XnRId

-62the five-region case. On the other hand, in the two-region case, the flux distribution itself differs from the reference reactor. The change in the eigenvalue for different nodal divisions is also investigated in order to examine the optimal number of nodes. Curve (a) in Figure 10 is the plot of eigenvalues with the change in the number of regions, calculated by the FOG code. It is noted from this figure that the eigenvalues for a multiregion reactor converge very rapidly to the value of the continuously distributed reactor. For more than four regions, this convergence is extremely good. Since only tne region-averaged neutron flux is used in the adiabatic approximation, and much of the information obtained by the diffusion code is thrown away in this approximation, it was investigated whether a simple, customarily used nodal solution for the diffusion equation might work as well. For this purpose a simplified diffusion code (SDC) was written in a similar way to the FOG code. Appendix D describes the code. Figure 11 is the results of the flux calculations obtained by SDC. The forty-region calculation by FOG is also plotted for reference. The deviations from the reference curve are much larger than those by the FOG calculation for tne same number of nodes. Curve (b) in Figure 10 shows the eigenvalues calculated by SDC. Although tne absolute deviation is small, many nodes are necessary to obtain the accuracy comparable to the FOG calculations. Figure 12 demonstrates that the convergence in the diffusion calculation becomes slower as the number of nodes increases. It is a generally observed phenomenon in calculations of large or loosely coupled

-63m 0 0 ~rr.rI LL W' z c r-i IIr Iy I _ 3nl 3/ _oi * 0 0 0 0 3nlVAN3913

-64Wr O 0 0 - ooo o o oX.P.p 4. p.,+i3 4..4 Co CD CD O C) C CD I dcr cd a d Pc d d (1) U ) U ) ) U U) O O O 0000 *~ r- H r rl.rH rH rl b 0D hD )hD bD bD bDO - 1 / 1 io oo C; X I>E <I / | I 0 I 1- ^ x -o o r o 8& F -. ~ Nn3N 0~0~~~~~~~~ 0 -W 0 0 ( IINn AV.LI8v) Xnfl- NoIn'3N

-65*d CM \s. o -H -4i 0 Cf) _\~c o w \ ~~~ — ^ ~0 i10H -, -< | 0 0 cd 0 0 0 0 0.H (r) l N _ dO (OGS) 3VJII I3ifldbO3

-66reactors. This result suggests that the computer time would be enormous, if tne time-dependent diffusion equation is solved by the differencedifferential method with sufficient accuracy. C. Feedback Characteristics The number of nodes must be sufficiently large to describe accurately the space-dependent feedback dynamics as well as the nuclear characteristics discussed in the preceeding section. In order to determine the optimal number of nodes for the feedback loop, the spacedependent neutron dynamics was excluded from the reactor model, and the neutron flux shape was assumed constant throughout the core. A computer program was written with a provision for varying the number of core regions. Using this program, the transient behavior of the sample reactor for a step change in reactivity has been calculated. Figure 13 shows the transient behavior of the reactor which has space-dependent feedback dynamics and uniform flux distribution. In these calculations, a constant void-reactivity coefficient of -0.001/1%-void was used. The heat transfer time constant of the fuel element was taken to be 5.0 seconds. The one-region reactor model has a triangular void distribution along the axis of the core. This gives a fairly good approximation at the steady state and also during the initial transient as the reactivity is increased in steps. At the outset or the transient, the boiling boundary moves toward the core inlet keeping the void distribution almost triangular. Later the boiling boundary remains approximately stationary and the growth of the total void fraction becomes more dependent

-67~~Xb,.~~~~~ 20-Region Reactor 2~80^^0'Qry~( 0 10-Region Reactor o C)s 5-Region Reactor >- "7 V 4-Region Reactor <'X Nt X 2-Region Reactor Z 60 3 1-Region Reactor 0 0 z + 1J 40 o + + 0 m 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 TIME, SEC 0.6 X 0.5 - rn V 0.4- _V |L+ Vx + + ~~O4~~ + -r X 12 x<160CM LL + + 0.5 1.0 1.4 TIME, SEC > + +v Figure 13. Behavior of Multiregion Reactors with Space-Independent Nuclear Characteristics. o,. ~'- + - o<x<40cH 0 0.5 1.0 1.4 TIME, $EC Nuclear Characteristics.

-68r- 0 0 014 q~~~180 0 ~ V/ 160 a: 100 a LL~- -- -_+ I ^^^/^^^ --- -+-120 0X~~~~ + 03'^. 0 0.5 1.0 1.4 TIME, SEC Figure 13. (Continued:)

-69on the increase of the void fraction at tne mid-core. The triangular approximation then becomes inadequate. Tne inadequacy of the one-region reactor model is clearly shown in Figure 13. Tne void fraction in the mid-core (80 < x < 120 cm) of the one-region model differs as much as 50 per cent from that of the twenty-region model. The results of the five-region model are almost identical to those of the twenty-region reactor model in all variables but the void fraction. For smaller number of regions the average void fraction of one region is not affected by the void fraction of the regions below this region for a longer period of time, since the void transit time across one region increases as the number of regions decreases. For the fiveregion model, this phenomenon is still noticeable, but the deviations from the twenty-region case are small. It is also noticed that in a model with less number of regions the void traction near the boiling boundary is somewhat less than that of the twenty-region model, because near the boiling boundary the void distribution in a model with less number of regions becomes more like triangular, whicn is not a very good approximation. D. Conclusion From these two investigations discussed in the foregoing sections, it is concluded that the division of the reactor into five regions will give a good representation of the space-dependent reactor. In addition to the five core regions, two reflector regions at the top and the bottom of the core are included in the present study. When the riser is used, the riser region may also be divided into regions to provide for a better simulation of the void transfer.

CHAPTER V SPACE-TIME BEHAVIOR OF BOILING WATER REACTOR A. General Remarks The space- and time-dependent behavior of a sample boiling water reactor is investigated by using the digital computer program discussed in Chapter III. The details of the sample reactor are described in Appendix C. Some physical properties of this reactor are list( II. This type of reactor was chosen since it is believed to b( the most typical boiling water power reactors. TABLE II PROPERTIES OF THE SAMPLE REACTOR Core Height 200 cm Core Diameter 200 cm Fuel Element UO0, 1.2 cm Water/U02 Volume Ratio 2.2 U235 Enrichment 1.5% Migration Area, 30% Void 99.2 cm2 Operating Pressure 1000 psig Design Power 200 MW(th) The perturbations given to the reactor are: (1) an increase of reactivity as a step function by a magnitude of 90 cents, (2) sinusoidal variations in reactivity with amplitude of 10 cents and -70

-71(3) sinusoidal variations in vertical acceleration with amplitude of 0.5g. No attempt is made to control the reactor while undergoing these perturbations. B. Stepwise Change in Reactivity The external reactivity is changed stepwise from zero to 90 cents at time zero. The reactor was assumed to be at the steady state before tne disturbance is introduced. The reactor was provided with a riser of the same length as the core height, and the initial boiling boundary was assumed at 80 cm from the core bottom. The subsequent reactor behavior is illustrated in Figures 14 and 15. The time step used for this calculation was 0.05 sec and the computer time for the 3 sec real time simulation was 3.46 minutes on the IBM-7090. From these figures, the following observations are made: 1. A large change in the flux shape is observed. Tnis is due to tne considerable change in the void distribution in tne reactor core. 2. The shape of the void distribution changed during the transient. The peak of the void distribution did not necessarily occur at tne core exit. This observation is in agreement with the analysis by Sanathanan. (29) 3. The general trend of the oscillatory behavior of the water velocity agrees with Fleck's result in that the inlet and the outlet velocities are out of phase by 180 degrees. It is observed by the present model that in most boiling regions the coolant has a similar behavior, while the coolant of the nonboiling region flows downward.

-72(.0 0 II 0, z /o X o+:lD a 0, m - q / \ - m I /" "" "" \'H'''' 11'''' / \ X 0 * * * * / \ S OOOHOi. I Mm Oo< X+^ I-) ( | --------------------- o 0 ( W** —^^ ) lH913H^ 0 310) < U) -H \V\U -H.H dC ~> _.____ \ o z 4 CM --- 00 oI> HH o

-73CJ_ _J \ J 2\ LL 0 8 \ 0 0 < / 80 L 0.2 0 2 3 TIME (sec) Figure 15. Behavior of a Natural Circulation Boiling Water Reactor to a Step Increase in Reactivity of 90 cents.

-74>.. 0x 80 z 7 70 60| 50200 / \\ - X200cm 160 - 120 a>, 1500 0 2 3 U Figure 15. (Continued LL500 2 3 TIME (sec) Figure 15. (Continued)

-75C. Sinusoidal Change in Reactivity A series of calculations simulating the reactivity oscillation test was performed. The reactor was operated with a riser of one-fifth the core height and at the steady state the boiling boundary was assumed at 50 cm from the core bottom. The external reactivity was fed sinusoidally to the reactor; Pext = 10 sin(2nft) (in cents) where the frequencies considered were f = 0.1, 0.2, 0,5, 1.0 and 2.0 cps. The time step for these runs was 0.0125 sec, and the ratio of the computer time to the real time ranged from 24.6 to 38,0. Figure 16 shows the response of the flux level as a function of frequency. It is noticed that both the amplitude and the phase are flat and the system is stable for the investigated frequency range, The phase curve shows that the flux led the reactivity at low frequencies. This phenomenon was observed in the actual experiment with VBWR (Vallacitos Boiling Water Reactor).(4) Figure 17 shows the void fraction behavior. For small frequencies, both f3 and f6 have similar characteristics; their phase lags are almost the same. The above statement does not hold tor large frequencies, since the void transient time becomes comparable to the period of oscillation. This phenomenon has been discussed by otner authors.(6Q41142) It is also observed that the amplitude of the void fraction variation is larger for the lower portion of the core. This may be attributed to the relatively high power level near the boiling boundary,

-76N(t) -N(0) N(O) 12'1~2 "AMPLITUDE 10 -60. PHASE LAG W 0 I i i I -60 0.1 0.2 0.5 1.0 2.0 f (cps) Figure 16. Frequency Response of Flux Level to External Reactivity. ext = 10 sin 2Tft (in Cents). ext

-77PHASE LAG 1.0 -' -90 0.8 -\ -150 f$ AMPLITUDE _ 0.6- -210 0.4- B e \ — 270 o 0.1 0.2 0.5 1.0 2.0 FREQUENCY (cps) Figure 17. Frequency Response of Void Fr action to External Reactivity. 0PX = 10 sin 2 0.4 - — 270 0 0.2 -330 04 0.1 0.2 0.5 1.0 2.0 FREQUENCY (cps) Figure 17. Frequency Response of Void Fraction to External Reactivity. P ext = 10 sin 2-f~t (in Cents).

-78In Figure 18, the detailed behavior of the reactor variables for the reactivity oscillation of 1 cps is shown. The following results are observed from these curves: 1. There is no difference among the reactor'lux responses in different regions of the reactor; no spatially dependent behavior of the neutron flux is observed. The observed space-independence agrees with the experimental result of the Dresden boiling water reactor.(2) 2. The variation o' the water velocity is in one direction in all the boiling regions, while in the non-boiling regions it is in the opposite direction. 3. The void traction variation is larger in the lower portion of the core where the power level is higher. The importance is large at this point, but the void-reactivity coefficient is small since the absolute magnitude of the void fraction is small near the boiling boundary. From these curves alone it cannot be concluded that the void-reactivity effect is more important in this region. If a single void-reactivity coefficient is to be used, a proper weighting is necessary to obtain the correct void effect on reactivity. D. Sinusoidal Change in Vertical Acceleration Another interesting example is the simulation of reactor behavior under varying acceleration field. A change in the vertical acceleration changes the core inlet water velocity through the momentum balance, Equation (2-33), since the bouyance term is directly proportional to g. It also affects the slip velocity ol the steam. In practice the variation in the vertical acceleration field is encountered in ship

-79I0 z 0 - ~ 05 1.o I.5 2.0 SEC -0 \.' t0.0I /SE 0.9 1 0.X0 1.5/2.0 SEC 0.9 1. \ i.0 0 0.5 1.0.I.1.0 SEC I. o SEC 0.9 1.0 0 0.5 0 15\ 2.0 SEC 0.9 Figure 18. Response of the Reactor to a Sinusoidal Change in Reactivity, Pext = 10 sin 2rt (in Cents).

-80I0 t i ~~;\ /2 SEC 0 1.002- 6W/VWs0 10*t 0.5 1.5 2 SEC 0.9981.002 - _W57VW 5. 0.5 1.5SEC 2,5 SEC 0.998- 1.002- ^_ VW4'VW4. 1.0 0 0.5 1 1.5 2 SEC 0.9981.002 - VWo7_ 1.0 0.5 I 1. 5 2 SEC 0.9981.002VW,2/VW, 0 0.5 1.5 2 SEC Figu8rre 180. ( Cont )SEC 0.998 - Figure 18. (Continued)

-810 0.5 15 2 SEC -1I0 0.001(I) SEC 0 0.5 1 1.5 2 SEC -0.001 0.001 C0 0.5 I 1.5 2 SEC -0.001 0.001AfA0 0.5 I 1.5 2 SEC 50.2 - Lnb SEC U 0 0.5 F 1.5 2SEC -0.001k 50.2 0 0 49.8 Figure 18. (Continued)

-82propulsion reactors. A previous investigation (3335) used a one-region reactor model, and it is interesting to observe the difference with the present model. The vertical acceleration was varied in the following way. g = go (1 + 0.5 sin 2rft) The reactor was at the same steady state as in the previous section. The coolant flow was again by natural convection. Figure 19 snows the behavior of the reactor when f = 0.5 From these curves, the following observations are made: 1. Unlike the reactivity oscillation test result, the amplitudes of the flux variation are not identical for different core regions. On the other hand, the phases are almost identical. 2. Unlike the case where the reactivity is perturbed, the coolant of both the boiling and the non-boiling regions moved in the same direction. This is understandable since the perturbed quantity is a body force. 3. The void Fraction variation shows a similar characteristics as in the previous section where the reactivity is perturbed; the variation is large near the boiling boundary. 4. The phase relations of the observed quantities are in general agreement with the previous investigation.(35) Although the space-dependence of the flux variation is observed, the effect of the space-dependent variation of the flux on the entire system behavior is not known from these curves. In order to investigate the space-dependent effect, another simulation calculation was compared.

-831.5 g/go I 0.9 WATER VELOCITY \ < SEC 0.100 " VOID FRACTION SEC 0.050 0.025 i=3 1.0 0 1 2 3 4 SEC Figure 19. Response of the Reactor to a in Vertical Acceleration, g=g(+5 sin t) 0.075 - o.025=2 0.75 ^^ 1.25

-84In the latter run, (1) the initial steady state neutron flux shape was maintained throughout the transient, and (2) the void-reactivity coefficient of the steady state flux distribution was used. In Figure 20, results of the two simulations are plotted simultaneously. In most cases it is not possible to distinguish the two curves. The agreement of the two simulations is attributable to the relatively small feedback effect, as well as to the fact that the external variation is a body force, and hence acts evenly on the entire system. Also plotted in Figure 20 are the output quantities obtained by the one-region model. The poor agreement is noticed. The comparison study was performed for other frequencies, i.e., f = 1.0 and 2.0 In all cases the agreement between the five-region adiabatic model and the five-region constant flux shape model is excellent, whereas the oneregion model yielded quite different results from the other two. The ratios of the amplitudes obtained by the one-region model to those by the five-region adiabatic model are shown in Figure 21. E. Discussion of Results Since only a few simulation calculations were'performed, it is not possible to draw a general conclusion as to the spatial effect on the boiling water reactor kinetics. However, the following can be said within the scope of the performed simulations: 1. When the external disturbances are large, the spacedependent effect due to the coupling between the neutron dynamics and the hydraulics becomes important. The simulation of the reactor under the step change in reactivity studied in this chapter is an example.

-85U) u) N, N 0 E-j H 0JI 0 I 0 0 0 0-0 <~ 0 0J ~ Q O C Nj _j r- in I-E- I E- HI I E-I00"-1 o \ y - CO CH cli 104-X~ QSl~~~~~~ ha 0 0 o - - - 5 0 >4 w w w' ~r~p;,~, u, cn u~h N~c~H 0) 0r 049 0r 0l M to H~CC 0~ 0 H 00)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~0)e 00n- c Ic 00F1 II /~ H~II IA 9 10 - - - 0~~

-86Iso~~~~~~~~~~~ S~ 160 C) o 140 \NONBOILING H LENGTH I 120o 100 - _j V- \ ^ — 0w in 80U_ fout o 60 - 0 v w~~~~ out 4020 0.5 I 2 f (cps) Figure 21. Ratio of Aplitude Variations Obtained by One-Region Model to Space-Dependent Model. g = g0(1+0.5 sin 2iftt).

-872. For small perturbations, a linearized model with constant flux shape may be sufficient, provided extreme care is taken in the feedback presentation. A few points where special attention is required are: (1) an exact steady state flux distribution must be provided, (2) the reactivity due to void must be evaluated with proper weighting processes, 3. The variations in vertical acceleration encountered at sea are of considerable magnitude and hence it becomes necessary to treat the system as non-linear. The space-dependence of the neutron flux variation is still negligible for the three cases studied in this chapter, 4. Compared to the simpler models the adiabatic model reveals more of the detailed characters of the reactor behavior. In some cases, however, the space-deperdent characteristic of the reactor is not important and the model with constant flux shape may serve as well. One (8) typical example of such a model is the one developed by Solberg.() In Table III a comparison is made among some representative models.

-88ed _ 0 Cd Cd - ho 4 *H *. H L 0^ 9% 0 40 (, P (^C CT4 t r Q H () 0a 0 C. (U O O 4 QQ)) uCQ' _ oI; ~ l at I ~I H c 04 (U l v O 4 1 H O q O qS I C)C>.d, I H IH 4I HOHy~d {~HtHY u OH H Io OH I0 4M ) CH 4 r i w 1 (1 p 0 4 0 4I E-ii 4 4p q I d0 I; ~ -1.~,' I I ~ O Rr0~+~ 0Q)H-HTQI CdQ)CdC)^c E, CdI 0 0 h-p 0 I 0.4 rQ 4- PA E @ OH I C) V 3' d Q * H Q, 4 * - o H -. — _ __ -—. —.-I..-.-.g.-. ^ 0 ~H 4; 14| S Co C |.,iC | ri ) o I I 2~ 4) HC.0) ~-0 d Cd a) 0O a0 0^ ( l 4-> l OH

CHAPTER VI CONCLUDING REMARKS A. Summary and Conclusion The conventional adiabatic approximation was successfully applied to the analysis of axially space-dependent boiling water reactor kinetics. In order to describe the space-dependent behavior, the reactor core was divided into a number of regions. The equations which describe the behavior of the steam void in each of the boiling regions have been derived. These equations were coupled with the nuclear equations through the adiabatic approximation. To solve these equations numerically a computer program was written and the space-time behavior of a boiling water reactor was investigated. 1. Feedback Dynamics and Formulation of the Adiabatic Model The reactor core was sectioned into axial regions and expressions for the void fraction and the water velocity for each region were obtained from the balance equations of energy and mass of two-phase flow. The void distribution in each boiling region was assumed trapezoidal or linear in between the lower and the upper boundaries of each region. With a similar assumption for the temperature distribution, an equation was derived for the temperature of the coolant in the non-boiling region. When the boiling boundary shifts from one region to another, numerical instability is possible. The unstable situation was avoided by using an extrapolation technique. As illustrated in Figure 13, the boiling boundary behaved identically in models with sufficiently large numbers -89

-90of regions. This means that the extrapolation was properly performed in determining the void fraction and the water temperature near the boiling boundary. A conventional method for the adiabatic approximation adopted by Henry and Curlee(20'2) for the analyses of space-dependent non-boiling water reactors was used to formulate a computer program. Since this technique is essentially a nodal approach to the solution of space-dependent reactor kinetic problems, the reactor core was divided into a number of regions. The number of space-nodes was determined as a compromise between the accuracy of information and the computer time. For the sample reactor studied, a five-region model yielded a satisfactory presentation of the nuclear and the feedback dynamics. Only a part of the storage locations of the IBM-7090 computer was sufficient for the five-region, onedimensional model. If two or three dimensional problems were intended, the computer capacity would limit the reactor space nodes. The regression program solved the difficult problem of taking the correlations of the flux shape and the reactivity to the void distribution. If it is required to have a more complex core division than the present model, it may become necessary to use a more advanced regression program, for example, a program which chooses the best possible transformations for the independent variables. ( 2. Space-Time Behavior of a Boiling Water Reactor In Chapter V the space-time behavior of the sample reactor was observed under three different types of external perturbations. They are (1) a step function change in reactivity, (2) sinusoidal changes in

-91reactivity and (3) sinusoidal changes in vertical acceleration. When the reactivity was added in steps of 90 cents to the sample boiling water reactor, a fairly large space-dependent excursion in the reactor variables was observed. The transient died out within a few seconds and a new steady state was reached. The sinusoidal reactivity was introduced into the reactor with an amplitude of ten cents. Negligible space-dependence was observed, This is mainly due to the fact that the shape of the void distribution was maintained almost constant for the observed frequency range (0.1-2.0 cps). For frequencies higher than those studied here, the shape of the void distribution may change during the oscillation, since the void transit time becomes longer compared with the period of oscillation. The relatively longer heat transfer time constant, however, will decrease the magnitude of the feedback effect. The observed space-independence agrees with experimental results of the Dresden boiling water reactor. 2) The amplitude of the vertical acceleration variation given to the reactor was 0.5 g. Small flux shape distortions were observed but the overall behavior of the reactor was not affected by the small distortions in the reactor power. The small space-dependent effect in spite of the fairly large perturbations is explained by the fact that the perturbed quantity is a body force acting uniformly on the entire coolant flow. 3. The Importance of Space-Dependence in Choosing Analytical Models The present investigation suggests the types of analytical models to be used in specific circumstances. In short, when the external

-92disturbance is large as in the case of large reactivity insertion, the space-dependent behavior of reactors becomes important and the adiabatic approximation provides an economical method of analyzing the situation. On the other hand, the space-dependence of neutron flux is less important when the external disturbances are small as in the case of simulating the reactivity oscillation test. In the latter case, however, the feedback characteristics are space-dependent, thus the constant flux-shape model may safely be utilized only if (1) the feedback model is space-dependent, (2) flux shape is correctly evaluated, and (3) the void-reactivity relation is obtained by a proper weighting process over the space. B. Future Work The present work dealt only with one-dimensional space problems in the axial direction. It is proposed that the work be extended to the radially space-dependent problems and further to the multidimensional problems. In particular, the radial behavior of a boiling water reactor is still to be investigated when used for marine propulsion, since for large reactors rolling in the sea has considerable effects on the behavior of the reactor and hence on its stability problems. The present investigation shows that the adiabatic approximation is a good method to handle such problems. The investigation of Chapter IV shows that the space-dependent feedback characteristics are more sensitive than the nuclear characteristics to the number of regions used. Better methods of presenting the feedback dynamics are worth studying. A better approximation may be found by assigning more nodes t tthe feedback than to the nuclear presentation.

-93As mentioned in Chapter I, the adiabatic approximation as applied in the present study is limited to the cases when the size of the reactor core is less than approximately 20 migration lengths. This limitation is due to the fact that the time rates of change of the delayed neutron precursors are neglected. It may be removed by modifying the adiabatic method as Henry and Willey indicated for the non-boiling water reactors. It is felt that more effort should be made in finding the methods of investigating the space-time problems of nuclear reactors. For instance, it may be possible to develop an adiabatic approximation in a "modal" sense. Another possible approach is a hybrid of the "nodal" and "modal" methods. The space-time synthesis seems to be a promising technique for solving three-dimensional problems. These and other new techniques need to be further explored before applying them to practical problems.

APPENDICES -94

APPENDIX A DERIVATION OF REACTOR KINETIC EQUATIONS FROM TRANSPORT THEORY Henry derived the reator kinetic equations from the transport theory. (5) The resulting time-dependent equations are exactly of the (18) same form as the classical kinetic equations. The time-dependent Boltzmann equation for the neutron flux density at time t, at position r, of neutrons with lethargy u, going into direction 2 is / a3 9(r n2 u t)= -[. ( r- t)it) v(u) at 21 (r 1 j (M12JSclnu> + ( r- l t Y-, AI JAI d - J 0... Q A, C(r) (u) + S(r_,lL:j K, —:)A wh ere.;:' (u.l i;' —; n - i.;d dQ = scattering frequency which is defined as the probability that a neutron having a lethargy u' and moving in a direction S2, will be scattered into du about u and dO about Q -95

-96S(r,u,_;t) = external neutron source ci(r;t) = delayed neutron precursor density of the i-th group gdi(u) = delayed neutron spectrum of the i-th group gf(u) = prompt fission neutron spectrum v(u) = speed of neutrons at lethargy u a, Ef, EZ = macroscopic absorption, fission, scattering cross section. = -delayed neutron fraction of the i-th group P = z Pi ki = delayed neutron precursor decay constant for the i-th group v(u) = expected number of neutrons born per fission by neutrons of lethargy u /v =eigenvalue, that can be so chosen as to keep the reactor critical. The equation for the delayed neutron precursor density is:c i ( _r,t) 1 | "(,C- L' i-t- 2(bl.l);EOr (F, t ri U t) a:t Jo ~- A ~Cj~(r;t) ~(A-2) In Equation (A-l) it is assumed that the reactor consists of only one scatterer and one fissionable material. To simplify these equations, a source free, steady state adjoint equation will be introduced.

-97o = QC71g( r,( u)_t[ZLa0(ft, I to(r L)] (d2,Qu) Ja$ ~ ~ A0 " rU ( -~ - vo ifw fo,(A-3 where 3t =,('-) + Z o P3 and the subscript zero indicates the steady state quantities. Let (r,,LaLkT, ) = (r L,,L;t) N (t) (A-4) Equation (A-4) does not imply the separability of flux into time- and space-dependent parts, and is a rigorous identity. The reduced fluxes are defined as (r,L', ) (,, t (A-5) (rLA) = | * (-r, ) A (A-6) J -,

-98Let d j Il U t t, ) du I r = ~ (A-7) avt u Lt and compute (ru ) < l ) i> p - N(t) (P;(f,t < 3A dZ td dI3r ti' r' A JlQ J where <1> and <3> means Equations (A-l) and (A-3). Some additional algebraic computations yield, dNJl) -P(__- ) L N(t) + A (t) (A-8) t A(t)tC C where F(t) I / ) I2(u/) ( ru )(, I t) r 9Lidu A (tJ-= F ((./iV v(u) 3 3=(if)

-99FM,2 (ru t)d= -- $-l * (-, r,_;t) r) u F(t)}3l 1 ^ I rrb, /i / (upu) (ru.t't) { Z (rt (ru';dut dL( r 01 ct = - LA q(u) =NLO) r J(i U t ) < 2'r U _ _) Q. In order to obtain the corresponding equation for the delayed neutron precursors, compute

-100which yields d C6it) = t3(t) 14(t) -; C; (t) (A-9) (A-9) Equations (A-8) and (A-9) have the same forms as the usual reactor kinetic equations. N(t), the time dependent part of the neutron flux, is called the time function and 0(r,_,u;t) or 0(r,u;t) the shape function. The latter is normalized in the sense of Equation (A-7).

APPENDIX B VESSEL PRESSURE DYNAMICS (24) In Brown's analysis4) of reactor vessel pressure dynamics, the reactor is divided into four parts as shown in Figure 22. They are: (1) subcooled water phase, (2) saturated water phase, (3) interphase, and (4) vapor phase. In the following, each of the four parts are explained. For each phase, balance equations for mass, volume and energy are given. From these equations the time rate of change of the system pressure is derived. Subcooled Water Phase Water is fed to this phase at a temperature below the saturation point. For simplicity, its mass is considered constant. After it is heated to the saturation point, the water moves into the saturated water phase. Saturated Water Phase Since it is assumed that there is no local change in the specific enthalpy of the saturated water, the net energy variation in the saturated water phase is zero. Inter-Phase This is a volume-less phase in which the phase change takes place. As the saturated water is heated, it changes into steam phase, forming bubbles in the coolant channel. If the system pressure changes, phase change occurs between the water and the steam phases (flashing and condensation). -101

-102SATURATED STEAM, M mII SATURATED STEAM s Oo o o^O o 0 0TURATED oo SATURATED 0 ~ 0 WATER, Mw n 0~ O / -II SATURATED WATER REACTOR CORE SUBCOOLED WATER I SUBCOOLED WATER Figure 22. Reactor Vessel Pressure Dynamics.

-103-. Va-apoi ~Phasel The saturated steam phase is produced from the inter-phase through addition of heat and fed from outside the system. Actually, there is no feed of steam from outside the vessel, but steam flows. out of the vessel. For simplicity, the sign of the feed-steam is taken as Figure 22 shows. Usually mIIn is negative. 1. Mass Balance Equations In the subcooled water phase, the mass is considered to be constant. Therefore 0l - - -an (B-l) In the saturated water phase, it is 1,7- W (n + mil (B-2) In the inter-phase, there is no mass. Hence 0' i., s. + Wa (B-3) The time reate of cha-,,ge in the vapor phase is svo =i- + a2t n (B-4) Thee surn of Equations (B-l), (B-2), (B-3) and (B-4) gives the entire sy,-tem mass balance. i-C +M P + YAl r (B-5).. 4-ivT +

-10io42. Volume Balance Equations The total volume of the reactor vessel, V, is constant, and it is assumed that the subcooled water phase also has a constant volume. Therefore, only the saturated water phase and the vapor phase are considered for the volume balance. In the saturated water phase, MW MW M a1v (s:) - _ _ __ __'p (B-6) \ W fVj (0L 2 p and in the vapor phase, ( M) M^ Ms r^)^ — ^ -t" (B-7) The system volume balance is the sum of Equations (B-6) and (B-7). \i O i/ M \ - / M \ V' 0 = -^ J^ (B-8) The mass and the volume balances are combined to give the pressure change in terms of steam and water properties and mass flow rates and stored mass. The volume balance equations give: if W n=V v, / 5 4- ^ i PS \ + + M M _ PS PW VV,

-105By using the mass balance equations, o=- (~li" -+ M, ~ ~G' 2 s p e p )rn ( s + rl, ) ( K - ) + ( I+ i +8 E) (B-9) es Pw Pr In Equation (B-9) hws is a dependent variable which can be eliminated by using the energy balance of the system. 3. Energy Balance Equations The rate of change of the total energy of a thermodynamic system is composed of the rates of change of the internal energy and the work done by the system, if there is any volume change, as given (45) by the first law of thermodynamics. In the subcooled water phase, X-I -= t ^ ( - h I v ) D -p P (B-o0) In the satarated water;ph- ase, ( M e ) M h -- j { ma{ I l A,; ) The left hand side can be rewxitten in terms of h, (M vX;- -+OWv — = h,'^ t ii m^ I;,J) (B-11) Pw

-io6The definition of enthalpy is used in deriving Equation (B-ll) from the previous equation. If the mass balance equation for the saturated phase is used, 0 = M/ ( p ) P (B-12) In the inter-phase, the energy balance equation is simply, Qw = - 1 W (B-13) qw2 s'm ws The energy balance for the vapor phase is O Jh), (B-14) The energy balance of the system is 7L,- +Q+ -V}P {{ ap -+lT,(^-) M+ MwS, ha (B-15) where Pi"P, P, P-s Now the dependent variable ms can be eliminated from Equations (B-9) and (B-15). The equation for the time rate of change of the system pressure becomes

-107(, ( pW ) t m + f) |j^+ ( m Mx) rp =-'-^ S "e, -— )( —-+ —-v- (B-16) his (, e, )( s~~ + *^ mop a V) PsZ 3| + P.2 aP In numerical evaluation of Equation (B-16), a conversion factor is needed to change the mechanical unit to the thermal unit. By using the common notation J for this conversion factor, V is replaced by V/J.

APPENDIX C SAMPLE BOILING WATER REACTOR The model of the reactor core investigated as a sample reactor is similar to the Dresden boiling water reactor,(3) which is composed of Zircaloy-clad uranium dioxide fuel-rods and a light water coolant-moderator. The U235 enrichment of the model reactor is 1.5 per cent, the same as the Dresden's first loading. The operating temperature of the reactor is 286~C (546~F), which corresponds approximately to 1000 psig at saturation. The water to fuel volume ratio is 2.2. Zirconium is the main cladding and structural material of the Dresden core. In the model reactor it is assumed the zirconium to fuel volume ratio is 0.5. In the calculations that follow, the moderator is assumed as a homogeneous mixture of water and zirconium. The fuel rod is 0.6 cm in radius, clad with 0.75 mm thick zirconium. The effective radius of the unit cell is 1.1541 cm. The composition of the core at the operating temperature is shown in Table IV. No steam void is assumed in these figures. Many investigators agree that the twoenergy group theory is sufficient for the boiling water reactor kinetics studies. Thus in the present study, a two energy group model is used in calculating the core constants. 1. Thermal Group Constants Thermal Cross Sections Since uranium 235 is a strong thermal neutron absorber, the thermal neutron spectrum is shifted from the Maxwell-Boltzmann distribution -108

-109TABLE IV COMPOSITION OF THE SAMPLE REACTOR CORE Material Density (g/cm3) Volume Ratio (%) Moles/cm3 U02 10.3 27.03 0.0382 H20 0.762 59.46 0.0412 Zr 6.5 13.51 0.0715 TABLE V CROSS SECTIONS FOR WIGNER-WILKINS SPECTRA Cross Void Fraction (%) Section (barns) 0 10 20 30 40 50 60 (1),,(a) ( v (a) 0.569 0.562 0.555 0.547 0.535 0.520 0.500 H sC 34.58 34.34 34,05 33.65 33.30 32.78 32,10 C25,(b) 361 356 350 343 336 326 313 a af 305 301 297 291 285 276 265 (a) Subsripts a, s, f indicate absorption, scattering and fission, respectively. (b) Superscript 25 refers to U235 The first digit is the atomic number minus 90, and the second digit is the last digit of the mass number.

-110towards higher neutron energy, or flux hardening occurs. Amster calculated the cross sections of various elements averaged over the (46) Wigner-Wilkins spectra under the presence of strong absorbers. These cross sections are interpolated for the operating temperature and for various void fractions. Table V shows the cross sections for H, U235 and a unit 1/v absorber. Macroscopic Absorption Cross Section The absorption cross section of the heterogeneous core is the volume and flux weighted average of the cross sections of the fuel and moderator. Letting E = absorption cross section of the entire core ~U = absorption cross section of the fuel aU saM = absorption cross section of the moderator VU = volume of the fuel V = volume of the moderator M = Ou= averaged thermal flux in the fuel M = averaged thermal flux in the moderator, then 7 = 21M M PM _ + ^VjEau f (C-l) VIAc +M VU qK Now since -U = (thermal disadvantage factor)-1 BM VM ZaM f VU 2aM 1-f

-111where f = thermal utilization Equation (C-l) becomes 1~ ziM (C-2) I - f + L + The calculation of the thermal utilization factor is described in a later section, By using the result obtained for f and the cross. sections of Table V, Ea is calculated. Figure 23 shows the absorption cross section for various void contents in water, Diffusion Area and Diffusion Coefficient The diffusion area L2 and the diffusion coefficient may be defined by the following equations. L1 = = (c-3) = D = (c-4) 3 1tr where tr is the transport cross section and is given by Ltr s ( r- o) (c-5) 1o = average of the cosine of the scattering angle in the center of mass system. With these definitions, the problem of obtaining these parameters reduces to that of cross sections,

-112I I I I 0.050 -0.4 -15 0.048 E E~ E 0. 046 - -0.3 -10 0.044 - L2 0.042 0.040 -0.2 -5 0. 038 0 10 20 30 40 50 60 VOID FRACTION, % Figure 23. Slow Group Constants.

-113Characteristic areas of media with holes differ from those of the media without holes due not only to the density difference but also to the streaming effect. This effect has been studied extensively by Behrens.(4748 49) For a medium whose mean free path is X and contains holes of hydraulic radius r, the characteristic area L2 is related to L, the characteristic area without holes, by the 0 following equation. I LA I \ (0 C 0 -F =cX r 2ogt (C-6) where a is the volume fraction of the holes. When the streaming effect is included, the transport cross section of boiling water for neutrons of energy E can be expressed as follows: tr \bo\\i9} waie ie _ V } bo/ii<:/ wct^er * ( i E ) = J; I * -r) bo; I,'n = 2 oiling (c-7) where?L (Eo)J is calculated by applying the Behrens' corrections to the known value of L (E). The diffusion coefficient of the fuel-moderator mixture is j) ~_D - -- __...._I (c-8) 3 ( Vovtui ewew'r i 2!11* tr ).

-114and the diffusion area is obtained from Equation (C-3), where in turn D and a are obtained by Equations (C-8) and (C-2). Figure 23 shows these slow group parameters. The radii of bubbles are assumed to be one half of the coolant channel thickness when the void fraction is unity and to decrease linearly with the void content. 2. Fast Group Constants Effective Fast Group Cross Sections The fast group constants can be calculated in a similar way to the slow group constants, if the effective cross sections are known for the fast group. Hellens obtained the effective cross sections to be used for this purpose. He fitted the two group age, fast removal cross section, and fast diffusion constant obtained by the MUFT II code* with the following "age diffusion" expressions. Df = 1 Ni (c-9) T: = \ 3 L Zt L Ctr M(} 5 )) (C-ll) L j where Ni is the atomic density of the i-th nuclei, and the fitted cross sections are shown in Table VI. * The MUFT code computes the energy distribution of neutrons having a given Fourier mode in an infinite medium. The output includes the fast constants for a variety of few group schemes.

-115TABLE VI EFFECTIVE FAST GROUP CROSS SECTIONS* Nuclei atr (barns) as (barns) Hydrogen 2.67 0.711 Oxygen 3.00 0.060 Uranium 9.10 0.264 Zirconium 4.00 0.000 * (Taken from Reference 50) Slowing Down Area The neutron Age or the slowing down area is calculated by Equation (C-ll), by using effective cross sections of Table VI. The streaming correction, though small, can also be made to the age. Figure 24 shows the corrected age. Fast Removal Cross Section The fast removal cross section is calculated by Equation (C-10). Figure 24 shows the calculated cross section as a function of void fraction. Fast Diffusion Coefficient The relationship between Dfast ft and T is D = By using Efast of Equation (C-10) and the stream-corrected T Dfat is calculated. Figumre 24 shows the fast diffusion coefficient.

-1160.025- y /-^Q -1150 fast _ 0.020 - 2.0 o 1LI I I. I 100 0.015 - -1.5 50 0* 10 20 30 40 50 60 VOID FRACTION (%) Figure 24. Fast Group Constants.

-1173. Calculation of Four Factors Fast Fission Effect, e For a homogeneous mixture of fuel and moderator, Hellens made a one-velocity fit for the fast fission factor using the result of (50) MUFT II calculation. 6-: = 2I-f- —'- (C-12) hnel C where 2,, c f' inel are macroscopic capture, fission and inelastic cross sections respectively. The fuel to water volume ratio was in the range of 0.25-1.0 and the constants v = 2.48, 08 = 0.14 were assumed. The microscopic cross sections for the fast effect are tabulated in Table VII. Fleishman and Soodak give the cross sections based on a threegroup model foar the fast fission factor evaluation of a cluster of fuel elements.(5l) In Table VII their values for Zr cross sections are listed. The fast fission factor for the homogeneous mixture of water, UO2 and Zr are calculated by Equation (C-12) and the cross section of Table VII. TABLE VII CROSS SECTIONS FOR FAST FISSION FACTOR CALCULATION Nuclei f(bar ) c(barns) ain 1(barns) Ref. Nuclei far ) ac cinel Hydrogen.00 0.00 1.38 50 Oxygen 0.00 0.024 0.171 50 Uranium 0.30 0.042 2.32 50 Zirconium 0.00 0.029 0.960 51..... J...,...

-118In heterogeneous reactors, the fast fission phenomenon is induced by the fast neutrons entering into the neighboring fuel elements without collision as well as by those that cause fast fission in the fuel element where they are born. Clark evaluated this effect for slab arrays and obtained the following expression for the collision probabil(ty52) ity. - U — bb _ TS b u Po(TSZb) = - a (- )(i — d (c-13) b I - e- b(l+TS)u a3 where b = thickness of fuel slab S = ratio of the thickness of moderator slab to that tM ku of fuel slab = t E tU ZM T = (tot) (tot)U = Etot = total macroscopic cross section for fast neutrons bt = width of the fuel slab in mean free paths = bU tM = width of the moderator slab in mean free paths If the medium is homogeneous PO = Po(TS,O) The.fast fission factor is considered to be dependent on the collision probability PO 0 For heterogeneous systems, the fast fission factor is obtained by the following relationship P0 (TS,Eb) ( 1)heterogeneous, fuel = )homogeneous p (TS,O) slab thickness b (c-14)

-119In media other than slab geometry, the heterogeneity effect can be estimated by using Equation (C-13) for a slab lattice of the same water to fuel ratio as the given lattice. Figure 25 shows the fast fission effect with this correction. Resonance Escape Probability, p The resonance escape probability is calculated in the usual way by Equation (C-15), p= _xp 1; M a 0) (,f) } (c-15) ( 6 ^MVM \s1 4 re S where NU = number of atoms of fuel tFe = slowing down power in the moderator = Au aM Au = resonance lethargy bandwidth = 3.0(53) f.(53 ) aM = o.46 cmi for H20 of p = 1.0 = 0.93 for H20 of p = 1.0 t[O/'M]res = resonance flux disadvantage factor Au [RI]eff = effective resonance integral = EU eff aU Nu The resonance flux disadvantage factor is, for the array of infinite cylinders, P = F'+ (E/-) va (C-16) ~}l) " ^/ EM

-120where F' = X1 r Io (x ru) 2 I C (X rv) E = CM(rM V ) II (crM) K<(MrU)+ IO(Mr,) K (aKd), r, I, ( K rt1)<o (YX ru)+ I | (r, nv) K, (XH i ) rUrM = radii of fuel rod and unit cell eM = 0.885 cm1 for water, p = 1.053 KU = 0.31 cm-1 for U2 p = 10.03) The effective resonance integral for water-UO2 array is given by Equation (C-17). (4) RI~( =RI ~00^) t 0^ t o,s8t+0, M )( - 0) (C-17) where U02 s1/2 s RI (G0) = 2.81 + 28.7 ( ), for 0.07 < M < 0.53 S/M = surface-mass ratio of the fuel rod = fuel temperature, ~K Q0 = 2930K The surface area of the fuel element must be corrected for the shadowing effect when the fuel array is tightly packed. Using (55) Bell's approximation, S S -u (c-18) corrected SU 1 + 4VM

-121where M = moderator macroscopic cross section. Table VIII shows the calculated resonance integral. The resonance escape probability is then p = exp ( - ) (C-19) where VM (I')F F', VU Nu (RI)e ff Figure 25 shows the resonance escape probability. TABLE VIII RESONANCE INTEGRAL AND RESONANCE ESCAPE PROBABILITY Void Frac- 0 10 20 30 40 50 60 tion (%) RI 16.144 15.977 15.779 15.543 15.256 14.897 14.437 p 0.8859 0.8752 0.8623 0.8463 0.8260 0.7992 0.7621 Regeneration Factor, r The regeneration factor is evaluated simply by au (c-2o) It is seen that r' changes very little with the void content. By taking v = 2.43, r = 1.604. (C-21)

-122r I''''I I I I -4 — rq (=1.604) 1.60 j -1.30 1.0 -1.20 0.90- -^ — ^^ ^^^ -1.10 0.80 \- 1.0 I I I I I, X 0 10 20 30 40 50 60 VOID FRACTION, % Figure 25. Four Factors and kc.

-123Thermal Utilization, f The thermal utilization factor, i.e., the ratio of thermal neutrons absorbed in the fuel element to all the thermal neutrons absorbed in the fuel and other materials, is calculated in the usual fashion(56) by using the diffusion approximation and is shown in Figure 25.

APPENDIX D SIMPLIFIED DIFFUSION CODE (SDC) The nodal approximation to the diffusion equation(57) is developed and programmed for the computer calculations. It is assumed that 1) neutron energy is one group, 2) the reactor is bare, 3) material is homogeneous in one region, 4) the reactor is one-dimensional slab reactor with infinite length, 5) in one region, flux changes linearly from the center to the boundary, 6) at the interface, neutron flux and the current are continuous, and 7) the flux vanishes at the outer boundaries. The diffusion equation is, then, Ad r7 U ) -, -\7+- --- + (D-1) where D = diffusion coefficient Z = absorption cross-section k = number of neutrons produced per absorption p = eigenvalue Figure 26 shows the slab reactor. -124

-1250 ",| Il <d" — ^ -+Z H H C P^ r'Oi< C\j + 0 0 +.-t...... H M-f r t _- t+ gt > — 4 0o H ~ H <:,-... p4 \F4 rr\J <] ^\r i, I^ \ 0

-126In the v-th region, Equation (D-l) becomes D 7V2 P- + kE = O (D-2) If Equation (D-2) is integrated in the v-th region, the first term becomes i JJ Dv V dV = D j Vn P d (D-3) w-th Regjon R-th Res on where V is taken normal to the surface. The second and the third n terms become, - fff ^dV = - f fvc v (D-4) frphkEdAd V =A k>3iv (D-5) The neutron current continuity at xv+1 requires 7DyV CV ---- v P 4tv l (D-6) In the right half of the v-th region, v<,Pt^ =V (D-7) 2 then the continuity condition becomes D AD _-__ -I_ (D-8) Z 2

-127Eliminating O~ v+l from Equations (D-7) and (D-8), 2DV _ = _ A. __ - ( V+I- vD) - I a), + ( )v Thus, in the v-th region, D Vn dS = (AJ (- )( + v-th Iv AV+( Region XDV-l (D P ( + 2t /^\v-/^ (( v-l *C/)v) (D-10) And the integration of Equation (D-2) reduces to Equation (D-ll). 2 (D) ~-,) ( i -k v i D\( ^-.1 (W L y-D 2 (s)( i ) + 2(') (ATh+l) -_.) v } - + l I. A I Ki,+ z D A i) ( <) ) = Q (D-ll) A^^~ (1:)')} [V l/ \l / +^

-128At the boundaries of the slab, neutron leakagaes are calculated from the current leaking out of the outer surfaces. At the Right Boundary J_ + - 0 (D-12) where AN =3A d \ s _ - C R, dx R AN Equation (D-12) becomes " 4 _ N D X (D-13) then (p, _ dp __ 2pDt Jt 4 6 (<)UI + + 4D, (D-14) Thus the integrated leakage term becomes J J( (IDN =A A M) (I N_ ANI I INN 2ADN (D-15) Att4Di N

-129and the integrated diffusion equation takes the following form in the N-th region. 4 4 _ _, 2 D- (DA- i k D - -+ -=? — __6) l^-^^^^t~i' ANl (A +4Dw)E kH kN At the Left Boundary J+ Id 4 x, Use of (D-18) in J- gives, PL 6 (, A, -D, - (D-17) 2 Therefore - ^ - Di 4, (D-18) Use of (D-18) in J- gives, T -- /^^i^} /,-t-&D (D-19) J- if- k[ d>! 1^ Al ^D

— 130and the integrated leakage becomes DVrn S 2 A- D2 + I D (I,Z ( 2o2 ) (D-20) 1;t Al AZ Reg on The diffusion equation then yields to -L^ + +1AA-J -'"(A D,) + )a),, +, kJ (n~~,+4lc/D,)k,D,4, i,),2 kf, k If it is defined that aMi* = 2 Di D,-' A — -:M -,,,D - 2., (J +4Dt)g, A,1 --— ~ -^-2D ki, (I +- -D, 2 A i - M -l, tM., = 2t3,,INA' ^ ^ ^ -

-131At~ - _,. - 2 D NA = - - -^ t- 2 N'NN k_ 4_ -__ __, - Ks __+ - (i +L9N)EN _M, V-I =2 A - _AA 0, then the original diffusion equation rlisteds to the following set of simultaneous equations. [A, A2 0 0... 0 0, 0 A32 A33 A, 0 | = 2 (D-22) o 0 0 ANN- NN JN) in a similar way as FOG. The program is listed on the following pages.

-132$COMPILE MAD, EXECUTE, PRINT OBJECT, PUNCH OBJECTs DUMP R SIMPLE DIFFUSION CODE D'N A(240,DIM)g B(240,DIM),D(15),DIM(2)tG(15),LAMBDA(15),M(15 1),S(15),SIG(15),NSIGF(15),FLUX1(15),FLUX2(15)9S1(15)~S2(15) INTEGER N,I,COUNT,CLOCK1,CLOCK2,TIMENNITERNONO1,N02,MAX INTEGER DAYTIM. EQUIVALENCE (DIM(2),NN) V'S DIM=2,1 R READ IN Lo Dl, DN, SIG1, SIGN. NSIGF19 NSIGFN RIA R READ IN NYPSILNoTHETA,MAX ANN R'A NN=N+1 R CALC DEL DEL=L/N P'S DEL R CALC D AND SIGMA T'H A1,FOR I=1,1GI.G.N D(I)=(DN*(I-.5)+Dl*(N-I+*5))/N NSIGF(I)=(NSIGFN*(I-.5)+NSIGF1*(N-I+.5))/N Al SIG(I)=(SIGN*(I-.5)+SIG1*(N-I+.5))/N P'S D(1).*.D(N), SIG(1).*.SIG(N), NSIGF(1)...NSIGF(N) R CALC M T'H A2,FOR I=1,l,I.G.(N-l) A2 M(I)=2.*)D(I)*D(I+1)/(DEL*(D(I)+D(I+1)) M(0)-2.*D(1)/(DEL+4.*D(1)) M(N)=2.*D(N)/(DEL+4.*D(N)) P'S M(O)...M(N) R CALC B T'H A3, FOR I=1I1,IGG((N*(N+1)) A3 B(I)=O. T'H A4, FOR I=191,I.G.N A4 B(II)=SIG(I)*DEL+M(I-1)+M(I) T'H A5, FOR I=11,PI.Ge(N-1) B( II+1)=-M(I) A5 B(I+1I)=-M(I) R SET INITIAL SOURCE T'H A6, FOR I=1II*.G.N A6 B( IN+1)=NSIGF(I) P'S B(1,1).*.B(N9N+1) PRINT COMMENT $1$ R START ITERATION PROCESS ITERNO=1 R RECORD TIME CLOCK1=DAYTIM.(0) COUNT=O R PUT B IN A ARBOR T'H A7, FOR I=191,I.G.(N*(N+1)) A7 A(I)=B(I) R STORE SOURCE T'H A8, FOR I=11,I.G.N A8 S(I)=A(I.N+1) R CALC FLUX X=GJRDT.(NN+1,A( ) DET) W'R X.E.1.O R CALC G AND LAMBDA T'H A9, FOR I=1,i,I.G.N G(I)=DEL*NSIGF(I)*A(I*(N+l)) A9 LAMBDA(I)=G(I)/S(I)

-133R CALC LAMBDA MAX AND MIN LMAX=LAMBDA (1) LMIN=LAMBDA( 1) T'H A10, FOR I1=11,I*G.(N-1) W'R LMAX.L.LAMBDA(I+1) LMAXLAMBDA (I+1) O'R LMIN.G.LAMBDA(I+1) LMIN=LAMBDA( I+1) A10 E'L R CALC EIGENVALUE LAMBDA=O0 T'H All, FOR I=1,l1.G.N All LAMBDA=LAMBDA+G(I) R STORE THE RESULTS IF EITHER CRITERION IS MET AND TERMINR ATE THE ITERATION AND PRINT THE RESULTS IF THE CRITERION R IS MET TWO TIMES IN RUN W'R COUNT.E.l NO1= I TERNO EIGN1=LAMBDA T'H Bl, FOR I=11,9I.G.N FLUX1(I)=A( I N+ ) 61 S1(I)=G(I) O'R COUNT.E.2 N02= I TERNO EIGN2=LAMBDA T'H B2. FOR I=1,lI.G.N FLUX2(I)=A(I N+1) 82 S2(I)=G(I) P'S N01, EIGN1, Sl(1).o. Sl(N), FLUX1(1)*..FLUX1(N) PRINT COMMENT $1$ P'S N02, EIGN2, S2(1)*o.S2(N), FLUX2(1)o..FLUX 2(N) PRINT COMMENT $85 R RECORD TIME CLOCK2=DAYTIM (0) TIME=CLOCK2-CLOCK1 PRINT COMMENT $OELAPSED TIME. IN 1/60 OF A SECOND.....$ P'S TIME PRINT COMMENT $1$ T'O ANN E'L R CHECK CONVERGENCE W'R(.ABS.((LMAX-LMIN)/LMAX).L.YPSILN).OR.(ITERNO.GE.(MAX-1)) COUNT=COUNT+1 O'E COUNT=O E'L R CALC NORMALIZED SOURCE9 G, AND PROCEED TO THE NEXT R ITERATION T'H A12, FOR I=11,,I.G.N A12 G(I)=G(I)/LAMBDA R CALC EXTRAPOLATED SOURCE T'H A13 FOR I=1,1,I.G.N G(I)=G(I)*(l.+THETA)-THETA*S(I) A13 B(I N+1)=G(I) ITERNO=ITERNO+1 T'O ARBOR E'L E'M

APPENDIX E MAIN PROGRAM LISTING $COMPILE MAD, EXECUTE, PRINT OBJECT, PUNCH OBJECT, DUMP R MP D'N Q(5),F(6),SR(6),VW(6)W(6),A(6 DM(1) ),DM(3),WW(6), 1FF(6),DEL( 10),BETA C(6,C(6),N(20) QSUM(6) tTW(6),X(25) 1Y(6C),Z(35),TEMP(35),QQ(5) FRACT(5),QQQ(5),FBAR(5),XBAR(5), 1LAMBDA(6),K(5 ),EFF(5),XC(6),VWBAR(5),B(240,DIM) DIM(3), 1FR(6),FRACTR(5) FRBAR (5) VWR 6),VWRBAR(5),SRR(6),XR(6) D'N IMAGE( 5000),PLOTY(25), NSCALE(5),FBARO(5) QQO(5) VWO(6) V'S Dt(t1)=2,0,3 V'S DIM=3,141,6,2C V' DAT. ( 10F7.5)*$ V'S OPT IN=$16I2*$ V'S NSCALE=,0,3,0,3 V'S LABEL=$ $ INTEGER I,CCUNT,SUREN,EQNO,BRNO, JBRNO1,IIJJIIMAX INTEGER BBNOBBNO1,WATCHCLOCK. INTEGER NSCALEFIG FIG=1 EQUIVALENCE (Y F(1 )),(Y(5),FR( ) ), (Y( 11),Q(1)), (Y(16),P) (Y(1 17),MTOT),(Y(18),TW(1 ) ),(Y(24) NBL),(Y(25),T),(Y(26),C),(Y(33) 1,QG(1)), (Y(38),QQQ(1 ) ) (Y(43),VW(1) ) (Y(49),QEFF(1 ) ) (Y(54) R 1HO) (Y(55),G) EQUIVALENCE (XR,XC(6)),(TSAT,TW(6)),(TINTW(1)),(FRACTRFRACT 1(5)) R AC=AC*10.P.-4 AND POWER=POWER*1O.P.+4 R'T DATA, B(00C0, )...B(1,5, 18) START R'T OPTION, N(1)...N(16) PRINT COMMENT 51$ PRINT CO1MMENT $ MP, AUG 28 $ PRINT COMEFNT S8$ R PRINT OPTIONS W'R N(1).EoO PRINT COMMENT $ CONSTANT SHAPE FUNCTION$ C'R N(1).E.1 PRINT COMMENT $ VARIABLE SHAPE FUNCTION$ E'L W'R N(2).E.O PRINT COMMENT $ ONE DELAYED NEUTRON GROUPE$ O'R N(2).E.1 PRINT COMMENT $ SIX DELAYED NEUTRON GROUPESS E'L W'R N(3).E.O PRINT COMMENT $ NO EXTERNAL REACTIVITY$ O'R N(3).E.1 PRINT COMMENT $ SINUSOIDAL EXTERNAL REACTIVITY$ O'R N(3).E.2 PRINT COMMENT $ RAMP EXTERNAL REACTIVITY$ O'R N(3).E.3 PRINT COMMENT $ STEP EXTERNAL REACTIVITY$ E'L W'R N(4).E.O PRINT COMMENT $ CONSTANT VOID-REACTIVITY COEFFICIENT$ O'R N(4).E.1 PRINT COMMENT $ EXACT VOID-REACTIVITY CORRELATION$ E'L W'R N(5).E.O PRINT COMMENT $ CONSTANT SLIP RATIO$ O'R N(5).E.1 PRINT COMMENT $ VARIABLE SLIP RATIO$ -134

-135E'L W'R N(6).E.O PRINT COMMENT $ CONSTANT VW(1)$ C'R N(6).E.1 PRINT COMMENT $ CONSTANT PUMP HEAD AND RISER HEIGHTS O'R N(6).E.2 PRINT COMMENT $ CONSTANT PUMP HEAD AND VARIABLE RISER HEIGHTS O'R N(6).E.3 PRINT COMMENT $ CONSTANT RISER HEIGHT AND VARIABLE PUMP HEADS E'L W'R N(7).E.O PRINT COMMENT $ CONSTANT TW(1)$ O'R N(7).E.1 PRINT COMMENT $ VARIABLE TW(1)$ E'L W'R N(9).E.O PRINT COMMENT $ CONSTANT BOILING BOUNDARY$ O'R N(9).E.1 PRINT COMMENT $ VARIABLE BOILING BOUNDARY$ E*L W'R N(10).E.O PRINT COMMENT $ CONSTANT PRESSURES O'R N(10).E.1 PRINT COMMENT $ VARIABLE PRESSURES E'L W'R N(11).E.O PRINT COMMENT $ RISER ONE REGIONS O'R N(11).E.1 PRINT COMMENT $ RISER FIVE REGIONS$ E'L W'R N(12).E.O PRINT COMMENT $ CONSTANT G$ O'R N(12).E.1 PRINT COMMENT $ SINUSOIDAL G$ E'L W'R N(13).E.O PRINT COMMENT $ CONSTANT STEPS O'R N(13).E.1 PRINT COMMENT $ VARIABLE STEPS E'L PRINT COMMENT $1$ P'S N(1)...N(16) R'A POWER=POWER*MW/200. PRINT COMMENT $1$ R INITIAL CONDITION CALCULATION SUBROUTINE TO FOLLOW R CALCULATE VW'S AND F'S W'R N(6).GE.1 II=l JJ=O K2=K(2)*RISER/5. K3=K(3)*(RISER+5. )/10. CIRCL CONTINUE E'L T'H CONVFOR I=1,1,I.G.6 R TAKE 1./VW CONV W(I)=1./VW(I) COUNT=O R CALC POWER IN THE CORRECT UNIT

-136BEGIN QP=(Q(1)+Q(2)+Q(3)+Q(4)+Q(5))/POWER T'H CO, FOR I=11,I.G.55 CO Q(I)=Q(I)/QP R CALC STEADY STATE VW'S AND F'S I=BRNO REPEAT W'R I.E.BRNO A( 11)=Q(I)*( l.-DELNB)/ARH+F( I )*SR(I)/ /W I) O'E A( l1 ) =Q(I)/ARH+F(I)*SR(I)/W(I) E'L A( 1 2)=-SR( I +1) A(13)=O. A( 2 1) A(1,1)+VW(1) A( 2 2)1.-SR(I+1) A(2t3)1l. W'R N(15).E.1 P'S A( 1 1)..A(2,3) E'L WW( I+1)W( I+1) FF(I+l)F( I+1) X=GJR. (2 3,AD) W'R X.E.1.0 W'R N(15).E.1 P'S A(1,2) 1./A(1 2),A(2,2),(A(2,2)+F( I))/2. E'L W( I+1)A(1,2) ( I+1)A(2,2) DEL(2*I-1)=.ABS.((W(1+1)-WW(I+1))/W(I+1 ) DEL(2*1 ).ABS. ((F(I+1)-FF(I+1) )/F(I+1)) O'E T'O START E'L 1=1+1 W'R I.LE.5. T'O REPEAT R STEADY STATE VALUES ARE CALCULATED FOR ALL FIVE NODES R ITERATE FOR THE NEW SLIP RATIO'S AND THE NEW SHAPE FUNCT R SLIP SUBROUTINE. CALCULATE SLIP RATIO G=GZRO W'R N(5).NE.O T'H SSLIP, FOR II1, 1,I.G.6 XC(I )RSW*F(I) *SR(I)/(l.-F(I)) SSLIP SR( I )=SLIP. (N(5),G,XC( I)) F'L R SHAPE SUBROUTINE. CORRELATE Q'S TO F'S R CALC FBAR, XBAR, VW, VWBAR VW(6 )1./W( 6 ) T'H SSFX9 FOR I1,1.I.G.5 VW( I )m1 /W(I) W'R I.E.BRNO VWBAR(I)m(VW(I+1)*(1.-DELNB)+VW(I )*(1+DELNB))/2. FBAR(I)FC (I+1*(1.-DELNB)/2. O'E FBAR( I ) (F( I )+F(I+1) )/2. VWBAR( I ) (VW( I)+VW(I+1) )/2 E'L SSFX XBAR( I )RSW*FBAR( I )*SR(I)/( 1-FBAR I)) W'R N(1).NE.O N=1 T'H SSX1. FOR I35,-1,I.L.2

-137T'H SSX1* FOR J=1, 1,J.G.3 X(N)=FBAR(I).P.J SSX1 N=N+1 T'H SSX2, FOR I=5,-1lI.L.3 T'H SSX2, FOR J=2, 1,J.G.3 X(N)=(FBAR(I)-FBAR(I-1)).P.J SSX2 N=N+1 T'H SSHAPE, FOR I=ll1I.G6.5 SSHAPE Q(I)=SHAPE.(I,XB) E'L R TEST CONVERGENCE OF F AND W FOR ALL POINTS W'R N(15).E.1 P'S DEL(2*BRNO-1)...DEL(10) E'L I=2*BRNO-1 C1 W'R DEL(I).LE.YPSILN W'R I.L.10 I=I+1 T'O C1 0'E COUNT=COUNT+1 W'R COUNT.E.1, T'O BEGIN E'L' E COUNT=O T'O BEGIN E'L R STEADY STATE VW'S AND F'S ARE CALCULATED FOR THE INITIAL R VALUES R ITERATE FOR NATURAL CIRCULATION W'R N(6).GE.1 R. CALC MDOT W'R N(6)E*E2 K2=K(2)*RISER/5. <3=K(3)*(RISER+5.)/10, E'L MDOT 1=MDOT MDOT =(VW( 1)P.2.-VW(6).P.2*( (1.-F(6))+RSW*F(6)*SR(6). P 2 )) 1/DELTA+G*RISER*F(6)-K 2 *VW(6).,P.18/(1.-F(6)).P.*2-K 3 *VW(1 1).P.1.8+G*HEAD T'H SSMD, FOR I=1l,,I.G.5 SSMD MDOT =MDOT -K(1)*VWBAR(I).P.1.8/(1.-FBAR(I)).P..2+G*FBAR(I) W'R JJ.NE.O T'O HALF O'R I.NE. 1AND.MDOT*MDOT1.L.O. JJ=1 HALF W'R N(6).E.1 DV=DV/2. O'R N(6).E.2 DR=DR/2. O'R N(6)*E.3 DH=DH/2. E'L E'L W'R MDOT.G.O. W'R N(6).E.1 VW( 1 )=VW( 1 )+DV O'R N(6).E.2 RI SERR I SER-DR

-138O'R N(6).E.3 HEAD=HEAD-DH E'L O'E W'R N(6).E.1 VW(1)=VW(1)-DV O'R N(6).E.2 RISER=RISER+DR O'R N(6).E.3 HEAD=HEAD+DH E'L E'L II=I1+1 W'R N(6).E.1 W'R VW(1).L.VMAX.AND*VW(1).G*VMIN.AND.II.L.IIMAX.AND.DV.G.DV 1MIN W'R N(15).E*1 P'S VW(1),MDOT E'L T'H NEWVW. FOR I=2,1,I.G.6 NEWVW VW(I)=VW(1)/(l.-F(I)) T'O CIRCL E'L O'R N(6).E.2 W'R RISER.L.RMAX.AND.RISER.G.RMIN.AND.II.L.IIMAX.AND.DR.G.DRM 1IN W'R N(15).E.1 P'S RISER, MDOT E'L T'O CIRCL E'L O'R N(6).E.3 W'R HEAD.L.HMAX.AND.HEAD.G.HMIN.AND.II.L.IIMAX.AND.DH.G.DHMIN W'R N(15).E.1 P'S HEAD, MDOT E'L T'O CIRCL E'L E'L W'R N(6).E.2, P'S RISER W'R N(6).E.3 HEADO=HEAD HEAD=HEADO+HEAD1 P'S HEADO, HEAD E'L W'R N(6).E.1 W'R DV.G.DVMIN, T'O FAUL O'R N(6).E.2 W'R DR.G.DRMIN. T'O FAUL O'R N(6).E.3 W'R DH.G.DHMIN, T'O FAUL E'L T'O PASS FAUL PRINT COMMENT $1INCREMENTED VALUE DID NOT CONVERGE$ T'O START PASS CONTINUE R CALC MTOT R MTOT=MTOT/(RHOW*DELTA) MTOT=(5.+RISER)*VW(1)+(1.-F(6))*ACR*VW(6)*RISER

-139MTOT=MTOT+VW(1)*(BRNO-l.+(1.+DELNB)* (5-F(BRNO+1)*(l.-SR(BRNO 1)*RSW)/6.))+VW(BRNO+1)*(l.-bELNB)*(.5-2.*F(BRNO+1)*(1.-SR(BRN 10+1)*RSW)/6.) T'H SSMT, FOR I=BRNO+1,1,IG*.5 SSMT MTOT =MTOT+VW(I)*(.5-(2.*F(I)+F(I+1))*(l.-SR(I)*RSW)/6.)+VW(I 1+1)*(.5-(F(I)+2.*F(I+1))*(l.-SR(I+1)*RSW)/6.) R CALC NEW VW(1) FOR CHECK VW(1)=MTOT-(l.-F(6))*ACR*VW(6)*RISER-VW(BRNO+1)*(1.-DELNB)*(. 15-2.*F(BRNO+I)*(1.-SR(BRNO+1)*RSW)/6.) T'H SSVW, FOR I=BRNO+1,1,I.G.5 SSVW VW(1)=VW(1)-VW(I)*(.5-(2.*F(I)+F(I+1))*(1.-SR(I)*RSW)/6.)-VW( 1I+1)*(.5-(F(I)+2.*F(I+1))*(l.-SR(1+1)*RSW)/6.) VW(1)=VW(1)/(4.+RISER+BRNO+(1.+DELNB)*( 5-F(BRNO+1)*(1.-SR(BR 1NO)*RSW)/6.)) P'S VW(1)..VW(6),MTOT E'L R OTHER INITIAL VALUES AND CONSTANTS W'R N(2).E.O C=BETA*T/(LSTAR*LAMBDA) R CALC C'S FOR N(2)=1. ALSO LAMC O'R N(2).E.1 LAMC=C0 T'H SSYO9 FOR I=11,I.G.6 C(I)=BETA(I)*T/(LSTAR*LAMBDA(I)) SSYO LAMC=LAMC+LAMBDA(I)*C(I) E'L R INITIAL POWER. Q=QQ=QQQ T'H SSY1, FOR I=1,19I.G.5 Q(I)=Q(I)/QP QEFF(I)=Q( I QQ(I)=Q(I) SSY1 QQQ(I)=Q(I) R CALC QSUM QSUM(1)= 0. T'H SSY2, FOR I=2,1,I.G.6 SSY2 QSUM(I)=QSUM(I-1)+Q(I-1) R CALC SUBCOOLING HSUB=(QSUM(BRNO)+DELNB*Q(BRNO))/(AW*VW(1)) P'S HSUB R N(9)=2 CASE IS NOT YET WRITTEN R CALC BOILING BOUNDARY DELNB=1. T'H SSY39 FOR I=1,1,DELNB.LE.1.-RASCH SSY3 DELNB=(HSUB*AW*VW(1)-QSUM(I))/Q(I) BRNO=I-1 W'R DELNB.GE.RASCH DELNB1=DELNB BBNO=BRNO O'E DELNB1=1.+DELNB BBNO=BRNO-1 E'L P'S BRNO, DELNB P'S BBNODELNB1 R CALC TOTAL NONBOILING LENGTH QBR=Q(BRNO)*(1.-DELNB) DELTAB=DELTA*(l1.-DELNB) NBL=(DELNB+BRNO-1.)*DELTA P'S NBL

-140R CALC WATER TEMPERATURE, TW T'H SSTW1, FOR I=5,-1I*L.BBNO SSTW1 TW(I)=TSAT W'R DELNB.G.O., TW(BRNO)=TSAT-Q(BRNO)*DELNB/(AW*VW(1)) T'H SSTW2, FOR I=BRNO-1,-1*I.L.1 SSTW2 TW(! )=TW(I+1 )-Q( I) / (AW*VW( 1) W'R N(7)*E.1, TINO=TIN R CALC RISER VARIABLES W'R N(11).E.1 R=RISEP/5. DELTAR=DELTA*R K2=K2*ACR.P.-1.8/5. T'H SSRR, FOR I=11,I.G.6 FR(I)=F(6) SRR( I )=SR(6) XR(I)=XC(6) SSRR VWR(I)=VW(6)*ACR T'H SSRR1, FOR I:=11,I.G.5 FRBAR(I)=F(6) FRACTR( I )=FRACTR SSRR1 VWRBAR(I )=VWR(1) P'S DELTAR E'L R CHECK THE INITIAL VALUES R CALC RHOZRO T'O SSRV(N(4)) SSRV(0) RHOZRO=-ALPHA*FBA( 1 )+FBAR( 2 )+FBAR( 3 )+FAR ( ( 4 )+FBAR ( 5) ) T'O SSREND SSRV(1) RHOZRO=SHAPE. (0X,B) SSREND CONTINUE R SS COMPUTATIONS ARE COMPLETED SURE=1 T'O S(1) BACK P'S Y(1)...Y(55),Z(1)...Z(25) W'R N(14).E.1, T'O START R PREPARE GRID FOR PLOT W'R N(16)*L.0 EXECUTE PLOT1. (NSCALE,25,10 11,10) EXECUTE PLOT2.(IMAGE, XMAX, XMIN, 25.90.) EXECUTE OMIT.(3.) XDEC=(XMAX-XMIN)/110. T'H BLANK, FOR BLANKX=XMINXDEC,BLANKX.G.XMAX T'H BLANK, FOR BLANKY=12.1*.1,BLANKY.G.12.9 BLANK EXECUTE PLOT3.($ $,BLANKX,BLANKY,1) NBLO=NBL PO=P QQO=T STEPO=STEP VWO (6) =VW(6) T'H INITIA, FOR I=1.1,I.G.5 FBARO( I )=FBAR( I ) QQO(I)=QQ(I) INITIA VWO(I)=VW(I) E'L SURE=0O TIME1=-SCALE WATCH=CLOCK (0) PRINT COMMENT $1$ R SOLVE DIFFERENTIAL EQUATIONS

-141EQNO=26+6*N(2) TIME=O. EXECUTE SETRKD.(EQNOY(1),Z(1),TEMPTIMESTEP) T'O S(2) RKDO T'O S(RKDEQO(O)) R CALC VARIABLES AND DERIVATIVES OF VARIABLES S(1) T'H RKZ1, FOR I=1-,lI*E.BRNO R CALC DERIVATIVES OF F'S RKZ1 Z(I)=O. W'R DELNB.G.O. Z(BRNO)=-2.*(F(BRNO+1)*(VW(BRNO+1)*SR(BRNO+1)-Z(19)/2.)-QBR/ 1ARH)/DELTAB O'E Z(BRNO)=-2.*(F(BRNO+1)*(VW(BRNO+1)*SR(BRNO+I))-QEFF(BRNO)/ 1ARH)/DELTA E'L T'H RKZ11, FOR I=BRNO+1, 1,I.G.5 RKZ11 Z(I)=-Z(I-1)-2**(F(I+1)*VW(I+1)*SR(1+1)-F(I)*VW(I)*SR(I)-QEFF 1(I)/ARH)/DELTA R CALC DERIVATIVES OF FR'S T'H RKZR, FOR I=1,1I.G.5 W'R N(11).E.l Z(I+5)=-Z(I+4)-2.*(FR( I+)*VWR (1+1)*SRR( +1 )-FR( I )*VWR( I )*SRR 1(I))/DELTAR O'E Z(I+5)=0O RKZR E'L R CALC DERIVATIVES OF Q'S T'H RKZ2,FOR I=1,1,I.G.5 RKZ2 Z(I+10)=(QQ(I)-Q(I))/TC R CALC DERIVATIVE OF P Z(16)=(QSUM(6)-QSUM(BRNO+1)+QBR-ARH*VW(6)*SR(6)*F(6))/PD R CALC DERIVATIVE OF MTOT R VW(1).L.O. CASE IS NOT COMPLETED W'R N(6)*GE.1 W'R VW(1).G.O. Z(17)=(VW(1).P.2.-VW(6).P.2.*((1.-F(6))+RSW*F(6)*SR(6).P.2*)) 1/DELTA-K 3 *VW(1).P.1.8-MDOT+G*HEAD T'H RKMD1, FOR I=1,1I.G6.5 RKMD1 Z(17)=Z(17)-K(1)*VWBAR(I).P.1.8/(1.-FBAR(I) ).P.2+G*FBAR(I) O'E Z(17)=(VW(1).P.2.-VW(6).P.2.*((1-F(6) )+RSW*F(6)*SR(6).P2.)) 1/DELTA+K 3 *.ABS.VW(1).P.1.8-MDOT+G*HEAD T'H RKMD2, FOR I=11,I.G.5 W'R VWBAR(I).L.O0 Z(17)=Z(17)+K(1)*.ABS.VWBAR(I).P.l.e/(l.-FBAR(I) )P..2+G*FBAR 1(I) O'E Z(17)=Z(17)-K(1)*VWBAR(I).P.1.8/(1.-FBAR(I)).P..2+G*FBAR(I) RKMD2 E'L E'L W'R N(11).E*O Z(17)=Z(17)+G*RISER*F(6)-K2*VW(6).P.1.8/(1.-F(6)).P.,2 O'R N(11).E.1 T'H RKMD, FOR I=1,1I.*G.5 RKMD Z(17)=Z(17)+G*R*FRBAR(I)-K2*VWRBAR(I).P.1.8/(1*-FRBAR(I) ).P. 12 E'L O'E

-142Z( 17)=0O E'L R CALC DERIVATIVES OF TW W'R N(7)*E.O Z(18)=0O O'R N(7).E.1 Z(18)=(TIN1-TINO)*EXP.(-TIME/TAU)/TAU E'L T'H RKZ21. FOR I=2,1,I.G.BBNO RKZ21 Z(I+17)=-Z( I+16)+2.*(Q(I-1)/AW-VW(1)*(TW( I)-TW( I-1 ) ))/DELTA T'H RKZ22, FOR I=BBNO+1 1,I.G.5 RKZ22 Z(I+17)=Z(23) W'R N(8).E.O Z(23)=0O O'R N(8).E.1 Z(23)=1.17*Z(16) E'L R CALC DERIVATIVE OF NBL W'R N(9).NE.O W'R DELNB1.L.1O.*RASCH Z(24)=-(2.*(Q(BBNO-BBNO (BBNO)*DELNB1)/AW-(1.+DELNB1)*DELTA*(Z 1(23)+Z(BBNO+16)))/(TSAT-TW(BBNO-1))+2.*VW(1) O'R DELNB1.L.l. Z(24)=-DELNB1*(2.*Q(BBNO)/AW-DELTA*(Z(23)+Z(BBNO+17)))/(TSAT1TW(BBNO))+2.*VW(1) O'E Z(24)=-(2.*(QBNO+Q(BBNO)+Q(BN )*(DELNB1-1.))/AW-DELNB*DELTA*(Z 1(23)+Z(BBNO+17)))/(TSAT-TW(BBNO))+2.*VW(1) E'L O'E Z(24)=0. E'L R CALC DERIVATIVE OF T R CALC DERIVATIVES OF C'S W'R N(2).E.O Z(25)=(RHO-BETA)*T/LSTAR+LAMBDA*C Z(26)=BETA*T/LSTAR-LAMBDA*C O'E Z(25) = (RHO-BETA) *T/LSTAR+LAMC Z(26)=0. T'H RKZ3, FOR I=1,1,I.G.6 RKZ3 Z(I+26)=BETA(I)*T/LSTAR-LAMBDA(I)*C(I) R DERIVATIVES ARE ALL CALCULATED R CALC C(O) FOR N(2)=1. ALSO LAMC LAMC=O. C=0. T'H RKY1, FOR I=1,19I.G.6 C=C+C( I) RKY1 LAMC=LAMC+LAMBDA(I)*C(I) LAMBDA=LAMC/C E'L R SHAPE SUBROUTINE* CORRELATE QQQ'S TO F'S R CALC FBAR, XBAR T'H RKFX, FOR I=11,I.G.55 FBAR(I)=(F(I)+F(I+1))/2. RKFX XBAR(I)=RSW*FBAR(I)*SR(I)/(1.-FBAR( I)) XR=F(6)*RSW*SR(6)/(1.-F(6)) W'R DELNB.G.O. FBAR(BRNO)=F(BRNO+1 ) * (l.-DELNB ) /2

-143XBAR(BRNO)=RSW*FBAR(BRNO)*SR(BRNO)*VW(BRNO)/VW(1) O'E FBAR(BRNO-1)=-F(BRNO)*DELNB/2. XBAR(BRNO-1)=RSW*FBAR(BRNO-1)*SR(BRNO-1)*VW(BRNO-1)/VW(1) E'L R CALC X-FUNCTIONS W'R N(1).NE.O N=1 T'H RKX1, FOR I=5.-1,I.L.2 T'H RKX1, FOR J=1 19J.G.3 X(N)=FBAR(I).P.J RKX1 N=N+1 T'H RKX2, FOR I=5,-1,I.L.3 T'H RKX2, FOR J=2, 1,J.G.3 X(N)=(FBAR(I)-FBAR(I-1)).P.J RKX2 N=N+1 E'L R CALC VW(1) W'R N(6).GE.1 VW(1)=MTOT-VW(BRNO+1)*(1.-DELNB)*(.5-2.*F(BRNO+1)*(l.-SR(BRNO 1+1)*RSW)/6.) T'H RKVW, FOR I=BRNO+1,1I.G.5 RKVW VW(1)=VW(1)-VW(I)*(.5-(2.*F(I)+F(I+1) )*(1.-SR(I)*RSW)/6.)-VW( I+1)*(.5-(F(I)+2.*F(I+1))*(1.-SR(I+1)*RSW)/6.) W'R N(11).E.O VW(1)=VW(1)-(1.-F(6))*ACR*VW(6)*RISER O'R N(11).E.1 T'H RKVW1, FOR I=1-1,I.G.5 RKVW1 VW(1)=VW(1)-(1.-FRBAR(I))*VWRBAR(I)*R E'L VW(1)=VW(1)/(4.+RISER+BRNO+ (1.+DELNB)* (.5-F(BRNO+1)*(.-SR(BR 1NO)*RSW)/6*)) E'L R SLIP SUBROUTINE. CALC SLIP RATIO W'R N(5).E.1 GG=G O'R N(5).E.2 GG=GZRO E'L T'H RKSP, FOR I=1,1,I.G.6 XC(I)=RSW*F(I)*SR(I)/(1.-F(I)) W'R N(11).E.1 XR(I)=RSW*FR(I)*SRR(I)/(1.-FR(I)) SRR(I)=SLIP.(N(5),GG,XR(I)) E'L RKSP SR(I)=SLIP.(N(5) GGXC(I)) T'H RKY2,FOR I=l,1,I.G.5 W'R N(1).NE.O.OR.N(4).E.1 QQQ(I)=SHAPE.(IX,B)/QP E'L FRACT(I)=1l+F(I+1)*(SR(I+1)-1.) R CALC NEW VW'S W'R DELNB.G.O. W'R I.L.BRNO VW(I+1)=VW(1) O'R I.E.BRNO VW(I+1)=(VW(I)+QBR/ARH)/FRACT(I) O'R I.G.BRNO VW(I+1)=FRACT(I-1)*VW(I)/FRACT(I)+QEFF(I)/(FRACT(I)*ARH)

-144E'L O'E W'R I.L.BRNO-1 VW(I+1)=VW(1) O'R I.E.BRNO-1 VW(I+1)=VW(1)/(1.-F(I+1)) O'E VW(I+1)=FRACT(I-l)*VW(I)/FRACT(I)+QEFF(I)/(FRACT(I)*ARH) E'L E'L VWBAR(I)=(VW(I.)+VW(I+1))/2. R CALC NEW Q'S R CALC QEFF W'R N(lO).E,1 QEFF(I)=Q(I)+(K(4)+K(5)*FBAR(I))*Z(16) O'E QEFF(I)=Q(I) E'L RKY2 QQ(I)=QQQ(I)*T W'R DELNB.G.O. VWBAR(BRNO)=VW( 1)*DELNB+(VW(1)+VW(BRNO+1))*(l.-DELNB)/2. E'L R CALC NEW VWR'S W'R N(11)*E.1 VWR(1)=VW(6)*ACR''H RKVWR, FOR I=1,1,I.G.5 FRACTR(I)=1.+FR(I+1)*(SRR(I+1)-1.) RKVWR VWR(I+1)=FRACTR(I-1)*VWR(I)/FRACTR(I) T'H RKVWB, FOR I=lglI.G.5 RKVWB VWRBAR(I)=(VWR(I)+VWR(I+1))/2. E'L R CALC NEW RHO W'R N(3).E.O REX=O. O'R N(3).Eo1 REX=RSIN*SIN.(6.2832*OMEGA*TIME) O'R N(3).E.2 REX=RRAMP*TIME O'R N(3).E*3 REX=RSTEP*(1-SURE) E'L W'R N(4).E.O RIN=-RHOZRO T'H RKRIN9 FOR I=1l,1I.G.5 RKRIN RIN=RIN-ALPHA*FBAR(I) O'E RIN=SHAPE (OXB)-RHOZRO E'L RHO=RIN+REX R CALC G W'R N(12).E.O GRAVE=1. O'R N(12).E.1 GRAVE=1.+GSIN*SIN.(6.2832*OMEGAG*TIME) E'L G=GZRO*GRAVE R CALC NEW TIN W'R N(7).E.1, TIN=TIN1+(TINO-TIN1)*EXP.(-TIME/TAU) R CALC NEW BOILING BOUNDARY

-145W'R N(9).NE.O BRNO1 BRNO BBN01=BBNO BRNO=NBL/DELTA+1 BBNO=BRNO DELNB=(NBL-(BRNO-1)*DELTA)/DELTA DELNB1=DELNB QBR =QEFF(BRNO)*(1.-DELNB) DELTAB=DELTA*(1.-DELNB) R THREE DIFFERENT CASES R CASE A. DELNB.G.1.-RASCH W'R DELNB.G.l.-RASCH F(BRNO+1)=F(BRNO+2)*DELTAB/(DELTA+DELTAB) BRNO=BRNO+1 DENBDELDELNB-1. QBR =QEFF(BRNO)*(l.-DELNB) DELTAB=DELTA*(1.-DELNB) R CASE B. DELNB.LE.1.-RASCH.AND.DELNB.G.RASCH O3R DELNB.G.RASCH W'R BRNO.L.BRNO1 F(BRNO+1)=F(BRNO+2)*DELTAB/(DELTA+DELTAB) O'E F(BRNO)=O. E'L W'R BBNO.G.BBNO1 TW(BBNO)=TSAT-(TSAT-TW(BBNO-1))*DELNB1/(1.+DELNB1) O'E TW(BBNO+1)=TSAT E'L W'R DELNB.G.(1.-5.*RASCH) WT=(DELNB-1.+5**RASCH)/(4.*RASCH) F(BRNO+1)=WT*F(BRNO+2)*DELTAB/(DELTA+DELTAB)+(1.-WT)*F(BRNO+1 1) E'L R CASE C. DELNB.LE.RASCH O'E TW(BBNO)=TSAT-(TSAT-TW(BBNO-1))*DELNB1/(1.+DELNB1) BBNO=BBNO-1 DELNB1=1 +DELNB E'L E'L W'R SURE*E.1, T'O BACK T'O RKDO S(2) CONTINUE R RECALCULATE TW(BRNO) AND F(BRNO) W'R DELNB1.G.1, TW(BRNO)=TSAT-(TSAT-TW(BBNO))*DELNB/(1.+DELNB) O'R DELNB.L,0O F(BRNO)=-F(BRNO+1)*DELNB/( 1.-DELNB ) O'E F(BRNO) =0 E'L R PRINT AND PLOT RESULTS W'R (TIME-TIME1).GE.SCALE P'S TIME, BRNO, BBNO, STEP W'R N(16).E.0 P'S Y(1)7..Y(55), Z(1)...Z(32) O'R N(16).E.1 P'S Y(1)...Y(55)

-146O'R N(16),E.-1 P'S Y(1)...Y(55) W'R N(3).G.O PLOTY(O)=REX/RNORM+24. E'L W'R N(12).G.O PLOTY(1)=(GRAVE-1.)/GNORM+24. E'L PLOTY(2)=ELOG.(T)*.86859+20. T'H PLOTF1, FOR I=1,1,I.G.5 PLOTY(I+2)=ELOG.(QQ(I)/1000.)*.86859+15. PLOTF1 PLOTY(I+15)=FBAR(I)*10. PLOTY(8)=P/5. PLOTY(9)=NBL/100.+13* T'H PLOTV1, FOR I=1,1,I.G.6 PLOTV1 PLOTY(I+9)=VW(I)/50.+6. O'R N(16).E.-2 P'S Y(1)...Y(55) W'R N(3).G.0 PLOTY(O)=REX/RNORM+24. EtL W'R N(12).G0O PLOTY(1)=(GRAVE-1.)/GNORM+24. E'L PLOTY(2)=(T-QQO)*QFACT/QQO+22. T'H PLOTF2, FOR I=11,9I.G.5 PLOTY(I+2)=(QQ(I)-QQO(I))*QFACT/QQO(I)+I+16. PLOTF2 PLOTY(I+15)=(FBAR(I)-FBARO(I))*FFACT+I PLOTY(8)=(P-PO)*PFACT/PO+15. PLOTY(9)= (NBL-NBLO)/NNORM+ 14. T'H PLOTV2, FOR I=1 19I.G.6 PLOTV2 PLOTY(I+9)=(VW(I)-VWC(I))*VFACT/VWO(I)+I+5. E'L W'R N(16).L.O PLOTX=TIME/TNORM W'R N(3).G*O EXECUTE PLOT3.($R$,PLOTX0PLOTY(0),1) E'L W'R N(12).G.O EXECUTE PLOT3.($G$,PLOTXPLOTY(l),1) E'L EXECUTE PLOT3.($*$,PLOTX9PLOTY(2),1) EXECUTE PLOT3.($+$,PLOTXPLOTY(8),1) EXECUTE PLOT3.($X$,PLTPLOTXPLOTY9)1) T'H PLT, FOR VALUES OF I=2,9,15 EXECUTE PLOT3.($5$,PPLTXOPLOTY(I+5) 1) EXECUTE PLOT3.($4$,PLOTXPLOTY(I+4),1) EXECUTE PLOT3.($3$,PLOTXPLOTY(I+3) 1) EXECUTE PLOT3.($2$,PLOTXPLOTY(I+2)91) PLT EXECUTE PLOT3.($1$,PLOTX.PLOTY(I+1),1) EXECUTE PLOT3.($6$,PLOTX,PLOTY(15)l1) E'L TIME1=TIME PRINT COMMENT $O******************************************** 1*************************$ E'L R CALC STEP W'R N(13).E.1.AND.TIME.G.TMIN

-147RATE=ABS.Z (25) W'R lO.*.ABS.Z(5) G.RATFE RATE=10.*.ABS*Z(5) W'R RATE.L..l STEP=STEPO*6. O'R RATE.L..2 STEP=STEPO*, O'R RATE.I.1. STEP=STEPO*4. O'R RATE.L..? STEP = STEP0 3 O'R RATE.L.10O STEP=STEPO*2. O'E STEP=STEPO E'L W'R DELNB.L.RASCH, STEP=STEPO E'L W'R TIME.G.TLIM WATCH=CLOCK.(0)-WATCH PRINT COMMENT $8$ P'S WATCH W'R N(16) L O PRINT COMMENT $1$ T'H FIGURE, FOR I=1,,II.G.FIG FIGURE EXECUTE PLOT4.(O0LABEL) PRINT COMMENT $15 E'L T'O START E'L T'O RKDO END E'M $COMPILE MAD, EXECUTE PRINT OBJECT, PUNCH OBJECT, DUMP EXTFRNAL FUNCTION ( K,X,B E'O SHAPE. INTEGER I: K, N W'R X(1).1 *.4 N=O O'E N= 1 E' L SFUN=B(N,K O) T'H LOOP, FOR I=1,1iI.G.18.OOP SFUN=SFUN+B(NK, I*X( I ) F'N SFUN E'N

-1.48List of Symbols in the Computer Program Symbol Definition and/or Implication A(I) Vector locations used in calculating the steady state water velocities and void fractions A(l,l) = Q(I)*(1.-DELNB)/ARH+F(I)*SR(I)/W(I) A(1,2) = -SR(I+l) A(l,3) = 0. A(2,l) = A(1,1)+VW() A(2,2) = l.-SR(I+l) A(2,3) = 1. ACR = AC/A ALPHA Linear void-reactivity coefficient ARH = AC.ps hv AW = AC.Pw B(I,JK) Vector locations of the coefficients used in the SHAPE. subroutine BBN0 Number of nodal point where the reference is made in the subcooled water temperature and non-boiling length calculations. Generally it corresponds to the region number in which the boiling boundary is found. BBNO1 BBNO at the preceeding computation time BETA(I) Delayed neutron fraction, Pi BLANKX, BLANKY Variables used for the grid set-up in the PL0T1. subroutine BRNP Number of node where the reference is made in the void fraction calculation. Generally it corresponds to the region number in which the boiling boundary is found.

-149BRN01 BRNO at the preceeding computation time C(I) Delayed neutron precursor densities, Ci C0UNT Number of succeeding iterations in which the convergence criterion is met DELNB Fraction of NBL in the BRN0-th region, 5nb. Could be negative in Case A.' (See Figure 5) DELNB1 Fraction of NBL in the BRN0-th region, 6nb. Could be greater than 1. in Case C. DEL(I) Relative differences of F(J) and W(J) in two succeeding calculations of the steady state iteration DELTA Core length of one region, A DELTAB Boiling length in the BRN0-th region, Ab DELTAR One-fifth of the riser length DH Increment of HEAD DHMIN Minimum value of DH DR Increment of RISER DRMIN Minimum value of DR DV Increment of VW(l) DVMIN Minimum value of DV EQN0 Number of differential equations F(I) Void fractions, fi FO(I) Void fractions at TIME = 0 FBAR(I) (F(I)+F(I+l))/2. FF(I) Void fractions at the preceeding calculation in the steady state iteration FIG Number of set of figures to be printed by the PL0T4. subroutine FN0RM Normalization factor for void fractions in the PL0T3. subroutine FRACT(I) = 1.+F(I+l)-(SR(I+l)-l. )

-150FRACTR(I) = 1.+FR(I+1)*(SRR(I+1)-l.) FR(I) Void fractions in the riser FRBAR(I) = (FR(I)+FR(I+1))/2. G Vertical acceleration, g GN0RM Normalization factor for g in the PL0T3. subroutine GRAVE Vertical acceleration relative to go GSIN Amplitude of sinusoidal variation in GRAVE GZR0 Acceleration due to gravity, go HEAD Pump head, in the unit of A, for the forced convection HEADO HEAD at TIME=O HEAD1 Additional HEAD after TIME=0 HMAX Maximum of HEAD HMIN Minimum of HEAD HSUB Subcooling, h ub II Number of iterations in the steady state calculation IIMAX Maximum of II IMAGE(I) Locations where the grid image is stored JJJJ J= 0 indicates the sign of M thas not changed. JJ = 1 indicates the sign was changed in the previous step. K(I) K(1)...K(3) are frictional factors when RISER=5., K(4) and K(5) are constants used to calculate QEFF(I) K2 = K(2)*RISER/5. K3 = K(3)*(RISER+5.)/10. LAMBDA(I) Precursor decay constants, ki 6 LAMC = E.iCi i=l LSTAR Equivalent to effective neutron lifetime, A

-151MID0T Derivative of momentum which must be made zero at the steady state MD0T1 MD0T at the preceeding calculation in the steady state iteration MT0T Total fluid momentum, Mtot N(I) Options for the program NBL Nonboiling length of the core NBLO NBL at TIME = 0 NBN0RM Normalization factor for non-boiling length in the PL0T3. subroutine NSCALE(I) Vector locations for the PL0T1. subroutine OMEGA Frequency of sinusoidal reactivity input 0MEGAG Frequency of sinusoidal vertical acceleration input P Pressure of the system, p PO P at TIME = O PD Denominator in the pressure equation, Equation (2-52) hS 6hw V S PS W P 5hs W vPs w I Pw} = M - + p- + hv tp2 p p PFACT Multiplier for P in the PL0T3. subroutine PL0TX Abscissa variable for the PL0T3. subroutine PL0TY Ordinate variable for the PL0T3. subroutine P0WER Total power output of the reactor QBR Effective power given to the boiling portion of the BRN0-th region QEFF(I) Effective power of the i-th region QFACT Multiplier for QQ(I) in the PL0T3. subroutine QP A normalization constant that gives Q(I) in the correct unit, = (Q(1)+ Q(2)+...+Q(5))/POWER

-152Q(I) Power given to coolant of the i-th region QQ(I) Power generated in fuel of the i-th region QQO(I) QQ(I) at TIME = O QQQ(I) Flux shape function of the i-th region QSUM(I) Power given to the coolant below the i-th node. R One-fifth of the riser height, in the unit of A RASCH A small fraction of A, enb RATE The larger of IZ(25)| and 10.*IZ(5)|, used in the STEP calculation REX External reactivity, pext RHO Total reactivity, p p ext + Pin RH0ZR0 RH0 at the steady state that must be compensated in the transient calculations, PO RIN Internal reactivity, causedby the steam void variation, in RISER Total riser height, in the unit of A RMAX Maximum value of RISER RMIN Minimum value of RISER RN0RM Normalization factor for RISER in the PLOT3. subroutine RRAMP Linear coefficient of xt when it is fed as rampfunction RSIN Magnitude of Pext that changes sinusoidally RSTEP Magnitude of pext that changes stepwise, p0 RSW = P /P SCALE Interval of TIME for printing out the results SRR(I) Slip ratios in the riser SR(I) Slip ratios in the core, yi STEP Step size of TIME for RKDEQ. subroutine STEPO Initially given STEP

-153SURE When 1, derivatives are calculated for initial check only. When 0, integrations of differential equations proceed. T Flux level, time function T(t) TAU Time constant for the inlet water temperature to change stepwise TC Fuel heat transfer time constant TEMP(I) Temporary locations to be used in RKDEQ. subroutine TIME Time in second after the perturbation is introduced, t TIME1 Preceeding TIME at which results are printed TIN Inlet water temperature, T in TINO TIN at TIME=O TIN1 TIN after TIME=O TLIM Maximum limit of TIME TN0RM Normalization factor of TIME for the PLIT3. subroutine TSAT Saturation temperature corresponding to P, Tsat T'W(I) Water temperature at the i-th node, Twi VMAX Maximum of VW(l) VMIN Minimum of VW(1) VW(I) Water speed at the i-th node, Vi VWO VW(1) at TIME = 0 VWBAR(I) = (VW(I)+VW(I+1))/2. VWN0RM Normalization factor of VW(I) in the PLIT3. subroutine VWR(I) Water speed at the i-th node of the riser VWRBAR(I) = (VWR(I)+VWR(I+1))/2. W(I) = 1./W(I) WT Weight factor to evaluate F(BRN0+l) WW(I) W(I) at the preceeding calculation in the steady state iteration

X(I) Power of void fractions and differences of void fractions of adjacent regions, used in the SHAPE. subroutine XBAR(I) = (XC(I)+ XC(I+1))/2. XC(I) Steam quality at the i-th node of the core, Xi XDEC Interval of PL0TX used to set-up the grid image XMAX Maximum of PL0TX XMIN Minimum of PL0TX XR(I) Steam quality at the i-th node of the riser YPSILN When all DEL(I) are less than YPSILN, convergence criterion is met for the steady state iteration. Y(I) Location of variables Y(o) = F(1) Y(5) = FR() Y(ll) = Q(l) Y(16) = P Y(17) = M0T Y(18) = TW(l) Y(24) = NBL Y(25) = T Y(26) = C Y(33) = QQ(l) Y(38) = QQQ(l) Y(43) = VW(l) Y(49) = QEFF(l) Y(54) = RH0 Y(55) = G Z(I) Locations of the time derivatives of variables with the corresponding subscripts

REFERENCES 1. J. A. Thie, Proceedings of the 1960 Idaho Conference on Reactor Kinetics, IDO-16791, pp. 107-117 (1962). 2. E. S. Beckjord, Trans. Am. Nucl. Soc., 3 (2), p. 433 (1960). 3. A. W. Kramer, Boiling Water Reactors, Addison-Wesley, 1958. 4. E. S. Beckjord, "Dynamic Analysis of Natural Circulation Boiling Water Reactors," ANL-5799 (1958). 5. J. A. DeShong, Jr., Nucleonics, 16 (6), pp. 68-72 (1958). 6. A. Ziya Akcasu, "Theoretical Feedback Analysis in Boiling Water Reactors," ANL-6221 (1960). 7. J. A. Fleck, Jr., J. Nucl. Energy, Part A, 11 (2/4), pp. 114-130 (1960). 8. K. 0. Solberg, "The'KJELLER MODEL' for the Dynamics of Coolant Channels in Boiling Water Reactors, Part I: Theory," KR-51 (1963). 9. Power Reactor Technology, 4 (2), p. 28 (March 1961). 10. R. L. Crowther and W. H. Wolf, Trans. Am. Nucl. Soc., 6 (2), p. 209 (1963). 11. R. L. Crowther and W. H. Wolf, Trans. Am. Nucl. Soc., 6 (2), p.210 (1963). 12. A. N. Nahavandi and R. F. Von Hollen, Nucl. Sci. Eng., 18 (3), pp. 335-350 (1964). 13. Albert Messiah, Quantum Mechanics, Vol. II, p. 740, John Wiley & Sons, New York, 1963. 14. H. Hurwitz, Jr., Nucleonics, I, pp. 61-67 (July 1949). 15. A. F. Henry, Nucl. Sci. Eng., 3 (1), pp. 52-70 (1958). 16. L. N. Ussachoff, Paper No. 656, Proc., 1955 Geneva Conference, Vol. 5, pp. 503-510, United Nations, 1956. 17. J. Lewins, J. Brit. Nucl. Energy Soc., 2 (3), pp. 326-334 (July 1963). 18. M. A. Schultz, Control of Nuclear Reactors and Power Plants, 2nd. ed., McGraw-Hill, 1961.-. -155

-15619. S. Kaplan and S. G. Margolis, Nucl. Sci. Eng., 7 (3), pp. 276-277 (1960). 20. A. F. Henry and N. J. Curlee, Nucl. Sci. Eng., 4 (6), pp. 727-744 (1958). 21. N. J. Curlee, Nucl. Sci. Eng., 6 (1), pp. 1-10 (1959). 22. S. Levy, J. Heat Transfer, Trans. ASME, 82, Series C, p. 245 (May, 1960). 23. R. W. Haywood et al., Proc. Inst. Mech. Eng., 175 (13), p. 669 (1961). 24. D. H. Brown, Paper No. 57-NESC-81, Paper Presented at the 2nd Nuclear Engineering and Science Conference, Philadelphia, Pa., (March 1957). 25. James E. Dallemand, "Stepwise Regression Program on the IBM 704," GMR-199, Research Laboratories, General Motor Corporation (Nov. 1958). 26. B. Arden, B. Galler and R. Graham, Michigan Algorithm Decoder, Computing Center, The University of Michigan, 1962. 27. T. Kanai, T. Kawai and R. Aoki, J. Atomic Energy Soc. Japan, 3 (3), pp. 6-16 (1961). 28. J. Randles, "Kinetics of Boiling Hydraulic Loops," AEEW-R 87 (1961). 29. C. K. Sanathanan, "Dynamic Analysis of Coolant Circulation in Boiling Water Nuclear Reactors," ANL-6847 (1964). 30. E. S. Beckjord and W. H. Harker, "Hercules I -- The Steady-State Calculation of Vertical Two-Phase Flow," GEAP-3261 (1959). 31. J. F. Marchaterre and M. Petrick, Nucl. Sci. Eng., 7 (6), pp. 525-532 (1960). 32. M. Muscettola, "Dynamic Model for a Boiling Water Reactor," AEEW-R 285 (1963). 33. G. L. West, Jr. and H. Nishihara, J. Joint Panel Nucl. Marine Prop., 6 (2), pp. 43-55 (Oct. 1962). 34. N. Zuber, Proceedings of the 1960 Idaho Conference on Reactor Kinetics, IDO-16791, pp. 189-208 (1962). 35. G. L. West, Jr. and H. Nishihara, Paper No. 63-WA-244, Paper Presented at the Winter Annual Meeting, Philadelphia, Pa., of the Am. Soc. Mech. Engrs. (Nov. 1963). 36. W. M. Rohsenow and H. Y. Choi, Heat, Mass, and Momentum Transfer, Prentice-Hall, 1961.

-15737. H. Schmidl, et al., "A Transfer Function Model of the HBWR Plant," HPR-5 (March 19g5). 38. H. P. Flatt, "The FOG One-Dimensional Neutron Diffusion Equation Codes," Unpublished, North American Aviation Corp., 1961. 39. S. Gill, Proc. Cambridge Phil. Soc., 47, pp. 96-108 (1951). 40. "VBWR Stability Test Report," GEAP-3971 (1963). 41. H. Christensen, "Power-to-Void Transfer Functions," ANL-6385 (1961). 42. A. Tapio Eurola, "On the Measurement of the Dynamic Properties on the Steam Void Fraction in Boiling Water Channels," ANL-6369 (1961). 43. F. H. Westervelt, Doctoral Thesis, The University of Michigan, 1960. 44. A. F. Henry and R. S. Willey, Trans. Am.Nucl. Soc., 2 (2), pp. 63-64 (1959). 45. See, for instance, Philip M. Morse, Thermal Physics, W. A. Benjamin, New York, 1962. 46. H. J. Amster, "A Compendium of Thermal Neutron Cross Sections Averaged over the Spectra of Wigner and Wilkins," WAPD-185 (Jan. 1958). 47. D. J. Behrens, "The Effect of Holes in a Reacting Material on the Passage of Neutrons, with Special Reference to the Critical Dimensions of a Reactor," AERE T/R 103 (1958). 48. D. J. Behrens, "Calculation of the Neutron Migration Length in an Infinite Lattice," AERE T/R 221 (1958). 49. D. J. Behrens, "The Migration Length of Neutrons in a Reactor," AERE T/R 877 (1956). 50. R. L. Hellens, "Neutron Slowing Down in Group Diffusion Theory," WAPD-114 (1956). 51. M. R. Fleishman and Ho Soodak, Nucl. Sci. Eng., 7 (3), pp. 217-227 (1960). 52. F. H. Clark, "Variation of Fast Fission Effect with Homogeneity in Infinite Slab Reactor," KAPL-743 (1952)o 53. A. M. Weinberg and E. P. Wigner, The Physical Theory of Neutron Chain Reactors, p. 669, The University of Chicago Press, 1958. 54. H. Soodak, ed., Reactor Handbook, 2nd ed., Volo 3, Part A, Physics, p. 225, Interscience Publishers, 1962.

UNIVERSITY OF MICHIGAN 3 9015 03483 0235 55. G. I. Bell, Nucl. Sci. Eng., 5 (2), pp. 138-139 (1959). 56. S. Glasstone and M. C. Edlund, The Elements of Nuclear Reactor Theory, D. Van Nostrand, 1952. 57. L. C. Just, C. N. Kelber, and N. F. Morehouse, Jr., "An Analog Computer Model of a Multiple Region Reactor," ANL-6482 (1962).