PHYS-E0562 Ydinenergiatekniikan jatkokurssi

Lecture 2 – Deterministic transport theory

Jaakko Leppänen (Lecturer), Ville Valtavirta (Assistant)
Department of Applied Physics
Aalto University, School of Science
Jaakko.Leppanen@vtt.fi

Mar. 2, 2016
Topics of this lecture
Introduction to the coupled problem:
I Non-linearity caused by coupling to thermal hydraulics and fuel depletion
I Scale of the coupled problem
I Reactor physics calculation chain

Transport theory:
I Neutron density, flux and current
I Formulation of the neutron transport equation

Physical properties of the transport problem:

I Fundamental and transient mode solutions
I Eigenvalue problem and criticality calculations
I Point kinetics approximation

Lecture 2: Deterministic transport theory

Mar. 2, 2016
2/74
What is expected to be known before the lecture
Reactor physics:
I Basic concepts: neutron flux, cross section, reaction rates
I Prompt and delayed neutrons

Misc:
I Basic stuff about differential equations
I Basic idea on how to solve neutron transport problems

Lecture 2: Deterministic transport theory

Mar. 2, 2016
3/74
Introduction to the coupled problem
The operation of a nuclear reactor is based on the self-sustaining chain reaction, carried on by
fission neutrons. Consequently, a major part of reactor analysis involves the transport and inter-
action physics of neutrons.

Solution to the neutron transport problem provides sufficient information on neutron-induced re-
action rate distributions within the reactor core, which can be used for calculating output, power
distribution, depletion and production rates of nuclides, etc.

There are major challenges in the solution of the transport problem:

I High level of heterogeneity, as neutrons are not uniformly distributed in the fuel
I Reaction probabilities (microscopic cross sections) are strongly dependent on neutron
energy
I Time-scales of prompt and delayed neutrons differ by several orders of magnitude ([µs] or
[ms] vs. [s])

Reactor analysis relies heavily on computational modeling, and any transport method used for
solving reactor physics problems needs to be able to cope with these challenges.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
4/74
Introduction to the coupled problem

Figure 1 : Neutron density distribution in reactor geometry. Left: BWR fuel assembly. Right: Small test
reactor core (ATR reactor, Idaho National Laboratory). When frozen in time, most of the neutrons in
thermal systems are found in large moderator regions. Density is considerably lower in fuel and
absorbers.2

2
The neutron density distribution is not to be confused with flux distribution, which also depends on neutron
speed, therefore emphasizing higher energies (see Fig. 2).

Lecture 2: Deterministic transport theory

Mar. 2, 2016
5/74
Introduction to the coupled problem
Even though microscopic cross sections are complicated functions of neutron energy, the trans-
port problem is linear by nature with the following assumptions:
I Neutrons do not interact with each other3
I Reaction probabilities are independent of reaction rates

The linearity assumption applies to criticality experiments and very low power research reactors.
For commercial power reactors, the reality is unfortunately different, because reaction probabilities
depend on:

1) Material temperatures and densities

2) Material compositions

This creates a non-linear problem, in which the neutronics solution is coupled to:

1) Heat transfer and coolant flow (via reactivity feedbacks)

2) Isotopic changes in fuel during the reactor operating cycle (via fuel depletion)

Solving the transport problem alone is not sufficient for the modeling of an operating nuclear
reactor – the calculation scheme needs to be able to handle also reactivity feedbacks and fuel
depletion
3 10 3
Neutron density in a power reactor is in the order of 10 1/cm , which can be compared to the density of
22 3
hydrogen atoms in water or uranium atoms in fuel (∼10 1/cm ).

Lecture 2: Deterministic transport theory

Mar. 2, 2016
6/74
Introduction to the coupled problem
The coupling between neutronics and thermal hydraulics is particularly strong in light water reac-
tors, because the reactor operates on thermal neutrons (∼80% of fissions), and the slowing-down
process is sensitive to operating conditions:
I Neutrons are born with a mean energy of 2 MeV, and they require ∼20 elastic collisions
with hydrogen nuclei to reach the thermal region (< 1 eV)
I If the neutron returns to fuel too early, it is more likely to be absorbed in the capture
resonances of 238 U than cause fission in 235 U
I Since coolant serves the purpose of moderator as well, any change in the reactor
operating state is immediately reflected in the slowing-down process

Reactivity feedbacks represent fast, almost immediate coupling between neutronics and reactor
operating conditions. In the long time-scale, similar non-linearity is induced by changes in fuel
composition by neutron irradiation:
235
I U is depleted and replaced by 239 Pu as the primary fissile isotope
I Non-fissile plutonium, minor actinides (Np, Am, Cm) and fission products are accumulated
in the fuel, increasing absorption
I Burnable absorber used for passive reactivity control is depleted

Nuclear fuel in light water reactors is loaded in the reactor core for the duration of the entire
operating cycle, which is typically 12 or 18 months. A single fuel assembly remains in the core for
several cycles, and is exposed to intense neutron irradiation for 3 to 4 years.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
7/74
Introduction to the coupled problem
In practice, the real challenge in reactor modeling is to obtain the solution to the coupled problem
by assuming that the problem can be linearized over some time interval and performing iterations
between the different solvers.

Obtaining high-fidelity solutions to neutronics, heat transfer, coolant flow and fuel depletion be-
comes a tremendous task:
I Typical LWR core contains 50,000 - 100,000 fuel rods, number of fuel pellets is counted in
millions
I Fuel temperature varies by hundreds of degrees, depending on local power4
I Significant variation in coolant density, especially in BWR’s
I Accurate simulation of fuel burnup would require tracking the concentrations of hundreds
of isotopes in millions of depletion zones

Calculations involving a single state-point without burnup are barely within the capabilities of mod-
ern super-computers.

Reactor design and safety analyses require covering various operating states, time-dependent
simulations, modeling of core depletion over multiple cycles, modeling of reactivity control, etc.,
which renders the direct approach unfeasible in practice.

4
Even the difference between surface and center-line temperature in a fuel pellet can exceed 400 K during
normal operation.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
8/74
Reactor physics calculation chain
Instead, the traditional approach to coupled problems is to divide the neutronics solution in parts:

Cross
I The scale of the modeled system is gradually
Cross section
section measurements
measurements
and
and nuclear
nuclear models
models increased, while simultaneously moving towards more
simplified description of physics
I Spatial resolution: nuclide-level → pin-level →
Isotopic
Isotopic micro-group
micro-group assembly-level → core-level
cross
cross sections
sections
I Energy resolution: continuous-energy → micro-group
structure → macro-group structure → few-group
Transport
Transport calculation
calculation at
at structure
fuel
fuel assembly
assembly level
level
I Transition from one stage to the next is carried out in
such way that local reaction rate balance is preserved

Homogenized
Homogenized few-group
few-group
I The details of the calculation chain depend on the
constants
constants modeled system and the methods used

The final stage involves a simplified description of the full-scale

Coupled
Coupled full-core
full-core system, which can be coupled to the solution of heat transfer
simulation
simulation
and coolant flow.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
9/74
Reactor physics calculation chain
From an engineering point of view, the typical LWR calculation chain can be divided into three
computational tasks:

(i) Production of nuclide-wise microscopic interaction data for neutron transport codes from
evaluated nuclear data
(ii) Solution of the local heterogeneous transport problem at fuel assembly level using
higher-order transport methods
(iii) Solution of the global homogeneous transport problem using diffusion theory

The second step is called spatial homogenization, and it provides the sufficient building blocks for
the full-scale simulation:
I Traditionally carried out using deterministic 2D lattice transport codes
I Input: isotopic micro-group reaction cross sections, detailed description of the geometry
at fuel assembly level
I Output: handful of macroscopic few-group constants representing the transport physics

Lecture 2: Deterministic transport theory

Mar. 2, 2016
10/74
Reactor physics calculation chain
Dividing the calculation chain into separate steps enables running coupled full-core simulations at
an acceptable computational cost without compromising the accuracy.

The cost of reducing the physical complexity by spatial homogenization is that all information on
reactivity feedback effects and dependence of reaction rates on fuel burnup is completely lost.

This information is recovered by repeating the procedure over and over again in such way that all
reactor operating conditions are covered:
I Assembly types
I Variation in thermal-hydraulic state variables: fuel temperature, coolant temperature and
density
I Reactivity control: soluble boron, insertion of control rods
I Fuel depletion

The result is a parametrized library of reactor-specific group constants, from which the data cor-
responding to local operating conditions can be obtained by interpolation.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
11/74
Reactor physics calculation chain
The neutronics solution in LWR core calculations is typically based on two-group nodal diffusion
methods, and coupled simulations can be used, for example, for:
I Steady state fuel cycle simulations, modeling normal reactor operation over one or
several cycles
I Dynamic simulations, modeling the behavior of the reactor core in transients

Fuel cycle simulations are needed for core design, to make sure that the reactor remains critical
and within the safety margins. Dynamic simulations are needed for safety analyses, and to study
the behavior in abnormal operating conditions.5

