Transport and diffusion theory boundary

conditions
Aljoša Gajšek
November 23, 2020

Abstract
Transport and diffusion equation boundary conditions with the fo-
cus on continuity of current and flux at material boundaries, definition
and derivation of partial flux and zero flux at extrapolated boundary.

1
Contents
1 Introduction 3
1.1 Revision . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Limitations of diffusion theory . . . . . . . . . . . . . . . . . . 5

2 Boundary conditions 6
2.1 Transport equation boundary conditions . . . . . . . . . . . . 6
2.2 Diffusion equation boundary conditions . . . . . . . . . . . . . 6
2.2.1 Boundaries at interfaces . . . . . . . . . . . . . . . . . 7
2.2.2 Infinitesimally thin source . . . . . . . . . . . . . . . . 8
2.3 Diffusion theory vacuum boundary condition . . . . . . . . . . 8
2.3.1 Partial flows . . . . . . . . . . . . . . . . . . . . . . . . 8
2.4 Extrapolated boundaries . . . . . . . . . . . . . . . . . . . . . 11

3 Conclusion 13

2
1 Introduction
Before we start lets ask ourselves why are boundary conditions so important?
The determination of neutron distribution inside a reactor is a central prob-
lem of nuclear reactor theory, as it determines the rate of nuclear reaction.
The general transport equation (5), or energy dependant diffusion approx-
imation (6) are both partial differential equations that have derivatives in
both space and time, so for solving any given problems appropriate bound-
ary and initial conditions are required.

This paper is meant for a seminar at course of physics of fission

reactors at Faculty of mathematics and physics, University of Ljubljana. I
will assume the reader is familiar with the neutron transport and diffusion
equation. I will also strive to use same notations as were used during the
lectures. If external reader is reading this I would suggest reading the script
of the subject [1] or alternatively book [2], which was also my main source
for this seminar.

1.1 Revision
For the sake of completeness I will however rewind some key concepts and
equation that will be needed for understanding the article.

ϕ is a angular neutron flux is a scalar and is a product of neutron

