ML18025B710

Method for Development of One-Dimensional Core Kinetics Data for RETRAN-02.

ML18025B710
Person / Time
Site:	Browns Ferry
Issue date:	09/21/1981
From:	Forkner S, Keys T, Winkler E TENNESSEE VALLEY AUTHORITY
To:
Shared Package
ML17276B803	List: ML17276B802 ML18025B710 ML20040D684
References
TVA-MDS-553, NUDOCS 8202020156
Download: ML18025B710 (55)

Text

TVA-MN-553 METHOD FOR DEVELOPMENT OF

'CORE. KINETICS 'NE-DIMENSIONAL DATA FOR RETRAN-02 S. L. Forkner E. N. &inkier T. A. Keys IIIIOTICE THE ATTACHED FILES ARE OFFICIAL RECORDS OF THE DIVISION OF DOCUMENT CONTROL. THEY HAVE BEEN CHARGED TO YOU FOR A LIMITED TIME PERIOD AND .

MUST BE RETURNED TO THE RECORDS FACILITY BRANCH 016. PLEASE DO NOT SEND DOCUMENTS CHARGED OUT THROUGH THE MAIL. REMOVAL OF ANY PAGE(S) FROM DOCUMENT FOR REPRODUCTION MUST BE REFERRED TO FILE PERSONNEL.

OBGIIQt + 50- Z. 5'7 COIttral ~ 8~O ~aO (55 DEADLINE RETURN DATE Date J=ZR=&Z Of DBCIIIIIBII REGNATORY DOGKH FRE Re%ran Task Force, ~

Nuclear Fuel Branch IIEMLA ORY HOCKEY FILE COPY RECORDS FACILITYBRANCH 8202020186 820i22 PDR ADOCK OS000259 P PDR

MDS-553, September -21, l98l METHOD FOR DEVELOPMENT OF ONE-DIMENSIONAL CORE KINETICS DATA FOR RETRAN-02 S. L. Forkner E. N. Winkler T. A. Keys G. W. Perry Nuclear Fuel Branch Tennessee Valley Authority

l. INTRODUCTION The capability to treat reactor core kinetics with a quasi-static approximation to the time dependent neutron diffusion equation .in one-dimensional (axial-slab) geometry is being added to the RETRAN.program (reference 1). In order to use the one-dimensional (1-D) kinetics capability, a library of cross section data (as a function of water density, fuel temperature and control state) which accurately represents the reactor in the axial-slab geometry is required.

For boiling water reactors (BWR's); the ma)ority of design and analysis calculations involving coupled neutronics and thermal-hydraulics's performed using three-dimensional (3-D) steady-state simulators. The .

3-D simulator often has the best representation of coupled neutronics and thermal-hydraulic phenomena that is practical for routine use. Generally the '3-D simulator has been extensively compared to actual plant operating data and normalized if necessary. Therefore, it is reasonable that the 1-D cross section data be based on 3-D simulator analyses..

The procedures presented in this report allow only one value of each cross section quantity for each 1-D model region to be obtained from a 3-D simulator case. Thus several 3-D simulator cases are required to obtain sufficient 1-D cross section data to define its dependence upon th' independent variables (water density, fuel temperature and control state in each 1-D model region). There are alternate procedures which do not require as many 3-D simulator cases; however, they have less potential for providing consistency between the 1-D model and the 3-D simulator.

There are several important elements in obtaining a suitable 1-D cross section library from the 3-D code including:

(a) determining cross sections for physics data available

~ in the 3-D code (b) collapsing out radial (x-y plane) cross section dependence by appropriate weighting (c) determining 1-D radial buckling for consistency with the 3-D code (d) transformation of independent variables between 3<<D and 1-D models (e) selection of appropriate 3-D cases for collapsing,to 1-D fitting of

'f) 1-D cross section data to RETRAN polynomial .

equa tions (g) verification of procedures The treatment of each of these elements is presented in later sections of this report.

2. DETERMINING CROSS SECTION DATA The TVA 3-D simulator CORE (reference 2) will be used as the basis for'development of the 1-D cross section data. Since the CORE code utilizes. lattice physics data based on a two-group model, the 1-D cross section data will also be based on the use of two prompt neutron energy groups. The two-group cross section data required by the 1-D model consists of the 19 quantities below'or each axial region:.

1. total delayed neutron fraction, 8

2. fast group absorption cross section, E 1 2

3. fast group transverse (radial) buckling, BIR

4. fast group diffusion coefficient, D>

5. fast group removal (slowing down) cross section, E r

6. fast group kappa times fission cross section, kZ

7. ,fast group Nu times fission cross section; vE

8. fast group average neutron volocity, Vl

9. fast group Nu (neutrons per fission), v

10. fast group microscopic. boron-l0 absorption cross section, a thermal group absorption cross section, E 2

12. thermal group transverse buckling, B2R

13. thermal group diffusion coef f icient, D2

14. thermal group microscopic Xenon absorption cross section, a x ~

15. thermal group Kappa times fission cross section, xEf2

16. thermal group Nu times fission cross section, vZf2

17. thermal group average neutron velocity, V2

18. thermal group Nu, v2

19. thermal group microscopic boron-10 absorption cross section, o 2

Thr kappa's have units of MeV, th~ ~

velocities ar< in cm/sec, and the remai.ning units'are standard.

The development of the fast and thermal group radial bucklings will be discussed in section 4. The fast and thermal group B1O absorption cross sections are not currently used. The remaining.l5 quantities are developed from'he physics data parameters available inside the CORE-code which are:

1. infinite multiplication factor, k~

2. age to thermal, v

3. fast absorption plus thermalization cross section, E 0

4. thermaliza t ion probabili. ty, p

5. fast fission factor, c 2

6. thermal diffusion area, L

7. thermal"absorption cross section, E a22

8. average energy release (MeV) per fission neutron, <e/v>

9. thermal group microscopic Xenon cross section, a x

10. to tal delayed neu tron frac tion,

11. fast group average neutron velocity, V h

12. thermal group average neutron velocity, V

13. average energy release (MeV) per fast group fission, ic>

14. average neutrons released per fast fission, v>

15. average neutrons released per thermal fission, v2 The first nine of these quantities are part of the standard lattice 1

physics data utilized in the CORE neutronics solution. The remaining six quantities are special input data added to CORE for use in obtaining the 1-D cross sections for RETRAN.

'I'I>> COIII.

~

Iatl. I>'( ~

.I>I>y>> I> ><>>>>>L I L I>>: may I>> ~;> I'<<>>> l lo>>>>l' l><'w>>la I

,values of some or all of, the following'independent variables.'a) exposure (b) currentwater density (implies moderator temperature)