The multi-stage calculation chain based on spatial homogenization and nodal diffusion methods
has the potential to produce very accurate results in full-scale 3D simulations at an acceptable
computational cost ...

5
Nuclear power plant simulations require modeling the reactor core as part of the coolant loop(s), including
pumps, valves, heat exchangers, etc. This type of simulations are carried out using system codes, in which the
reactor is modeled using similar methods as in core simulations, or using simplified methods, such as 1D neutronics
or point-kinetics approximation.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
12/74
Reactor physics calculation chain
... but the approach also has its limitations:
I Homogenized group constants are case-specific, i.e. they must be generated separately
for each problem
I Different reactor types (LWR, SFR, HTGR) may require different methods both in
assembly and core-level calculations
I Traditional deterministic lattice transport codes usually rely on a 2D solution, which is not
sufficient for describing 3D effects encountered in some advanced reactor types
I Dividing the full-scale transport problem into isolated sub-problems creates artificial
interfaces, which affect the quatlity of the solution

In general, there is a multitude of deterministic transport codes and methods developed for dif-
ferent purposes, and no universal calculation scheme that could handle all applications with the
same level of reliability. The two-group nodal diffusion methods covered in this course apply
mainly to conventional LWR analyses.

The results are always subject to uncertainties and errors in both methodology and input data,
and the calculation scheme can be validated by comparison to experimental data and high-fidelity
methods.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
13/74
Reactor physics calculation chain
The remaining lectures in this course cover the various stages of the multi-stage calculation
scheme and method used for obtaining the solution to the coupled problem:
I This and the following lecture cover the deterministic and stochastic (Monte Carlo)
approach to transport theory
I Diffusion theory, which forms the basis of nodal diffusion methods, is derived from
transport theory in Lecture 4
I Lecture 5 deals with fuel depletion and the solution of the burnup problem
I Lecture 6 covers the basics of heat transfer and coolant flow
I Lectures 7 and 8 cover spatial homogenization and nodal diffusion methods, i.e. the
solution of the full-scale transport problem
I Lectures 9 and 10 cover the methods used thermal hydraulics and fuel behaviour
modeling
I Lecture 11 covers reactor dynamics (time-dependent solution of the coupled problem)
I Lecture 12 covers design and simulation of reactor operating cycle

Lecture 2: Deterministic transport theory

Mar. 2, 2016
14/74
Basic concepts of transport theory
The main objective of neutron transport calculation is to solve the nuclear reaction rate distribution
within the reactor core. All solution methods rely on transport theory, which is based on relatively
simple conservation laws and a few mathematical definitions.

Transport theory of neutrons shares some similarities with other transport problems, for example,
in fluid dynamics and plasma physics, with a few significant differences:

1) The interaction probabilities are very strongly dependent on neutron energy

2) The forces acting in collisions have very limited range, and the neutron path is broken into
a “random walk” from one collision to the next
3) Neutron-neutron interactions can be ignored

As discussed later on, neutron transport can be considered a linear problem under certain ap-
proximations.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
15/74
Basic concepts of transport theory: cross sections
The probabilities of neutron-induced reactions are characterized by cross sections. As discussed
in Lecture 1, the total macroscopic cross section, Σ, gives the total interaction probability per
traveled path length:
dP = Σ(r, E)ds (1)

Scattering reactions are associated with a change in neutron energy and direction of motion,
and they are described by double-differential cross sections, which combine the probability of
the scattering event to distributions of energy transfer and scattering angle. The probability of a
0
neutron scattering from direction Ω̂ to Ω̂ and energy E to E 0 per traveled path length is given
by:
0 0 0 0
dP = Σs (r, Ω̂ → Ω̂ , E → E )dsdΩ dE (2)

The change in direction is often written using the

0
scattering cosine µ = cos θ = Ω̂ · Ω̂ , i.e. the cosine of
′
dΩ̂ ϕ
dΩ̂ θ
the angle θ between the direction vectors before and after
the collision (the second angle, ϕ, describing the rotation
of the scattered direction vector around the incident
direction is called the azimuthal angle).

The double-differential scattering cross section is then written as Σs (r, µ, E → E 0 ). This is a

valid approximation in an isotropic medium, in which case the scattering rate does not depend on
the absolute direction of the incident neutron.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
16/74
Basic concepts of transport theory: cross sections
Macroscopic cross section is related to microscopic cross section, σ , by:

Σ(r, E) = N (r)σ(E) (3)

where N is the nuclide density. Microscopic cross section can therefore be interpreted as the
interaction probability per traveled path length and nuclide density. An alternative interpretation is
that this parameter represents the effective cross-sectional area of the nucleus.

Microscopic cross sections are nuclide-specific constants, which depend on the reaction type and
neutron energy. The total cross section can be written as the sum over partial reaction modes:

σ(E) = σcapt (E) + σfiss (E) + σela (E) + . . . (4)

The same applies to macroscopic cross sections:

Σ(r, E) = Σcapt (r, E) + Σfiss (r, E) + Σela (r, E) + . . .

  (5)
= N (r) σcapt (E) + σfiss (E) + σela (E) + . . .

And if the medium consists of multiple nuclides:

X X
Σ(r, E) = Σi (r, E) = Ni (r)σi (E) (6)
i i

where the summation is carried over all the constituent nuclides.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
17/74
Basic concepts of transport theory: neutron population
All deterministic transport methods treat neutrons as a collective population of independent parti-
cles, with the help of three concepts:
I Neutron density
I Neutron flux
I Neutron current (or current density)

Neutron density and flux are best understood as a two closely related density-like functions in the
six-dimensional phase-space, depending on the position and momentum variables, although the
momentum variable is usually replaced by energy and direction of motion.
The direction of motion Ω̂ depends on two angular variables,
z, Ωz
and it can be written using the three Cartesian direction
vectors as:
Ω̂
Ωx = sin η cos ϑ η
Ωy = sin η sin ϑ (7)
Ωz = cos η ϑ x, Ωx

where η ∈ [0, π] is the polar and ϑ ∈ [0, 2π] is the azimuthal y, Ωy

angle.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
18/74
Basic concepts of transport theory: neutron density
Angular neutron density, n(r, Ω̂, E), gives the number of neutrons inside an infinitesimal volume
element dV at position r:
n(r, Ω̂, E)dV dΩ̂dE (8)
The infinitesimal element in the angular space is the differential solid angle dΩ̂ about the direction
of motion Ω̂. The energy is distributed within interval dE about E .
z

dΩ̂
dz
dy
dx
x

The density of neutrons moving in all directions, also called the scalar density, is given by:
Z
n(r, E) = n(r, Ω̂, E)dΩ̂ (9)
4π
and the density of neutrons moving in all directions at all energies:
Z Z
n(r) = n(r, Ω̂, E)dΩ̂dE (10)
4π E

Lecture 2: Deterministic transport theory

Mar. 2, 2016
19/74
Basic concepts of transport theory: neutron flux
Neutron flux, or more precisely, angular neutron flux, has no straightforward physical interpreta-
tion, but it is defined as the angular neutron density multiplied by speed:
ψ(r, Ω̂, E) = vn(r, Ω̂, E) (11)

The difference between the two concepts can be understood by considering a randomly dis-
tributed population of neutrons passing through a volume at different energies:
I If a snapshot, frozen in time, is taken of the population, it shows the distribution of
neutrons in different positions, energies and directions of motion – the angular neutron
density in the six-dimensional phase space.
I If another snapshot is taken a short time interval dt later, each neutron has moved forward
by distance ds = vdt. The faster the neutrons are moving, the longer the combined
distance covered by the entire population. This distance will be later related to the integral
of the neutron flux and the number of neutron-induced reactions within the volume.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
20/74
Basic concepts of transport theory: neutron flux

Figure 2 : Neutron density and flux distribution in a BWR fuel assembly integrated over energy and
direction. Left: Density distribution peaks in the internal and external moderator channels where thermal
neutrons are collected. Right: Flux distribution is clearly more uniform, because the neutron speed
multiplier (ψ = vn) emphasizes the contribution of high-energy neutrons with longer mean-free-paths
and more uniform distribution over the geometry.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
21/74
Basic concepts of transport theory: neutron flux
The neutron flux is not a physical, measurable quantity,6 and its practical significance actually
results from a simple connection to reaction rates. The rate of reaction x in the infinitesimal phase-
space element can be written as the product of angular flux and the corresponding macroscopic
cross section:
dRx = Σx (r, E)ψ(r, Ω̂, E)dV dΩ̂dE (12)

When integrated over variables, the result is a physical reaction rate:

Z Z Z
Rx = Σx (r, E)ψ(r, Ω̂, E)dV dΩ̂dE (13)
V Ω̂ E

that gives the number of reactions x per second, induced by neutrons within the domain of inte-
gration. For example, the total capture rate in volume V is given by:
Z Z Z ∞
Σγ (r, E)ψ(r, Ω̂, E)dV dΩ̂dE (14)
V 4π 0