density n and their speed v
~ t) = vn(~r, E, Ω,
ϕ(~r, E, Ω, ~ t) (1)
~j is angular current density, it is angular neutron flux, defined as
a current of magnitude ϕ and direction Ω̂

~ t) = Ω̂ϕ(~r, E, Ω,
~j(~r, E, Ω, ~ t) (2)
Neutron flux φ is sum of angular neutron fluxes over the complete
spacial angle

3
 Z
φ(~r, E, t) = ~ t)
dΩ̂ϕ(~r, E, Ω, (3)
4π
~ It is sum of all angular current densi-
Neutron current density is a vector J.
ties, and therefore represents the current of neutrons.
Z Z
~ r, E, t) =
J(~ ~ t) =
dΩ̂~j(~r, E, Ω, ~ t)
dΩ̂Ω̂ϕ(~r, E, Ω, (4)
4π 4π

Based on the power of Ω̂ in the equations (3,4), we call them first and second
moment of ϕ.
The most fundamental equation for describing the neutron popu-
lation is transport equation (5). It is however very difficult to solve, as it has
derivatives both in space and time, and integration over angle and energy of
scattered neutrons.

~ t)
1 ∂ϕ(~r, E, Ω, ~ · Ωϕ(~
+∇ ~ r, E, Ω,
~ t) + ΣT (~r, E)ϕ(~r, E, Ω, ~ t) =
vZ ∂t Z (5)
∞
~ →Ω
dE 0 Σs (E 0 → E, Ω ~ 0 )ϕ(~r, E 0 , Ω
~ 0 , t)dΩ
~ + s(~r, E, t)
0

Where s is source of neutrons and Σs scattering cross section.

The equation can be simplified by integration over angle Ω̂, this yield neutron
continuity equation (6). Loosing angular dependence creates new unknown
in the form of J~

1 ∂φ ~ r, E, t) + ΣT (~r, E)φ(~r, E, t) =
− ∇ · J(~
vZ∂t
∞ (6)
dE 0 Σs (E 0 → E)φ(~r, E 0 , t) + S(~r, E, t).
0

Assuming stationary conditions, P1 approximation and adding

equations (4,3). Allows us to simplify this equation even more into diffu-
sion equation

−∇ ~ · D(~r, E)∇φ(~
~ r, E) + Σt (~r, E)φ =
(7)
Z
S0 (~r, E) + Σ0s (~r, E 0 → E)φ(~r, E 0 )dE 0 .

4
 Since the diffusion equation was derived from less general transport
equation we can predict, there might be similar boundary conditions for both.

1.2 Limitations of diffusion theory

We mostly treat neutron motion as a diffusion process. In principle we as-
sume that neutrons are travelling from regions of high neutron density to
regions of low neutron density. Similarly as we would study the temperature
profile, or mixing of gasses. While those two can be very precise, that is
because of very frequent collisions, that makes particles move in a unpre-
dictable random directions (think about, Brownian motion that microscop-
ically describes diffusion), this is not always the case for neutron diffusion.

Scattering cross sections for neutrons are very small and in orders of
magnitude of ≈ 10−24 cm2 , mean free path for fast neutrons x = Σ1t = N1·σt ≈
cm. The other problem is that some critical components of nuclear reactor
are of similar dimensions as mean free path, such as nuclear fuel elements,
and therefore cannot be adequately described by the neutron diffusion. The
third problem is anisotropy of real nuclear reactors. The regions close to
source or boundary can have high angular dependences as seen from figure 1

Figure 1: Angular anisotropies in special cases. From [2]

5
2 Boundary conditions
2.1 Transport equation boundary conditions
First we will start with transport equation (5).
The conditions are dependant on a particular problem. Lets first
look at the vacuum boundary condition. We must also assume convex shape
of our reactor so that if the neutrons leave the reactor, there is no way they
will scatter back into it, because in vacuum there is nothing to scatter of
see Figure 2. The appropriate boundary condition would be zero neutrons
entering the medium from the outside

n(~
rs , E, Ω̂, t) = 0 if Ω̂ · es < 0. (8)
Or in terms of angular neutron flux:

ϕ(~
rs , E, Ω̂, t) = 0 if Ω̂ · es < 0, all rs on S. (9)

Figure 2: Reentry and non reentry surface. From [3]

2.2 Diffusion equation boundary conditions

The neutron flux φ is a real continuous function, so φ(+~r) = φ(−~r) ∀ ~r .
This might not be true for some useful mathematical models such as point
source, but this wont be a problem at our boundary. The neutron flux
boundary condition is therefore

6
 φ1 (~
rs , E, t) = φ2 (~
rs , E, t), (10)
where 1 , 2 are regions on each side of boundary and rs is a position vector
of boundary.

2.2.1 Boundaries at interfaces

Since the diffusion equation is just a simplification of more exact transport
equation, we will try to derive the appropriate boundary conditions from it.
Lets imagine the boundary between two regions of different cross sections.
The appropriate boundary conditions are

ϕ1 (~ rs , Ω̂, t) ∀ Ω̂
rs , Ω̂, t) = ϕ2 (~ (11)
where ϕ1,2 are angular fluxes for regions 1 and 2. This might not look to clear
at first but lets imagine, that we zoom in a lot directly at the bounder. On
that scale even small mean free paths might look big, and nothing would ac-
tually happen at the moment of passing the border. This condition satisfies
the conversation of neutrons in every single angle, but isn’t true for diffu-
sion theory. In diffusion theory we assume instant and isotropic scatterings.
We can only ensure that angular moments of equation (11) are satisfied. In
derivation of diffusion theory we only used first and second moments of an-
gular flux (φ(~r, t); J(~r, t)), our boundary conditions can only hope to satisfy
those two
Z Z
dΩ̂ϕ1 (~
rs , Ω̂, t) = rs , Ω̂, t) ⇒ φ1 (~
dΩ̂ϕ2 (~ rs , t) = φ2 (~
rs , t) (12)
4π 4π
Z Z
dΩ̂Ω̂ϕ1 (~
rs , Ω̂, t) = rs , Ω̂, t) ⇒ J~1 (~
dΩ̂Ω̂ϕ2 (~ rs , t) = J~2 (~
rs , t). (13)
4π 4π
We can substitute J in equation (13) and get the following bound-
ary conditions:

φ1 (~
rs , t) = φ2 (~
rs , t)
(14)
−D1 ∇φ1 (~
rs , t) = −D2 ∇φ2 (~ rs , t)
Or in other words as many neutron as pass the surface on one side,
must exit it on another. We can happily see that our logical thinking from
the start of subchapter was also satisfied in equation (12)

7
2.2.2 Infinitesimally thin source
Sometimes we might find it interesting to observe thin source of neutrons. In
between two different regions. How to set our boundary conditions in that
case?

Figure 3: Currents on each side of infinitesimally thin source

As seen from the picture 3, we have a thin source of neutrons

between two regions. φ must be analytical function, while current of the
source is J~ = ( S2 , 0, 0). The external current is J~0 = (J0 , 0, 0). The boundary
conditions are therefore:
φ1 (~
rs , t) = φ2 (~
rs , t)

eˆs · J~2 (~
rs , t) − eˆs · J~1 (~
rs , t) = (15)
S0 (t) S0 (t)
(eˆs · J~0 (~
rs , t) + ) − (eˆs · J~0 (~
rs , t) − ) = S(t)
2 2

2.3 Diffusion theory vacuum boundary condition

2.3.1 Partial flows
J(r,t)from equation (4) is number of neutrons travelling over specific area in a
unit of time. It is what we in usually call current. Here we will call it neutron

8
current density. The units of both neutron flux and current density are iden-
tical, but J as a vector quantity describes number of neutrons travelling over
a surface oriented in some direction. While φ simply describes all the neu-
trons that pass given surface, regardless of orientation. Although they seem
very similar and they are closely related, there is no simple relation between
φ and J. φ is useful, for describing reactor rates inside the nuclear reactor,
and J for neutron flow or leakage (for example from surface of a reactor).

For that purpose we will also introduce partial current densities,

J± . They correspond to the rate at which neutrons flow over area into positive
direction (from left to right) J+ and in negative direction (from right to left)
J− . From equation (4) and looking at the Figure 4 it is clear that:
Z ∞ Z
~
J± (r, t) = dE dΩ̂~j(r, E, Ω̂, t) (16)
0 ±2π

±2π is just another notation of performing angular integration in the direc-

tion of surface normal +2π and into opposite direction −2π

Figure 4: Current densities over area dA. From [2]

When we discussed Transport equation the boundary condition

was a logical statement that there can be no incoming neutrons at a vacuum
boundary.

rs , Ω̂, t) = 0, for Ω̂ · dS < 0, all r~s on S.

