Argonne National Laboratory
Reduced Enrichment for Research and Test Reactors
Nuclear Engineering Division at Argonne
U.S. Department of Energy

RERTR Publications:
Analysis Methods for Thermal Research and Test 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,

yII(0,m) = yI(0,m), m >0

yII(t,m) = yIII(t,m), m<0

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). Thus,

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) and yn(t,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 ). Thus,

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 the equation

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


Hence, , , , .

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 becomes


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 5.

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 coefficients.

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 b(P5.).

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. Thus,


fn = C cosh kxn

fn+1 = C cosh k(xn + h)

fn-1 = C cosh k(xn - h) .

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 becomes

Deff .

For greatest accuracy the effective diffusion parameters should be evaluated using the spectrum-weighted blackness coefficients <a(P5)> and <b(P5)>.

2016 RERTR Meeting

The 2016 International RERTR Meeting (RERTR-2016) will take place in Belgium. Stay tuned for further details.

2015 RERTR Meeting

The 2015 International RERTR Meeting (RERTR-2015) took place in Seoul, Korea on Oct. 11-14, 2015.
For more information visit RERTR-2015.


ARGONNE NATIONAL LABORATORY, Nuclear Engineering Division, RERTR Department
9700 South Cass Ave., Argonne, IL 60439-4814
A U.S. Department of Energy laboratory managed by UChicago Argonne, LLC

Last modified on July 29, 2008 11:33 +0200