and the total number of fission neutrons produced by thermal neutron-induced fission:
Z Z Z E
1
νΣf (r, E)ψ(r, Ω̂, E)dV dΩ̂dE (15)
V 4π 0

where E1 is the boundary between thermal and fast energy group.

6
This interpretation is subject to argument, but what cannot be argued is that flux can only be measured via
physical reaction rates.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
22/74
Basic concepts of transport theory: neutron flux
For scattering reactions the product of angular flux and the double-differential scattering cross
section describes the differential reaction rate density, i.e. the rate at which neutrons are scattering
0
from direction Ω̂ to Ω̂ and energy E to E 0 :7
0 0 0 0
dRs = Σs (r, Ω̂ → Ω̂ , E → E )ψ(r, Ω̂, E)dV dΩ̂dEdΩ̂ dE (16)

Or if the change in direction is written using the scattering cosine:

0 0
dRs = Σs (r, µ, E → E )ψ(r, Ω̂, E)dV dΩ̂dEdµdE (17)

Similar to total reaction rate, the differential reaction rate can be integrated over variables, scat-
tering angles and emission energies. The total scattering rate is given by:
Z Z Z ∞Z 1 Z ∞
0 0
Rs = Σs (r, µ, E → E )ψ(r, Ω̂, E)dV dΩ̂dEdµdE (18)
V 4π 0 −1 0

Or similarly for the down-scattering rate, i.e. scattering from fast to thermal group:
Z Z Z ∞ Z 1 Z E1
0 0
Rs12 = Σs (r, µ, E → E )ψ(r, Ω̂, E)dV dΩ̂dEdµdE (19)
V 4π E1 −1 0

where E1 is the group boundary.

7
To be precise, dΩ̂ about Ω̂, etc.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
23/74
Basic concepts of transport theory: neutron flux
The neutron direction of motion plays no role in capture and fission events,8 which means that the
reaction rate can be integrated over direction of motion:
Z Z 
dRx = Σx (r, E)ψ(r, Ω̂, E)dV dΩ̂dE = Σx (r, E) ψ(r, Ω̂, E)dΩ̂ dV dE
4π 4π (20)
= Σx (r, E)φ(r, E)dV dE

where the scalar neutron flux is defined similar to the scalar neutron density (9):
Z
φ(r, E) = ψ(r, Ω̂, E)dΩ̂ = vn(r, E) (21)
4π

The scalar flux can be interpreted as the distance covered by neutrons with energy dE about E
inside an infinitesimal volume element dV at r in time interval dt:

ds = φ(r, E)dV dEdt (22)

When integrated over volume, energy and time, the result is the total combined distance cov-
ered by the neutrons, as suggested by the previous example illustrating the differences between
neutron density and flux.

8
The incident neutron is lost in capture and fission neutrons are typically assumed to be emitted isotropically,
which is a very good approximation in the energy range of fission reactors.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
24/74
Basic concepts of transport theory: neutron flux
The connection to reaction rates is seen by recalling that the macroscopic cross section Σ is
defined as the neutron interaction probability per traveled path length.

This connection reflects the fact that the rate of neutron-induced reactions is increased if:

1) Neutron density is increased, increasing the number of neutrons that may collide with the
nuclides in the medium
2) Neutron velocity is increased, increasing the path length traveled within the given time
interval and the number of chances of colliding with the nuclides

The same connection applies in the opposite case, when the neutron flux (density, velocity or
both) is decreased.

The reaction rate also naturally depends on the macroscopic cross section:

Σ(r, E) = N (r)σ(E) (23)

which in turn depends on:

1) The nuclide density N , which affects the number of nuclides encountered by the neutron
within its path
2) The microscopic cross section σ , representing the interaction probability between the
neutron and a single target nucleus

Lecture 2: Deterministic transport theory

Mar. 2, 2016
25/74
Basic concepts of transport theory: neutron flux

Figure 3 : Neutron flux is connected to reaction rates by multiplication with macroscopic cross sections
and integration over variables. Reaction rates therefore also depend on flux spectrum. This emphasizes
the contribution of low-energy neutrons, as cross sections increase towards lower energies. Left: Total
collision rate distribution in the previous BWR fuel assembly. Right: Total fission rate (“hot” shades) and
thermal flux (“cold” shades). Dark and bright colors indicate low and high values, respectively. Since more
than 80% of fissions are induced by thermal neutrons, the fission rate distribution follows the distribution of
thermal flux.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
26/74
Basic concepts of transport theory: neutron flux
Angular and scalar neutron flux are both density-like functions describing the neutron population,
and their differences become apparent in how they are primarily used:
I Angular flux, ψ(r, Ω̂, E), is used for determining the rate at which neutrons are moving
from one phase-space position to another by streaming or scattering
I Scalar flux, φ(r, E), is used for determining the rate of reactions where the incident
neutron is lost and secondary neutrons, if emitted, have no preferential direction.

It is important to note that, even though (angular) neutron density and flux depend on the direction
of motion, they are both scalar functions. The direction is a variable, similar to position and energy,
and it provides two of the six phase-space coordinates.9

9
Since human senses are limited to three geometric dimensions, it is not possible to visualize angular or scalar
flux accurately, but if you must, think about density (scalar field) instead of flow (vector field).

Lecture 2: Deterministic transport theory

Mar. 2, 2016
27/74
Basic concepts of transport theory: neutron current
The vector equivalent of angular flux is the angular current density:10

j(r, Ω̂, E) = vn(r, Ω̂, E)

(24)
= Ω̂ψ(r, Ω̂, E)

where v = Ω̂v is the neutron velocity. The angular current density describes the rate at which
neutrons with energy dE about E and direction dΩ̂ about Ω̂ pass through an infinitesimal surface
element dS located at r:
dL = j(r, Ω̂, E) · dS (25)

When the angular current density is integrated over the full solid angle, the result is another vector
quantity called neutron current density:
Z
J(r, E) = j(r, Ω̂, E)dΩ̂ (26)
4π

which gives the net rate of neutrons with energy dE about E passing through an infinitesimal
surface element dS located at r:
dL = J(r, E) · dS (27)

10
The terminology used in the literature for different flux and current quantities is not unambiguous. It should also
be noted that the concept of flux in neutron transport theory is not equivalent with the mathematical definition of flux,
used correctly, for example, in electromagnetic field theory.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
28/74
Basic concepts of transport theory: neutron current

Figure 4 : Illustration of angular and neutron current densities as vector fields in a BWR fuel assembly.
Left: Angular current density of thermal neutrons traveling in a 30◦ direction with respect to the positive
x-axis. Right: Neutron current density, representing the net flow of thermal neutrons. Distribution of
thermal neutron density is plotted in the background.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
29/74
Basic concepts of transport theory: neutron current
One of the crucial differences between the flux and the current quantities is seen by comparing
Eqs. (21) and (26). If the neutron density is completely isotropic (all directions equally likely), then:
Z
J(r, E) = j(r, Ω̂, E)dΩ̂ = 0
4π
Z (28)
φ(r, E) = ψ(r, Ω̂, E)dΩ̂ > 0
4π

Whereas the neutron flux can be associated with the total distance covered by neutrons in a
volume and the rate of neutron-induced reactions, the current density can be associated with the
flow of neutrons over boundary surfaces.

Neutron current, i.e. the rate at which neutrons are crossing surface S is given by surface integral:
Z Z Z Z Z Z  
J = j(r, Ω̂, E) · dSdΩ̂dE = j(r, Ω̂, E) · û dΩ̂dSdE (29)
Ω̂ S E Ω̂ S E

where û is the surface normal.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
30/74
Basic concepts of transport theory: neutron current
Similar to integral reaction rates, neutron current can also be divided into components with respect
to space, energy and direction. If surface S defines the boundary of a volume, the rate at which
neutrons are streaming in is given by the inward current
Z 0 Z Z  
−
J = j(r, Ω̂, E) · û dµdSdE (30)
−1 S E

and the rate at which neutrons are streaming out by the outward current:
Z 1Z Z  
+
J = j(r, Ω̂, E) · û dµdSdE (31)
0 S E

where µ = Ω̂ · û is the cosine between the neutron direction of motion and the surface normal.11
The values of inward and outward currents are negative and positive, respectively.

Integral of the neutron current density yields the net current:

Z Z  
Jnet = J(r, E) · û dSdE (32)
S E

The value is positive for outward and negative for inward flow. The net current is related to the
partial currents (30) and (31) by:
− +
Jnet = J + J (33)

11
Assuming that û points in the outward direction.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
31/74
Basic concepts of transport theory: summary
Angular neutron flux, ψ , can be thought of as a mathematical density-like function in the six-
dimensional phase space that connects the behavior of the neutron population to physical reaction
rates: Z Z Z
Rx = Σx (r, Ω̂, E)ψ(r, Ω̂, E)dV dΩ̂dE (34)
V Ω̂ E

The directional dependence is important in scattering reactions, described by double-differential

