Nuclear Criticality Safety Engineer Training

Module 6 1

Diffusion and Transport Theory - Part I

LESSON OBJECTIVES - PART I

To introduce the concepts of radiation transport, the neutron balance equation, approximations to
the solution of the transport equation including the diffusion theory approximation, and energy
group structure.

REFERENCES

There are many textbooks and references for radiation transport and diffusion theory. Two
widely used references are listed below, and the material in this module is taken from both of
them.

1) Duderstadt, J.J. and Hamilton, L.J., Nuclear Reactor Analysis, John Wiley and Sons,
Inc., New York (1976).

2) Lamarsh. J.R., Nuclear Reactor Theory, Addison-Wesley, Reading, MA (1972).

NEUTRON TRANSPORT

In the introduction to Chapter 4, Duderstadt has some well-worded comments about neutron
transport theory. At the risk of taking quotes out of context, he states:

"... the so-called neutron transport equation ... usually strikes far more terror in the
hearts of fledgling nuclear engineers who are intimidated by its frightening
reputation within the nuclear reactor community. Neutron transport theory has
come to be associated with a hideous plethora of impenetrable mathematics,
unwieldy formulas, and (eventually) the expenditure of enormous amounts of
money on computer number-crunching."

"... As in other areas of physical analysis, we will find that the usual maxim applies:
that the “brute force” numerical approach (i.e., discretize everything in sight and
slap it on a computer) is conceptually the simplest approach to understand and

1
Developed for the U. S. Department of Energy Nuclear Criticality Safety Program by T.
G. Williamson, Ph.D., Westinghouse Safety Management Solutions, Inc., in conjunction with the
DOE Criticality Safety Support Group.

NCSET Module 6 Diffusion Theory 1 of 18

 computationally the most expensive calculation to perform. The more elegant
approximate methods require far less computational effort but far more in the way
of mental gymnastics in order to understand the significance and reliability of their
predictions."

Solution of the transport problem begins by writing a neutron balance equation. This process is
similar to balancing a checkbook, with entries for income and spending. The balance equation
approach is used in many academic disciplines. In fluid mechanics a mass and momentum balance
is written, resulting in the Navier-Stokes equations; in electromagnetic theory an energy balance is
written, resulting in Maxwell’s equations; in quantum mechanics an energy balance is written,
resulting in Schroedinger’s equation; in transportation a people moving balance equation is
written to derive equations for highway design; in economics balances of goods and services and
cash are written which result in government policies. In many of these disciplines the balance
equation is written at the beginning of an academic course, and the rest of the term is spent in
solving the equations, often with numerous approximations. The same is true for the neutron
transport equation.

In general, one can define a neutron density N( r ,t) such that N( r ,t)d3r is the expected number
of neutrons in the small volume d3r at a time t. This can be generalized to define the angular
neutron density that defines the state of a population of neutrons with a differential volume having
energy E and velocity v at time t. That is,

n( r ,E, Ω ,t) = expected number of neutrons in d3r about r , with energy dE about E,
moving in direction Ω in solid angle d Ω at time t

The bold, over-lined symbols represent vectors, and Ω is the unit vector of the velocity, v / | v |.
The angular neutron density has units of neutrons/(cm3 eV ster) assuming energy units of eV are
used. With this definition, the total neutron interaction rate within d3r is given by

F( r ,E, Ω , t) = v Σ(E) n( r ,E, Ω ,t)

= Σ(E) ϕ( r ,E, Ω ,t)

where ϕ( r ,E, Ω ,t) is the angular neutron flux which has units of neutrons/(cm2·s·eV·ster)
assuming energy units of eV are used. Note that n, ϕ and F are all functions of seven scalar
variables: x, y, z, θ, φ, E and t.

Two more definitions are necessary before writing the neutron balance equation. The local
neutron source term (e.g., from fission events) is defined as s( r ,E, Ω ,t). Neutrons can also
change state due to scattering events which cause them to change energy and velocity. This term
in the balance is represented by

 
 ∫ v' Σ s ( E' → E, Ω ' → Ω )n(r, E' , Ω ', t )d 3 r  dEd Ω
 V 

