ML20011E255

Draft Arrotta:Advanced Rapid Reactor Operational Transient Analysis Computer Code Documentation Package,Vol 1:Theory & Numerical Analysis

ML20011E255
Person / Time
Site:	McGuire
Issue date:	04/30/1986
From:	Eisenhart L S. LEVY, INC.
To:
Shared Package
ML19310C691	List: ML20011E252 ML20011E253 ML20011E254 ML20011E255
References
NUDOCS 9002120372
Download: ML20011E255 (136)

Text

{{#Wiki_filter:- M. 5V. fy4 ,C 7 ~~l ~' i p4 R Y.J31.g..; [ L.L ' g...:g .l' $5lls- ? 5 L 4 t COMPUTER CODE MANUAL System Documentation ARROTTA 1 , s r s, s Advanced Rapid Reactor 0 s 1 Operational Transient Analysis L Computer Code Documentation Package l -. Volume 1: Theory and Numerical Analysis '.,.4 April, 1986 ARROTTA And-Tne ARROTTA Computer Code Documentation Package M Were Prepared Under Research Project 1936-6 L-Prepared by: S. Levy Incorporated 3425 South Bascom Avenue Campbell, California 95008 Principal. Investigator: Laurance D. Eisenhart L 7, b%; -ARROTTA was originally initiated as ANITA under Research Project 1936-4 to Brookhaven National-Laboratory f ~ Prepared For:- -Electric Power Research Institute.. Inc. l-3412 Hillview Avenue y Palo Alto, California 94304 EPRI Project Manager: Walter. J. Eich Analysis Methods and Verification Program ,^- Nuclear Power Division 1. s e m; n-a 9002120372 900129 f@ L PDR ADOCK 05000369 M [ P PDC

e T ,I; l l f: i .i i NOTICE I This report was prepared by the organization (s) named below as an account of j - work sponsored by tne Electric Power Research Institute Inc. (EPRI).- Neither EPRI, members of EPRI, the organization (s) named'below, nor any person acting on i their. behalf: (a) makes any warranty or representation, express or implied, with respect to'the accuracy, completeness, or usefulness of the information con-U tained in -this report, or that the use of any information, apparatus, method, or 4 process. disclosed in this report may not inf ringe-privately owned rights; or (b)- j assumes any-liabilities with respect to the use of,-or for damages resulting from the use;of, any information, apparatus, method, or process disclosed in ' this report. Prepared by: S. Levy Incorporated Campbell. California i l f

l_ ..j ' 3 ,' ji' "T %*yy!NN ?' TJi W "~N"h"Np/ 54t?"e*%g d l '. : g i', a { \\ 1 4 -I^ } i ..j 1 j k -li .n "},. ABSTRACT ( d i fi. * .a - ) tk .-I e .,1 4 i Yt 4 p j i / ? c Wil t e i k 5

{4.

mn 4,, y n 6 )i 1 k F 6 4s '.i4 D (- h 4 h ~{,. f-f ?* 2 e s.- '6 L n . h.. f, 4x 'y y. b C-g . r E y' n D . :i .h b 't'.. c digtSAM A %f g 4',jg'g,_

e.., ,, - a g :), ' .i ' l.h 1 - t ( x q .s ACKNOWLEDGEMENTS' .c ') b .e t P t ~ t ..:f 1 IA 4 -i 1 ? l 9 i s \\ i 3 4 l' I ? . a $,l pe hP, '~g 3 ?C 1 y; }e g. 1 e. t [ 1 y, t. r t f;'.I i. 2 g i L. E . 1 i 1 s a 1 9 ,h, '. 'Y'N*Tk Aby&sNka_va>ganu;& ' _v _.,___.___i_,_-_____ 'i

0 4 ')'t

$I iJ.;

g, TABLE OF CONTENTS 'Section Page-r

1.-INTRODUCTION-l-1

2. Neutronics Model 2-1 i

2.l Introduction 21 2.2 Xenon lodine Formalism 2 - 2.3 Cross'Section Formalism 2-7 -2.4 Control: Rod Movement. 10-2.5. Time Behavior. 2-11 m 2.5.1 Time Dif f erencing 2-11 2.5.2 Derivative Estimation 2-11 2.5.3 Transient Flux / Leakage' Extrapolations 2-12 2.6 Analytic Nodalization Methodology 2-12 2.f.1 Derivation of-Spatial Coupling 2-12 2.6.2 Discontinuity Factors in a-Easier Setting 2-18 1 2.6.3 Formal Leakage-Equations 2-20' 2.6.4 Suitable Algebraic Form for the Leakage & Flux Operators 2-21 2.6.5' Transient Nodal Equations 2-25 2.7. Numerical Analysis Considerations 2-28 2.7.111ntroduction and Definitions 2-28 l ~ 2.7.2 Leakage Approximations 2-29 2.7.2.1 Introduction 2-29 2.7.2.2 Development of' the Leakage. Equations 2-32 2.7.2.3 Explicit Leakage Iterative Schemes 2-34 2.7.2.4 Implicit Leakage Tratment Schemes 2-37 2.7.2.5 Leakage Computational Strategy-2-40 2.7.3 Outer Iteration Numerical Analysis .2-42 2.7.3.1 Introduction 2-42 1 2.7.3.2-Proolem Formulation 2-42 ,c 2.7~.3.3 Chebyshev Acceleration Scheme 2-44 2.7.3.4 Chebyshev Parameter Estimation 2-48 2.7.4 Inner Iteration Algorithm 2-51 2.7.4.1 Introduction 2-51 1 2.7.4.2 Iterative Methodology. 2-52 l 2.7.4.3 S.0.R. parameter Estimation 2-53 j 2.7.5 Reactor nonlinear Iterations'and Searchs 2-55 2.7.5.1 Introduction 2-55 l 2.7.5.2 Search Procedures 2-57 2.7.5.3 Reactor Iteration Procedures 2-58 2.7.5.4 Acceleration of the Reactor Iteration Procedures 2-59 j i j

ht L t y? 1 b TABLE-OF CONTENTS (Continued)' 4 ' Section-- Page-

3. Fluid Dynamic.Model-3-1 a;

3.r Introduction 3-1 3.2-Two-Phase Flow Equations 3-4 / 3.2.1 Single-Phase liquid. Flow (Tc2' < Ts) 3-6

3.2.2 Single-Phase Liquid Flow (Tc2 > Ts) 3-9 3.2.3 Two-Phase Flow Subcooled Boiling 3 3.2.4 Two-Phase Flow Bulk Boiling _

3-12 3.2.5 Single-Phase' Flow Vapor Flow 3-13 3.3 Numerical Algorithms 3-15 3.3.1 Transient Algorithm 3 1 g' 3.3.2 Steady-State Algorithm 3-22 -j 3.4 Boundary Conditions 3-26 q 3.5 Special Features _ 3-26

4. Heat Conduction Model 4-1 4.1 Introduction.

4-1

j 4.2 Steady-State Heat Conduction 4-2

]

4.2.1 Pellet Region and Gap 4-2 4.2.2 Clad Region 4-10 4.3' Transient Heat Conduction 4 4.3.1: Average Clad Temperature 4-18 1

4.3.2 Average Fuel-Temperature 4-21 4.3.3 Clad Inner Surface Temperature _ 4-23 . 4.3.4 Clad WallETemperature 4-25 .4.4~ Clad Surface: Heat Flux 4-27 4.4.1 Steady State Clad Surf ace Heat Flux 4-27 'i 4.4.2 Transient Clad-Surface Heat Flux-4-28 4.4.3 Special Features 4-29 4.5 Fuel Rod Material Properties 4-30. 5.. References; 5-1 l Appendix: Derivation and Evaluation of Spatial Coupling Equations A-1 3 (.

LIST OF FIGURES Figure Page 1-1 Chiculational Flow Diagram for ARROTTA 1 l 2-1 Mesh Example 2-1B 3-2 A Discretized Flow Channel 3-15 4-1 Typical Fuel Element Temperature Distribution 4-3 4-2 Thermal-hydraulic Cell 47 4 j o

cr b i 2 1: g ~ LIST OF TABLES - -Table-Page 2-1 Maximum Overrelaxation Factor vs. Iter 6 tion Index 2-3 +. t i.A - i

Ii s Section 1 ARROTTA - GENERAL DESCRIPTION AND CAPABILITIES ARROTTA - Advanced Rapid Reactor Operational Transient Analyzer - is a three - dimensional XYZ neutronics computer program for both static and transient appli-cations that includes the effects of-the thermal hydraulic feedback mechanisms. ARR0TTA neutronics are based-on the Analytic Nodalization Method as developed for QUANDRY (J_) in EPRl_RP 1936-1. The thermal hydraulics model.in ARROTTA is taken directly from the BEAGL program (2) as developed under EPRI RP 1761-18. ARR0TTA can solve steady state simulations, xenon (and samarium) transients, and kinetics transients. The-analytical nodalization method will generate an exact solution to the neutron diffusion equations if a) the node is considered as a homogeneous entity and b) the shape of the transverse leakage function is considered to be of a known quadratic form. The equations in the limit of small mesh revert to the same limit as the-standard flux centered finite difference equations. The time derivative-term is' solved by a standard finite differencing while the time dependent cross section behavior is handled either implicitly or -semi-implicitly. Normal zero flux and zero current boundary conditions are allowed - independently on all six sides of the problem. ARR0TTA uses a full two group model for the diffusion equations and up to six groups of delayed neutron pre-cursors. Initial conditions that may be calculated include a fixed input state-point, a. criticality search based on critical boron concentration, and a search to criticality. based on a rod bank position. In all cases these solutions are nonlinear since they are explicit functions of the thermal-hydraulic conditions which are coupled to the pointwise cross section values which in turn are coupled to the pointwise power. Also, the nonlinear effects of the xenon and iodine. poison (as well as the samarium and prometheium) may also be handled on a pointwise basis. Collectively, all of the above feedback mechanisms (nonlineari-( 1-1 ~ ~ -

ties) are nandled on the reactor iteration level during eigenvalue problems, in actual transient calculational sequences, the possible feedback mechanisms) are delayed by decoupling the nonlinearity between successive timesteps. The thermal-hydraulic model is comprised of a fluid dynamics model and a heat transfer model. The fluid dynamics is described by. an inhomogeneous, non-i equilibrium, two phase flow that is based on energy splitting between the pha-ses. One-dimensional (axial) conservation equations are solved by a marching technique for multiple parallel flow channels taking into account six possible flow regimes. Consequently, flow reversal and cross flow phenomenon are outside i the scope-of tnis fluid model. The relative flow split amongst the channels may be specified by the user but is time invariant. The system pressure, inlet 4 coolant temperature, and inlet flow rates may all be specified as time depen-dent. The heat conduction model is based upon spatially averaged, time-dependent equations for the averege pellet temperature. Given the properties of the media, the gap and the inner and outer clad wall temperatures are then determined. In the time-dependent case these heat conduction equations are solved by a semi-implicit formulation. Notable Code Characteristics: 3 Neutron diffusion equations solved in full two group rigor over a rectangular geometry limited to 40 x 40 x 40 nodes (64,000 total)- ARROTTA will generate initial conditions the include either a control rod group search or a soluble boron search Thermal-hydraulic feedback mesh solved at up to 64,000 nodes Three separate automatic levels of approximation to the neutron leakage operator Time discrezitation of the neutronics is either fully implicit or semi-implicit Fixed or variable timestep size algorithm 1-2

b 4 Automatic adaptive computation of acceleration parameters e Flexible restart procedures

• ~

Flexible spatial editing routines Assembly Discontinuity Factors allows greater accuracy for a . fixed number of neutronic nodes Comprehensive and computationally efficient cross section representation formulation Flexible and accurate control rod representation algorithm A simplified calculation flow diagram for ARROTTA is shown in Figure 1-1. 4 -Detailed expla' nations of the calculational modules of ARROTTA are contained in Volume 1, the Theory and Numerics Manual, of the ARROTTA Computer Code Documentation Package. ARROTTA is written completely in FORTRAN and is currently available on both CDC and IBM equipment at the UCCEL data center. ARR0TTA requires the EPRI-LIBRARY (3) set of scratch disk drivers. Even though some of these routines are machine dependent, EPRI-LIBRARY is being maintained as part of the expected software base. The normal ~ implementation of this code is to use the FORTRAN-4 compiler on the CDC system (4) and the FORTRAN-H EXTENDED and ENHANCE 0 compiler using the' Automatic Precision increase (AUT0 DBL) feature (AD option 15331) on-the IBM system (5). The ARROTTA Computer Code Documentation Package consists of five volumes. Volume '1, the Tneory and Numerics Manual, discusses the general equations and approxi-mations, geometry representations, constitutive relationships, and numerical analysis, both theoretical and observable, techniques and procedures employed in-the various models in ARROTTA. Volume 2, the User's Manual, provides detailed descriptions and examples for using ARR0TTA including complete descriptions of the input data that is required and cogent descriptions of the output data possible. Volume 3, the Programmer's Manual, provides the complete implemen-1-3

t V l s1 tation procedures for the ARROTTA computer code and describes the general cooing cf ARROTTA_ including functional descriptions of all' ARROTTA subroutines ano-functions,-FTB scratch file storage, and COMMON variables. Volume 4, the Applications Manual, and Volume 5, the Modeling Guidelines Manual, will be ] completed at a later date. 1 This volume present the mathematical theory and modeling of' both the neutronics and the coupled thermal hydraulics phenomenon contained in ARROTTA. In addition, the computational problems of the solution numerical analysis are explored in 1 - order that the reader of this volume might have an appreciation of the underlying procedures and techniques that h6ve been used. Chapter.2 describes the Analytical Nodal Method as implemented in ARROTTA including a sketch of the derivation of both the static and transient neutronic equations. Other equations pertinent to the neutron-formalism are also documented. Additionally, emphasis f is' placed on the numerical analysis techniques used especially the adaptive choosing of estimates for the needed acceleration parameters. Chapter 3_ explores. -l the fluid dynamics model, an in-homogenious, non-equilibrium, two-phase flow model model based on energy splitting between the'11guid and vapor phases. One-dimensional conservation egyations for the mass and energy.are solved in the i axial direction under th'e hypothesis of no crossiflow by a marching technique._ Chapter 4. discusses the fuel rod heat conduction model. It is based on the spa-tially averaged, time dependent heat conduction equations for the average pellet ] temperature and the clad inner and outer surface temperatures. These equations i are solved semi-implicitly. i 4 i l 1-4 ~ u.

l _e. >, e o. I ru 1 1 I ... i.. u s..,. n.. I i.. i Q.

a t,i w..

. i.,,, i i t.. x l.n = s i . an,i .., ;.q C1 p.i i.u n,.. t ui. l... c. t,i un 1 [ i i ,/ v ,t, i ..t. c sec u... nr. i I use.i.t. =t 4. i..nin... J t/ {' l L14 tit 8 Wit l.litt' ] t tel l V tt ItWD) I W t. hebttteill l f 1P.mittet llevt99 mitt l l j.Due' Di 'ibettet l { a.* tittetbwt1De ] l = v l$1087SW1.ThtN1byt'.wittl l l1MMitM TMN1 M'h'tM l { 0btP9httl IO90h... l %? ) l wette.1 steep f tie. If toewestec j g l ti....i.. uit ie ti e j n....,..... m... t,....,,....J 't,..,,,,, 1-i u,..,.t,..,,,,. ..,,.,t. t g tei toi i *e n h/ \\/ I W Olbrium 8ence feest.tt. Seisct i t. Ste9 Lengt* l if a, .t. n u l i,....) c.i..i. u....u , I t,i = g/ %/ tC l A 1-5

Section 2 ARROTTA NEUTRONICS MODELS

2.1 INTRODUCTION

The original starting point for the development of the ARROTTA equations is the standard two group three-dimensional neutron dif fusion equation in both its transient (Eq. 2.1-1) and static (Eq. 2.1-2) forms. All arguments that assess the suitability of using these equations will for the sake of brevity u ll simply be deemed to be valid for the purposes of this exposition. v'I I = 7'