cross sections, which combine the probability of the scattering event to the distributions of energy
transfer and scattering angle.

In reactions where the incident neutron is lost or secondary neutrons are emitted isotropically,
angular flux can be replaced by the scalar flux, which is obtained by integration over full solid
angle: Z
φ(r, E) = ψ(r, Ω̂, E)dΩ̂ (35)
4π

The integral reaction rate is then reduced to:

Z Z
Rx = Σx (r, E)φ(r, E)dV dE (36)
V E

The integration of flux over the phase space variables and time yields the total combined distance
traveled by the neutrons. This is intuitively associated with reaction rate by recalling that the
macroscopic cross section is defined as the interaction probability per path length traveled by the
neutron.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
32/74
Basic concepts of transport theory: summary
The vector equivalent of angular flux is the angular current density:

j(r, Ω̂, E) = Ω̂ψ(r, Ω̂, E) (37)

which integrated over full solid angle yields the neutron current density:
Z
J(r, E) = j(r, Ω̂, E)dΩ̂ (38)
4π

It should be noted that, even though angular flux and angular current density are related by (37),
similar connection does not apply for scalar flux and neutron current density:

J(r, E) 6= Ω̂φ(r, E) (39)

The current densities are associated with the rates at which neutrons cross the boundaries of a
specified volume. When (37) is integrated over a surface, direction, energy and time, the result is
the net number of neutrons that have passed through the surface.

Reaction rate integrals appear in the source and removal terms of the neutron transport equation.
The streaming term can be written using the neutron current density, and in diffusion theory it is
associated to flux gradient by an approximation known as Fick’s law.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
33/74
Neutron transport equation
General neutron transport theory forms the basis of all deterministic transport calculation meth-
ods, including diffusion theory, which is used in full-scale reactor simulator calculations. The
theory is based on the following fundamental assumptions:

1) Neutrons can be treated as independent particles, traveling in straight lines between

collisions
2) Neutron-neutron interactions can be ignored
3) Reaction cross sections are independent of flux and constant in time

With these assumptions the neutron transport equation describes a linear problem.12 Additional
approximations include:

4) All materials are isotropic, i.e. differential scattering cross sections depend only on
scattering angles and not on neutron direction of motion
5) Fission neutrons are emitted isotropically
6) Fission spectrum is independent of incident neutron energy

Approximations 5 and 6 are valid at energies relevant for fission reactor applications.

12
As pointed out earlier, assumption 3) does not hold in coupled problems describing an operating nuclear reactor
subject to reactivity feedbacks and fuel burnup, in which case the transport problem actually becomes non-linear.
The solution to the coupled problem is obtained by iteration between the different solvers, assuming that the
transport problem can be linearized over sufficiently short time intervals.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
34/74
Neutron transport equation
The transport equation is formulated simply by considering the particle density balance inside an
infinitesimal element of the six-dimensional phase space:13
I Position: dV about r
I Direction: dΩ̂ about Ω̂
I Energy: dE about E

The equation can be written as:

1 ∂
ψ(r, Ω̂, E, t) + Ω̂ · ∇ψ(r, Ω̂, E, t) + Σ(r, E)ψ(r, Ω̂, E, t) = q(r, Ω̂, E, t) (40)
v ∂t
| {z } | {z } | {z }
(A) (B) (C)

where q is the general source term and:

(A) is the time-rate of change in neutron density inside the phase-space element14
(B) is the streaming term, i.e. the net rate at which neutrons moving in direction Ω̂ with
energy E are entering and leaving the volume element (see next slide)
(C) is the removal term, i.e. the rate at which neutrons are removed from the phase-space
element by absorption or scattering to new energy and direction
13
See figure on slide 19.
14
1 ∂ψ ∂n
Results from the relation between flux and neutron density: ψ = vn =⇒ =
v ∂t ∂t
Lecture 2: Deterministic transport theory
Mar. 2, 2016
35/74
Neutron transport equation
The form of terms (A) and (C) is easily understood, but the streaming term is less straightforward.
It was stated earlier that the rate at which neutrons with energy E and direction Ω̂ pass through
an infinitesimal surface element dS located at r is given by the angular current density:

dL = j(r, Ω̂, E, t) · dS (41)

If S is the complete boundary of the infinitesimal volume element dV , then the net rate at which
neutrons are flowing in or out is given by surface integral:
I
j(r, Ω̂, E, t) · dS (42)
S

Since angular current density is a vector field, divergence theorem yields:

I Z
j(r, Ω̂, E, t) · dS = ∇ · j(r, Ω̂, E, t)dV (43)
S V

Angular current density can be written using angular flux j = Ω̂ψ , which results in:
Z Z
∇ · Ω̂ψ(r, Ω̂, E, t)dV = Ω̂ · ∇ψ(r, Ω̂, E, t)dV (44)
V V

showing the final form of the streaming term integrated over volume.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
36/74
Neutron transport equation
The source term can be written as the sum of external, scattering and fission source:

q(r, Ω̂, E, t) = Q(r, Ω̂, E, t) + S(r, Ω̂, E, t) + F (r, Ω̂, E, t) (45)

The external source term, Q, is independent of flux and it gives the rate at which neutrons are
emitted into the phase-space element by external sources.

The scattering source can be written using the double-differential scattering cross section, by
integration over all incident energies and directions:15
Z Z ∞
0 0 0 0 0 0
S(r, Ω̂, E, t) = Σs (r, Ω̂ → Ω̂, E → E)ψ(r, Ω̂ , E , t)dΩ̂ dE (46)
4π 0

Since fission neutrons are emitted isotropically, the fission source term can be written using the
scalar flux: Z ∞
χ (E) 0 0 0
F (r, Ω̂, E, t) = νΣf (r, E )φ(r, E , t)dE (47)
4π 0

where νΣf is the fission neutron production cross section and χ(E) is the fission spectrum, i.e.
the probability that the energy of the emitted neutron falls on interval dE about E . Factor 1/4π
results from the source isotropy.
15
It is assumed here that the double-differential scattering cross section is written as the sum over all nuclides and
scattering modes. Neutron-multiplying (n,2n), (n,3n), etc. reactions are accounted for in the distribution of
secondary energies and scattering angles (normalization to 2 or 3 instead of 1).

Lecture 2: Deterministic transport theory

Mar. 2, 2016
37/74
Neutron transport equation
The transport equation describes accurately the balance in angular neutron density with the help
of angular flux. Similar equation can be derived for the scalar density and flux, by integrating (40)
over the full solid angle.

The integration of the time-derivative term and the removal term is straightforward:
Z Z
1 ∂ 1 ∂ 1 ∂
ψ(r, Ω̂, E, t)dΩ̂ = ψ(r, Ω̂, E, t)dΩ̂ = φ(r, E, t) (48)
4π v ∂t v ∂t 4π v ∂t

and
Z Z
Σ(r, E)ψ(r, Ω̂, E, t)dΩ̂ = Σ(r, E) ψ(r, Ω̂, E, t)dΩ̂ = Σ(r, E)φ(r, E, t) (49)
4π 4π

The fission source term (47) was already written using the scalar flux and integration over full solid
angle cancels the 1/4π factor:
Z
χ (E)
Z ∞ Z ∞
0 0 0 0 0 0
νΣf (r, E )φ(r, E , t)dE dΩ̂ = χ (E) νΣf (r, E )φ(r, E , t)dE (50)
4π 4π 0 0

The external source is independent of the flux and its angular dependence depends on the source
type. The integration of the two remaining terms is less straightforward.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
38/74
Neutron transport equation
The integration of the scattering source (46) is written as:
Z "Z Z ∞
#
0 0 0 0 0 0
Σs (r, Ω̂ → Ω̂, E → E)ψ(r, Ω̂ , E , t)dΩ̂ dE dΩ̂ (51)
4π 4π 0

Since the integration is carried over the full solid angle, and because the double-differential scat-
0
tering cross section depends only on the angle between directions Ω̂ and Ω̂, the integral over
the double-differential scattering cross section can be written using the scattering cosine:
"Z #
Z Z ∞ 1
0 0 0 0 0
Σs (r, µ, E → E)dµ ψ(r, Ω̂ , E , t)dΩ̂ dE (52)
4π 0 −1

The term in brackets describes the total scattering probability from energy E 0 to E , and can be
written as: Z Z ∞
0 0 0 0 0
Σs (r, E → E)ψ(r, Ω̂ , E , t)dΩ̂ dE
4π 0
Z ∞
"Z # (53)
0 0 0 0 0
= Σs (r, E → E) ψ(r, Ω̂ , E , t)dΩ̂ dE
0 4π

The term in brackets is the scalar flux, and the final form can be written as:
Z ∞
0 0 0
Σs (r, E → E)φ(r, E , t)dE (54)
0

Lecture 2: Deterministic transport theory

Mar. 2, 2016
39/74
Neutron transport equation
The final term in (40) left to be integrated is the streaming term, which can also be written using
the angular current density (24):