(c) exposure averaged water densi'ty (d) current fraction of control rod insertion r

(e) 'exposure averaged control fraction (f) fuel temperature (g) Xenon concentration (h) soluble boron concentration The general form of the relationships betw'een the independent variables C

and the lattice physics data is discussed in section 5 of reference 2.

The details of these relationships are given in the input description for each CORE code version. Using these relationships, each, of the 15 CORE lattice physics quantities can be evaluated at each of the 3-D nodes for a given case utilizing the final converged values of the independent variables.

Comparing the lists of data required by RETRAN and that available

'in CORE indicates that the quantities. o , 8> V , V , 'v , v , and Z x 1 2 1 2 a2 appear on bo th lis ts and thus no special manipula tions are required to obtain these items for each of the 3-D nodes. The remaining RETRAN cross sections are obtained from the relationships below:

21 Z ~ pE' 0

22 E al =E 0 (1 p) 2.3 D = E T 1 o

2 2.4 D = Z L a

2.5 vZf1,,~ Z 0

k~(c-1)/c 2.6 vZf2 = E a22k'/(cp) 2.7 eEf1 E k~(K /v )(c-1)/s 2.8 MZf1 ~ v k Z o fl Ef1~(E a22 E r)

Utilizing equations 2.1 through 2.8 allows the remaining 1-D kinetics cross sections quantities to be evaluated for each of the nodes in the CORE case.

(,'l(OIIH !II(( I'll)N (OUI IAI'tiIN(

Using the procedures described in the previous section the RETRAN cross section data is obtained for each node (ijk) in the 3-D simulator case. The cross sections at all nodes (ij) in a given axial plane (k)

.are weighted to obtain the axially dependent cross sections.

3'l Z k

Z Z IP' Z W The appropriate weighting function (M) for each cross Section quantity (g) is cons true ted from the neutron dis tribution calcula ted by CORE. The cross sections can, on option, be weighted by either flux qr the product of flux and an approximate adjoint flux.

The CORE neutronics (described in section 3 of reference 2) do not directly calculate the neutron flux. Instead, CORE computes a quantity called the "nodal'leakage" (L ). CORE also utilizes an absorption ijk defined such that Aiijkk Li.k operator (Ai.k) ijk is the rate of fast neutron absorption (and thermalization) in the node.

The fast neutron flux (actually the product of. flux and nodal volume) can then be calculated as:

f

'Aijk Lijk'"oijk where Z 0

is defined in section 2. As described in reference 3, the fast adjoint flux, thermal flux and adjoint can be approximated by:

3.3 ijk

~t ~ ~f Z r /Za2) ijk ijk 1+L2B2 .ijk

(

t* ~ fA ijk ijk keff (i+L+II2) ijk

/

2 Where k ff is eff the 3-D code effective multiplication factor and B is defined .by:

(k /k ff)-1 Xjk ~v+ L flak The cross sections are collapsed using equation 3.1 with' defined as either the flux or product of flux and ad)oint flux for each energy group. Diffusion coefficients and group average velocities are inverse weighted, e.g.,

3.7 D k E Wijk/E {8j)k/D i)k)

The delayed neutron fraction is weighted by:

f .X t k fl f2 if'jk fl $2 ijk i.e., the delayed neutron fraction is source rate (importance-source on option) weighted at each level.

The procedure described above produces 1-D cross sec'tions for each axial level (typically 24) in the 3-D simulator case. As an option the cross sections can be further collapsed to obtain one value of each cross section to use for N consecutive 3-D simulator axial. levels {i.e.,

only 1/N as many 1-D neutronic regions as axial levels in the 3-D simulator case).

,4. 4 DETERHINING RADIAL BUCKLING The collapsing procedures outlined in section 3 produce a set of two-group cross sections, neutrons per fission, average neutron veloc-ities and a delayed neutron fraction as.a functi'on of axial level for each 3-D simulator case. This leaves the axially dependent radial'uckling for each group to be determined.

4.1 Fast Grou Radial Bucklin The axial variation of the fast group radial buckling is readily determined from the converged neutron distribution calculated by CORE.

The leakage of fast neutrons through all node faces on the active core surface is treated by'lbedo boundary conditions. The leakage of neutrons through all the vertical node faces on the core surface at each level is summed:

.4.1 LEAKAGE> = E.

boundary

'< k (1-a< )

k o

where J is the outward directed current (readily obtained from. the k

nodal leakage, L, and coupling coefficients) for node surfaces on the boundary at level k and a< is the corresponding nodal" surface albedo.

k The total fast neutron removals (absorption plus thermalization) can be found for level k by summing nodal values.

4.2 REHOVALSk Z Ai kLi k The fast group radial buckling is then obtained by:

4.3 (

2

)