NCSET Module 6 Diffusion Theory 2 of 18

The neutron balance equation is now formed by requiring that the sum of the time derivative of
the neutron density, the production terms and the loss terms will be zero, giving

∂ n( r , E, Ω , t )
+ vΩ • ∇ n(r , E, Ω , t ) + vΣ t n(r , E, Ω , t )
∂t
∞ (1)
= ∫ dΩ '∫ dE' v' Σ s (E' → E, Ω ' → Ω )n(r , E', Ω ', t ) + s(r , E, Ω , t )
4π 0

This is it - the neutron transport equation! (Duderstadt Eqn 4-40)

Equation (1) balances the losses and gains of neutrons at a point in phase space which is described
by the neutron location (x, y, z), the direction of motion (θ, φ), the energy (E) and the time (t).
Neutrons move in straight lines with constant speed until they interact with a nucleus. When an
interaction occurs at location (x, y, z), the neutron energy and direction may be changed. In the
balance equation, the terms on the left side of the equality represent particles moving out of the
phase space element and those on the right side are those moving into the element. The fives
terms in the equation represent:

• the first term on the left is the time rate of change of the neutron population;
• the second term on the left is the spatial divergence of the neutron population;
• the third term on the left is the neutron interaction term;
• the first term on the right accounts for neutrons scattered from energy E’ and direction
Ω ' into the phase space (E, Ω );
• the second term on the right accounts for neutron sources, such as fission.

This transport equation is:

• a linear equation in the unknown independent variable n;

• contains seven independent variables;
• has derivatives in space and time;
• contains integrals over energy and angle; and
• requires boundary conditions and initial conditions.

The equation is most often written in terms of the angular flux (Duderstadt Eqn 4-43)

∂ϕ ( r , E, Ω , t )
+ Ω • ∇ ϕ ( r , E, Ω , t ) + Σ t ϕ ( r , E, Ω , t)
v∂ t
∞ (2)
= ∫ dΩ '∫ dE' Σ s (E' → E, Ω ' → Ω )ϕ ( r , E' , Ω ' , t ) + s( r , E, Ω , t )
4π 0

NCSET Module 6 Diffusion Theory 3 of 18

If the only source term is fission then the source term becomes

∞
χ ( E)
s f ( r , E, Ω , t ) = dΩ ' ∫ dE' ν ( E')Σ f ( E')ϕ ( r , E, Ω , t )
4π 4∫π
(3)
0

where χ(E) is the energy distribution of the fission neutrons (assuming they are all prompt
neutrons for now). All we have to do now is solve this equation. To do that we must first
discuss some approximations.