Ω̂ · ∇ψ(r, Ω̂, E) = ∇ · Ω̂ψ(r, Ω̂, E) = ∇ · j(r, Ω̂, E) (55)

Integration over the full solid angle yields

Z Z
∇ · j(r, Ω̂, E)dΩ̂ = ∇ · j(r, Ω̂, E)dΩ̂ = ∇ · J(r, E) (56)
4π 4π

where J(r, E) is the neutron current density (26).

When the results of Eqs. (48)–(56) are collected and external source term is omitted for conve-
nience, the result can be written as:
1 ∂
φ(r, E, t) + ∇ · J(r, E) + Σ(r, E)φ(r, E, t)
v ∂t
Z ∞
" # (57)
0 0 0 0 0
= Σs (r, E → E)φ(r, E , t) + χ (E)νΣf (r, E )φ(r, E , t) dE
0

This form of transport equation is also known as the neutron continuity equation, and it forms the
starting point for the derivation of diffusion theory, discussed in Lecture 4. It should be noted that
no approximations were made in the derivation of (57) from the original transport equation (40).

Lecture 2: Deterministic transport theory

Mar. 2, 2016
40/74
Neutron transport equation
Even though the transport equation is linear, it is extremely difficult to solve because:

1) The geometries are complex and heterogeneous

2) Cross sections are complicated functions of neutron energy
3) The angular dependence of streaming and scattering source term are difficult to handle

In practice, the integrals over space, direction and energy cannot be resolved while holding on to
the continuous dependence on phase-space variables.

For practical reasons, most deterministic transport methods therefore rely on at least three ap-
proximations:

1) The geometry is discretized into a number of homogeneous material regions

2) The continuous energy dependence of cross sections is condensed into a number of
discrete energy groups
3) The angular dependence of double-differential scattering cross sections is represented by
functional expansions, the directional dependence of flux is represented by functional
expansions or discrete directions

The first two approximations are common to all deterministic transport methods, the treatment
of angular dependence is what differentiates the methods from each other (Sn , Pn , method of
characteristics, diffusion theory, etc.).

Lecture 2: Deterministic transport theory

Mar. 2, 2016
41/74
Neutron transport equation
The solution by applying diffusion approximation is left for Lecture 4, and the remainder of this
lecture is devoted to looking at the different physical aspects of the transport problem from the
viewpoint of deterministic transport theory.16

From here on, the transport equation is written for the sake of simplicity without the phase space
variables:
1 ∂ψ
+ Ω̂ · ∇ψ + Σψ = Q + F + S (58)
v ∂t
It is assumed that the accurate solution is available, even though this is clearly not the case in
reality.

Integration of (58) over variables preserves the neutron balance, and turns streaming term into
total leakage rate: Z Z Z I Z
Ω̂ · ∇ψ dV dΩ̂dE = J · dSdE (59)
V Ω̂ E S E

Integration of the removal term yields the total reaction rate, and the absorption rate is given by:
Z Z Z  
Σψ − S dV dΩ̂dE (60)
V Ω̂ E

Integration of fission and external source terms yields the corresponding total source rates.

16
In Lecture 3 the perspective is switched to that of the individual neutron with the introduction of Monte Carlo
neutron transport simulation.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
42/74
Equilibrium state and criticality
When no difference is made between prompt and delayed neutrons, the criticality state of the
reactor can be related to whether the chain reaction is below, at or above the self-sustaining
condition:17
I When the system is in sub-critical state, the chain reaction is not self-sustained, and the
fission chains die out after a number of generations
I When the system is in critical state, the chain reaction is self-sustained, and each
terminated fission chain is replaced, on the average, by exactly one new chain
I When the system is in super-critical state, the chain reaction is self-sustained, and each
terminated fission chain is replaced, on the average, by more than one new chain

In the absence of external sources, it is easy to see that the time-dependence of neutron popula-
tion takes the form of:
ωt
n(t) ∼ n0 e (61)

where ω depends on the source and loss terms. Steady-state condition (ω = 0) always implies
criticality, i.e. that the chain reaction is in a self-sustained state, and fission source rate equals the
sum of absorption and leakage rates.

The absolute population size is not fixed, which results from the fact that the system is character-
ized by a homogeneous differential equation.
17
In contrast to the description given in Lecture 1, the fission chains here are assumed to consist of both prompt
and delayed neutrons.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
43/74
Equilibrium state and criticality
In the presence of constant external source, time-dependence of neutron population takes the
form of:  
ωt Q ωt
n(t) ∼ n0 e − 1−e (62)
ω
where ω 6= 0. If the system is sub-critical (ω < 0), the asymptotic solution is constant in time:
Q
lim n(t) = − (63)
t→∞ ω
reflecting the fact that source neutrons are multiplied in fission chains that are finite in length.
In fact, a steady state solution in the presence of external source can only exist if the system is
sub-critical.
In the special case of criticality, the transport equation with external source is reduced to:

1 ∂ψ
=Q (64)
v ∂t
as the streaming, removal, and fission and scattering source terms are summed to zero. The
time-dependence of neutron population is then of form:

n(t) ∼ Qt (65)

which is also understood by recalling that in critical state each neutron replicates itself (on the
average) without additional multiplication. In other words, the external source initiates new non-
branching infinitely long fission chains at a constant rate.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
44/74
Equilibrium state and criticality

Sub-critical Sub-critical
Critical Critical
Super-critical Super-critical
Neutron population

Neutron population
Time Time

Figure 5 : Neutron population as function of time in sub-critical, critical and super-critical states. Left:
Equilibrium in the absence of external sources implies criticality. When the system is sub- or super-critical,
the population is exponentially decreasing or increasing, respectively. Right: Equilibrium in the presence
of external sources may only exist if the system is sub-critical (and the source is constant). The growth
rate is linear at criticality and exponential when the system is super-critical.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
45/74
Fundamental and transient flux modes
Since the transport problem is characterized by a linear equation, its solution can be written as a
linear combination of different flux modes:

ψ = ψ0 + ψ1 + ψ2 + . . . (66)

where ψ0 is called the fundamental mode and the remaining modes the transient modes. The
different flux modes ψn are associated with the eigenfunctions of the transport equation and time
constants ωn , such that:
ω0 > ω1 > ω2 > . . . (67)

Since the time dependence of flux is characterized by an exponential amplitude function, and:
ω0 t ω1 t ω2 t
e >e >e > ... (68)

it is easy to see that the fundamental mode is the asymptotic solution to the transport equation,
which persists after the transient modes, excited by changes in the physical conditions, fade away.

Time constants ωn can be associated to eigenvalues kn , in such way that:

k0 > k1 > k2 > . . . (69)

Eigenvalue k0 corresponding to the fundamental flux mode is interpreted as the multiplication

factor.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
46/74
Fundamental and transient flux modes

8 ψ0 ψ1 ψ2 ψ3 1.0

6
0.8
4
Neutron flux (A.U)

Eigenvalue
0.6
0

−2 0.4
−4

−6 0.2

−8
0.0
0 100 200 300 400 0 10 20 30 40 50
Axial coordinate (cm) Eigenmode

Figure 6 : Left: First four axial flux modes in a PWR full core calculation (amplitudes not to scale). The
solutions are asymmetric, because of axial profiling of burnable absorber. Total flux and fundamental flux
mode ψ0 are always positive, but transient modes may have negative components. Right: First 50
eigenvalues. As the plot shows, k0 > k1 > k2 > . . . . Each eigenvalue is associated with a time
constant, and the transient modes fade away as t → ∞.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
47/74
Fundamental and transient flux modes

Figure 7 : First 15 radial flux modes in the PWR core calculation. The fundamental flux mode (top left)
corresponds to the flux shape in steady-state condition. Calculated by Monte Carlo simulation using the
fission matrix method. See also Lecture2_anim1.gif for an animation on the convergence of flux into
fundamental mode in a fast criticality experiement (Pu-Flattop).

Lecture 2: Deterministic transport theory

Mar. 2, 2016
48/74
Criticality eigenvalue problem
Several applications in reactor physics involve the solution of the steady-state transport equation,
even though there is no balance between the physical source and loss rates. An example of such
task is criticality calculation, which essentially means determining how far the system is from the
self-sustaining state.18

In such case, the transport equation is usually written in the criticality- or k-eigenvalue form:

1
Ω̂ · ∇ψ + Σψ = F +S (70)
k
i.e. time-dependence and external source are dropped, and balance between source and loss
rates is attained by scaling the fission term with constant 1/k.

In theory the solution consists of an infinite number of eigenfunctions ψ0 , ψ1 , ψ2 , . . . with as-

sociated eigenvalues k0 , k1 , k2 , . . . , but in practice it is the fundamental flux mode ψ0 which is
taken as the solution to (70).

The fact that the higher flux modes can be dropped is based on the idea that the solution repre-
sents a steady-state system. As discussed earlier, this is the asymptotic solution that persists as
t → ∞ and the higher modes fade away.