(LEAKAGE k E + Z k Dl k

Where Eal, E r and Dl are the weighted values of fast group absorption cross section, removal cross section and diffusion coefficient for level k as defined in section 3.

Since the CORE code does not do, a detailed tracking of thermal neutrons (the 1.5 group, model is normally used for BWR's), the thermal g'oup radial buckling is not readily obtained from a CORE soluti'on.

2 Therefore, another technique is needed to determine (B )

4.2 Thermal Grou Radial Bucklin Xt is desired that the 1-D calculation utilizing, the .collapsed cross sec'tions reproduce the axial power shape and multiplication fa'ctor of the 3-D simulator case for which the cross sections were collapsed.

A method for ensuring this consistency while determining the axial variation in the thermal group radial buckling has been devised. The technique (termed the "Power-Buckling" iteration algorithm) is based on a simple variation of Crowther's power-control iteration (reference 4).

The original power-control algorithm determined the axial control material distribution to obtain a specified power shape and k-effective for a fixed set of cross sections. The current procedure determines the axial variation in the transverse (radial) buckling of the thermal group yields consistency between the collapsed 1-D cross sections'and 'hich fast group buckling) and the axial power shape and k-effective of'he 3-D simulator case from which they were obtained. This procedure does not require that a specific weighting function be used in collapsing the cross sections.

Development oF the Power-Buckling algorithm begins with the two-group diffusion equation in slab geometry; 4.4 VD (Z)Vf (Z)+[D (Z)B (Z)+E (/) ]$ (7.) =,

[vE (Z)$ (Z)+vE (Z)$ (/) ]/k 4.5 -VD2(Z) V~2 4)+[D2(l 82R(Z)+E2(Z) ]$ 2(Z) = E (Z)$ 1(Z)

Define the axial source shape:

( ).=- [ fl( )+ f2( )02( )/yl( )]yl(Z)

Define the axial power shape:

4. 7 =

~

P (Z) [ <E 1(Z)+4 E 2 (Z) ~2(Z) /~l(Z).']~1(Z) f f Define the power to source conversion factor, f(Z):

)+vE 2(Z) 42 (Z) /$ (Z) 4.8 f (Z).: f1 ( f 1 fl f2 2 ~l

,Dropping the notation for Z dependence and using 4.8 the fast'roup diffusion equation becomes:

1 1 11R 11 0 Since the values of k and 0

P (from 3-D code) are known along with all the collapsed cross sections and B,2 the only unknown in equation 4.9 is the effect of ($ 2/$ ) on f.. Hy making an approximat'ion for ($ 2.1

/$ )

equation 4.9 can be solved as a "fixed source" problem for the fast group flux.

Equating the right hand sides of equations 4.4 and 4.9 results in the following relationship between the fluxes and axial power shape:

4.10 g2

= (fP-vEfl~l)/vEf2

\

thus having solved equation 4.9 for 41, equation 4.10 can bc solved for $ 2 and both $ 1 and $ 2 will be consistent with the 3-D power shape (P).

The required thermal group radial buckling can be found'by solving equation 4.5:-

4.11 B 2 r-1 2 2 2 2R D2~2 D2 D The last term in 4.11 can be recognized:as the axial buckling so that:

412 B 2 =

2R D242 1 2 2 D2 B

2Z

'2 2 By estimating B2 equation 4.12 can be solved for an estimate of B (after equations 4.9 and 4.10 have been solved), then equation 4.5 can be solved for a more exact approximation of B 2

. Then f is recalculated by equation 4e8 based on the calculated (f '/$ ), and the solution of the equations 4.9, 4.10, 4.12,, and 4.5 repeated until converged.

In summary the power-buckling algorithm is the procedure below:

given:. all axially varying cross sections, fast group radial buckling, axial power shape, k-effective and 1-D finite difference mesh.

Determine.'he. thermal group axially varying radial buckling which yields consistency for the given data.

Procedure:

I. Make first estimate of:

( 2Z) ("/ )

(42/0 ) = ~ /~

II. Calculate iteration n+1 estimate of f(z) using iteration n estimate of ($ 2/$ 1) in equation 4.8.

III. Solve finite difference form of equation 4..9 for '4 n+1 using

~

fn+1 n+1 IV. Compute 42 using $ .,

and fn+1 in equation 4.10.

V. Use (B 2

),n. ~

$ , and $

n+1 in'.,equation 4.12 to obt'ain (B )

2"

• n+1 VI. Use (B ) and $ in finite difference form of equation n+1 4.5 .to obtain'$2 n+1 n+1

• VII. Use $ , $ , and (B ) in equation 4.12 to obtain (B2 )n+1 2Z VIII. Test co'nvergence of $ , $ , and B 2

. 'If not adequately converged repeat steps II through VII.

The values of (B2 ) determined by the above procedure guarantee that the collapsed cross sections will reproduce the k-effective and axial power shape of the 3-D'imulator case from which .they', were collapsed when used in a 1-D diffusion theory solution with the same finite difference mesh as used in the power-buckling algorithm.

5. TRANSFORMATION OF INDEPENDENT VARIABLES The models utilized for calculation of water density.and .fuel temperature in CORE hive significant differences. from those employed in RETRAN. Thus for identical (steady state) core conditions (i.e.,

power level, flow rate, pressure, inlet enthalpy and power distribution) c CORE and RETRAN s

will not necessarily arrive at the same axial distribution of, fuel temperature and water density. It should be noted that even if the models were identical the ayers e of the nodal values of fuel'ern eratures and water densitites at an axial level (computed by CORE) would.not necessarily agree with the values calculated for an ~avera e node (by RETRAN). Thus in order for the 1-D model to reliably reproduce the 3-D simulator results it is necessary to perform a transformation of the independent variables (to RETRAN 1-D model values) used to fit the collapsed 1-D cross sections rather than utilizing the CORE value.

5.1 Fuel Tem erature Transformation The "RETRAN equivalent" model for fuel temperature utilizes a single cylindrical fuel rod temperature solution for each axial level.

The average rod surface heat flux and temperature are applied as .

boundary conditions with, the fuel rod gap conductance used as a parameter.

The fuel rod average surface heat flux at any level (q s ) is obtained, from the core power level, axial power shape, number of fuel rods and outer diameter. The fuel rod surface temperature is calculated near the core inlet by the Dittus-Boelter equation:

qs 51 TDB +Tb s hDB 5.2 hDB 0.023(kk h (Re) 'Pr)

Os8 Os4

where Tb is the bulk water temperature at the axial level obtained by an energy balance on the flow channel. The remaining quantities are defined at the end of this section.

Higher into the core the fuel rod surface temperature is calculated by the Thorn equation:

T + 0.072 exp(-P/1260) q 5.3 T s

T sat s The transition from Dittus-Boelter to Thorn equation is made at the first T DB level where T s

is smaller than Ts The fuel rod is divided into a mesh of concentric nodes as used by the actual RETRAN model (typically 10 intervals in the pellet, 2 in the gap, and 4 in the clad). An analytic solution for the temperature increase from the outside of each node (T0 ) to the inside (Ti) is used (reference 5).

5.4 T ~ T + [Q(R -R )/4 + (~qQR -q R )Rn(R /R )]/k The analytic solution assumes a constant conductivity at the average interval value k = (k + k0 )/2 and.since the conductivity depends upon the temperature T an iteration is required. The heat source distribution (Q) is assumed consistent with the RETRAN model (normally flat in the pellet and zero everywhere else) ~ The conductivity of the pellet and clad are interpolated from the same tabulated data used in the actual RETRAN model and theinput value of gap conductanc'e is converted to a conductivity by multiplying by the gap width. The calculated temperatures at points in the pellet are, volume averaged to obtain the fuel temperature used in fitting the collapsed cross sections.

Nomenclature D

h

= volume equivalent hydraulic diameter k = volume liquid conductivity

h D IL

=- heat transfer coopt ic i( nt I'rom l)it tus-f)<>el tor q<<ation c

Re Reynolds- number Pr = Prandtl number qs -"average heat flux at fuel rod outer surface P = volume fluid pressure T

s

= fuel rod outer surface temperature T

b

= bulk liquid temperature I

T sat

= saturation temperature at pressure P k = average thermal conductivity in a fuel rod node Q

= fuel rod node heat generation rate q heat flux at the inside of a fuel rod node R & R o

inner and outer fuel rod node radii T & T =

0 fuel rod node inner and outer temperatures 5.2 Water Densit Transformation The "RETRAN equivalent" model for obtaining the water density at each axial level utilizes a single control volume at each level (as in the actual RETRAN,system model). The steady-state form of the RETRAN energy, momentum and continuity equations are solved for each control volume (or junction between volumes). The control volumes are used as a vertical stack with each control volume in the stack representing all the CORE nodes corresponding to, it's axial position. On option two CORE axial levels can be represented as one control volume. The solution procedure utilized by the "RETRAN equivalent" model for water density is presented in appendix A.

The RETRAN-02 thermal-hydraulic equations assume thermodynamic equilibrium between the liquid and vapor phases and thus cannot account for the effects of subcooled boiling. However, a profile fit subcooled void model is used to correct the densitites calculated by the thermal-hydraulic

models for use in determining the cross sections for the 1-D kinetics equations. To ensure consistency between the water densities used to develop the fitted cross sections and those used to evaluate them, the profile fit subcooled void correction is also made to the water density obtained from the "RETRAN equivalent" models as described in appendix A.

5.3 Control Distribution Transformation In the RETRAN l-D kinetics model it is .assumed that the effect of the control rods on cross sections can'be represented by a linear u

interpolation between cross sections for the all rods withdrawn (Z )

x c

.and all rods inserted (E x ) configurations.

The interpolation factor'Cf k, control fraction) for each cross section (x) at each level'(k) is provided as input to RETRAN, but at any xk'uantity point in time the values at each level depend only on the distance the rod bank has moved. That is, the Cf xk distribution can shift axially with time, but its shape is not changed. For purposes of representing a scram where all of the control rods move at the same speed this representation appears adequate since the relative locations of control rod tips axially would remain constant. Thus, the RETRAN representation should be adequate if the radial (X-Y plane) flux shape around control rod tips changes relatively little as the position of the tips moves.

Although RETRAN allows a different interpol'ation factor for each cross section, there does not appear to be any advantage to using a different factor for cross sections in the same energy group. Therefore, only fast f t (Cfk) and thermal (Cfk) group interpolation factors are used in the TVA model.

Given the collapsed cross sections for fully uncont'rolled (superscript u),

fully controlled (superscript c) and initial control distribution (superscript i) states for each axial level the appropriate interpolation between the uncontrolled and controlled values to reproduce the 3-D simulator axial power shape and k-effective for the initial control state can be determined by'he power-control iteration (reference 4). The power-control iteration algorithm is presented below and is very similar to the Power-Buckling algorithm presented in section .4e2 Given: The finite difference mesh, all axially varying cross sections'Z, u c x x Zx ), the initial control state axial power shape Z, (P) and k-effective (k0 ) ~ ~

I e

Determine: The axial variation in 'the interpolation factor between uncontrolled and controlled cross sections for fast and thermal groups which yields consistency for the given data.

Procedure I. Compute thermalization probability for the initial control state.

5.6 pk [E /(E .1+E ]k

~

fast group control distribution to preserve'k. 'nd E u(l-p i)-E u p

i 5.7 II. Make first estimates of:

2Z) k 5.9 ~2 ~l k r a2 k

5. 10 k k

III. Compute interpolated cross sections and power-to-source ratio using estimated control and flux ratio distributions.

5.11 E [E + Cfk(E -E ]k k

5.12 E [E + Cfk(E -E k k 5.13 f = ([v f + f2(42/yl)]/[ fl + f2'(42/41)] k IV. Solve finite differenced fast group diffusion equations:

5 14 V Dl V/1 + (E + E')$ 1 fP/k V. Compute thermal flux consistent with power:

5.15 ~2 (fP-vZflgl)/vEf2

VJ.. C<>>>>(><>L<. >>< w <':I l>>><<I <<. >I ll>< r>>>>> I <<'>>>Lr<>I <IlsI rll><>l I<>>>:

5.16 Cfk = ([E $ 1-E 2/2-D2BPZ$ 2]/[(E 2-E 2)f2])

VII. Solve finite differenced thermal group diffusion equation and compute" revised flux ratio ($ 2/$ 1).

E

5. 17 VD2V$ + E = E 2 2)2 VIII. Compute revised .estimate of thermal group axial buckling:

5.18 (B2Z k [(E $ 1/D2$ 2)-(E 2/D2)]

IX. Test convergence of buckling and control distribution if not sufficiently converged, repeat steps III through VIII.

6. SELECTION OF 3-D SIMJLATOR CASES A suggested set of standard CORE cases to utilize for producing collapsed 1-9 cross sections is shown in table 1. The subscript "a"'s'used to denote cases with. the initial rod pattern (control state) to be. used in the transient while "b" denotes completely uncontrolled (ARO) cases and "c" de'notes fully controlled, (ARI) cases. If the initial control state is all rods out then the "a". and "b" cases are the same. All 12 case conditions (36 total cases) utilize the xenon distribution calculated for the transient initial state or base case (0 a

). At each case condj.tion (0:through ll) the "b" and "c" cases utilize 'lways the. same 3-D simulator water density. (U-array) and fuel temperature (TF array) distributions as were used for the corresponding "a" case. Thus the only differences in the three cases for a given condition are due to control configurations.

The "loadline" conditions in table 1 indicate that the reactor flow, pressure and inlet temperature are calculated for the stated power level assuming the reactor reaches the power level from the base condition via a loadline. A loadline has a prescribed relationship for t'e steady-state power, flow, pressure and inlet temperature and involves no control rod movements; The "mixed" case conditions in table 1 refer to 3-D simulator cases with either the U or TF arrays (or both) taken from previously run cases which may not be at the same conditions as the current case (for example the TF array may be taken from a case at a lower power level than that used to compute the U-array). The use of "mixed" cases reduces any coherence between U an TF values for the various conditions and allows them to be truly independent variables.

Table 1 SUGGESTED CORE CASES FOR KINETICS DATA GENERATION Cases Conditions 0 a,b,c Base point (time zero power, flow, pressure, inlet temperature-and control rod configuration for the transient) 1 a,b,c Same conditions as base point except reactor pressure increased (50-60 psi).

2 a,b,c Same conditions as base point except reactor pres'sure increased by a large amount (120-140 psi).

3 a,b,c Loadline conditions for base point power increased by a small amount (approximately 10-percent NBR).

4 a,b,c , Loadline conditions for base point power decreased by a small amount (approximately 15-percent NBR) 5 a,b,c Loadline conditions for base point power increased by a large amount (approximately 20-percent NBR).

6 a,b,c Loadline conditions for base point power decreased by a large amount (approximately 35-percent NBR).

7 a,b,c Mixed case with U array from condition 1 and TF array from condition 0.

8 a,b,c Mixed case with U array from condition 2 and TF array .from condition 5.

9 a,b,c Mixed case with U array from condition 0 and TF array from C condition 6.

10 a,b,c Mixed case with U array from condition 5 and TF array from condition 6.

ll a,b,c- Mixed case with condition 5.

U array from condition 6 and TF array from

The suggested conditions in table 1 are only 'one 'set of many which would provide'n adequate basis for developing the RETRAN cross section.

The conditions in table 1 are only guidelines and mi'nor modifications are acceptable. The conditions in table 1 have been tested for base cases r

from approximately 60 percent of rated power to design po~er conditions and yield an accurate representation of the interaction of water density, fuel temperature and control over the range of values normally encountered during iJ the critical portion of pressurization transients.

For purposes of illustrating the suggested cases, table 2 shows the approximate values of power, flow, pressure and inlet temperature for a base condition at design power for the Browne Ferry units-r Since water density and fuel temperature arrays are not calculated for cases at conditions 7 through l (these cases use arrays calculated at condition 0 to 6), ~

the power level, flow, pressure and inlet temperature used are irrelevant.

~ ~

Table 2 TYPICAL CASE CONDITIONS FOR KINETICS DATA QENERATION FOR BROWNS FERRY UNITS WITH DESIQN POWER BASE POINT Thermal Total Reactor Inlet Condition Power Flow pressure tomp ..

Nabab er ~HW ~Mlb hr ~sia) ~F Descri tion 0 3441 102.5 1055.0 528.0 Base point 3441 102,5 1110.0 528.0 Pressure increase 3441 102,5 1200.0 528.0 Largo pressure increase 3 3770 115.8 1064.6 '29.8 Loadline small deorease 2947 75.9 1037 .9 525.4 Loadline small deoreaso 4100 129.2 1075.4 531.5 Loadline large inorease

'288 55.9 1028.2 520.9 Loadline large docrease

7. I'fT'I'1NC I<1?'I'RAN CROSS Sl'.(,"I'TON I'OI.YNOHIAI.S Ln RETRAN each of the required kinetics parameters for each axial neutronics region is represented as a polynomial. The polynomial independent variables are:

7.1 X =,' (U-U )/Uo 1

7.2 X = MTF ~TF 2 0 where U and TF are the water density and fuel temperature assigned to the ne'utronics region 'and the subscript "0" indicates the time zero or'ase value.

The cross sections and other parameters are then represented as:

N M N M 7.3 E = (1-Cf )[E E CnmPj, X2 ]+(Cf )[E E C nmX1 X2 ]

n=l m=1 n=l m 1 where the.,superscript u and c denote f'ully uncontrolled and fu11y controlled states; the orders of the polynomials in water density and fuel temperature "

are M-1 and N-l, respectively, and thus there are M+N fitting coefficients at each control state which are denoted by C . Typically, values of M=3 mn'nd N=2 or 3 are employed for normal pressurization transients.. The Cf values represent the control fraction to be used with cross section Eg.

The suggested cases in chapter=6 provide 36 (24 for an all.rods out initial state) values of each cross section in each region which,. are utilized to determine the fitting coefficients. Associated with each value of the cross sections are values of the* three independent variables (water density fuel temperature, and control fraction) determined by'he methods described in chapter 5. Thus for each cross section in each neutronic region, equation 7.3 is utilized for each of the cases to obtain a set of 36 equations in which the only unknowns are the fitting coefficients. Since

there are only 12-18 unknown coefficients, they can be determined by the least squares method. The 36 equations are arranged in matrix form and a standard (reference 9) least squares matrix inversion routine used to determine the fitting coefficients.

The above method simultaneously fits the cross sections for all control-states thus minimizing the error.

8. IMPLEMENTATION The calculation of the nodal cross sections described in chapter 2 is performed in a special Kinetics data Interpolation and Logical Reduction (KILR) subroutine added to the CORE code. On user option the KILR subroutine.

can be executed after any 3-D simulator calculation. It-also performs the cross section collapsing'o 1-D described in chapter 3 and determines the fast group radial buckling as described in section 4.1.

Each CORE case for which the KILR option is activated will result. in Collapsed Output for RETRAN being Produced, Summarized and Edited.(a "CORPSE" ).

The CORPSE contains the collapsed cross sections, the conditions of the 3-D simulator case (e.g., power, flow, pressure, etc.) and results of the simulator calculation (e.g., effective multiplication factor, axial power, water density and fuel temperature profiles). Since several CORPSEs are required to produce the RETRAN cross section polynomials, each CORE/KILR case places the resulting CORPSE into a MORGUE file (Multiple Outputs for RETRAN Grouped for United Evaluation).

The MORGUE file is processed by the MORTICIAN program (Multiple Operations,,

on RETRAN Transient Input Cross section Including Analysis of Neutronics),

which has several user options. The power-buckling iteration described in section 4.2 and the calculation of the "RETRAN equivalent" model water densities and fuel temperatures'escribed in chapter 5 are performed hy the EMBALM option (Equivalent Model and Buckling Algorithm results Loaded into MORGUE).

Determination of the control state interpolation factors by the Power-control iteration described in section,5.3 is performed by the CASKIT option (Control Algorithm ecifying Kinetics data Interpolation during Transients).

~S

CASKIT can also perform static 1-D diffusion theory calculations for k-effective, axial flux, and power shapes for movements of the controls rods as a bank from the base position (no feedback is used so only the effect of control on the cross sections is simulated).

The determination of the cross section fitting coefficients as described in chapter 7 is performed by MORTICIAN's FUNERAL option (Fits Used in Neutronics Evaluations by RETRAN Arrived at and Listed). The fitting coefficients along with the arrays defining the control interpolation factor's are placed in a RICH (RETRAN Input Cross, section Hierarchy) file. The FUNERAL option displays a summary of the root-mean-squared fitting error for each cross section in'ach region and on option provides a comparison of fitted and actual value at each data point.

To facilitate verification of the methods used in KILR and MORTICIAN options EMBALM, CASKIT, and FUNERAL, the GRAVE option (Generate Results of static Analyses Verifying Evaluated cross section fits) was added to MORTICIAN. GRAVE uti.lizes the "RETRAN Equivalent" model water densities and fuel temperatures for each case along with the RICH file to determine the 1-D cross sections corresponding to each of the CORE/KILR cases. The cross sections evaluated from the fitting coefficients are utilized in static 1-D diffusion theory solutions and the results compared .to the .actual 3-D simulator results stored in the MORGUE file. GRAVE also performs ad)oint flux diffusion theory calculations for the base condition cases allowing point model effective delay neutron fraction and prompt neutron lifetimes to be evaluated.

9 ~ VERIFICATION The purpose of the methods described in 'this report is to produce a set of cross section fits which allows the neutronics behavior of the reactor to be accurately reproduced when utilized in RETRAN's 1-D kinetics calculations.

Since the 3-D static simulator is the best representation of the reactor neutronics available for routine use and has been extensively compared to actual operating data, it is appropriate that we attempt to produce collapsed kinetics data which will allow the 1-D solution in RETRAN-to very closely with the 3-D simulator. If close agreement between I'gree reactivity and. power shapes between the 3-D and 1-D solutions is achieved, then the main uncertainty in the kinetics reactivity calculation will be from the uncertainties in the 3-D simulator. Thus the best available reactivity calculation will be used and any improvements in the 3-D simulator will be reflected in the kinetics calculations.

To ascertain that close agreement between the 3-D and 1-D calculations is achieved, the results of the CASKIT and GRAVE calculations of HORTICIAN

're used along with actual RETRAN results. The agreement between the 1-D diffusion solution performed by GRAVE using the fitted cross sections and the 3-D simulator cases is shown in table 3 which shows the average and standard deviation of differences in axial power shape and reactivity for the'10 cases used to develop the cross section fits at each of three control states. From the results in table 3, it is evident that the collapsing and fitting procedures produces good agreement between the 1-D and 3-D calculations.

To further test the cross section fits and the transformation of variables between CORE and RETRAN, a series of RETRAN cases was run where a perturbation was introduced (20 psi pressure increase, Q.5 o F inlet

Table 3 CO)PARISON OF 1-D TO 3-D CALCULATIONS k-eff Diff (pcm) e Arial Peak Power Diff(%)  % Eeror in Reactivity Change~a Conditions ~Av . S. Dev ~Av S .Dev ~Av S, Dev.

BF3 BOC1 all rods out 2~ 11. -0.305 1.170 1.10 0.61 all rods in -4 10. 0.001 0.292 ~ -0.13 0.91 actual rods 14. 0.614 1.302 = 2.23 - 3.94 PB2 TTl

-all rods out -63. 37. 0.337 1.220 6.75 5.93 all rods in 40. 14. -1.376 1,116 .2.31 actual rods -58. 20. -1.040 2.954 -2.94 2.44 PB2 TI2 all rods out 23 ~ 42. 0.002 1.860 -3.45 11.92 .

all rods in 5. 12. 0.083 0.111 -2.66 0.46 actual rods -3 14. -0.157 0.542 -1.49 2.33 PB2 TX2 all rods out '-43 110. -0.460 3.47 0.49 14.40 all rods in 10. 17. 0.045 0.518 0.49 3.79 actual rods -8. 13. -0.146 0.686 2.73 3.88 epcm ~ 105 1n(k2/kl)

~eReactivity change is the 'absolute value sum of component (Doppler and void) changes from the base condition at some control state.

temperature increase or 15-percent flow increase) and the transient calculation performed until the reactor reached steady state with a new power level. The final steady state power, flow, inlet temperature and pressure were then utilized in CORE to calculate the k-effective and power shape at the final conditions. By comparing the k-effective from CORE for the base case and the final conditions case, the total reactivity difference between CORE and RETRAN were determined (note both RETRAN states are critical). Utilizing the cross sections developed for conditions before the start of Peach Bottom u'nit 2, turbine trip teat 2 (EOC-2), the results of the steady state to steady-state comparisons were:

Reactivity, Error (% of Reactivity Difference cm Added b Perturbation Pressure Increase 5.7 1.6 Inlet Temperature Increase -8. 6 19.2 Core Flow Increase 7e 7 -2.16 The larger percent error in reactivity for the inlet temperature increase is primarily due to the small size of the perturbation (less than 45 pcm) since the actual error in reactivity was essentially the same for all three cases.

The axial power shapes from RETRAN and'CORE were also compared for the final steady-state conditions and were in good agreement (maximum relative power difference at. any axial region was less than 4 percent).

RETRAN's method of representing control rods was tasted to confirm its ability to predict the worth of the control rods accurately when moved, as a bank. The CASKIT option in MORTICIAN was used to calculate the static reactivity change with control bank movement (CASKIT was compared to actual

RETRAN runs, and its results were shown to be essentially identical for static rod movements). The CASKIT results were compaxed'o CORE code static 3-D calculations with control'rods moved as a bank while maintaining all other variables (xenon, water density, fuel temperature, etc.) at those for the base condition. Thus these CORE and CASKIT calculations show the effect of rod movements only. Figures 1 through 10 show the results of',these calculations utilizing cross sections developed for conditions during startup testing of Browns Ferry unit 3. Comparisons for rod movements for an initial state with several partially inserted control rods and from an initial all rods out state axe shown.

The RETRAN procedure of interpolating between fully controlled and uncontrolled cross sections to obtain cross sections for the actual xod position is shown by figures 6-10 to be very accurate for rod movements from an initial all rods out state (the worth of inserting contol rods five feet was underpredicted by the 1-D method by l.l-percent). Since bank movement of control rods from an all rods out state keeps all control rod tips at the same level as the rods move in so that essentially all nodes in a plane have the same control fraction the RETRAN linear interpolation would be expected to be accurate. Rod movements from initial rod configurations with many partially inserted control rods are not predicted as accurately. Figure 4 shows that an 8-percent underprediction of the worth of control rod movement.

to five feet compared to the 3-D simulation for a typical initial rod pattern.

Similar comparisons have been made for other reactors and conditions, and they confirm that the 1-D method accurately predicts the static worth of rod movements from an all rods out configuration and increasingly underpredicts the worth of rod movements up to five feet as the amount of initially, partially inserted control rods is increased.

Figure 1 TYI'I CRL ROD CONF I B -RODS MOVED OFT.

BF3CY.1 933NNO/l1TU RRTEO 2.0 30 K"-1.00297 CC 10 K-"1.00297 OC 1 ~ 5

)

QJ cr. 1 .0 0 ~ 5 CL 0 ~ 0 2- 4 8 10 12 14 16 18 20 22 24 RX IRL NQOE

Figure 2 TYPICRL ROD CONFIG '-RODS NOVED 3FT Q

IJJ BF3CY1 933080/t1TU RATEO C) 2.0 30 K=0.98957 (Z 10 K=0 ~ 99131 /

OC

/I 1 ~ S 1

)

UJ

/

I I

cc 1.0 4J 0 ~ 5 (Z

0.0 4 6 8 10 12 14 16 18 .- 20 22 2,4 RXIRL NODE

Figure 3 TYPICRL RQD CONE IB -RODS NQVED SET BF3t.Y1 933t1hlo/l1TU RFITED

/

2-0 90 K=0.97787

/

/

CL 10 K"0-97966

-/

OC 1-S hJ I

1.0 Lal CL O-S CC 0-0 16" 18 20 22 24 2 4 6 8 10 12 14 .

RX IRL NODE

0 Figure 4 TYI' CRL IN IT IRL .RQO'QNF I BLIRRT ION .

3-0

,BF3CY1 933MWD/tlTU RRTEO Z ~5 SO CRLC 10 CRLC Z-0 1 ~ 5 i ~ 0 0 ~5 0 ~0 0 Z '3 . 4 5 ROO 'OVEMENT ( F T.) FROl1 IN I T I RL

Figure 5

,TYP ICRL INITIRL ROD CONF IGlJRRTI ON 3.0 SF3CY1 933NhJD/MTU RRTED 3D CRLC Z.S 10 CflLC 2.0 1'5 '

~ 0 0.. 1 . -2. 3 4 5 RQD NQVEHENT t I=T ) FRQN IN I T I RL

Figure 6 IN I T IRL: RRQ CQNF I G -RQDS NQ'Ij'ED OFT BF3CY1. 933t1MO/I1TU RRTEO 2.0 30 K=1.03342 CL 10 K "1.033CZ OC 1 ~ 5

)

I cc 1.0 CL 0 5

~

CC 0 ~0 2 4 6 8 10 12 14 16 18 20 22 24 RX I RL NQOE

Figure 7 INI TIRL RRQ CONF I G.-RQOS NQVED 3FT UJ BF3CY1 933l1NO/f1TU RRTED C) 2.0 30 K"1 o02789 Ct: 10 N=l.02812 OC 1.5 UJ l

.CC 1 -0

/

/

Q /

0 ~ S CI:

0 0

~

2 4 6 8 10 1.2 14 16 18 20 22 24 RX IRL NODE 0

0 Figure 8 I NIT IRL RRQ CQNF ID.-RQDS NQVEQ SFT.

CL UJ BF3CYl 933N4IO/t1TU RRTED

.C3 2.0 30 N=l.02087 CC 10 8=1.02101

/

I 0<

1 ~ S /

/

LU /

I cc 1-0 Q

0 ~5 CC 0 0

~

2 4 6 8 10 12 14 16 18 20 22 24 RXIRL NODE

Figure 9 I N I T I RL RRQ CQNF I DURRT I QN bJ LLJ BF3CY1 933t1MO/t1TU RFITEO QL 30 CRLC

'10 CRLC

)

4J bJ LL 4

LLJ I

CZ CQ 0

0 3 4- 5 .6 . 7 RQO l1QVEt1ENT ( FT ) FROW IN I T l RL

0 Figure l0

.INI T IRL RRQ CQNF IDURRT IQN CL UJ 4.0 C3 BF3CY1 933t1WO/MTU RRTEO

. CL 3-S 30 CfllC CC 10 CRLC CC 3 0

~

2 5

~

2 0

~

1.5 1 ~ 0 0 2 3 4 5 6 7 ROO t1QVEHENT (FTi FRQl1 INITIRL

The reason for the underprediction of bank movement worths for initial configurations with partially inserted control rods is probably due to changes in the planar flux distribution near the rod tips as they are inserted. In

'a BWR, insertion of the control bank moves the control rod tips to levels with increasingly higher void fraction. Since the flux depression caused by a control rod is reduced at higher voids, the control interpolation factor (for a partially controlled region) will be too low since it was based, on the flux depression caued by the partially inserted control rods at a plane with lower void. content.

REl'L'RENCES

l. "RETRAN-02 A Program for One-Dimensional Transient Thermal-Hydraulic Analysis of Complex Fluid Flow Systems," EPRI NP-1850-CCM, Volumes I, II and III, May 1980.

2. S. L. Forkner, G. H. Meriwether and T. DE Beu, "Three-Dimensional LWR Core Simulation Methods," TVA-TR78-03A, 1978.

3.. Hsiang-Shou Cheng, et al., "A Dynamic Analysis of BWR Scram ,

Reactivity Characteristics," .BNL-NUREG-50584,,December 1976,

4. R. L. Crowther, "Extensions of the Power-Control Method for Solution of Large Inhomogeneous Reactor Problems," Proceedings of the Conference on the Application of Computing Methods to Reactor Problems, May 17-19, 1965. ANL-7050.

M. E. Garrett, "Corrections to Subroutine TEMPF in Program FCPS,"

internal memorandum CAS-360, October 13, 1978.

'VA

6. "RETRAN A Program for One-Dimensional Transient Thermal-Hydraulic Analysis of Complex Fluid Flow Systems, " EPRI NP-408, Volume 1:

Equations and Numerics, January 1977.

7. Personal Communication, J. H. McFadden, Energy Inc., July 1980.

8. "An Evaluation of State of the Art Two Velocity Two-Phase Flow Models and their Applicability to Nuclear Reactor Transient Analysis,"

EPRI NP-143, Volumes 1, 2, and 3, February 1976,

9. G. W. Westly and J. A. Watts, "The Computing Technology Center Numerical Analysis Library," CTC-39, Union Carbide Corporation Nuclear Division, 1970 Subroutine ALSQ, pages 383-389.

Appendix A Solution of Steady-State Versions of RETRAN Thermal-Hydraulic Equations Since steady-state is assumed all volume and junction total. flows are identical and the continuity equation solution is not required.

The implementation of the steady-state form of the remaining RETRAN equations is presented in the following sections.

A.l Momentum E uation The RETRAN dynamic slip momentum equation (reference 6) is solved to obtain the core pressure distribution. For volumes k-1 and k connected by junction the pressures are obtained by:

.1.A P = < APFk 1 APL Pkk-1 Q APHkk-1

"'k) where:

2.A bPHk Pk Zk/144

=

k k f~2PW 144 g k 'k"k 4.A APL = 4 k W /(144 g p A )

APM [W2/(144 g A2p ) ]f1 + V2 (1-u)p Qp A2/W2]

c k SL R gk The single-phase friction factor (f) is determined by interpolation on tabulated data based'on the volume Reynolds number. The friction factor

table was developed using the RETRAN equat1ons. The two-phase multiplier

($ 2 ) is also obtained by table lookup based on the volume pressure, mass velocity, and quality. The two-phase multiplier tables are reduced versions of the RETRAN tables for the Baroczy correlation.

Initial estimates of all volume and junction pressures are made using guessed densitites, void fractions, and qualities. All pressures are updated after the complete set of, fluid conditions has been updated.

The momentum and gravity head terms at the inlet and outlet are subtracted from the pressure difference between the upper and lower plenum volumes to obtain the core pressure drop to compare to the CORE code value. The difference in the two core pressure drops is used to adjust the outlet loss coefficient to obtain agreement. The pressure in the lower plenum is adjusted to force the core mid-plane junction pressure to be the same as the "reactor" pressure used in CORE.

A.2 Ener E uation The RETRAN energy equation at steady-state can be solved for the vapor flow rate at junction j to obtain:

6A W ~ [(qA)s s j W+Wj-1U j-1 +Wj-1Ug j-1 -WJ )/(U -U) with 7.A = h + 4 V~ /(778 g ) +

k/778 U Q Z c

8.A l18 - hg + 'z V~

J gj

/(778 gc ) + '. k /778 where (q A )

S S is the total heat added between the midpoints of volumes k-1 and k. Solution of equation 6.A is started by using

h j-1 equal to the inlet enthalpy of the CORE case (and W g

1=0). If the vapor flow rate calculated by 6.A is not positive then the fluid is still subcooled so W 8 is set to zero and the equation solved for h .

The junction void fraction is obtained from the relationship:

9.A (x. [W P~/(W (P<-P ) + P (W-VSL(l-a) P~A))]j n n-1 iteration where a and e refer to the current and previous estimates of void fraction.

The volume void fraction is calculated by:

k 1-Xf Pf f (

1 g

) H+1 k

ll.A (Xf) = Q j+1 +

(W.~l W )/W 12.A H = [Q(V 1 + V ) (VSL k] [Q(V 1

+ V.) ]

where the volume slip velocity (VSL ) k is obtained as described in the next section. The volume average density is calculated as k PR k(PE Pg)k.

For two-phase flow the liquid and vapor phase may be flowing with different velocities. The difference between the liquid and vapor velocities is called the slip velocity and RETRAN utilizes a differential equation to describe the slip velocity behavior. The slip velocity depends on the forces acting upon each phase and interphase forces.

I II<.' <<I << . ~ 'l'O'll' ll<'< <' I y I'<<>><'I I <<<<<< <

il' I<< ~ I<r<<::<<<< < >.<<':<I I <

~

<< I <<<<<I f l<<w regime. Based on references 6 and 7, the steady-state form of the RETRAN slip equation can be written as:

1 1 k 1 k 13.A VSL

~(

pp p

)

j ('k-1 'k) g k-1 k k 1 1

( p<(l-a) ) AgII, Bga p a g

where A and B are terms in the representation of interphase g g friction and are functions of the phase densi'ties, hydraulic diameter, phase velocities, viscosities and wall roughne'ss (references 7 and 8).

The volume slip velocities are obtained by averaging the junction values 14.A (V L)k 0.5(V + V )

The energy and slip equations are iterated for a junction before pro-ceeding to the next junction. After all junctions (and volumes) have been solved the pressure distribution is updated and the iteration repeated. The iterations are terminated when all pressures and void have adequately converged.

A.4 Profile Fit Subcooled Void Model A subcooled void model developed by Zolotar and Lellouche (reference 9) is used in RETRAN to calculate the water density profile used in the neutronics model. This void model relates the "true" flow quality to the equilibri'um quality, where X

e

- X [1 eD

- tanh (1 - Xe /XeD )]

1 X [1 tahn (1 Xe /XeD )

eD

X Q

= (h-h)/h f fg X "true" flow quality X

e

= flow equilibrium quality (neglecting subcooled boiling)

X eD

= flow equilibrium quality at the bubble detachment point Note: Boiling occurs if and only if Xe XeD The equilibrium quality at the detachment point is given by:

X = -C Z/h eD p fg

where, A =,4B (HDB + HN) 2H2 (H ~<i'H )+8q H (H +H )

B WD WD DB 2A Re0 ~ 662Pr k/D HDB CDB Re Pr. " k/D~

HB

= exp(P/630)/5.184(10 ~)

= 0.033 c + 0.013 CDB

= 0.2 DHy/(4R d)

CHN c = fraction of area available for flow = fL A

fL + Arod R

rod heated god radius (ft)

A fL

= flow area (ft2)

N rod Number of rods 2

wD DB sat RD B wD sat H (Tw - T )2 H (T T )2 B sat HN sat RD Calculation Procedure

1. For a given channel obtain X and e q (wall heat flux) at the center of each axial mesh.

Starting at the bottom of the channel:

2. Evaluate X eD for each node using q w

3. If Xe < X eD'.

, there is no voiding. Proceed to next node (go to step 2)

Using Xe and X eD evaluate X for rod

5. Evaluate a as given in the next section

6. Proceed to the next node at step 3 Void Relationshi Using the flow quality, X, calculated previously, the void fraction is calculated from C X+ (1-X)p /p<~ +p V /G 0 g gg3 with G total flow (lb /ft2-hr) then C

0

= L(a,P)/[k0 + (1-k )a ]

0 r = (1 + 1.57 p /p~)/1-kl) ko ~ k + (1 1

k)1 (p /p )0.25 g

~ 4

'j V C GL [

r (pq-p R

Pa g

2

) og g c

'1-a

>0 25 0)

~

sinO sin O ~ 1 for upflow of mixture

-1 for downflow of mixture o = surface tension (lb/ft) g, g ~ 4.16 x 10 ft/hr P nodal pressure P =

c critical pressure

- c c

= 1 e 1 0 L(e,P) 1-e

= 4P2/[P(P -P) ]

1 G P kl mill (kl, Kl) k =

1 1/[1.0 + exp (-Re ~

10 )]

k1 = 0.80 a

0

= a , where n = iteration number n-1'G

= 1.41 Calculation Procedure

1. Start at bottom of channel, utilize X as calculated previously and gues a

0 for each node.

2. If a 0 0.0, set L(a,P) = 1.0 temporarily

3. Calculate C and V and evaluate a gk

4. If )a-n o )

> 0.001 set a a 0 and go to step 3.

5.'f (a-a (

< 0.001 proceed to'ext node.

This void distribution is used to calculate the water density for each volume for which cross sections are evaluated. The void and water density distributions from the profile fit/drift flux calculation arc used only for the evaluation of cross sections for the 1-D neutronics and do not effect the system thermal-hydraulics.