First a comment on the equation. Since the equation balances the neutron losses against the
gains, it can only be solved for a critical configuration, that is, when the neutron production
equals the neutron losses. What is normally done to get around this dilemma is to introduce an
eigenvalue, λ, as a number that modifies the number of neutrons per fissions. That is, replace
ν(E') with ν/λ and solve for the value of λ which will balance the equation. This eigenvalue is the
one mentioned in the mathematical definition of the effective multiplication factor given in
ANSI/ANS- 8.1.

TIME DEPENDENCE

For one approximation, consider the steady-state solution. That is, drop the time dependence by
eliminating the term which includes the partial derivative with time. This is a reasonable
approximation for our purposes since the majority of criticality problems are steady state. Cases
in which the time dependence is important include nuclear excursions (which should never
happen) and the analysis of subcritical pulsed-neutron source experiments.

ANGULAR DEPENDENCE

Approximations to simplify the angular dependence include neglecting any angular dependence,
the discrete ordinates approximation, polynomial expansion, and ray tracing.

For the diffusion equation the angular dependence is neglected, which can be a serious
approximation in special cases. We have seen that neutron scattering by hydrogen is not isotropic
in the laboratory, so ignoring the angular dependence will not treat that class of problems
correctly, and that certainly includes a large number of problems. Also, the neutron flux density
might not be isotropic near a boundary between two dissimilar materials. An example is a
scattering material next to a non-scattering material, such as a vacuum. The characteristic of a
vacuum is that there are no return neutrons, so in the scattering material near the boundary the
majority of the neutrons are moving toward the boundary. At this point in space an isotropic
approximation for the flux density might not be a good approximation.

NCSET Module 6 Diffusion Theory 4 of 18

In spite of these comments, diffusion theory with no angular dependence and with corrections to
the scattering cross sections is a powerful tool which gives reasonable results, particularly in
large, well-thermalized systems.

For the discrete ordinates approximation, the cosine of the polar angle, µ = cosθ, is divided into a
finite number of intervals and the angular distribution is assumed constant over each interval. The
discrete ordinates approximation has led to a large family of computer codes, commonly called
the Sn approximations, such as DTF, ANISN, XSDRN, DOT, TWOTRAN, TWODANT,
THREEDANT and TORT. The higher the level of approximation, S4, S8, etc., the better the
approximation. The Sn codes are most useful to the criticality specialist for calculating simple
geometries in one dimension, e.g., doing parametric studies for spheres, long cylinders or large
slabs.

Polynomial expansion solutions of the transport equation are not used as “production codes",
although scattering cross sections are represented as expansions in Legendre polynomials in many
codes.

The ray-tracing approximation has led to the family of codes under the name of Monte Carlo,
such as KENO, MORSE and MCNP, and some shielding point-kernel integration codes, such as
QUAD.

ENERGY DEPENDENCE

The approximation for the energy dependence is to establish an energy group structure with a
finite number of groups and assume that the cross sections are constant over each energy group.
The energy variable spans many decades from fission energy to thermal energy and over these
decades the cross sections may vary over wide ranges. There are some natural energy groupings,
such as the fission energy range, the slowing down range, and the thermal group. However, we
might want more energy detail in important ranges, for example: near important resonances such
as the 238U capture resonances for slowing down problems; more detail in the fission energy
region if we are concerned with fast fission effects; or more detail in the thermal group for a
predominately thermal system. The energy groups are normally ordered from g = 1 in the highest
energy group to g = G in the lowest energy group, so as a neutron slows down it increases in
group number.

Once a group structure has been selected the group-averaged cross sections must be calculated:

Eg −1
∫E g
dE Σ a ( E ) ϕ ( E )
Σ ag ≡ Eg −1 (4)
∫E g
dEϕ ( E )

NCSET Module 6 Diffusion Theory 5 of 18

Note the dependence of these group cross sections on the neutron spectrum, which kind of leads
to a circular argument. You need the spectrum before you can calculate the cross sections, but
you can't find the spectrum without knowing the cross sections. One choice is to iterate and
calculate the group averaged cross sections for each problem, which can be a time consuming
process. Another choice is to select a representative spectrum and tabulate cross sections for that
spectrum. Many criticality analyses are done for well-thermalized systems in which the neutron
slowing down spectrum is dominated by hydrogen scattering. Cross sections that are calculated
for a hydrogen (water) slowing-down spectrum are used in many analyses. Because the neutron
slowing-down spectrum used for the cross section generation may be different from the spectrum
in the problem at hand, it is important that a proper bias be established for the calculation.

As an example of group structure, here is the 27-group structure in the ORNL SCALE system,
which is often used in criticality studies.

Upper Upper Upper

Group Energy Group Energy Group Energy
1 20 MeV 8 100 keV 19 1 eV
2 6.434 MeV 9 17 keV 20 0.8 eV
3 3 MeV 10 3 keV 21 0.4 eV
4 1.85 MeV 11 0.55 keV 22 0.325 eV
5 1.4 MeV 12 100 eV 23 0.225 eV
6 0.9 MeV 13 30 eV 24 0.1 eV
7 0.4 MeV 14 10 eV 25 0.05 eV
15 3.05 eV 26 0.03 eV
16 1.77 eV 27 0.01 eV
17 1.3 eV
18 1.13 eV

This group structure has seven groups in the fast neutron energy range above 100 keV, eleven
groups in the intermediate range between 100 keV and 1 eV, and nine groups in the thermal
neutron energy range.

Neutrons can be transferred from an energy group (transferred into a group) by:

• Elastic scattering (up and down)

• Inelastic scattering
• Fission
• (n,2n) reactions.

NCSET Module 6 Diffusion Theory 6 of 18

Some calculation tools do not require a group structure and cross section weighting. Those
codes, e.g., some Monte Carlo codes, use cross sections commonly referred to as continuous
energy cross sections.

SUMMARY - PART I

A neutron balance equation, commonly called the neutron transport equation, has been written.
The equation is an integral-differential equation in seven variables. Approximations to simplify the
time variable, the angular variables and the energy variable have been discussed.

NCSET Module 6 Diffusion Theory 7 of 18

 Diffusion and Transport Theory - Part II

LESSON OBJECTIVES - PART II

To investigate analytic solutions to the one-speed diffusion approximation to the neutron

transport equation and introduce the concepts of geometric and material buckling. At the same
time, illustrate the limitations of the one speed diffusion theory approximation.

DIFFUSION EQUATION

In Part I we wrote an equation based on a particle balance and noted that this equation in all of its
full-blown glory is very difficult to solve exactly, and even if it could be solved exactly, the
solution could provide more information than we might want. We then proceeded to make
approximations to the equation that removed the time dependence, the angular dependence and
the energy dependence. That left an equation in the three spatial variables, which can often be
quickly reduced to one spatial variable.

In Lamarsh (Reference 2) starting on page 123 there is a neutron balance development which
leads to an equation he calls the equation of continuity. This equation relates the neutron flux
density, which is in the source and interaction terms, to the neutron current density, which is in
the leakage term. He then proceeds to relate the flux and the current, ϕ and J, with Fick's law and
says that this leads to the diffusion approximation. The following assumptions are made:

a) the medium is infinite;