18
As will be discussed in Lecture 7, this type of problems are also encountered in spatial homogenization, in which
the neutronics solution is obtained in some sub-domain of the full-scale geometry, and the accurate description of
boundary leakage currents is not available.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
49/74
Criticality eigenvalue problem
Eigenvalue k0 corresponding to the fundamental flux mode is also called the effective multiplica-
tion factor, keff , which is used to define the criticality condition as:

< 1 ⇒ system is sub-critical

keff = 1 ⇒ system is critical (71)

> 1 ⇒ system is super-critical


If the integral source and loss rates19 are known, keff can be associated to neutron balance:

F
keff = (72)
T −S+L

where

F is the integral fission source rate

T is the integral total reaction rate
S is the integral scattering source rate
L is the integral leakage rate
The integral scattering source rate includes neutron production from multiplying (n,2n), (n,3n), etc.
reactions.
19
integrals over space, direction and energy

Lecture 2: Deterministic transport theory

Mar. 2, 2016
50/74
Criticality eigenvalue problem
It is important to note that using the criticality eigenvalue form of the transport equation as a
representation of the physical transport problem is an approximation. Consequently, the solution
to (70) is not the the solution to the physical problem, but rather to a modified problem in which:
I The fission source term is scaled up when keff < 1
I The fission source term is scaled down when keff > 1

The result is that the contribution of fission source on flux distribution and spectrum is correspond-
ingly over- (keff > 1) or under-estimated (keff < 1).

The root cause of this issue is not in the way the transport problem is solved,20 but rather in its
formulation: a time-dependent system is forced into a steady-state condition by adjusting one of
the source terms. There is no solution to this problem, which can have a significant impact in the
results when the system is far from criticality.

In infinite-lattice calculations performed for the purpose of group constant generation the distortion
in flux spectrum is often adjusted using leakage corrections, which account for the contribution of
inward or outward net current, which in a realistic steady-state system would balance the source
and loss rates.

20
The same bias appears in the Monte Carlo method when the transport simulation is run in criticality source
mode, as discussed in Lecture 3.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
51/74
Other eigenvalue forms of the steady-state transport equation
The k-eigenvalue form is just one approach to approximating the system with a steady-state
transport equation. Similar scaling to attain balance between source and loss rates could be
accomplished using the collision eigenvalue, c, which adjusts the total source term:

1
Ω̂ · ∇ψ + Σψ = (F + S) (73)
c
The solution to this equation is again an approximation to the physical problem, with its biases
when c 6= 1.

PWR full-core reactor simulator calculations are often run to obtain critical boron concentration,
which can be formally thought of as an eigenvalue problem, in which the contribution to boron
absorption is separated from the removal term and adjusted to attain balance between source
and loss rate:
0
Ω̂ · ∇ψ + (Σ + NB σB )ψ = F + S (74)

where Σ0 is the total cross section excluding absorption to coolant boron, and NB and σB are
the atomic density and microscopic absorption cross section of 10 B, respectively.

Another example of a similar approach is introduced in Lecture 7 along with leakage corrections,
in which the balance between terms is attained by critical buckling search.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
52/74
Other eigenvalue forms of the steady-state transport equation
Another interesting approach to obtain a steady-state form of the transport equation is to assume
that the time-dependence of flux can be fully separated from the rest of the variables. Time-
dependence takes an exponential form and the separation of variables can be written as:
αt
ψ(r, Ω̂, E, t) = ψ(r, Ω̂, E)e . (75)

The substitution of this into the time-dependent transport equation (40), yields for the time-
derivative term:
1 ∂ αt α αt
ψ(r, Ω̂, E)e = ψ(r, Ω̂, E)e (76)
v ∂t v
When the exponential multiplier is canceled in all terms, the transport equation can be written as:
α
ψ + Ω̂ · ∇ψ + Σψ = F + S (77)
v
where α is the eigenvalue that can be adjusted to attain balance between the source and loss
rates.21 The criticality condition can be written correspondingly as:

< 0 ⇒ system is sub-critical

α = 0 ⇒ system is critical (78)

> 0 ⇒ system is super-critical


21
To be precise, α is the eigenvalue corresponding to the fundamental flux mode.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
53/74
Other eigenvalue forms of the steady-state transport equation
The additional multiplier α/v has the same units as macroscopic cross section, and it can be
interpreted as an absorption (α > 0) or production (α < 0) reaction with cross section inversely
proportional to neutron velocity. The reaction can also be interpreted as time-absorption or -
production, as seen by writing the probability per traveled path length:
α
dP = ds = αdt (79)
v

It should be noted that the α-eigenvalue form makes no adjustments to any of the physical source
or loss terms, and it can be shown that the solution to (77) corresponds to the asymptotic (funda-
mental mode) solution of the time-dependent transport equation.

The shape and spectrum of the flux is adjusted by the time-absorption / production term, which
can be understood by considering slow neutrons, for which the adjustment is most significant:
I When the system is super-critical, the slowest neutrons cannot keep up with the
exponentially increasing neutron population, and their contribution is scaled down by
time-absorption
I When the system is sub-critical, it is the slow neutrons that are dominating the
exponentially decreasing population, and their contribution is scaled up by time-production

The differences between the k- and α-eigenvalue solutions to the transport equation are best
illustrated by an example.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
54/74
Comparison of k - and α-eigenvalue mode
Example:Comparison of k - and α-eigenvalue mode

Consider the k- and α-eigenvalue solutions to transport equation in an infinite 20 cm thick

uranium plate with 30% 235 U enrichment:
1) The k-eigenvalue solution gives keff = 1.23, suggesting that the system is (prompt)
super-critical, and the balance between the source and loss rates is attained by adjusting
the fission term by factor 1/1.23 = 0.81.
2) The α-eigenvalue solution gives α = 1.04·107 1/s, which also means that the system is
super-critical. The balance between source and loss rates is obtained by adding a
time-absorption term with effective macroscopic cross section 1.04·107 /v 1/cm.

In the first case, the fission source is artificially adjusted to obtain a solution to the modified
transport equation.
In the second case, the solution to the asymptotic fundamental flux mode is obtained by
separating time-dependence from the rest of the variables, without modifications in the original
problem.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
55/74
Comparison of k - and α-eigenvalue mode
Example:Comparison of k - and α-eigenvalue mode

The differences are emphasized when the uranium plate is surrounded by an infinite water
reflector. The thermalization of neutrons has a major effect in the k-eigenvalue solution, in which
keff is increased to 1.46. The α-eigenvalue for the reflected geometry is 1.23·107 1/s.
The most dramatic differences between the two methods are seen in the corresponding flux
solutions. The k-eigenvalue solution completely ignores the time-dependence of flux, which
emphasizes the contribution of thermal neutrons returning to the fissile material with high
probability to cause new fissions.
In reality, the exponentially growing chain reaction is maintained by fast neutrons alone, and the
flux grows several orders of magnitude during the time it takes for a thermal neutron to return
from the reflector.
The α-eigenvalue flux represents the asymptotic solution to the time-dependent problem, and
reflects the correct shape of the exponentially growing flux. The contribution of thermal neutrons
is scaled down by time-absorption, which is inversely proportional to neutron speed.
The time-dependent behavior of the systems is illustrated in animations Lecture2_anim2.gif
and Lecture2_anim3.gif. It is seen that in the given time-scale, neutrons scattered to low
energy appear to be stopped in the reflector, unable to contribute in the continuation of the chain
reaction.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
56/74
Comparison of k - and α-eigenvalue mode
Example:Comparison of k - and α-eigenvalue mode
k-eigenvalue k-eigenvalue
0.03 α-eigenvalue 0.03 α-eigenvalue
Normalized flux

Normalized flux
0.02 0.02

0.01 0.01

0.00 0.00
−30 −20 −10 0 10 20 30 −30 −20 −10 0 10 20 30
X-coordinate (cm) X-coordinate (cm)

Figure 8 : k- and α-eigenvalue solutions to transport equation in an infinite slab geometry (20 cm thick
uranium plate with 30% 235 U enrichment). Left: unreflected system, Right: infinite water reflector. The
unreflected system is prompt super-critical, with keff = 1.23 given by the k-eigenvalue calculation. The
reflector increases the value to keff = 1.46. Both eigenvalue modes result in very similar flux shape for the
unreflected system, but the thermalization of neutrons in the reflected geometry is clearly emphasized in
the k-eigenvalue calculation. In reality, fast neutrons are sufficient for maintaining an exponentially-
growing chain reaction, and the fact that thermal neutrons follow far behind is not properly represented.
Calculated by k- and α-eigenvalue Monte Carlo simulation.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
57/74
Comparison of k - and α-eigenvalue mode
The previous example is an extreme case, presented to demonstrate the fact that the solution of
the k-eigenvalue form of the transport equation can lead to completely false conclusions on the
physical behavior of the system.

The k-eigenvalue method is still the most commonly used approach to solving reactor physics
problems, and when close to criticality, the flux solution can be a good approximation on what the
solution would look like if the system was critical.

