The One Speed Diffusion Theory

Model
Dr. Khurrum Saleem Chaudri
P&D Wing, D-Block
Reference Book: Nuclear Reactor Analysis
by
James J. Duderstadt
Louis J. Hamilton
Why One Speed??
 One Speed = Single Average Kinetic Energy
 One speed model allows preliminary design estimates and is much
easier to solve
 Any model characterizing all of the neutrons in a reactor by a
single speed (hence energy) and treating their transport from
point to point as a diffusion process cannot be expected to yield
accurate quantitative estimates!!!
 If the cross sections appearing in this theory are chosen properly,
one can use the one-speed diffusion model to make preliminary
design estimates
 Also this will provide as a basis to solve and analyze more
sophisticated models (multigroup diffusion), cause of being
identical
Example of Cross-Section
Heuristic Derivation of One Speed
Diffusion Equation
 Heuristic: Refers to experience-based techniques for problem
solving, learning and discovery that give a solution which is not
guaranteed to be optimal
 Advantage of one speed?? Greatly simplified mathematical study of
nuclear reactor behavior
 Point of suspicion for using one speed theory?? Neutron energy
range 10-3 – 107 eV
 But, we will later see that if we regard one speed approximation as
an average description and choose the appropriate average x-
sections, we can actually get quantitative description of nuclear
reactor
 Most energy dependent theories (MG diffusion theory) are solved
by performing a sequence of one-speed calculations for each
successive energy group
Heuristic Derivation of One Speed
Diffusion Equation (Contd…)
 Let us characterize the neutron distribution in the reactor by