b) the medium is uniform, so that all cross sections are constants, independent of position;
c) there are no neutron sources in the medium;
d) scattering is isotropic in the laboratory system;
e) the neutron flux is a slowly varying function of position;
f) the neutron flux is not a function of time.

The author then spends several pages discussing the significance of each of these assumptions and
that discussion is left for your reading enjoyment.

The one-speed diffusion equation with no time dependence and with the caveat that the diffusion
coefficient has no spatial dependence is

− D∇ 2 φ ( r ) + Σ a φ ( r ) = S( r ) (5)

where the diffusion coefficient, D, in this approximation is

Σs
D= (6)
3Σ 2t
with the subscripts s and t indicating scattering and total cross sections, respectively.

NCSET Module 6 Diffusion Theory 8 of 18

The diffusion equation is a second-order differential equation with solutions characteristic of each
coordinate system. The problem also requires suitable boundary conditions. Since the steady-
state diffusion equation is a partial differential equation, there are an infinite number of functions
which satisfy the equation. However, for a real physical situation there is only one function which
correctly represents the flux. To find this function we place boundary condition restrictions on
the problem.

BOUNDARY CONDITIONS

At the outer boundary of a region of fissile material we expect a condition that no neutrons enter
from the outside world. One way to state this condition is to impose a requirement that the
current of incoming neutrons is zero at the boundary; that is

J ( rb , Ω _ ) = 0

where r b is the system boundary and Ω_ indicates the direction into the surface. Physically this
condition requires that there be no material (i.e., a vacuum) beyond the system boundaries, but we
know that diffusion theory is not valid near a surface over which the material properties change.
To get around this problem, assume a boundary condition at the surface of the medium of the
form
1 dϕ 1
= −
ϕ dn d

where the derivative is the normal derivative of the flux and d is the extrapolation distance. For a
planar free surface it is shown by transport theory that the extrapolation distance is related to the
transport mean free path by

d = 0.71 λ tr

where λ tr = 1/Σ tr . The transport cross section, Σ tr , is a correction to the scattering cross section
to account for some of the inadequacies of diffusion theory. Then the vacuum boundary
condition becomes a requirement that the solution of the diffusion equation vanishes beyond a free
surface at the extrapolation distance,

ϕ ( rb + d ) = 0

At interior interfaces the boundary conditions are that the flux and the current across an interface
A-B be continuous, that is,