It is important to note that the modeled system itself is often an approximation of physical reality,
for example, an infinite lattice of identical fuel assemblies. In such case, the asymptotic shape and
spectrum of the time-dependent flux (solution to α-eigenvalue problem) may not be any closer
to physical reality than the solution to the adjusted transport equation (solution to k-eigenvalue
problem).

There are also differences in the numerical algorithms used for solving different eigenvalue prob-
lems. α-eigenvalue problems can be subject to numerical instabilities when the system is well
below criticality, and k-eigenvalue problems are easier to solve using iterative techniques.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
58/74
Point-kinetics approximation: prompt neutrons
The simplest way to account for time dependence in deterministic transport theory is the so-called
point kinetics approximation, in which the flux is approximated by the fundamental mode solution,
assuming separation of variables:

ψ(r, Ω̂, E, t) ≈ ψ0 (r, Ω̂, E)T0 (t) (80)

where ψ0 is time-independent form function and T0 is the time-dependent amplitude function.

When the sub-indexes are dropped for convenience, the flux amplitude takes the form:
ωt
T (t) = T (0)e (81)

In the absence of delayed neutrons, the coefficient in the exponential can be written as:
ρ
ω= (82)
Λ
where reactivity ρ represents the branching of the fission chains and the amount of multiplication
in the system, and generation time Λ can be interpreted as the average time taken for a single
prompt neutron to reproduce itself by fission.
Coefficient ω is the inverse period, related to the reactor period T :
1
T = (83)
ω
which is time taken for the neutron population to grow by factor 2.7 (positive reactivity) or decrease
by factor 0.37 (negative reactivity).

Lecture 2: Deterministic transport theory

Mar. 2, 2016
59/74
Point-kinetics approximation: prompt neutrons
Reactivity can be written using the eigenvalue corresponding to the fundamental mode solution,
i.e. the effective multiplication factor, keff :

k0 − 1 keff − 1
ρ= = (84)
k0 keff

In other words, reactivity gives the fractional deviation of the system from criticality.22

In some text books, generation time is related to another time constant, known as the prompt
removal lifetime
τr = Λkeff (85)

which represents the effective mean time between the generation and the removal of a neutron
from the system by new fission.23

The emission of delayed neutrons complicates the solution of the transport equation, by coupling
fission rate to the decay of delayed neutron precursors. Time dependence of flux still follows an
exponential form, but the simple relation between ω and ρ in Eq. (82) is lost.

22 −5
One of the common units for reactivity is pcm (per cent mille), 1 pcm = 10 and 1000 pcm = 1%.
23
Since this time interval has the physical interpretation as the average lifetime of the emitted fission neutron, τr
is often called the prompt neutron lifetime.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
60/74
Point-kinetics approximation: delayed neutrons
As discussed in Lecture 1, delayed neutron emission results from the radioactive decay of certain
fission products, known as delayed neutron precursors. For example:

87 β− 87 ∗ 86 1
35 Br −−−−→ 36 Kr −→ 36 Kr + 0n (86)
55.7s

The decay of the excited state of 87 Kr is practically instantaneous, and the delay between fission
and neutron emission depends on the half-life of the precursor isotope 87 Br.

There are hundreds of fission products that act as delayed neutron precursors,24 but most of them
are either short lived or have very low fission yields. The half-lives of significant precursor isotopes
range from hundreds of milliseconds to almost a minute (for 87 Br above).

In practice, all precursor chains are not handled separately, but the reactions are instead divided
into a number of representative precursor groups with different lifetimes and neutron yields. The
most conventional representation uses six groups, with half-lives ranging from about 0.02 to 55
seconds.25

The physical total delayed neutron fraction β depends on the yields of prompt and delayed neu-
trons, which depend on isotope and neutron energy.

24
Fission yield data typically consists of about 1000 fission product isotopes and isomeric states.
25
The six-group representation is used in the ENDF/B and JENDL evaluated nuclear data files. The European
JEFF-3.1 file and later evaluations are based on eight precursor groups.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
61/74
Point-kinetics approximation: delayed neutrons

Table 1 : Ten most abundant delayed neutron precursor isotopes in thermal fission of 235 U and 239 Pu and
fast fission of 238 U. The yields are cumulative, and they refer to average nuclide production per fission.

235 239 238

Nuclide Half-life (s) U yield Pu yield U yield
87
Br 55.70 0.021 0.007 0.016
137
I 24.51 0.036 0.023 0.056
88
Br 16.50 0.018 0.005 0.019
138
I 6.46 0.015 0.007 0.040
93
Rb 5.80 0.035 0.017 0.046
89
Br 4.37 0.014 0.003 0.021
97
Y 3.75 0.021 0.012 0.032
94
Rb 2.70 0.015 0.007 0.034
98m
Y 2.00 0.020 0.019 0.026
143
Cs 1.79 0.016 0.006 0.035
141
Xe 1.73 0.016 0.005 0.034
99
Y 1.48 0.019 0.013 0.047
97m
Y 1.17 0.028 0.025 0.023
97
Sr 0.43 0.017 0.007 0.033
95
Rb 0.38 0.007 0.003 0.022

Lecture 2: Deterministic transport theory

Mar. 2, 2016
62/74
Point-kinetics approximation: delayed neutrons

Table 2 : Half-lives (in seconds), decay constants and relative yields of delayed neutron precursor groups
in 235 U, 238 U and 239 Pu fission. The physical total delayed neutron fraction depends on neutron energy.

235 238 239

U U Pu
j T1/2 λj ν j /ν d T1/2 λj ν j /ν d T1/2 λj ν j /ν d
1 54.51 0.013 0.038 52.38 0.013 0.013 53.75 0.013 0.038
2 21.84 0.032 0.213 21.58 0.032 0.137 22.29 0.032 0.280
3 6.00 0.116 0.188 5.00 0.139 0.162 5.19 0.134 0.216
4 2.23 0.311 0.407 1.93 0.359 0.388 2.09 0.332 0.328
5 0.50 1.398 0.128 0.49 1.406 0.225 0.55 1.263 0.103
6 0.18 3.872 0.026 0.17 4.030 0.075 0.22 3.209 0.035

Lecture 2: Deterministic transport theory

Mar. 2, 2016
63/74
Point-kinetics approximation: delayed neutrons

0.05 233
U 2.0 233
U
235 235
U 1.8 U
238 238
U U

Delayed neutron fraction (%)

0.04 239 1.6 239
Pu Pu
Delayed neutron yield

241 241
Am 1.4 Am
243 243
Am Am
0.03 1.2
1.0
0.02 0.8
0.6
0.01 0.4
0.2
0.00 −9 0.0 −9
10 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 101 10 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 101
Neutron energy (MeV) Neutron energy (MeV)

Figure 9 : Physical total delayed neutron yield (left) and fraction (right) of selected actinides as function of
neutron energy. The fraction also depends on the prompt neutron yield, which varies from nuclide to
nuclide and increases practically linearly as function of neutron energy.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
64/74
Point-kinetics approximation: delayed neutrons
The emission of delayed neutrons changes the fission source (47) into:
Z ∞
χ (E) 0 0 0
F (r, Ω̂, E, t) =(1 − β) νΣf (r, E )φ(r, E , t)dE
4π 0
Jd  (87)
χj (E)
X 
+ λj Cj (r, t)
j
4π

where β is the total delayed neutron fraction26 and the second term gives the emission rate in Jd
precursor groups. The delayed neutron source depends on precursor concentration Cj , which is
coupled to decay equations:
Z ∞
∂ 0 0 0
Cj (r, t) = βj νΣf (r, E )φ(r, E , t)dE − λj Cj (r, t) (88)
∂t 0

where λj is the decay constant and βj is the delayed neutron fraction in precursor group j . It is
also assumed that the precursor nuclides are not transmuted by neutron interactions.

As pointed out in Lecture 1, the concentrations of delayed neutron precursors follow changes in
fission power with a considerable delay. Consequently, the momentary precursor concentrations
in Eq. (87) depend on the immediate operating history.

26
For now, β is written without dependence on E .

Lecture 2: Deterministic transport theory

Mar. 2, 2016
65/74
Point-kinetics approximation: delayed neutrons
When the reactor has operated at constant power for a sufficiently long time, the precursor con-
centrations are in equilibrium, and
∂
Cj (r, t) = 0 (89)
∂t
The saturation concentrations are then given by:
Z ∞
1 0 0 0
Cj (r) = βj νΣf (r, E )φ(r, E )dE (90)
λj 0

which yields for the fission source term:

Z ∞
χ (E) 0 0 0
F (r, Ω̂, E) =(1 − β) νΣf (r, E )φ(r, E )dE
4π 0
Jd  Z ∞ (91)
χj (E)

0 0 0
X
+ βj νΣf (r, E )φ(r, E )dE
j
4π 0

Fraction (1 - β ) of fission neutrons is then emitted as prompt with spectrum χ and fraction β as
delayed with spectra χj , suggesting that in steady-state operation the impact of delayed neutrons
is seen only in the energy distribution of emitted fission neutrons.27