ϕ(~ (17)
Again diffusion will only be able to approximate this boundary condition.
We can see from equation (17) that we only have a boundary condition over
half a surface of solid angle (only for incoming neutrons). Lets try to find
the boundary condition by integrating the above equation:

9
 Z Z
rs , Ω̂, t)Ω̂ · ês dΩ̂ =
ϕ(~ ~j(~
rs , Ω̂, t) · ês dΩ̂ = J− (~
rs , t) = 0 (18)
−2π −2π

Even this boundary condition can only be approximated since dif-

fusion theory cannot yield the exact form for J± . We will again resort to the
P1 approximation, in which
1 ~ r, t))
(φ(~r, t) + 3Ω̂ · J(~
ϕ(~r, Ω̂, t) = (19)
4π
We will also need a couple others equations:

Z Z Z
4π
J~ = −D∇φ, dΩ = 4π, ΩdΩ = 0, Ω2 dΩ = .
4π 4π 4π 3
The first equation comes from continuity equation, the others three we have
discussed during the lectures.
Z
J± (~r, t) = ês · Ω̂ϕ(~r, Ω̂, t)dΩ̂ (20)
±2π

Z
1 ~ Ω̂
J− = ês · Ω̂( (φ + 3Ω̂ · J))d (21)
−2π 4π
Lets rewrite the integral of the sum as the sum of integrals. If we remember
φ, J were both derived by integration over angle, and they are both angle
independent so we can easily write them outside the integral.
Z Z
1 3 ~
J− = φ ês · Ω̂dΩ̂ + J · ~es Ω̂2 dΩ̂ (22)
4π −2π 4π −2π

On the right side we have known integral (over half a sphere), on the left we
can rewrite omega by definition:

dΩ̂ = sin(θ)dθdς (23)

where θ and ς are azimuth and polar angles. We can choose (rotate) our
coordinate system such that the scalar product will be equal to cosθ

~es · Ω̂ = cosθ (24)

10
and for our half a sphere (both angles go from 0 to π) our equation becomes:
Z Z Z
~es · Ω̂dΩ = cos(θ)sin(θ)dθdς (25)
π π

We solve it by well known trick sinxdx = d(cosx). We must also

change the borders to equal our new variable.

Z Z π Z 1
1
~es · Ω̂dΩ = cos(θ)d(cos(θ))dς = 2π =π (26)
0 0 2
the partial solution to right side of equation is
Z
4π 1
Ω2 dΩ = . (27)
2π 3 2

If we combine both parts and write J~ = −D∇φ we get the final

solution as:
1 D
rs , t) + eˆs · ∇φ(~
J− (~r, t) = φ(~ rs , t) = 0 (28)
4 2

2.4 Extrapolated boundaries

Lets apply those boundary conditions to one dimensional geometry, so that
d
ês = 1 and ∇ = dx
1 D dφ
J− (xs ) = φ(xs ) + |x = 0 (29)
4 2 dx s
dφ φ(xs )
|xs = − (30)
dx 2D
Equation (30) gives us the gradient of flux near the border of
medium. We can use it to extrapolate point of zero flux.
dφ
φ(x) = φxs + ∆x = 0 (31)
dx

∆x = 2D (32)

11
 2
x
e = xs + (33)
3Σtr
Where x e is zero flux extrapolated boundary. In the last step we
have used D = 3Σ1tr . Where Σtr is a transport cross section, and is a cor-
rection of total cross section for anisotropic scattering. More on that can be
found in [1] (3.39).

The vacuum boundary condition

J~− (xs ) = 0, (34)

can therefore be replaced, by somewhat easier

φ(e
xs ) = 0. (35)

We must however, ask ourselves about the validity of that result.

We can imagine, as in figure 3 neutron flux near the source of neutrons. If we
imagine a source in the shape of a giant plane the neutron flux won’t even
change regardless of position (in reality neutrons undergo beta decay and
transform into lighter protons, with half life of t 1 ≈ 15min). Our boundary
2
for the outside world is nothing more than a source of neutrons, and their flux
should not change as drastically near the boundary. It can however change
due to geometry of the problem according to the Gauss law. Comparison
between diffusion flux and measured can be seen on figure 5

12
 Figure 5: Neutron flux in diffusion theory and measured values. From [2]

More advanced transport theory calculations of extrapolated bound-

ary result in a bit different value, for plane geometry:
1
x
e = xs + 0.7104 (36)
Σtr
Extrapolated value is in reality also dependant on the geometry of
the surface, but the curvature can be neglected unless its radius is comparable
to the mean free path.

3 Conclusion
Neutron diffusion equation, will be very important tool for nuclear reactor
analysis. It works very well for big isotropic systems, but we must be careful
when using it, in the areas where the irregularities might be comparable
to the mean free path, and near the boundary of the medium. In those
cases transport equation is a better tool, but in real cases the Monte Carlo

13
method will be used. When applying the neutron diffusion equation we can
make a problem easier by applying zero flux at extrapolated border boundary
condition. But we must be careful, not to take this value too seriously it is
not the real flux but just different boundary condition, as in real cases the
neutron flux does not fall to 0!

References
[1] Andrej Trkov, Luka Snoj and Matjaž Ravnik, REAKTORKSA IN RA-
DIACIJSKA FIZIKA Študijsko gradivo, Ljubljana, 2.10.2013).

[2] James J. Duderstadt and Louis J. Hamilton NUCLEAR REACTOR

ANALYSIS, Department of Nuclear Engineering The University of
Michigan, JOHN WILEY & SONS, Ann Arbor Michigan, 1976

[3] Estimating the volume of a 3D image using the alpha shape algorithm
in R, URL: https://towardsdatascience.com/estimating-the-volume-
of-a-3d-image-using-the-alpha-shape-algorithm-in-r-61169ac00020, date
published 11.9.2009

14