the neutron density N(r,t) = (#/volume) at position r at
time (t)
 Let us define flux ø = vN(r,t) (#/cm2-sec)
 Flux can be multiplied by corresponding macroscopic cross-
section to get certain reaction rates e.g. ∑f(r)ø(r,t) = rate of
fission reaction per cm3 at point r
 Continuity Equation = balance equation i.e. time rate of
change of number of neutrons in a given volume (∂n/∂t) is
simply the difference between the rate of production and rate
of disappearance (absorption + leakage) for the said volume
Heuristic Derivation of One Speed
Diffusion Equation (Contd…)
 Let us consider arbitrary volume “V” of surface area “S”
located anywhere within the reactor
 Control volume = mathematical abstraction for volume
under consideration
 Total number of neutron in “V” at time “t”
3 1
V V r v  (r , t )
3
d rN ( r , t ) d

 Time rate of change of number of neutrons in “V”

d  3 1  3 1  ( r , t )
 d r  ( r , t )   r
 d
dt V v  V v t
Heuristic Derivation of One Speed
Diffusion Equation (Contd…)
 Partial derivative as flux is a function of position and time
and derivative is with respect to “time” only
 Time rate of change of neutrons in volume “V” = Production
in “V” – Absorption in “V” – Net Leakage from “V”
 Because this is a heuristic derivation so let us write different
terms
Production in “V” =  d rS (r , t ) where S(r,t) = neutron
3
 V

source density basically based on (fission + external source)

 Absorption rate density at any point “r” is ∑(r)ø(r,t). Hence
total rate of absorption in “V” =  d 3ra (r ) (r , t )
V
Heuristic Derivation of One Speed
Diffusion Equation (Contd…)
 Neutron leakage is defined by using neutron current density, then
net rate at which neutrons pass out through a small surface
element dS at rs is J(rs,t).dS. Hence total net leakage through the
surface of “V” is just =  dS .J (r , t )
S

 Production and absorption are volume integrals but leakage is

defined using surface integral to combine them together, we can
use Gauss’s theorem to convert surface integral into volume
integral
 dS .J (r , t )   r.J (r , t )
3
d
S V

 Gauss’s theorem is basically “sum of all sources minus the sum of

all sinks gives the net flow out of a region” or “divergence
represents epansion or compression of J i.e. net flow in or out
Heuristic Derivation of One Speed
Diffusion Equation (Contd…)
 Substituting all the terms in basic continuity equation, we get
3  1  
V d r 
 v t
 S   a  . J 

0

 This equation is for any arbitrary volume so we can simply

write this equation as
1 
 S   a  . J  0
v t
or
1 
 .J   a  S
v t
Heuristic Derivation of One Speed
Diffusion Equation (Contd…)
 This equation contains two unknowns i.e flux and current
and we have only one equation
 To close the system we need another equation relating flux
and current
 There is no exact relationship between flux and current so
we have to approximate
 Borrowing from other fields of physics (gaseous diffusion),
current density is approximately proportional to the negative
spatial gradient of the density (in our case flux) i.e. particles
will tend to flow from regions of high to low density at a rate
proportional to the negative density gradient
Heuristic Derivation of One Speed
Diffusion Equation (Contd…)
 This relationship is known as Fick’s Law i.e.
J (r , t )   D(r ) (r , t )
 Here D(r) = Diffusion coefficient
 When using Fick’s law, we have to remember the underlying
assumptions which come with it. These assumptions are
 Medium is infinite (i.e. not valid near boundaries)
 Medium is uniform (not strictly true as reactor is
heterogeneous)
 Neutron flux varies slowly in position (Not valid near sources
or sinks)
 Isotropic scattering in Lab Coordinate System (LCS)
Heuristic Derivation of One Speed
Diffusion Equation (Contd…)
 Chapter 4 deriving diffusion approximation, we got

 
1
D 1  3  t   o  s 
3tr  

 Fick’s law is valid provided it was used to describe

 Neutron flux several mean free path (mfp) away from
boundaries or isolated sources
 The medium was only weakly absorbing and
 Provided the neutron current was changing slowly on a time
scale comparable to the meantime between neutron-nuclei
collisions
Heuristic Derivation of One Speed
Diffusion Equation (Contd…)
 Putting Fick’s Law in general continuity equation we derived
earlier, we arrive at the one-speed neutron diffusion equation
i.e.
1 
 .D(r )   a (r , t )  S (r , t )
v t
Initial and Boundary Conditions
 In order to solve diffusion equation, we need suitable initial
and boundary conditions
 These IC and BC came from neutron transport theory
(Chapter 4), we will use physical argument (heuristics) here
 This is also a reason of good results we get from diffusion
theory although so much approximations and limitations
 Initial Condition:  (r,0)  o (r) for all r (i.e. steady state flux)
 Boundary conditions are a bit more complicated and depend
on the type of physical system under study
Vacuum Boundary
 At the outside boundary of a reactor, one would like to
construct a boundary condition corresponding to the fact
that no neutrons can enter the reactor through this surface
from outside (Fick’s Law not applicable near boundary)
 This BC has implicit assumption that the reactor is
surrounded by an infinitely large vacuous region (in real life
reactor is surrounded by air, concrete, steel etc.)
 It is frequently convenient to assume that the reflection of
neutrons back into the reactor from such materials is
negligible so that non-reentrant BC applies
Vacuum Boundary (Contd…)
 Diffusion theory is incapable of exactly representing a non-
reentrant BC
 The closest one can come to is demand inward partial
current to vanish on boundary (clearly an approximation as
expression of J- already is approximate)
Vacuum Boundary (Contd…)
 This leads to d =0.7104 tr
 Since D = diffusion coefficient = tr/3 and D ≈1 cm
 Hence d = .7104 x 3 ≈ 2.1 cm << reactor size (meters)
 Hence ø(surface) ≈ 0
 This means that one can simply ignore the extrapolation
length and simply assume that the flux vanish on the true
boundary
Interfaces (Material Discontinuities)
 The structure of nuclear reactor core is extremely complex,
containing fuel, structural material, coolant, control regions
 So if homogeneous medium is not possible (∑(r)), at least we
can use “sectional uniform” medium and hence have cross
sections that change abruptly across an interface separating
regions of differing material compositions and solve diffusion
equation for each material separately. But what at interfaces??
 Physically fictitious but mathematically expedient
convenience of including an infinitesimally thin neutron
absorber or source at the interface
   
es . J ( rs )  J ( rs )  S
ˆ
 
Other Types of Boundary Conditions
 Flux must be non-negative, real and finite (will help for cases
when medium is infinite or sources are present)
 Occasionally use symmetry boundary condition to rule out
physically irrelevant solutions of diffusion equation e.g. for
one dimensional slab with center at origin ø(-x)=-ø(x) is
discarded)
Summary of the One-Speed Diffusion
Model
 Initially we will use to describe the neutron population in a
nuclear reactor by using neutron diffusion equation
1 
 .D( r )   a ( r , t )  S ( r , t )
v t
 Along with Initial Condition  (r, 0)  o (r ) for all r
 And boundary conditions:
 Free surface:  (rs , t )  0 or J  (r s , t )  0 
 Interface: flux and normal component of current should be
continuous across interface
 0≤ø(r,t)<∞ (except in the neighborhood of localized sources)
Summary of the One-Speed Diffusion
Model (Contd…)
 In many cases, we will deal with situations for which the
medium in which the neutrons are diffusing is uniform or
homogenous such that D and ∑a do not depend on position
 Then the one-speed diffusion equation simplifies to
1 
 D 2   a (r , t )  S (r , t )
v t
 Laplacian operator depend upon co-ordinate system being
used
 We will frequently consider situations in which the flux is
not a function of time, hence steady state reactor
D   a (r )  S (r )
2
Summary of the One-Speed Diffusion
Model (Contd…)
 This equation is known as the Helmoholtz equation
a S (r )
 2 (r )   (r )  
D D
1 S (r )
  2 (r )  2  (r )  
L D

 Here “L = Neutron Diffusion Length √(D/∑a)

 Neutron diffusion length is essentially a measure of how far
the neutrons will diffuse from a source before they are
absorbed
Elementary Solutions of the Diffusion
Equation
 Apply equation (A) to study diffusion of neutron from a
steady-state source in a non-multiplying medium
 Plane Source in an Infinite Medium
 Infinitely wide plane source located at the origin of an infinite
(no boundary), homogenous medium
 Source is assumed to be emitting neutrons isotropically at a rate
S0 (n’s/cm2-sec)
 Since both the source plane and the medium are of infinite
extent, the neutron flux can only be a function of distance “x”
from the source plane
Plane Source in an Infinite Medium
1 S (r )
 Hence  2 ( r )   ( r )  
L2 D
 Becomes d 2 1 So
  ( x )     x
dx 2 L2 D

 Where source is modeled by a Dirac delta function

 If we restrict x ≠ 0, the source term disappears to give us
d 2 1
2
 2  ( x)  0  for x  0
dx L
 Our approach is to solve this homogenous equation for x ≠ 0
and then use a BC at x = 0 to “fix up” these solutions
Plane Source in an Infinite Medium
 We can obtain this BC directly by integrating from
x 0  x  0
 Across the source plane and then taking the limit to find
d d 
D |  D |  J x (0 )  J x (0 )
dx dx
 So
 Similar to interface BC, if we use the symmetry of the
geometry to assert that
J x (0  )  J x (0  )  J  0 
S0
 lim J
x 0
 x 
2
Plane Source in an Infinite Medium
 This equation implies that net neutron current at the origin
on either side must be just half of total source strength as
source is isotropic
 As we have second order derivative term, we will need
another BC
 We will use the BC of finite flux as x approaches ∞
d 2 1
2
 2
  x  0 for x  0
dx L
with BC ' s
d S0
 a  xlim
 0
D 
dx 2
 b  lim   x   
x 
Plane Source in an Infinite Medium
 The solution for x<0 will be deduced using symmetry
condition
 General solution of equation then becomes

 x   x 
  x
 A exp     B exp  
 L   L 
Applying BC  b 
 B  0
While BC  a  gives
d DA  x 
D  exp   
dx L  L 
Applying lim
x 0

DA S0 S0 L
  A 
L 2 2D
Plane Source in an Infinite Medium
 Hence our solution becomes
S0 L  x 
  x  exp    x  0
2D  L
and by symmetry, we can say
S0 L  x 
  x  exp   x  0
2D  L
 Hence the neutron flux falls off exponentially as we move
away from the source plane with a characteristic decay length
of “L= √(D/∑a)” so larger “L” means smaller absorption
cross section so less neutron flux attenuation as we move into
the medium
Plane Source in an Infinite Medium
 Magnitude of flux is proportional to the source i.e. if S = 2S0
then ø=2 ø0
 This solution provides a reasonable description of neutron
flux in a medium (provided absorption cross section is not
too high) while this solution is not valid within certain mfp of
source plane itself
Point Source in an Infinite Medium
 Isotropic (no angle dependence) point source emitting S0
n’s/sec at the origin of an infinite medium
 Diffusion equation in spherical coordinates reduces to

1 d 2 d 1
2
r  2  r   0 for r  0
r dr dr L
with BC ' s
 a  lim
r 0
4 r 2 J  r   S0
 b  lim
r 
 r   
Solution will be

 r   A

exp  r 
L B
exp  r 
L
r r
Point Source in an Infinite Medium
 Applying BC (b) gives B=0
 Applying BC (a)
d
J r    D r   D
dr

 D
d 

A

exp  r
L



dr  r 
 

 DA 

 exp  r
L  
exp  r
L 
 
 r 2
Lr 
 

S0  lim DA 

 exp  r
L 

exp  r
L   
r 0  r2 Lr 
 
S0
S0  4 DA  A 
4 D
Point Source in an Infinite Medium
Hence

 r  
S0 exp  r  L 
4 rD
 An interesting application of this result is to calculate the
man-square distance to absorption in a non-multiplying
medium
 Let us say that number of neutrons absorbed between “r” and
“r+dr” is just  S exp   r L  
 4 r dr    
0 2
 4 rD  a

 
or in terms of probability in "dr"

p  r  dr 
D
r


exp  r
L  dr
  
 a 

r
p  r  dr  2 exp  r
L

L
dr 
Point Source in an Infinite Medium
 We can then calculate

r 2
  dr r 2 p  r 
0

r3   r 
  dr exp    dr
 L 
2
0 L
1   r 
2 0
 r exp  
3
 dr
L  L 
Applying integration by parts
r2  6 L2
or
1
L 2
 r2
6
Point Source in an Infinite Medium
 Hence neutron diffusion length “L” has the interesting
physical interpretation as being (1/√6) of the rms distance to
absorption

 We have calculated rms distance from source to the point of

absorption, not the total path length traveled by the neutron
Point Source in an Infinite Medium
 For graphite L=59 cm
 <r2>=(√6)L=144 cm
 Mfp of thermal neutron = 2.5 cm and average number of
collisions suffered by the neutron before absorption is 1500
then its apparent average path length = 1500 x 2.5 > 3700
cm which is considerably larger than (√6)L=144 cm
1 2
 L 2
 r can be used to define neutron diffusion length in
6
situations in which diffusion theory would not apply (i.e.
strongly anisotropic scattering or large absorption)
 Then first determine flux ø(r) resulting from isotropic point
source using transport theory and then you can calculate “L”
 d r r  r 
3 2
1 2
L2  r 
 d r  r 
3
6
Finite Slab Geometries
 Isotropic plane source imbedded at the center of a slab of
non-multiplying material of width “a” surrounded on both
sides by vacuum
 Diffusion equation becomes

d 2 1
   x  0 for x  0
dx 2 L2
with BC ' s
d S0
 a  xlim
0

D 
dx 2
 a 
b     0
 2 
Solution will be
 x   x 
  x  A exp     B exp  
 L   L 
Finite Slab Geometries
 Applying BC (b)


a  a   a 
   0  A exp     B exp  
2  2L   2L 
 a 
 B   A exp   
 2L 
  a  x  a 
   x   A  exp     exp  
  L  L 
Applying BC  a 
  x  1  x  a   1  S0
 AD  exp       exp     
  L  L  L  L  2
AD   x  x  a  S
lim  exp     exp     0
x 0 L   L  L  2
AD   a   S0
1  exp  L    2
L   
1
S L   a  
 A 0 1  exp  L  
2D   
Finite Slab Geometries
1
S L   a     x   x  a 
  x  0 1  exp  L    exp   L   exp  
2D        L  
  x   a  x 
exp     exp   
S0 L   L  L 
  
2D  
 a 
1  exp  

  L  

From symmetry
  x   a  x    a 
 exp     exp      exp  
  x 
S0 L   L   L  .  2L 
2D   a     a 
 1  exp     exp  

  L  
   2L 
a2 x 
Sinh  
  x 
S0 L  2 L 
2D Cosh a 
2L 
Finite Slab Geometries
 The sketch of solution reveals somewhat similar result to
infinite medium except dropping off more rapidly near the
boundaries due to neutron leakage (solution not valid within
certain mfp of source and boundary)
Finite Slab Geometries (Multi Region
Problem)
 A variant on the problem involves replacing the vacuum by a
material of different composition than the slab itself
 General procedure for attacking such multi-region problems
is to seek solutions of the diffusion equation characterizing
each region and then match these solutions using interface
BC
 We can use geometrical symmetry to solve for 0<x<∞
 In region (1), we will seek a solution to
d 21 1 a
2
 2
1  x   0 0 x  for region 1
dx L1 2
Similarly
d 22 1 a
2
 2
1  x   0  x  for region 2 
dx L2 2
 Finite Slab Geometries (Multi Region
Problem)
Boundary Conditions
S0
 a  xlim
 0
J 1  x  
2
 b  2  x   
In addition to these, we have interface conditions
 c  1  a 2   2  a 2 
 d  J1  a 2   2 
J2 a 
Using our earlier work as a guide, we can seek general solutions
 x   x 
1  x   A1 cosh    B 1 sinh   in region 1
 1
L  1
L
 x   x 
2  x   A2 exp     B 2 exp   in region  2 
 L 2   2 
L
Finite Slab Geometries (Multi Region
Problem)
 After applying the boundary conditions and simplifying, we
get
Finite Slab Geometries (Multi Region
Problem)
 The comparison of the solution with bare slab shows that flux in
the central region falls off somewhat more slowly when the
vacuum is replaced by a diffusing medium
 This is due to the neutrons scattering back into the slab that would
have otherwise been lost to vacuum
 Such materials are used to reduce neutron leakage are known as
reflectors (Large ∑scatter and small ∑abs) e.g. water (PWR),
graphite (HTGR)
 A concept very closely related to that of neutron reflectors is the
reflection coefficient or albedo, defined as the ratio between the
current out of reflecting region to the current into the reflecting
region
Finite Slab Geometries (Multi Region
Problem)
J out
 
J in

 Suppose we want to attach a reflecting slab of thickness “a” to

a reactor core
 If we are given a current density Jin = J+ entering the
reflecting slab surface, which we locate at x=0 for our
convenience, then we can solve the diffusion equation
characterizing the reflecting region
d 2 1
2
 2  ( x)  0 for 0  x  a
dx L
subject to boundary conditions J +  0   J in ,   a   0
Finite Slab Geometries (Multi Region
Problem)
 We can solve for the flux and then calculate the albedo which
will come out to be
 2D   a
1   coth  
   L   L
 2D   a
1   coth  
 L   L
Finite Slab Geometries (Multi Region
Problem)
2D
1
    L
2D
1
L
 For example, the albedos characterizing infinitely thick
reflectors of graphite, water and heavy water are 0.93, 0.82
and 0.97 respectively
 The albedo can be used to replace the detailed solution in the
reflecting region by an equivalent BC at the edge of core
J
using   out

J
in
Finite Slab Geometries (Multi Region
Problem)
 That is, if we use our earlier definitions of the partial current
densities J±, then the effective BC on the flux in the region
x<0 is just
1 d 1  1 
D |x  0    
 dx 2  1  

 Used in this manner, the albedo becomes a very useful device

for obtaining BC for reactor core calculations
Numerical Methods for Solving the
Neutron Diffusion Equation
 Thus far we have confined our attention to neutron diffusion in
homogenous (or region wise homogenous) media since in this case
the one speed equation could be solved analytically
 However in any realistic reactor calculation the heterogeneous
nature of the core must be taken into account
 One, not only, must consider non uniformities corresponding to
fuel pellets, cladding material, moderator, coolant, control
element but also spatial variations in fuel and coolant densities due
to non-uniform core power densities and temperature
distributions as well
 So analytical solution is out of picture and more preferable path of
numerical solution should be chosen
Numerical Methods for Solving the
Neutron Diffusion Equation
 General procedure is to rewrite the differential diffusion
equation in finite difference form and then solve the resulting
system of difference equations on a computer
 Suppose, we wish to solve  D d      x   S  x 
2

2 a
dx
 Subject to BC’s characterizing a finite slab width “a” i.e
ø(0)=ø(a)=0 (Ignoring extrapolation length for convenience)
 Discretize spatial variable “x” by choosing a set of N+1
a
discrete points equally spaced (for convenience) a distance  
N
Numerical Methods for Solving the
Neutron Diffusion Equation
 Re-write diffusion equation at each of these discrete points
xi, but to do so we need an approximation for d2ø/dx2
 Suppose we Taylor expand flux at xi±1 in-terms of its value at
point xi
d  2 d 2
i 1    xi 1   i   |i  2
|i
dx 2 dx
d  2 d 2
i 1    xi 1   i   |i  2
|i
dx 2 dx
 Adding these two expression, we get

d 2 i 1  i 1  2i

dx 2
2
Numerical Methods for Solving the
Neutron Diffusion Equation
 Sufficiently small values of step size is chosen (in order to
insure slow variation of flux) (a solution (kind of) near source
or sink for limitation of diffusion theory)
 This three point central difference formula should be a
reasonable approximation to the values of d2ø/dx2 at the
point xi
 Using this difference formula to re-write diffusion equation
   2i  i 1 
 D  i 1    ai  Si when i  1, 2,.....(N-1)
  2

Re arranging the equation
D  2D  D
  i 1   2   
a  i  i 1  Si
 2
    2

ai ,i 1i 1  ai ,ii  ai ,i 1i 1  Si when i  1, 2,.....(N-1)

Numerical Methods for Solving the
Neutron Diffusion Equation
 Diffusion equation has been reduced to a set of (N-1)
algebraic equation equations for N+1 unknowns i.e.
(ø0,ø1,ø2,……..,øN)
 Using two boundary conditions N-1+2 = N+1 solutions can
be obtained
Derivation of Difference Equations for
One Dimensional Diffusion Problems
 We will now consider the more general form of the one
dimensional diffusion equation in plane geometry
d d
 D  x   a  x   x   S  x 
dx dx
 This equation will be subject to BC or interface conditions
that we will leave arbitrary for the moment
 Rarely composition varies in a continuous way from point to
point
 Rather system properties are assumed to be essentially
uniform in various sub-regions of the reactor core
(represented by spatially averaged or “homogenized”
properties in each sub-region)
Derivation of Difference Equations for
One Dimensional Diffusion Problems
 First step is to set up discrete spatial mesh with non-uniform
mesh spacing

 Various methods to generate difference equation

representation of this one dimensional equation (we used
Taylor series expansion in general case)
 A more common scheme is to integrate the original
differential equation over an arbitrary mesh interval and then
to suitably approximate these integrals using simple means
values or difference formulas
Derivation of Difference Equations for
One Dimensional Diffusion Problems
 Simplest scheme is approximate the integral by value of
integrand at midpoint times the integration interval
 i 1
xi   i  i 1 
 2
dx a  x    x    aii  2 
2 
i
xi 
2  
 i 1
xi   i  i 1 
 2
dxS
S  x 
i  
2 
i
xi 
2  2 
The derivative term requires a bit more work
 i 1
xi  d d  x  d   x  xi  2i 1
 i
2
dx D  x  D  x | i
xi  dx dx dx xi 
2 2

To handle the derivative terms, ue can use a simple

two point difference formula
d  x    i
| i 1  i 1
dx xi 
2
 i 1
d  x    i 1
| i  i
dx xi 
2
i
Derivation of Difference Equations for
One Dimensional Diffusion Problems
 Furthermore, we will use centered average for “D”
  i 1  1
D  xi     Di 1  Di   Di ,i 1
 2  2
   1
D  xi  i    Di 1  Di   Di ,i 1
 2  2
Hence diffusion term can now be written as
d d xi  2i1
 i 1 
xi  d
 i
2
dx D  x  D  x | i
xi  dx dx dx xi  2
2

Di ,i 1  Di ,i 1 Di ,i 1  Di ,i 1
 i 1     i  i 1
i   i 1 i   i 1
Derivation of Difference Equations for
One Dimensional Diffusion Problems
 Combining relevant terms, we arrive at a set of difference
equations very similar to earlier results i.e.

ai ,i 1i 1  ai ,ii  ai ,i 1i 1  Si

where
 Di  Di 1  1
ai ,i 1    
  i   i   i 1
 Di 1  Di Di 1  Di  1
ai ,i   ai    
  i 1  i   i   i 1
 Di 1  Di  1
ai ,i 1    
  i 1   i   i 1
Derivation of Difference Equations for
One Dimensional Diffusion Problems
 Hence, once again, we have arrived at a set of N-1 three
point difference equations for N+1 unknowns discretized
fluxes
 In the particular case in which mesh size is constant and the
coefficients D(x) and ∑a(x) do not depend on x, we again get
same result as derived in previous section by Taylor series
expansion
 Coefficients aij depend only on a single subscript j, however
double scripting is useful for we will rewrite these algebraic
equation as matrix equations
Derivation of Difference Equations for
One Dimensional Diffusion Problems
 Final task is to append to these equations two additional
equations taking into account the BC
 We could simply use the vacuum extrapolated BC i.e.
ø0=øN=0 (taking care to place the mesh points x0 and xN on
these extrapolated boundaries)
 More general BC (such as non-reentrant current) can be
developed by taking the final two difference equations in the
set as
a0,00  a0,10  S0
aN , N 1N 1  aN , N N  S N
Derivation of Difference Equations for
One Dimensional Diffusion Problems
 The coefficients aij will depend upon the scheme used to
derive the difference equations!!!
 Fortunately if the mesh spacing is small, these differences will
be insignificant in actual calculations
 Since the spatial variation of the flux is essentially
characterized by diffusion length L, one generally chooses a
mesh spacing which is less than neutron diffusion length “L”
 Similar 3 point difference equations will also arise in
curvilinear geometries with one dimensional symmetry
 For convenience, we will assume region wise uniform
properties
Derivation of Difference Equations for
One Dimensional Diffusion Problems
 For cylindrical co-ordinates
 d 2 1 d 
D  2
    a  r   S  r 
 dr r dr 
Spherical co-ordinates
 d 2 2 d 
D  2
    a  r   S  r 
 dr r dr 
Hence difference equations come out to be
ai ,i 1i 1  ai ,ii  ai ,i 1i 1  Si
Where now
D  c 
a i,i-1 =- 1 
2 
 2i  1 

2D
a i,i =  a
2
D  c 
a i,i+1 =- 2 1 
  2i  1 

where c=0,1,2 for plane, cylindrical or spherical geometries
Derivation of Difference Equations for
One Dimensional Diffusion Problems
 For vacuum BC ø0=øN=0 and S0=SN=0
 Solve for i=1….i=N-1, eliminate ø0,øN
 For cylindrical geometry, at rN, øN=0
 At origin, we use symmetry to imply
d
|r  0  0  0  1 , S 0  S1
dr
 Again we ignore i=0 and i=N
 Very similar considerations hold for the case of spherical
geometry
Solution of 3-point Difference
Equations
 Suppose we have developed an appropriate set of difference
equations which we must now solve for discretized fluxes øi
Solution of 3-point Difference
Equations
 Where A is an (N-1 x N-1) matrix, and ø and S are (N-1)
dimensional column vectors
 A is the finite difference diffusion matrix which is tridiagonal
 This situation arises only for such ID geometries in which 3
Point difference equations are used
 In 2 and 3-D geometries, more complicated matrix
structures are encountered
  A1S
 Such tridiagonal matrices can be inverted directly using
Gaussian elimination (“forward elimination-backward
substitution” method)
Solution of 3-point Difference
Equations
 The general scheme is to subtract ai,i-1/ai-1,i-1 times the (i-1)th
row from ith row to eliminate the ai,i-1 element in the ith row
 At each step, the ith row is then divided by its diagonal
element and the procedure is continued
 For example
Solution of 3-point Difference
Equations
 Such that we eventually arrive at a system of equations of the
form
Solution of 3-point Difference
Equations
 Gaussian elimination used as it is not easy to find the inverse of A
 The Gaussian elimination of Forward Elimination and Backward
Substitution can be used to directly solve the difference equations
 This scheme is particularly important since it frequently appears
as an integral part of the iterative methods used in 2 and 3D
diffusion problems
 For this reason it is useful to formalize Gaussian elimination a bit
by noting that what we have in fact accomplished by forward
elimination is the factorization of the matrix A into a product of a
lower “L” and upper “U” triangular matrix
Solution of 3-point Difference
Equations
Derivation of Multidimensional
Difference Equations
 Most detailed neutron diffusion calculations characterizing
nuclear reactors require either 2D or 3D treatments
 Such details are particularly important in studying power
profiles in large reactors subject to non-uniform fuel loading
and depletion
 Hence we must consider the numerical solution of the more
general diffusion equation
.D  r    a  r   r   S  r 
 Once again the geometry of interest is discretized into a
mesh of cells (rectangular grid)
Derivation of Multidimensional
Difference Equations
 A 2 and 3D dimensional grid example is shown here

 Like last time, integrate diffusion equation over the spatial

volume of a given cell, using this to define the spatially
averaged cell properties
Derivation of Multidimensional
Difference Equations

 Here the sum (in gradient term or leakage term) is taken

over adjacent mesh point neighbors j=1,….,J where J=2, 4
or 6 in 1-,2-, or 3-D cartesian geometries
Derivation of Multidimensional
Difference Equations
 Where the mesh coupling coefficients lij are determined by
the particular mesh geometry and finite difference scheme
one chooses
 For example in Cartesian coordinates using essentially the
approximation schemes discussed earlier i.e.
lij  Dij / ij2
 Difference equation representing diffusion equation then
takes the form
J
Dij  J Dij 
  j    2   a
 i  Si
 ij 
 j 1  ij
2
j 1 
where "i" runs over all mesh points
Iterative Solution of Multidimensional
Difference Equations
 For solution, first task is to cast the set of equations into
matrix form
 This requires assigning a single index to each mesh point.
Why???
 For example, in a 2D mesh we could label mesh points as
 (i,j) k = i+(j-1)(N-1)
 For a 2D problem, the matrix structure takes the following
form in this case, for a five by four mesh point array
Iterative Solution of Multidimensional
Difference Equations

 Notice that the tri-diagonal form of 1-D is now augmented

by two additional side band diagonals (as 2D case yields a five
point difference formula. For 3D, we get 7 point difference
formula)
Iterative Solution of Multidimensional
Difference Equations
 Now let us consider how we might solve such systems i.e.
invert such matrices (A-1S=ø)
 Since Gaussian elimination can be applied to any matrix
(formally, at least) we might first consider applying this
technique to obtain a direct inversion
 For 1-D case, the FE sweep reduces original tri-diagonal
matrix to a form with
Iterative Solution of Multidimensional
Difference Equations
 For a two dimensional problem, similar forward sweep
results in filling of all zero entries between the main and
outer diagonal

 Problem with this approach is the requirement of

considerably more computer memory to allow for such
direct inversion of matrix
 Also such a direct algorithm is quite complicated to program
Jacobi-Richardson (Point-Jacobi)
Scheme

 “D” can be easily inverted

 Here enters the “iterative” philosophy in picture
 Flux “ø” is present on both left and right hand side
 Guess the RHS and then calculate LHS
 Process goes on until convergence
 Draw back of this method is its slow convergence
Jacobi-Richardson (Point-Jacobi)
Scheme
 Why slow convergence??
 Only a small chunk of the original matrix (its diagonal) is
being inverted on each iteration hence takes a lot of iterations
to invert the matrix completely
 For better understanding of scheme, let us write it out
explicitly in terms of the algebraic system
Jacobi-Richardson (Point-Jacobi)
Scheme
 Or it can be written as

 It can be seen that Jacobi scheme does not use all of the
available information during each iteration
 If the equations are solved in sequence from i=1 to i=N, as
they would be on a computer, then the solution of the first
equation yields (m+1)st guess of ø1 while for the same
iteration ø2 uses mth guess of ø1
 A better solution is to use the guess as soon as it is generated
Gauss-Seidel (Successive Relaxation)
Method
 In this method, the system of equations in each iteration is
solved as
Gauss-Seidel (Successive Relaxation)
Method
 This can be rewritten in matrix form by decomposing “A”
into the sum of an upper and lower triangular matrix
Successive Over Relaxation (SOR)
Method
 It is possible to accelerate the convergence of the iteration
scheme even further by introducing an acceleration
parameter to extrapolate the iterative flux estimate (SOR)
 The first step is to computer the Gauss-Seidel estimate,
which we label as (m+1/2)th guess for øi

 And then (m+1)st value for øi is calculated as a linear

combination of (m+1/2) and the previous SOR iteration
Successive Over Relaxation (SOR)
Method
 i.e.
i m1  i m1/2  1    im when   1  2 
 Iterative algorithm for each element can then be written as
Comparison of Methods
One Speed Diffusion Theory of a
Nuclear Reactor
 Thus far we have studied neutron diffusion in non-
multiplying media with one-speed diffusion theory
 Now we wish to apply this theory to nuclear reactors with
fissile material present
 Life cycle of neutron in a nuclear reactor (thermal) can be
shown as
One Speed Diffusion Theory of a
Nuclear Reactor
 For a fast reactor, neutron life cycle is a little modified
The Fission Source Term
 Let us assume that diffusion, absorption and fission all occur
at the same energy
 Fission source can then be given as S f  r , t    f   r , t 
 If this is the only source of neutrons in the reactor (prompt
neutrons), then the appropriate diffusion equation becomes
1 
 .D   a  r , t    f   r , t 
v t
 Here we have multiple components in macroscopic
absorption cross section i.e.
 a   amoderator  structure
a   coolant
a   fuel
a

when
fuel
a   fuel
   fuel
f
The Time-Dependent Slab Reactor
 Uniform slab of fissile material characterized by ∑a ,∑tr ,∑f
 This very simplified and unrealistic case helps clear many
important concepts
 Appropriate mathematical equation for this case is
1   2
D   a  x, t    f   x, t 
v t x2

with initial condition   x, 0   o  x   o   x 

a   a 
and boundary condition   , t      , t   0
2   2 
 As compared to time-independent case, now we have partial
differential equation
 One way to solve it is the separation of variables
The Time-Dependent Slab Reactor
  x, t     x  T t 
 Substituting this form back in diffusion equation and dividing
by   x  T t  we would get
1 dT   d 2

 D       x    constant  -
T dt   dx 2 f a

 LHS (function of time) is equal to RHS (function of space),
hence both must be equal to some constant which is
unknown for the time being
dT
 -T  t 
dt
d 2 
D 2
  f   a   x  -   x
dx v
The Time Dependent Slab Reactor
 Time dependent equation can be solved easily
T  t   T  0  e  t
 T(0) is the initial value to be determined later
 To solve the space dependent part
d 2   
       a   x   0
D f
v 
2
dx
a  a
with boundary condition          0
2  2
which is of the form of eigen value problem
d 2
 n  x   0
2
2
B
dx
a  a
with  n     n     0
2  2
The Time Dependent Slab Reactor
 As per the set boundary condition, we are only interested in
symmetrical solution
Eigenfunction  n  x   cos Bn x
 n 
2

Eigenvalues Bn2    with n=1,3,5...

 a 
 Using this analogy on equation at hand, we have

  va  vDBn2  v f  n when n=1,3,5...

 These values are known as time eigen values of the equation
 General solution (both time and space combined) becomes

n x
  x, t    An exp(nt ) cos
n a
odd
The Time-Dependent Slab Reactor
 Using orthogonality, we find
2 a2 n x
An   a dx0  x  cos
a 2 a
 Thus we found that the flux (for any symmetrical initial
distribution) can be represented by a superposition of modes
weighted by an exponential factor
n
n  v a  vDBn2  v f while Bn 
a
 Initially we assumed solution to be variable but then both solutions
were superimposed to yield a non-separable function of space and
time
 Hence separation of variables does certainly not imply a separable
solution
Long Time Behavior
 Interestingly , there is one very important situation in which
separability will occur, that involving the neutron flux for
very long times
Notice B  B  ...B   n / a 
2 2 2 2
1 3 n

Hence 1  3  5  7  ...
 This means that the modes corresponding to larger n decay
out more rapidly in time
 If we wait long enough, then only the fundamental mode
remains
  x, t  A1 exp  1t  cos B1x as t  
The Time-Dependent Slab Reactor
 This implies that regardless of the initial shape of flux, it will
decay into the fundamental mode shape
 Coefficient of fundamental mode A1 is non-zero as

2 a
 x
A1   2
a dx0  x  cos
a  a
2
 For sufficiently large value of fission cross section, -lamba
will become positive corresponding to an exponentially
growing flux
 However the same argument will hold since 1  3  5 ...
 Hence regardless of flux grows or decays, it will eventually
approach a fundamental cosine distribution
The Time-Dependent Slab Reactor
 It is customary to refer to this value of fundamental mode as
 
2

B1     Bg  Geometric Buckling
2 2

a
 This curvature is used since Bn2 is a measure of the curvature
of the mode shape i.e.
1 d 2 n
B 2

n
 n dx 2
 Larger buckling means larger curvature showing larger
leakage
 Hence we might expect that the mode with least curvature
will persist in time the longest
The Time-Dependent Slab Reactor
The Criticality Condition
 Criticality condition is about making the flux in reactor time-
independent (in the absence of sources other than fission)
 If we write the general solution of the flux as

  x, t   A1 exp  1t  cos B1 x   An exp  nt  cos Bn x
n 3
odd

 Requirement for a time-independent flux gives us that

1  0     a   f    DB12
Rewriting the criticality condition using the B12  Bg2
 f   a
 Bg2
D
It has become customary to refer to
 f   a
 Bm2  material buckling
D
The Criticality Condition
 Material buckling depends only on material composition
while geometrical buckling depends on core geometry
 Hence our criticality condition can be written as
 (Material composition) Bm2= Bg2(core geometry)
 To obtain criticality we adjust material composition along
with reactor size
Bm2  Bg2  1  0  supercritical
Bm2  Bg2  1  0  critical
Bm2  Bg2  1  0  subcritical
The Criticality Condition
  f /  a 
1   a 1  L B  
12 2


 2 g 2
 1 L Bg 
 f  Fuel a
Fuel
 f
  f  k
a aFuel
a
Also
Rate of neutron absorption
Rate of neutron absorption plus leakage

 d 3 r  a
 V

 d 3 r  a   d 3 rD 2
V V

a  d 3 r 1
= V

a V
d 3 r  DBg2  d 3 r
V
1  L2 Bg2 
 neutron non-leakage probability  PNL
The Criticality Condition
 Then

1  1 
 PNL  
 a  1  L Bg 
2 2
  a 
 l  Neutron lifetime in a finite reactor
Thus
 f /  a
k= fPNL 
1  L2 Bg2
Thus we can identify our fundamental time
eigenvalue as just the inverse of the reactor period
k 1 1
-1  
l T
The Criticality Condition for More
General Bare Geometries
 Only quantity characteristic of the reactor size or geometry
appearing in k or PNL is the geometric buckling
 Geometric buckling varies for different shape of geometries
 For multi-region reactors expressions become very
complicated
 Consider a bare reactor in critical state (time independent) of
uniform composition surrounded with vacuum boundary
conditions
 D    a  r    f   r 
2

with   rs   0
The Criticality Condition for More
General Bare Geometries
Or
  f   a 
     r   0
2

 D 
  f /  a  1 
 2    r   0
 D / a 
 k 1 
 2    2   r   0
 L 
Again equation is identical to
 2 n  Bn2 n r   0 when  n  rs   0
 Again like slab reactor we get eigen functions and eigen
values with only the lowest eigen value surviving
The Criticality Condition for More
General Bare Geometries
 Reactor criticality condition will be similar to slab reactor
i.e.
  f   a 
B 
2
m   B1
2
 Bg
2

 D 
Thus we can continue to use
PNL = 1  L2 B 
2 1
g and k   f /  a  1  L2 B
2 1
g 
for more general bare geometries provided we identify
the geometric buckling Bg2 as the smallest eigenvalue
B12 of the Helmohotz equation
 2  Bg2  r   0
The Criticality Condition for More
General Bare Geometries
 As discussed earlier, we don’t get any information about flux
level
 If the usable energy produced per fission event is wf, then the
thermal energy deposited in the core per unit volume per
second is just given in terms of the fission reaction rate
density as q  r   wf  f   r 
 This is just the local thermal power density at position r in
the core
 Total power generated by the core is just the integral of the
power density over the core volume i.e. P  V d 3rwf  f   r 
Example: A Right Circular Cylindrical
Core
 DO IT YOURSELF
 It will lead you to

 v0 r  z 
  r, z   AJ 0   cos  
 R   H 
when
  
2 2
 v0 
B 
2
g   
 R  H 
and
3.63P
A
wf  f V
Geometric Buckling and Critical Flux
Profiles for Some Common Geometries
Reflected Reactor Geometries
 To illustrate the complications that arise with multiregion
core geometries, let us consult slab reactor again
Reflected Reactor Geometries
 We will proceed directly to examine the time-independent
diffusion equations that must be satisfied by the fundamental
mode flux shape
 We will seek a solution for each region separately
d 2 C
  Ca   f  
a
Core: D C
2
C C
 x  0 0 x
dx 2
R d 
2 R
a a
Reflector: D   R
a  R
 x  0  x b
dx 2 2 2
with boundary conditions

C 
a R  a 
a      
2 2
a a
b JC    JR  
2 2

R  
a
c  b  0
2 
Reflected Reactor Geometries
 This problem will have no solution unless we choose the
proper combination of core composition and size
 We should anticipate that a criticality condition relating these
core characteristics would emerge in the course of our
analysis. General solution for the core is
 C  x   AC cos Bm
C
x
Sin term discarded due to prevalent symmetry
Also
 Cf   C
B C2
m  a

DC
While for reflector
 a 
  b  x 
 R  x   AR sinh  2 R 
 L 
 
D 
1/ 2
where LR  R
/ a
R
Reflected Reactor Geometries
 Applying the interface boundary conditions to find
B a
C
 b 
AC cos  m   A R sinh  R 
 2  L 
 BmC
a DR  b 
C C
D B A sin 
m
C
  R
A R
cosh  R 
 2  L L 
dividing both expressions
 Bm
C
a DR  b 
C
D B C
m tan    coth  R 
 2  LR L 
 This equation represents a relation between reactor
composition and size that must be satisfied if a solution to the
steady state diffusion equation to exist
 Hence this is just the criticality condition for this particular
geometry
Reflected Reactor Geometries
 This doesn’t look anything like Bm2  Bg2
 This, in fact, is a transcendental equation i.e. no explicit
solution for the critical size or composition
 Numerical or graphical techniques must be used
 Rewriting equation to draw a graph
 BmC a   BmC a  DRa  b 
  tan   C R
coth  R 
 2   2  2D L L 
 If we plot the LHS against  BC
m a / 2
, there will be many
intersections but we will choose only the lowest value of  B C
m a / 2
Reflected Reactor Geometries

 From the graph we notice that

BmC a   
2

  
2
or BmC
2 2 a
 The width required for criticality is somewhat smaller after
reflector addition
Reflected Reactor Geometries
 It is conventional to define the difference between bare and
reflected core dimensions as the reflector savings
   a(bare)  a(reflected )
 For the case at hand, reflector savings can be written as
1  D C BmC LR  b 
  C tan 
1
tanh  R 
Bm  DR  L  
For a thick reflector b LR
DC R
  R L
D
 Which is a measure of maximum reflector savings that can be
realized
Reflected Reactor Geometries
 Reflectors serve another function besides reducing neutron
leakage
 They flatten the flux and power distribution in reactor core
 One speed diffusion theory is not adequate to describe this
effect which results in a peaking of thermal flux in the
reflector region, as shown in figure
 This need multigroup treatment
Numerical Criticality Searches