ϕA = ϕB
and

NCSET Module 6 Diffusion Theory 9 of 18

 (J A ) n = (J B ) n .
Other conditions which seem obvious but are nontrivial are that the solution to the diffusion
equation must be real, nonnegative, and single-valued in those regions where the equation applies.

DIFFUSION EQUATION IN A MULTIPLYING MEDIUM

Now include fissioning material in the system and investigate solutions in cartesian coordinates for
a one-dimensional slab of thickness a. Now include time dependence also. The fission source is
given by
S f ( x, t ) = vΣ f ϕ ( x, t ) .

If this is the only source term, the diffusion equation, with time dependence, becomes

1 ∂ϕ ∂ 2ϕ
− D 2 + Σ a ϕ (x, t ) = ν Σ f ϕ ( x, t ) (7)
v ∂t ∂x

with the initial condition that

ϕ ( x,0) = ϕ 0 ( x) = ϕ 0 ( − x) (symmetric)

and boundary conditions

 a′   a′ 
ϕ  , t = ϕ  − , t = 0
 2   2 

where a' = a + 2d, which includes the extrapolation distance.

Solve this partial differential equation by using the separation of variables method to look for a
solution of the form

ϕ ( x, t ) = Ψ ( x) T( t ) .

Substitute this into the diffusion equation and divide by Ψ(x)T(t) to find

1 dT v  d 2 Ψ 
=  D 2 + (ν Σ f − Σ a )Ψ ( x)  = constant ≡ − λ
T dt Ψ  dx 

where λ is a constant which is not known at this time.

NCSET Module 6 Diffusion Theory 10 of 18

The partial differential equation has now been separated into two ordinary differential equations:

dT
= − λ T( t )
dt
and
d2Ψ λ
2 + ( ν Σ f − Σ a )Ψ x = −
D ( ) Ψ ( x) .
dx v

The time dependent equation is easily solved:

T( t ) = T( 0) e − λt , (8)

leaving the spatially-dependent equation to be solved as shown below.

EIGENFUNCTION EXPANSIONS

Consider a homogenous problem similar to the spatial diffusion equation,

d2Ψ
2 + B 2 Ψ ( x) = 0 (9)
dx

with the boundary conditions

 a′   a′ 
Ψ   = 0 and Ψ  −  = 0
 2  2

where B is an arbitrary parameter. The general solution is

Ψ ( x) = A 1 cos Bx + A 2 sin Bx .

The boundary conditions require that

 Ba ′   − Ba ′ 
A 1 cos   = 0 and A1 cos   =0
 2   2 
and
 Ba ′   − Ba ′ 
A 2 sin   = 0 and A 2 sin   = 0.
 2   2 

NCSET Module 6 Diffusion Theory 11 of 18

One solution is to let A1 = A2 = 0, which doesn’t tell us a lot about the behavior of the system.
Another approach is to place conditions on the constant B; for example

nπ
B = Bn ≡ ( n = 1,3,5...)
a

with A2 = 0 gives
nπ x
Ψ n ( x) = A n cos ( n = 1,3,5...) (10)
a′

We also can write another set of solutions for A1 = 0 and

nπ
B = Bn = ( n = 2,4,6 ...)
a′
and
nπ x
ϕ n ( x) = A n sin ( n = 2,4,6 ...)
2

This is an eigenvalue problem with solutions which have many interesting properties, and a lot of
time could be spent playing with them. However, for now we will only consider the fundamental
mode, that is n = l.

 π x
Ψ ( x) = A 1 cos   (11)
 a′ 

If we identify this eigenvalue problem with the space dependent diffusion equation, then

1λ 
B12 =  + νΣ f − Σ a 
D v 
or
λ 1 ≡ vΣ a + vDB12 − vνΣ f

The fundamental solution of the one-speed diffusion equation in cartesian coordinates for a slab
with thickness a' (including the extrapolation distance) is

πx
ϕ ( x, t ) = A 1 exp(− λ 1 t ) cos . (12)
a′

This solution satisfies the boundary conditions, and the values of A1 can be determined from the
initial condition.

