ML19270F613
| ML19270F613 | |
| Person / Time | |
|---|---|
| Site: | Browns Ferry  |
| Issue date: | 04/15/1978 |
| From: | Beu T, Darnell B, Perry G TENNESSEE VALLEY AUTHORITY |
| To: | |
| Shared Package | |
| ML19270F611 | List: |
| References | |
| TVA-TR78-02, TVA-TR78-2, NUDOCS 7902280330 | |
| Download: ML19270F613 (150) | |
Text
TVA-TR78-02 METIIODS FOR TIIE LATTICE PHYSICS ANALYSIS OF INR's B. L. Darnell T. D. Beu G. W. Perry 94 Topical Report Tennessee Valley Authority April 15, 1978
TABLE OF CONTENTS Pagg ABSTRACT . . . . . . . . . . . . . . . . iv ACKNOWLEDGMENT . . . . . . . . . . . . iv
- 1. INTRODUCTION . . . . . . . . . . . . 1-1
- 2. CROSS SECTION GENERATION . . . . . . . . . 2-1 2.1 Cross Sect ion Libraries . . . . . . . . . . 2-1 2.1.1 Epithermal Data . . . . . . . 2-1 2.1.2 Thermal Data . . . . . . . . . . . . . . 2-2 2.2 Dancoff Factors . . . . . . . . . . . . . 2-5 2.2.1 Theory . . . . . . . . . . 2-5 2.2.2 PWR Method . . . . . . . . . . . . 2-7 2.2.3 BWR Method . . . . . . . . . . . 2-8 2.3 Resonance Calculations . . . . . . . 2-11 2.3.1 Resolved Resonances . . . . . . . . 2-12 2.3.2 Unresolved Resonances . . . . . . 2-18 2.4 Epithermal Spectrum Calculations . . . . . . . . . 2-22 2.5 Thermal Spectrum Calculation . . . . . . 2-27 2.5.1 Theory . . . . . . 2-27 2.5.2 Method of Soluticn . . . . . . 2-31 2.5.3 Geometry . . . . . . . . 2-34 2.5.4 Thermal Cross Sections . . . 2-39 2.6 Effective Control Rod Parameters . . . . . 2-44 2.6.1 Control Blades . . . . . . . . . . . 2-44 2.6.2 RCC Rods . . . . . . 2-55 2.7 Fission Product Treatment . . . . 2-65
- 3. TWO-DIMENSIONAL FLUX AND POWER DISTRIBUTION . . . 3-1 3.1 Two-Dimensional Multigroup Diffusion Equations . . 3-1 3.2 Method of Solution . . . . . . . 3-7 3.3 Relative Power Distribution . . . 3-9 3.4 Mesh Description . . . 3-10 i
TABLE OF CONTENTS (Cent.)
O Page
- 4. DEPLETION CALCULATIONS . . . . . . . . . . . . . . . . . . 4-1 4.1 Isotopic Depletion Equations . . . . . . . . . . 4-1 4.2 Method of Solution . . . . . . . . . . . . . . . . 4-2 4.3 !ormalization . . . . . . . . . . . . . . . . . . . . 4-7 O
O il
LIST OF FIGURES Figure Title Page 2-1 Geometry for Calculation of Moderator Blackness . . 2-6 2-2 Fuel Pin in a Square Lattice 2-10 2-3 Coordinates for the Integral Transport Equation . . . 2-28 2-4 Definition of Currents . . . . . 2-32 2-5 BWR Lattice Configuration . 2-35 2-6 Cylindricized Fuel Cell . . . . . . 2-36 2-7 Super Cell Model . . . . . . 2-36 2-8 Cylindricized Assembly (8x8 BWR) . . . 2-38 2-9 Two-Dimensional Mesh for a Super Cell . . 2-42 2-10 Geometry of Poison Pins in a Control Blade . 2-45 2-11 RCC Rod Cell . . 2-56 2-12 Blackness of Cylindrical Absorbers . . 2-59 2-13 Transport Correction Factor for Cylindrical Absorbers 2-60 2-14 R versus R r RCC Rods . 2-64 0 a O 3-1 Mesh for the Five-Point Difference Equation . 3-2 3-2 Model for a BWR Assembly . . . . . 3-11 3-3 Model for 1/4 PWR Assembly . 3-12 4-1 Depletion Chains Represented . 4-3 ili
ABSTRACT The methods currently in use at TVA for the lattice physics analysis of LWR's are described in this report. A chapter is devoted to each of the major calculations: the generation of neutron spectra and cross sections, the two-dimensional fev group diffusion theory computation of the flux and power distribution in the assembly, and the change in isotopic concentrations with exposure.
ACKNOWLEDGt1ENT Appreciation is extended to Dr. S. R. Specker for his considerable contributions to the development and application of the methods presented in this report. Other contributors include C. H. Chen, S. L. Forkner, and D. L. liut s on .
O iv
- 1. INTRODUCTION This topical report describes the lattice physics methods in use at the Tennessee Valley Authority. These methods are embodied in a computer program called LATTICE, which TVA has been actively developing for several years. LATTICE consists of several components which perform the various calculations: the computation of the neutron spectra and cross sections, the two-dimensional diffusion theory calculation, the fission product treatment and depletion. The only calculation that is done in a separate program is the generation of macroscopic cross sections for control rods. Although LATTICE evolved from the United Nuclear Corporation code LOCALUX (Reference 1-1), which was a modi-fication of LASER (Reference 1-2), there is now very little similarity to either of these codes.
LATTICE requires as input easily obtained design information such as dimensions, densities, loadings, and number densities. Relatively few number densities are required, since the code calculates fuel and moderator number densitie:5. The output consists of data required by the TVA reactor simulator code, such as the assembly infinite multiplication factor, the two-group constants, the energy release per fission neutron, the xenon thermal microscopic absorption cross section, xenon and control reactivity worths, and the fraction of fissions occurring in uranium. Also printed are many other quantities of interest, such as isotopic concentrations, the fast and thermal spectra, reaction rates, and the assembly local fission power distribution. The methods which TVA has implemented to calculate these parameters are documented in this report, 1-1
CIIAPTER 1 REFERENCES
- 1. R. R. Calabro and J. M. Fiscella, "LOCALUX A Depletion Program for Multizoned Reactor Assemblies," UNC-5262, 1970.
- 2. C. G. Poncelet, " LASER - A Depletion Program for Lattice Calculations Based on MUFT and TIIERMOS," WCAP-6073, 1966.
O O
1-2
- 2. CRCSS SECTION GENERATION 2.1 CROSS SECTION LIBRARIES 2.1.1 Epithermal Dat..
HRG3 (Reference 2-1) is incorporated into the LATTICE code to calculate the epithermal neutron spectrum and the spectrum averaged epithermal cross sections. The epithermal energy range from 10 MeV to the thermal cutoff of 1.855 eV is divided into 62 energy groups, each of lethargy width 0.25. Following is a list of the materials present in the current HRG library.
Table 2-1 List of Materials in the HRG Library MATERIAL _ID SOURCE ENDF/B MATN0 U-235 12 III 1157 U-236 67 IV 1163 U-238 13 IV 1262 Pu-239 14 III 1159 Pu-240 77 III 1105 Pu-241 78 IV 1266 Pu-242 17 IV 1161 Gd-155 241 III 1252 Gd-157 243 III 1253 Xe-135 93 IV 1294 Sm-149 94 IV 1027 Fission Product 47 HRG ----
H-1 1 III 1148 0-16 4 III 1134 B-10 20 HRG ----
Zr-2 10 IV 1284 SS-304 11 HRG ----
The III or IV in the source column indicates that ENDF/B-III or ENDF/B-IV data is used. The data was processed using the SUPERTOG code (Reference 2-2), modified to write its output in the CAM-I format used by HRG3. The last column gives the ENDF/B material numbers that were 2-1
used. HRG in the source column signifs that TVA has not updated the cross section data for that particular material from that included in the HRG library obtained f rom the Argonne Code Center in 1972. The origin of this data is discussed in Reference 2-1.
There are resonance parameters in the HRG library for U-235, U-238, and Pu-239 for a maximum of 99 resonances each. Other resonances for these materials and all resonances for other materials have been treated assuming infinite dilution with Doppler broadening at 300 K. Resonance contributions are added to the non-resonance cross sections of the material.
Two other comments are in order at this point. First, the data for Zircaloy-2 is also used for Zircaloy-4. Second, cross sections for Gd-154 and Gd-156, while not included in either the epithermal or the thermal library, are programed into LAT'I1CE for use in the calculation of acti-vation rates for the gadolinium depletion chain. Since the macroscopic cross sections of these gadolinium isot< ces are small, they are not included in the spectrum calculations.
2.1.2 Thermal Data Thermal neutron spectra are calculated for 35 groups below 1.855 eV using a variation of THERMOS (Reference 2-3). Table 2-2 shows the mate-rials present in the current thermal cross section library. The material with ID 55555 is a lumped fission product with a 1/v 2200 m/sec cross section of 50 barns.
Scattering kernels were generated for hydrogen bound in water using S(a,p) data (ENDF/B material number 1002). A special version of the XLACS code (Reference 2-4), developed for TVA by ORNL, was used to process the data. For atoms heavier than hydrogen, an analytic free gas or diffusive model is used, as in FLANGE-II (Reference 2-5). The third column 2-2
Table 2-2 List of Materials in the THERMOS Library MATERIAL ID KERNELS UPDATED U-235 23592 0 Yes U-zJ6 23692 0 No U-238 23892 0 Yes Pu-239 23994 0 Yes Pu-240 24094 0 Yes Pu-241 24194 0 Yes Pu-242 24294 0 Yes Gd-155 20 0 Yes Gd-157 2713 0 Yes Xe-135 13554 0 Yes Sm-149 14962 0 Yes Fission Product 55555 0 No H 18 6 Yes 0 168 7 Yes B-10 115 6 No Zr 9140 6 No SS-304 30400 6 No Heavy Scatterer 92940 7 Yes indicates the number of scattering kernels present for that material.
For all materials with no kernels, the scattering kernel for the heavy scatterer, ID 92940, is use6. These kernels were generated using U-238 data and the TVA version of XLACS. A seven in the third column indicates that scattering kernels are present for temperatures of 80, 212, 390, 500, 550, 600, and 1160 degrees F; a six indicates that kernels are present for the first six temperatures listed.
Table 2-2 also indicates whether TVA has updated the cross-section data of the particular material from that in the origin- 1 LASER library to the ENDF/B-III data. We have made comparisons of the B-10, SS-304, and Zr-2 2200 m/sec absorption cross sections in our library to more authoritative values. For Zircaloy, our value is 3 percent higher than the ENDF/B-IV cross s-ction. Co.:garison values for the two other mate-rials were directly obtained from, or calculated from values given in, 2-3
EPRI-221 (Reference 2-6). The results were that for B-10, there was a 2b nut of 3837b dif ference, and for SS-304, our value was 3 percent higher than the iluoted value. These differences cause little error in the calculations and thus are acceptable.
O O
2-4
2.2 DANCOFF FACTORS An important aspect of LWR lattice physics methods is the calcula-tion of Dancoff factors which are used in evaluating the resonance inte-grals of the fuel. Sauer's method (Reference 2-7) for an infinite uniform array of pins is suitable for analyzing PWR's. However, a more detailed calculation is necessary for use in analyzing BWR's due to the presence of water gaps and water holes. Rod dependent Dancoff factors are calculated for BWR's using the DANCOFF JUNIOR routine of Gellings and Sauer (Reference 2-8).
2.2.1 Theory The Dancoff factor accounts for the reduced flux at a fuel rod surface caused by neutren absorption in the lattice surrounding the rod.
Since the blackness p of a region is the probability that a neutron entering the region is absorbed in that region, the Dancoff factor is equivalent to the blackness of the surrounding " moderator" region.
Referring to Figure 2-1, let n be a unit vector directed into the moderator, perpendicular to a fuel rod surface element dS, and let g be the direction of a chord in the moderator of length 2. Then (g n) is the cosine of the angle betweer a and 9 The following three assumptions are made:
a) the fuel rods absorb all incident neutrons, b) the neutron source distribution in the moderator has no angular or spatial dependence, and c) the resonances are sufficiently narrow that any scattering collision in the moderator re:. oves the neutron from the energy range of the resonance.
2-5
O 9
I
! i i
u s,u , l j
i
.' ~-
os .
j p
20 ly 4:~g 'l
/
l e <x .
I FIGURE 2-1: GEOMETRY FOR CALCULATION OF MODERATOR BLACKNESS O
2-6 4
Then the blackness can be written as (Reference 2-8):
(Q n)e" dgdS
=1 '
(0 n) > 0 (2-1)
(0 n)dOdS where f is a function of S and O. Define the chord distributieri function f(f) so that f(f)df is the probability that a chord is of length between E and (f+df):
(0 n)dQdS
//f 2 =
f(f)df = - - (0 n) > 0
-~
(2-2)
(0 n)dgdS where f is the length of the chord of interest (Reference 2-9).
Rewriting Equation (2-1) using Equation (2-2) yields the general expression for the blackness and the Dancoff factor:
C=0= f(f) (1-e~ ) d2 (2-3) 2.2.2 PWR Method Sauer (Reference 2-7) approximated the chord distribution function in the moderator as a simple exponential distribution, characterized by the mean chord length E, :.hifted by a distance L:
1 0 2<t f(f)df =
, exp -
\f-t/
dh
\f-t 2>t Substituting this distribution into Equation (2-3) gives the following equation for the Dancoff factor:
-TIE C=1- ---
(2-4) 1+(1-I)lf 2-7
where I=
fI+ -1 - 0.08 (for a square lattice),
V = the volume fraction of the moderator, V = the volume fraction of the fuel, and 7
I = the total macroscopic cross section of the moderator.
The moderator is defined as tne homogenization of all non-fuel regions.
2.2.3 BWR Method Referring to Figure 2-1, Gellings and Sauer (Reference 2-8) make the following substitutions into Equatic (2-2):
dS = R d$ dz, 0 $ $ $ 2n, -a< z<=
dQ = du sino dQ = du sino
~
(
0$o$n
/I [2 T -1<,<1 O
(0 n) = cos$ sinu = [1 T2 sino 2= $Nsino The first collision blackness becomes e n 2n +1 dz sin o du d$ exp -
g d(
_ y _
-m 0 0 -1 m n 2n +1 dz sin o du d$ d(
-= 0 0 -1 or 2n +1 E=1- Ki *
-2 3 0 -1 0
2-8
For a n/4 section, n/4 +1 (f)=1- fK 3 s) 4 (2-6) 0 -1 where Ki 3* "Y "" " "
exp(-x) (x + 6.399407x + 5.066719) 40.6366198x + 1.6211389 (x2 + 6.766116x + 5.066719)
I = the moderator total macroscopic cross section, s($,&) = the mean chord length in the moderator between the reference rod arnt the nearest rod in the lattice,
$ = the angle around the rod, 0 $ $ $ n/4 , and
& = the perpendicular distance from the line & to the point where the chord leaves the surface of the reference rod (Figure 2-2).
Thus the following information is required for the BhR lattice:
a) the distance from the reference rod to the water gaps for each n/4 direction, b) a library of the lattice configuration within a n/4 angle, and c) the rod type in each rod location.
Once the blackness is calculated for each n/4 direction, the rod Dancof f factor can be calculated as the sum of the blacknesses. The Dancoff calculation is repeated for each different rod type. For both 7x7 and 8x8 BbH fuel assemblies, we have assumed that the Dancoff factor for any rod more than two rod positions from a water gap is not affected by the gap. The water hole in the 8x8 design is also represented.
Af ter rod-by-rod Dancof f factors are calculated, the assembly Dancof f factor is found by averaging the rod-by-rod valurs.
2-9
O
, e d . 0 E
. e, U i
%_l) j *
.__.____ _ _____ _____ N a
l
/
%dL$3 l l
l L _ - _. __
_..-_____________.J
. e ..
4 .D - - - >
0
- l ()
2.3 RESONANCE CALCUATIONS The method used for the resonance calculations is identical to that used in the IIRG3 code. The derivation presented here is adapted from the itRG3 manual (Reference 2-10).
The contribution of a single resonance to the cross section for process x, absorption or fission, of the fine energy grotp g in an homogenized cell is:
E g-1 o*(E)$(E)dE E
o = g (2-7) x,g E ,3
$ cell (E)dF.
E g
where $(E) is the flux spectrum in the absorber region, and $ ell (E) is the flux spectrum of the homogenized cell. It is assumed that the resonances are isolated and therefore the total resonance contribution to the cross section in an energy group will be the sum of the contribu-tions of the individual resonances within that group.
In practice, it is much simpler to calculate the total number of x processes by a resonance:
m I =
x ox(E)$(E)dE (2-8) 0 and then find an approximate value for the fraction of x processes which occur in group g:
E 8-l f
x,8 1
I x( (' (2-9) x E g
2-11
Thus, Equation (2-7) can be rewritten as I .f
- *'8 a =
(2-10) x,g Q 8
where Q is the integral of the homogenized cell flux over group.g.
The calculation of the quantities in Equation (2-10) for a resolved resonance will now be considered followed by a description of the modifi-cations made in the method when it is applied to the unresolved region.
2.3.1 Resolved Resonances The resonance integral I x is evaluated using a modification of the method of Adler, Ilinman, and Nordheim (AHN method, Reference 2-11) that includes an intermediate resonance approximation for both the resonance absorber and an admixed moderator. The result is F
0 I =
E 0
- A '"A K where f is the level width for process x, E is0 the resonance energy, and n
0
= 2.6035x10 6 _F g (2-12) 0 The parameters F and F are the neutron scattering width and the total width of the resonance, respectively, and g is the statistical spin factor.
The Doppler broadening constant & is t=h y (2-13) where T is the temperature in eV and A is the atomic weight in neutron mass units. The parameters p d nd I d re determined by the expressions (K-A)o" + Ao F p P KA
=
o I. + AI, (2-14)
O a n 2-12
T Kg (K-A)o, + Aa p d (2- N shere K and A are the intermediate resonance pa rameters , N is the lumped nuclear density of the absorber, E is its equivalent mean chord length E = S(1-C) (2-16) where V and S are the volume and surface area of the absorber region, and C is the Dancoff factor. The parameter o is the off-resonance scattering cross section per absorber nuclide of all the atoms in the absorber region:
o 8' heavy
=a + (2-17)
N where o, is the scattering cross section per absorber nuclide of the admixed light atoms I
s, light o = (2-18) m N The functions J((,pd) nd L(Ig,(,pd) re defined as m
T( ,x) dx (2-19)
(b'OKA) " TY ,x[) + p 0
f'
- L(Igr>&,Dgg) ij,(g P 0
where x=f(E-E} 0 The Doppler function T is m
T((,x) = 5- dY exp
- h (x y)2 (2-22)
}47t 3 f2
-m 2-13
and P0 is the average escape probability from the absorber region. The J and L functions are evaluated through the use of tables.
The method by which Equation (2-11) was obtained is a. follows.
The AllN method uses the Narrow Resonance (NR) approximation for narrow resonances and the Narrow Resonance intinite flass (NRiff) approximat ion for wide resonances. For absorbing regions containing an admixed modera-Lor, the AllN equations f or these two approximations can be combineil into a single equation by letting A = 1 for NR or A = 0 for NRIN. This suggests that an intermediate resonance (IR) approximation, which assigns A a value between 0 and 1, would give improved results. The IR method used is that of Ishiguro (Reference 2-12), who also extended the method to include admixed moderators through the introduction of addit ional pa rameters.
Iloweve r , only a single admixed moderator (oxygen), characterized by the pa ramet e r K , is assumed to be present in the absorber (fuel) anil higher order refinements of Ishiguro are neglecteil. Equation (2-11) is the result when K is introduced.
The IR parameters are evaluated from expressions which depend on the resonance integral itselt. Iloweve r , it is assumed that this depend-ence is suf ficiently weak that a simplified version of the resonance integral can be used for this purpose. It is also assumed that unbroadened cross sections can be used to evaluate the IR parameters. Using the Wigner rational approximat ion to the escape probability P0 '
G p - "
0 0 + o (2-23) t e where a, = 1/(EN) is an effective cross section for escape, and expressing I'0 in terms f the quantities in Equation (2-11) yields 2-14
Dx g J(&,Ed) + L(Id ' t' ' Od) = p{ J((,ph)
(2-24) where
, gg (K-A)o,+ Ao +o e F (2-25)
0 n When unbroadened cross sections are used with Equation (2-24) in evaluating the resonance integral I and the IR conditions are applied, one finds that A and K are solutions of
~I tan X "A A=1- (2-26)
~I tan Y d
K= 1- y ,
(2-28)
KA (A + 1)'4 f I+ y
^ E 0 KA (2-29) m O KA (A + 1)2 7 [\ 1 +}/
where A, is the admixed moderator mass in units of the neutron mass and all other quantities are as previously defined. Equation (2-11) is then evaluated by use of these K and A. Note that the AHN assignments of each resonance to either the NR or NRIM approximation are used only in starting the iterative solution of the transcendental equation for A.
The Wigner rational approximation and unbroadened cross sections are also used in evaluating f in Equation (2-9). As f was defined, it is the ratio of the two integrals, each with the same inte-grand but with different limits. Reference 2-13 shows that the integrand inEquation(2-9)isthatofJ((,ph):
2-15
0x(E) Q(E) = P((x)
P({,x) +
Assume that the resonance cross sections are unbroadened, i.e., that T + 0 implies that ( + m [see Equation (2-13)]. Then the integral in Equation (2-22) can be evaluated:
I T(&,x) = (2-30) 2 1+x Substituting the above two equations into Equation (2-9) yields F *g-1 1
j 1+p'A+p7 K KA x
g f
_ (2-31) x,g pm 1
1+pyy+p}x Z J-m where 2
x =-
g F (E
g
-E) 0 (2-32)
These integrals can be evaluated analytically to give
'x,g
- 5 1 i
'^" 1 g-1
[ O ,A K
\ 1 + p})
}b '
-l
[ O ,A g\1+p]]
K
- f. 1 (2 '3)8 This expression gives a contribution from the resonance int gral ta every energy group. Since the contribution is very small for groups far from the resonance, the following cutoff procedure is used. Two energies equally distant from the resonance energy E 0 re < e emine such that 90 percent of the resonance integral is between them. The resonance integral is then allocated according to Equation (2-33) to the energy groups O
2-16
containing the cutoff energies and all intermediate groups. The contribu-tion of the resonance integral to groups outside the cutoff energies is included in the bounding groups.
Note that the value of f depends on the geometry and the scattering properties of the admixed moderator through {,butbecauseunbroadened cross sections are used, f is not dependent on the temperature.
The final quantity in Equation (2-10) remaining to be evaluated is the group flux $ . An equivalent two region model is assumed. For convenience the absorber region will be referred to as " fuel", the other region as " moderator," and the homogenized system as " cell." It is assumed that the cell flux is the volume weighted average of the two region fluxes:
V V
$cgig(E) =
y $ mod ( )
V D fuel (2%
cell cell In the resonance calculation, it has been assumed that $ mod ) I4' Using this assumption with the Wigner rational approximation, the fuel flux becomes 1
0 A (2-35)
@ fuel (E) = E T((,x)+p}
From the homogenization of the nuclide concentrations it follows that fuel fuel cell cell.
With these assumptions, E
$g & cell (
E g
2-17
E g-1 V V
$g = dE E
V en $m d(E) + \,g ; $ fuel(E) -
g N f n cel x,8 (2-36)
= 0.25 2E N O fuel / p, (1 + p, )
2.3.2 Unresolved Resonances In the unresolved range, Equation (2-10) still applies in principle but is modified in detail. For the unresolved resonances the following assumptions are made:
a) the narrow resonance approximation (i.e., K = 1, A = 1),
b) the Wigner rational approximation, c) spatial sel f-shielding of the resonance is neglected, d) the absorption and fission widths ox each resonance are equal to their averages in the unresolved range (i.e., f = <F > and f = <F >), and a a f f 0
e) the reduced neutron widths, I = f /E , ire distributed n n 0 according to t he l' ort er-Thomas d ist ribution "y/2 P(y) -
b where y=I /<f >
As a consequence of the first three assumptions, the resonance integral I""#"*
x can be calculated once the parameters of a resonance are known, using the approximation given by Equation (2-24) in Equation (2-11) with K and A = 1. can be set From Equation (2-36), &g is 0.25; and f x,g t_ a 1 in the fine group in which the resonance appears and 0 otherwise. Then, from Equatioe (2-10), the contribution of a typical unresolved resonance about E 0 t the fine group cross sections is 4I nres(EO )'
2-18
The unresolved resonance integral is found from Equation (2-11) using the above five assumptions:
a <V >
lunres(E) x = E J(E) (2-37) 0 where 3(E) is the J function averaged over the Porter-Thomas distribution:
m J[((y), g)(y)le ~Y!
Y '
3(E) =
J2ny 5 0 and "m,eff U
- 0 11 p e An adaptation of the GAM-I procedure is used to evaluate 3(E)
(Re ference 2-13) . A change of variables is made from y to the Doppler broadening constant:
( = B + Gy (2-38) where
<F>
B=
- ^
2 E 0
0 o= <r,(?)
so that the integral becomes
=
3(E) = -
J[(, 3(()) exp - (- * (2-39)
Eid [E B
Note that the integrand is singular at ( = B. The singularity can be removed by dividing the range of integration into two parts: from B to B + A, and from 8 + A to m, where A is a small number arbitrarily 2-19
set to 0.00095. Then the contribution of the first range of integration to the integral can be neglected. In addition, little error is introduced when the upper limit of iategration is changed from m to 1 since the integrand is very small for f, > 1. Therefore, the integral in Equation (2-39) is evaluated from B + A to 1. Simpson's rule is used to perform the integration.
The total contribution of unresolved resonances to the cross section in a fine group is E
g-1 unres unres @
o = 41 (2-40) x,g x g E
g where 5 is the average level spacing. Then the factor dE/5 in the above equation is the probability that a resonance occurs in the energy range E to E + dE. Substituting Equation (2-37) into Equation (2-40) yields E
unres <f> E~I o = 40 3(E)NE (2-41) x,g m,eff g E
g This integral is evaluated by a method that differs from the GAM-I method and is based on the fact that over the unresolved range, In 3 can be accurately approximated by a quadratic equation in the lethargy u.
To evaluate the coefficients in this quadratic, 3 must be calculated at only three values; two of these, u and u3, are chosen to be the 3
lethargies corresponding to the upper and lower energy 1imits, respec-Lively, of the unresolved range for the nuclide and the third, u i" 2'
chosen midway between them. Expressing the quadratic in terms of the mean lethargy u 0
1 Lpr n f fine group g WiMi is in W O
2-20
unresolved range and substituting into the integral in Equation (2-41) yielils "g-1 "g-1 J(u)du = J O "*P "("~"O}
+
("~"O "
Il 11 M E Assuming that the quadratic term in the exponent can be expanded, "g-1 3(u)du = 3 0 (" - "
8
+h (1-d + d2 )e' - (1 + d + d2 )e" (2-42) where a = 2u 0
b= 2 (In 3 - 2 In 32 + 1" 33) 3 c=f2 (u2 * "3)I" 1 - 2(u) +u)" 3 2 * ("I * "2)1" 3 d = 6 a/2 In 3 I" 0 * [u u 3)
"O~ "2)("0" "3)1" 1 - 2(uo - "1)("0- "3 2
+
(u0 - "l)("0" "2} " '3 3 = 3(u i )
j and 6 is the lethargy width of that portion of fine group g which is in the unresolved range.
2-21
2.4 EPITIIERMAL SPECTRUM CALCULATION The method used for computing the epithermal spectrum is the homogeneous B method used in the code llRG3. The derivation which is presented here is adapted from the HRG3 manual (Reference 2-14).
The steady-state Boltzmann transport equation can be written as 0 Vf(r,E,9) + q(r,E)f(r,E,0)
S(r,E)
=
1,( r ,E '4E ,g '49) f ( r , E ',g ')dE 'dg ' + { (2-43) where f(r,E,0) is the angular flux and the source is assumed to be isotropic.
To determine the energy dependence of the angular flux, the spatial dependence of Equation (2-43) must be removed. The spatial dependence appears explicitly only in the leakage term. Since leakage of epithermal neutrons is weakly dependent on the reactor's geometry, slab geometry can be used to approximate the first term of the transport equation.
Representing the spatial variable by z, Equation (2-43) becomes a
p g, f(z ,E ,0) + 1 ( z , E ) f (z ,E ,0)
= fl,(z ,E ' *E ,O '40) f(z ,E ',9)dE 'dg ' + hS(z,E) (2-44)
Now introduce the Legendre polynominal expansions for the flux and the scattering kernel:
2E+1 f(z,E,9) = 1 g - Y g(z,E)Pg(p) (2-45) f=0 I g(z ,E '+E ,g '4]) = I 4 Isf(z,E'4E)Pg(pg ) (2-46)
E=0 0
2-22
where p is the cosine of the angle between 9 and the z axis and pg = 9' g is the cosine of the scattering angle in laboratory coordi-nates. Substituting these expansions in Equation (2-44) yields 22+1 3 I g P g(p) p g $p (z,E) + 1 (z,E)$p(z,E) f=0 m
=1 Pg (p) Ist(z ,E ,' +E)$g(z ,E ')dE ' + S (2,E)620 (2-0) f=0 where 6 s .i Kr necker delta equal to one when f = 0 and egeal to zero 20 otherwise The Fourier transform is now taken:
2 P()
p -ikp 6 p(k,E) + 1 g(k,E) $g(k,E) f=0
= - - Pg (p) I gg( k ,E '-*E)$p(k ,E ')dE ' + 5(k ,E)6f0 (2-48) f=0 where w
h(k,E)=
g e * $g(z,E)dz
-m JHd m
5(k,E) = e
- S(z,E)dz
-m Note that k is equal to B, the square root of the geometric buckling.
Also, we have assumed that the cross sections are independent of z and that the flux approaches zero as z approaches plus or minus infinity.
Equation (2-48) is now multiplied by I (E)P (p) t 1 (E)-ikp 2-23
and integrated over all p, obtaining 2 -
1 (E)Q (k,E) = I (2f+1)Ag ,(o) 1,p (E ' +E)Q p (k ,E ')dE ' + 5(k ,E)620 f=0 m = 0,1, (2-49) where i t(E) o= g-and 1
Pg (p)P (p)
Ag ,(o) = gb o-U--- dp (2-50)
-1
=
o P (o) Q (o) j = min (f,m) n = max (1,m)
The nth order Legendre f metions of the seco..u 'ind Q (o) are defined Qo(o)= fin Q(o)=foin 3
-1
'I he B) equations are obtained fr<n Equation (2-49) by assuming that 1,g(E'4FT = 0 for f > 1. After some simplifications the following equations are obtaine<!:
kJ(E) + 1 (E)Q(E) = Is0,E '4E)$(E ')dE ' S(E) (2-51) kQ(E) + a(E)1 (E)J(E) : 1,3 (E '4E)J (E ')dE ' (2-52) where
-A 00f ")
a(E) = - g -
(2-53) 30 1-A 00 "
O 2-24
$(E)=6(k,E) (2-54) 0 J(E)=i@;(k,E) (2-55)
S(E) = 5(k,E) (2-56)
Integrating Equations (2-51) and (2-52) over each fine energy group results in the multigroup equations kJ +1 =1 I (2-57) a t,8 $g s'<g "U(E* '4E*) $* , + S
- 1 3 k$g + a Ig t,gJg = ,1 131(E 8 '4E8)J8. (2-58) g g where E
g-1
=
$ $(E)dE (2-59) 8 E g
g-1 J = J(E)dE (2-60)
E E
g g-1
=
S S(E)dE (2-61)
E E
g g-1 3
I t,g
=-
Q I t(E)$(E)dE (2-62)
E E
g Additional assumptions made were that the possible differences in I due to its different energy weightings in Equations (2-57) and (2-58) could be neglected, 1,3(E'4E) could be flux weighted rather than current weighted, and that a can be evaluated from the fine group equivalent. of Equation (2-53).
2-25
Equations (2-57) and (2-58) are applied to the homogenized assembly, hence the name " homogenized B method" The cross sections required to 3
selve the equations are averaged over a Cranberg fission spectrum with a 1/E tail for the iine groups in the last energy range in the resonance range, cont ribut ions to the line g r o u p r r o t: ', sections are calculated and ave ra y,vd as dese r i hed in the previous section. A lactor dependent on the assembly hydrogen-to-heavy metal ratio (ll/ft) is used to reduce the lumped fission product epithermal cross sections calculated by CINDER, Section 2.7, to account for resonance shiclling effects (Reference 2-15). A second tactor, obtained from Monte Carlo studies, accounts f or enhanced f ast fission of U-238 due to heterogeneity elfects in the lattice, and is also a f unct ion of II/fl. The flux spectrum $ can now be calculated by solving Equations (2-57) and (2-58) simultaneously beginning with t he highest energy gronp. Broad group parameters are obtained by averaging fine group values over the calculated spectrum for use in the CINDER fission product treatment and the assembly two-dimensional diffusion theory calculation.
O 2-26
2.5 THERMAL SPECTRUM CALCULATION The thermal spectrum is obtained by using a nodal integral transport method (Reference 2-16) to solve the one-dimensional steady-state multi-group integral t ransport equation. This method is used with isotropic scattering in a heterogeneous system in cylindrical geometry. The method consists of constructing a set of coupled equations for neutron currents in terms of loca! escape and transmission probabilities and first-collision sources.
2.5.1 Theory The steady-state, energy-dependent integral transport equation in the coordinate system shown in Figure 2-3 can be written as
$(r,E,g) = d r ' K( r '-+r ,E ,0) [Il( r ',E ,0) + S( r ',E ,0) ] (2-63) where
$(r,E,0) = the angular neutron flux, K(r'q ,E,9) = the first-flight transport kernel, S(r',E,0) = any external source, and ll(r",E,9) = the neutron emission density m
= dQ' dE ' I s( r ',E '+E ,g '1) $( r ',E ',g ')
4n 0 m
+
dQ 'f dE ' ol (g r ',E ',g ') X(r ',E '4E ,0 ' +0)
Ln 0
$(r',E',9') (2-64) 2-27
O I
Y a
~,
y AA1 di ,5 E
- f dx
>a f_ g' ~
v
> x 0
Y FIGURE 2-3', COORDINATES FOR THE INTEGRAL TRANSPORT EQUATION O
2-28
where 1 = the scattering kernel, x = the fission spectrum, and A = the eigenvalue of the system.
The nodal integral transport method t rans f orms the convent. i ona l integral transport equation into an equivalent neutron current equation.
The partial current at position r is defined as J(r,E,AD) = di2 g Q(r,E,g) (2-65)
Substituting Equation (2-63) for the llux and defining the first collison source as Q(r',E,9) = II(r',E,9) + S(r',E,9) (2-66) yields J ( r , E , oi)) = dr' y dil il K( r '+ r , E ,12)Q( r ', E ,12)
(2-67)
The integral over space is divided into two parts: an integral over the local region of interest and another over the neighboring regions. Define the local escape probability, local I AQ' O O (E I' ,Q)Q([ ,E,9)
P(r,E,a0) =
local I 4n'"2 9 Q(r ,E,9) the local average volume source, l '"l 'f 4 "'"' f' '}
@ E) -
E dg
, local J4n 2-29
the local transmission probability from neighbor n, dr' dQ 0 K(r'+r,E,0)Q(r',E,0)
T ' '
- n n
dr' y da 0 Q(r',E,9) and the part ial current of neighbor n:
1 (E) : dr' du O Q(r',E,0)
Then Equation (2-67) can be written as N
J(r,E,AQ) = P(r,E,aR)y(E)aV(r) + I Tn (r,E,ag)J n(E)
(2-68) 3 where AV(r) is ti.e volume of the local region.
Equation (2-68) is the neutron balance equation, expressed in terms of partial currents and a volume source, applied to the local region at r under consideration. The calculation of the point-to point transport kernel is now reduced to a calculation of a set of local escape and transmission probabilities, and the iterative solution of the fluxes is reduced to that of partial currents.
Once the converged set of partial currents is obtained, the neutron fluxes can be calculated using the law of neutron conservation applied to a given region of interest:
Collision rate = Volume source Net entering current I t(r,E)Q(r,E)aV(r) = Q(E)aV(r) + Jin(E' ') (
where the angular dependence has been removed by integration over all solid angles.
2-30
2.5.2 Method of Solution Equations (2-64), (2-68), and (2-69) form the theoretical basis of the nodal integral transport approximation. They are solved numerically in a modified TIIERMOS (Reference 2-17) with the method of cosine currents.
To be compatible with TliERMOS, the energy variable has been changed to velocity. The velocity space is divided into I intervals of width Av g and r space into K regions of volume AV . Define k
F ki Nf k'Vi)
S g = S(r_k'Vi}
li ki "(Ik '"i}
Q ki
- 9(Ek '"i)
= 11 . + S .
kt kt A uniform and isotropic volume source Q g is assumed for each region.
Also, it is assumed that t.he scattering kernel and fission spectrum are known so that. }{g can be calculated from Equation (2-64). The external source S g is known, thus Q g can be determined.
T.'ie currents J(r,E) can be defined in cylindrical geometry as shown in Figure 2-4. Then Equation (2-68) can be approximated as 5
Jkt. = P kt . Qkt . AV k + T kt? J+k-1,1. + Tkt? J'kt. (2-70)
~
J~ P*. (2-71) kt Q kt . AVk+TtJ.
. =
k-1,1 kt kt where for region k and velocity v g, P
g
= the fraction of the uniform volume source which escapes through the outer surface, P* = the fraction of the uniform volume source which escapes through the inner surface, 2-31
O
/
/ /
- /
/ , 'ni cioN j/
i j /
/ /
,/
/
/
l _
i d 'k .;
J' ,
,i /
/
/ / /
/
J*k w /. /
Ju
/ /
' i
/
,! i l
/
/
/
MESH POIN T k-l / ,e
/ /. '
/ / /
/ /
/
/
/
,/
MESH POINT k FIGURE 2-4: DEFINITION OF CURRENTS O
T = the fraction of the cosine-shaped current on the inner surface which passes ihrough the region and escapes through the outer surface,
'l[} - t ue f raction of the cosine-shaped current on the outer surface which passes through the region and escapes through the outer surface, and T = the f raction of a cesine-shaped current on the outer surface which passes through the region and escapes through the inner surface.
The above escape and transmission probabilities are interrelated:
T =u ( -72) k T?=1-p.P.
kt i
kt kt (2-73) 5 T"ki?=1 -u k 1 + pki . (Pki . -Pki) (2-74) with a = R /R k-1 (2-75)
O ki ki k k (2-76) where l
g = the total cross section of region k at velocity v j, R k = the outer radius of region k, and Sk = the outer surface area of region k.
Thus, only the escape probabilities P g and P need to be calculated.
As a result of the assumption of cosine currents, the escape probabilities are available from a correlation dependent on o and the region optical thickness.
2-33
The nodal flux is approximated from Equation (2-69) as:
F ki I Vp 9ki # k * ( k-1,i ~ k-1,1) ~ f ki ~ ki} (~ )
Since the volume source is a function of the flux, the above equation must be solved iteratively. Using an initial flux estimate, the volume source is calculated. The inward currents for each region are determined from Equation (2-71) beginning with the outermost region. Then the outward currents for each region are found using Equation (2-70) beginning from the innermost region. A zero net current boundary condition is satisfied by setting J+g = J~g at the outer suriace. The flux spectrum is determined using Fquation (2-77) and the process repeated until convergence. The rate of convergence is improved using over-relaxation:
F iI = Fd . + 6(F l ki ki ki
-Fl)ki (2-78) where j is the iteration counter and 6 is the over-relaxation parameter.
2.5.3 Geometry In the spectrum calculations, the fuel pins are grouped according to their spectral environment. For a typical 8x8 BWR assambly, four fueled spectrum regions are used with a fifth spectrum regiot containing the channel and water gap surrounding the assembly (Figure '-5). The flow area inside the chacriel not included in the unit cell areas defined by the lattice pitch is adoci to the moderator area associated with the outer ring of pins. Any wate. ._ a re included in spectrum region one.
Each fuel or poison pin and its surr3unding clad and moderator is represented by a cylindrical unit cell of infinite length with the moderator radius chosen to conserve the total area of the cell (Figure 2-6).
An average fuel cell is determined for each fueled spectrum region, 2- M
_____3 l 5 i r@@@@@@@@%
l @@@@@@@@
l @@DDD@DD l @@@@DODD l @@@@@@@@
l @@@@'@@@
C fD D @\DD @ @,
\s \ \
-FUE L PIN ZlRCALOY CHANNEL
' CONTROL BLADE WATER ROD BURNABLE POISON PIN FIGURE 2-5: BWR LATTICE CONFIGURATION 2-15
O O
%n =
MODERATOR FIGURE 2-6: CYLINDRIClZED FUEL CELL O
Ro i
J Rm POISON ODERATOR c // HOMOGENIZED FUEL, CLAD &
MODE R ATOR C , WATER RING FIGURE 2-7; SUPER CELL MODEL O
2-36
ignoring any poison cells, for which the spectrum will be calculated.
This cell is typically divided into five rings in the fuel, a single ring in the clad, and six rings in the maderator. For poison cells, more spatial rings are included in the poison region, and a coupling region of homogenized fuel, clad, and moderator is added around the outside of the cell to create a super cell model of the poison cell and the surrounding lattice (Figure 2-7). The coupling region is needed to obtain the correct thermal spectrum in the poison pin; omission of the coupling region would represent an infinite array of poison cel. The composition and cross sections of the coupling region are obtofn' from the fuel cell calculation for the spectrum region in which the poison , resides. The coupling region radius corresponds to five times the lattice cell pitch. Outside the coupling region, a water ring is added with the thickness chosen to match the average neutron velocity in the coupling region to that in the corresponding spectrum region obtained from the assembly spectrum calculation described below.
A space-dependent spectrum is determined using a one-dimensional cylindrical representation of the assembly to account for the interaction of spectrum regions. The assembly is divided into concentric rings, each corresponding to a spectrum region (Figure 2-8). The inner rings are the fueled spectrum regions ardered by spectrum region number with region one in the center. The outer ring corresponds to the homogenized channel and water gap and has a sign;ficant sof tening effect on the other spectrum regions. The ring thicknesces are chosen to conserve the volume fractions from the actual assembly configuration. Each ring is subdivided into several n.esh intervals with the number proportional to its thickness.
2-37
e R1 R3
\ j e FIGURE 2-8; CYLINDRIClZED ASSEMBLY (8x8 BWR) e 2-38
2.5.4 Therma _1 Cross Secti_ons Space dependent thermal spectrum calculations are made for the average fuel cell in each spectrum region. Effective microscopic cross sections in each energy group are determined by spatially aver-aging over the cell flux for use in the corrtsponding homogenized assembly region. For spect rum regions containing poison pins, the fuel cell cross sections ar- used in the coupling region of the super cell. The spectrum from the super cell calculation is used to average the cross sections in the super cell for the corresponding assembly region. Reaction rates from the fuel and super cells are conserved by modifying the assembly cross sections as follows.
For nuc1ide i, the celi reaction rate is X (E) = o.(E) N.(r)$(r,E)dr C 1 1 (2-79) where the integral is over the cell volume, c.(E) = the THERMOS library cross section for nuclide i, 1
N. (r) = the concentration of nuclide i, and 1
$(r,E) = the cell flux spectrum.
Expressing the nuclide concentration and the flux as volume averaged quantities, Equation (2-79) becomes X'(E) e = o.(E)N'$c(E)V i c d*(E) cc (2-80) where Nf= N (r)dr g
'c
$ (E) = y-1 $(r,E)dr C
c V = the volume of the lattice cell, and df(E)=thecelldisadvantage factor for nuclide i 2-39
l
N g (r)$(r,E)dr O d (E) - I Ne $((E)V c Since nuclide reaction rates in the cell are conserved in the assembly region, X'(E)
A
= X'(E) e or o'g(E)N'g$ 4 (E)V 3 = 0;(E)d'(E)N'$ (E)V c (2-81) where A denotes the assembly region. Thus the effective cross section for an assembly region is i
N? $(E)V '
J*(E)
= 0.(E)d*(E)
(2-82)
N'gQ A(E)VA Assuming that the flux integral is the same for the cell and the correspond-ing assembly region, Equation (2-82) becomes
. . N' oA *(E) = 0.(E)d*(E) i c 4 (2-83)
N*
A or, in terms of non-averaged lattice cell parameters:
N.(r)$(r,E)dr 1
c oA*(E) = o.(E) i (2-84)
N,'g $(r,E)dr c
The effective thermal cross sections determined from Equation (2-84) are used to calculate the thermal neutron spectrum in the assembly. Then the broad thermal group cross sections for the two-dimensional diffusion theory calculation can be found. The unit cells are modeled as homogenized regions in both the assembly spectrum calculation and optionally in the O
2-40
diffusion theory calculation; hence, the broad group cross sections for each nuclide in the cell are the fine group cross sections for the corresponding assembly region averaged over the flux spectrum:
E th .
o#g'(E)$,g(E)dE th 0 (2-85) i E th QA (E)dE where E g is the thermal group cutoff energy.
Since diffusion theory significantly underestimates the flux depression in a highly absorbing rod, the poison cells are not homogenized in the two-dimensional dif fusion theory calculation of the assembly. Each poison cell is modeled as a square poison region surrounded by its homogenized clad and moderator. This requires a procedure similar to that above to adjust the poison region cross sections to agree with the transport theory calculation.
No correction is made to the clad-moderator cross sections. Therefore, a fixed source, one group, two-dimensional diffusion theory calculation is performed for a square representation of the super cell (Figure 2-9). Then the diffusion theory absorption cross section of the poison region is adjusted to conserve the fraction of absorptions in the poison region determined from the transport theory calculation for the super cell:
Vi 6 Vi p a,p p = p a,p 5 n (2-86)
Vi 6 c a,t c Vic a,c hc where Vp = the volume of region R, l
a,R = the average macroscopic thermal absorption cross section of region R 2-41
O WATER RING I
! N N HOMOCENIZED
! \ FUEL, CL AD AND MODERATOR I ~
, COUPLING RECl0N i
t /
/
jI, l /
CLAD AND -
WMAMf MODERATOR f l POISON I%)
1-- ,L I
l ALL BOUNDARIES HAVE MIRROR SYMMETRY t
FIGURE 2-9: TWO-DIMENSIONAL MESH FOR A SUPER CELL O
2-42
E h
I N.(r)o.(r,E)drdE y _ i 0 R a,R E th dE dr 0 R
~)g = the average thermal flux of region R F
l Q(r,E)drdE 0 R E
h dE dr 0 R The subscript p refers to the poison region, c refers to the super cell, and the superscripts D and T refer to diffusion theory and trans-port theory results, respect ively.
Solving for the dif fusion theory absorption cross section in the poison region yields i",p a = TD i a,p (2-87) where
-D -D -T 1 $ $
T -
(2-88)
D I
a,c $c 0 p This correction is applied to all thermal cross sections for the poison region:
ogh = T a th 1,p D i Q-89) where o is the cross section of the nuclide in the homogenized cell, given by Equation (2-85). Eiguations (2-87) and (2-88) must he solved iteratively because of the interaction of the dif fusion theory cross sections and flux, but few iterations a re required for convergence.
2-43
2.6 EFFECTIVE CONTROL ROD PARAMETERS The effective few group parameters for use in the control rod regions of the two-dimensional assembly calculations are determined using the c-de RODWORTil (Reference 2-18), modi f ied by the T(nnessee Valley Authority. Different methoils are used to determine the effective few group parameters depending upon the geometry of the control rods; for cruciform rods consisting of a linear array of cylindrical poison pins, the method devised by flichelini (Reference 2-19) is used, while the method of Sha (Reference 2-20) is used for rod cluster control (RCC) rods.
2.6.1 Control Blades
?!ichelini's mrchod (Reference 2-19) combines P3 blackness theory with geometrical considerations of tbc lifferent paths that neutrons tan take through a control blade filled with cylindrical poison pins.
There are two possible paths of neutron penetration across the blade:
through the slit due to the clad of the poison pins and through the poison itself (Figure 2-10). This differentiation allows the poison to be first considered black in order to calculate the transmission through the clad, then the additional transmission through the grey poison region can be evaluated.
First, it must be determined whether the poison i., grey or black to neutrons in each energy group. To do this, the optical thickness 21,(R a) of the poison region is estimated using fine group cross sections from a built-in library. If the optical thickness is greater than 7, the poison is considered black, and the t ransmission of the neutrons through the blade must be due entirely to the slit represented O
2-44
0 X V,
I W 3 I
I I
i
-p-la ',4--
I I I
I I I i 0
1
~
l i A l
___4____ - __ __
_e____ _ _ .____L
/
FIGURE 2-10: GEOMETRY OF POISON PINS IN A CONTROL BLADE 2-45
by the clad. Due to the small thickness of the clad, the probability that a neutron crossing the blade will undergo a collision in the clad is very small. Thus it is assumed that the fraction of the incoming neutrons which successfully cross the wall is, at any abscissa x, related to the angle (y} + Y2 ) between the two trajectories tangential to the cylindrical poison regions (Figure 2-10). This angle is given by sin y (x) = ( **} ( +* + ' ~ (
(2-90) 3 2 (R+x) + R sin y,,(x) = (R-x)[(R-x)2 + 2aR?
- a ]b - R(R-a)
' ', ~ 9 I }
' 2 (R-x)" + R where a is the thickness and R is the outer radius of the clad. Assuming a diffusive outside medium, the angular distributioi of the incoming neutrons is T.1 (u) = -WL ($.1 + 3J. coso)
I (2-92) where Jj are the incoming net currents, $g are the fluxes on the i = A,B opposite sides of the blade, and u is the angle between the neutron trajectory and the norm.il to the blade surface If the 0-axis is thosen parallel to the axis of the poison pin, the angular distribution becomes 1
P (0,$) = g (&f + 3Ji sine cos&),
g (2-93) where the Q-axis is normal to the interfaces between the pins.
Consider now the net currents entering the two sides at the point of abscissa x (Figure 2-10). The balance between the incoming and outgoing neutrons is, in spherical coordinates, O
2-46
n y (x) 3 2
J (x) 3 =f$(x)+fJ(x)-
3 3 Ig sin 0 do cos$ d$
(2-94) n y (x) 3
- [J g 0
3 sin 0 e
-y2 (*)
2 cos 4 do n y (x) 3 J(x)=j$g(x)+y 3
J (*
B
~
A O -y2(*)
(2-95) n y3 (x) 3 2 3
A sin 0 do cos $ d$
0 -y2 (*)
where $g and 3g (i = A,B) are mean values along x.
Defining two transmission probabilities as V3 (x) p (x) E f g cos$dQ=f(siny(x)+siny2(x)] (2-96)
-y2 (*)
v (x) 2 2
'Y (x) + y (x) l 2 sin 2y3 (x) + sin 2y2(*)
p2 (x) c s Q dQ = - ~+
2
-y2 (*) 4 (2-97) then Equations (2-94) and (2-95) become J(x)=fQ(x)+
3 A JA(*) ~ BP i(X) - 3 P B 2(* (
JB(*) B*
- B
~
P AI ~ i A P2(* (
To clear Equations (2-98) and (2-99) of the spatial dependence, the mean values of the currents, fluxes, and the transmission probabilities are introduced:
2-47
h: i P1 (x)dx , 6=
2 P2 (x)dx (2-100) 0 0 where Ax is the distance at which p3(x) and p2(x) v nish. Now combining Equations (2-98) and (2-99) yields the P gblackness theory parameters a and p:
A+ B 1 I- I a= _
= -
6A+ B I*
2 (2-101)
A-B I*
- p. =
_ 1 I 6
A' B - 2 The expresrions for p 3and U have 2 been approximately evaluated as P
10 1 Pi R (2-102)
- P20I2 p = - ~ ~
R The first factor in these formulas is the transmission probability evaluated from Equations (2-96) and (2-97) at the point x = 0 under the assumption that siny j (0) ~ y g(0), i.e., a/R << 1:
P 10 I ~ 2(R+a) (2-103)
P.' i R~ ~ 2(R+a) (2-104)
The second factor is an estimated shape factor of the function p g(xt in the interval 0 < x 5 Ax, 0.167
- f. = 12 + 1+w.a/R , i = 1,2 i
(2-105) 1 At the point x = 0 these functions have a maximum and then gradually decrease and reach zero at Ax. Michelini reports that the values for w) and w are 2.33 and 3.70, respectively, for a/R < 0.3.
2-48
The last factor, ax=Hla(2R-alf2 R-a (2-106) is, as defined for Equation (2-100), the distance at which p (r.) 3
= p (x) 2
=0 (the distance at which the slit is no longer seen by the neutrons), and is determined from Equations (2-90) and (2-91) by setting y3 (Ax) = -y 2(O*)*
In practice, the validity of the assumption used in the evaluation of the first and second factors of Equation (2-102) is confirmed since the clad thickness is not large with respect to its outer radius. Thus, the trans-mission probabilities for the case of black poison in the pins are:
^ " ( - ) 3+
- p =K i = 1,2 (2-107) i i 1 2(R+a) R(R-a) 2 1+w.a/R 1
where K, = 1 K2 = 4/n w g = 2.33 w = 3.70 2
If the poison is grey, i.e , the optical thickness is less than 7, then a fruction of the neutrons which impinge upon the poison region is transmitted through the pin. If one sets p Md = (P; + P2)/2, and normalize to one entering neutron, then pblack represeits the neutrons passing through the slits and the (1 p black) remaining neutr ns must have trajectories which cross the poison material. If p, represents the attenuation suffered by the neutrons crossing the poison, then the total transmission nrobabili;y across a blade with grey poison in the pins is P=Pblack + (1 Pblack) P o (2-108)
Since neutrons may pass through more than one poison region, the calculation of p is a formidable task. Therefore, p is approximately 9
evaluated based on the fact that the probability p,(x), i.e., the 2-49
attenuation suffered at any point, does not vary significantly along x.
To explain this fact consider the neutrons crossing the poison material as being of one of two groups: those which cross only one poison region and those which cross more than one.
The first group of neutrons cross only one cylindrical poison region with the cylinder seen at different angles by the neutrons. Consider the mean attenuation suffered by the current of neutrons which crosses a poison region starting isotropically from a point on the x axis at the distance [R + (R+x) ]b from the pin axis (Figure 2-10):
n c 2
sin 0 de cos(y+Q)exp'-1,c(0,Q)]dQ p(x) - n (2-109) e sin 0 de cos(y+Q) dQ 0 -c wl.e re
. -1 R-a .-1 R+x c = sin ' y = sin 2
[R +(R+x)2)\ [R +(R+x) ]
and h
2 c(0,Q) = - [R2 ,(g, )2] sin Q 2 (2-110) siaO
{(R-a)2 is the segment of trajectory withu. the poison.
Expanding the exponential up to the second derivat.ive and performing the integrations results in the following expression which is valid for low absorhing materials (1 d < 0.3):
a -
p = l-1 d + 2/3 (1 d) a a (2-111) where d = 2(R-a) is the diameter of the poison region. Note that Equation (2-111) does not show any dependence on x. Thus, the mean attenuation p remains unchanged within the poison, neglecting the high-order derivative effects not appearing in Equation (2-111).
2-50
The second group of neutrons suffers an attenuation p which produces a small decrease in the total attenuation p . Since the sum of the 9
path lengths within the various poison regions is about the same for any given trajectory, p is practically the same at any point in the poison.
Numerical estimations of p and related averages on the populations of the two neutron groups gives consistently for a wide range of a/R ratios p9 2 p [0.84 + 0.16 exp(-0.61,d)] (2-112) which is valid under the same restrictions as Equation (2-111).
Since the restriction on p is due only to the truncation of the power expansion, one suspects that a more general law exists whose power expansion approaches Equation (?-111). The following expression has been found valid in general:
p9 2 exp(-I d)[0.84 + 0.16 exp(-0.61 d)] (2-113)
It should be noted that owing to the uncertainties of the spectral averages and to the approximations needed in the derivation of p , the 9
differentiation of the transmission probability into two components, p 3
and p 2, I ses its significance. llence, the standard blackness parameters reduce to a 1. I ~ P 21+p (2-114) 1 1 +p 0_51
~
p where p is given by Equation (2-108).
The group average blackness parameters are obtained by weighting a(E) and p(E) by the flux spectrum. For the single thermal broad group, the flux spectrum is chosen such that the average parameters' yield the same absorption rate obtained from the multigroup calculation. The spectrum required is the asymptotic flux in the medium surrounding the 2-51
control rod (Reference 2-21). This is the same as the flux spectrum in the region with the control rod withdrawn, which is calculated by LATTICE for input to RODWORTH. Thus, for the thermal group, th a(E)$,(E)dE 0
"th E g
Q,(E)dE O
(2-115)
E g
E(E)Q,(E)dE il g =b -- -
Q,(E)dE O
where E g = 1.855 eV is the cutoff energy of the thermal group.
In the epithermal energy ranges, the flux spectrum used for averaging is that at the control rod surface (Reference 2-21):
F'g+ 1 a(E) Q(E,R)dE E
E Ir =
g E ,3 Q(E,R)dE E
8 (2-116) g+1 (E) $(E,R)dE E
g p =
g Eg ,3 Q(E,R)dE E
g where E is the lower energy limit of broad group g, and R denotes the control rod surface (Figure 2-10). Diffusion theory is used to determine 9
2-52
the flux distribution in the region external to the control rod, assuming a uniform spatial source distribution (Reference 2-21). The flux dis-tribution is evaluated at the control rod surface and substituted into Equation (2-116) to determine the average blackness parameters.
The diffusion equation for any energy group in the region around the control rod is D(E)V 2 Q(E,x) - l (E)$(E,x) + S(E) = 0 (2-117) a with the boundary conditions Vc = 0 at x = 11
=a at x=R where 11 is the distance to the external region boundary. Note that for nonthermal energy groups, I includes bot.h absorption and removal from the group due to scattering.
For a control blade, Equation (2-117) yields a A1 (*) ~ A 1(R)~
+ liK A., (R )
$(E,x)={q - - U3-(2-118) 3 A_.
1 DK A.,(R) where d = 1/D a A,(x) = sinh (Kx) - roth(Kil) cosh (Kx)
A (x) = cosh (Kx) - roth(Kil) sinh (Kx)
When evaluated at the control rod surface, Equation (2-118) becomes
$(E,R)=f- g-- (2-119) j - d su) .
2-53
For practical problems, KR is small and Kll is large, so A (R) t - coth (Kll) E - ]
3 A2 (R) 2 1 Also assuming that I z1 ,
DK 1 = J31a/1tr ? [3-Substituting yields I) I Q(E,R) - (2-120) a( ) ,
1 + /5 o(E)
The factor in brackets accounts for the flux depression at the control rod surface relative to the asymptotic flux S/I,. Equation (2-120) can be rewritten as Q(E,R) = Q ,(E) -
(2-121) 1 + /5 a(E) where Q,is the flux spectrum calculated by LATTICE and input to RODWORTil.
Therefore, for the epithermal groups, Eg+1 a(E)Q,(E)dE
=
g I + /5 a(E) g,g)
E g+1 Q ,(E)dE b
g 1 + d5 a(E)
A similar derivation using p instead of a as the boundary condition at the rod surface for the dif fusion equation yields EE+I (E)Q,(E)dE g = __3 1+J5_fl(E)- 2-R3) g EE *I Q,(E)dE g 1 + d5 (E) 2-54
The macroscopic absorption cross section and the diffusion coefficient of the blade for each energy group are calculated for zero and one internal mesh point by the equations:
a) Zero internal mesh points a
=2U (2-124)
D= (h-5) where h, the mesh thickness, is the diameter of the poison pins.
b) One internal mesh point I,=f6-[h(h-5)
(2-125)
D =h /h(h - 5) where h is the poison pin radius.
2.6.2 C R_C Rods Sha's method for determining the effective few group parameters for RCC rods (Reference 2-20) is based on applying the diffusion equation in the region exterior to the control rod subject to transport theory boundary conditions at the surface of the rod.
For cylindrical absorbing elements thac are uniformly distributed throughout the core (or follow a definite pattern), the problem may be treated by considering a control red cell as shown in Figure 2-11. The control rod with radius 20is located at the center of the cell ard the lattice area associated with the control rod is replated by the concentric annulus with radius R g. In the following derivation the control rod is assumed to be infinitely long and thus the problem under consideration is reduced to one dimension.
2-55
O FIGURE 2-11: RCC ROD CELL 9
2-56
The neutron balance through the control rod surface is J.in(E) - Jout(E) = p(E) J.in(E) (2-126) where J g(E) = the neutron current per unit surface area in the energy interval E to (E+dE) entering the surface of the control rod, Jg (E) = the neutron current per unit surface area in the energy interval E to (E+dE) leaving the surface of the control rod, and (E) = the incident neutron capture probability for neutrons with energy E (blackness).
Assuming diffusion theory to be valid in the region external to the rod, the net current at the surface of the control rod is D(E)V$(E,R 0 ) = J.in(E) - Jout(E) = (E) J.in(E)
(2-127)
The entering current Jg(E) can be expressed as J (E)=f$(E,R)*i O ( ' 0) (~
Substituting Equation (2-128) into Equaticn (2-127) yields
$(E,R D(E)V$(E,R '
0 4
-0 + D(E) 2 ' 0)
Therefore, the ratio of neutron current to flux at the surface of the rod is D(S)V$(E,RO ) p(E) a(E) 2 - (E,R ) 2-p(d O
In crder to calculate a(E), the blackness of the control rod p(E) must be evaluated. Blackness is a difficult quantity to calculate, and is a function of the control rod material, the size of the rod, and the neutron current distributica entering the surface cf the rod.
2-57
However, has been calculated extensively (References 2-22 and 2-23) and is shown in Figure 2-12. The energy dependence of the blackness is implicit in the control rod cross sections used in reading the plc,t.
For the thermal energy range, represented by one broad group, Maxwellian averaged cross sections are used to determine p since studies indicate that the use of hardened constants leads to poorer agreement between theoretical and experimental results. Therefore, only an average value for the blackness is obtained, and the thermal group current-to-flux ratio at the rod surface is calculated from Equation (2-129):
O th
~
ag = (2-130) 2(2-E,dt)
It should be pointed out that Equation (2-129) was derived based on cylindrical geometry and diffusion theory. lloweve r , these results will be used in the two-dimensional calculation of the assembly with square geometry and transport corrected cross sections. To account for these differences, Equation (2-130) is modified as follows:
- O th og = C C (2-131)
T 3 2(2-p,g)
The transport correction factor C accounts for the er1cet of predicting T
the extrapolation distance with diffusion theery instead of transport theory. Since a is inversely proportional to the extrapolation distance, T
(
tr)!0 (2-132) where 6 is the transport theory extrapolation distance and A is the trans-port mean f ree path of the mediu.a surrounding the control rod. Figure 2-13 is a plot of 6/A as a function of R 0! tr' 2-58
1s/1't 1.0 i
!~ '~~ _ o b ~ ~ _ o.f
/ '
2 '__. ; _ _' ~ ~~ o ,
~
i ' :
[' ~. - ~~~ ~ ! O]
j
~~~ ~ ,'_ %_ .a 0.5 l
0.8 ~ - ' .
, ,_ i
- -to.s '
~~ 40.7 0.6 ) w- >
/ //
/
,i - %,
D l /f i l j ///,,;/ ,
s .4 , ~
l +~~
$ I
- /
,///f ~ , ; ~~-
_:[9 i
/, ', , >
i
- 0. 2 ~
IIV// /m,vl
,;=**
___ +~~. ~
t i
' ~_._^' L ' ' 4 I
i 'r/ /:' ' ,
/l / >
0.- i 1
O 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 Sa Ro FIGURE 2 32; gLACgyggg LINDRICAL ABSORBERS
O 1.4 p
1.3 E
E 1.2 e
LL 4 /A tr e
1.1 1'I 2
o
'e D 1.0 1.0
.c o
E 0.9385 E 0.9 6 0.9 C -
7 -
.h N
3 '
S 0g % 0.8 N
la h 0.7104 0.7 07 C
m 0.6 0.6 0.5 0.5 0 3 3 5 g
- A tr. Radius og 9 in L* nits of Mean Free p,,9 FIGung g_,3 swart CORRECTION FACTOR FOR CYLINDRICAL AssoRBEas
- -60
The factor C is to correct for the change in leakage from cylindrical 3
to square geometry. It is assumed that the leakage from the rod should be conserved and the neutron current would remain constant after the trans-f ormat ion o f geometry. Therefore, C is de rado of de actual rod sudace 3
area to the surface area of a square rod 2nR C =
3 4t 0.8862 where t is the length of the side of the equivalent square rod. Thus, for the thermal energy range, th E
g = 0.4431 _
C T (2-133) 2-S th In the epithermal energy range.. a(E) is calculated from Equatica (2-129) where p(E) is obtained from Figure 2-12 using control rod cross sections from a built-in multigroup library. The group average a is determined by weighting by the flux spectrum at the rod surface. The required flux spectrum is derived in the same way as ,for the control blade.
$(E,R * $m }
a A1 (R O }
O
~ Dil A2 (RO.
)
where, for a cylindrical rod:
Kg (KRI )
O )
- I (KR ) 0(* 0}
- 0(* 0)
Aj (R g g K3 (KRg )
A2(R0 =
I (KR ) I 1(KRO ) - 1(KRO )
g and K = 1 /D is evaluated for the region surrounding the control rod.
For practical cases, KR0 is small and KR 3is large, so 2-61
A3(RO ) K0(KRO )
A2(RO ) ~ 1(KR0 )
Also assuming that I, 2 I yields Q(E,RO ) * @m(E) (2-135) g(KR]
0
- 1 + [3 a(E) K 1(KR0) -
Note that, unlike the control blade, the flux spectrum at the RCC control rod surface is dependent on the external region properties. Thus, E8 +I a(E)Q,(E)dE K0(KR 0 I*
E "- ) K (KRO )
3 a = E,. (2-136) g c.
8 Q,(E)dE K0 (* 0)
E 1 + [i a(E) K 3(KRO )
Now treating the control rod as a diffusive medium instead of an internal boundary, the flux and current distributions in the rod as derived from diffusion theory are:
$(E,r) = $,(E) I O(Kr) (2-137)
J(E,r) = Q,(E) KD 1 (Kr) (2-138)
By definition, the blackness parameter is l y(KRO) 10 (KRO )
I
= [1 D I)(D R 0/
~
R 0 O e
2-62
To determine the control rod absorption cross sections, its diffusion coefficients are assumed to be unity (the reactivity worth of the control rod is inse:nsitive to the diffusion coefficients). Then Equation ('-139) can be solved for I, of each group using the appropriate 5: Equation (2-133) for the thermal group and Equation (2-136) for the epithermal groups. Figure 2-14 is a plot of R as a funcdon of R 0 a 0 determined by this method (Reference 2-18).
2-63
O 1.8 16 12 --
s a"
~.
08 04 #
0 0 02 0.4 06 08 R&
FIGURE 2-14: R,' Ea VERSUS R,& FOR RCC RODS 2-64
2.7 FISSION PRODUCT TREATMENT A modified form of the CINDER point depletion routine is used to calculate the thermal absorption cross section for the lumped fission product. The calculation includes 69 chains consisting of .,3 nuclides.
Since CINDER is well documented (Reference 2-24), only the modifications to the routine will be discussed.
Average epithermal and thermal fluxes plus fission and absorption cross sections of the fuel nuclides are required to calculate the con-centrations of nuclides in the fission product chains. The fluxes are assembly averaged values obtained from the two-dimensional few-group diffusion theory calculation. The epithermal microscopic cross sections are determined by averaging assembly trans p rt theory multigroup values over the broad energy groups. The thermal microscopic cross sections are assembly averaged values determined from Equation (2-85).
The CINDER library contains the thermal cross sections of the fission product chain nuclides at a reference velocity of 2200 m/sec. Ilowever, the thermal cross sections used in LATTICE are average values, not 2200 m/sec values. Therefore, the CINDER library cross sections n.ust be modified to appropriately averaged cross sections for the CINDER calculation to be compatible with LATTICE.
The U-235 thermal absorption cross section o,23s(v ) calculated by LATTICE has the case-dependent rpectrum and fuel flux depression effects incorporated into it. Thus, the ratio of o,23s(v ) to the 2200 m/sec value o,23s(y0) represents the cross section reduction due to these effects:
o 235(v )
F=" (2-140) o,23s(y0}
2-65
This factor can not be applied directly to correct the CINDER cross sections since U-235 is a non-1/v absorber while CINDER assumes 1/v behavior of the nuclides in the fission product chains it represents.
To calculste a correction for the non-1/v behavior of the U-235 cross sections, a LATTICE depletion case was run with a trace amount of B-10, a 1/v absorber, in the fuel pins. The ratio of the effective B-10 thermal absorption cross section to the 2200 m/sec value edited at each depletion step is the total correction factor to be applied to the CINDER cross sections. Comparison of this factor and Equation (2-140) for the same case indicates that Equation (2-140) gives results 10 percent too low.
Therefore, the correction factor is adjusted to 023s(y )
^
F = 1.10 (2-141) 023s(y0) a The above comparison was based on 2200 m/sec cross sections of 694 b for U-235 and 3837 b for B-10.
CINDER uses Equation (2-141) to modify its library cros. sections for its calculations, and uses it again to convert the lumped fission product thermal cross section back to the ';30 m/sec reference velocity for use by the rest of LATTICE since computation of the spectra that determines the thermal velocity is repeated for each depletion step. In addition, the epithermal cross section calculated by CINDER is too large because the effect of resonance interaction was neglected (Reference 2-15). A cor rec-tion factor dependent on H/M is used to reduce the fission product cross section in the epithermal spectrum calculation (Section 2.4).
O 2-66
CllAPTER 2 REFERENCES
- 1. J. L. Carter, "lIRG3: A Code for Calculating the Slowing-Down Spectrum in the P yand B Approximation,"
3 BNWL-1432, 1970.
- 2. "SUPERTOG - Data Generator: Fine Group Constants and P Scattering Matrices from ENDF/B," ORNL-TM-2679, 1973.
- 3. 11. C. Honeck, " THERMOS - A 1hermalization Transport Theory Code for Reactor Lattice Calculatinns," BNL-5826, 1961.
- 4. N. M. Greene, et. al., "XLACS: A Program to Produce Weighted Multigroup Neutron Cross Sections from ENDF/B," ORNL-TM-3646, 1972.
- 5. H. C. Honeci and D. R. Finch, " FLANGE-Il - A Code to Process Thermal Neutron Dats from an ENDF/B Tape," DP-1278, 1971.
- 6. B. R. Leor.a rd , J r. , et. al., " Cross-Section Standardization for Thermal l'ower Reactors," EPRI-221, 1975.
- 7. A. Sauer, " Approximate Escape Probabilities," Nucl. Sci. & Eng.,
Vol. 16, pp. 519-335, 1963.
- 8. H. Gellings and A. Sauer, " Description of Dancoff Junior Program,"
ORNL-tr-4158, 1963.
- 9. K. M. Case, F. deHoffmann, and G. Placzek, Introduction to the Theory of Neutron Diffusion, Volume 1, Los Alamos Scientific Laboratory, Los Alamos, New Mexicc, 1953.
- 10. J. L. Carter, "HRG3: A Code for Calculating the Slowing-Down Spectrum in the Py or B y Approximation, Appendix B," BNWL - 1432, 1970.
- 11. F. T. Adler, G. W. Hinman, and L. W. Nordheim, "The Quantitative Evaluation of Resonance Integrals," paper No. P/1988, Proceedings.
of the Second United Nations International Conference on the Peaceful Uses of Atomic Energy, Vol. 16, pp. 155-171, United Nations, Geneva 1958.
- 12. Y. Ishiguro, " Exact Treatment of the Resonance Absorption of Neutrons of Intermediate Energy," Nuc. Sci. & Eng., Vol. 32, pp. 422-425, 1968.
- 13. J. Adir and K. D. Lathrop, " Theory of Methods Used in the GGC-3 Multigroup Cross Section Code," GA-7156, 1967.
- 14. J. L Carter, "IIRG1: A Code For Calculating the Slowing-Down Spectrum in the P, or B) Approximation, Appendix A," BNWL - 1432, 1970.
- 15. R. L Ifellens, "The Physics of PWR Reactors," CONF-720901, pp. 3-22, 1972.
- 16. H. S. Cheng, C. T. McDaniel, and A. Leonard, "A Nodal Integral Transport Method for Calculation of Two-Dimensional Power Distribution in Non-Uniform Lattices," CONF - 710302, pp. 613-678, 1971.
2-67
CHAPTER 2 REFERENCEE (continued)
- 17. II. C. Honeck, " THERMOS: A Thermalization Transport Theory Code for Reactor Lattice Calculations," BNL-5826, 1961.
- 18. P. Lemberg, "The RODWORTH Code - A Blackness Theory Calculation of Ef fective Few Group Dif fusion Theory Parameters for Cruciform and RCC Shaped Control Rods," UNC-5228, 1968.
- 19. M. Michelini, "Neutran Transmission Probabilities Across Control Blades Filled with Round Tubes: Formulation and Accuracy," Nucl.
Sci. & Eng., Vol. 42, pp. 162-170, 1970.
- 20. W. Sha, "An Analysis of Reactivity Worth of the Rod Cluster Control Elements (RCC) and Local Water Hole Power Density Peaking,"
WCAP-3269-47, 1965.
- 21. A. F. IIenry, "A Theoretical Method for Deterre.,ng the Worth of Control Rods," WAPD-218, 1959.
- 22. G. W. Stuart, " Multiple Scattering of Neutrons," Nucl. Sci. & Eng.,
Vol. 2, pp. 617-625, 1957.
- 23. G. W. Stuart and R. W. Woodruf f, " Method of Successive Generation "
Nucl. Sci, & Eng. Vol. 3, pp. 339-373, 1958.
- 24. T. R. England, " CINDER - A One Point Depletion and Fission Product Program," WAPD-TM-334, 1962.
O 2-68
- 3. TWO-DIMENSIONAL FLUX AND POWER DISTRIBUTIONS In order to obtain a more accurate calculation of the reactivity of the assembly and the flux distribution for the depletion study, a two-dimensional few group diffusion theory calculation for the assembly is performed. Then the rod-by-rod relative power distribution is calculated from the flux distribution.
3.1 TWO-DD1ENSIONAL MULTIGROUP DIFFUSION EQUATIONS A numerical solution must be found ts the diffusion equation for each broad energy group:
-V D(r)V$(r) + 1(r)$(r) = S(r) (3-1) where the space dependent parameters D = the ditfusion coefficient, 1 = the total macroscopic removal cross section,
$ = the flux distribution, and S = the source distribution.
The removal cross sections in the nonthermal energy groups include both absorpti,n and scatter out of the group. In two-dimensional ectangular geometry, the dif fusion equation is D + I$ - S = 0 (3-2)
Te determine the flux distribut.on, the assembly is divided into an array of mesh points, and Equation (3-2) is integrated over the region asscciated with each mesh point (the dashed region in Figure 3-1):
D -
D + I$ - S dx dy = 0 (3-3) 3-1
O i,J + :
R t i i.j ) R(i,j) g ---
f--q e hy.j l
4 3 l hy g I
I
'-' I he i1 ,d i + i,;
i 1 2 Y;I-I a
L____ _ ___ _ .J h x , i-i b hx,, C h y,j.I 2 2 R (i-i.j - s ) R (i .j -l)
, i p-l hx,i-l hx,e FIGURE 3-1: MESH FOR THE FIVE-POINT DIFFERENCE EQUATION O
3-2
The first two terms can be written as D ~ ~
Dg dx dy =
-a a 3
-L.D.. $+ 1 dx dy (3-4)
IJ IJ \ 3x2 Dy2 /
assuming each region R(i,j) is homogeneous and therefore has constant properties. Transform the surface integrals to line integrals using Green's Theorem:
-g dx dy = (Pdx + Qdy)
Using Q = S$ $
ax '
P = Dy yields
, o b+b dx dy =
D
__Q dx + E* dy\ (3-5)
(8x Oy / c Y /
Substituting into Equation (3-4) and integrating around the dashed contour in each region gives 9 ,
1
- 2 D..
i,j '3 fb
\Dx
+ By - / dx dy = D.1- 1, J - l 2
a 0Y
$ dx + D.
l '3 -l b
S$dx Y
d a f
-D. .
1, j - 1
$ dy -
dx D..
tj Ox dy + D.. S$ dx (3-6) tj Y d e
+ D. U dx - D. $ dy - D. U dy 3*
1 - 1, j. Dy 1-1,j Dx i - 1, j.-1 g h Note that the integrals along the contours between regions cancel due to the continuity of current at the material interfaces. Approximating the integrals in Equation (3-6) as 3-3
B A
B 3 dy = Ay A
with a backward dif ference for the average derivatives, Equation (3-6) becomes
- I D + dx dy = (3-7) i,j ij Dx Oy "ij Di +1,j " i - 1,j Di - 1,j + ij D ij C
ij Di ,j+1 C i,j-1 k ,j-1 i
where a = - 1-- [D h \
ij 2h ( i,j-1 y,j-1 + Dijh y,j)
I h j+ Dg,3 h ij - 2h ij ,)
+ C
'ij ij i-1,j ij i ,,j - 1 The last two terms of Equation (3-3) are (lQ -S) dx dy = ; 1 $ dx dy - S dx dy (1-8)
The integrals can be approximated by
$ dx dy = $ AxAy S dx dy = 5 AxAy where $ and 5 represent averages in region R(i,j). Then Equation (3-8) can be written as 2 l ij $ dx dy - S dx dy =
b'$
i
+ f (3-9) 1,J O
/
3-4
where
- 'i ij Y'j + i'j~l Y'j ~ l b 'j -
i 2 2 2 .
h . I . 1, j. h . I. h .
, x,1-1 1- y,j , 1 - 1, j.-1 y,J-l 2 2 2 h . 5.lj. h .
5.l ij-. l h .
y,j-1 g = _ x,1 y,j ,
ij 2 2 2 x,i-1 fi-1,j y,j , i- 1, j - 1 y,j-1 2 \ 2 2 /
The five point difference form of Equation (3-3) is the sum of Equations (3-7) and (3-9):
+ + '
"ij i i+1,j ^i-1,j @i-1,j ij @ij +C ij i,j+1
- 'i , j - 1 *i ,j - 1 + fij =0 (3-10) wh e re b'.' . = b . . + b ' .
lj lj IJ Equation (3-10) is solved for each group, but th- form of the source is dependent on the group. If all fission neutrons are born in the highest energy group, and if only downscattering between groups is allowed, then the average source for group g in each region can be written as f G l
I ol E $6. g=1 (3-lla) l g=1 ij #3 ggij -_h IE~ $$ g/1 (3-1lb) r i
gj lj where G = the number of energy groups, k = the multiplication factor of the assembly, eI E = the macroscopic nu-fission cross section for group g IJ in region R(i,j), and 3-5
Ig-1 r..
= the macroscopic removal cross section out of group g-1 IJ into group g in region R(i,j).
Examination of Equation (3-10) and its associated coefficients indicates that parameters for a non-existent region external to the system boundary are required to calculate the flux at the boundary.
Mirror (zero-current) boundary conditions are assumed:
U =0 U =0 Dx g '
ay g x y hence the properties of the imaginary region must be chosen so that the above conditions are satisfied when Equation (3-10) is solved for points on the boundary. The simplest choice for such an external medium is a vacuum. Then all required external region material properties for the calculation of the coefficients in Equation (3-10) on the boundary are zero. The coefficient of the flux outside the boundary becomes zero, eliminating the outside flux from the equation. The re fo re , the external flux can be set equal to zero and the coefficient need not be calculated.
The other coefficients now depend only on regions within the boundary.
Since each term in a coefficient is the product of a region material property and that region's mesh spacing, the external region mesh spacing can also be set equal to zero. Thus any material properties may be used for the imaginary region with the same results.
O 3-6
3.2 MET 110D OF S01.UTION Since the fission source, Equation (3-11a), is a function cf the flux, a two part iteration scheme is necessary to solve the coupled set of equations formed from Equation (3-10). The fission source is held constant during inner iterations on the flux and then recalculated using the converged flux in the outer iteration to begin the process again.
To begin the calculation, an initial guess of k and the flux for each group is required. The fission source is calculated and the inner iterations for the highest energy group are performed. Then the inner iterations are repeated for the lower energy groups using the slowing down source, Equation (3-11b). A new estimate of k is made:
I 018 $f.lJ k i LbE il X IE $f.
i,j ,g "ij 13 to rec.ilculate the fission source. The process is repeated until both k and the flux converge.
In the inner iterations, the alternating direction implicit method (Reference 3-1) is used to solve the matrix of mesh equations obtained from Equation (3-10). First the equations are solved in the i direc ion using the results from the previous iteration to get an inter--diate value of the flux. Then this intermediate flux is used to solve the equations in the j direction to get the new flux and complete the iteration.
In the nth iteration, for the i direction sweeps, Equation (3-10) can be rewritten in tridiagonal form by combining the constant terms:
3-7
a.
-b 1- 1, j. $"i- 1, j. + b'.'
tj
. 69I tj
+ a..
lj $"i+1,j. + g?I ij
=0 (3-13)
O where n-1 n-1 E
ij C
i,j-1 h ,j-1 + Cij $n-1 i i,j+1 ij and $" is the intermediate flux. Equation (3-13) is solved by Gaussian elimination for each value of j. Substituting the intermediate flux into Equation (3-10) gives the equation for the j direction-c.
1,j-1 $91,j-1 + b'.' - 49. + c.. $9 .+1 + g?I* = 0 (3-14) tj tj tj 1,j ij where g?I' ij = a.1-1,j. $"i-1,j. +
a..
ij
'? +
4"1,j&1 f..
lj Equation (3-14) is solved for each value of i to determine the new flux estimate Q" Between iteratiors, the value of the flux is mou fied using linear acceleration in order to improve convergence:
$" = $n-1 + 6 @" - Q"'I) (3-15) where 6 is the over-relaxation pa rt meter.
O
., . e
4
% 44 c)S+p %g=A
_ E E-,,O~
TEST TARGET (MT-3) s$
1.0 g m iga P "3 Olip=
cw l,l bU bb
.8 j l.25 1.4 1.6
< 6" =g MICROCOPY RESOLUTION TEST CHART
- r % //p
%ef f, %y 5rDp 44+W-
+>
4r" S.
,. S O
<h'&
'V
& f/p' "hi;>'hk NNN eg
' #4bh,
(>x+ > _ E E_ _
TEST TARGET (MT-3) 1.0 g m alg d
u [fl$$M lM 1.8 1.25 j IA 1.6
< 6~ >
MICROCOPY RESOLUTION TEST CHART 4*#% '# ++4
'hVk
y,,,
%,,?O
3.3 RELATIVE POWER DISTRIBUTION The rod-by-rod relative power distribution in the assembly can be calculated af ter the two-dimensional flux distribution has been determined as outlined above. Neglecting the contribution of the gamma ray source, the rod power is proportional to the fission rate in the cell-G N IK 8 8 P.=V. I (3-16) 1 1 g=1 n=1 " 1'ni $*
where Vg = the volume of cell i, K = the energy release per fission in nuclide n, IE = the cell macroscopic fission cross section of nuclide n, nl in group g, and
$f = the cell average flux in group g.
The relative power is obtained by dividing by the average rod power of the assembly:
P.
P. =
1 I (3-17) 11 I .
P.
I t=1 where I is the number of fuel rads in the assembly.
3-9
3.4 MESH DESCRIPTION The model used for a typical BWR assembly with the control rod inserted is shown in Figure 3-2. Each region boundary corresponds to a mesh line; additional mesh lines are included to improve the accuracy of the calculation.
All lattice locations are treated as separate regions because of the differences in spectrum and composition. Cells containing burnable poison are divided into a discrete pellet region surrounded by homoge-nized moderator and cladding. In non-burnable poison cells, the fuel, cladding, and moderator are homogenized into one region with two additional mesh lines within the region. As an option, the fuel cells can be modeled the same way as burnable poison cells. Extra in-channel water not associated with the unit cells is included in the cells adjacent to the channel.
The channel is homogenized into the water gaps outside the assembly.
Each of these regions has two equally spaced internal mesh lines except in the control rod slot when the control blade is present. Then the two mesh lines are spaced to coincide with the blade structure (Figure 3-2).
The symmetry of a typical PWR assembly allows the use of a one-quarter assembly model (Figure 3-3), but the option to model a full assembly is retained. The burnable poison cells are treated the same way as for a BWR, and the control rods are modeled the same way as burnable poison rods. The homogenized fuel cells have one internal mesh line instead of two except when a burnable poison or corErol rod cell is in the same row and/or col..mn. Since the control rods are within the lattice, the water gap between the assemblies is narrow and requires.no internal mesh lines.
3-10 -
I 2 3
4 fl.
5 7
y i i ij iii i i ii gi ii. g i i .;
l is e i i O 5 to 15 20 25 30 1 -CONTROL BLADE STRUCTURE 2- ACTIVE CONTROL BLADE 3- HOMOGENIZED EX-CHANNEL WATER AND CHANNEL 4 -HOMOGENIZED FUEL CELL 5-BURNABLE PolSON 6-HOMOGENIZED MODERATOR AND CLADDING 7-WATER HOLE FIGURE 3-2: MODEL FOR BWR ASSEMBLY 3-11
O 1 -
1 3 2
N punamagemm 6
m O
R I I ii i. i, , , . i, . . . i .
0 5 to 15 20 1-WATER CAP 2-HOMOGENilED FUEL CELL 3-8URNABLE POISON OR CONTROL ROD 4-HOMOGENIZEE MODERATOR AND CLADDING FIGURE 3-3: MODEL FOR M PWR AS$EMBLY O
'l-12
CHAPTER 3 REFERENCE
- 1. D. W. Peaceman and H.11. Rachford, Jr. , "The Numerical Solution of Parabolic and Elliptic Differential Equations,"
J. of Soc. Indust. Appl. Math., Vol. 3, pp. 28-41, 1955.
3-13
e
- 4. DEPLETION CALCULATIONS The lattice physics methods which have previously been described apply only to a fixed point in time. However, during operation of the reactor, the fuel concentrations in the lattice change and the con-centrations of strongly absorbing fission products increase. The changes in the nuclei concentrations are accounted for by solving the depletion equations. The basic assumptions used in solving these equations are that the power density and the activation rates (absorp-tion, fission and capture) remain constant during the time step.
4.1 ISOTOPIC DEPLETION EQUATIONS The depletion equations are a coupled set of first order differential equations which describe the variation with time of the number densities of the nuclei in the lattice. The time rate of change in the concentration of nuclide i at position r can be expressed in terms of " source strengths", s. . ,
IJ anel an " absorption strength", a j.
M S N.(r,t) = 1 s..(r)N.(r,t) - a r)N.(r,t) i= 1,2,. .,M dt i ij j i i j=3 (4-1) where N j (r,t) is the concentration of nuclide i at time t and position r and M is the maximum number of nuclei.
The source strength of a nuclide i from nuclide j is the sum of contributions from all possible mechanisms, i.e., fission, decay, and neutron-capture. Thus the general expression for the source strength may be written as s..(r)
IJ
= Y.J+1.A((r)+A..+A?(r)
J IJ J (4-2)
, 4-1
where O
v = the yield .: nuclide i from the fissioning of nuclide j, Ak = 2 O Ef 8 J 8 j$ , the fission activation rate, 0 = the microscopic fission cross section of nuclide j in group g,
$E = the ilux in group g, A. . = the decay constant for production of nuclide i from the
- 3 decay of nuclide j, and AS=A'-Ak,thecaptureactivationratewhereA is an absorption J J J f
. J rate analogous to A..
J The absorption strength of a nuclide i is the sum of its absorption rate and decay constant:
a.(r) 1
= A"(r) + A.1 1
(4-3)
A schematic diagram of the depletion chains that are used is given in Figure 4-1. Note that the short lived and low cross section stages are neglected. Also only the fission products with large macroscopic absorp-tion cross sections are included explicity while the rest are lumped into a single pseudo (ission product (see Section 2.7).
4.2 METHOD OF SOLUTION To solve the depletion equations, Equation (4-1), the assumption is made that the activation rates (and thus the absorption and source strengths) are constant through the time step. Therefore the solution to Equation (4-1) is
-a.t -a.t a.*M N.(t) = N.(0)e +e e I s..N.(t)dt (4-4) 1 1 JO j=1 13 3 0
4-2
- 1. liepletion Chain for U-235 235 g ,
236 g 92 n,y 92
- 2. Depletion Chain for U-238 238 g 92 N
s _
n,p 239 Np 93 s p ,t =0 239
\ 241 242 94 " 240'"
I Pu 94 Pu n,y 94 n,y 94 n,y D'
o
- 3. Depletion Chain for Gd-154 154 155 156 I Gd : Gd : Gd : Gd 64 n,y 64 n,y 64 n,y 64
- 4. Fission Product Chains 235 236 240 241 92 '
92 ' 238"'
92 239 94 "' 94 "' 94 "' 94 242'"
u y u 135 p 149 p Lumped 37 61 Fission Product (CINDER) 0 u
149 Sm 62 FIGURE 4-1: DEPLETION CIIAINS REPRESENTED 4-3
where for simplicity the spatial dependence has be.en neglected. In practice, each fuel pin is treated as a separate depletable region which may have a different isotopic content.
The solution to the depletion equations may be divided into two parts: the homogeneous part as indicated by the first term of the right-hand side of Equation (4-4) and the inhomogeneous part which represents the contribution of the sources as the second term of the right-hand side of Equation (4-4). Since the depletion equations are linear in the source terms, the contributio. from each source s..N. may IJ J be calculated independently of the other sources.
For the first nuclide in a depletion chain there are no sources; thus the homogeneous part of Equation (4-4) is the complete solution:
N'(t) = N;(0)e- l' (4-5)
In continuing to solve the equations, use is made of depletion functions O in a manner similar to that described in Reference 4-1. Defining the first depletion function as
~
F)(x) = e * (4-6)
Equation (4-5) may be written as N (t) = N (0) F (a y) 3 3 7 t (4-7)
For the second nuclide N in depletion chain with nuclide N) as 2
source, the homogeneous part of Equation (4-4) is N *(t) = N 2(0)*F1(a2 t) (4-8) and the inhomogeneous part due to the contribution of the source s 21N g(t) is O
4-4
. -a t baI N'"(t) = c 0" *21 1 (I)(IT
-a t a T 2
=e e s 0 21 1( 1("1 }
-a g t -a 2'
= Ny (0) s g -
a 2 - "I F (a t) - F (a 3 g 7 2t.)
= N (0) sg t -
(4-9)
,2 , ,1 Defining the second depletion function as F (x ) - F 3(x2) 7 3 F2 (*1'*2)
- x -x (4-10) 2 1 Equation (4-9) becomes Nf"(t) = Nj (0) s 21 t F2 ("It,a 2) t (4-11)
Therefore the complete solution for the second nuclide in a chain with a single source s N is 21 1 N (t) = N *(t) + N "(t)
= N2 (0) F ( 2 ) + N;(0)gs tF(
1 t
2 l t, a ,, t ) (4-12)
For M different sources s N the general solution is 2j M
N2 (t ) = N2 (0) F;(a.j t) + lE N (0) s t F2 ("jL'"2 ) t (4-13)
J:1 Consider now the third nuclide N 3
na ep e n cli n wi i e secon(1 nuclide N as the source 2 Using Equation (4-12) with Equation (4-4) yields
-a t 3
h I N
3 3( }' 1( 3 ** do
- 3 8 32 2(I I
=N( 3 1("3
+
2(0) s32 2("3 ' 2}+
Ng (0) s g t s3I't F 3(*3 ' 2 l') ( "I )
4-5
" " I "E"" " """ "' " I"" " "*
where F3 (*1'*2'*3}
Likewise it can be shown that for the k-th nuclide in a depletion chain with the (k-1)-th nuclide as source, the solution to the depletion equation is (14eference 4-1) k k-1 Nk (t) = I .
N.(0) j F k+1 j.(a k t,. ,,a 1
t)
- II (s..t) ij (4-15)
J:1 1J For cases where there are more than one source t.crm between some of the nuclei in the chain, new chains have to be created following the new sources, and the contributions from all the branches of the different chains are added to get the general solution for each nuclide as was done for k = 2 (Equation (4-13)).
In evaluating the depletion functions use is made of the recurrence relation (Reference 4-1)
Fk-1(*k *3'*l) - F k-1(*k *3'*2 (4-16)
F k(
- k '
- 1',- l ' ' *l) = -
x 2
-x 1
f or k > 2.
It should be noted that in evaluating the depletion functions, the numerator and denominator of one of the functions may become small at the same time, thereby resulting in loss of information. To combat this possibility a series expansion for the exponentials is used:
- n e =1 n=0 h
Using this expansion the second and third depletion functions become m (x,,-x )n-1 3
F2 (*2'*l) = jF2 (x n=1 )1- -
(4-17)
=
1 1 n-1 3 *3'*2'*l) = F;("3 x
~
( 3 *1
- ~
- 3 *2)n-1 2 '*1 4-6
Also, in order to avoid other numerical problems, the fission products which have large absorption cross sections or decay constants are treated as though their concentrations were in steady state. The steady state condition implies that the time rate of change in the nuclide density is zero:
gfN(t)=0 f
(4-19)
Then Equation (4-1) can be solved using Equat ion (4-3)
M I s..N.(t)
IJ J J"I N.(t) = (4-20)
A^ + A.
1 1 1he depletion equations are solved at the end of each time step.
Due to the high absorption rate in the barnable poison pins, a predictor-corrector method is used to obtain a better estimate of the burnable poison number densities. For each time step the thermal neutron spectrum and flux distribution calculations are repeated. However, because of the slow change in the fast spectrum with exposure, the fast spectrum and flux distribution calculation is repeated only when the exposure interval from the previous fast calculation reaches a predetermined value.
4.3 NORMALIZATION The depletion equations are solved for the fuel or super cell representing each spectrum region in the one-dimensional transport calculations to account for spectrum effects and the radial dependence of the nuclide concentrations across the cell. Each pin is depleted separately in the two-dimensional calculation of the assembly to deter-mine average cell concentrations which include the effect of interaction 4-7
of the pin with other assembly regions. For the results of these calcula-tions to be consistent, they must be normalized to one another. Since the two-dimensional depletion more accurately models the conditions in the pins, the one-dimensional cell results are normalized to agree with the assembly calculation.
The normalization of the one-dimensional results is in three parts.
First the activation rates for each cell are normalized to the average activation rates calculated for the time step by two-dimensional diffusion theory for the pins in the spectrum region represented by the cell. Then the depletion equations are solved.
Further normalization of the activatior ates is required since the powe r density P0 is assumed to be constant. The power density is defined as I
P 2 0
I K.N.(t)A[.
ii i (4-21) gg where M = the number of fissionable nuclei, F
K = the energy released per fission of nuclide i, and Af=Io 8 f
QE , the fission activation rate for nuclide i.
Since the nuclide concentrations changed and the activation rates were assumed constant during the depletion step, the power density is no longer the same as it was before the depletion. Therefore, the absorp-tion and fission activation rates are normalized by the ratio of the power before depletion P0 and the power after depletion P :
P A (t) = A (0) S i i k = a,f (4-23)
P' O
4-8
This completes the normalization of the activation rates; however, the nuclide concentrations must also be normalized. This is done in a similar manner as the first normalization of the cell activation rates.
The nuclide concentrations from the assembly depletion are assumed to be correct and the cell concentrations are normalized so that they match the average values for the pins in the spectrum region represented by the cell. Once these normalizations have been completed, the two sets of depletion calculations are consistent.
4-9
Cl{ APTER 4 REFERENCE O
- 1. II. Siewers, "An Analytical Method for Solving Depletion Equations,"
Atomkernenergie (ATKE) Bd. 27, Lfg. 1, pp. 30-34, 1976.
O L
O 4-10