• Diqgq+y*I(1-6)[v{fi 3+v{f2 8)

(2.1-la) 4 2 + A C d d del 84 4> 1 2=?*D Y'2*It2 0 '2

• 1 'l (2'I*IDI v

2 R subject to the initial conditions'given as 4 1*Y*l[vl i V41+[g1 1+v[f2'2) (2.1-26) -9 ' D fi $ 2 Y'2+It2*2*1*1-(2.1-2b) -v + D R where the terminology is standard and is repeated here for clarity. g is the neutron group index,- g = 1 or 2 0 is the total number of delayed neutron groups - d = 1.....D $ is the scalar neutron flux in group g [/cm cm sec) Cd is the delayed neutron precursor density for delay group d [/cm em em) D is the group g diffusion coefficient [cm) g 2-1

t I i {t is the total matroscopic cross section for group g {fg is the fission macroscopic cross section for group g [R is the macroscopic removal cross section for group 1 to 2 transfer l x is the prompt neutron fission sper,trum xy is the delayed neutron fission spectrum for delayed group d i is the mean number of neutrons emitted per fission in group g vg Y is the critical system eigenvalue A is the decay constant for delay group d [/sec) d B is the fractional yield of decay neutron precursors in group d per d fission g is the neutron velocity for group g [cm/sec) v The above equations are subject to either zero flux or zero current external boundary conditions on each of the six problem boundaries independently. Although it is possible to implement other external boundary conditions, currently ARROTTA only considers these two types. At all material interf aces the continuity of the flux and the normal component of the current is assumed. Since this analysis will always assume that there are exactly two groups and that the group cutoff is nowhere near the lowest probable fission energy even for the delayed fission, all fission source spectra cccurring in the above defining equations will be replaced by the following obvious thermal reactor spectrum: [1,03. All coefficients (cross section functions) in the above equations are shown as functions of both space and time. However, in the actual reactor problems of interest they can also be functions of a variety of feedback variables usually based on the calculations of various thermal hydraulic phenomenon. These genera-lizations will be given later in the section on ARROTTA cross section modeling. For ease of presentation only, it will be assumed that the above diffusion equations are linear. l l The precursor equations for each delayed neutron group and given by Eq. 2.1-3 [ with the corresponding initial conditions specified by Eq. 2.1-4 2-2 4 ~.

!~ BC d=1... D (2.1 3) I"I}81*"If2 '23 **d C Y' 6 f d*g d 6

1. )

d=1,...,0 (2.1 4 ) d*y I"If1'l+"If2 C 2 d The derivation will begin classically by assuming that the reactor volume has been partitioned into a three dimensional regular array of homogeneous right rectangular parallelepipeds with the following grid structure: X={x(i): x(0) < x(1) <... < x(i) <... < x(l) } Y = {y(j): y(0) < y(1) <... < y(j ) <... < y(J) } Z={2(t): z(0)<2(1)<...<2(k)<...<z(K)} An arbitrary node (i.j,k) in the reactor volume represents the space defined by the following inequality: x(1-1) < x < x(i) and y(j-1) < y < y(j) and z(k-1) < z < z(k) The nodal volume is defined similarly as i v(i.j,k) = [ r(i) - x(i-1) ) * [ y(j) - y(j-1) ) * [ 2(k) - 2(k-1) ) or alternatively as h (i)

• h (j)
• h x

y r When referring to an arbitrary direction the generalized terminology - u where h (t)

• u(t) - u(1-1) will be used. The v and w directions will be the other two u

directions in this terminology when appropriate, in the ensuing derivations, we will work with the static diffusion equation f or convenience unless the time dependent form is needed for a particular reason. The net currents (J - a two group vector) on the f aces of an arbitrary nooe i i -(1,j,k) as a function of the two transverse directions can be represented by the defining Eq. 2.1-5 o Jui,j,k (v.w) = -[D,j,k) d[$(ug,v,w)) / du (2.1-5) i ) l 2-3 l l

p l then the correspending integral over the $Urf ace of the V-w f ace will give the proper expression for the f ace averaged current J (i,j,k) as Eq. 2.1-6 i J (I'd' ) ' h /dw b ')' u i'j*k u h (j) h (k) y I In order to derive the nodal dif fusion equation, we must first integrate the original two group diffusion over an arbitrary node (i j,k) obtaining Eq. 2.1-7 i-for the static case. h (j) h,(k) [J, -Jx ) + h,(i) h,(k) [J -J 3 (2.1-7) y y + h,(1) h (j) [J,k } - J,k)+yi,j,k){T j 4 T* =i V(i,j,k)xv{f & and 4 is the volume averaged nodal flux as given below as Eq. 2.1-8: + t(i.j,k)= /dxJdyJdz[4(x,y.2)) ( 2.1-8) V(1,j,k) Thus' the net leakage L (i,j,k) can be presented as the difference of the two u face currents in the same direction as Eq. 2.1-9 L (1,j,k) = J (i,j,k) - J (1-1.j,k) ( 2.14) u u u Thus rewriting the static equation we have Eq. 2.1-10 h h l, + h,h hyy{(= h*h[h*v{f hhb +hhL + &+hhhx y z I 'I y2x x2y R gy (2.1-10) while the corresponding transient equation with time still a continuous variable yields Eq. 2.1 11 2-4

y t v {h,hl + h,h,L7+hhh (E

• 1 ~

"1)#} I2*}~lI} y 2 T R f f + hhh [A X xyg ddd d=1 2.2; XENON IODINE FORMALISM The fission product poisoning equation treatment in ARROTTA is patterned directly from the SIMULATE-E (J,) program. The ARROTTA cross section model f or-mally expects that both the iodine-135 to xenon'-135 and promethium-149 to samariam-149 explicit nodewise poisoning are calculated in ARROTTA. However, only the primary chain is present for the Pm and Sm. As with the SIMULATE-E program the user may wish to define an effective direct yield for Pm-149 so that i the indirect yields are approximately taken into account. Conceptually, two separate and distinct sittple two-isotope linear chain differential equations have been installed in ARROTTA. If the problem definition contains these isoto-pes, then ARROTTA just assumes that the contributions to the macroscopic cross sections have been excluded from the cross section input and are to be carried explicitly. The initial steady state calculation of an ARROTTA sequence will produce the equilibrium concentrations of the isotopes, If the user desires a realistic simalation of a reactor configuration for a transient analysis, then the proper initialization methodology includes the pointwise equilibrium xenon and samarium poison distribution, in this case the equilibrium concer.trations that are calculated on the initial steady state are never et anged. They are just used as a transient independent fixed poison. In the classic xenon analysis scenario af ter the initial equilibrium is found, the calculation of the time dependent shape assumes that the fission (appear-ance) rate and the absorption (disappearance) rate can be assumed to be very slowly varying relative to the time step size being used to follow the event. This hypothesis can be checked only by taking a smaller time step for the same event and assessing the calculated difference. 2-5 l t

~ p Since the two f amilies of ordinary differential equations are exactly of the same form, they will be given only once with a mother nuclide m (1135 or Pm-149) and a daughter nuclide d (Xe-135 or Sm-149). It should bc noted that if the user wants to track the Be-La chain, then it is easy to use the correspond-ing Ba and La parameters in place of the Sm and Pm parameters in the ARRDTTA cross section model input. The basic differential equations for the mother nuclide is given by Eq. 2.21 and the daughter nuclide d by Eq. 2.2-2. l dN(m)/dt Y(m) A(m) N(m) (2.2-1) i = and dN(d)/dt Y(d) A(d) N(d) + A(m) N(m) (2.2-2) = where the disappearance rate for nuclide k (m or d) is given by Eq. 2.2-3 A(k) = A(k) e,3 e(1) + o2

• (2) 52*2*3)

+ a and the yield rate for nuclide k is given by Eq. 2.2-4 Y(k) y(k) [ If3 4(1) + If2 'I2) ) I2'2*A) = where ((g) are the absolute fluxes [n/bn-sec) for g=1,2 A(k) is the decay constant [/sec) for nuclide k N(k) is the number density of nuclide k [ atoms /bn-cm) y(k) is the fission yield fraction for nuclide k The equilibrium solutions (dN(k)/dt = 0) to the above equations 2.2-1,2 are given by Eq. 2.2-5 and 2.2-6, respectively. Neq(m) Y(m) / A(m) (2.2-5) = 2-6 4

I. and Y(d) ) / A(d) (2.2-6) N,q(d) { Y(m)

• A(m) / A(m)

+ = Even though simplifications exist for particular mother-daughter pairs, the entire solution has been programmed in ARROTTA for completeness and analogy to SIMULATE-E. The time dependent solution proceeds given a set of nuclide concentrations at the beginning of the time step of length At where all the parameters are eva-luated at the beginning of the time step. The new concentration of the mother nuclide, N(m, At) it given by Eq. 2.2-7 while the new concentration of the daughter nuclide, N(d At) is given by E4 2.2-8. 4 [1 - e(m))

• Y(m) / A(m)

(2.2-7) N(m,0)

• e(m)

N(m,At) + = and N(d,0)

• e(d) + ci * [1 - e(d)) + c2 * [e(m)-e(d))

(2.2-8) N(d At) = where exp ( - A(k)

• At ) for k = m or d e(k)

= Y(d) / A(d) _ [A(m)

• Y(m) ] / [ A(m)
• A(c) ]

c1 + = and X(m)

• N(m.0)

A(m)

• Y(m) c2 A(d) - A(m)

A(m) * (A(d)-A(m)) = 2.3 CROSS SECTION FORMALISM The cross section formalism in the ARROTTA program is patterned af ter the Brookhaven National Laboratory experience with the two-dimensional kinetics codes that culminated with EPRI-BEAGL, Although the form of the model is 2-7

somewhat arbitrary, it does, in actuality, simply account for most of the variation of the cross sections throughout a very large class of nonlinear tran-sients without having an unduly complicated cross section feedback represen-tation as given by Eq. 2.3-1. t -{ =R [ A + B * (X - X') + C * (X - X')**2 ) l + (1. - R) * ( A' + B' * (X - X' ) + C' * (X - X' )**2 ) [d-{ / d T,) (T,- T,i) (2.3-1) + (d { / d /T ] (sort (T ) - sqrt(T )) { + g g f Function (N, N,,, N, Ngg, pm} + N b g where R is a monotonically increasing function rf the fraction of the node that. is controlled and R(0) = 0 and R(1) = 1. The basic variation is a quadratic one in void for a BWR and in water density for a PWR. Two separate and distinct sets of quadratic coefficients exist for the rodded and the unrodded configurations simultaneously. Consequently, the j existence of control rods does nothing to complicate the composition overlay structure. Also, t'he moderator temperature feedback is defined to be a linear. phenomenon. The Doppler feedback is treated in the fast group only and is assumed to be linear in the square root of the absolute temperature. The Doppler and moderator feedback mechanisms are assumed to be independent of the existence j of control rods in a given composition. This representation is true separately for each of the nine cross sections including the kappa sigma fission cross sec-tion. Thus, no relation is assumed between the value of the average number of neutrons per fission (nu) and the average energy release per fission (kappa). Additionally, for the absorption cross sections and for the removal cross sec-tion only a composition dependent boron-10 cross section triplet is available for problems that contain an explicit boron representation for PWR's. Also, a composition dependent pair of absorption cross sections are available for xenon, iodine, samarium, and promethium in XISP problems only. 2-8

i \\ l Thus, to model this variation fully in ARROTTA, 54 numbers are needed to mocel theifunctional behavior for the X-variation of the moderator - six f or each of nine cross section types, Another nine coefficients are needed for mooerator temperature feedback, and Doppler will account for another five entries. Since f I all of these entries require knowledge of the reference feedback variables, another three entries are required for the reference fuel and moderator tem-peratures and the reference value of the X variable. Therefore, with no special treatments 71 input values are needed to fully specify the allowed macroscopic cross section variation in an ARROTTA composition. Basically, there only exist three types of kinetics parameters that are required for each ARR0TTA composition - inverse neutron group velocities, precursor group yield fractions, and precursor group decay rates. Again, it is assumed that these parameters are independent of control rods as well as feedback conditions. The assumption in ARROTTA is that each of these parameters is permitted to be composition dependent, in the usual kinetics problem with six delayed precursor groups there are fourteen kinetics parameters to be input per composition. Xenon (and samarium) " transients" can also be ar.alyzed with ARROTTA. Actually, ~ only the isotopic equations for xenon, iodine, samarium, and promethium are solved in a time dependent f ashion. Both the neutronics and the thermal-hydraulics are solve as a series of steady state problems. The total input required is precisely 18 numbers per composition since it is assumed in ARROTTA that all properties in a XISP problem are composition dependent. Two group nu's are needed to get the correct explicit fission product production. The isotopic yields and decay rates for each of the four nuclides of interest are required. Also the f ast and thermal absorption cress section for each of the four nuclides of interest even though a common assumption has been to ignore the fast contri-bution. Also the correct pointwise poison distributions for use in the initial con-ditions for kinetics analyses can be calculated as part of the ARROTTA initiali-zation procedure. This insures a more realistic description of the conditions 2-9

that exist in a reactor at the start of the event. The equilibrium xenon an:1 samarium distributions are found as part of the eigenvalue feedback iteration. Then, these pointwise distributions remain constant during the transient. i 2.4 CONTROL ROD FORMALISM Since the representation in ARROTTA is implicitly based upon an effective homo-genization of the assembly, the control rod inventory in a core is described relative to the assembly mesh. Consequently, even if the solution is to be found over a more refined computational mesh the definition of the control rods do not. First, the individual rods are defined to be part of rod groups. These group definitions are totally arbitrary reif tive to any actual plant and control system. They are meant to f acilitate the user input. Second, the actual posi-tions at the initial conditions are specified; and, third, the time behavior of each rod group is then specified. There are four possible control rod configurations that are possible in ARROTTA: a. top (PWR) or bottom (BWR) rod insertion, and b. rod is completely contained within the given assembly position (e.g. RCC " finger rodlets" or homogenized rod representation) or a cruciform type rod that is inserted between the actual assemblies. Control rod groups are defined by an initial position from the bottom of the problem and a series of pairs of " time, velocity" pairs that govern the movement _ of each given rod group by an assumed linear velocity equation as given by Eq. 2.4-1 for the time interval t(j-1) < t < t(j) while the velo-city is constant in that time interval. Care is taken if a given com-putational timestep goes over to defining velocity intervals to make the appropriate adjustments to Eq. 2.4-1. z(t) = z(t(j-1)) v(j-1) *[t t(j-1) ] (2.4-1) + l 2-10 j l

1 i 1 l t 2.5. TIME BEHAVIOR j ) I 2.5.1. TIME DIFFERENCIN3 \\ In the simulation of the time dependent field equations (the neutron diffusion equations) and ordinary differential equations (the precursor equations) encoun-tered in ARROTTA, standard finite differencing techniques will be utilized. The standard first order expression for the approximation of the derivative is used as given in Eq. 2.5-1. df/dt ) = [ f(t(j)) - f(t(j-1) ] / [ t(j) - t(j-1) ) (2.5-1) Also a function or combination of functions that is time dependent is to be eva-luated by the " theta" method for time integration as given by Eq. 2.5-2 6 f(t(j)) + ( 1 - 0 ) f(t(j-1)) (2.5-2) f (t) = where'f(t) is usually a time dependent reaction rate. Note that 0 = 1/2 is central differencing and 6 = 1 is fully implicit differencing. B = 0 which is fully explicit is not allowed due to mathematical stability considerations. ~2.5.2 DERIVATIVE ESTIMATION In the derivation of the quadratic leakage equations it is necessdry to estimate the derivative of the flux and precursor densities in terms of the given func-tion itself in order to be able to complete the formalism of the Analytic Nodal Method, it must be emphasized that this estimator is only used in the analytic leakage coupling integration. It is, therefore, deemed appropriate to estimate J the derivative as a frequency proportional to the original function as given by' Eq. 2.5-3 below, f(t(j)) (2.5-3) df(t(j))/dt = w where the w-estimate comes from the preceding time steps by Eq. 2.5-4 In [ f(t(j)) / f(t(j-1) ] / [ t(j) - t(j-1) ] (2.5-4) w = 2-11

V I, For_the first time step when the transient has just been initiated, the assumed f requencies are zero. Sufficient information is stored on the restart file to allow ARROTTA to continue the problem with the proper frequency information. 2.5.3 TRANSIENT FLUX / LEAKAGE EXTRAPOLATION At each time step in the transient calculation, ARROTTA must decide on the flux and leakage distributions to use for the initial guess to begin the iterative procedure. Usually the converged distributions f rom the just calculated time step are a reasonable initial guess. However, given that we have already calcu-lated the pointwise f requencies associated with these distributions it is possible to do much better. The initial fluxes at each node in each group are extrapolated by the current estimate of the f requency of the node in that group as shown by Eq._ 2.5-5. ((t(j+1)) = ((t(j))'* exp [ w * (t(j+1) - t(j) 3 (2.5-5) The corresponding leakages in the case of the explicit leakage formulation and the leakage source in the case of the implicit formulation are extrapolated similarly using the f requency for the flux directly as given in Eq. 2.5-6 L(t(j+1)) = L(t(j))* exp [ 'w * (t(j+1) - t(j) ] (2.5-6) 2.6. ANALYTIC NODALIZATION METHODOLOGY 2.6.1 DERIVATION OF SPATIAL COUPLING If we can reduce the defining equations 2.1-10 or 2.1-11 as applicable to one direction at a time by integrating over the other two transverse directions, it will be_ possible to derive a set of generalized one-dirensional dif fusion - equations that essentially contain the same amount of information but, hope-fully, in a more solvable or approximable form. Integrating over y and w to obtain the equation for u yields Eq. 2.6-1 for the static equation 2-12

1 m s i l l f [.: 2. 2 .2 -h h,D (U(u)-D(< $>+<'2(>}+hh

( (u) y y

Bu By Sw a xv[g( (u) (2.6-1) j

• Y hh y

l-l where i l t (u) = h h, ff((u,v,w)dvdw y y and (d 4 / dv2> =ff ((u,v,w)dvdw 2 dv Defining the continuously u-dependent v-leakage by the expression 2.6-2 h,L (u) = -D < d2 (/dv2> (2.6-2) y 4 integrating the above expression now in the direction of prime interest u gives Eq. 2.6 4 L (t) = f L (u) / h (2.6-3) l J (1+1) - J (t) = y y u y y Thus giving the nodal f ace-averaged v-directed net leakage. The total, transverse leakage to the direction u per unit u divided by the f acial area.h h, is given y as Su(u) by Eq. 2.6-4 L,(u) / h, (2.6-4) S (u) = L (u) /'h + u y y Combining, _the one-dimensional diffusion equation for the U-direction now has a familiar form given as Eq. 2.6-5 t 2-13 i

.a L e i L i- .2 ~ t (u) + ((T ~ T X"E)'u(u)=-S(u) (2.6-5) f u y Similarly, for the transient case, the equation becomes Eq. 2.6-6 i 2 I1'0)X"I)'u(ut)+ f h (u(u t) -D..it(u,t)+([T*T f u D A xd d(u t) = S (ui ) C t + d u (2.6-6) t However, here we make two assumptions that are needed for the time dependent equations in order that these parabolic partial differential operators may be assumed to be behaving pseudoasymptotically. We must state that the flux and the precursor derivatives can at least for a given time step be approximated by a frequency-like multiple of the original unknown as previously intimated ~by Eq. 2.5-3. This gives in particular Eq. 2.6-7 for the group fluxes and Eq. 2.6-8 for the precursor concentrations. d4 (2.6-6 ~ dt 9 '9 9 = 1,2 and dC e

• C (t)

(2.6-8) dt d d d = 1,...,0 but since the defining rate equation for the precursor contentration is-given by Eq. 2.1-3 r Y' 8d"1$*ACdd" o d = 1,...,0 (2.3-1) f t this implies that the precursor concentrations can be eliminated from this ana-lytical form since modest algebra yields Eq. 2.6-9. 2-la

6 C (u,t) = Y'I Sd " f * (U't) I I*d

• c) d=

1,...,0 (2.6-9) d Therefore, the one dimensional transient based diffusion equation is given by Eq.'2.6-10 6d 31 v[f} (2.6-10) -D ((ut)+{[+w V~I - [(1-6) x + [ g p p Su d "d* d Y $u(u,t) = - S (u t) u Now the terms in the static and transient have exactly the same form, alDeit somewhat more complicated expressions at a given time t in the transient, if we remember that the node has been chosen to be homogeneous, there' ore giving constants in the given node for all the coefficients in the differential equation. If there existed a form for the function S (u) (or S (u,t) at a given u u time t) we would be able analytically to integrate this equatica. The approxima-tion originally due to Finneman (8_) is to assume that the spatial dependence of each transverse leakage function can be represented by quadratic polynomials as given by Eq. 2.6-11 whose coefficients are determined by the average transverse leakages from the given nodes and its two nearest neighbors. Thus, L[i] (x).. l[i-13 [i-13 (x) ((i] p[i] (x) ((i+13 [i+13 (x) (2.6-11) p p and the polynomials p(x) have the integral properties such that the integral of the polynomial over the center interval is the mesh width and the integral of the polynomial over either end interval is zero, i The transverse leakage source function S (u) across a given node in a given u direction u has become by ansatz of the form of Eq. 2.6-12 where the values that appear in Eq. 2.6-12 are average nodal net rates. 2-15 t

U 4' 3; E i 1-1 t+1 g + [Tu,3 - Tu ) p (u) + [Tu,g - Tu ) p (u) (2.6 12) S b) = u g g g g g g u The expansion quadratic polynomials can be expressed by Eqs. 2.6-13 and 2.6-14 w:lere the two functions are simply referred to as p (s) and p'(s) where s is the normalized distance within a given node in the u-direction, i.e. s = [u - u(t)) ~ / h (A)* u + b* 5 + c s (2.6-13) p,(s) = a ano p"(s) = a' + b* s + c' s (2.6-14) f 2 i The coeffittents appearing in Ecs. 2.6-13 and 2.6-14 are a function of the mesh regularity in the particular u-direction only. Let h a h (1), g = h (1-1), and f u u = h (1+1) for convenience, then given the temporary d of Eq. 2.6-15, the needed u coef ficients are given by Eqs. 2.6-16 to 2.6-21. d o h'f g / 6 + h (f 9. f 9 ) / 2 + h [f g/6 + f f /2 + fg /6) 3 2 2 22 32,f9)/6 (2.6-15) 23 +h[f9 then 4 3 2 2 3 d a' = h f g / 6 + h gf /3+h9f /6 (2.6-16) d b* = 2 h*f g / 3 - h g f / 3 = h g f3 3 2 2 /3 (2.6-17) 1 4 3 2 -d c' = h f 9'/2 + h g f /3 (2.6-16) 3 2 23 d a+ + -h f 9 /6-h9f/6 (2.6-19) 4 23 d b* * -h f g / 3 + h 9 f / 3 (2.6 20) and 2-16

d c' s hfg/2+hff/2 (2.6 21) Conceptually, the derivation proceeds first to find particular solutions for each of the constant, linear, and quadratic transverse leakage source terms sub. ject to the zero boundary condition and to combine them with the homogeneous solution that was subject to a fixed value of the surf ace flux and current boun-dary conditions. The analytic solution is then integrated over the entire node and a relationship between node averaged fluxes and face averaged currents is obtained in terms of the surface flux at the right of the given node. This pro-cess is repeated but now going f rom the same surf ace to the end of the lef t side in order to obtain another equation relating node averaged fluxes and f ace-everaged net currents in terms of the same surfece flux. Since both equations are expressed in terms of the common surface flux it can be eliminated f rom con-sideration since it is not an unknown in the nodal balance equation. Conse-quently. a relationship now exist between node-averaged fluxes, f ace averaged net leakages and the interface current. Doing this same procedure at the next nodal interf ace will allow a similar set of relations. Taking the difference of these two will finally allow one to obtain the desired equation relating node-averaged fluxes and f ace-averaged leakages. Thus, through a long and involved tedious process given in Appendix A, the interface surf ace flux can be be elimi-nated between the adjacent nodes giving a resultant expression that relates the leakage, L, ih the direction of interest, u, to the three closest average fluxes, t, and the five closest transverse direction leakages, S. At first this derivation can be thought to assume that this surf ace flux is equal on both sides of the nodal interf ace. However, in general this is not true. Only the fine mesh heterogeneous flux is actually continuous. Advanced homogenization techniques dues to Koebke (9,) show that the Koebke equivalence theory does not require the nodal surf ace fluxes to be the same on either side of the interf ace since one side of the interface could contain a multiplier that forces the con-tinuity. In the Smith (3,H) modification to this theory where there appears to be a factor on either side of the interface we shall define " discontinuity fac-tors" f*(l) and f*(i) with Eq. 2.6-22 such that they trivially force the con-tinuity by 2-17

L l thet(i) = f*(i) * (*(i) = f (i) * * (i) (2.6 22) in ARROTTA it shall be presumed that these discontinuity f actors have been a priori computed by the user and are available es input. i. 2.6.2 Dist0NTINVITY FACTOR $ IN A EASIER $ETTING The concept and use of the discontinuity factors in the coupling equations in the context of the Analytic Nodal Method leads to very complicated equations. However, the concepts can be equally gleaned from a simpler system such as the mesh centerec finite difference equations. The face average flux will be expressea in terms of its coupling with and without the added complexity of the discontinuity f actors in a two dimensional setting for ease as given in Fig.21 below. hx(1-1) hx(i) i-1,j i,j t x(i-1) x(i) x(i+1) Fig. 2-1 Mesh Example The equations for the f ace-averaged current on the x(i) f ace are given by Eq. 2.6 23 and 2.6-24 without discontinuity factors J<x(i)> = -2 DN ) [4(i) - t<x(i)> ] (2.6-23) h (i) x 2-18

i and J<x(i-1)> = -2 DN-1) [4<x(1-1)>-I(1-1)) (2.6-24) h(M) l x However, with the addition of the discontinuity f actors these equations become generalized to 2.6-25 and 2.6-26. -2 0(i) = g7(g),,<x(g), j g-(i) ) (2.6-25) Jcx(i)> h (i) x i and -2 0(1-1) g,<x(g)37f+(1-1)-I(1-1)) (2.6-26) J<x(i)> = ) h (i-1) x In order to simplify the following the ensuing algebra let a and b be defined by Eq. 2.6-27 and 2.6-28, respectively, t -2 D(i) / h (i) (2.6-27) a = x and -2 D(1-1) / h (i-1) (2.6-28) b- = x then the ' resultant ef fective interf ace flux can be written as Eq. 2.6-29 a

• T(i) + b T(i-1)

(2.6-29) $<x(i)> = b / f*(1-1) + a / f*(i) while the actual f ace. average interf ace current is given by Eq. 2.6-30 f*(i)*I(i)+f*(i+1)*I(1-1) (2.6-30) J<x(i)) = + f"(i) / a l f*(i-1) / b Whereas the discontinuity f actors in the setting of the Analytic Nodal Method of the ARROTTA code represent a more difficult algebraic formalism, the type of 2-19

t i i generalized nodal coupling given above contains the essence of the modeling in ARROTTA. 2.6.3 FORMAL LEAKAGE EQUATIONS t i Folding together the relationships thusfer concerning the average leakage in a given node we find that the leakage is related to the three closest fluxes and the five closest transverse leakage sources in the given direction sought. This formalism yields an equation of the form of Eq. 2.6-31. t+1 1+2 L (t) = F II)*#IA)

• tat-2 u(t)
• S (1)

(.6-31) I b u tel-1 u u Given that the nodal averaged fluxes are availcble, we can now rewrite all of the leakage relationships as a three equation system in the unknown leakages. The leakage operators can be rewritten as the matrix given by Eq. 2.6 32. -l h G h O l F y x 2 x x x y -1 -1 h G -1 h, G, L F 4 (2.6-32) = - x y y y -1 -1 h G, h, G x 2 2 2 where each G is a pentadiagonal operator in the v-direction and each F is a y y tridiagonal operator in the V-direction. It must be emphasized that in both the stttic and transient applications the form of the above matrix equation is exactly the same; but, some of the actual terms are not. 2-20

1 4 I T t 2.6.4 SUITABLE ALGEBRAIC FORM FOR THE LEAKAGE AND FLUX OPERATORS The Analytic Nodal Method analytically solves a series of one-dimensional dif-fusion problems in the process of obtaining the actual three-dimensional problem of interest. These one dimensional problems have been reduced to a statement of a leakage balance equation in the three-dimensional problem of interest. This section will explore the structure of the spatial coupling equations in order that the more fundamental equation behavior will become apparent. This section assumes a f amiliarity with Appendix A - The Detailed Derivation of the Spatial Coupling Equations. Therefore, for each of tne three perpendicular directions (u=x,y, and 2) 7 submatrices have been defined as a function of the cross sec-tions in the node and the mesh length in the direction of interest. Dif ferences in the cross sections in the submetrices dependent on whether a static or tran-sient in being solved have already been addressed. The submatrix indices will be dropped in this section for the sake of clarity, but the D and E submatrices will still carry the + and - distinctions. Additionally, for each face of the node let the discontinuity f actors defined as a two by two diagonal matrix where fi(-) and fi(+) will denote the left and right faces of the ith node the the current direction of interest u. Repeating Eq. 2.6-31 as a starting form, the leakage equations are given by 1+1 1+2

• 1 O (A)
• S (1)

(2.6-32) L (t) = t F (*)*'(*) u u u u t=t-1 tot-2 The general forms of the F's and G's will now be given as Eqs. 2.6-33-40. Momentarily consider only those nodes that are neither boundary nodes nor adja-cent to boundary nodes. Therefore, the following regular form f or the coef-ficients are realized. It should be noted that when discontinuity factors are not being used these equations simplify somewhat since all the f-matrices become the two by two identity operator. 2-21

(2.6 n, (2.6 34) B1 R g g.1 (2.6 35) i+ t t1 }B 1 p& L + 1+1 f 1+1 ~* -(R (2.6 36) i i t F g i+1 f g,1 (2.6 37 ) R y t+1 p 1 O1 g R g 03 (2.6 33) t+2 + g i L~1 ~ t 0 R +1 L +T3 (2 6 39) t t+1 1 [l 1 G" 3R t t + U +1 t [Jt g L, 1 ,~ E t Gg 3 Rg (2.6 40) + T,1 g [lt t+1 p Gg O. u, R. I, J. K, nd a 0,E+. E ,f+,n Eqs. 2.6 1 lastly. saatrices 6 nd a defined. Thes A. B, C, D+,as give K g,y L+2 R g,y be nts rms yet matrice coefficie = to known L te nsion 2.6 47 (2.6 41) interinediateof the e pa terms akage through x T are give the kn Eqs. 2.6 41 's the in le n own N' with ed with ing nd 2.6 21. Thus giv ~ a ombin .f c 3'

• hro gh

~+ f'g A3 u l ~ m%y~

~ j 1 i I i 1 D' 1 - c' 1 E* 1 3 (2.6-42) I f't-1 [ -a' 1 ~ D -b 0 = 1-1- t-l 1-1-1 1-1 L I i f{ [-a[C -b'gj - c'g g] (2.6-43) D E' U = g g fi [(1-a)i-at)Ct-(bi+bi)Di-(c + c't) E'g (2.6-44) Tt = 1 1 f - c* lf) M-N ] t.1 [ -a I C -b ,l l, ,l t,l L,I L,l t t-1[-(1-ai.1-a{.1)C-1+(b{.1+b.1)0{.1+(c{.1+ci.1)E{1 ~ J f t g t (2.6-46) and [a'gg+b{D*g +c{E{] (2.6-47) g f *g C K = P Now all is well defined except for the impact of the external boundary con-ditions. Let us now assume that the line under consideration has N nodes. A zero current boundary condition at the lef t implies that the matrix R1 = 0; similar-ly, at the right of the line it implies that R +1 = 0. A zero flux boundary N = 0 or A 1 = 0, respectively, i condition is represented equally simply by A g o Therefore, the impact of the boundary conditions on the flux coupling terms (Eqs. 2.6-33 to 2.6-35) of the original Eq. 2.6-31 is very straightforward. However, the impact of the boundary conditions.on the transverse leakage coupling (Eqs. 2.6-36 to 2.6-40) is not so straightforward since the coupling i is pentadiagonal rather than tridiagonal. Not only are the bourdary nodes spe-cial but the nodes that are adjacent to the boundary also need special treat-ment due to the quadratic shape assumption inherent in the transverse leakage source term. Consequently, the eight cases (flux or current condition with the j bouncary node or adjacent to the boundary node and lef t or right end) will now be exhibited for the total leakage source coupling 'Z' of 2.6-31. These 2-23

'a L equations for the-case of the zero current are derived assuming a mirror sym-metry of the flux iSape outside of the domain of interest. For the cases of the zero flux condition the derivation assumes that the flux is zero everywhere out. side the domain of interest rather than assuming a negative symetry through a zero value on the problem outer surf ace.

1. Zero current at lef t boundary I = R E01*U1+

2 U)3 + R bI +T)S2*NE5223 l 2 1 1 2

2. Zero current at right boundary 3[J,g+V)b-1-R3 [In-1 + i *K3b Z

-R Q -l N R gn b = g N n n N n

3. Zero flux at lef t boundary

[R (J +V ) - R (1 + T )) S + [R (1 +T) RK)S +RKS 2 = 2 1 2 1 1 1 1 2 1 2 11 2 223

4. Zero flut at right boundary 2

= -R Q S +[R 0 - R (J +U)) S N N-1 N-2 N+1 N N N-1 N N-1 +[RN+1(J + U ) - R (1 +T)) S N N N N-1 N N S. Zero current adjacent to lef t boundary [RQ - R (Q + J + U ) ) S + [ R (J + U ) - R (! +T ) ) S Z = 3 +[R(I+T.)-RK3S + RKS

6. Zero current adjacer.t to right boundary 2

= -R 0 S + [R0 R (J +U ))S 2-24

g~7%= L N-1) ) $ N1(! +[R(J +0) +T R N-1 N2 N N-1 N R K 3$ +[R(1 +T +K) i-N N-1 N N N-1 N-1 N l 7; Zero flux adjacent to lef t boundary Some as the general case Eq. 2.6 31 with the matrix that multiplies S set te g tero.

8. Zero flux adjacent to right boundary Same as the general case Eq. 2.6-31 with the matrix that multiplies S.g set g

to Zero. lP 2.6.5 TRANSIENT NDDAL E0'JAT10NS The detailed exposition of the transient nodal ecuations can nov be presented. The tine and space diffusion equations are given by Eq. 2.1-1 with the precursor equations by Eq. 2.1-3 which are repeated here for completeness in matrix form. V"If1*7*DV4-[T'*T (I*6) X V t' (2'I~II ~ et D +fy XC Addd 4 and BCd Y'I Bd"I'*ACcd, d=1,...,D (2.1-31 f 5t Let an arbitrary series of time increments at(n) be given. The system of governing equations can be then written much more compactly as Eq. 2.6-4B for the flures and Eq. 2.6-49 for the precursor concentrations. 2 25

n; '). i. i n !i i' D $(t) = f(t)t(t) - L(t) + ]) A C (t) (2.6 43) V~ od i _ 8t os 6 L 8 C (t) = M (t)t(t) - C(d(t) d=1,...,0 (2.6-49) d d I where = h 3,(t) + h, G (t) 3 L It) f L(t) y y u I 8 It is certainly poss ble to introduce a variety of " theta" differencing for-y malism - the most pnilosophically straightforward is to allow a different single value of theta f or each neutron group and each delayed precursor f amily where this single value is inferred by some other auxiliary condition. This will in general yield a 2

• D by 2 + D nonsymmetric matrix in order to calculate an l

optimum value for etch dif ferential equation. This was demonstrated by Porsching and Henry (R) not.to be a practical approach. Consequently, we will adopt the E common " numerical folklore" of having at most two values, one for the fluxes, O and one for the precursors, E.. Using the thete differencing formalism as f p given previously by Eq. 2.5-2 when applied to a product and the standard dif. ference approximation for the first derivative as given by Eq. 2.5-1, the tiuxes and precursor concentrations necome Eq. 2.6-50 and 2.6 51, respectively. i 2-26

4 t { t"*1 - ("} = e F"*I "*1 + (1.e )F"e" (2.6-50) V*I I t g g At - n - (gl - (1 - 0 ) L f 9 D n+1 n I'pd II*'p)C} + A C d d and j n+1 "*1 + (1 - e ) H t" (2.6 51) n+1 n n I {C -C} = eM ( f o d o fe 4t n n+1 n I'pd

• II*'p)C}

d=1,....D C A d o The precursor equation 2.6-51 can be rearranged in order to separate the new concentration as given by Eq. 2.6 52. k n+1 n C = (1 + A ot 'p) $El*AAt Il * 'p)) C (2.6-52) o d n d n d n+1 n + AinE'fd ' II ~ 'f N t ]} d=1....D M d The consequences of having the form 2.6-52 imply that during the flux solution of a given time step, the C"*1 term can be replaced by its constituents as given in 2.6-52. Also, this means that the precursor concentration update is computationally a generalized edit that occurs after the fixed flux calculation even though it formally appears to be intrinsically in the discreZetized tran-sient equation. The actual fixed flux equations are given by Eq. 2.6-53 for the two energy groups encountered. The leakage equations are exactly of the same form as in the eigenvalue problem. 2-27 m

f-Mh sq jr 4 , :1 g bYi y# s D -A At"6 ~ n"1 [l. y-1 --e' f{F"*I

{

]),n+1~ t d D

At d ' I I

• d 'I 'p n

n r, . + e ' L + F. 1,--(),9)g n n f D A at"e n r [" I. y-I -(l.-e){F"-[ d P g Mb d At d=1 1+A at *p n d n D 'x n i +1 C y dal 1+ A at 0 N d-d np. (2.6-53)f 2.7 ARROTTA' NUMERICAL ~ ANALYSIS CONSIDERATIONS i -2.7.1. INTRODUCTION AND DEFINITIONS-Tne.ARR0TTA neutron iterative processes will be discussed in this chapter.: The' } - j basic iterative concepts will be~ introduced here for clarity. The following four-1 concepts of. iterative approximations as ~ used in ARR0TTA will be discussed.in-- this chapter: g [ a. leakage approximations andLiterations l < b;- neutron inner ' iterations c.- neutron outer iterations ~ d. nonlinear or reactor iterations. AnLeigenvalue calculation in ARROTTA'can simultaneously entail all four of the above iterative processes. Any time step of a transient:calculatian, however, q only entails the first two' levels of iteration, i .6.' Y The' leakage oporator involves the leakage (difference of adjacent currents in n ' given coordinate direction) as it affects the basic flux equations. In ARROTTA, an inner iteration simultaneously solves the two group neutron diffusion equation with a fixed term on the right hand side. The neutron outer iteration 2-28 4 i 4 i

~ y 1- ?n-1 solves allows' the fission source expression on the right hand side to converge to the critical statepoint with its associated critical eigenvalue given a constant set of material properties. The reactor or nonlinear iteration allows. -the material properties to change due either to feedback or changes in the control' state of the reactor (boron or control rod search) as well as the changes in-the neutron coefficient matrices that result from the explicit use of the system eigenvalue. .2.7.2 ARROTTA LEAKAGE APPROXIMATIONS ~2.7.2.1 Introduction 1 The basic flux and lea'kage equations of ARR0TTA for the entire problem that is j to be solved can be symbolically given by the matrix equation 2.7-1. A B flux Source (2.7-1) F C leakage 4 0.0 The flux equations contain two unknowns per node (one for each group) and the -leakage equations contain 6 unknowns per node in the 3-dimensional setting (one -for each group and each orthogonol direction). In an eigenvalue problem, the l ' source term in the flux equation is the fission source term, while in a tran- .sient solution process the source to the flux equation is the aggregate of all the previous timestep terms plus the precursor term. There are three levels of treatment of the leakage source to the flux equation present in the ARROTTA program. The conventional approach is to add the leakage equation as another level in the iterative processes and explicitly iterate. L ARROTTA uses-a line Gauss-Seidel formulation for this process while the original -QUANDRY: code uses a point Jacobi formulation. The iterative process (whether it .be Gauss-Seidel or Jacobi)'is basically simple because of the exceedingly simple structure of the matrix C of Eq. 2.7-1. The C matrix structure is simple because 2-29 i ___ n _

g l i the leakage equations are algebraic relationships between each other rather than ] differential equation relations as was shown in the previous section on the neutron equation derivation. To understand the genesis of the other leakage approximations schemes in ARROTTA it is first necessary to look at the matrix C of Eq. 2.7-1. If the leakages are ordered as all X then all Y then all Z components, the matrix C = 1 - K where 1 l is the identity matrix and K has nonzero entries far from the diagonals (at least 2

• number of nodes). The actual ordering of the leakage vector within ARROTTA is somewhat more complicated. All leakages within a given row come together ordered as all X-leakages for that row, followed by all the Y-leakages, then followed by all _the Z-leakages. Then all rows within a given plane; and finally all possible planes. This ordering is partially determined by the data structures within ARR0TTA rather than solely by the numerical algorithms used in the program. For the purposes of the analysis exposition of this class of methods, however, it is appropriate to think of the simpler leakage ordering.

All of the solution methods assume that the spectral radius of K is less than unity. This'is precisely the same assumption that allows one to use either of the above explicit iterative treatments. The second equation in the matrix Eq. 2.7-1 can then be rewritten as Eq. 2.7-2 [ l - K ] leak - F

• flux (2.7-2)

= This is solvable directly if the inverse of [I - K] is computable. It is well known that the Neumann series expansion for [1 - K] converges if and only if the spectral radius of K is less than unity. Consequently, since we have assumed this valid for the explicit method we will assume that it still is valid. This type of methodology is dubbed an " implicit" approximation type to distinguish it from the. iterative " explicit" class. The true solution of Eq. 2.7-2 is given by Eq. 2.7-3 l 2-30 l L 1 l'

r: u {l - K)'l = 1 + K'+ K2+K3 +... (2.7-3) y The implicit method in ARROTTA is BALUNS (Brookhaven Approximate Leakage Using Neumann Series) which has been shown to be comparable in accuracy with the explicit treatments. However, the computational complexity (operations count) of the BALUNS methodology is f ar less than the explicit formulation. The essence of the BALUNS approximation is the assumption given by Eq. 2.7-4 which is the first order Neumann series approximation. [1 - K]'1 = 1 + K (2.7-4) Experimentally, we have coserved that the implicit treatment and the explicit treatment are comparable in accuracy; however, the implicit leakage updating is much faster. Consequently, internally in ARROTTA even when a user requests the explicit leakage treatment, until certain partial convergence criteria are met, the leakage.updatirg is implicit regardless of the input requested. It is possible.to construct yet another implicit method albeit somewhat dege- - nerate. -The zeroth order Neumann series approximation is trivi311y given by Eq. 2.7-5 '[I - K]"I 1 (2.7-5) This approximation is much faster than either of the other two treatments but-lacks the essential-information of the second equation of the matrix 2.7-1. Consequently, all ARR0TTA calculations (that are restarts of partially converged previous problems) begin with this " truncated" mode of approximation until the l proper level of partial convergence has been achieved. Then, if the user has requested either the implicit or the' explicit treatment, ARROTTA automatically switches to the implicit mode either forever if the implicit mode was requested or'until the next convergence plateau is reached if the explicit mode was requested. 2-31

~ ( 'p J, 4 2.7.2.2-Development Of Th'e' Leakage Ecuations The analytic nodalization operators for the ARROTTA neutrorics can te expressed as the following four-balance equations when given as a general matrix operator Eq. 2.7-6. ~ h,h T' ~ ~ "kM & or S~ W h n I' hhl yg g 7 F -1 h *IG h,~1G L 0 x y x x x .F h *IG, h "IG L 0 y x 7 y y F h *1G h ~IG -1 L 0 [ 7 x 3 y 7 7 where h iS=the' vector of meshes in the U-direction u G is the pentadiagonal leakage operator in the u-direction u L is the leakage variable in _the u-direction u W is'the neutron spectral operator of the two group diffusion equation and M is the two group fission source operator. The iterative properties of the inner iteration-flux solvers that were used-in-the early experimental QUANDRY work, however,-did not converge well since tne submatrix W lacked diagonal dominance. Consequently, the simple transformation Lof. the first equation of Eq. 2.7-6 to Eq. 2.7-7 by combining judiciously the given four equations as given below did wonders for the convergence behavior. Ed. d. J. J ) [$, L, L. L F = k M $ or S-(2.7-7) x y z x y 2 2-32

h. [.,cp R where 9 h,h F, J =: W +:hhf +hhF + yy xgy J = h,G h G,. + y x y ' + h G, O hG y" j zx c r and J =hG + hG q z x 'I Given the fluxes from some flux iteration treatment, we now have a three equation. system in the ' unknown leakages. Consequently, let us concentrate on the leakage formulation for_ the moment. The leakage operators can be rewritten 'as the three by three system'given'by Eq. 2.7-8. l l h),* 1'G h *1G L F 4 l -1 x y x 7 x F (2.7-8) h *1G l_ h G L =- x y z y 7 j

4 h ~lG

' h -IG .-I L, F 4 x g y y 3 At this point in the exposition, a notational change is in order so that.the j t . matrix properties themselves will be highlighted rather than the equation struc- ] ture f rom' the dif ferential operator. Let P, Q, R, 5, T, and V be defined in term of the operators of Eq. 2.7-8 such we get the simpler looking Eq. 2.7-9. I ~! -P - -Q L F '* ~ '~ ~ x x F 4 (2.7-9)' -R-1 -S L = y y -T -U -I L F, 4 2 iz 2-33

2.; W %1 :J.~ l i 1

lf' the matrix of Eq. 2.7-9 were computationally feasible to invert, then the answer would formally be given by' Eq. 2.710.

L .l. -P -Q~ "I F ~ e x x -R 'I- -S F 4 -(2.7-10) .L = y y L -T -U 1 F 4 2 z ~ We will-investigate the ramifications of two separate-classes of approximation techniques for 'oDtaining the inverse of the the set of equations given sym. l bolically by Eq. 2.7-10. First, explicit iterative. updating of the system will- - be' advanced both in the: form of Jacobi and Gauss-Seidel techniques. Second, a new class of approximations, dubbed implicit, based on the. von Neuman series - expansion will be presented. Then the various relationships between them.will 7 be explored. - 2.7.2.3' Explicit Leakage Iterative Schemes 4 Classical iterative techniques of Jacobi and Gauss-Seidel will now be strate-gized for.the system of equations' given before as' Eq. -2.7-9. The three leakage - equations may be written in the following Jacobi equation nomenclature as given below in Eq. 2.7-11a-c for the' iteration index [k). L [k+1)l, p ,p [k] g [k] -(2.7.11a) j L,[k+1]. F,.-g q[k]. 3 q[k] (2.7.11b) L k+1) = F, 4 + T L [k] + U L [k] (2.7.11c) z X y Alternatively, in matrix operator notation we get Eq. 2.7-12 2-34 P t

--_7_-_-_-___-_-_ z. ? g L[k+1]' p j - --[k] ~ '~ 0 P :Q L x x x y f L, (2.7-l') F 4- +- R 0 S L = t.2-F, 4 .T U 0 L, ' The matrix' operator of'Eq. 2.7-12 is really the basic matrix operator needed to study all of the iterative schemes that we will consider. Let 0 P. Q . Klbe defined as K = R0S Then one can hope for a convergent iterative TU0l process'if and only.if the spectral radius of the matrix operator K is less than ~ ' unity'-(p(K) < 1). ' The corresponding Gauss-Seidel scheme for the solution of the leakage equations 'is,given as Eq. 2.7-13a-c p t [k] g g [k]- (2.7-13a). L [k+1], px, x Ek+I) + SL (2.7-13b) L Ek+13 = F e + RL j y x z T l [k+1] y ( [k+1]- (2.7-13c)- L [k+1].= F 4 + 7 2 x ~ However, the corresponding Gauss-Seidel matrix operators given as Eq. -2.7-14 is somewhat more complicated.

ll - K') L[k+1],

p4, gn t[k] (2.7-14) 0 0 0 0 P Q 0 0 S R 0 0 and K" where-K' = = T U 0 0'O 0 2-35

.31_f i j Thus, the Gauss-Seidel operator G5 implicit in Eq. 2.7-14 can be expressed as Eq. 2.7-15. i GS' = ( 1 - K' )*I K" (2.7-15) - Again-theoretical?y, we need the spectral radius (p(GS)) to be less than unity.. -for convergence and fully expect that p(GS) < p(K) for a better asymptotic rate. The. explicit form of the operator GS' can be gleaned by substituting Eq. 2.7-13a - into' Eq. 2.7-13b and then substituting both of them into Eq. 2.7-13c, After . rearrangement this.sigebra yields a new system for the GS operator Eq. 2.7-16. l t-'[k+1) px,-' L -' [k] x X L =[I + H) F 4 + GS L (2.7-16) y y L F, 4 L, 2 where 0 -0 0 -0 P Q H-= R 00 and GS = 0 RP RQ+5 T+0R U -0 0 TP+0RP TQ+URQ+US It is then natural to assess the error vector for both the Jacobi and Gauss- - Seidel-methodologies given the same initial guess leakage vector L(*). Since we have assumed that the spectral radius (p(K) < 1) to be convergent ~, we know that ~ formally the inverse leakage operator is given below by Eq. 2.7-17 a von Neuman series expansion q '[ ! - K )~1 = 1 -+ K + K2+.K3 +... (2.7-17) 2-36 l

Consequently, the true soluti_on L[true) is given by Eq. 2.7-18 ^ 2 3 L[true), p , gp, g F+ K F+ (2.7-18)- Assuming now that the starting vector is the same, let us compare the results for both one, e(J-1), and two, e(J-2), Jacobi iterates with the result for one Gauss-Seidel iterate, e(GS-1), based on the relative dif ferences f rom Eq. 2.7-18 explicitly keeping all term up through the third order. The resultant algebra yields the following Eq. 2.7-19.a-c e(J-1)'= K (F - LE*)) 2 3 (2.7-19a)- + K F + K F 2 (p, t[*)) 3 ~ , gF (2.7-19b)- e(J-2) = K e(GS-1) = (K - H ) F - GS LE*3 + K F + K F (2.7-19c) 2 3 2.7.2.4 ' Implicit Leakaoe Treatment Schemes An alternate approach for updating the net leakages has been proposed by Trapp (14) and the computational complexity analysis later cone by his associate Kime @). This approacn involves using a truncated Neumann series expansion as an approximation theory tool for the net leakage equations known as the BALUNS approximation - Brookhaven Approximation to Leakage Updating with Neumann Series. _Thus, the exact expression given in the previous section as Eq. 2.7-17 when truncated at any finite integer N-gives an Nth order appro_ximation to the sought. inverse. Again, the pivotal mathematical necessity is that the spectral radius'of the operator K be convergent, i.e. p(K) < 1. This is the initial assumption for the' justification of the Jacobi iterative technique. This allows the computation of the leakages L, L, and L to be considered just x y Z an intermediate quantity as one is computing the new leakage term of the source for the right hand side of the flux equation. Eq. 2.7-20 given below specifies for an unknown N the general form of this approximation theoretic tool. 2-37

0 M 1 '-' [new) p,-- t X X L ={1+K +K2 +... +' KN} F,, (2.7-20) y l F 4 z 7 From a consideration of computational complexity only the zeroth and-first order' - approximations will be seriously analyzed. Higher order approximations have f ar to high an operation' count to be seriously considered as a programming alter-native. However, it is still possible to use the theoretical results from the higher order approximations to enlighten our understanding of the BALUNS pro-cess. The zeroth and first order BALUNS schemes are given as Eq. 2.7-21 and 2.7-22. L -'[new) 1 0 0-' F e ~' X X L 0 y - 1 0 F 4 (2 7-21) = y L 0 0 1 F 3 z' and t -'[new) i p. gT p 4 X X L '= R 1 S F (2.7-22) y y y L T U 1 F 7 7 Given the' true theoretical salution-d the problem from Eq. 2.7-18, it is tri-vial to write the error. in th0 BALUh3 approximation for any given order N as e(B-N). Then the remaining higher order terms of the infinite series in the operator K'as given in Eq. 2.7-23. 2-38

I 4 K +2 p..,.. for N = 0,1.... (2.7-23) N e(B-N) =K*I F + N 3; ~ Before proceeding with the leakage treatment numerical analysis procedure as implemented inl ARROTTA, we must take note of an interesting equations result invol_ving the use of the BALUNS approximation in the flux equation. Recall that - we have used two forms of the neutron diffusion equations _either the first equations of the matrix of Eq. 2.7-6 or the weighted sum 'of the flux and leakage equations given as-Eq. 2.7-7 where Eqs. 2.706 and 2.7-7 are reproduced here for convenience. [W, h h 1, - h h l, h h 1 ) [,, L, L, L7 [ = k M, or S .(2.7-6a) x7 xz y7 x y and M, or S (2.7-D [J. d, J, J 3 b*l ' b ' 'z = x y z x y where J W + _hhF +hhF + hhF = yzx xgy xy7 J + hG _h G = x 7y j7 J hG + h G, = y gx x and-J hG + hG = g yx xy For any arbitrary value of N, an Nth order Neumann approximation when applied to Eq. 2.7-7 exactly yields the same flux equation as an (N+1)st order Neumann approximation applied _to Eq. 2.7-6. The proof of this claim come directly from a substitution. This gives added impetus to use _Eq. 2.7-7 for the total solution algorithms. Because of this result, the following two particular BALUNS approxi- 'i mations impact on the flux Eq. 2. -7 6 should be especially noted - the zeroth and first order approximations. 2-39

The zeroth order BALUNS expression is simply Fu'*l for u=x,y, or z-since u inverse (1. - K) is approximated by the identity 1. When substituted into Eq. 2.7-6a' we get a special lower order approximation called the truncated leakage i . option in ARROTTA. Its accuracy is minimal but its computational compisxity is very low, making it highly attractive for initial iterations until a decent level of'overall iterative convergence has been achieved. Thus, the truncated approximation can be given simply as Eq. 2.7-24 below [ W + h h f.+ h h F +hhF ) [,] = k M, or S ( 2.7-24 ) ygx xgy xy7 The first order BALUNS expression truncated the inverse approximation as (1+K). When the first order BALUNS is substituted into Eq. 2.7-6a we get an approxima-tion denoted in ARROTTA as the implicit leakage option. Its accuracy is linear in the operator K and' compares favorably with the' explicit techniques presented in the previous section. The implicit approximation to the flux equation is pre-sented as Eq. 2.7-25. 1 [W+hhFygx+hhF +hhF,) [9] = [k M, or S] - L' (2.7-25) xzy xy where L' = [J, J, J 1 [Fx,, F v, F $[ x y 7 y 7 t 2.7.2.5 ARROTTA Leakage Computational Strategy r The_ eigenvalue and transient fixec source calculations that comprise the ARROTTA calculational modules must of necessity have very different strategies for the efficient. updating of the leakage source to the flux equations whether or not the explicit directional leakages are being separately updated. This section will_ outline the interaction between the leakage evaluation process and the other iterative processes that directly affect it. 2-40

ARR0TTA eigenvalue calculations currently use a complicated leakage com-putational strategy. 'The three schemes discussed in the previous sections are ordered from one to three as truncated, implicit, and explicit, respectively. Also, the problem convergence as measured by the outer iteration convergence criterion is logarithmically split into two_ intermediate levels of " partial-convergence" relative to the user's final eigenvalue criterion. No matter what the user has requested for the leakage scheme, initially the trun-cated approximation is in force. When the convergence reaches the first inter-mediate level and the user has requested either the implicit or the explicit leakage treatment, ARROTTA switches to using the implicit treatment. Similarly, when the con-vergence reaches the second intermediate level and the user has requested the explicit leakage treatment, Gauss-Seidel iterations are done. In both the implicit and explicit cases, the actual leakage updating procedure only occurs every other fission source outer iteration. Thus, with the ARROTTA static iterative procedure there exists a natural progression of increasing accuracy with regard to the calculation of the leakage source. ARROTTA transient calculations use a very different scheme for inferring the leakage source to the fixed source neutron diffusion equations at each time step of a transient calculation. Once, a particular leakage mode is enabled only leakage calculations in that mode are done. When a time step begins the leakage source from the previous time step is extrapolated by the flux level omega made exponential extrapolation based on the results of the previous two time steps at that given. space point and energy group. In the truncated mode there is no leakage source to approximate. In the implicit mode the extrapolation occurs for the total leakage source at the given point and group. In the explicit mode the extrapolation occurs individually for each directional component of the leakage. In both the implicit and the explicit computational modes the user may choose to have the leakages corrected once during the process of the fixed flux iteration for that timestep. This updating occurs after the first pass of con-current flux iterations has occurred. 2-41

^ w -i t e ~ it' must be noted that there is no current option available to' update the leaka-ges less frequently than once' per timestep. The only other restriction on the f transient leakage option is that the leakage option used to produce the H corresponding initial leakages must contain sufficient data - e.g. the transient ^ can not be 'done in an explicit mode if the initiating eigenvalue calculation was - .done in either the. implicit or truncated mode. '2.7.3 ARR0TTA OUTER ITERATION NUMERICAL ANALYSIS 2~7.3.1 Introduction-4 The neutron eigenvalue problem in ARROTTA is a fission source outer iteration problem. '.This problem occurs during the initialization phase. of a kinetics problem and at each time step in a xenon transient. The basic iterative.proce-dure solves the system for the largest eigenvalue-and its associated eigenvector of all. positive components - the fundamental mode of the flux. The basic process is a power iterative method.- After a description of the power iterations, we will discuss the Chebyshev acceleration technique and then the esf.imation of the parameters needed for optimum implementation of the Chebyshev strategy. t 2.7.3.2 Problem Formulation 7 The matrix formulation in the explicit and implicit modes formally can be repre-sented by the matrix Eqs. 2.7-26 and 2.7-27, respectively J$ =kM& - [J, J, J ] [L, L, L ]T (2.7-26) x [ J + J' ] 4 =kM4 (2.7-27) where -J = W + hhF + hhF + hhF yyx xyy xyz J'.= -[h h l, h h l, h h l ) [ l + x 3 [F, F, F 3T yy xg xy x y z dx= hG + hG gy y3 2-42

1 h G, ') J hG + = y gx x J hG + hG = z yx xj ] .and-k is the reciprocal of the system eigenvalue M is the integrated fission source operator d 4 is the flux eigenvector corresponding to the eigenvalue k -In-both cases one is' interested in obtaining a nontrivial positive component eigen-pair solution to the above equations.. The implicit eigenvalue problem is of th'e classical form {A 4 = k M &} while the explicit eigenvalue problem also contains the extra explicit leakage term. In the implicit case the eigenvalue sought is from the_ operator [J + J']'1 M while in the explicit one it is f rom J"I - M. Given an initial starting pair k[0] and $[0], subsequent pairs k EA3 and EA3 4 are generated at lesst formally by inverting the operators on the left hand sides. of Eqs. 2.7-26 and 2.7-27. These inverse are not actually computed but are approximated by the inner iteration process that is described elsewhere in the text. For the purposes of the exposition, it is sufficient to assume that . this. inversion process'is exact. Computationally, one need only satisfy oneself-that the error comitted by the inner iteration procedure is mathematically small enough not to contaminate the eigentode calculation of the outer process. Thus, the eigenvector estimates generated are respectively given by Eq. -2.7-28 and 2.7-29 for the explicit and implicit formulations and the eigenvalue by Eq. 2.7-30 4Ei+13 = J-l k EA3 M $ *3 - J~1 [J J J,)[y y y] (2.7-28) E p p I ,4El+13 = [ J + J' J-l kEA3.M 4EA3 (2.7-29) kEA*13 = < M $ E A+13 > / < M 4 [ A3 > (2.7-30) where < M. > represents a Euclidean norm based on the fission operator integrated over the entire problem volume. The power method is guaranteed to 2-43

q .y l converge to the largest eigenvalue if it has a single multiplicity,given that the inner iteration has been accomplished without error. 2.7.3.3.. Chebyshev Acceleration Scheme Chebyshev acceleration is a particular exa'mple of a semi-iterative method with respect to the original power method of Eq. 2.7-26 or 2.7-27 (H) and (H). Conceptually, one is seeking an algebraic formula to combine the iterates such that the overall convergence rate is enhanced. Given the two previous iterates-(E13 and (EA'13.and the current unextrapolated iterate $'EA+13 it is desired to find t'wo parameters aEA+13 EA+13 and b such that the linear combination as given .by Eq. 2.7-31 will increase the convergence rate. 4[i+13 4 [t+13 + a[t+1] (4 [t+1],,[t].) + bEi+1) (4EA3 - (E*"I) ) (2.7-31) Notice'imediately that if the parameters a'and b are always zero, then this degenerates into the original iterative strategy, One of the choices for the computation of these " acceleration" parameters is the Chebyshev polynomial for-malism. First we shall explore those properties that have classically led mathe-maticians'to pick Chebyshev theory, and then we will discuss the difficulties to ascertain whether or not the proper estimation of the Chebyshev parameterization i has occurred in ARROTTA. ) The converger.ce rate of the original power methoo basically depends on how well separated the sougnt eigen-Dair (k(1), x(1)) is from all the other eigenstates (k(i), x(t) for i = 2...., N) of the matrix operator. Assuming the eigenvec-- tors form a basis for the vector space we can write the eigenvector estimate at l anyiteration[t)'asacompleteexpansionintermsofthebasisvectors{x(i), i =1,...,N } ' as shown i n Eq. 2.7-32 l N i EA3 = x(1) + {' c(i) x(i) (2.7-32) ) x .i =2 i where the c(i) are scalars and the eigenvectors have been ordered such that x(1) is the desired fundamental mode. Similarly, the absolute error after [1] itera-2-44

-e ' n'; e -I l l tionsLis--just the sumation of _ the higher order modes as given by Eq. 2.7-33. I N. { c(i)x(i)- (2.7-33) } ~ i=2 e-.=- L 'Doing one more iteration, the 1+1st, the expansion of.xEI+1) can now be expressed y Eq. 2.7-34 y N EA+13='{ktrue (1)/kEI)}x(1)+((ktrue(i)/kE13l_c(i)x(i) (2.7-34)- x i=2 ,kbl*13,... EA) If we now assume that the sequence of eigenvalue estimators k is sufficiently close to the true fundamental eigenvalue k[true] then after r more power iterations the expansion Eo 2.7-34 can be approximated by Eq. 2.7-35 N x[t+r]=x(1)+{-{ktrue(i)/ktrue(3)}r c(i) x(i) (2.7-35) i=2 <f .Thus, yielding Eq. 2.7-36 for' the error expansf on - N' ,[t+r],.{ {ktrue(i)/ktrue(1jgr c(i) x(i) (2.7-36) i=2 However, the eigenvalues have been assumed ordered such that l _ k(i) l. for i > 1 (2.7-37) 'k(1) h[ The dominance ratio d is defined as the c3ximum over all higher modes given by Eq. 2.7-37. Therefore, after r power iterations, every error mode is multiplied j by a number no larger then d r61 sed to the r power. Those eigenmodes that are [fo smaller have been better annihilated than the largest mode which has only been 2-45 Ei-In

reduced by this factor. This is basically a sketch of the elemen s that show that the power method eventually converges albeit slowly, j The r iterations that started with the expansion x(* of Eq. 2.7-32 can be viewed as repeatedly applying the original matrix operator of Eq. 2.7-26 or 2.7-27 that was normalized by the true system fundamental ei0envalue. The natural mathematical extension is to ask is there an rth degree matrix polyno-mial operator whose needed properties would include Eq. 2.7-38a-b7 Q (1) = 1 (2.7-38a) p [l k(i)l / k(1) 3r (2.7-38b) Or [k(i)/k(1)] The efficiency condition of Eq. 2.7-38b is the crux of the acceleration problem.- Assuming that we know both a maximum, d - the dominance ratio, and a minimum-ratio c Eq. 2.7-38b need only be applied over the interval 0 < c (l k(i) l / k(1) < d < 1 for i > 1. The ARR0TTA code actually allows this eigenvalue ratio to lie'in an' ellipse pro-vided that-the real axis is the major axis rather than the minor. axis. However, that generalization will not be covered here. The polynomial of degree r having the least maximum modulus'over the range [c,d] of the real line that also satisfies that boundary condition Eq. 2.7-38a is pre-cisely the normalized Chebyshev polynomial as given by Eq. 2.7-39 T [(2y-d-c)/(d-c)) p P (y) for r > 0 (2.7-39) = p Tr [(2-d-c)/(d-c)] where cos -( r arccos w) if l wl < 1 T (w) r cosh ( r arccosh w) otherwise = l 2-46

p w ,s V p Q s

l 'g Chebyshev polynomials-are trivial to generate since the following initial' con-h' (diti,ons' 2.7-40a-b and the recurrence relation 2.7440c well define the sequence

~computationally. T (w)' = 1 (2.7-40a): O .T (w) = w . (2.7-40b)- 1 and-T (w) =-2~w T r-1(w) - T r-2(w) for r > 1 (2.7-40c) r 'Using the-normaliz' ed Chebyshev= polynomials P (y) af ter completing an initial-[1] r power iterations ofLthe original system and assuming that we know the parameters d and c, it.is easy, to generate a sequence of parameters to use in Eq. 2.7 p tfor.the Chebyshev sequence as given below in Eq. 2.7-41 -) E A+11 : = 2 / ( 2 -'d - c ) a l 4 bel +1E=~0 (2.7-41) a[g41) =' 4 T,1 (q) / ( (d - c) T (q) ) for r > 1 j p r o cand bE A+1) V JTr-2 (q) / T I9) i r a q . where-1 ( 2 - d - c')- / (d-c) /- ~ .q = Therefore, with a minimal introduction to Chebyshev theory, it has been sketched

thatt the semi-iterative Chebyshev polynomial acceleration gives a f aster asymp-totic rate of convergence for the original set of eigenvalue equations 2.7-26 or 2.7-27.

It is yet' unanswered how to estimate these values and how to decide 'what order of polynomial is proper. These subjects will be approached adaptively in the next section. I 2-47

/ 1 2.7.3.4 Chebyshev Parameter Estimation In what has preceded it has been assumed that we had perfect knowledge of the dominance ratio d., This is, of course, never true in pra'cti,ce. We must be able to estimate a reasonable value for the parameter d (Jl,18). Also, since we do ] not know the dominance ratio exactly, the order of the Chebyshev polynomial cycle to'use is another parameter to be determined. Thus, it is our respon-sibility to define the overall. Chebyshev strategy that is to be employed. . During power iterations,.the dominance ratio can be estimated by examining the - i ratio of the successive norms of residuals as given by Eq. 2.7-42a below ll s *3 - s '~l) ll / ll E'~l) - ([1-2] ll (2.7-42a) E E Q(r.) = or alternatively, elements of just the computed fission source, fs, from suc-cessive iterations could be used as in Eq. 2.7-42b Q(1) = - ll fsEA3 - fsEA'13 ll / ll fsEA~l)- EA~23 ll (2.7-42b) fs t This provides an initial starting estimate d0 for the dominance ratio d. The assumption is made that the minimum value c is zero and that a low order . Chebyshev is best - e.g. third order. Estimates are made af ter every Chebyshev cycle that allow the program to compare the actual observed convergence rate with the theoretical rate that is achievable if that estimate of the don.inance were truly correct. In order to use the Chebyshev polynomials before the actual values of d and c are exactly known, approximate values are used as is given by Eq. 2.7-43 , r EI2 'd -c0)/(d ~CO)) T 7 0 0 for e > 0 (2.7' M Tp [(2-d ~CO)/(d -c0)3 0 0 If the iterate x[t+r] can be approximated by simply the fundamental and the second mode, which is inherently true eventually in the iterative process, then 2-48

f Fjt p I x[1+r] can be given as Eq. 2.7-44 . e ii x[1+r) = P(1;r.d ) x(1) + P(d;r,d ) c(2) x(2) (2.7-44)~ 0 O s Since.the Chebyshev polynomial P(1;r,d ) = 1, the residual vector can be O approximated by Eq. 2.7-45 i .l,[1+r+1), 4[t+r) [ = l (d - 1) P(d;r d ) c(2) x(2) l (2.7-45) 0 ~ and'the residual ~ quotient becomes expressible by the ratio Eq. 2.7-46 Q(1+r))=j P(d;r.d ) / P(d;r-1,d )l 2.7-46) O 0 However,'. based _on the' initial condition 2.7-40a, P(y;0,d ) = 0 for 'alI y. Tnis 0 a causes the total ~ residual for the r iterations of the Chebyshev cycle given by-l

the' product of the r residuals of the form'of Eq.- 2.7-46 for the iterations 1+1 to 1+rLto be given.by the simple expression Eq.

2.7-47 l Q[r] = ll P (d; r, d )! 52' '4 I O .Therefore, Eq. 2.7-47 is a computable measure of the dominance ratio that can .be calculated.by' solving the expression for the value of d. It is now possible ] .to compare the computed residual of Eq. 2.7-47 with P(d ;r,d ) and make:an O O l. intelligent. decision of how to proceed. o 3 E > P (d ; r, d ) Case 1:- 1>Q O 0 .The actual value d is greater than the estimated value d. Consequently, the 0 L theoretical convergence acceleration is not being obtained. Given Eqs. 2.7-43 i and:2.7-46, Eq. 2.7-47 comes from simple substitution Eq. 2.7-48 l 2-49 L l 1 ~

T [(2d-d 'CO)/(d 'CO) I 0(7) O O p for-r > 0 (2.7-48)- =- t r [(2-d 'C )/Id 'C )) T O O O O Using P(d ;r,d ) =.1 / T [ (2 - d0~CO) / (d0~CO) ) implies that 0 O r E3 T [ (2_- do~C)/ Ido*C))=Q / P(d II'd ) r O O O O It is now possible to use Eq. 2.7-49 in order to obtain the new estimate for the dominance ratio d. d=[d'C3{ cosh-(1 [0[r] 3 [d 'C 3 O 0 O 0 P(d ;r,d ) ) + [d +c 3 } (2.7-49)' arctosh 2-r O O O 0

it

lf the minimum of c is replaced by the estimate zero, then a slightly simpler 0

form for:the new dominance is given by Eq.- 2.7-50 d= { cosh-(1 arccosh E ) +1} (2.7-50) '2 r P(d ;r,d 0 0 Case 2: Q[r]< P (d ; r, d ) (# 1) O O Thi s' t 3 case implies that the true dominance ratio d is less than d. ' Consequent-0 ly,-tne formulae given by Eqs. 2.7-49 and 2.7-50 above are modified to use the trigonometric cosine and inverse cosine rather than the hyperbolic ones, respec-tively.; 1: t Case 3: 1<QEP3 [ In this predicament there is no error reduction. Some of our assumptions must be L invalid. The current strategy within ARROTTA realizes that the adaptive process based on the Chebyshev residuals has broken down. Power iterations are restarted until a robust estimator is again gleaned from the process. 2-50 ^ 1s

q a t t 2.7.4-ARROTTA INNER ITr. RATION ALGORITHM . 2.7.4.1 LIntroduction ~ .The. inner iteration problem is given symbolically by Eq. 2.7-51 where F repre-sents the particular form of the two group neutron diffusion operator currently being solved, & is the unknown two group flux vector, and q represents the fixed' source. F4=q (2.7-51) In ARROTTA the operator; F is a block matrix partitioned by plane and row so that the ' smallest iterative unit that is considered is an entire line for both neutron. groups simultaneously by use of the natural integer ordering rather than-a P-cyclic ordering such as the " red-black". It is entirely possible to have a different block size. However, the mathematical ordering of the unknowns when .i taken together with the ARROTTA scratch disk data structure as discussed elsewhere in this report, make this block size a natural. It should also be. noted that any block size smaller than this will increase the demands on the iterative process while lowering the asymptotic estimates of convergence rates. Thus, since there is no difficulty in the direct line. inversion for this blocking containing a maximum of 40 tridiagonal elements each of which is a 2 x 2 matrix, no other direct inversion blocking structure was seriously considered. This is done by especially written block forward elimination back substitution (BFEBS) tnat utilizes the known structure of the block diagonal component-of F, denoted as D. F is now denoted by the partition given as Eq. 2.7-52 where L 'and V are respectively block lower and block upper subsets of the original operator F. F=D-L-U (2.7-52) 2-51

A l 2.7'.4.2 Iterative Methodology It is important-to note that this matrix structure leads one to first think of I the Jacobi iterative technique as symbolically given by Eq. 2.7-53. J= D'l (' L + U) (2.7-53) - However, the two classic reasons for avoiding the Jacobi technique are again valid _here: first,.a slow asymptotic rate of convergence, and second, an unne-cessarily-high storage requirement. i : Initially, the' inner iterations are solved by using a classic Gauss-Seidel pro-cedure given in operator form by Eq. 2.7-54 (1 - D-1 L )~1 D'l U -(2.7-54) G = where I is the proper size identity matrix. Given alternatively in Eq. 2.7-55 in terms of the iterative unknowns {P(n) : n= 1,...,N}forNthenumberofY-nodestimes2-nodesandthesuperscriptkisthe iteration index. Note th6t in-the form of Eq. 2.7-55 the blocking is nota-tionally simply the direct inversion-of D. D P(k+1]., y p[k+1], g p[k] + q (2.7-55) It is well' known that in the symmetrizable positive definite case, the spectral radius of G can be estimated from the residual error norms that are available during the course of iterating. The PDQ (J9_) series of finite difference programs utilizes this feature to calculate a priori spectral radii estimators for.each neutron group prior to the actual beginning of the iterative process. l However, it-is also possible by adaptive techniques to dynamically estimate the spectral radius during the iterative process itself. This is rigorously true for the positive definite symmetrizable case alluded to above. Also, experimen-tally it seems quite reasonable to use the same error residual indicators in the nonsymmetrizable case to ' guesstimate' such parameters bearing in mind that if 2-52 m.___

t ,+ i u

c i

4 . these. estimators seem unreasonable then one must restart the estimating proce- ' dure while continuing to iterate. Whenever a firm estimator of the spectral radius of the process exists, one can then infer an ' optimum'. value of an 7 overrelaxation parameter w(opt) that can be used until'a new estimato becomes available. This successive overrelaxation procedure (SOR) yields a matrix operator of the 'orm Eq. 2.7-56 with the needed property that when w = 1.0, Eq. 2.7-56' degenerates into Eq. 2.7-54 as expected. S =-(1 - w 0*l) -~1 (w D'l V4 (1-w) 1 ) (2.7-56) 1 . and-the constant vector q is replaced by h as given in Eq. 2.7-57 h * (1 - w D'I L)'1 (2.7-57) 4 When given in explicit form of the unknowns we have Eq. 2.7-58 for the iterative' index [k]. ([k+1], 3.,(k] + h (2.7-58) i 2.7.4.3 Sor Parameter Estimation 'It.is only necessary that we are able to estimate.the spectral radius of the associated Jacobi iterative matrix J (p(J)).of Eq. 2.7-53 or alternatively the spectral radius of the Gauss-Seidel, operator G (sqrt(p(G)) = p(J)) of Eq. 2.7-54'in order to proceed with an acceptable value for w in the S operator. Given p(G), w(opt) is calculable from the well known relationship given as Eq. 2.7-59.. sqrt (1 - p(G)) 3 (2.7-59) L w(opt) =.2/[1 + It is possible to estimate the required spectral radius either by'doing a priori iterations to determine it or while doing less than optimal iterations to con-tinue-to update the estimators of the spectral radius. The first scheme con-sists of. solving the Gauss-Seidel operator of Eq. 2.7-54 when the valuation of l 2-53

i the constant vector q is zero everywhere and.the initial guess for the unknown vector is~ an arbitrary nonzero vector. We know that the true solution of this degenerate system of equations is the zero vector. However, the rate at which the iterative solution procedes toward the zero vector is a reasonable measure ~ i , of the spectral radius of Gauss-Seidel operator G as expressed by Eq. 2.7-60 for iteration k for a suitable Euclidean norm, 3 ~3 (G) = ll x ll / ll x ll (2.7-60) p where x[k] is the pseudo flux of the degenerate iterative problem. I i However, the time spent on computing the spectral radius this way does not j directly help the computation of any problem of interest. _ Consequently, more i interest has been shown in those methods that are truly adaptive (20). Although j such mathematical parameter estimation formally requires the matrix operator to be symmetric positive definite, the ARROTTA operator seems to produce fine j iterative acceleration behavior when each " point" of the matrix is the 2 by 2 nonsymmetric matrix operator derived from the two group neutron diffusion equation. It is obvious that the Gauss-Seidel operator _ is just the proper - starting guess for the optimized SOR operator S of Eq. 2.7-56. The iterations proceed with u = 1.0 as the initial w value. Let the iterative difference resi- - dual as measured in a convenient norm be given by Eq. 2.7-61. l J r[k] , g g x[k] x[k-1) g j f g j x[k-1), x[k-23 ll (2.7-61) Let m be the number of iterations since we last changed parameter estimations. Then it is time to again estimate the spectral radius of the iterative operator if. and only if the following criteria are all satisfied: m > NUMIN = 5 4 i 2-54

[ w ,.(, 3 i e i i s and d ,4 m ( 1 - w)**I < Psp = 0.5 1 f. and .-10*R'spl<r[k-'1)-r[k]<Rsp=.0001 h 4 . where the def ault values for the parameters are given above. If the previous group'of iterations were Gauss-Seidel, then the spectral radius _ estimate comes .'directly f rom Eq. 2.7-57. However, if we have already heen iterating with a value of the overrelaxation parameter w that is less than the optimum, then the l inferred next value in~ the w-sequence is available as Eq. 2.7-62 from Eq. l 2.7-59 based on-the following estimate for the spectrel radius, p = (r[k], g, 3) / (w

• sqrt (r[k)) )

(2.7-62) i

I In order to ensure that the estimates will always approach w(opt) from below, j

upper bourds *re -imposed:on the ettimates of w as a function of estimation

. to 8.. At each stage w < tES3 where t[sl is given in table 2.1. j number 1 l Table 2.1. Maximum Overrelaxation vs Iteration Index 1 1 2 3 4 5 6' 7 8 >8 4 1.6 1.8 1.9 1.95 1.975-1.985 1.99 1.995 1.995. 1 -i 2.7.5 REACTOR N0h! LINEAR ITERATIONS AND SEARCHES 3 r -4 2.7.5.1 Introduction The-eigenvalue problem sought in ARR0TTA has really a quite ccmplicated set of - dependencies - mathematically much more difficult than the fixea source tran-g sient' evaluation. This is true partially due to the innate difficulty of the- ' feedback mechanisms presented by the thermal hydraulic and xenon interactions -with the neutron flux distribution. However, this is aggravated by allowing 2-55 l A

\\ t i searching of either control bank position or soluble boron level to achieve a desired eigenvalue during the iterative process. Implicitly, it is assumed that j most of the changes in the valuation of the search variable will occur early in the iterative process. Consequently, as the iterative procedure homes in on the desired critical state with a reasonable set of intermediate distributions, one needs to know both what are satisf actory convergtnce tests for the varying para-meters of interest and what procedures are avcilable to enhance the convergence rate of the process. Care must also be exercised in how to terminate this pro-cess in order that the resultant restart filec be well defined relative to the ARROTTA reactor iteration process. Given the current value, the search algorithms will prodne the next estimate of position or concentration based on the current value and the current estimate of the reactivity worth per unit change under the assumption of linearity in the small. This reactivity change estimate implicitly includes the effects of all other feedback mechanisms as well as the present state of convergence rather than the true converged state. Therefore, *.hese estimators are evaluated each time through the nonlinear iterative cycle. No computing effort is wasted by this since the feedback evaluations have already forced the pointwise cross sec-tions to be regenerated. The exact mathematical nature of the nonlinear feedback problem is certainly not fully understood by exact analysis. However, it has been observed that an adap-tively generated one term acceleration similar in form to the successive overre-laxation (50R) formalism can improve the observed convergence rate once the gross inadequacies of the original input guess h6ve been removed. The procedure as given in ARROTTA contains indicators that tell the program when and how to update this parameter if the necessary indicators are robust. Otherwise, no ~ acceleration is attempted until these indicators again produce reasonable esti-mates of the nonlinear process. 2-55

i i -2.7.5.2 Search Procedures The eigenvalue solution routine in. ARROTTA will allow either the position in either a BWR or a PWR of a single control rod group defined arbitrarily by the user or the soluble boron poison in a PWR to very to achieve criticality - e.g. t A'= 1.0 where A' is the inverse of k-ef fective. Let n be the reactor iteration index. Then given the current value of the search variable q[n-1), the present iterate of the eigenvalue [n-1), and the current estimate of the reactivity change per unit change of the search variable (d/dq)[n-1), the next value of the search variable q[n] can be inferred f rom Eq. 2.7-63 as given below since that relationship can be approximated by a simple linear expansion. A'. A[n-1) {q[n], q[n-1)} (dA/dq)E"'I) (2.7-63) This yields the next estimate of the search variable as given by Eq. 2.7-64 {x,,x[n-1)} /(dA/dq)E"*1) (2.7-64) q[n] q[n-1) eigenvalue A "3 that will allow E The iterative procedure will now produce a ner a the reactivity derivative to be updated by Eq. 2.7-65. (dA/dq)E"3 = (A "3 - A "~l) ) /-( q "3 - q "*13 ) (2.7-65) E E E E This procedure assumes that an estimate for the initial reactivity derivative is known, if this'information is known to the user it should be input as the RHSHGS (rho Search Guess) variable in the real input control series. While it is possible to have such-information it is not necessary. If no derivative infor-mation exists on the first search variable modification, an arbitrary ten per-cent perturbation will be made. The sign of the change is determined by the-initial calculated eigenvalue [1]. The only specific modifications to the above rules are that the initial boron change must be at least 5 PPM. Also, at all times in the searching the control group rod position must be within the problem definition or the boron poicon must be nonnegative. 2-57 i 6 ~

2.7.5.3 Reactor Iteration Procedures The reactor iteration eigenvalue problem addresses the entire functional depen-dency typified by the circular logic of the nonlinearity; i.e. a cross section distribution implies a power distribution which in turn implies a feedback variable distribution which in turn implies a new cross section distribution. It is implicitly assumed that the convergence of this total process can be assessed by just monitoring the convergence of the pointwise power distribution. The actual ARROTTA reactor iteration ptoceeds in the following circular manner af ter an initial set of conditions are assumed for the starting guess: l 1. infer spatial xenon and samarium poisoning unless xenon tr:3.sient 2. generate pointwise cross sections 3. generate flux distribution 4 infer power distribution 5. check for power distribution convergence 6. calculate thermal hydraulic feedback variables 7. exit if converged or iteration limit is reached 8. Update search variable if needed The reactor iteration is considered converged when the largest power difference at all of the power producing nodes of the thermal-hydraulic mesh has converged when the criteria given as Eq. 2.7-66 has been reduced below a prescribed input value. cEh3=MaximumP)0lP")-PEh*I)l / PEh) E (2.7 66) nodes It should be noted that in order to have consistency among all possible problems that could be executed including restarts, one additional thermal-hydraulics calculation must be executed even when it is known that the iterative process is about to terminate. 1 2-58

I 2.7.5.4 Acceleration of The Reactor iteration Procedures Each pass or entrance to the thermal hydraulic routines of ARROTTA uses the { current thermal hydraulic mesh power to generate a new calculation of the feed-I back variables, it is reasonable, therefore, to ask the question: is there a way to combine the results of the last power calculation with one or more previous I ones in order to accelerate the convergence process ? Thtse operators are f ar too complicated to attempt to find a closed form mathematict) solution to this problem. However, it is possible to give a reasonable anthe.atical form to the problem by mocking the situation that was posed in the adaptivo ton.putation of the successive overrelaxation parameter w as was doQ in MEKibD Q}), Let the general reactor iteration operator be denoted by fi Then the conven-tional SOR scheme is schematically expressed as Eq. ?.7-67 pin], p(n-1].,(n) {p p(n-1), p[n-1]} (2.7-67) Obviously when w ") = 1, this operator degenerates into the ur'lginsi power l iteration. Assuming that SOR theory is operable, it is nectraary to find a value of w that minimizes the spectral radius of the total Operhtor. The spectral radius of this process is assumed to be given by the L2 norft (square root of the sum of the squares) of the iteration nodal power difference:. as Oven by Eq. 2.7-6B x(n],jgp[n].p(n-1]jj j g g p(n-1), p(n-2]ll (2.7-68) Then the new w(n) is given by Eq. 2.7-69 sqrt(1-u"))} (2.7-69) E wE")=2 /{1 + where u(n), {x(n). y[n 1} / {sqrt (A ") )

• wEh'I)}

E The above criteria are only considered valid when the estimator for the spectral radius u ") has stabilized to a robust estimator as defined in relationships of b Eq. 2.7-70 2-59 1 y

i r l t l l (0.1 - ) e "3 > l xl"3 - xl"'13{ / xE"3 (2.7 70a) and 1 > xC"3 > l (e "3 - 1)2 (2.7-700) if the above are not valid, tnen ARRoTTA just continues with the current set of parameters. The only exception is that if the estimator of the pseudospectral radius of Eq. 2.7-6B is greater than unity, then the knowledge of the operator is highly suspect. Consequently, the SOR parameter w "3 is reduced by thirty E percent but not allowed to become less than unity. I i 2-60 1

i F u i s Section 3 [ t FLUID DYNAMIC MODEL i i 3.1 ' INTRODUCTION 'The fluid dyn6mic model of ARROTTA can be characterized as an inhomogeneous. non-equilibrium, two-phase flow model based on the concept of energy split 1 between the liquid and vapor phases (22). Basically, it is a four-equation mooel: one mass balance equation for the liquid-vapor mixture, two energy balance equations for the liquid and vapor phases, and an algebraic slip rela-tionship.- The momentum balance is not considered because it is decoupled from the 1ergy equations and because the inlet mass flux is specified for the core y as a boundary condition. 1 As'we shall see later, the three original conservation equations are transformed into'two separate mass balance equations for the liquid and vapor phases involving a vapor generation term and one erargy balance equation for the liquid ' phase. Tnis final set of equations is actually solved in ARROTTA. l. The fluid flow in a BWR is a collection of single-channel parallel flows. In a PWR this can be assumed as a first approximation. Hence for both types of reac-tors we consider the conservation equations reduced to one dimension, in what follows, we shall start with a general set of conservation equations for a one-i dimensional two-phase flow, including the momentum balance. We shall then show which terms 'are neglected under what cor.ditions, leading to the set of conser-i vation equations that are actually solved in ARROTTA. L / The general conservation equations for a constant cross section, one-dimensional two-phase flow are of the following form: 3-1 I+ M

7, m h# ~ L, l t j h Conservation of Mixture Mass i i + [(1-a)pg + ap 3 + [(1-a)p Yt t * "#g#g) = 0 (3.1-18) l Conservation of Vapor Energy (' (ap H ) + (ap H v ) - a +q ( 3.1-lb ) = g ggg Conservation of Liouid Enercy 3 [(1-a)p H ) + [(1-a)p H v ) - (1 9 ) =qyg + q,g (3.1-Ic) gg ggg Conservation of Mixture Momentum I [(l*a)p v.+ op v 3 + h [(1-a)c v vg g g + ap v v 3 ( 3.1-ld ) gg gg ggg P .+ [(1-a)og + se ) g + F, = 0 + g in the above equations, g is the volumetric heat source per unit volume per i yg Unit time for tr.e vapor phase and g is-that of the liquid phase. These volu-yg met.'ic sources 'are the direct (prompt) moderator heating due to gamma ray atte. nuation and neutron slowing down in the coolant. 9 The net convective ne6t source per unit volume per unit time from clad wall to tne coolant is represented by q for vapor and by q,g for liquid. These sur-gg fece heet sources depend on the heat transfer mechanisms considered in each of the flow regimes assumed. -The wall-fraction loss is represented by F, which depends on a friction factor usually given by some empirical correlation. L 3-2 i

.s-g i r i t . The terms involving the pressure P represent the energy production due to the ~ rate of pressure chenge and the pressure gradient (note that D/Dt = 8/8t + vat'). The current version of ARROTTA assumes that the pressure is constant in space and that the rate of pressure change is sufficiently small so that 9P/8t is negligible relative to the other terms in the energy equations (2,2,). Consequently, $=0 (3.1 2) Dt and the momentum balance equation becomes decoupled from the mass and energy balance equations. However, the system pressure is specified as a function of time and the thermodynamic properties are reevaluated at each time step with the new pressure. Neglecting DP/Dt.is a good approximation for normal operating conditions and for

most anticipated transients including turbine trip events in BWRs if the rate of pressure increase is less than 100 psi /sec. For a rapidly depressurizing tran-sient such as 1.0CA, eq. (3.1-2) may not be acceptable.

The usual boundary condition in a light water reactor is that the pressure drops across all flow channels be equal. This calls for pressure drop calculations 'from the momentum balance equation; the calculated pressure drops are'then used to adjust the inlet mass flow distribution to achieve the equal pressure drop ~ boundary condition. ARROTTA prescribes the inlet mass flow distribution as a function of flow channel by an input array derived from auxiliary calculations external to ARROTTA. L lt is important to note that in writing down the conservation equations (Eqs. 3.1-1) the following assumptions have been implicitly made: (i) Viscous dissipation is negligible. I (ii) Mechanical energy exchange between phases is negligible. (iii).Covariance terms are negligible. 3-3 1 1

f" I i In the conservation equations, a is defined as the area fraction for vapor com-monly called the void fraction. The vapor density.p is uniform over cross-g sectional area. Finally, Hg is the density weighted area-averaged enthalpy of vapor.. All the other symbols are defined in the Nomenclature. .) 1 3.2 TWO-PHASE FLOW EQUATIONS ) For ease of writing, let us define the following quantities: 1 F n (1-a)p ( 3.2-la ) g g ) F E ao (3.2-lb) ] g i G n (1-a)p v ( 3.2-Ic ) g gg G ( 3.2-1d )- l g aaegvg We shall ignore the momentum balance equation as discussed before. The mass and energy balance equations (Eqs. 3.1-la, 3.1.lb and 3.1.1c) become: h(F+F)+h(G+G)=0 (3.2-2a) g 9 g g h(FH)+ (G H ) = q +q (3.2-2b) f gg gg h(IH)+h(GH)=q ( 3. 2-2c ) tg gg +q The volumetric heat sources are: 4yg = 0 (3.2-3) i 9vi = Y 02f (3.2-4) 4 3-4 f

~ ( t is the f raction of total. fission energy that is' directly deposited in where Y2 the coolant (moderator) and O is the total heat source per unit volume per unit f time due to fission. I Five flow regimes are considered in ARROTTA. They are defined as follows: 1. Single. phase liquid-flow, Tc2 < T3 2. Single-phase liquid flow, Tc2 > T3 3. Two-phase flow Tc2 > T ; Tg<T3 (subcooled boiling) D .4 Two-phase flow, Tc2 > T ; Tt =T3 (bulk boiling) 3 5. Single-phase vapor flow. Tg>T 3 where. Tc2 = clad outer surf ace (wall) temperature, T = bulk coolant (liquid phase) temperature, g i T = vapor temperature, g T saturation temperature at system pressure, a g i critical clad wall temperature at point of void detachment. T = D We now discuss the conservation equations and the heat transfer mechanisms in each of the five flow regimes. All correlations to be presented in the following are expressed in British units except for time which is _in seconds rather than hours. 3-5

rg,. J 6 i l I.' 3.2.1 Single-Phase Liquid Flow: Tc2 < Ts The fluid dynamic equations (Eqs. 3.2-2a through 3.2-2c), upon setting a = 0 -reduce to o-L t nF, + 301=0 (3.2-5) r nt 42 r I h (Fg g) + (G H ) = Y 02f+h H (3.2-6) [ tg p l' where the surf ace heat source d[l is represented in terms of the Dittus-Boelter heat transfer coefficient hDB IE) l d((=D h Tc2 - T ) (3.2-7) DB g eq k 0 0.4 DB = 0.023 Re.8 pp g t) (3.2-B) h eq Here D is the equivalent diameter of the coolant channel, Re is the liquid eq a Reynolds number, Pr the liquid Prandt1 number, and k is the liquid thermal con-g ductivity., They are given by the following correlations and definitions: G Re = D,q fF. ( ) (3.2-9) t k = 1.60497x10'# (ak0 + akl *k + ak2 *k + ak3 *k) (3.2-10)- g Pr = Cpg ug/kg (3.2-11) 3-6 r

i I where = 4.0 x 10*3 H (3'2*l2} 'xg L a0 = 0.57373B622 (3.2 13a) k a = 0.2536103551 (3.2-13b) g3 e ak2 = -0.14546B269 ( 3.2-13c ) a3 = 0.01387472485 (3.2 13c) k 445 (T +208) = 1.622208 x 10-5 10 (1+4.009409x10-8 (3.2-14) v g ~ ~ i (P-P,g) (T -89) 3 gat = 3206.2'- 6483.29 (705.4-T ) [1-2.9x10-4 (3.2 15) P g (705.4 T)]/(1012-T ) g The coefficient f is a f actor which adjusts the nominal core flow rate (input p value) to account for the effect of equal pressure drop but unequal power pro-duction in each thermal-hydraulic channel and must be specified by the user. The liquid specific heat Cp will be defined in the following discussion. The liquid (water) temperature (T ) in this flow regime is determined from the g liquid enthalpy (H ), which is calculated from the fluid dynamic equations (Eqs. g 3.2-5 and 3.2-6) by' a Newton-Raphson iteration scheme: (") .T (n+1) = T ("I + s - a 1 1 ( 3.2-16 ) i i C (") PL b H(n)=b)+b T (n), 3 b4 g 2 g (802-T (D)) (T (n)-750) 1 g g 3-7

p C (3) = b + (T (")-750)# (3.2-18) EA (802-T (")) t g where the above correlations for H and C are taken from Jordan Q4) g pg bg = - 47.62 + 0.01618P (3.2-196) 2 = 0.96566 - 4.984x10-6P (3.2-19b) b l i b3 = - 8843.4 4.8.2445P (3.2-190) b4 = - 20392. + 6.6112P (3.2-19d) and P.is the system pressure in psia. The superscript (n) is the iteration index. The iteration is continued until the following criterion is satisfied in each thermal hydraulic region: T ("*II - T (n) 1 i < 10*7 (3.2 20) T ("*I) - Tg+1 g where Tgn is the inlet water temperature. Note that there are three unknowns in Eqs. 3.2-4a and 3.2-4b; namely, G, H and g g F. The additional equation needed to solve for the three unknowns is supplied g by the equation of state describing the liquid density o for the pressurized g water (24). (Note that F, = p in this flow regime): g g = 92.924 + 5.761x10*4 P+ (3.2-21) o (39440.2 + 1.6386P)/(Hg - 1377.35 - 0.035704P) 3-8 i i

7 3.2.2 Single-Phase Liquid Flow: Tc2 > Ts In this flow regime. the dynamic equations Ecs, 3.2 5 and 3.2-6 still apply, but i-the surface heat source is now given by f c2-T ) + hg (Tc2-T,)) (3.2-22)' a : * (D ) IhDB (T s g eq where B = 0.05358 e /630 (Tc2-T ) (3.2-23)- P h 3 l-which is the Thom (,25,) version of the Jens Lottes correlation (26,). g h in Eq. 2 3.2-23 is still in BTU /sec-f t ,.7, 3.2.3 Two-Phase Flow (Subcooled Boiling): Tc2 > T, It < Ts D 'This subcooled boiling regime is an important part of the two-phase flow model in ARROTTA. The fluid dynamic equations (Eqs. 3.2-2) apply in this flow regime. However,'some manipulation is done to reduce the dynamic equations so as to f acilitate their solution. In this flow regime TD is the clad temperature at which net voidage-is initiated. Special attention is given to determining where inside the axial mesh subcooled boiling actually starts and linear weighting of heat' transfer coefficients is used in each subregion of the critical mesh zone. When entering saturated boiling, a similar procedure is used to locate the saturated boundary. Equation 3.2-2a can be separated into phasic form as 8Fd + 3G9*T (3.2-248) 8t 82 9 8F gg g + J = -r (3.2-245) 8t az 9 3-9 6 I

s -, ' where 'g is the vapor generation term. If we assume that H does not vary - [ meaningf ully during the transients that are meant to be considered, then Eq. (3.2-20) yields H99 $9 (3.2 25) Adding the energy equations, the total heat source is given by .a. q = Agg + g, (3.2-266) g g and we take (see section 4.4) 1 .u.. q, =-(D4) hDB(Tc2-T,) + h (T 2-T ) ) = qqg + gg (3.2-26b) B g eq where B and NB stand for boiling and nonboiling. We postulate now that '"g = (AB '91)/Hfg (3.2 27) 4 where the condensing heat flux is taken as: h (T - T,) (3.2-2B) qq 0 = j g eq i' and hj is described by the Hancox-Nicoll correlation (27) k' h = 0.4R,0.662. p (3.2-29) k j p eq Hence. Asg = H (qB' ~91)/Nfg (3.2-30) g 3-10

p, -~ b l l-i provides the energy split. The boiling heat flux is derived in Section 4.4 as 3 98 = (34 ) (h g + 7 os),(Tc2 - T,) (3.2-31) h eq r i ~ and the 'nonboiling (convective) heat flux as: l que = (g2)hDB (Tc2 + T, - 2T,) (3.2-30 eq Using Eqs. 3.2-26a, -30 and.31, we have -Hqfg)/Hf (3.2-33) Ast*QNB + (H qq and the liquid energy equation becomes (F,H,) + h (G,H,) = v2 f + h (3.2-34) 0 t In summary,'the fluid dynamics'of the subcooled boiling regime is described by i, f-Eqs. 3.2-24 and 3.2-34. Note that this flow regime has seven unknowns; namely, F,.. F, G,, G, n, H, and H. We need four additional equations to solve for g 9 g i .the seven unknowns. The needed four-auxiliary equations are: (3.2435) H =Hf+Hfg g 8at' (3.2-36) .F =ao 9 9 F, = (1-a) o,, (3.2-37)- a SK'I (3.2-38) =W a Here pt. is given by Eq. 3.2-21, S is the slip ratio defined as the ratio of vapor velocity'to liquid velcity, and Eq. 3.2-38 is based on the Bankoff-Jones correlation (28, 29) for the slip: 3-11

y d** S=- (3.2-39)- where R K =K, + (1-K ) n (3.2-40a) K, = 0.71 +.29 ( ) (3.2-40b) c R = 3.448275 - (1.875x10~4 - 5.85x10*7 )P (3.2 40c) P Pc = critical pressure = 3206.2 psia (3.2-40d) P Properties pertaining to the liquid phase are those of pressurized water as given by Eqs. 3.2-8 through 3.2 21. The vapor properties are described in i Section 3.2.5 (Eqs. 3.2-49 through 3.2-57). 3.2.4 Two Phase Flow (Bulk Boiling): Tc2 > Ts, Te = Ts When the liquid temperature reaches saturation, we have bulk boiling two-phase flow. The fluid dynamic equations, Eqs. 3.2-29, 3.2-31 and 3.2-34, as well as the auxiliary equations Eqs. 3.2-35 through 3.2-40, still apply in this flow regime except that H, = H, T,, = T,.and 7 3 sat F,, = (1-a) o (3.2-41) In this case the surface heat flux source term becomes . f.. i"q+9. *(h)(hg+hDB)(Tc2 - T,) (3.2-42) s eq [ 3-12

i-

3

i..^

l I 3.2.5 Single Phase vapor Flow: T9>T3 l In this flow regime the voic f raction reaches unity. It is assumed that the vapor temperature (T ) is at saturation, and all vapor properties are held g constant at saturation. However, the clad and fuel temperatures may continue to change. The next few equations show that ARROTTA contains nearly all the necessary equipment to calculate superheated vapor characteristics. In the present ver-sion of the code, the actual calculatin of those characteristics has not been implemented. The program currently does not track the two phase / single-phase t interf ace inside a volume, but it does allow such an interf ace to form and either rise or fall in time. The surface heatsource& is given in terms of a heat transfer coefficient h DB h=(D IhDB (Teg-T) (3.2-43) g ec The coefficient hDB employs the Dittus-Boelter correlation for superheated steam, that is, Pr.4 ( 9 ) (3.2-44) 0 DB = 0.023 Re.8 0 h eq G p(9) (3.2-45) Re = D,q f 95 Pr =l.a -b g[exp[-c(T*-T)3-1)l (3.2-46) g g g 3 e F T *

• 1/2 (Tc2 + T ) ; T * > T (3.2-47) g 3

g 3 L l 3-13 l-

2 k = 1.6x10' i (17.6 + 0.00587 T + 1.04x10'4 T - 4.51x10' T3) g ge ge ge 2 + (103.51 + 0.4198 T - 2.771 x 10-5 Tge) ( gs ) gg .2 ( gs )2 ) ( 14 + [ 2.1482x10 T (3.2 48) Tge = (T * - 32)/1.8 (3.2-49) g I o = (3.2 50) [ + 0.5682 exp(5.757x10'# O.9) U * - T,) ] I P g ogat P Here P is the average system pressure in psia and gs = 4.916x10-6 + 1.5195x10-8T * + 3.8x10*7 u p g g3 2 + 1.168x10-8,g + 2.823x10'll 3 o (3.2-51) g3 In the above equations all the properties are those of superheated steam and T

• g represents the average superheated steam temperature.

If P < 130 psia the saturated vapor density is evaluated as follows: olat=[ + 0.0028529 ) P.948 (3.2-52) 0 In Eq.-3.2-46 the coefficients a, b and c are given by g g g g = 1.01 - (.066017 (.0058825 .00046958 y) y) y ( 3.2-53a ) a where y = [(T3 - 150.)/150 32 3-14 i l

I f v i 1.01 - 9.6 x 10-6 T -a (3.2-530) b a g 3 g 3 - 5.43 x 10*3 ( 1 - )3 ( 3.2-Sk ) c = g 3.3 NUMERICAL ALGORITHMS i in this section we shall describe the numerical algorithms employed in ARROTTA for solving the conservation equations discussed in the previous section. The algorithms are developed by means of a conventional finite-difference approxima-tion in space and time. We shall consider only the more general case of two- + i phase flow regimes. The airorithms developed for this general case can be readily reduced to the single-phase flow regimes. i 3.3.1 Transient Algorithms The conservation equations to be solved for the two-phase flow regimes are of the forms: nFY + nG9=" (3.3-1) at mz 9

1. F, + nG, = - t (3.3-2) g h(FH)*

I0 H,,) = 0,, (3.3-3) tt t where

07. " Y 02f+9 (3'3-4)

Consider a vertical flow channel that is discretized into a number of thermal-hydraulic (TH) regions as shown in Figure 3-1. The indices (m-1), m. (m+1) etc., represent the center points of each TH region; whereas, the indices (m-1/2), Im+1/2) etc., represent the interf aces between adjacent TH regions. 3-15

.l \\ I i j k r t i ir - o m+1 0 2 +1 m Gm+1/2 b

m+ 1/2 - -

[dg')/m (7 } 0 ), 9 .g az, 2 (s (S.\\ mi m,H,G f a m m st /m i n u m-l/2 -- S G m-1/2 o m-l AZm-l 1 ____JL-l 'l n Figure 3 1 A Discretired Flow Channel l 3-16 1 ~-~~-' ~ ~ i

i I We wish to find the solution of the conservation equations'in each TH region at [, contiguous discrete time points: t ' t, t ' ' t,... such that t) = 0 l 2 j t),3 + at). Applying a semi-implicit finite dif ference approximation to Eas. 3.3 1, 3.3-2 and 3.3.3 yields-j j j j F,- F,-1 G,,3fp - G,,3fp, {j (3.3-5) g g g g P at)- 62, I j j-1 j j ,Otm+1/2

• Otm-1/2, pj Itm-Fim g4

{ gm r at 62, j L j j~ j-1 j-1 3 3 3 j tm+1/2 tm+1/2 ~ Otm-1/2 gm,3fp, j H H I h

• I H,.G tm tm tm g

3,3_y) D at 42, j where gm

• m i

) (3' *00) ( f9 j j .aij 1 .aij-1 j j O,=Y02 fm + 9stm +Am -r,H, (3.3-8h) L g i g g th Define the following average quantities for the m TH region at time step j: d G = 1/2 (Ggm+1/2 + Ogm-1/2) (3'3*9) gm d d G,* 1/2 (Gim+1/2

• O m-1/2)

(3.3-10) g GfmHd d 3 d d tm = 1/2 (G im-1/2 im-1/2) (3'3*11) H im+1/2 im+1/2 + G H 3-17 e

k.' e I' i These definitions are approximations that are basic to the numerical algorithms developed for ARROTTA. Their validity may be subject to question in the case of flow oscillation, flow reversal or flow stagnation. For this reason, flow { reversal and stagnation are disallowed. Solving Eqs. 3.3 9, 3.3-10 and 3.3 11 [ for the quantities-at the upper interf ace (m+1/2) we obtain J Ogm+1/2

• 20
• O (3.3-12) gm m 1/2 l

1m+1/2=2Gfm*O (3.3-13) 0 m-1/2 Gfm+1/2fm+1/2*20 m" m *ON (3'3'I4) H m-1/2 m-1/2 Substituting Eqs. 3.3-12 througn 3.3-14 into Eqs. 3.3 5 through 3.3-7, then collecting the still unknown quantities at time step j on the left-hand side, we obtain gm ) = F f + ( 3)G 0 d [F i gm-1/2

• b O

gm Atj (3.3-15) gm.* 67 g m m E Ifm + 24tm)*

• I

) m-1/2 + f 24t 5 AZ AZ m j (3.3-16) At m m 2At j-1 j-1 ( ) j j j tm-1/2 tm-1/2

• b.m"j im im m

H,= g (3,3,37) 2At Fj tm + (67 )O m m Note that the quantities on the RHS of the above equations are known at time step j if we march from core inlet where boundary conditions are specified and 3-18 4

m 2st d+K)=Kd[1+(g 3)v,J (3.3-29) d 3 (0Zm - 1 + Oy 9 m Thus dividing Eq. 3.3-27 by Eq. 3.3-29, we obtain the algorithm for af: j j l E a*3 .(3.3-30) Vm 'm 03 -1+05 + Km d d Since K, depends' on a, as shown in Eqs. 3.2-40, a simple successive substitu-d d tion.is used to update the a, and hence K,. The interaction is terminated d when the absolute value of the f ractional change in a, becomes less than 10-7 s After the converged a3 is obtained, the vapor and liquid mass fluxes at the m current time step j are calculated f rom Eqs. 3.3-18 and 3.3-19 as follows: e j j 'j (Ovm * "m) Dgm (3,3,33) G 2At (0 3) m j j j (0 * ~ l + *m) P

• 2 t

G (3.3 32) -2At3) ( AZ, d is obtained from Eq. 3.2-21 for the pressurized water using the Here p tm reference pressure at the current time step. The vapor pd is assumed to be gm -the saturation value evaluated at the reference pressure. The vapor and liquid densities (Fd,and Fjg,) at the current time step j are then g calculated from Eqs. 3.3-20 and 3.3-21, respectively. The vapor enthalpy Hj is also assumed to be the saturation value at the gm current time step. The liquid enthalpy H3, is calculated f rom Eq. 3.3-17, g 3-21

e if we solve the heat concuttion equation first to obtain the surface heat L fluxes, or contain the unknown value of'Td This equation is then solved g,. iteratively using the thermodynamic relationship between H and T. g g Define the following quantities: p-I.: 26t 3 [f + A2 Ogm ) Qj gm m Vm (3.3-18) E

• bm i,

i 24t 3 [ Ff, + IA2,) Gf,) olms (3.3-19) Pkm Note that the O's are actually calculated f rom the RHS of Eq. 3.3-15 and 3.3-16. - and are theref ore known before both the F's and G's at step j are determined. Recall the definitions: d F (3'3*20) gm = o, pgm Ff,=(1-o,)ofm (3.3-21) d d d G v, (3.3-22) a o o gm gm gm g Gfm*II'"m)'m (3.3-23) Using the above definitions we can reduce the O's to p: I 3-19 l

e i r ( 2tt i 5 d of, o[ [ 1 + v,) (3.3-24) ,7 g m 2tt L* Qj,=(1-of)[1+ v{, J' (3.3-25) ,7* 3 $1nce the Q's are known up to this point (calculated from the RHS of Eqs. 3.3-15 i and 3.3-16), we can determine the void fraction o, if the velocities vd,and d g vd,are known. We employ the Bankoff-Jones correlation for slip to obtain the g I M as follows: m j j 1 - 0* V Sfn 95 (3.3-26) = v{m Kk-ad where Kf is given-by Eqs. 3.2-40. Note that d d Of, Kf of rf [ 1 + - v,) (3.3 27) g m (Q 1+Of'm+K)*I 5 ) E "5 V + (l** ) V 3* e Zm gm im j (1-o ) )

• K (3'3'23) 2tt

=(- k)V E "m

• m gm gj m

m Now, using the Bankoff-Jones correlation (Eq. 3.3-26) in the above equation, we obtain 3-20

l. i d Note that the denominator is simply Q 2m im which is known even before the 0 is determined. m The bulk coolant temperature T0 tm is calculated from Eq. 3.2-16 using the just i calculated Hd, for the liquid enthalpy. Since it is a nonlinear equation in g d,, the Newton Raphson iteration is used to obtain the final Converged Td T g g, a s follows: J if, (n+1) Tf,(n) + H** - H (n) g (3 3-33) ] where n denotes the iteration index, H (n) is evaluated from Ea 3.2-17 and g Cpg (n) is evaluated from Eq. 3.2-18. The iteration is terminated when 6 j j T,(n+1) - T I im ") < 10-7 g Th(n+1) - Tg +1 Here T is the core inlet coolant temperature. in It is instructive to note that the use of slip for the void calculation is also basic to the hydraulic model of ARROTTA. Since the concept of slip loses its physical meaning in the case of flow reversal or flow stagnation, they are not permitted in the hydraulic calculation. 3.3.2 Steady-State Algorithms - At steady _ state the time derivative terms in the conservation Eqs. 3.31 and 3.3-3 vanish. Thus we have BG g {3.3-35) 3-22

E t i l. g. h(GH)=O (3.3-35) gg g Let G denote the total flow rate per unit area (mass flux). G is known f rom g g the core inlet boundary condition. To conserve total mass, we must have Gg+G

• O (3.3-37) t t

These are the conservative equations to oe solved for steacy state. Referring to Fig. 31 we apply a finite-difference approximation in space to Eqs. 3.3-35 and 3.3-36 to obtain 02 (3'3-38) Ogm+1/2

• Ogm-1/2 + gm m

and tm A2, (3.3-39) Him+1/2

• Otm-1/2 im-1/2 + O O

H tm+1/2 Applying the definition Eq. 3.3 9 to Eq. 3.3-38 to eliminate Ggm+1/2' "' ODt'i" A2 (3.3-40) I 2) G,=Ggm-1/2

• gm g

Thus, we can calculate the average vapor mass flux for each TH region marching f rom the core inlet where the initt masG flux is specified by the boundary con-ditions. The liquid mass flux G, is then calculated from the total mass balance Eq. g 3.3-37: Ga=Gg-G, (3.3-41) g Applying the definition Eq. 3.311 to Eq. 3.3-39 to eliminate the quantities at the upper interf ace cf the TH region m, we obtain 3-23

i {' 7 .) j LZ 1 Ntm

• g iG Him-1/2.+ Otm(2)3 (3.3 42) tm-1/2 a

i Thus,. once again marching from the core inlet, we can calculate the liquid. enthalpy H for each node by Eq. 3.3-42. im Once the liquid enthalpy is obtained, the liquid temperature T, is computed g from H,'by means of the Newton-Raphson iteration scheme (Eq. 3.3-33); namely, g H - H (n) g g T,(n+1) = T,(n) + g g Cpg (n)

• ~

where n is the iteration index, H (n) and Cpg (n) are the Jordan correlations for g -liquid enthalpy and specific heat in terms of the liquid temperature T, and g reference pressure as given by Eqs. 3.3-18 and 3.2-18, respectively. The itera-tion is' terminated when Tim (n+1) - Tim ("I j < 10-7 (3.3-44) T,(n+1) - Tin + 1 g After the-mass fluxes (G, and G,) are calculated, the void f raction-a, is g g calculated as follows. First, define the ratio i G O S gm tm, "m m p -(3,3 45)' m*G,p I ~ "m g gm .where p is calculated from Eq. 3.2-21 for the pressurized water in the two-im . phase flow regimes and set to the saturation value in the vapor flow regime. 1 gm. is assumed to be the saturation density (psat). Note that we have o applied the definitions for the mass fluxes as given by Eos. 3.3-22 and 3.3-23. S -'is.the slip ratio defined as v,/v,, m g g Next, define'the quantity l 3-24

g[{Q <\\ y s ., q .i f ((__9 y r y. ;y - ] te a m lI. 4 a g i R* a*S* + -Y2m'.E (3.3-46)< s' 1 + R, 1.- o,;+ g S, ] ' we have.used Eq.'3.3-45 for R,.. Now, applying the Bankoff-Jones correla-j 'm -tion for the: slip' ratio,- 1 - o,_ S := '(3.3-47) m g "m ; m In Eq.'3.3-46',.we obtain the algorithm for a,' m"Y E ( *

• b}

a - zm m j The parameter K, as given by _Eq.- 3.2-40 depends on the void _ f raction a,.

Thus,

[ ~ an iteration-between o and _K,.is required to obtain the final a,. A simple m successive ~ substitution is used to carry 'out the-iteration. The iteration is ~ -terminated'when the absolute value-of the fractional change in a,7 ecomes less { b .than'10~7 pa -The liquid and vapor densities are then obtained from the converged o as g L follows: lu l I (3.3-49) .__ Fim "!Il~~ % - Pim ilr -F P (3.3-50) gmi" m gm Only F,'is calculated-in the code. F,is not explicitly calculated, nor. g 9 _printsd. l The entire thermal-hydraulic calculation consists of two parts: the hydraulic . calculation as described above and the thermal calculation (solution of the heat conduction equation). The two calculations are intimately related.. This will 'be. discussed after we describe the heat conduction model and its solution algorithms in.Section 4. 3-25 i

1 i i 3.4 BOUNDARY CONDITIONS The bouncary conditions in ARROTTA are inlet flow rate per unit area (mass flux), system pressure, and inlet coolant temperature. All of them are spe-cified as a function of time. They are input via tables. The relative flow into each channel is also specified. This distribution remains fixed while the total flow varies during the transient. 3.5 SPECIAL FEATURES Basically, ARROTTA contains a nonequilibrium thermal model. However, through use of the input array (card series 20 and 21) described in the user's manual.. various dif ferent thermal models can be called into existence. We list them here and suggest the user refer to the user's manual for more detail. 1. HEM model (no subcooled boiling, no slip) 2. HEM model with slip 3. Subcooled boiling but no slip 4.- Constant slip 5. Moderator properties (temperature, void fraction, etc.) do not change during a transient j 6. Any combination of the above. The value of these options lies in comparison calculations which can be used to quantify the margin lost in performing a conservative licensing calculation. 3-26

~ k 4 Section 4 HEAT CONDUCTION MODEL

4.1 INTRODUCTION

The heat conduction model in ARROTTA uses an analytical approach to obtain the average fuel (pellet) temperature, average gap temperature, average clad tem-perature, clad inner surf ace temperature, and clad outer surf ace temperature (H). Because of the analytical treatment, solution of the heat conduction equation is f ree ~ of the usual mesh ~ spacing effect inherent in the finite-difference approximation commonly used. ' The heat conduction modl involves the following assumptions: (i) The volumetric heat source is either piecewise constant in space, or any integrable function of. space (in the pellet only). (ii) The thermal properties of the fuel rod are constant in cach subregion; namely, pellet, gap and clad. They are a function of the average temperature in each subregion and are re-evaluated whenever a new average temperature is obtsined. (iii) The heat conduction in the axial direction is negligible. T_he time-dependent heat conduction equation for an infinitely long cylindrical - fuel roc is given by pC =1 (rk ) + q"' (4.1-11 p where T is tha temperature, p is the density, C is the specific neat at p constant pressure, k is the thermal conductivity 6nd q"' is the heat generation rate per unit volume. Convective heat transfer occurs at the clad outer surface between the clad wall and coolant. It ente rs into the heat conduction model as ? boundary condition. 4-1

., y W& 4 .s Q_. l

q p1,

'4.2 STEADf-STATE HEAT CONDUCTION a;;, y z At steady state the time d6rivative term vanishes, thus Eq. A.1-1 reduces to k i 15-- (rk S) ' = q" ' _ 4.k-ll ( r dr ~ dr By assuming that: k is constant in each subregion of interest and that q"' is either a constant or some integrable function of r, we can solve Eq. 4.2-1 ana-l lyticallyi to obtain the _ average fuel (pellet) temperature, average gap tem-c perature average clad temperature, clad inner surface temperature, and the clad. ' wall tempereture. 4.2.1' Pellet Region and Gap q T Consider a typical fuel element as sketched.in Fig. 4-1. The heat generation in' l the pelist is w:11-represented by a parabolic distribution at beginning of life j or at-moderate exposures. At high exposures a higher order polynomial represen-tation may be necessary. In this section we consider a parabolic distribution 'and a Bessel function distribution; however. the derivation can be applied to any integrable function. s .i ..For a. parabolic distribution the heat generation _in the pellet is writtea + q"'(r)=(b)[1+ 1)nl (4.2.2) P V R2 - 2 f Lwhere Q is the total heat due to nucle?r fission in-the pellet, V = =R 2 c AZ = wR2 f (62 =L 1"for convenience), and n is the parabolic di:tribution parameter to be .specified; usually by doing an auxiliary neutronic calculation. The permissible t crange is 0 < n < 2. When.n = 0 the heat generation rate is uniform across the pellet. ( 4-2 y.'

x; .x c.: e i t-y . c a 4* \\ cOObb cL40 G49 p - pELLN bG i Bt t c N

r r

2,! R l I l/ /"s w l Tfi, I l .w. e-l l I o F-l l 4 If2 -l -w CL. = l 2-Tei l =wr i I Tc2 i l I I

r O

R,. R R 2 3 RADIUS Figure 4-1 Typical Fuel Element Temperature Distribution t3

[N y e-LE 7 (N .gz g fi.

r2 7-dsing:.the subscript f. to denote the pellet, we have the steady-state heat con-

~ duction equation-for the pellet: N (rk ; dr ). = ( R2) [1 + ( 2

1) n)'

(4.2-3) I C r dr n R 2 L t f Itu' ltiplying the above equation by rdr, - then integrating from. zero. to r, the. point of-interest, we obtain dTg ',. Oe , 3) y 3 k I dr 2rR2 2~ M f f y Integrating again from r = Q to r, we obtain 2 f 3 - (4s k ) (I--) [(1

3) +

C T-(r) =T 4R{ r ] (4.2-5) f R 2 f f At-(tha pellet surface (r = R ), we have-f .Tf2 = Tf1 - (4 k '~ ~ f-and-

t

.,dTf Oc f (dr )r = R .k (4.2 7) = f 2nRf c. Before we derive.:the formula for the average fuel' (pellet) temperature, we will U need information for the gap between' the pellet and clad. This is obtained by, solving the conduction equation in the-gap: ,E I L l l 4-4 ,y l '.

g 1. 1 s .p -c. > y

g. > -

.,l

t'. 1 9

r.dr-(rk9 cr );=.0 (4.2-8) since the heat generation is essentially Zero.. Integrating Eq. 4.2-8 fron'R - to f r, we obtain' f dT ' dT -1 9.- A = 9 (dr )r=R (4.2-9) 'rk Rk 'I f j dr. The-RHS of the above equation can be determined by imposing the continuity of -heat flux at interface r = Rf and Eq. 4.2-7: .R kfg( 9)p, (4.2-10) Rkff( )p,g = =. .Thus, Eq 4.2-9 leads.to' a b a --(25k ) r Oc 1 -(4.2-11) dr g_ -Integrating Eq. 4.2-11 f rom r = R f to r, we obtain-s T.(r),=.Tf2 - ( q) in:( ) (4.2-12)- g At the clad inner surf ace' (r = Rf + AG), we-impose T (Rf + AG) =-T ~ c1 'h'"C' I g Tf2 = Tc3 +12k} "( + }' (** ~ } g and b) Oc (4.2-14) =- dr 2nk (R +AG) r=R g G_ g f Since we need the average temperature in the pellet and gap to evalaute the ' thermal properties, we now derive these quantities. 4-5

n i t e, 4

p. ;

4'- -1 ? i L: Define the average fuel (pellet) temperatures as a volume average:. c. t k g r T (r)'rdr f g_ H 'Tf= (4.2-15) Rf l rdr. 0 We' now ' substitute Eq. 4.2-5 into Eq. 4.2-15 to obtain f f - Oc T T 1 angf (1-n/3) (4.2-16) which, upon using Eq. 4.2-6 for Tf1 and Eq. 4.2-13 for Tf2, leads to (T -_Tc 1 ) ' * 'oiI (1 ) + l-A" (l & )) (4'2-17) f-4k k R 9 f ' If we-define the overall (composite) heat transfer coefficient V as follows, 1=-I (1 U) + 1-in (1 + $ ) (4.2-18)- U 4k 6 k R f g f ~ Eq;-4.2-17 leads to--

Qc = U (T p

f-Tc1): (4.2-19) - Since AG/Rf u l',1 we may approximate in (1+$) = AG AG . g.._ R 2R 2(R +AG) (4.2-20)- f f f F l, 4-6 l 4 s

p y l n' 7.4 ,F Hence, Eq. 4.2-18 reduces to g 'AG AG-1, 1 (3. n) + (4.2-21)' + U 4k 6 2k R 2k (R +AG) f gf g f a

which.is used in ARROTTA.

Another common representation of the heat generation in the pellet is.. based on [ the diffusion approximation for ne'utrons,which yields o i Y (r) = ( C L{e 1 (e " )/ 13(e)) (4.2-22) SRl) 2 R . t f i

.where~10 and 13.are the modified Bessel functions of zeroth and first. order and e is :a function'.of the neutronic properties 'of the fuel pellet (30) (speci-fically 'the pellet radius divided ' by the thermal neutron. di f f usion length).

Following the same procedure as before, we obtain, in this case, j I (0) 2 AG AG

1. 1 O

(4.2-23) - t .k e -[ g, ] U -0 2k R 2k (R +aG) f gf g f F Assume'.that -the -reactor core. is composed of a - collection of. repeating. unit I thermal-hydraulic (TH) cells. ~ The unit TH cell consists of' a typical -fuel rod and'.its associated coolant (moderator) as depicted in Fig. 4-2. The heat source

Q 'in the pellet is represented by-e.

f(wRf) (4.2-24) 'Q

• (1-Y 'T ) C Q c

1 2 f is the direct.(prompt) ' where Y is..the direct heating fraction for clad, Y2 i L moderator heating fraction for the coolant, C is-the ratio of the volume of'the f unit' TH cell to.that of the' fuel pellet and Q is the total nuclear heat per f (unit' volume: of the unit TH cell, obtained from the neutronic calculation (see Eq'. 2.2-20). ^ i s t' l 9 ..--. -a r

n s 6 - Using the definition of O as given by Eq.-4.2-24, we obtain the formula for the-I 4 c average fuel temperature' from Eq. 4.2-19: 2 Tr=Tc3 + (1-Yg-Y ) 22, ( ) (4*2'25} 2 Thus, if the clad inner surface temperature T is known, the average fuel tem-gg .perature can be readily determined. t mm tac AG f a, i ( /

1. 3 PITCM (p) '

q r Figure 4 2 Thermal-hydraulic Cell 4-8 f

..a; hl',i.; # lH.. f (, - e,, dp 6 -{- 2 ll k - We 'now cerive the. formula for the-average: gap temperature. ~ Define the average'

-c gap temperature as-follo~ws

-

t ( (,

Jff+lG(r)'rde f

g MG T-(4.2-26) ,j '9 'f AG-Rf+ R rdr e(, f-Using Eqs. 4.2-12, 4.2-13 and 4.2-19 in Eq. 4.2-26, we obtain ef ter some manipu-if

1ations, pw L

69 T=Tel + ji, 2 in -(1 + R )[ (T U-f f-Td) ( '2'2y) g 2k (1 4 a, 32,) r g R f q If-we define 2in(1+M) R U I ~ S.E [1 .j -(4.2-28) g: (1 +

E 1

2k f t'he average gap temperature can be expressed as 1 - T = e T, + (1-s) Tc1 B 2-20 g . No%. -1 -(1:+ $ j ~- 1 = 2 ( $ ) [1 + 1/2 (66 'jj (4.2-30) 2 R R R f f f 1 -h 2 i N

f~ e 1 ' a n~d - .in(1+ 06)=1(R0 );[1 - 1/2 ( R j)- (4.2-31)- x R. f f p 4 - since AG/R u.1.- The parameter 8 then ' reduces to f 3 U (R60 ) ._(4.2-32) l 8= 2k 1 g f. -Note that'B <'l~. In summary, the average-fuel temperature..T is calculated by Eq. 4.2 f rom f the" clad 'inneri surf ace temperature Te3 which is calculated prior to T. The f average. gap temperature T is then. calculated by Eq. 4.2-29 using Eq. 4.2-32 for' g ~B.- Since k is correlated to the average gap temperature and U depends on k, f g g J cas:given by Eq. 4.2-32 is a function of T. Hence, an iteration'between-8 and-g I is! required to obtain the final T. A simple' successive substitution is used g g -to. perform' the iteration until the absolute value of the fractional change in LT 'becomes 'ess than 10-5, l l g It 'should' be ; mentioned that k is correlated.'to the average fuel temperature f iT'f and that AG takes into account ithe thermal: expansion between the pellet and s

clad !(i.e.,-- AG depends on the. average clad-temperature T and T ).

In the ~ c f Tf; calculation',-the old values of T 'and T are used.. The correlations for k, c f g AG and kg will be presented in Section 4.5. 4.2.2 Clad Region . y } We now2 consider the clad region where we are interested in the clad inner sur-i O f ace ' temperature T el, the clad wall temperature Tc2, and the average clad tem-jo Perature T as'shown in Fig. 4-1. The steady-state heat conduction equation-in e .the clad region with a constant heat source is given by s 4-10 i ,--...i--.---i"'=="==' ' ' '

m .ml. n C, kh p di .). .g c ( ' ~33} dr.-(k;;drJ l C f C r m Where 'Y)- is th'e direct, clad heating f raction and C is the ' ratio.of the unit -. thermal-hydraulic cell volume to the total-clad volume; that is 2 2 P P (4.2-34) C =. - e 2 3(R -R ) -45AC(R +4G+6C) f 3 2 .l tor,weobtain '. Integrating Eq.'4.2-33-from r = R2 dT dT y)C Qc f (r -R) (4.2-35) e 2 ' 2 2g ( cJr=R R~k - rk 2-c* 2 WeJnow imposelthe continuity of the heat flux at the clad inner surf ace}(r=R ) 'I 2 to determine the~ first term in Eq. 4.2-35: ir dT dT k : (or Ir=R 2 g(or -R c r=R = 1 U (Tf-Tc1) (4.2-36) k 2 2 2 1 .Here Eqs. 4.2-14 and-4.2-19'have been used. Equation 4.2-35 then reduces to l l Y)C0Rf]-() J dT TC0-cf c Icf 9 r+[k (Tf .T 3) - r) (. 4.2.37) 'l = c dr 2K 2k. e c e. i 2 to r, we obtain j Integra' ting again from r =.R l lCfd TC0 0 I 3cf 2 T (r) = (Tc3 + 4k ~4k e_ c c 4-11 A

r: V:% ? i; $f. L ~r (I - Tcl)

• f 2k R

2 At the clad well (r=R ), we have T (R ) = T 3 e 3 c2. Hence Eq. 4.2-38 leads to t Tcl = Tc2*( ) k- -2d"( ( ~ 1~ 2)2 3 2 c 2 2-(4.2-39) Here.we~ have used Eq. 4.2-25 for T. Thus, the clad inner surface temperature f T can be determined if T is known. e1 c2 ~ I The -clad wall temperature T is calculated-prior to' Teg as follows. N rs t,- c2 determine the clad wall ~h' eat flux q'" from Eq. 4.2-37 by the definition. di e 9s"'Kc (p )r=R {4'2'40) 3 1 which, upon.using Eq. 4.2-25 for T and Eq. 4.2-34 for C, lead to r c I N[=2nR Op (1-Yg) (4.2 41) The surface heat flux *q'[ is related to the bulk coolant temperature T via the g heat transfer coefficient h : g kh = hw(Tc2 - T ) (4.2-42) g 4-12

r: e y q[ T' X ~ 3: Here h, depends onl the' flow regimes-defined in Sect. 3.2.. For single-phase liquid - flow, hy = hDB, the Dittus-Boelter heat transfer ' coefficient.- _Thus, i using Eq'.' 4.2-41 for q", we obtain ? O-J1"2)h Tc2 = Tg+.g f 0B ' where T is ' determined from the _ hydraulic calculation.in the coolant channel as. ) g - d' escribed in Sect.~ 3.2. r For-Tc2. ) T (which includes part of the single phase and al of the two phase 3 ~ flow regimes), tne total surf ace heat flux q" consists of the' boiling. heat flux' and-tne non-boiling heat flux; that'is, c2 - T ) = (hg +- hDB) (Tc2 ~ 5) (*' ~ }

• hg (Tc2 - T ) + hDB(T g

.s 3 +hDB(T -T) s g where--h - is the Thom correlation for the' boiling heat transfer', and T is the-g 3 saturation-temperature. If we define i E;( - hDB (T - T ). (4'.2 45) ~ 3 g we have from Eq. 4.2-44 4 i=(n+ nog)-(Tc2 T ) (4.2.46) g s whence-B Tc2 = Ts+ (4.2-47) (hb+h0B) H 4-13

^ ~ p Ik i ? ~ w i r < Wq What= is,actually done'in the code to compute T is the following. First,' d" is-gg calculated by Eq. 4.2-45.- Recall that'tne Thom correlation hB is proportional to the wall superheat, ATg =; TC2 - T, as.'shown in Eq. 3.2-23; namely,- 3 hB.= FB (Tc2 = T )2-(4*2'#0) s-where b.= 0.5358 e /630 (4.2-49) P F ~ Thus, substituting Eq.14.2-4B into Eq. '4-2.46 for h gives rise to a quadratic g equation..in.the wall superheat.aT : ~ g i I fg (Tc2 = T )2<+ hDB (Tc2.- T ) - qB = 0 _(4.2.50) ~ 3 3 The positive. root-of-this quadratic equation determines the! wall ~ superheat;- . therefore DDB + 4 F -qB hDB -: T =T 3 b O lWe now -derive t' e formula for the average clad temperature that -is needed to ~ h evaluate the clad -thermal properties. Define the average clad temperature as' follows: 'h i ~ c 4 14 w

by ;; ,Q b. i 1 R '/ 3Te (r) rer m R Tg E_2 (4.2-52) R l 3 lg rdr m 2-y Substituting Ea.. 4.2-38 into Eq. - 4.2-52 and carrying out integrations, we-obtain after some manipulation. k R2 R2 i 1;1 t-2 1 3 1

)+1,3t

). (4.2-53) + c R,g 2 in ' (- 3 ). 2 2 2in(3) 2 2 R -R -3 2 R 3 2 R 2_ 2 2 2 Yg c f [- (R +R )'* (R-R)) C0 3 2 in(b)_ 8k c .R l 2 . Note that-R.= R2 + 2aC (see Figure 4-1). Since 3 An( ).=(2aC)(1 ) (4.2-54) t 2 2 k '. and-R ). = 2 R (2aCJ (1 + ) (4.2-55) (R 2 es m, 4-15 r 4 )

1 1 we have-2 2 (l+ )~ ((1 )~4( )) 2 2 (R,-Rj 2 (-'(R+Rj+ } (4.2-55) 3 2 2 R in(2) ( 1.,4.t,_) C l 2 R 2 which tends to-vanish since ( AC/R ) = 0. Thus, the average clad temperature f 2 reduces to -4 I=a Tel.+ a2 T (4.2-57) c 3 c2 where R22 1 a s-4- 1 2 2 (R -R ) 2in( j -l l 4 R2 3 1 a2 ' (4.2-59) ~ 2 2 R (R -R) 2 in (-.3_; 3 2 R l 2-and a1 + a2

• l' Since AC/R2 << l. it can be shown that a3 reduces to 2R

+ AC 2 a3 = 4 (R2 + AC) (4.2-60) which-is simply the ratio of the volume of the inner clad region to that of the total clad region. (Recall that the clad region is split into 2 subregions, each having a thickness of aC as shown in Fig. 4-1.) Equation 4.2-60 is 4-16

p,. W, __ G

15 pl

.,n: !} r actually usedjin the code.: The quantity.a is then calculated cy the relation 2 a2 = 1 - a.. y z = Using. Eq., 412-34 for C and Eq. 4.2-59 for. a, we can reduce Eq. 4.2-39 :for c 2 -T to tne following simple expression: c3 2 Q 2) .2-61L 1 Tc3 r Tc2'(2,) )( 92"1 8 2 c J. - Here' we have u'3ed the_ approximation in( ) = in (1 + 2aC;, 260 (4.2-62) 2 2 2 i since 2aC/R2"I'- ..In_ summary, the steady-state solution of the heat conduction equation is reduced to -that:'of o the average fuel temperature. ' average gap temperature.- They are. calculated: by_ Eq. 4.2-25 for the average fuel temperature, Eq. 4.2-29 for' the l average gap temperature, Eq.-4.2-57 for the. average. clad temperature, Eq. 4.2-611 ' forEthe cladfinner surf ace temperature, and Eq. 4.2-43 or 4.2-51 for the. clad ' wall temperature. l+ '4.3 ~ TRANSIENT. HEAT CONDUCTION The transient heat conduction equation is given-by Eq. 4.'l-1 at any point of i g i< -interest in the. fuel element. To be consistent with the steady-state solution, we wish: to derive the. transient heat conduction equation for the average fuel j- . temperature,- average clad temperature, clad' inner surf ace temperature, and the - clad - walls temperature. Since the heat capacity in the gap is nearly zero, the-time-dependent' heat conduction equation of the average gap temperature need not be t solved. Instead, it-is derived from the time-dependent average fuel tem-perature and clad inner surface temperature using Eq. 4.2-29. 4-17 p^

1 i J4.3.1 Average Fuel Temperature To derive the. time-dependent heat conduction equation for. the average fuel tem-perature, we integrate Eq. 4.1-1 over the entire fuel pellet, then divide each term by the volume of the pellet to obtain ( R BT R BT R 1 jf f g f 3 y 3 f,,, o Pi St. V o ar f Br V o ~ V f (4.3-1) f f where R-R{ f Vf=f rdr = (4.3-2) .lf we assume that the' heat capacity (o Cpf) is constant in space and use the f definitions of. the average heat. source and the average fuel temperature, Eq. L4_.2-15, the above equation becomes t f C =( )R kf( I)r=R <4" ( 4 ~' 3-3 ) P pf f f f The average heat source term is given by i <4 > = (1-Y1 - Y ) C'f Q (4.3-4)' 2 f ^ The term invol'ving the slope at r = R ' n Eq. 4.3-3 is treated as follows. Note i f 'that the time derivative term of the heat conduction equation is nearly zero in ~ ' the gap (where the density and specific heat are vanishingly small), and that -there is no heat generation in the gap. Thus, we have for the gap, L 4-18

~ ^ er 'Is 0: i di 1 d -(rk9 dr ) = (4.3-5)' 9 r.. dr to r = R, we obtain-Integrating from r = Rf 2 dT 'dT fg( jr=R rRk2g( Jr=R (4.3-6). Rk f 7 1 4 I Forlthe continuity of heat flux at'r = R, we must have f dT dT dT 2g( 9)r=R (4*3-7) -Rkf f _( jr=R[=Rkfg( 9)pg =Rk 2 l u c. Since - AG, the' gap width, is very small, it is a good approximation to evaluate the slope at-r = R ;as -follows: 2 T-T (dT (*~}

• ~

dr =R AG/2 2 Hence.. .'dT - 2R Rkff( I)r=R k f-Tci) (4'3~9)'

• ~

p as g g

We wish to express T in terms of the average fuel temperature T. This is done

~ g f by. defining T for the transient as follows: g T(t)=6(t)-T(t)+[1-6(t))Tc3(t) (4.3-10) g f .The above definition is formally _ exact if-we calculate the weighting factor 6(t) as a1 function of. time. Thus, Eq. 4.3-9 becomes 4-19 m.

1 1 s ;q 1 v s i 3 dT -2R .Rkf 9 '( Jr=R k B(t)'(T - Td) I#'3*III { f_ f g The. transient-heat conduction equation Eq. 4.3-3, for T then reduces to f dT 4k R I I C 9 .of pf = (1-Y -Y )LC 0ff-( 2 ) 8 (I - Td) (4.3-12) 3 2 f f - i .To -calculate the time-dependent 8 rigorously is difficult. Instead..we assume that the steady-state expression (Eq. 4.2-32) for 8 is applicable in the.tran-sient, i.e., B(t)L= ( ) ( 4 ~.3-13 ) ' ] 9-f l Note that 8 is a function 'of time since both U and k are updated at'each time-g step. Using this approximate.B and-further approximating R /R = 1 (since R2* l 2 f Rf _ + AG/Rf u-1), we obtain the final time-dependent heat conduction equation for T that is used in ARROTTA-f 1 dT. of pf f = (1-Y -Y ) C Qf f * [ U (Tf -Tc3 J .(4.3-14) C 2 ~~ 3 2 f We' are interested in thegnumerical solution of the above' equation at the con J tinguous : discrete time point: t, t, .... t).3, t,..., where tj = + At). 0 3 j applying a' semi-implicit time differencing-to Eq. 4.3-14, we obtain - i 4-20 ? ___;__.__.___m__._._____._m- ^'

1 4 i 4 t [ r, ., y y, fM j-Tj.1 j j U*1(T ' d) (4.3-15) d 1 j.1 f-.(}. je f-gy d -f s L d SolYing for Tf yields the. numerical algorithm: j- - w{-l%*1+ wf Q{ + w{'-1 Th l T.. (4.3-16): p

,j-1 g -1 1.

j w 3- - where wd~lE(ofpC jjd-1: (4.3-17a) g d w (1~T 'Y ) C 'at (4.3-17b) 1 2 1 2 f j 9-1 E 2..U -1 d at/Rf '(4.3-17c) 3 3 Note'that the heat capacity (of pf) and the coefficient U are evaluateo with the - C ! previous ' temperatures at time step (j-1). It is' also important that if i Tf reaches the meldng or-boiling temperature the' calculation' becomes one for fuel enthalpy' rather than temperature. 4.3.2fAverageCladTemperature j.- . Before.'we derive the time-dependent heat conduction equations for Te1 and Tc2' .we ; shall first derive the transient heat conduction equation' for the average j clad temperature T (t). This is done by multiplying Eq. 4.1-1 by the element of c volbme'rdr, ' integrating the resultant equation from r =R ~ p to r = R, then 3 dividing each item by the total volume of the clad region. The result is 4-21 t

O., p Lr ie dY [ dT . dT '#c pg og' = [Rk3g(dr Ir=R.-R2 cl r [.[ C 2 2 (g. 2) 3 d r=R2 (4.3-18) s- 'l s-- where = Y C 0. e 3gf M. W - The, boundary condition at the clad wall (r = R ) specifies that 3 s dT e c -(dr )r=R =.R .Rk 3 h,-(Tc2 - T ) (4.3-20) t 3 where h,..is the heat ? transfer. coefficient that varies with the flow-regimes . defined in Sect. 3.2, and T is the bulk coolant temperature from the hydraulic g i . calculation. To evaluate the term involving the slope at r = R in.Eq. 4.3-18 we apply the-2 o . continuity'of heat fl'ux at r = R2 and use Eq. 4.3-7: i dT dT dT dT Ek2 c (dr 3r=R =Rk2g( Jr=R =Rkfg(- j=R =Rkff-( Jr=R , (4.3-21). r 2-2 f f ' Applying l the. same approximation.regarding the 8(t), Eq. 4.3-13, to the above continuity condition, we obtain y dT Rk2 c (dr 3r=R (TI~TCl J (*~ } "~ 2 This is an approximation whose implication was discussed in Sect. 4.3.1. ?. u 4 22 .e - ~ L

.x < s g w .7 i =d f ,7 Using Eqs.._4.3-20 and 4.3 in Eq. 4.3-18, we obtainf the time-dependent heat. . conduction equation. for the average clad temperature i (t). J g .t

W 6 - T )) _ R'3 hf(Tc2 - T )-

dT f c g (4.3-23) e 'Y 0 Q pc e PC dt ' 3 C f + 2AC(R +6G+aC) 24C(R +aG+AC). --= f f 'where we have:Used the relationship, j 2 2 '(R - R ) = 4AC (Rf +1aG + AC) (4.3-24) 2 It should be mentioned that Eq. 4.3-23 is not solved in the code. Rather, we wish-to derive f rom Eq.14.3-23 two separate transient. heat conduction equations LJ forf Tei(t) and 4.Td2(t), - respectively. This is described in the following sections. I L 4'.3.31 Clad ~1nner Surface Temperature We now make t a 'very basiciassumption th'at the average' clad temperature can be represented at each time by a' linear combination of-the clad inner surf ace tem-o as given by'Eq. 4.2-57, i.e., 'peratureTf3andthecladwalltemperatureTc2 a \\ Y(t)L=ai Tc3 (t) + a2 Tc2 (t) (4.3-25) e whereLa +'a = 1. As we have shown before, this is valid at steady state., It. g 2 is', however, an. approximation in the transient. It may be argued that, for most reactor transients of. interest with. ARROTTA, the form of the clad temperature [ profile is expected to stay close enough to the steady state' form so that-Eq. l L ,4'.3-25 is. a-good approximation. - Adding the term.k (Tc3 - Tc2)/(2AC)2 on botn sides of Eq. 4.3-23 and applying c l' the definition (Eq. 4.3-25) for T .we obtain c l-4-23 )} ,:l ll _ '

~ y j 1 s U(T - T )) dT k e3 f_ c PC

• 1 c pg dt

-aYCQy3cf-2aC(R +6G+6C) koC) f l dT R h,(Tc2-Ti) k ep. 3 e C +eYCQ = - a2'c pg dt p 3 g f 24C(R +aG+6C) C f 4( AC) (4.3-26) The LHS of Eq. 4.3-26 is. approximately the heat balance -for ~ the inner clad region; whereas, the RHS is' that of the. outer clad region.- Therefore, we assume that both sides of Eq. 4.3-26 are identically zero. Setting the LHS of Eq. - 4.3-26 equal to zero, we obtain the time-dependent heat conduction equation for Tel(t): i dT U (Tg - T.3) k P C '1 cC = a)Y C Q3 c f +: 26C(R +aG+AC)" e pg dt f 4(AC) Applying a semi-implicit time differencing:to the above equation, we obtain the I numerical algorithm for Teg: 9 'I j~l d d + e "I Tf+e4 c2 d d d d-I d T T 1 c1 + e Qf. 3 e Tc1- + (4.3-28)- 1 ed*1 + eJ~l + '8 -1 3 1 3 4 where-j-1 l E ag (o Cc pg)d-l (4.3-29a) e j E a Y c at (4.3-29b) 2 1ic j 4-24 :y.,

erm. u a @',.9 1 (l)A' l '} gi, s i t s j.1-U 4t) I ' '3 (4.3-29c) 2AC(R +4G+4C) f d .k' At jal: 0 s (4*3*29d)' N C 4 2 4(oC Note that', 'if one wishes. to obtain a fully implicit solution for' all the tem = [ - peratures, an iteration is.necessary to update successively the thermal proper. -ties that' depend on.the average temperatures. - This 'is not done in ARROTTA. In the code the. thermal-hyoraulic time step is the same as the'neutronic' tine step- ' hich is usually an. ceder of magnitude-smaller than would be necessary..f or-the w For suen a.small Lt-it is therefore a goo'd approximation thermal-hyd rauli cs. j t'o calcul'ation'the temperature with.the semi-implicit scheme using coefificients j ~ b ~ based on the properties at the preceding time step. o 4.3.4 Clad Wall Temperature

f To obtain the time-dependentlheat conduction. equation for the. clad : wall tem-

. perature we-set tha RHS of-Eq. A.3-16 equal to zero: } ~di k R hw (Tgp - T,) 3 c2 c (4.3-30)~ .(Tc1 - Tc2j 2AC(R +6G+4C) apC =aYCO f + 4(oC)2 2 C PC dt-21C f N' t. Again opplying -the: semi-implicit' time 'differencing to the above equation, we obtain the numerical' algorithm for the clad wall ' temperature: 4 T jl+ 6f"'3 f + 6j~I f3 + 6f Tf d E -j 6 0 T e (4.3-31) T.c2 = 6 -1 + 6 -l + 6d 0 d 1 3 4 r 4-25 b

1 i where J j*1 2.('c pe)j.} 6 54 C (4.3-32a) j 1 j 6 IaTC At (4.3 32b) 2 2Ic j j.1 kf'I At 3 6 a (4,3 32c) 3 4(40), t d R At j h 3 3 6 a (4,3 32d) 4 2AC(R +4G+4C) g Note that by substituting Eo. 4.316 and Eq. 4.3-2B into Eq. 4.3-31 it is .possible to represent the outer clad temperature in terms 'of the unknown liquid coolant, temperature (a.ll other cuantities being known from the previous-time. step). this may be referred to as the resolved outer clad temperature equation. In. solving the system of equations, the outer clad temperature appears as a driving function in the surface heat flux (Eq. 4.41 et sec). By substituting the resolved outer clad temperature equation into the heat flux relations, we produce.new relations which formally contain only the liquid temperature as the unknown. Thus, the hydraulic chiculation is formally uncoupled from-the fuel rod calculation ynd af ter the completion of the hydraulic calcul'ation T 15 c2 back calculated from Tg and then in sequence Tc3 and Tg are found. After Teg and Tc2 are c61culated at each time step, the average clad temperature at time step j is calculated by Eq. 4.3-25; namely, l 4-26 b- ,~, r .w

[ ~ i 1 4 Tf = a. T ) + (1-a)) T (4'3'33} g 2 where a3 is given by Eq. 4.2-60. The accuracy in the clad temperature calculations depends on the definition for a3 as well as that for B(t) as in Eq. 4.3-10. The validity of the approxima. tions for ag and B(t) depends on how well behaved the temperature profiles are in different parts of the fuel rod. In both the steady state and trcnsient the average gap temperature is used to I obtain the overall heat transfer coefficient. The gap temperature at time step j is legged, hence what is used is (1-8 *l) T (4.3-34) 3 B *1 T l T g f 4.4 CLAD SURFACE HEAT FLUX I The clad heat flux is calculated from the outer Clad surf ace temperature (Tc2) l and the bulk Coolant temperature (T ) using various heat transfer coefficients g which depend on the flow regimes defined in section 3.2. 4.4.1: Steacy-State Clad Surface Heat Flux -In steady state,.all heat lee.as the fuel and enters the coolant. Except for the subcooled boiling. reg'me, the surf ace heat flux is calculatable directly from the fission heat and the fuel and clad temperatures need not be calculated until after the hydraulic calculation. In the subcooled boiling regime the clad outer temperature (Tc2) determines the split in the total surf ace heat flux be-tween that portion which heats the liquid and that portion which produces vapor. In this regime we must define: Convective heat flux: 4-27

W L I i n 1 L 4 i DB (T*g - T ) (4.4-1) { q a h g g i Boiling heat flux: q[g a hg (Tc2 - T ) + hDB (Tc2 - T*g) (4.4-2) g 5 where T*g is the film temperature ((Tc2 + T )/2). g The 6humption here is that that portion of the convective heat flux (hDB ' (Tc2-T )) entering the film above the film temperature produces voidage g while the rest produces convective heating of the liquid. '4.4.2 Transient Clad $Urface Heat Flux The calculation of the transient clad surf ace heat flux is more involved. it is calculated from the clad well temperature (Tc2) and bulk coolant temperature (T ) using various heat transfer coefficients, which depend on the flow regimes g defined in Sect. 3.2. For single-phase liquid flow where T is less,than T, the saturation C2 g temperature, 4"*hDB {Ic2 - T ) (4.4-3) g where h is the Dittus-Boelter heat transfer coefficient given by Eq. 3.2-8. DB For the' single phase liquid flow where Tc2 > T, 3 4"= hog (Teg-T)+hg (Tc2 - T ) (4.4 4) s g g where hg is the Thom heat transfer coefficient given by Eq. 3.2-23. j 1 l 4 28 1 7

y, L i

k. :.,

1 t For the two-phase subcooled boiling regime where T 2>13 g g and T < T, Eq. (4.4 4) is still valid, but q, is taken as the sum of q"g ' and q g which are defined as in Eq. (4.4 2, -2). e For the two phase bulk (scturated) boiling regime where Tc2 > T, and T =T. g 3 q'g = (hDB + b ) (Tc2-T) (4.4-5) g g For the single-phase vapor flow regime where Tg>T, 6[=h IIc2-T) (4.4-6) DB g where h is the Dittus-Boelter heat transfer coefficient for superheated steam DB as given by Eq. 3.2-48. 4.4.3 Special Features In steady state the code will produce an error message and stop if the void f raction reaches unity or if any fuel reaches a fully vaporized. state. In a transient, however, several special heet transfer and flow states are allowed: 1. Fully vaporized flow is allowed, but only the ef fects on the metal temperatures are accounted for. 2. A user-supplied DNB correlation is permitted as Subroutine DNBR (Param) with a transmitted parameter list. heat flux through the clad QDNB = i coolant temperature TCO = outer clad temperature TR3 = coolant enthalpy TFV1 = vapor flow rate (G ) TFV2 = 9 void fraction ALFA = 4-29

i i t liquid flow rate (G ) TFVA = g TFV5 liquid density = PRS pressure = inlet enthplpy HIN a saturatian enthalpy HSAT = HFG heat of vaporization l = DEO hydraulic diameter = TH(4) = fuel channel length REG flow region number j = ANS answer = 0.if below DNB a = 1.0 if above DNB ? If ANS = 1.0, the channel is in DNB. If we are in flow boiling, we calculate a transition heat transfer coefficient based on a linear interpolation - in void fraction - between what would be the coefficient if we were not in DNB and the single phase vapor coefficient. Radiative effects are not accounted for. Alternatively, if we are in single phase flow, a subcooled heat transfer coef-ficient is calculated based on steam flow and a Prandtl number exponent of 1.23. 3. If the void fraction has reached unity or the clad temperature is above the Leidenfrost temperature, a single phase vapor coefficient is calculated. Here we define (in 'F): TLON = T + 500 (1 - e-0.022 8 ) g where P is the pressure in psia. 4.5 FUEL ROD MATERIAL PROPERTIES This section presents the material properties of the pellet, gap and clad that are used in ARR0TTA. The properties are expressed in British units in the following. 4-30

.s' t. .For Pellet (UO II 2 3*9 8I + 6.02366X10'3(0.46+T /1000)8 (4.5 1) kf= g (0.69261+T /1000) f where k is the pellet thermal conductivity in BTV/sec-f t 'T and T is the g f average fuel temperature in 'F. W 3 pf

• Of0 [1 - ( y 10 ) T /1000)

(4.5-2) f L where pf is the fuel (UO ) d'"SitY i" ID /ft * 'fo is the density at room tem. 2 m perature, and g is the pellet thermal expansion coefficient ('F'l). The specific heat C is a piecewise constant function which depends on the pg average fuel (pellet) temperature (T ). The fuel enthalpy in turn is a f piecewise linear function of the fuel temperature: H(T ) = Cpf (Tf-T) (4.5-3) f n where H(T ) is in BTU /lbm, Cpf.in BTU /lb, *F ang Tf and To.(a ref erence tera-f perature) are in 'F. Table 4-1 gives the specific heat and reference tem-L perature needed to calculate the fuel enthalpy in five temperature ranges up to the melting point of UP2 (5247'F). It also' lists the enthalpy at tho upper limit of each range. The heat of. fusion is taken as 126 BTU /lb, and the molten specific heat is 0.12 BTh/lb *F. The boiling point is taken as m v TFB0ll = 7800 + 860 log 10 (P/1500) and the heat of vaporization as 541 BTU /lb,, 4-31 1

p

p.,

b. h For Gap: 'k = 1.314x10'3 (T +460)0.665 (4,$,4). g g Table 4-1 i i FUEL HEAT CAPACITY AND ENTHALPY Maximum Maximum Specific Reference Temperature Enthalpy H(Tf) Heat C f, Temperature Range T, 'F BTU /lb cal /gm BTV/lb,p*F . T, 'F f m o 1 1500. 105.0 58.3 0.07332 68. 2. 3092. 228.8 127.1 0.07776 149. 5 3 4036. 322.4 179.1 0.09915

785, i

14 4591. 397.6 220.9 0.1355 1657. 5 5147. 477.0 265.0 0.1428 1807. where Tg is the average gap temeprature in 'F and'k is the gap thermal conduc-- g tivity in BTV/sec ft *F. The gap thickness AG in ft is given by AG = AGO * (*T R)Tf + (oTC R) Tc + AG, (4.5-5) f where a is the pellet thermal expansion coefficient (*F'1), s is the clad i TC thermal expansion coefficient ('F'1), AG is the nominal gap thickness (ft) and O . 6Ga is the effect of f asperities i.e., the deviation from nominal dimensions (ft). The user also has the option of using a constant gap conductance spe-cified in the input in lieu of having the code calculate kg and 4G. 4-32 -ww -.-g -,w a + -s.

i i l For Clad (2irca11oy-21: j 8.99 + 0.0346 T + 7.25X10~0-T C C (4.5-6) k = g 1 + 0.005 T e I I ' s the clad thermal conductivity in BTV/sec-ft*F and T is the average where k i e e i clad temperature in 'F. ) (o c ) = 27.8 (1 - 1.08x10-5 I ) (1 + 3.515x10-' Y ) (4.5-7) c pe c e t where T is the average clad temperature in 'F and o C is the clad heat capa-e c pg 3 city tri BTV/f t

• F.

i h k l I l t 4-33 l* b

I t f l + l i i i i Section 5 REFERENCES l 1. K. S. Smith, "An Analytic Nodal Method for Solving the 2-Group, Multi-dimensional Static and Transient Neutron Diffusion Equations," Nuc. Eng. Thesis, Dept. of Nuc. Eng., MIT, Cambridge, MA., (February,1979). 2. D. J. Diamond, H. S. Cheng, and L. D. Eisenhart, "BEAGL-01, A computer Code for Calculating Rapid Core Transients, Volume 1-Modeling," EPRI NP-3243-CCM Voi.1. Electric Power.Research Institute (1983). 3. J. A. McClure and R. L. Hatteberg, "EPRI Environmental Library Computer Code Manual," Part II, Chapter 9C, Advanced Recycle Methodology Program Computer Code Manuals, CCM-3, to be issued by Electric Power Research Institute (1985). 4 FORTRAN VERSION 5 REFERENCE MANUAL, Publication Number 60497800, Control Data Corporation. 5. IBM System /360 and System /370 FORTRAN IV Language, Publication Number l GC28-6515-lD, International Business Machines Corporation. 1 6. UPDATE REFERENCE

MANUAL, Publication Number

60497800, Control Data i.

Corporation. 7 D. M. Ver Planck, W. R.

Cobb, R.

S. Borland and P. L. Bersteegen, '" SIMULATE-E: A Nodal Core Analysis Code for Light Water Reactors," EPRI NP-2792-CCM, Electric Power Research Institute (1983). '8. F. Bennewitz, H. Finneman, M. R. Wagner, " Higher Order Corrections on Nodal Reactor Calculations," Trans. Amer. Nuc. Soc. 12,250(1975). 9. K. Koebke, "A New Approach to Homogenization and Group Condensation" 1 AEA Technical committee Meeting on Homogenization Methods in Reactor Physics, l l -Lugano, Switzerland (November 13-lb, 1v/5). (- 10. K. S.

Smith, A.

F. Henry, and R. A. Loretz, "The Determination of Homogenized Parameters for Coarse Mesh Nodal Analysis," in 1980 Advances in l Reactor Physics and Shielding, Sun Valley, Idaho (Sept. 14-19, 1978). 11. K. S. Smith, " spatial Homogenization Methods for Light Water Reactor Analysis," Ph.D. Thesis, Dept. of Nuc. Eng., MIT, Cambridge, MA., (June, t l 1980). 5-1

e 1 12. T. A. Porsching and A. F. Henry, "Some Numerical Methods for Problems in i Reactor Kinetics," in Numerical Solution of Field Problems in Continuum Physics, American Katnematical 50ctety (1y/U). 13. R. S. Varga, " Matrix Iterative Analysis," Prentice-Hall (1962). 14 G. E. Trapp, personal communication, Brookhaven National Laboratory (1983). 15. A. M. Kime, " Algorithmic Analysis Comparing Gauss-Seidel and Neumann Series Approximation in Evaluating Smith-Henry Neutron Diffusion Equation Model," M. S. Thesis, Dept. of Com. Sci., WVU, Morgantown, W. V., (Feb.1985). 16. L. A. Hageman, "The Chebyshev Polynomial Method of Iteration," WAPD-TM-592 (1966). 17. L. A. Hageman and R. B.

Kellogg,

" Estimating Optimum Acceleration Parameters f or' use in the Successive Overrelaxation and the Chebyshev Polynomial Methods of iteration," WAPD-TM-592 (1966). 18. L. A. Hageman and R. B. Kellogg, " Estimating Optimum Overrelaxation Parameters," Mathematics of Computation, Vol. 22, No. 101, 1968, pp. 60-68. 19. C. J. Pfeifer and C. J. Spitz, "PDQ-8 Reference Manual," WAPD-TM-1266, (1978). 20. L. A. Hageman and D. M. Young, " Applied Iterative Methods," Academic Press (1981). 21. A. L. Aronson, H. S. Cheng. D. J. Diamond and M. S. Lu, "MEKIN-B: The BNL Yersion of the LWR Core Dynamics Code MEKIN," BNL-NUREG-2807 Brookhaven National Laboratory (1980). 22. G. S. Lellouche, "A Model for Predicting Two-Phase Flow," NBL-18628 (RP-1032), Brookhaven National Laboratory (1974). 23. F. W. Dittus and L. M. K. Boelter, University of California, Berkeley, Pub. M.,2,433(1930). 24 W. B. Jordan, " Fits to TherpoJynamic Properties of Water," KAPL-M-6734 Knolls Atomic Power Laboratory (1967). 25. J. R. S. Thom, W. M. Walker, T. A. Fallon, and G. F. Reising, " Boiling in Subcooled Water During Flow up Heated Tubes of Annuli," Proc. Inst. Mech. Eng.,180, 226 (1965-1966). n -26. W. H. Jens and P. A. Lottes, " Analysis of Heat Transfer, Burnout, Pressure Drop, and Density Data for High Pressure Water," ANL-4627, Argonne National Laboratory (1951). 5-2

l 27. W. T. Hancox and W. B. Nicoll, "A General Technique for the Preciction of Void Distributions in Non-Steady, Two-Phase Forced Convection,' Int. J. Heat Mass Transfer, 14, 1377-1394 (1971). ~~ 28. S. G. Bankoff, "A Variable Density Single-Fluid Model for Two-Phase Bubble Flow," Int.d.MultiphaseFlow. 2,, 79 (1975). l 29. A. B. Jones and D. G. Dwight. " Hydrodynamic Stability of a Boiling Channel," KAPL-2208, Knolls Atomic Power Laboratory (1962). 30. R. T. Lahey, Jr. and F. J. Moody, The Thermal-Hydraulics of a toiling Water f Nuclear Reactor, American Nuclear Society, La Grange Park, IL (1977). t i i F 1 5-3

p-- i m ks. A n- \\P"END1x: Derivation and Evaluatiori of Spatial Coupling Equations i (To be submitted st a later date) h, Y t 9 i b P i f f f 6 4 /}}