NCSET Module 6 Diffusion Theory 12 of 18

It is customary to refer to the value of B12 for the fundamental mode as the geometric buckling

2
 π
B12 =  ≡ B2g ,
 a′ 

where the factor π/a' is characteristic of cartesian coordinates.

The next step is to make the flux distribution in the reactor time independent, i.e., to make the
neutron balance steady state. Define the state of criticality as follows:

criticality : when a time-independent neutron flux can be sustained in the reactor

with no sources other than fission.

The requirement for a time-independent solution is that the fundamental eigenvalue vanish.

λ 1 = 0 = v(Σ a − νΣ f ) + vDB12

The solution becomes

ϕ ( x, t) → A1 cos B1 x ≠ function of time

with the criticality condition

νΣ f − Σ a
B12 = B2g = = B 2m . (13)
D

That is, the condition for criticality is that the material buckling, Bm2, which depends on the
material properties (cross sections and atom densities) equals the geometric buckling Bg2, which
depends on the geometry (slab width, a', in one-dimension cartesian coordinates). To achieve a
critical configuration we must adjust either the size or the core composition.

Some expressions for the geometric buckling of various geometries are listed below (in each case,
λ is the extrapolation distance and γ is the axial extrapolation distance for the cylinder.)

π2
Infinite slab of width a: B2g = (14)
( a + 2λ ) 2

π2 π2 π2
Cuboid with sides a, b, and c: Bg = + + (15)
( a + 2λ ) 2 ( b + 2λ ) 2 ( c + 2λ ) 2

NCSET Module 6 Diffusion Theory 13 of 18

 π2
Sphere of radius r: B2g = (16)
(r + λ ) 2

π2 ( 2.405) 2
Finite cylinder with radius r and height h: B2g = + (17)
(h + 2γ ) 2 (r + λ ) 2

For some criticality hand calculations it is a common practice to compare geometric shapes of the
same buckling, for example for a given material a critical sphere and a critical cube will have the
same geometric buckling.

SUMMARY

Beginning with the neutron diffusion equation in one dimension for one energy group we have a
derived a criticality relation between the size of the system and its material properties, Bg2 = Bm2.
Along the way we have looked at approximations and boundary conditions used to solve the
diffusion equation for simple geometries, and laid the ground work for more general, more
complicated solutions of the diffusion problem.

NCSET Module 6 Diffusion Theory 14 of 18

PROBLEMS

1. Figure 7 of LA-10860-MS (Critical Dimensions of Systems Containing 235U, 239Pu and 233U,
1986 Revision) shows that the effective extrapolation distance for a cylinder with h/d = 1 is about
2 cm for 94% enriched uranium. Figure 6 of the same reference shows that the ratio of the
spherical extrapolation distance to the cylindrical extrapolation distance is nearly 1 for cylinders
with h/d = 1. Table 29 of LA-10860-MS gives the critical mass of a bare sphere of U(94) as
49.12 kg at a density of 18.74 g/cm3.

Calculate the radius and geometric buckling of the critical sphere. Calculate the radius of a
cylinder of U(94) with the same buckling. Assume that the extrapolation distance applies to a
cube of the same material and calculate the length of the sides of such a cube. Calculate the
volume (V), surface area (S) and 4V/S for each geometry and compare the results.

2. A mixture of 235U and water is contained in a sphere with radius 14.20 cm. The atom densities
(atom/barn-cm) and thermal neutron cross sections (barn) are listed in the following table.

σf σa σs ν
235
U 1.2811E-04 582.6 680.9 14.3 2.425
H 6.6559E-02 0.333 20.49
O 3.3279E-02 0.00054 3.76

For this mixture calculate the thermal diffusion coefficient, D, and the material buckling, Bm2.
Equate the material buckling to the geometric buckling and calculate the radius of a critical bare
sphere. Calculate the 235U concentration in this mixture and the mass of uranium in this sphere.
Compare this value with the 235U critical mass curve (see Module 5). Note that this problem does
not give a “good” solution and points to a weakness of diffusion theory. This is discussed below
in the solution to the problem.

NCSET Module 6 Diffusion Theory 15 of 18

PROBLEM SOLUTIONS

1. Since the volume of the sphere is simply the mass divided by the density, V = 2621.1 cm3. The
radius can then be calculated as 8.553 cm. Since the ratio of the spherical to cylindrical
extrapolations distances is 1, and the spherical extrapolation distance is given as 2 cm, Equation
16 gives
π2
B2g = = 0.08862 cm −2
(8.553 + 2 ) 2

With this value of the buckling, and assuming that the axial and radial extrapolation distances for
the cylinder are the same, Equation 17 can be used to find the cylinder radius by inserting the
buckling calculated above and λ = 2 cm, then solving for r = 7.649 cm.

Since no information is given about the cube extrapolation distance, assume it is also 2 cm.
Equation 15 can then be used to determine the length of a side of the cube, a = 14.279 cm, or a
half-length of 7.139 cm.

The calculated values of V, S and 4V/S are given in the following table.

Dimension (cm) Volume (cm3) Area (cm2) 4V/S (cm)

Sphere 8.553 2621.1 919.3 11.4
Cylinder 7.649 2812.4 1103.0 10.2
Cube 7.139 2911.2 1223.3 9.52

The sphere has radius 8.55 cm, the cylinder a radius of 7.65 cm and the cube a half-width of 7.14
cm and the best surface-to-volume ratios are in the order sphere, cylinder and cube.

2. Using the diffusion theory approximation as discussed on Page 8 of this module, the diffusion
coefficient, D, is
Σs
D=
3Σ 2t
with the subscripts s and t indicating scattering and total cross sections respectively. Using the
data given in the problem, the macroscopic cross sections are

Σ s = 1.491 cm -1
Σ t = 1.600 cm -1
giving
D = 0.194 cm.

NCSET Module 6 Diffusion Theory 16 of 18

From Page 13 of this module, the material buckling is given by

νΣ f − Σ a
B2m = .
D

Using the data given in the problem,

νΣ f = 0.1810 cm -1
Σ a = 0.1094 cm -1 .

Combining these values with D calculated above,

Bm2 = 0.3689 cm -2
.

For a bare sphere of radius r, Equation 16 gives the geometric buckling as

π2
B2g = .
(r + λ ) 2
Equating the geometric and material bucklings, then solving for (r + λ) gives

(r + λ) = 5.173 cm,

the extrapolated radius of a critical bare sphere of this uranium water mixture. Using the given
atom density of uranium, the mass density of 235U in the sphere is calculated to be

ρ(235U) = 0.050 g/cm3 .

Next using the extrapolated critical radius to calculate the volume, the mass of uranium is found
to be
m(235U) = 29.0 g .

We know that this is a ridiculously low answer. One standard reference gives a subcritical mass
limit for this mixture of about 850 g (TID-7016 Rev. 2, Fig. 2.1). The problem is in the diffusion
approximation in which the hydrogen scattering is assumed to be isotropic.

A better approach is to redefine the diffusion coefficient with the transport approximation (see,
e.g., Duderstadt, p. 136):
1
D=
3Σ tr
where
Σ tr = Σ t − µ 0 Σ s
= Σ a + Σ s (1 − µ 0 )

NCSET Module 6 Diffusion Theory 17 of 18

with
Σa = Σt − Σs
and
2
µ0 =
3A
where µ 0 is the cosine of the mean scattering angle, approximated by 2/3A. Calculating µ 0 for
each component of the mixture and combining these values with the cross section data given in
the problem,

Σ tr = 0.6929 cm -1

giving

D = 0.481 cm.

Repeating the calculations as above, this value of D gives

Bm2 = 0.1488 cm -2

and the extrapolated critical sphere radius is 8.144 cm. The mass of 235U in this sphere is 113.1 g.
While this value is better than the first one calculated, it is still too far from reality to be of much
use. This example points out a weakness in applying one-group diffusion theory to real problems.
The criticality safety specialist must know the limitations of any method used for safety
calculations.

NCSET Module 6 Diffusion Theory 18 of 18