27
It is pointed out later that this is not exactly the case, as the time constants depend on the importance of
neutrons for the continuation of the fission chain.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
66/74
Point-kinetics approximation: delayed neutrons
In the point-kinetics approximation it is assumed that the spatial distribution of delayed neutron
precursor concentrations is separable from the time-dependence, which results from the similar
separability for flux and, consequently, fission rate.
The relation between inverse period and reactivity depends on delayed neutron constants βj and
λj . The relation is not easily resolved in simple closed form, but it can be written as:
Jd  
X ωβj
ρ = ωΛ + (92)
j=1
ω + λj

This is known as the inhour equation. When the number of precursor groups Jd is six, the time-
dependent amplitude function can be written using its seven roots (ω0 , ω1 , . . . ω6 )28 as:
6
X ωj t
T (t) = T (0) Aj e (93)
j=0

where constants Aj are:

6  
βi
X
Λ+
i=1
ωj + λi
Aj = (94)
6  
1 X βi λi
Λ+
1 − ρ i=1 (ωj + λi )2
28
Not to be confused with the ωn ’s of the different eigenmodes.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
67/74
Point-kinetics approximation: delayed neutrons
Components corresponding to smaller ωj decay faster, and the time dependence becomes dom-
inated by the solution determined by the largest root ω0 29
ω0 t
T (t) ' T (0)A0 e (95)

The corresponding period T0 = 1/ω0 is known as the asymptotic or stable reactor period, and it
can be measured experimentally by monitoring the rate of change in fission power.

The relation between (92) and (82) is seen by setting β = 0. This corresponds to the prompt
super-critical state,30 in which ρ  β .

The relation between ρ, ω and the time constants is illustrated in Fig. 10. The transfer to prompt
super-critical state is shown as a steep drop in the asymptotic period. It is also seen that:

lim ω(ρ) = −λ1 (96)

ρ→−∞

which means that the decrease rate of neutron population is limited by the slowest decaying
precursor group.

29
When ρ > 0, only ω0 is positive, leading to exponential growth. When ρ < 0, ω0 is the least negative root,
corresponding to slowest exponential decay.
30
Another common practice is to measure reactivity in units of delayed neutron fraction, or dollars. When reactivity
is above 1$, the system is in prompt super-critical state.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
68/74
Point-kinetics approximation: delayed neutrons

Λ = 0.1 ms (negative reactivity)

102
Λ = 0.1 ms
Λ = 0.1 µs
101

Asymptotic period (s)

ω = −λ1
0 ρ=β
Reactivity

100
β→0

10−2

β→0
10−4

10−6
−λ6 −λ5 −λ4 −λ3 −λ2 −λ1 0 10 100 1000 10000
Inverse period Reactivity (pcm)

Figure 10 : Relation between ρ, ω and the time constants. Left: Illustration of the 7 roots
ω0 , ω1 . . . , ω6 . of Eq. (92) with fictitious delayed neutron constants. The roots are found at the crossing
points between the curves and a horizontal line corresponding to ρ. Right: Asymptotic reactor period
corresponding to the largest root as function of reactivity. The top curve shows the negative of reactor
period (decaying amplitude function) as function of negative of reactivity.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
69/74
Point-kinetics approximation: summary
The point kinetics approximation relies on separation of variables, and essentially assumes that
the time-dependence of the system can be approximated by a single form function and the asso-
ciated amplitude function, and that all higher-order (transient) flux modes can be discarded.

The method is best applied to transients caused by uniform changes in the operating conditions,
or in small or otherwise closely coupled systems, for example:
I Criticality experiments
I Small research reactors
I Fast reactors

The approximation breaks down especially in localized transients in neutronically large systems,
such as PWR control rod ejections, which require more elaborate methods and time-space re-
solved solution.

LWR transient analyses are typically performed using reactor dynamics codes with time-dependent
one- or three-dimensional neutronics. This topic is covered in Lecture 11.

It should also be noted that the inhour equation (92) describes the dynamic behavior of the reactor
after a step change in reactivity. In reality, both positive and negative reactivity insertions take time,
and the results are not accurate when ρ is changing.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
70/74
Calculation of time constants
The delayed neutron fractions in the coupled fission source (87) and precursor concentration (88)
equations was written without energy dependence. This is not because β is independent of neu-
tron energy, but rather because the use of physical delayed fractions leads to inaccurate results
for the dynamic reactor response.
The reason is that some source neutrons are more likely to introduce long fission chains than
others, because they are born further away from boundaries or strong absorbers, in a region of
higher reactivity, or simply at a different energy. The point-kinetics approximation fails to capture
these effects, even though the separation variables is otherwise well justified.
In order to account for the fact that fission chains that are terminated early contribute less to the
asymptotic neutron population than those forming longer chains, the reactor time constants must
be calculated by importance-weighting. In deterministic transport theory this implies using the
adjoint neutron flux ψ † :
Z Z Z Z Z
† 0 0 0 0 0
ψ (r, Ω̂ , E )Bj (r, E, E ) ψ(r, Ω̂, E)dV dΩ̂dEdΩ̂ dE
V Z Ω̂Z EZ Ω̂0Z E 0
βeff ,j = Z
† 0 0 0 0 0
ψ (r, Ω̂ , E )F(r, E, E ) ψ(r, Ω̂, E)dV dEdΩ̂dEdΩ̂ dE
V Ω̂ E Ω̂0 E0
(97)
where F is the fission operator:
0 0
F(r, E, E ) = χ(E )νΣf (r, E) , (98)

and Bj is the delayed neutron emission source operator:

0 0
Bj (r, E, E ) = χd,j (E )νd,j Σf (r, E) (99)

Lecture 2: Deterministic transport theory

Mar. 2, 2016
71/74
Calculation of time constants
The adjoint flux is a measure of neutron importance, and it is the solution to the adjoint trans-
port equation, obtained by inverting the scattering and fission operators.31 Similar adjoint flux
weighting can be written from the calculation of prompt generation time:
Z Z Z Z Z
† 0 0 1 0 0
ψ (r, Ω̂ , E ) ψ(r, Ω̂, E)dV dΩ̂dEdΩ̂ dE
0 0 v
Λeff = Z Z VZ Ω̂Z E Z Ω̂ E (100)
† 0 0 0 0 0
ψ (r, Ω̂ , E )F(r, E, E ) ψ(r, Ω̂, E)dV dΩ̂dEdΩ̂ dE
V Ω̂ E Ω̂0 E0

The importance-weighted β and Λ calculated from (97) and (100) are called the effective delayed
neutron fraction and effective generation time, respectively.

Neutron importance and adjoint flux play a major role in sensitivity and uncertainty analysis, and
it also has use in shielding applications and variance reduction techniques used in Monte Carlo
simulation. These topics, however, are beyond the scope of this course.

31
In Monte Carlo simulation, the adjoint transport problem can be intuitively understood as tracking the neutron
histories backwards, for example, starting from a fission event and following the path to its origin. Keeping track of
the neutron’s descendants and calculating their contribution to the asymptotic population forms the basis of the
iterated fission probability (IFP) method, used by Monte Carlo codes to calculate importance-weighted time
constants.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
72/74
Summary of main topics
The modeling of an operating nuclear reactor is a complicated task, involving the couple solution
of neutronics, heat transfer from fuel to coolant, heat removal from the core by coolant flow and
fuel depletion.

Reactivity feedbacks and fuel burnup turn linear transport problem into a non-linear coupled prob-
lem. The solution is obtained by iterating between different solvers. The calculations are based
on a multi-stage scheme, in which the scale of the problem is gradually increased while simulta-
neously moving towards simplified physics.

All deterministic transport methods are based on the solution to the transport equation, which is
formulated based neutron balance in the six-dimensional phase space. The transport problem
cannot be solved without approximations.

Transport calculations can be roughly divided into time-dependent and eigenvalue calculations.
The most commonly used method is the k-eigenvalue criticality calculation, in which the balance
between neutron source and loss terms is obtained by artificially scaling the number of emitted
fission neutrons.

The simplest approach to time-dependent problems is the point-kinetics approximation, in which

the flux is approximated by the fundamental mode solution and separated into space-, direction
and energy-dependent form function and time-dependent amplitude function. The approxima-
tion is commonly used in small or otherwise closely-coupled systems, but the simulation of LWR
transients requires more elaborated time-space resolved solution.

Lecture 2: Deterministic transport theory

Mar. 2, 2016
73/74
Topics of next lecture
The next lecture (9.3.2016) introduces the Monte Carlo method, which provides an alternative
insight into the transport problem from the viewpoint of an individual neutron

Specific topics include:

I Basics of Monte Carlo calculation, probability and statistics
I CSG geometry model used in Monte Carlo codes
I Particle tracking routines
I Interaction physics
I Simulating the neutron population
I Collection of results

Lecture 2: Deterministic transport theory

Mar. 2, 2016
74/74