Analysis Methods for Thermal Research and Test Reactors
COMPUTING CONTROL ROD WORTHS
IN THERMAL RESEARCH REACTORS
Blackness Coefficients and Effective Diffusion Parameters in Slab Geometry
Mesh-dependent effective diffusion parameters are often used to calculate control rod worths within the limitations of diffusion theory. For slab absorbers these effective diffusion parameters are functions of a pair of blackness coefficients. This appendix illustrates a method for calculating the blackness coefficients in the P1 approximation with and without scattering in the absorber slab and discusses how the method may be extended to the P3 and P5 approximations. Equations are derived for the effective diffusion parameters as functions of the blackness coefficients and the mesh interval size within the absorber.
A.1 Matching Conditions on the Absorber Surfaces
Consider the three-region slab configuration shown below.
The absorber region (Region II) extends from x=0 to x=t and is assumed to be a source-free region which scatters neutrons isotropically. Since the angular fluxes y(x,m) incident on Region II from Regions I and III must be continuous at the boundaries,
where is the cosine of the angle between the flux direction and the normal
to the slab. The boundary fluxes in Regions I and III may be expanded into
a power series over the range of m (-1 to 1).
Maynard (Ref. 1) has shown that the angular flux continuity requirement
leads to the matching conditions
if the absorber (Region II) is source-free and scatters neutrons isotropically.
It turns out that either the even or odd moments can be matched. The odd
moments are usually chosen. Thus, only odd values of m are used in these
matching equations with mmax = L. Rmn and Tmn
are reflection and transmission coefficients defined by the equations
Rmn and Tmn are the reflected and transmitted contributions
to the outgoing mth moments due to the incoming flux. They are
defined to be positive and are functions only of the absorber thickness
(t) and its macroscopic scattering and absorption
cross sections (Ss and Sa).
A.2 Evaluation of the Reflection and Transmission Coefficients
For the purpose of evaluating Rmn and Tmn it is
convenient to assume that Regions I and III are voids and that a mn
(n=0,1,2,...) source distribution is incident on the absorber slab from
the left. For this source distribution the monoenergetic one-dimensional
Boltzmann transport equation is solved for the surface angular fluxes yn(0,m)
Using the one-dimensional option of the TWODANT code8 (i.e. ONEDANT),
an angular quadrature of 24 (S24), and double PN quadrature
constants, Rmn and Tmn can be readily integrated by
the Gauss-Legendre quadrature method. Thus,
where N is the angular quadrature order (SN) and Wi
are the required Gauss-Legendre weights which are normalized so that in
the range from i=1 to i=N/2 they sum to unity. The abscissas (mi)
and weights (Wi) are given below for N=24.
For the special case of a pure absorber slab (Ss
= 0) Rmn is zero and Tmn can be easily evaluated.
For this case the transmitted angular flux is the product of the incident
flux mn and the probability of a neutron
passing through the slab without absorption (e-Sat/m
where Em+n+2 (Sat)
is the exponential integral of order m+n+2.
A.3 Evaluation of the Blackness Coefficients
In the regions outside the control blade the angular fluxes satisfy the
monoenergetic, one-dimensional, time-independent Boltzmann transport. To
obtain analytical expressions for the blackness coefficients the angular
fluxes on the surfaces of the absorber are expanded into a series which
in the PL approximation becomes
where yn(x) is the nth spherical
harmonic moment and where Pn(m) is
the nth Legendre polynomial. The spherical harmonic moments are given by
Although the higher order moments are more complicated, y0(x)
and y1(x) are just the neutron flux
and neutron current, respectively.
In the P1 approximation (L=1=m) the angular fluxes on the
left-hand and right-hand surfaces of the absorber plate become
Note that B1 is negative because the blackness coefficients
are defined in terms of currents into the absorber slab. Adding the
Maynard matching equations for the L = 1 case gives
Using the above values for (A0 + B0) and (A1
- B1) this equation becomes
The blackness coefficient in the P1 approximation is therefore
Using the same procedure but subtracting the Maynard matching equations
gives the blackness coefficient.
Note that the b-equation can be obtained from
the a-equation by simply changing the sign of
the transmission coefficients. This is a general result which is also valid
in the P3 and P5 approximations.
As discussed above, for the case of zero scattering the reflection coefficients
vanish and the transmission coefficients reduce to exponential integrals.
Thus, the no-scattering P1 approximation for the blackness coefficients
Recall that improved values are obtained for these blackness coefficients
by using the modified forms (a0m and
b0m) in which the (1/2) coefficient
is replaced with the value 0.4692.
These P1 approximations for the blackness coefficients show
that they are functions only of the properties of the absorber slab (Sa,
Ss, and t).
This is also true for the higher order approximations. Although the algebra
is tedious, these same methods can be used to find the blackness coefficients
in the P3 and P5 approximations. However, for these
cases the spherical harmonic moments yn(x)
need to be evaluated on the surfaces of the control slab for n=2,3,4 and
To calculate these higher-order spherical harmonics it is assumed that
the control blade is surrounded by an homogenized fuel zone. For this medium
the monoenergetic, one-dimensional, time-independent Boltzmann equation
is multiplied by (2n + 1)Pn(m) and
integrated over all m. Using the recurrence relation
(n + 1)Pn+1(m) + nPn-1(m)
= (2n + 1)mPn(m)
this integration leads to a set of first-order coupled differential equations
for the spherical harmonics which, in the PL approximation, are
subject to the requirement that yn(x)
= 0 for n>L. Solutions to these differential equations are of the form
yn(x) = gn(n) enx/l
where l is the total mean free path of neutrons
in the homogenized fuel regions outside the absorber slab. Substituting
this solution into the differential equations for the spherical harmonics
leads to a recurrence relation for the gn(n)'s
with g0(n) defined to be unity. In
the PL approximation yL+1(x)
= 0 which requires gL+1(n) = 0. This
condition leads to the allowed values of n which
are positive and which are (L-1)/2 in number. Recall that in the PL
approximation only odd values of L are used. Since the blackness coefficients
are functions only of the properties of the absorber slab, it is necessary
to assume that the surrounding fuel regions are infinite in extent (no leakage)
with either isotropic or linear anisotropic scattering.
Following these procedures Ref. 3 shows that
g0(n) = 1
g1(n) = 0
g2(n) = -1/2
g3(n) = 5/(6n)
g4(n) = 3/8 - 35/(24n2)
g5(n) = -[9g4(n) + 4 n g3(n)]/(5n)
g6(n) = -[11g5(n)
+ 5 n g4(n)]/(6n)
In the P3 approximation there is one allowed n-value
since (L-1)/2 = 1 and this value is obtained from the requirement that g4
= 0. Thus, n2(P3) = (35)1/2/3.
Similarly, in the P5 approximation n2(P5)
= 1.2252109 and n3(P5)
= 3.2029453. These values of n and gn(n)
determine the higher order spherical harmonics evaluated on the absorber
slab surfaces. Thus,
Spherical Harmonics in the P3 Approximation
y0l = fl + a2
y0r = fr + b2
y1l = Jl
y1r = - Jr
y2l = g2(n2) a2
y2r = g2(n2) b2
y3l = g3(n2) a2
y3r = - g3(n2) b2
Spherical Harmonics in the P5 Approximation
y0l = fl + a2 + a3
y0r = fr + b2 + b3
y1l = Jl
y1r = - Jr
y2l = g2(n2) a2 + g2(n3) a3
y2r = g2(n2) b2 + g2(n3) b3
y3l = g3(n2) a2 + g3(n3) a3
y3r = - g3(n2) b2 - g3(n3) b3
y4l = g4(n2) a2 + g4(n3) a3
y4r = g4(n2) b2 + g4(n3) b3
y5l = g5(n2) a2 + g5(n3) a3
y5r = - g5(n2) b2 - g5(n3) b3
where a2, a3, b2, and b3
are arbitrary constants which are eliminated in the evaluation of the blackness
With these values for the spherical harmonics the blackness coefficients
can be calculated in the P3 and P5 approximations
following the same approach as was used above for the P1 approximation.
As before, the continuity requirement of the angular fluxes at the surfaces
of the absorber slab determines the expansion coefficients An
and Bn in terms of the spherical harmonics ynl
and ynr. To determine the blackness
coefficient a in the P3 approximation,
for example, Maynard's matching equations are added together for m=1 and
for m=3. From these two equations the constant (a2 + b2)
is eliminated and from the resulting equation a(P3)
is determined. By subtracting the equations and eliminating the constant
(a2 - b2), b(P3)
is obtained. b(P3) can also be found
by simply changing the sign of all the transmission coefficients (Tmn)
in the final expression for a(P3).
Although the algebra is very tedious, the same procedure is used to find
a(P5) and b(P5).
For this case Maynard's two matching equations are added together for m=1,
m=3, and m=5. From these three equations the constants (a2 +
b2) and (a3 + b3) are eliminated and the
resulting equation solved for a(P5).
Changing the signs of the reflection coefficients in this equation determines
For more details regarding these procedures see Ref. 3. As for the P1
case, the P3 and P5 blackness coefficients are functions
only of the properties of the absorber slab (t,
Sa, and Ss).
A.4 Evaluation of the Effective Diffusion Parameters
The effective diffusion parameters are chosen so as to preserve the current-to-flux
ratios on the surfaces of the control slab as given by the blackness coefficients.
Since these effective diffusion parameters are to be used in a finite difference
solution of the diffusion equation, they contain an explicit dependence
on the mesh interval size, h. This allows the use of a very coarse mesh
in the absorber for diffusion calculations. Equations will be derived for
the effective values of D and Sa for
use in those diffusion codes, such as DIF3D4, which evaluate
fluxes at the center of mesh intervals. The same methods may be used to
determine the effective diffusion parameters for codes which evaluate fluxes
on mesh interval boundaries. These techniques for determining mesh-dependent
effective diffusion parameters as functions of a
and b were first proposed by E. M. Gelbard.17
Consider the following diagram.
It is convenient to assume that the same material extends to regions
outside the absorber slab of thickness t. Since
a and b depend only
on the properties inside the slab, this assumption leads to no loss in generality.
If the flux varies linearly from the center to the edge of the mesh cell,
fl = (f-1
+ f1)/2 and Jl = (D/h)(f-1
- f1) .
The solution to the diffusion equation within the source-free absorber
consists of a symmetric (cosh kx) part and an asymmetric (sinh kx) part
with respect to the centerline. For the symmetric solution, fl
= fr , Jl = Jr
, f1 = C cosh k(t
- h)/2, and f-1 = C cosh k(t
+ h)/2 so that after some manipulation
Similarly, for the asymmetric solution, Jl = -Jr,
fl = -fr,
f1 = A sinh k(t
- h)/2, f-1 = A sinh k(t
+ h)/2 so that
The ratio of these equations gives
from which it follows that
The expression for the effective diffusion coefficient is obtained by
adding the equations for a and b.
An expression for the effective macroscopic absorption cross section
is obtained by writing the diffusion equation for the control slab in finite
difference form and solving for Sa.
fn = C cosh kxn
fn+1 = C cosh k(xn + h)
fn-1 = C cosh k(xn -
These equations for k, D, and Sa
determine the effective diffusion parameters in terms of the blackness coefficients
(a and b) and the
mesh interval size h for the case of mesh-centered fluxes.
The same procedures can be used to find the effective diffusion parameters
for diffusion codes which evaluate fluxes on mesh boundaries. The results
for k and Sa are the same as those
given above. However, the expression for the effective diffusion coefficient
For greatest accuracy the effective diffusion parameters should be evaluated using the spectrum-weighted blackness coefficients <a(P5)> and <b(P5)>.