ML19270F614
| ML19270F614 | |
| Person / Time | |
|---|---|
| Site: | Browns Ferry  |
| Issue date: | 06/01/1978 |
| From: | Beu T, Forkner S, Meriwether G TENNESSEE VALLEY AUTHORITY |
| To: | |
| Shared Package | |
| ML19270F611 | List: |
| References | |
| TVA-TR78-03, TVA-TR78-3, NUDOCS 7902280339 | |
| Download: ML19270F614 (80) | |
Text
TVA-TR78-03 TlIREE-DIMENSIONAL. LWR CORE SIMULATION METHODS S. L. Forkner G. H. Meriwether T. D. Beu 203 3
1W Topical Report Tennessee Valley Authority June 1, 1978 4
TAHl.E OF CONTENTS Page iv ABSTRACT ACKNOWLEDGMENT iv 1.
INTRODUCTION 1-1 2.
CORE DESCRIPTION 2-1 2.1 Numbering System 2-1 2.2 Physical and NucIcar Fuel Types 2-3 2.3 Orifice and Thermal-Hydraulic Types 2-3 3.
NEUTRON TRANSPORT AND CONSERVATION ?!0DELS 3-1 3.1 Neutron Energy Dependence 3-1 3.2 Neutron Tracking.
3-2 3.2.1 General Neutronics Equation 3-2 3.2.2 Fast Neutron Tracking 3-4 3.2.3 Thermal Neutron Tracking 3-5 3.3 Boundary Conditions 3-7 3.4 Evaluation of Coupling Coefficients 3-8 3.4.1 Determination of a 39 3.4.2 Determination of p and u 3-10 3.5 Iteration Strategy
. 3-13 4.
TilERMAL-HYDRAULIC MODELS 4-1 4.1 Flow Distribution 4-1 4.2
!! eat Source Distribution 4-4 4.2.1 In-channel Distribution 4-4 4.2.2 Bypass Region Distribution.
4-5 4.3 Coolant Conditions 4-5 4.3.1 Nonboiling Regime 4-6 4.3.2 Wall Voidage Regime 4-7 4.3.3 Subcooled Boiling Regime 4-8 4.3.4 Bulk Boiling Regime 4-12 4.3.5 Bypass Region Effects 4-13 4.4 Effective fuel Temperature.
4-14 i
TABLE OF CONTENTS (Cont.)
Page 5.
LATTICE PliYSICS DATA REPRESENTATION 5-1 5.1 Lattice Physics Table Sets 5-1 5.2 Corrections to Nodal Lattice Physics Parameters 5-3
- 5. 2.1 Doppler Effect 5-3 5.2.2 Moderator Temperature 5-4 5.2.3 Water Density llistory 5-5 5.2.4 Control ilistory 5-5 5.2.5 Xenon Concentration 5-6 5.2.6 Soluble. Boron Concentration 5-6 5.2.7 Combination of Correction Factors 5-7 6.
REACTIVITY DEPLETION 6-1 6-1 6.1 Forward Step Method 6.2 Power-Exposure Technique 6-3 7.
AUXILIARY CALCULATIONS 7-1 7.1 Thermal Limits 7-1 7.2 Simulated In-Core Detector Readings 7-3 7.2.1 Traversing In-Core Probes 7-3 7.2.2 Fixed Monitors 7-4 7.3 Special Fuel Management Features 7-5 7.3.1 Fuel Isotopics 7-5 7.3.2 Relative Control Rod Worths 7-5 7.3.3 Target Control Rod Patterns 7-5 7.3.4 Fuel Reload Patterns 7-7 7.3.5 Additional Capabilities 7-9 APPENDIX. Tll\\S THERMAL-ilYDRAULICS CODE A-1 A.1 Coolant Enthalpy Distribution A-1 A.2 Void Fraction and Quality A-2 A.3 Flow and Pressure Drop A-4 O
m
LIST OF FIGURES Figure Title Page 2-1 Typical LWR Planar Core Map.
2-2 2-2 Physical and Nuclear Fuel Types 2-4 3-1 Notation for Adjacent Nodes in the Neutronics Equations 3-11 3-2 Power-Leakage Iterations..
3-14 6-1 Forward Step Depletion 6-2 6-2 Power-Exposure Technique 6-4 111
ABSTRACT The CORE code was developed by 17A to simulate LWR cores over the range of condit ions encountered during operat ion and f or those analyzed during design and safety studies. The core is modeled as a three-dimensional array of rectangular nodes with coupling between neutronic and thermal-hydraulic phenomena in combination with operational para-meters.
The code can be used for detailed ccre analysis, fuel management, and optimization of operating strategies. The methods used in CORE are described in this report.
ACKNOWLEDGHENT The following people made significant contributions to the develop-ment and implementation of the methods described in this report:
T. L. Hayslett, T. F. Newton, G. W. Perry, and also H. A. Stewart who collaborated on the Appendix.
iv
1.
INTRODUCTION 1)ue to spatial variation of the moderator density, exposure, and fuel composition in a light water reactor (LWR), a three-dimensional steady-state reactor core simulator with thermal-hydraulic feedback is needed for core design analysis, fuel management, and optimization of operating strategies. Simalation of a large LWR requires modeling cf the inter-relationships of nuclear and thermal-hydraulic phenomena in combination with operational variables. The Tennessee Valley Authority (TVA) has developed the CORE code to simulate LWR cores over the range of conditions encountered during operation and for those analyzed during design and safety studies. Due to the large number of phenomena to be represented, a compromise in the detail of some of the methods used must be made to achieve a feasible computational algorithm. The methods used in CORE are described in this report.
The reactor core is represented as a three-dimensional array of rectangular nodes, each having homogeneous internal properties. The fission source distribution is determined by a nodal coupling technique based on the TRILUX method (Reference 1-1).
The neutron energy depend-ence can be represented by a modified one group model or by a model which explicitly accounts for thermal neutron interchange between nodes (Reference 1-2).
Nodal neutronic properties are determined by multi-variable Lagrangian interpolatica of input table sets of lattice physics data obtained from TVA's LATTICE code (Reference 1-3).
Nodal water density, exposure, fuel temperature, control fraction, soluble boron content, and history effects are includea.
In addition, either transient or equilibrium xenon effects can be represented.
1-1
The thermal-hydraulic model determines the flow distribution to the individual channels and a t.ngle leakage path by equalizing the pressure drops across all flow paths. The flow quality in the channel is deter-mined by Ahmad's method (Reference 1-4) such that both subcooled and bulk boiling are included in a consistent manner.
Several void quality relationships are available through options. The effect of direct heating of the in-channel and bypass water regions is also included.
Thermal limits are determined and compared to specified values. Simulated in-core detector readings can be calculated as an option.
Depletion of the core for a cycle can be performed in a series of steps or by an optional power-exposure technique (Reference 1-5).
The power-exposure algorithm calculates the Haling distribution (Reference 1-6) for use in guiding plant operations and for fuel management studies. Special fuel management options are available for the estimation of relative control rod worths, tabulation of fuel isotopic concentrations, and development of control rod patterns or fuel shuffling patterns.
The CORE code is programed in FORTRAN IV except for a few IBM 370 assembler language subroutines. Flexible dimensions allow dynamic allocation of computer storage, and a free-form input routine is used.
CORE has many output edits available and extensive restart and case stacking capabilities.
O 1-2
CIIAPTER 1 REFERENCES 1.
L. Goldstein, F. Nakache, and A. Veras, " Calculation of Fuel-Cycle Burnup and Power Distribution of Dresden-1 Reactor with the TRILUX Fuel Management Program," Trans. Am. Nucl. Soc., Vol. 10, p. 300, 1967.
2.
S. L. Forkner, "Two-Group Nodal Core Simulator Based on TRILUX Style Coupling," Trans. Am. Nutl. Soc., Vol. 23, p. 579, 1976.
3.
B. L.
Darnell, T. D. Beu, and G. W. Perry, " Methods for the Lattice Physics Analysis of LWR's," TVA-TR78-02, 1978.
4.
S. Y. Ahmad, " Axial Distribution of Bulk Temperature and Void Fraction in a lleated Channel with Inlet Subcooling," J. IIcat Tran.,
Vol. 92, pp. 595-609, 1970.
5.
R. L. Crowther, "Burnup Analysis of Large Boiling-Water Reactors,"
Proceedings of a Panel on Fuel Burnup Predictions in Thermal Reactors, IAEA, Vienna, 1968.
6.
R. K. IIaling, " Operating Strategy for Maintaining an Optimum Power Distribution Throughout Life," TID-7672, 1963.
1-3
2.
CORE DESCRIPTION The reactor core is represented as a three-dimensional mesh of rectan-gular nodes with equal spacing in the horizontal (X,Y) direction. The axial (Z) spacing can vary from level to level if desired. The horizontal dimensions of a node normally correspond to a fuel assembly or a flow channel while the vertical dimension is specified as input. Mirror and rotational symmetry options are available for representing half or quarter core geometry. The maximum number of nodes allowed is limited only by computer resources because of the use of flexible dimensioning.
2.1 NUMBERING SYSTEM Figure 2-1 illustrates the horizontal arrangement and numbering sys'em for a typical (BWR) core.
The minimum values of the X and Y
-dinates occur at the bottom left of the figure. Both X and Y values are incremented by two at each channel location and are paired X-odd, Y-even.
Locations between channels are paired X-even, Y-odd.
Therefore, if control rods and in-core detectors are located between assemblies, as in a BWR, their locations are numbered with opposite parity f rom the channel locations. This system eliminates confusion in describing loca-tions in the core If control rods and in-core detectors are associated with a single channel, as in a PWR, their locations are identified by the channel X,Y coordinates. The minimum Z value is at the bottom of the core and is incremented by one for each level up the core.
The X,Y,Z coordinates are denoted by i,j,k indices for use in CORE, each of which is incremented by one for sucecssive positions in their respective directions. The minimum i and k coordinates occur at the same positions as the minimum values of X and Z, respectively. The minimum value of j is at the same location as the maximmn Y coordinate.
2-1
l+1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 j
i 8
i i i i i
i i
i i i
i i
- 37. - - - - -l j ! -l J l
' I i
I i 36 i i 1
I l
~
~8---'--r-"-i-----I 35 - - - - - ' J L '-.;- 1 34 I i i--'-'--"-l---~~-2 33 - - - - - '
32 8
8 i
- - - - ' ' ' ' - - - - - - 3 31 - - - - -
1 30
~'----------4 29 - -
i i 28 g__
'-'-----5 26
6 25-----
24
---_-7 22 i-l-
23 - - - - -
=
8 21_ - _ -..
4.
20
g 19 - - - - -
l_
18
- - - - 10 17 - - - -
16
- - - - 11 3 3 _ _ __ _ _
14
- - - - 12 13 - - - - -
12 j j _ _ _ _ _l _
=
r
- - - - - - 13 10 1
1 i ~ - - - 14 09 -- - - - ' '
08
-_.--------15 07 - - - -
0 5 - - - - -l' l
', l - -
2-i 06
' - - - - - - - - 16 i
i 04
-~~-II 03 - - - -
L-i i
i i
02 i
i i
i i
., - - - - 18 01 - - - -- ' '
8 i
i i
i i
i i i I i i
i i i i i i i
i i
i i
i t
l I
I I I
l l
l g
g l
l l
l 1
I l
I Y
a l i I i l
l l
l l
l l
1 i
I 3
00 02 04 06 08 10 12 14 16 18 20 22 24 26 28 30 32 34 36 X+
01 03 05 07 09 11 13 15 17 19 21 23 25 27 29 31 33 35 h
CONTROL ROD LOCATION e
REAL OR PSUEDO DETECTOR LOCATION FIGURE 2-1: TYPICAL LWR PLANM CORE MAP O
2-2
2.2 PilYSICAL AND NUCLEAR FUEL TYPES Fuel assemblies with different physical dimensions or axial varia-tien in nuclear properties are modeled as dif ferent " physical fuel types."
The physical fuel type number of each channel is defined by the input array NFT...
A zero value of NFT.. indicates that the channel at loca-IJ IJ tion i,j lies outside the problem boundary.
The axial variation of lattice physics characteristics within physical fuel types is represented by defining " nuclear fuel types." Nuclear fuel type numbers are input in the array NFTK where the subscript m is the k
physical fuel type number NFT...
All channels of the same physical fuel 1J type have the same arrangement of nuclear fuel types (Figure 2-2).
Together, the NFT; and NFTK arr ys define the problem boundaries k
and the physical and nuci_ar p operties of each node.
2.3 ORIFICE AND TIIERMAL-liYDRAULIC TYPES Since scme cores have orifices built into the core support plate to improve the flow distribution, an "ori fice type" is defined for each location on the plate. Orifice type numbers are input in the NOT..
IJ array.
Since the type of orifice is not dependent on the fuel assembly, the NOT.. array is fixed irrespective of fuel assembly movements.
IJ Channels with different flow and heat transfer characteristics are defined by the possible combinations of orifice and physical fuel types and are referred to as different " thermal-hydraulic types."
The thermal-hydraulic type number of carh channel is input in t he IOT, a rray where the subscript m is the physical fuel type number NFT.. and n is the IJ ori fice type number NOT...
IJ 2-3
O PHYSICAL FUEL TYPES y3
/
2 /
f2f1/
1
,s',,,
4 1
k 2
4 NUCLEAR FUEL 2
TYPES s'
i 3
3 3
/
3 FIGURE 2-2:
PHYSICAL AND NUCLEAR FUEL TYPES O
3.
NEUTRON TRANSPORT AND CONSERVATION F10DELS 3.1 NEUTRON ENERGY DEPENDENCE The neutronic properties of each node are determined by assuming that the node is homogeneous and by utilizing the results of a detailed lat-tice physics analysis for its nuclear fuel type. These analyses (Refer-ence 3-1) account for heterogeneity effects and the spatial and energy distribution of neutrons within the node. The detailed calculations are used to determine effective homogeneous cross sections collapsed to two broad energy groups. Use of these cross sections amounts to assuming that the neutron distribution within the node is separable from the overall distribution in the reactor core.
This assumption is generally good for LWR's and allows the calculations to be separated into the detailed lattice phv.cs analyses of the nodes and the determination of the core-wide n-atron balance and distribr. tion.
In LWR's, the mean free path of thermal neutrons is small compared to the nodal dimensions. Thus, relatively few thermal neutrons escape from the node in which they are thermalized. This causes the overall distribution of neutrons in the core to be relatively insensitive to the thermal neutron tracking model. Many core simulators neglect the effects of thermal neutron interchange between nodes and use a modified one group model in which thermal neut ons are confined to the node in which they are produced. However, CORE explicitly accounts for thermal neutron interchange using a simpler model than that used for fas't neutrons.
The modified one group model is available as an option.
In cases where the assumption of homogeneous nodes may not be adequate to track thermal neutrons accurately, corrections for heterogenity 3-1
effects can be used. These corrections account for the variation in thermal ne'itron interchange across different faces of the node due tu the location of water gaps, or a control blade for INR's.
3.2 NEUTRON TRACKING The overall neutron distribution in the core is calcult.ted by con-sidering the neutronic coupling of the nodes as treated in the TRILUX model (Reference 3-2) which can be shown to be a synthetic kernel approximation to the response m _rix method (Refere ice 3-3).
Howev",
the implementation of the model is not identical to that in TRIU 3.2.1 General Neutronics Eqt.iati_on The derivation et the neutronics equations (Reference 3-4) hegin with the response matrix equation for an arbitrary energy group:
J..
= P..S.
+ 2 R
..J IJ IJ 1 n nij ni
( a'- 1 )
where J..
= the current directed from node i to node j,
IJ S.
the source in node i, r
1 P..
= the probability that source neutron will escape from a
IJ node i to node j, R.. = the response caused at the face of node i adjacent to ntj node j due te neutrons impinging upon node i f ro.a node n, and the summation is over all nodes adjacant to node i. There is a maximum of six adjacent nodes since diagonally adjacent nodes are not included.
The synt' :t ic kernel approach allows the escape and response m.itrix elements to be written as a combination of simpler probabilities. One set of probabilit.ies which can be used :o approximate P.. and R.. is the IJ Hlj 3-2
set of coupling coef ficients defined in the development of the TRII.UX equations:
o
= the probability that a source neutron in node i is absorbed f
or thermalized without escaping from node i.
The probability of escaping fr)m node i is (1 - o );
f
- r.. = the probability ' hat a neutron escaping from node i enters IJ the adjacent node j. Thus, the summation of r.. over the ij six adjacent nodes is equal to one;
- p.. = the probability that a neutron entering node j from node i 2J is directly absorbed in node j; p'.. = the probability that a neut ron entering node j from node i J
is directly reflected back to node i; p..
= 1 - p..
p..,
the probability that a neutron entering IJ IJ LJ node j f rom node i behaus like a neutron born in node j,
i.e.,
the fraction p..o.
is absorbed without escaping IJ J node j and p.. (1 - c. ) escape r. ode j ; and IJ J
../p.., which depends only on properties of node j.
a.
=
J IJ tj Using the above coupling coefficients, P.. can be written as tj P.. = (1 - c.'...
(3-2)
IJ l'
IJ and the response matrix elements are R.. = p.(1 - o.)r.. + p.6.
(3-3) nij ni i tj ni nj where 6 is the Kronecker delta which is equal to one when n = j and equal nj to zero otherwise.
Thus, the responses at the five nonincident faces of the node are equal and the response at the incident face is larger by the direct reflection probability p.. Rewriting Equation (3-1) yields nt 3-3
J.. = (1 - c. )r.j S + 1 p.J.
+p J..
(3-41 IJ 1
i i
n n1 n1 ji ja 3.2.2 Fast Neutron Tracking Equation (3-4) is the starting point for the development. of the CORE neutronics equations for the fast group. The implementation of the neutronics equations is simplified and calculational requirementa are greatly reduced by elimination of the directed nodal currents J..
from lj the neutron balance equation. The eigenfunction L.,
called the nodal 1
leakage, is used to replace the currents and is equivalent to the total current out of node i neglecting che direct reflection. terms:
L. = (1 - o.)
S. + 1 p.J (3-5) 1 1
1 n ' n1 ni Substituting Equation (3-5) into Equation (3-4) yields the neut.ron balance equation in terms of L..'
O I C..L.
+ S.
. J3 J 1
L; = j y _, g (3-6)
IJ l
J where r..(1 - p'..)
C.. = -(*1 - p..p..)
]
J-(3-7) lj 2.1 J i u.
+ 0.
A =
(3-8) i (1 + a.m./keff)(1 - c.)
1 1 1
I S. = k,gg m.A.L.
(3-9)
I t 1 1 The quantity C..
is a current operator defined such that IJ J..
J..
= C..L.
C..L.
(3-10)
IJ J1 IJ 1 J2 J O
3-4
The absorption operator A.
is defined so that A.1,. represents the fast i
1 1 neutron absorption plus thermalization rate in the node. The para ieter m.
is the nodal multiplication factor, i.e.,
the ratio of fast neutron I
production to fast absorption and thermalization.
In the modified on--
group model, the neutrons are assumed to be absorbed in the node in which they are thermalized, so m.
is identical to the nodal infinite multi-1 plication factor k,.
Iloweve r, in the two group treatment, m is corrected g
explicitly to account for net thermal neutron interchange between nodes.
represents the The eigenvalue for the problem is denoted by keff, and Sg fission source of fast neutrons.
3.2.3 Thermal Neutron Tracking Since the modified one group model assumes no interchange of thermal neutrons between nodes, the above equations are all that are required to determine the neutron distribution in the core.
However, the net exchange of thermal neutrons between nodes can have a significant effect if there is a substantial difference in the magnitude of the fast neutron source or in thermalization properties between adjacent nodes.
A thermal neutron tracking model can be developed by the same method used above for the fast group (Reference 3-4).
The response matrix equation, Equation (3-1), is written for the thermal group with the fission source S.
replaced by a thermalization 1
source Q..
1 Q. = p.A.L.
(3-11) 1 1 1 1 where p is the thermalization (or resonance escape) probability. The t
t coupling coefficients from THILUX are used to determine P..
and R..
using IJ IJ thermal group properties. Reference 3-4 shows that the only significant 3-5
response to a current of thermal neutrons entering a node is at the incident face, indicating that p..
for thermal neutrons can be assumed IJ to be zero.
Then the thermal current from node i to node j is given by JI.=(1-ch)r..Q.+pk.Jk.
(3-12) lj i
tj i ji ji Using Equation (3-12) in a thermal group neutron balance with the definitions of m; and k, yields
_qf.
= k. + A.L.
lE Ck.Q.(1-ak)-Ch.Q(1-ch)
(3-13) m.
et ji a j
ij i
i i
i t j
where % is the number of neutrons produced in fission per thermal neutron absorbed in node i, assuming that the node is homogeneous. The curient operators are evaluated from Equation (3-7) using thermal group coupling coefficients. The above equation represents the effect that thermal neutron interchange between nodes has on the fast group neut ronics. The thermal group neut ronics is coupled to the fast group through the thermalization source Q..
1 Equation (3-13) adequately uescribes the effect of therm:.1 neutron exchange between homogeneous nodes. However, the thermal current across an interface will vary depending on the orientation of the interface with raspect to heterogeneities within the node.
This effect is especially important in BWR's where fuel assemblies are separated by water gaps to accommodate the cruciform control blade located on two sides of each assembly. The effect of heterogeneities on the nodal multiplication factor can be accounted for by modifying the thermal current operators used in Equation (3-13).
O 3-6
Due to the uncertainty in evaluating the heterogeneity corrections, this effect is neglected in cu'crsat practice The modified one group model is used for 14'R's since the two-group method shows no advantage in the accuracy of results. The two group model is normally applied to PWR's.
3.3 BOUNDARY CONDITIONS The neutronics equation, Equation (3-6), requires coupling coeffi-cients and the leakage L for all six nodes surrounding each node. For nodes on the troblem boundary, one or more adjacent nodes will Jie outside the boundary and trus are undefined. The parameters for these imaginary nodes are determined in two ways.
The first of these methods is used for a symmetry boundary condition.
In this case, the node outside the boundary is assumed to have the same n'eutronic properties and leakage as the node located symmetrically to it within the problem boundary. The second method is to apply an albedo condition at the boundary. Only albedo boundary conditions are allowed in the vertical direction, but eithar albedo, symmetry, or a combination of conditions can be used on the aorizontal boi daries.
The albedo boundary condition is defined as J..
= B.J..
(3-14)
J1 1 IJ where J..
= the current from node j outside the problem boundary J1 to node i on the boundary, B
= the albedo for nertrons leaving node i toward node j, and f
J.. = the current from node i to node j.
1J 3-7
Since the albedo is equivalent to the direct reflection probability of the imaginary node outside the boundary, the current operator ft-such nodes becomes C.. = -(
- I ( 1 - B. )
r..
(3-15)
IJ 1 - B P.-)
t IJ and Equation (3-6) can be solved for node i on the boundary. Albedos for the fast and thermal groups are input for each channel in the core.
fhe albedo is set equal to zero for nodes which have all adjacent nodes within the problem boundary.
3.4 EVALUATION OF COUPLING COEFFICIENTS In order to implement the neutronics models described in the previous sections, relationships must be established between the coupling coeffi-cients and the lattice physics properties and dimensions of the node, The relationships developed here are similar to ti.ose given in Reference 3-5.
The general form of these relationships is derived from reciprocity, blacknesr. theory, and analytic diffusion theory solutions to simple problems.
In order to improve accuracy, the relationships are modi fied by inserting arbitrary constants at points where approximations might have affected the basic analytic equations. Values for the constants found by minimizing the deviation of the core multiplication factor were and nodal powers for a series of problems analyzed using CORE from the results obtained by fine mesh, two group, finite-difference diffusion theory calculations. There are other procedures for determining the nodal coupling coefficients; a comparison of several different methods is presented in Reference 3-6.
O 3-8
3.4.1 Determination of a The nonleakage probability o for node i is related to the blackness f
b by the approximate formula j
- b. = 1.E.(1 - c.)
(3-16) i 1 1 1
where i = the mean chord length of the node, and I = the absorption plus thermalization cross section for fast neutrons or the absorption cross section for the thermal group.
The ratio $ of the flux at the node surface to the average flux in the 1
node is also related to blackness:
2-b.*
Q. = 1.1.
(3-17) i i1 2b j Combining Equations (3-16) and (3-17) yields g (2$*. + 1.1.) - 2 g
(3-18) o =
g)(2$f+I2) g1 where g was introduced to account for higher order chord length terms 3
neglected in Equation (3-16).
The values of $ and g) are calculated from 1
$ = 1 + lx C; + C M +
K [C3* 4 M)]
(3-19)
M+
K (C5* 6 2
8 =C7+C8 (3-20) 1 where K = 1/D, D = the diffusion coefficient, M = the relative difference between the nodal mean chord length and the mean chord length of a cubic node of equal volume, and C = the optimized coefficients.
3-9
3.4.2 Determination o_f p d a In order to satisfy the reciprocity theorem, the following relation-ship must exist between the coupling coefficients:
a.
+ 0.
= p*..
F..
(3-21) r..
'3 l
'3
- 2. 5.(1-0.)
1 1 1
The parameter F..
is the fraction of the total surface area of node i IJ shared with node j.
In practice, r..
is assumed equal to F..
since the IJ 1.1 remaining factor is approximately equal to one.
A second relationship between the coupling coefficients is obtained by comparing the nodal neutronics equation to an analytic diffusion theory solution for a large homogeneous multiplying medium. The neutronics equation using the notation of Figure 3-1 can be written as O
2r 2r y
y 1+p i
X+ ' 'X-) #
1+p
( 'i -
Y+ -
Y-X y
2r (m. - 1)(a. + o.)
1+p (L. - L + - LZ-) -
(3-22) t Z
Z (1 + a.m.')(1 - c.)
1 1 1
where m; = m /kd f, the ratio of the nodal multiplication factor to the g
ef fective multiplication factor of the problem.
The asymptotic diffusion theory solution about node i in the R direction is I'R +
R-i R
+
cos(B where AR is the nodal spacing and B is the buckling in the R direction.
R 2
The cosine is approximated as a series neglecting terms higher than BR and is substituted into Equation (3-22) to obtain 3-10
\\
>z+
/
Y+
G-----
i
/
x+
x-
/
Y-Z-
FIGURE 3-1: NOTATION FOR ADJACENT NODES IN THE NEUTRONICS EQUATIONS 3-11
r (AX)2 r (AY)~
r (AZ)2 2
X y
7 y,y g
(3-24) 1 + p) 1+p 1+p T+o.
Y Z
t I + um.
I where M is the migration area of the node.
The critical condition m.*
=1 (3-25) 1 + M.28.
I 1
was used to eliminate the buckling from Equation (3-24).
Equations (3-21) and (3-24) are combined using the relationship 1 =
p.. + p..
+ p..
13 13 IJ
=
p..
+ p..(1 + a.)
(3-26)
IJ ij i
to arrive at a relationship between p and a:
1+a 1-P
-A-Rj (3 W )
P
~]~[Rj Rj 1+a 3
I t am j with CD C
9 10 PRj *5
- lg (3-28)
J J J where C and C are p mize ee icien s.
ecm n n
qu ns g
10 (3-21) and (3-24) also gives the following quadratic equation for a:
a(1+a) + b(1+a) + c = 0 (3-29) where a, b, and c are dependent on P ' * '
and nodal dimensions.
R In the thermal group, the coupling coefficient p..
is zero and the IJ relationship between p' and o :
parameter o is not used.
This leads to a g
O 3-12
p!. = 1 g 1.I!(1 - ci)
(3-30)
IJ 21 1 1
where the superscript t denotes the thermal group. The parameter g I'
2 a correction factor to account for the approximation made in Equation (3-16):
33+(C12 + D )
13\\
(3-31)
=C g3 E
where the C re optimized coefficients.
n Only four of the six coupling coefficients originally defined are r
,p
, a, and a. Of these, required to solve the neutronics equation:
g g
r.. is approximated as being proportional to nodal dimensions. Therefore, IJ only three coefficients are dependent on lattice physics parameters and are determined from the equations listed in the following table.
Table 3-1 Equations for Coupling Coefficients Fast _ Group Thermal Group o.
3-18 3-18 L
p..
3-27 3-30 IJ a.
3 c' not used 1
3.5 ITERATION STRATEGY The neutronics iterations form the innermost loop of the power-leakage iteration strategy shown in Figure 3-2.
The thermal-hydraulics calculations performed in the outer iteration provide the nodal operating conditions required to evaluate the lattice physics parameters at each 3-13
INI'UT ESTIMATE k gg, POWER DISTRIBUTION, AND PRESSURE DROP o
CALCULATE WATER ADJUST PRESSURE DENSITY AND FLOW DISTRIBUTIONS DROP n
if EVALUATE NODAL I;EUTRONIC PARAMETERS o
CALCULATE SOURCE AND LEAL \\GE DISTRIBUTIONS n
LFAKAGE
^
NO DISTRillUTION AND l.EAKAGE DISTRIBliTION k,g g CONVERGED?
AND k 77 h
LEAKAGE LOOP YES u
CALCULATE POWER DISTRIBUTION o
POWER
^
NO DISTRIBUTION POWER CONVERGED?
DISTRIBUTION POWER LOOP
\\
YES OUTPUT O
FIGURE 3-2:
POWER-LEAKACl: ITERA110NS l-14
node.
Thus the nodal coupling coefficients, multiplication factor mg, and the current and absorption operators C.. and A. in Equation (3-6) tj i
remain constant during any set of Icakage iterations. However, they are updated in the outer loop due to thermal-hydraulic effects.
A successive over-relaxation procedure is used to update the leakage distribution between iterations. The sequence of calculations in the inner loop is shown below.
S"i' A.I"-1 (3-32)
=
I **
k"'ff e
3 6
3"' C..L" + 39~I C..L'."I
+
L'.'=3{I
.I 3'3 J'3 (3-33) i IC..+A.
lj 1
J t
L.
LS~I (3-34)
A" -
1 n-1 1
g n-1 n L" = L"-
+6 3
(3-35) 1 1
1 where the superscript n is the iteration counter and 6 is the over-relaxation parameter. The leakage distribution for the first set of inner iterations is estimated as either a flat or input distribution.
The initial value of k is set equal to one as is the normalization df factor A.
As k and the leakage distribution converge, the value of eff A will approach unity. The normalization factor is determined for each iteration as the Rayleigh quotient of successive L estimates:
g I (L;)2 A" = '
(3-36)
I L7 L"'I 1
1 3-15
For each iteration, the effective multiplic:.uon factor of the core is determined as follows:
l I k"~ff S"~
. e 2
4 k"-
(3-37) 1 C.. L" ~ I + I A. L" ~
B
'd
- i k"ff = 1 (k"~ff + A"k )
(3-38) e 2
e The parameter k is the altiplication factor obtained by a neutron bal-ance using the leakage distribution from iteration n-l.
The first term of the denominator of Equation (3-37) represents the leakage out of the core; thus, the summation is over only those nodal interfaces with albedo boundary conditions. Equation (3-38) is used to obtain a better estimate of the multiplication factor f or each iterat ion; when the calculation and k"gf will be equal.
converges, k The leakage iterations, are repeated unt il both h nd the ler! gr eff distribution converge withir. input criteria or until a specified maximum number of iterations have been performed.
In a two group calculation, the nodal multiplication factors are then recalculated from Equation (3-13).
The fission power distribution is determined from the source distribution:
P. = E
.S.
(3-39) 1 R1 1 where E is the energy released per fission neutron produced, obtained g
from the lattice physics calculations. The new P and two group m are j
g passed to the outer loop for the thermal-hydraulics calculations and the process repeated until the flow and power distributions and kd f " '"
converged.
t 3-16
CHAPTER 3 REFERENCES 1.
B. L. Darnell, T. D. Beu, and G. W. Perry, " Methods for the Lattice Ehysics Analysis of LWR's," TVA-TR78-02, 1978.
2.
L. Goldstein, F. Nakache, and A. Veras, " Calculation of Fuel-Cycle Burnup and Power Distribution of Dresden-1 Reactor with the TRILUX Fuel Management Program," Trans. Am. Nucl. Soc., Vol. 10, p. 300, 1967.
'. Weiss, "Some Basic Properties of the Response Matrix Equations,"
3.
t Nucl. Sci. & Eng., Vol. 63, pp. 457-492, 1977.
4.
S. L. Forkner, "Two-Group Nodal Core himulator Based on TRILUX Style Coupling," Trans. Am. Nucl. Soc., Vol. 23, p. 579, 1976.
5.
Z. Weiss, " Nodal Equations Derived from Invariant Imbedding Theory,"
Nucl. Sci. & Eng., Vol. 48, pp. 233-247, 1972.
6.
A. Ancona, et al
" Nodal Coupling by Response Matrix Principles,"
Nucl. Sci. & Eng., Vol. 64, pp. 405-417, 1977.
3-17
4.
THERMAL-HYDRAULIC MODELS Because of the dependence of the nodal lattice physics data on nodal operating parameters, the thermal-hydraulic model of the core has a significant effect on the calculation of the reactivity and power distribution. Another important neutronic coupling to the thermal-hydraulic calculation is through the Dappler effect on reactivity with its dependence on fuel temperature. Thus, the thermal-hydraulic analysis may be considered of equal importance as the neutronics analysis in a core simulator.
Formulas from the 1967 ASME Steam Tables (Reference 4-1) are used to evaluate the required water and steam properties as a function of system pressure and inlet enthalpi l.
4.1 FLOW DISTRIBUTION For a specified total core flow rate, the flow rate through each channel and a single leakage path is determined by equalizing the pressure drops across all flow paths. The pressure drop is determined by contri-butions from elevation, friction, acceleration, and local losses which are influenced by the channel thermal-hydraulic type, coolant density, and crud thickness. The flow through the channel is also a function of the total channel power, the axial power shape, the inlet subcooling, and the reactor pressure.
To avoid performing repetitive calculations of pressure drop and flow for the large number of channels in a LWR core, the detailed thermal-hydraulic calculations are performed externally to the CORE code.
The THAS single-channel code developed by TVA does the thermal-hydraulic calculations for each channel thermal-hydraulic type using a specified 4-1
symmetric axial power Jistribution, inlet subcooling and reactor pressure.
These calculations are repeated for several channel powers at varying flow rates to determine the pressure drop.
The results are combined into two-dimensional tables of flow rate versus channel power and pressure drop for use in CORE,
'lhe methods used in TIIAS are described in the Appendix.
The flow in individual channels is determined by Lagrangian interpola-tion in the base table for the channel therrnal-hydraulic tyi e.
The base channel flow W is then corrected for increased flow resistance due to crud o
buildup, and for variations of the reactor pressure, inlet subcooling, and axial power shape from the assumed table values. TilAS calculations are made to determine the effect that varying the above parameters has on the iIow rate The actual channel flow W is least squares fit using the TilAS results to a quadratic equation for use in the code:
9 aW' + bW
+c=0 (4-1) c c
The coefficients in Equat. ion (4-1) are Co+f C1 + X.in(C2 + x. C )I in 3 J M-2)
P in P
11 b=I-(C4+C5 fV) 11--+C X.
(4-3) p 6 in P in fg c
c=p
-p
- W (b
+aW)
(4-4) in in o o oo where the o denotes quantities evaluated at base conditions, f = a crud thickness factor, x.
= the inlet quality, in p.
= the inlet density, In 4-2
11, = the heat of vaporization, P = the channel power, C =
the fitting constants, and V = a channel power shape parameter K
Z I P OZ l ~ ~k k
k f
_ k=1 (4-5)
K I P aZ q
k k=1 with P = the relative power in node k, AZ = the length of the node, k
Z = the distance from the bottom of the active core k
to the center of node k, L = the total length of the fuel rods, and K = the number of axial levels in the core.
The core leakage flow is computed by setting the pressure drop in the bypass region equal to the sum of the elevation head and bypass region entrance losses. For the purpose of calculating the leakage flow rate, it is assumed that the coolant mixes uniformly in the bypass region; thus, all bypass regions are lumped together.
In determining the channel and leakage flow rates, the value of the core pressure drop is assumed to be known. Normally, the pressure drop is not known a priori and must be considered as a function of the flow rate as well as the parameters affecting the flow rate.
The pressure drop is usually estimated initially and adjusted between thermal-hydraulic iterations by forcing the sum of the calculated channel and leakage flow rates to equal the specified value of the total core flow.
4-3
4.2 IIEAT SOURCE DISTRIBUTION The heat addition to the coolant flowing through the core is the sum of three components. Most of the heat deposited in the coolant is transferred from the fuel through the clad. Smaller fractions are added directly to in-channel and bypass water by neutron t hermalizat. ion or gamma heating in the coolant or in structural material in the channel and then transferred to the coolant. The heat flux distribution is proportional to the fission power distribution calculated by the neutronics niodel.
The direct heating components are determined by the energy deposition distribu-tion for fast neutrons and gamma rays which tends to be smoother than the fission power distribution. The distribution of fast neutrons and gammas in the core is estimated by
$ jk y + (1-y)P..k (4-6)
=
i tj where y is an input parameter which is equal to one for a uniform distribu-tion and zero for a distribution ident.ical to the relative fission power distribution P..k.
tj 4.2.1 In-channel Distribution The contribution of direct heating of in-channel water to the total core heat rate is d
CF P Q... d..k (4-7)
Q
=
ij k d n lj K IJ where Y = Ow framon of core power from M rect lwating of g
in-channel water, 4-4
<l
= the direct heating energy deposition distribution, which ij k is proportional to the local water density, and C = a proportionality constant.
Thus, the in-channel heat addition distribution is the sum of its components:
Q'ij k = q"ij k ij k + Qij k S
(4-8) is the heat flux out of the clad, and 5 is the surface area where q';'p of the clad in the node.
4.2.2 Bypass Region Distribution The contribution to the core heat rate from direct heating of the bypass region water is F
h (4-9) ij k b n ij k where F is the fraction of power from direct heating of the bypass region b
water, and depends on whether the adjacent nodes are controlled or uncon-trolled. Since the variation of water density in the bypass region is small, its effect on the energy deposition distribution is neglected, and d..k is lj set equal to one.
In addition, the fission power P used in Equation (4-6) to evaluate Q in the bypass region must be a weighted average of the fis-sion powers of the nodes surrounding the bypass region between fuel assemblies.
4.3 COOLANT CONDITIONS The neutronics culculations are coupled to the thermal-hydraulic analysis through the dependence of lattice physics parameters on local operating conditions. Especially important is the effect water density 4-5
has on reactivity and the power distribution. The bulk coolant temperature in the channel is also required to calculate the effective fuel temperature (Section 4.4) for the evaluation of the Doppler effect on reactivity.
In normal operation of a INH, subcooled coolant enters at the bottom of the core and it, heated as it moves up the channel. This can lead to a maximum of four flow regimes with characteristic coolant conditions:
a) the nonhoiling regime where the coolant remains subcooled, b) the wall voidage regime where steam bubbles form but do not detach from the fuel rod surface, c) the subcooled boiling regime where the steam bubbles detach and may recondense in the subcooled liquid phase, and d) the bulk boiling regime where the liquid is at the saturation temperature and the steam bubbles are not recondensed.
The following sections present the methods used to calculate the bulk coolant temperature and density in each flow regime. The effect of steam formation in the bypass region between fuel assemblies in BWR's is also discussed.
4.3.1 Nonboiling Regime In ti.is flow regime, the coolant is subcooled and the surface temperature of the fuel rods is below the saturation temperature. No vapor is formed and thermal equilibrium exists. Thus, the average bulk coolant temperature in the node is calculated by an energy balance Q
T = T; +
(4-10) b fp where 4-6
T = the average bulk coolant temperature in the previous node, g
Q = the heat addition to the in-channel coolant between the midpoints of the nodes, G = the mass flux through the node, A = the flow area of the node, and g
C = the specific heat of the coolant.
p The coolant density is calculated from the 1967 ASME formulas as a function of bulk coolant temperature and system pressure.
4.3.2 Wall Voidage Regime In the wall voidage regime, steam bubbles are formed on the heating surface, but do not detach from the wall. Heat transfer from the wall is enhanced by agitation of the coolant layer adjacent to the wall caused by the growth and subsequent collapse of bubbles on the surface. There is a negligible net generation of vapor; thus, all the heat goes to raise the bulk liquid temperature and thermal equilibrium can be assumed.
The bulk coolant temperature is calculated from Equation (4-10) as for the nonboiling flow region. The coolant density is computed from the equation U=p
- a(p
-p)
(4-11) g f
where p = the density of saturated water, f
p = the density of saturated steam, and a = the void fraction in the node.
The void fraction is calculated from Griffith's model discussed in Reference 4-2:
4-7
a=6h (4-12)
'f where P is the heated perimeter and 6 is the bubble layer thickness given by h
q" - h(T -T ) Pr k w b bb 6=
(4-13)
- 1. 07is (T - T )
g b
The single phase heat transfer coefficient is h = 0.030 Re Pr (4~l4) b b
and the wall temperature is calculated from the Jens and Lottes equation:
-6)II' 60 (q"x10 T
=T
+
(4-15) e /900 P
w s
Other parameters appearing in Equations (4-13), (4-14), and (4-15) are a) the surface heat flux q",
b) the bulk coolant thermal conductivity k 'b c) the bulk coolant Praniti asi. Reynolds numbers Pr and Re b
b' d) the saturation temperature T '
s e) the equivalent diameter of the channel D '
e f) the reactor pressure P-4.3.3 Subcooled Boiling Regime The subcooled boiling regime is defined as that between the onset of bubble detachment and the beginning of bulk boiling. Vapor bubbles detach from the heating surface and can eventually recondense in-stream in the subcooled liquid.
Thus, heat transfer is enhanced due to the mixing of the steam bubbles in the subcooled liquid before collapse. Since vapor and subcooled liquid exist simultaneously, the assumption of thermal equilibrium O
4-8
is not applicable and the calculation of coolant conditions becomes complex.
Ilowever, the bulk temperature and density can be determined if the quality is known and a slip model for the liquid and vapor phases is specified. A modified version of Ahmad's method (ib.crence 4-3) is used to calculate coolant conditions in the subcooled boiling regime.
In subcooled boiling, the total heat flux is divided into two components q"=q"+q{
(4-16) where q'jistheportionthat raises the temperature of the liquid phase andq{isthepartwhichproducesboiling. Ahmad assumed that q}= HAT (4-17) g where AT = T - T is ne local subcooling, and develop d a correlauon g
b for the subcooled boiling liquid phase heat transfer coefficient:
I!
k 1/2 1/3 11.
11 hp = 2.44 Re P#
b b
it il
(~
where li is the inlet enthalpy and li is the enthalpy of the saturated g
g liquid. The axial level in the channel where bubble detachment occurs can be determined by finding the point Z where the subcooling calculated using d
Equation (4-10) is less than the subcooling required for detachment AT d
obtained from Equation (4-17).
Cases where AT exceeds the inlet subcooling d
AT signify that subcooled boiling occurs from the start of the heated in length, and AT is set equal to AT with Z set equ 1 to zero.
d g
d Ahmad wrote an energy balance on s length of heated channel in subcooled boiling using several simplifying assumptions to obtain a differential equation for the bulk coolant temperature in ditensionless form:
4-9
0 )=-
3+fC AT (4-19)
C 2
dZ The dimensionless subcooling AT is the ratio of the local subcooling to the subcooling at bubble detachment. The dimensionless distance Z is Z-Z d
M-20 Z =
7sb where Z is the axial distance up the core and Z is the total significant sb boiling length up the channel f rom Z. The heat transfer parameter C 3
1 and the condensation parameter C2 "'"
hPZ g h ab
(~
C :
3 GA C fp Kil Z
_2 fg_ sl>
(4-22)
C2 3 GA C g
where K is a condensation constant proportional to the square of the heat flux. For practical problems, the condensation of steam bubbles in sub-cooled boiling is very slow due to the high flow rate through the channel; thus, the condensation parameter C is usually very small and the second 2
term in Equation (4-19) can normally be neglected.
Assuming a uniform heat flux, Equation (4-19) can be solved to give
- 3/2)}
(4-23)
AT = exp
- (C Z
+CZ 3
2 The heat flux in a LWR normally will not be uniform over the entire channel, but the assumption of a uniform heat flux inside each node allows the heat flux to vary axially in steps and is consistent with the neutronics model.
O 4-10
Thus, the bulk coolant temperature in each node can be determined from Equation (4-23) using the outlet conditions of the lower node to evaluate the required parameters:
T
=T
- AT AT (4-24) b s
d The average coolant density in subcooled boiling is found using Equation (4-11) with the void fraction X
a=
(4-25)
P E
X+P S(1-X) g where S is the slip ratio determined below. The flow quality obtained from Ahmad's method is AZ - 1 + AT X=
g, B + AT where Q
A = GA C AT
( 27) fp d 11 B=
(4-28) p d The quantity Q in Equation (4-27) is the total heat addition to the in-c channel coolant, including the direct heating contributions, between Z and d
the position of interest, and can be calculated using Equation (4-8).
- Thus, the effect of nonuniform heat flux is also included in the calculation of quality.
It should be noted that the Levy equation used in TilAS for flow quality (Reference 4-4) is an approximation to Equation (4-26) for the special case where ATd < ATin 2
and in-stream bubble condensation is negligible (C = 0).
4-11
Development of a general model for the slip ratio is difficult due to its dependence on the flow, the physical properties of the coolant, and the channel geometry. Most slip modele are hased upon empirical correlations of data taken under a liinit ed range of conditions. Several slip models are available as opt. ions, but the CISE model (Reference 4-5) is normally used. The CISE slip ratio is obtained as follows:
ah Y = 1-u (4-29) h Y
'l = FC }l
-C
~
3 l
1.0 T<0 S=
(4-31) 1.0+C,[l' T>0 where o is the void f raction assuraing homogeneous flow [S = 1 in Equation (4-25)] and the C's are empi rical parameters
-0.51 p
-0.08 g
C = 0.0273 Re Wb -
(4-32) 3 p
-0.19 p
.22 C4 = 1.578 Re
-g (4-33) p 8
where Wb is the Weber number.
4.3.4 Bulk Boiling Regime The bulk boiling flow regime begins where the bulk coolant reaches the saturation temperature. The total heat addition is utilized in the production of vapor; thus, the bulk coolant temperature remains constant at the saturation temperature The average in-channel coolant density is determined using the same method as for the subcooled boiling r egime 4-12
4.3.5 B pass Region Ef fects y
The bypass water between fuel assemblies in a INR can experience boil-ing which may significantly affect the lattice physics properties of the adjacent nodes through a reduction in moderator density. This effect is accounted for by using an ef fective nodal water 6ensity in the neutronics calculations:
U
=U+a (p -p)C7+
(4-34) eff g
where U is the in-cha:mel water density, a is de vo M fracMon in de bp bypass region asso;iate.. with the node, and the C's are fitting constants which depend on the channel thermal-hydraulic type.
Data for the fits is obtained from detailed lattice physics calculations by varying in-channel and hypass water densities.
The void fraction in the bypass region is calculated using Equation (4-25) with the CISE slip model.
The required bypass quality is deter-mined from an energy balance assuming thermal equilibrium between t'te liquid and vapor phases:
Ob7 + lig 7
- 11 (4-35)
X
- ~
bp II fg where Q is the total heat addit. ion rate to the bypass water up to the mid-b point of the node.
The flow rate W in the bypass region surrounding each channel is assumed to be the same and thus can be calculated by dividing the total leakage flow by the number of channels in the core.
Since the rate of heat addition in the bypass region is low, the assumption of thermal equilibrium does not lead to significant error.
4-13
4.4 EFFECTIVE FUEL TEMPERATURE Tae neu' ronic a..d thermal-hydraulic analyses are coupled through the dependence of the Doppler reactivity effect on fuel temperature.
Since the Doppler correction is applied on a nodal basis, a method of calculating the effective fuel temperature in each node is required.
By using a detailed fuel rod thermal-mechanical code, the temperature distribution in individual pin segments can be determined as a function of power and exposure for each fuel design. The ef fective fuel temperature for pin i is calculated as a weighted average of the volume averaged and surface temperatures of the fuel:
(1-a)T*
(4-36)
T* = aT*
+
f avg sur where a is the weighting factor, typically about 0.85 (Reference 4-6).
To arrive at the proper value of the nodal effective fuel temperature, the variation of power and exposure from pin to pin within the node must be taken into account. The individual pin powers and exposures can be determined from the nodal average values with the aid of the pin-by pin relative power distribution calculated by TVA's I.ATTICE physics code (Reference 4-7).
Since exposure is directly proportional to the time integral of power, the pin-by pin relative exposure distribution in the node can be obtained by integrating t he power distributions from a series of LATTICE calculations at different nodal average exposures.
Thus, the pin powers and exposures can be calculated by multiplying the nodal aver-age value by the appropriate distribution. The individual pin effective f uel temperatures are calculated f rom Equation (4-36) and averaged to obtain the nodal value of the fuel temperature.
O 4-14
The results of the above calculations are reduced to the following equation for use in the core simulator:
Tg=Tb + q' C3 + E(C2+C3 4* (5+ 6 )J
('~
)
where T = the bulk coolant temperature in the node, b
q' = the nodal average linear heat generation rate, E = the nodal average exposure, and C = fitting coefficients.
4-15
CllAPTER 4 REFERENCES 1.
C. A. Meyer, et al., " Thermodynamic and Transport Properties of Steam," American Society of Meenanical Engineers, 1967.
2.
G. W. Maurer, "A Method of Predicting Steady-State Boiling Vapor Fractions in Reactor Coolant Channels," WAPD-BT-18, 1960.
3.
S. Y. Ahmad, " Axial Distribution of Bulk Temperature and Void Fraction in a IIeated Channel With Inlet Subcooling," J. lieat Tran., Vol. 92, pp. 595-609, 1970.
4.
S. Levy, " Forced Convection Subcooled Boiling - Prediction of Vapor Volumetric Fraction," GEAP-5157, 1966.
5.
A. Prenoli, et al., "An Empirical Correlation for Evaluating Two-Phase Mixture Density Under Adiabatic Conditions," An Informal Contribution to the Private Meeting, !! cat Transfer Laboratory, Grenoble, France, 1970.
6.
D. II. Risher, Jr., "An Evaluation of the Rod Ejection Accident in W PWR's Using Spatial Kinetics Methods," WCAP-7588, 1975.
7.
B. L. Darnell, T. D. Heu, and G. W. Perry, " Methods for the Lattice Physics Analysis of LWR's," TVA-TR78-02, 1978.
O O
4-16
5.
LATTICE PllYSICS DATA REPRESENTATION The neutronics and thermal-hydraulic analyses are coupled through the dependence of the lattice physics parameters on the operating condi-tions at each node. The lattice physics data required for the neutronics calculations is obtained from detailed analyses of each nuclear fuel type represented by a single fuel assembly and its surrounding water (Refer-ence 5-1).
These analyses are performed parametrically as a function of in-channel water density, assembly averaga exposure, and exposure-history dependent in-channel water density, with the control rod either withdrawn or inserted during the analysis. The spatially averaged two group lattice physics parameters are represented by polynomial fits and by Lagrangian interpolation of tabulated values for each nuclear fuel type.
5.1 LATTICE PHYSICS TABLE SETS Most of the lattice physics parameters are functions of the in-channel water density U, exposure E, and the control fraction CF (the fraction of the node height a control rod is inserted). Since the dependence on these variables tends to be nonseparable, the basic method for represent-ing lattice physics data is in sets of tables at specific values of U, E, and CF.
Multi-variable Lagrang,ian interpolation is used to obtain nodal parameters at the required values of U, E, and CF.
The following para-meters are represented by the U-E-CF table sets:
a) the infinite multiplication factor k,,
b) the fast fission factor c, c) the thermalization or resonance escape probability p, d) the age-to-thermal I,
5-1
e) the fast absorption plus removal cross section 1,
9 f) the thermal absorption cross section I,,
g) the thermal diffusion area L,
h) the xenon thermal microscopic absorption cross section X'
i) the xenon reactivity worth p '
X j) the reactivity worth p "E
C inserted versus depletion with the rod withdrawn, k) the energy release per iission neutron produced E 'R
- 1) the fissile and fertile fuel isotopic concentrations, and m) the reactivity ef fect of water density history Ull for BWR's or the soluble boron concentration PPM for PWR's.
Due to variation in operating parameters, the in-channel water density will not be constant nor can control rod positions he maintained throughout the cycle. These ef fects are accounted for by defining exposure history dependent variables. The water density history is defined as Ull(E) = E U(E)
(5-1) where U is the exposure averaged in-channel water density E
1 tl(E) =
E O
The control history Cli is defined similarly using the exposure averaged control fraction CF.
For lattice physics parameters a-k above, the table sets assume that the node was depleted at the instantaneous in-channel water density O
5-2
with the control rod withdrawn throughout the depletion, i.e., U = U and CF = 0.
Depletion at a constant soluble boron concentration PPM, is also assumed for PWR's.
Corrections made for deviation of U, CF, and PPM froni the base values are discussed in Section 5.2.
Since the nodal multiplication factor for BWR's is very sensitive to the water density history, the correction to k,due to deviation of U from U is included in the lattice physics tables for greater accuracy (item m).
For PWR's, the correction to k,due to variation of PPM from PPM is tabulated as 9
Fuel isotopic concentrations (item 1) are item m for the same reason.
related to history dependent rather than instantaneous values of water density and control fraction. Thus, the isotopics are interpolated f rom the tables as a function of exposure using U and CF instead.of U and CF.
5.2 CORRECTIONS TO NODAL LATTICE PHYSICS PARAMETERS The lattice physics table sets assume base values or fixed relation-si' i ps for nodal independent variables that are not as closely coupled as are U, E, and CF.
These variables include the fuel temperature, the moderator temperature, the equilibrium xenon concentration corresponding to the base power density, and the soluble baron concentration. Many of the lattice physics parameters are not sensitive to deviation of such variables from th assumed base values.
In these cases, the parameter is determined solely by interpolation in the table sets.
Other parameters, most notably the nodal multiplication factor k,, must be evaluated by correction of the interpolated values to account for the variation of independent variables from values assumed in t'e tables.
a 5.2.1 Doppler Effect Deviation of the nodal fuel temperature T from the base value T 7
affects the nodal multiplication factor k, through Doppler broadening 5-3
of resonances in the fuel. This effect is accounted for by a correction to the base k,from the lattice physics tables:
CD "D
f ~
~
where the Doppler coefficient a I8 "E*"
D and E.
In addition to correcting k,, C is also used to make correcuons D
to other parameters affected by changes in resonance absorption, includ-ing the age-to-thermal, the absorption plus removal cross section, and the thermalization probability. The correct ion factor for l and p is g
simply the ratio of the Doppler corrected k, to the base k,, the correction factor for I is the inverse ratio.
5.2.2 Moderator Temperature e
The lattice physics tables assume a fixed relationship between the moderator temperature and the in-channel water density at a base reactor pressure. Subcooled boiling or variations irom the nominal reactor pressure can cause the nodal moderator temperature calculated in the thermal-hydraulic analysis to differ from the base temperature. The correction to the base k,due to this effect is the factor:
C = 1+ f (T
-T)
(5-4)
M M m m
where T, = t he nodal mode rat or temperat u re,
= tte base moderator temperature, and T
4 f
=a fitting function which is ilependent on U, E, CF, T,
g and the soluble boron concent rat '
.n.
O 5-4 s
5.2.3 Water Density Ifistory The lattice physics tables for items a-k assume that the exposure averaged and instantaneous water densities are equal. This assumption may be inadequate for EWR's where the nodal in-channel water density can vary greatly depending on operating conditions. Thus, for BWR's the effect of water density history on the nodal k, is accounted for by multiplying the base k,by a correction factor obtained by interpolation in lattice physics tables (item m).
These tables contain the ratio of k,at a specified U to the base k, for three values of U.
Quadratic interpolation between the three values determines the water density history effect for any nodal value of U.
Since PWR's do not experience as large a variation in water density, the correction factor for the nodal k, is expressed as a polynomial fit:
CF (5-5)
C
= 1 + E(U-U) fmil(1-CF) + fGI2 gg The fitting functions f and f re dependent n E and (U-U).
Ulfl Ulf2 5.2.4 Control }{istory The lattice physics table sets assume that the node is depleted with the control rod always fully withdrawn from the node.
To account for the effect of previous control rod insertions, the nodal multiplica-tion factor is multiplied by a factor dependent on the exposure averaged control fraction:
C
=1+p CF (5-6)
C C
where p is the reactivity worth if depleted in a fully controlled C
condition (CF = 1) as opposed to the base case of CF = 0, and is determined by interpolation in a table set (item j).
5-5
5.2.5 Xenon Concentration The poisoning ef fect of the highly absorbing fission product Xe-135 must be included in determining the nodal k,.
The base k, interpolated from the tables (item a) incorporates the effect of the equilibrium xenon concentration at the base power density assumed in the table. The effect of the actual nodal xenon concentration on k,is accounted for by a correction factor related to the fractional variation f rom the base concentration:
( X-X \\)
(5-7)
C
- l~f P
X X X X
o o
where X = the noda! xenon concentration, X = the equilibrium xenon concentrath n corresponding to the 9
base power density,
" ttvity worth of the base xenon concentration, p
- me n g
o interpolated from the tables (item i), and f
=a fitting function that corrects p f r dev ati n f the X
Xo soluble boron concentration from the base value.
The actual xenon concentration can be input from a previous case or set equal to zero, or can be calculated as the equilibrium value for the actual power density in the node or from the transient iodine and xenon equations.
5.2.6 Soluble Boron Concentration The lattice physics table sets for items a-k for PWR's assume a soluble boron concentration ppt 1 equal to the base value PPfl. Since g
the nodal k,is sensitive to the soluble boron concentration, a 5-6
co r rect ion factor to account for variations of PPM from PPM, is inter-polated from the table set for item m for PWR's as was done for the water density history effect for BWR's.
Soluble boron is not used in IWR's except in an auxiliary shutdown system; t hu s. _ the dependence of BWR nodal k,on PPM is expressed as a polynomial correction factor to the base k,.
C
= 1+fppg (NM - WM,)
(5 W ppg The fitting function f depends on U, E, CF, PPM, and the relative ppg xenon concentration.
Since the soluble boron affects the spectrum in the node, the cross sections are also affected by a deviation of PPM from PPM. The xenon 9
cross section is modified as follows to account for the spectral shift:
(5-9) where K is a fitting constant. The effect on other cross sections is lumped into a correction to the nodal macroscopic thermal absorption cross section descr; bed in the following subsection.
5.2.7 Combination of Correction Factors The nodal multiplication factor k,is calculated from the base value k" by combining the correction factors discussed in the previous subsections:
k, = k (1 + C }
(5-10)
D M UH C X PPM Since both xenon and soluble boron compete with other materials in the node for thermal neutrons, the thermal group constants I, and L must be corrected for deviation of xenon and soluble boron concent rations from 5-7
the base values. The correction Iactors for these two parameters are proportional to the relative change in k, caused by the variation of both xenon anti soluble boron.
O O
5-8
CilAPTER 5 REFERENCE 1.
B. L. Dat' ell, T. D. Beu, and G. W. Perry, " Methods for the Lattice Physics Analysis of LWl's," TVA-TR78-02, 1978.
5-9
6.
REACTIVITY DEPI.ETION 6.I FORWARD STEP METii0D The forward step method is used to simulate depletion with core operat-ing conditions that change during the cycle. The core average exposure is incremented in steps during which the operating parameters (control rod pattern, core power, flow rate, etc.) experience only minor variations.
Figure 6-1 diagrams the sequence of calculations for the forward step method of depletion.
Using representative operating parameters for the step and the nodal fuel properties (exposure E, water density history Of and control history Cli) calculated at the end of the previous step, the core multiplication factor k and relative power distribution P are determined by the power-g leakage iterations. A restart file containing the power distribution, nodal fuel properties and other information is output for use in later cases. The fuel properties can be updated for the next step:
I E" k = E".k + P9.k AE" (6-1) ij tj tj W
Uti?U = Ull?.k + U".k P".k AE" Vf (6-2) tjk tj ij tj I
CIf?
= CII". k + CF". k P".k AE" W (6-3)
IJk tj tj tj ij k where U
he eHecuve nodal water denshy, ij k CF..k = the control fraction of the node, IJ W
= the mass of heavy metal per unit length in the node, and Vf = the core average mass of heavy metal per unit length.
6-1
O I NPl!T o
POWER-LEAKAGE ITERATIONS (FIGURE 3-2)
If WRITE RESTART FILE v
UPDATE NO FUEL PROPERTIES TO END OF STEP?
YES UPDATE NODAL FUEL PROPERTIES (E, Ull, Cil)
ANOTilER YES EXPOSURI:
STEP?
NO
,g ANOTilER YES CA;iE?
NO u
STOP O
FIGURE 6-1:
FORWARD STEP DEPLETION 6-2
In this manner, the fuel can be depleted step by step through the cycle.
6.2 POWER-EXPOSURE TECilNIQUE An alternate procedure for the depletion of a core through a cycle is the Power-Exposure Technique (Reference 6-1).
This technique (PET) directly calculates the unique, consistent power and exposure distributions at the The method can also be used end-of-cycle (EOC) given the desired EOC keff.
to step from a known initial state to another exposure level with a specified k
The power distribution determined by PET is that of the Ifaling eff.
Principle (Reference 6-2) which, if maintained, gives minimum power peaking throughout the cycle. The llaling distribution is useful for many survey calculations and in the development of target control rod patterns (Section 7.3.1).
The sequence of the PET calculations is shown in Figure 6-2.
For an inputk;OC F
and a known initial state, the exposure step AE and the lialing fg power distribution E are estimated. The EOC or final state nodal fuel properties are updated for iteration n+1:
E" k = E!9 + E9.k AE" WT (6-4) lj tjk tj ij k UlI9
= Ull'. 9 + U". k 57.k AE" h"f (6-5) tj k tjk tj tj ij k Cll" k = Cll!9 + CF".k P".k AE" hT (6-6)
IJ tj k tj tj where the superscript BOC indicates the beginning-of-cycle or initial state values and the other parameters are defined in the previous section. The nodal fuel properties are used in the power-leakage iteration which 6-3
INPUT
~
o SPECIFY k eff GUESS AE, P ijk o
ESTIMATE EOC 1
NODAL FUEL PROPERTIES
-v
' ' ~ ' ' ^ "
INTR /.PolATI:
11ERATIONS p
lj k (FIGURE 3-2) o ADJUST
!:0 n+1 EOC,,
AL eff eff' YES
/
"+
NO
-P
-P
=
ij k ij k *
/
YES WRITE RESTART FILE o
AN0111ER YES CASE?
NO STOP FIGURE 6-2 : l'0WER-EXPOSURE 11.CilNIQU E 6-4
re-estimates k and P..k.
The procedure is repeated with an adjusted eff tj AE and extrapolated P..k until the power distribution converges and the lj calculated k is equal to the desired value.
g.f 6-5
CHAPTER 6 REFERENCES 1.
R. L. Crowther, "Burnup Analysis of Large Boiling Water Reactors,"
Proceedings of a Panel on Fuel Burnup Predictions in Thermal Reactors, IAEA, Vienna, 1968.
2.
R. K. Ifaling, " Operating Strategy for Maintaining an Optimum Power Distribution Throughout Life," TID-7672, 1963.
O O
6-6
7.
AUXlLIARY CAI.CULATIONS 7.1 THERMAL LIMITS Core thermal performance characteristics can be determined using the converged power distribution (Section 3.5).
The thermal limits evaluated consist of peaking factors, linear heat generation rates (LHGR's), and margins with respect to transition boiling. Of course, the core simu-lator is not intended for monitoring thermal limits, a function of the plant process computer, but for comparison and predictive purposes.
In order to evaluate thermal margin, it is useful to cxamine the following peaking factors for each node or channel:
a) the nodal local peaking factor F is the ratio of the power of the g
highest powered rod segment in the node to the average rod segment power in that node, b) the channel peaking factor F is the ratio of the power of the channel to the average channel power in the core, c) the nodal axial peaking factor F, is the ratio of the power of the node to the average nodal power in its channel, d) the nodal gross peaking factor F =FF is the ratio of the power g
ca of the node to the average nodal power in the core, and e) the nodal total peaking factor F =FF is the ratio of the power g
of the highest powered rod segment in the node to the average rod segment power in the core.
The local peaking factor is evaluated using input data developed from the detailed lattice physics analysis for each nuclear fuel type. The repre-sentation reflects the influence of the configuration of the fuel and control rods around the node as well as the nodal exposure and water density. The 7-1
other peaking factors are determined from the relative power distri-bution. All the peaking factors are edited for the nodal positions carresponding to maximum F, maximum F, and minimum transition boiling margin.
Utilizing the peaking factors, core thermal power, and assembly design data, the rod segment LHGR and the average planar LHGR can be calculated and compared to limiting values.
LWR's are restricted in operation to conditions that will not lead to the transition frem nucleate to film boiling. Two indicators of thermal margin with respect to transition boiling are used:
the departure from nucleate boiling ratio (DNBR) for PWR's and the critical power ratio (CPR) fo r BWR ' s. Minimum values for these parameters are determined and compared to limiting values in order to evaluate the transition boiling margin.
The DNBR is the ratio of the DNB heat flux to the operating heat flux of the highest powered rod segment in a particular node.
The DNB heat flux is that required to cause transition boiling at any location in the node and is calculated as an option using the W-3 correlation (Reference 7-1).
The CPR is the ratio of the channel power required to initiate transi-tion boiling at any location in the channel to the operating channel power.
Generally, either the critical heat flux or quality at the point of transition is related to power to create a critical pocer correlation. The correlations used inclede CISE (Refarence 7-2), XN-2 (Reference 7-3) and an assembly power limit method (Reference 7-4).
The CISE formulation uses the critical quality and the XN-2 correlation calculates the c ritical heat flux for con-version to the critical power.
The assembly power limit method is similar to that used in TVA's Browns Ferry Nuclear Plant process computer and in 1-2
the auxilisry method described in Reference 7-4.
A precalculated base critical power, provided as a function of nuclear fuel type, channel flow and exposure, is adjusted to account for the effects of axial peaking, total core flow rate, and inlet subcooling.
7.2 SIMULATED IN-CORE DETECTOR READINGS INR's incorporate fission chamber detectors at fixed radial locat$ons as axially traversing in-core probes (TIP) or fixed monit,rs to measure operating power distributions. These detectors are located in guide thimbles in the bypass region between fuel assemblies for BWR's and within the assembly itself for PWR's.
7.2.1 Traversing In-Core Probes Since the detectors are fission chambers, the response at any location is proportional to the fission rate of the detector at that location:
( -1) fk ifk ifk I
where R
= the TIP reading at radial location 2 and axial level k, h
P
= the power density in node i adjacent to the detector at ak position 2,k;
= the ratio of the fission rate in the detector to the Ek power density of node i, determined from a correla-tion of detailed lattice physics data as a function of the nodal water density, exposure, and control fraction for each nuclear fuel type, V = the volume of the detector, C = a proportionality constant, and 7-3
the summation is over the I nodes directly adjacent to the detector. Thus the possible variation in tr, clear fuel types around the detector is taken into account.
Simulated TIP readings are made for each real or pseudo detector location (Figure 2-1) at each axial level.
Then the readings are normalized to those for a specified detector thimble:
R I<f k =c [
(7-2) n Rm where n
R
= the normalized TIP reading at position E,k; g
R
= the average reading for the specified detector thimble which is located at the radial position m, and C
= the input normalization value of E _
7.2.2 Fixed rionitors Af ter the normalized TIP readings are calculated, simulated signals for iixed in-core monitors can he generated for BWR's.
Local power range monitors (LPRf!'s) are located at fixed axial levels throughout the core in order to continuously monitor the power distribution. Average power range monitors (APRf!'s) combine signals from several LPRtt's to estimate the total power generated in the core Rod block monitors (RBfl's) are used to restrict control rod withdrawal based on the power change in the region affected by the rod.
The required RBtl signal is obtained by averaging the signals of the LPRF1's adjacent to the control rod.
Thus only LPRr! signals need to be evaluated.
The LPRfl signals are determined using the input detector locations and the calculated TIP readings.
I f t he LPRr1 is located on an axial level boundary, then it s signal is the.iverage of the TIP reailings for t he nodes 7-4
above and below the boundary. Otherwise, the LPRM reading is set equal to the TIP reading at the detector location. The APRM and RBM signals can be calculated as the average of the appropriate LPRM signals.
7.3 SPECIAL FUEL MANAGEMENT FEATURES 7.3.1 Fuel Isotopics s
The fissile and fertile fuel isotopic concentrations can be determined
(
1 if the fuel properties (exposure, water density history, and control history) are available from a previous CORE run or from the plant process computer.
Nodal values of E, U and CF are used to interpolate the isotopics from input table sets described in Section 5.1.
The mass of each of the following materials is tabulated on a nodal basis:
U-235, U-236, U-238, Pu-239, Pu-240, Pu-241, Pu-242, fissile Pu and total Pu.
Summations of nodal values are made to obtain the isotopic content by physical fuel type, by assembly, in specific core regions, or in the entire core.
7.3.2 Relative Control Rod Worths In order to calculate the core shutdown margin, it is necessary to determine the highest worth control rod.
A simple perturbation technique is used to estimate the worth of each control rod relative to the other rods by a weighted change in multiplication factor of the affected aodes when the control rod is withdrawn. The shutdown margin is then determined by performing detailed calculations using the core simulator for several of the rods having the highest relative worth.
7.3.3 Target Control Rod Patterns Target control rod patterns are used to guide plant operations. As an aid in developing these patterns, CORE contains an algorithm (Refer-7-5) which estimates the rod pattern requited to approximate a ence 7-5
specified power distribut. ion, normally the Ifaling distribution calculated with the power-exposure technique (Section 6.2).
This option is valid only for those reactors that use control rods for power shaping, i.e.,
BWR's.
The algorithm is based on the following guidelines:
a) the control rod pattern should minimize the deviation of the k,
e distribution from that associated with the desired power shape; b) the control rods are grouped into two basic checkerboard patterns with half the rods in each pattern. Only one pattern can b used at any time, Each pattern is subdivided into diagonally alternat-ing deep and shallow rods; c) the control rods are withdrawn in a sequence such that the distances between inserted rods and the distances between inserted rods and the core boundary are optimized to reduce interaction between control rods; and d) all symmetrically located rods are inserted the same distance into the core in order to maintain core symmetry.
The control fraction required in each node ijk to obtain the desired d
power distribution for a specified core k is determined from gg d
k k"
"ff E
k
=
=
ij k kp ijk CF..k -
(7-3) 1J k"
- k" m
e ij k ij k where k" = the nodal multiplication factor with the control rod withdrawn, k[=thenodalk,withthecontrolrodinserted, O
7-6
P = the nodal k, associated with the desired power distribution, k
and kI
= the core k associated with the power distribution.
7 gg The number of control rods at axial level k corresponding to the required control distribution is estimated by I
J N =1 I I
CF..k (7-4) k 4
ij i=7 j=7 where the division by 4 is necessary since each control rod controls four fuel assemblies in a INR.
Since control rods are inserted from the bottom of the core, the number of rods inserted at level k cannot be larger than at level k-1.
The dis-tribution calculated by Equation (7-4) may not satisfy this requirement.
Therefore, a pattern search procedure (Reference 7-6) is used to find the monotonically decreasing distribution with the minimum least squared devia-tion from N. The revised distribution is used to establish the number of control rods inserted at each axial level. The individual rod groups to be inserted are selected based on the guidelines presented earlier. Control rod patterns developed by this method can be evaluated by performing a normal simulator calculation.
7.3.4 Fuel Reload Patterns An algorithm is available to assist in the development of a fuel reload pattern for a cycle. The algorithm uses random assignment and exhaustive enumeration of exchanges in fuel assembly positions to generate trial reload patterns. Quarter core symmetry with both automatic and user-specified shuffling instructions is used to reduce the number of feasible patterns to 7-7
be searched. The deviation from a desired power distribution is employed as the objective function which is estimated af ter each assembly movement based upon the local distribuon of assembly infinite multiplication factors. The desired distribution is that which yields maximum reactivity within a constraint on power peaking and is obtained from an auxiliary pattern search routine (Reference 7-6).
The user-specified shuffling instructions restrict the possible assign-ments of the fuel assemblies:
a) discharge specified assemblies, b) insert fresh assemblies into specified locations, c) move selected assemblies into specified locations, and d) restrict the assignment of selected assemblies to within a set of specified locations.
O The automatic instructions have lower priority than those specified by the user and, in order of implementation, include a) discharge the fuel assemblies with the lowest ictivities to make up the total required number of assemblies to be discharged, b) randomly insert the required number of fresh fuel assemblics into available locations, and c) randomly assign assemblies to available locations.
The above instructions are incorporated in an assignment table that controls the movements of the assemblies. The table is dimensioned for a quarter core and indicates the permissible assignments for each it I assembly.
If th" problem i:, a hal f or full core case,.issemblies with similar tuel properties to those in the modeled quadrant are assigned to symmetric locations in the rema ining quadrants.
7-8
The number of fuel assemblies to be shuffled simultaneously cars be varied from two to the total in the quarter core.
Reload patterns are generated until a limit on the number of assembly movements is reached. The pattern with the minimum va;ue of the overall objective function is designated as the best pattern. This optimization tecanique does not guarantee that a global optimum pattern will be reached. However, a good fuel reload pattern can be obtained through the use of a good initial fuel pattern and a logical s.et of shuffling instructions. The ganerated pattern can be refined by the user as necessary and evaluated by a detailed simulator calculation.
7.3.5 Additional Capabilities The CORE code also has the capability to automa;ically adjust a specified parameter during the power-leakage iterations to achieve an input value of the effective multiplication factor. Parameters available for adjustment include the tot.al core flow, thermal power, inlet subcooling, and soluble boron ecntentration.
In order to calculate reactivity coefficients and adiabatic reactivity worths, various thermal-hydraulic propcities can be held constant at values f rom a previous case while performing the power-leakage iterations. The combinations of parameters include a) the water density, moderator temperature, and effective fuel temperature, b) the water density and moderat or temperature, or c) the effective fuel temperature.
"'le nodal xenon concentration can also - ueld constant with any of the above combinations.
7-9
CHAPTER 7 REFERENCES 1.
L. S. Tong, " Prediction of Departure from Nucleate Boiling for an Axially Non-Uniform lieat Flux Distribution," J. Nucl. Energy, Vol. 21, pp. 241-248, 1967.
2.
L S. Tong, " Boiling Crisis and Critical lleat Flux," TID-25887,1972.
3.
K. Galbraith and J. Jacch, "The XN-2 Critical Power Correlation,"
XN-75-34-A, 1975.
4.
V. W. Mills, "TBAR Method for LPRM Calibration and Core Performanca Evaluation for Browns Ferry Unit 3 Cycle 1," NEDC-21366, 19/6.
5.
T. L. llayslett, Jr., "MYROD, A Computer Program for Predicting Control Rod Patterns for a Boiling Water Reactor," Master's Thesis, Department of Nuclear Engineering, Mississippi State University, 1976.
6.
R. Ilooke and T. A. Jeeves, " Direct Search Solution of Numerical and Sta i stica l l'roblems," J. Assoc, for Computing Mach., Vol. 8, pp. 219-229 1961.
O O
7-10
APPENDIX TilAS TlfERMAL-llVDRAULICS CODE The TilAS code developed by TVA performs the thermal-hydr.iulic analysis of a single heated channel for the range of operating conditions in LWR cores. The equations of conservation of energy and momentum are applied to control volumes defined by dividing the channel into a number of egnal length segments. Nonboiling, subcooled boiling and bulk boiling conditions are treated. Experimentally based correlations are used to calculate the flow quality, the steani void f raction, the distributed friction loss coefficient, and local pressure losses. Water and steam properties are determined from the 1967 ASME formulas (Reference A-1) except for a least squares fit of the surface tension.
A.1 COOLANT ENTIIALPY DISTRIBUTION Conservation of energy in each control volume implies FP ll
= 11; +Cj (A-1) where 11
= the coolant enthalpy at the volume outlet, II. = the coolant enthalpy at the volume inlet, 1 11 W = the mass flow rate through the channel, P = the channel power, F = the fraction of power in the control volume, and C = a conversion factor.
The above eiguation is valid for all flow regimes.
/
A-1
A.2 VOID FRACT!.h AND QUALITY The f r a c t i c,..
of steam voids in the coolant is calculated from X
a=
(A-2) 9 E
X+O S(1-X) f where X = the flow quality, S = the slip ratio evaluated using the CISE model discussed in Section 4.3.3, p = the density of saturated steam, and p = the density of saturated water.
g Flow quality is defined as the ratio of the flow rate of steam to the total flow rate in the channel:
8 X=
(A-3)
The flow quality in the bulk boiling regime is equal to the thermodynamic quality X,,
H - li g (A-4)
X=X
- H
-f e
g f
where H = the local coolant enthalpy, li = the enthalpy of saturated water, and g
li = the enthalpy of saturated steam.
In subrooled boiling, thermodynamic equilibrium does not exist and the flow quality is calculated from the f.evy correlation (Reference A-2).
Levy's method begins by finding the subcooling at bubble detachment AT
=T
-T as a functi n 1 a dimensionless distante Y cha ract erist ic d
b g
A-2
of the bubble layer thickness:
I I
[ AT,p - W YpB B
AT
- SQ Prg + In I + Pr
-1 55Y I
^~
AT
=
g B
d g+In(1+SPr)+ fin Y 3 0
( AT
- SQ Pr g
B The dimensionless distance Y is derived from a force balance on a bubble:
B 0.015 [0g D p (A-6)
Y =
B p
ce2 g
where pg = the viscosity of water, p = the density of water, g
o = the surface tension of the bubble, and D = the equivalent diameter of the channel.
The single phase subcooling AT, is given by
=b (A-7)
AT h,p sp where q" is the heat flux out of the surface and h, is the single phase heat transfer coefficient calculated from the Dittus-Boelter correlation h
= 0.023
- Reg" Pr,.4 (A-8) o The thermal conductivity of water is denoted by k, Re and Pr are the g
g g
Reynolds and Frandt1 numbers, respectively, evaluated at the bulk liquid temperature. The parameter Q in Equation (A-5) is defined as 9"
Q=
(A-9)
C./t gpg where C is the specific heat of water. The wall shear stress is P
A-3
G I
=f (A-10)
W 8x P c g where G is the mass flux through the cont rol volume and f is the Darcy-Weisbach friction factor f = aRe
+c (A-II) with the constants a, b, and c supplied by the user.
Then the quality at hubble detachment is calculated C AT (A-12) x =-
d fg and the subcooled flow quality is interpolated by 0
X $Xd p
(A-12)
X=
l X
exp ( {
'g g
l X
>X X, - Xd e
d 9
A.3 FLOW AND PRESSURE DROP Conservation of momentum in the control volume is expressed as a pressure loss between the inlet and outlet to the volume. The total pressure drop across a control volume is the sum of components from f riction, elevation, acceleration, and local ef fects:
AP = APf + APh + AP + AP (A-14) a s
At the user's option, the flow can be fixed to determine the pressure drop, or can be changed iteratively until the channel pressure drop reaches a preset value.
The f riction pressure drop is calculated f rom the Darcy-Weisbach equation (A-15)
AP = f4 - 2g p f
fD p
A-4
where I = the Darcy-Welsbach single phase friction factor,
$ = the two phase multipler of the f riction f actor, and 7
AZ = the length of the control volume.
The two phase friction multiplier is interpolated at
(.4ction of mass flux and quality in tables developed from Baroczy's correlation (Reference A-3).
The elevation pressure drop is calculated from ap + (1-a)p A7. 8-(A-16)
AP =
h g
2 gc where g is the acceleration of gravity.
The pressure loss due to acceleration of the coolant as it is heated is 2
AP =-
op v 2 + (1-a)p v g ',
gg Ef out a
(A-17) op v 2 + (1-a)p v 1
2 E
g8 EI in c
where
= the velocity of the vapor phase v
=5 G a
v = the velocity of the liquid phase f
= (1-)() g (1-a)
The local pressure losses are caused by changes in flow area along the channel, such as spacers:
2 AP. = K$
(A-18) i f 2g P c1 A-5
The single phase local loss coefficient K is an input parameter, and the two phase locallossmultiplier$fusedisamodifiedRomiemultiplier (Reference A-4):
(1-X)2 2p L f 1
(A-19) 4 = (1-a,), u,p 2
where ah (A-20) a =
C. + a (1-C.}
t h
i The parameter a is the void fraction assuming homogeneous flow, and C is h
f an empirical constant which depends on the flow restriction. Note that this pressure drop is omitted for control volumes in which no spacer is located.
The inlet pressure to the first control volume is determined by subtracting the inlet pressure loss from the input lower plenum pressure.
Then the above calculations are performed successively for each cont rol volume using the outlet conditions of the lower volume as the inlet condi-tions for the next. The channel outlet pressure is obtained by subtracting the exit pressure loss from the outlet pressure of the uppermost control volume.
O A-6
APPENDIX REFERENCES 1.
C. A. Meyer, et al., " Thermodynamic and Transport Properties of Steam," American Society of Mechanical Engineers, 1967.
2.
S. Levy, " Forced Convection Subcooled Boiling - Prediction of Vapor Volumetric Fraction," GEAP-5157, 1966.
3.
C. J. Baroczy, "A Systematic Correlation for Two-Phase Pressuie Drop," Chem. Eng. Prog. Symp. Ser. *-
64, Vol. 62, 1964.
4.
J. A. Wooley, "Three-Dimensional BWR Core 5 mulator," NEDO-20953A, 1977.
A-7