PROPAGATION OF NEUTRON WAVES
THROUGH HETEROGENEOUS MULTIPLYING
AND NONMULTIPLYING MEDIA
By
VICTOR RALPH CAIN
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMIENUTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
August, 1965
To Bliz
ACENOWL;EDGMENTS
The author wishes to express his appreciation to his supervisory
committee for their assistance. Special appreciation is due to the
committee chairman, Dr. Rafael Perez, for his unusual patience and
constant encouragement and guidance.
Acknowledgment must also be rendered to the Department of Nuclear
Engineering and to the Computing Center of the University of Florida for
their combined financial support of the computational work for this
dissertation.
Finally, the author wishes to acknowledge his debt to the Oak
Ridge National Laboratory and to the U.S. Atomic Energy Commission for
the assignment which made this work possible.
TABLE OF CONTENTS
PAGE
iii
v
vii
1
4
21
30
58
79
88
101
104
ACKNOWLEDGNENTS . . . . . . . . . . . .
LIST OF FIGURES . . . . . . . . . . . .
ABSTRACT .. . . .. . .. . . . . .. . ..
INTRODUCTION .. . . . . . . .. .. . . .. .
CHAPTER I. USE OF A ONEGROUP ONEDIMENSIGNAL GREEN'S
FUNCTION IN HETEROGENEOUS MEDIA .........
II. USE OF A ONEGROUP TWODLKENSIONAL GREEN'S
FUNCTION IN HETEROGENEOUS MEDIA .........
III. USE OF AN AGEDIFFUSION TWODLMENSIONAL GREEN'S
FUNCTION IN HETEROGENEOUS MEDIA .........
IV. CALCULATIONAL RESULTS ...............
V. SUMMARY AND CONCLUSIONS ..............
:ES
A. The Complex Arithmetic Version of GINV .......
B. Digital Computer Techniques for the Calculation
of the Error Function with a Complex
Argument .. . . . . .. . . . ..
C. Galanin's Thermal Constant for a Sinusoidally
Varying Flux .. . .. . . . . . ..
D. Physical Constants Used in Sample Calculations ...
:ES .. . .. . . . . . . ..   
APPENDIC
REFERENCE
LIST OF FIGURES
FIGURE
PAGE
1. Geometry for the OneGroup, OneDimension Green's Function
Calculation ............................................. 5
2. OneGroup Flux Amplitude vs Distance Along the Z Axis for
One Cadmium Foil in Gra~phite, at 100 eps, with 1, 2, 4, and
10 Spatial Modes .................. ................... ....... 66
3. OneGroup Flux Amplitude vs Distance Along thle Z Axis for
One Cadmium Foil in Graphite, at 1000 cps, with 1, 2, 4,
and 10 Spatial Modes ................... .................. ... 67
4. OneGroup Flux Phase Lag vs_ Distance Along the Z Axis for
One Cadmium Foil in Graphite, at 100 eps, with 1, 2, 4, and
10 Spatial Modes .......................................... 68
5. OneGroup Flux Phase Lag vs Distance Along the Z Axis
for One Cadmium Foil in Graphite, atl1000 cps, with 1, 2,
4, and 10 Spatial Modes .................. .................. 69
6. Flux Amplitude vs Distance Along the Z Axis for Two
Cadmian or Uranium Foils in Graphite, at 200, 500, and
800 eps, with One Spatial Mode ................... ...****** 71
7. Flux Amplitude vs Distance Along the Z Axis for Two
Cadmium or Uranium Foils in Graphite, at 200, 500, and
800 eps, with Six Spatial Modes ......................****** 72
8. Flux Phase Lag vs Distance Along the Z Axis for Two
Cadmiun or Uranium Foils in Graphite, at 200, 500,
and 800 eps, with One and Six Spatial Modes ................ 73
9. OneGroup Flux Amplitude vs Source Frequency for Two
Cadmium Foils in Graphite, at z = 2, 8, 8.89, 15, and
30 em, with Six Spatial Modes .............................. 74
FIGURE PAGE
10. OneGroup Flux Phase Lag vs Source Frequency for Two
Cadnium Foils in Gratphite, at z = 2, 8, 8.89, 15, and
3O cm, with Six Spatial Modes .............................. 75
11. AgeDiffusion Flux Amplitude vs Distance Along the Z Axis
for Two Uranium Foils in D2O, at 200 cps, with One and
Six Spatial Modes ................... ................... .... 76
12. Components of AgeDiffusion Flux Amplitude vs Distance
Along the Z Axis for Two Uranium Foils in D20, at
200 cps, with Six Spatial Modes ............................ 78
13. Regions of Complex Plane Requiring Different Techniques
for Calculating the Error Function ......................... 97
Abstract of Dissertation Presented to the Graduate Council
in Partial Fulfillment of the Requirements for the
Degree of Doctor of Philosophy
PROPAGATION OF NEUTRON WAVES THROUGH HETEROGENEOUS
MULTIPLYING AND NONNUTLTIPLYING MEDIA
By
Victor Ralph Cain
August 1965
Chairman: Dr. Rafael B. Perez
Maj or Department: Nuclear Engineering
Prediction of the behavior of neutron waves in heterogeneous media
is performed using onegroup diffusion theory and agediffusion theory. The
onegroup theory is used in two different, finite geometries which allow
reduction to a onedimensional problem and to a twodimensional problem.
These cases result in diffusion kernels, or Green's functions, for the two
finite configurations.
The theory is then extended to include Fermiage, or continuous
slowingdown, theory for higher energy neutrons. The approach is similar
to the FeinbergGalanin heterogeneous reactor theory except that it is a~p
plied to a finite geometry and includes time dependence. A finite diffu
sion kernel is obtained which is similar to the results of the simpler
calculations. In addition, a finitemedium, Fermiage kernel results which
describes the behavior of the slowingdown neutrons in the finite geometry.
Both the onegroup and the agediffusion developments are used to
calculate numerically sample configurations which are suitable for experi
mental verification. These include rectangular assemblies of graphite and
heavy water which have poison or fuel rods inserted.
In order to demonstrate better the improvements of this work over
the FeinbergGalanin theory, an extension to critical assemblies is made.
This offers the possibility of doing criticality calculations for physically
small assemblies which are not amenable to homogenization techniques.
Using the agediffusion theory results, it is also shown that data
from neutron wave experiments may be processed in such a way as to give
experimental measurements of both the diffusion and the slowingdown
kernel.
v111
INTRODUCTION
In 1960 a program was initiated at the University of Florida in the
field of neutron wave propagation. The initial effort was in the determina
tion of diffusion and thermalization parameters in moderating media. It was
shown by Perez and Uhrig (1)1 that these parameters were related to the
damping coefficient and phase shift per unit length of the neutron wave.
The work of Hartley (2) and Booth (3) confirmed these predictions, and
lately Perez and Booth (4) were able to perform accurate measurements of
thermalization parameters in graphite using the neutron wave technique.
The above mentioned studies were performed in homogeneous neutron moderat
ing media. Some work has been done by Booth (5), in heterogeneous media,
and by Denning et al. (6) in tworegion moderating media. Denning et al.
were able to show a definite reflected wave from the interface between the
two media which interferes with the incident wave.
The goal of the present work is to predict, as accurately as is fea
sible, the propagation of neutron waves in heterogeneous multiplying and
nounultiplying media. The incentive to perform this work is twofold. From
the purely academic viewpoint, the propagation of waves through periodic
structures is of sufficient interest in the fields of quantum mechanics,
sound theory, solidstate physics, and others that extension of the
lUnderlined numbers in parentheses refer to corresponding numbers
in the list of references.
previously mentioned work in heterogeneous media becomes de s irab le. That
neutron wave propagation theory is applicable to these other fields may
perhaps seem reasonable by pointing out that the neutron wave propagation
we are concerned with here is a macroscopic phenomenon, qualitatively
similar to sound wave propagation. One simply replaces gas density with
neutron density. Quantitatively, the neutron wave propagation is the more
difficult problem since nuclear reactors are highly dispersive and absorp
tive media in which disturbances of the neutron field are quickly damped.
Secondly, there is a strong incentive from the practical viewpoint
of gaining understanding about heterogeneous reactors. Traditionally,
homogenization techniques based on the WignerSeitz unitcell theory have
been used to design heterogeneous reactors, and it is only recently that
theories accounting for the interactions between fuel plates are evolving
from the works of Galanin ( ), Feinberg (8), and others.
The study of neutron wave propagation through heterogeneous media
involves analytical techniques and assumptions similar to the developments
used in heterogeneous reactor theory. Hence studies of this type offer the
possibilities of experimental tests of the various theories in a way which
enlarges the realm of application of purely static experiments. Also, any
disturbance of the neutron field in a heterogeneous reactor can be analyzed
as a superposition of neutron waves, thus yielding information on possible
shortrange instabilities for a particular reactor design.
Theoretical studies in this field are quite complicated because of
the interplay of the effects of the spatial heterogeneities produced by the
fuel plates and the energy distribution of the neutron population. Before
meaningful but costly experiments can be performed, the theoretical founda
tions have to be developed, which is the goal set forth in this disserta
tion.
At first sight, the obvious approach to the theoretical work is
simply to modify the FeinbergGalanin theory to include time dependence so
that it can be applied to subcritical assemblies. In practice, this is not
satisfactory since most reactors are rather large assemblies. Because of
this the present heterogeneous reactor theory assumes that the moderating
medium in which the fuel is embedded is infinite spatially. Practical
experimental assemblies for neutron wave propagation experiments, however,
tend to be quite small in comparison. Consequently, it was necessary to
start from the basic equations including the finiteness of the feasible
experimental assemblies.
CHAPTER I
USE OF A ONEGROUP ONEDIMENSIONAL GREEN'S
FUNCTION IN HETEROGENEOUS MEDIAL
The :problem to be considered is the calculation of the transmission
of a thermal neutron wave, according to onegroup diffusion theory, in a
simple heterogeneous geometry. The geometry chosen is such that it permits
a practical experimental setup and also results in an analytical form reduc
ible virtually to a onedimensional problem.
Figure 1 shows the geometry of the problem. The semiinfinite
parallelepiped of moderating material extends from x = a to x = a, y = b
to y = b, and z = 0 to z = oo. The transverse dimensions a and b are assumed
to include the diffusion theory extrapolation distance. The perturbation
in the medium is taken to be in the fonn of very thin sheets, the kth sheet
extending from x = xo to x = xo and y = b to y = b and of thickness h,
centered around z = zk. The system is driven by a sinusoidally modulated
source of thermal neutrons placed at the face z = 0.
The timedependent, onegroup diffusion equation describing this
system is
1 at(ryt) = U [Z + 2 (r)3 ] (r, t) (1.1)
v at a a
1The solution obtained in this chapter is an extension of work
originally performed by R. S. Booth.
Fig. 1. Geometry for the OneGroup, OneDimension Green's Function
Calculation.
where
v = the neutron speed,
Q(rjt) = the spatial and timedependent neutron flux,
r = a position vector,
t = time,
D = the diffusion coefficient (units of length),
ye = the La~placian operator,
La = the macroscopic cross section of the moderating media,
6Ca = the spacedependent macroscopic cross section of the
absorbing sheets.
The spacedependent perturbation may be written as
Ba(r) = BlaR( u x xo u(x x )
k=1
u z z k 3 ,(1.2)
where Blak the macroscopic cross section of the kth foil, u(y) = the
unit step function, and M = the total number of foils. If h is taken to be
sufficiently small, the first two terms of the Taylor's series expansion
for the unit step function in z may be used in Eq. (1.2), reducing it to
Bla(r) = h lu(x + xg) (x x)] a~,k 8(z z) (1.s)
where 6(y) is the Dirac delta function.
It should be noted that if the sheets contain fissionable material,
8a,k is replaced by 82a,k Ff~k, where v is the neutron multiplicity
and 82f~k is the macroscopic fission cross section of the _kth foil.
The sinusoidally modulated neutron source is assumed to be isotropic
with a Maxwellian energy distribution and can be represented by
S(r~t) = So + Re [Slr)(= e iwt (1.4)
where i =9 J, w is the angular frequency of the modulation and So is the
steady component of the source present at z = 0. If the neutron flux is
assumed to be separable in space and time, the oscillating part of the flux
can be represented as (1)
Q(r t) = Re[r (1.5)
Equation (1.5) may be introduced into Eq. (1.1), resulting in
[ + 6Cz (r) 3 Q, v (r) (1.6)
or, in slightly different notation,
D Q = iw L(x,z) e(r) (1.7)
where
a a vCa (1.8)
Do vD (1.9)
GL(x,z) a 6Ca(r)/D .
(1.10)
8
The boundary conditions will be taken to be
m(r) )
x=+a
= i (r) = 0
zCD
(1.11)
In addition, at the source plane,
iwt
e
(1.12)
ae(r t) 1
D = Re [S(r) I
dD z 1 2
z=o z=o
Using Eq. (1.5), this reduces to
ae(r) 1
 D =S(r) .
az 2
z=o z=o
(1.13)
The spatial flux, Q(r), may be represented by an expansion in the
complete set of transverse eigenfunctions defined by
(9 + B a) 4 m(x1y) = 0
and by the boundary conditions at x = +a and y = +b plus the conditions of
symmetry around x = 0 and y = 0. Normalized, these eigenfunctions and their
related eigenvalues are
1 (2R 1)nx (2m 1)xJy
4 m(x'y) = cos cos
Am~ 2a 2b
(1.15)
and
(1.16)
wJhere R and m may take on all integral values.
= (r)
y=+b
(1.14)
(20 i+ 1 (2m 1)Jrn
The expansion of Q(r) is given by
Q(r) = Qm(z) bQm(x,y) .
,m
(1.17)
Substituting Eq. (1.17) into Eq. (1.7) gives
al + iw
Do i im,(x,y) Q~m(z) = BL(x,z) Q~m(x'y) ~Qm(z) .(.8
a~m a~m
dy (xy) esutsin
b
a
a
a,m
4 m(z) ,
(1.19)
where
a, + iw
Do
p B
pq9 1pq
(1.20)
and
BLpqim(z)' S ax dy, *F (x~y) 6L(x~z) Qm(x,y)
a b
(1.21)
Before continuing further with Eq. (1.19), note that the two boundary
conditions in z have yet to be applied. In order to make use of the source
boundary condition, S(r) must be expanded in the same transverse eigenfunc
tions. That is,
Am
S(r)
z=o
(1.22)
whe re
a
Sem; =ax
a
"m(x'y).
It should be pointed out that this expansion, as well as the flux expansion,
assumes symmetry around the x and y axes. If one uses a source containing
asymmnetric components, oQm(x~y) must also include sine tents. Inserting the
expansions, Eqs. (1.17) and (1.22), into the source boundary condition,
Eq. (1.13), gives
de (z) S
am am
 D (1.24)
dz 2
z=o
The boundary condition at z = OD becomes
lim
z me(z) = O (1.25)
reduction of Eq. (1.19), substitute Eqs. (1.3)
The y integration may be performed immediately,
delta. Also performing the x integration gives
Now returning to the
anid (1.10) into Eq. (1.21).
yielding 8qm, the Kronecker
the result
M
6L (z) =h 8 8Fp 82~~ 8(z z )
pqim D qm p ~
k=1
where 6Fp is given by
b
Sdy S(r) Z=
b
Now define the Green's function G q(z,S), satisfying
ac; 2~ pqpq(ZS D_( z
(1.32)
and the homogeneous forms of the boundary conditions satisfied by 4 q(z).
Multiply Eq. (1.50) by G q(z,S) and Eq. (1.92) by Op (5). Subtract the
co
two results and operate with d6 to obtain
d i
o
(G(z,S) ddg(5) 8(5 G(2cc ) d"i F(i) G(zi) = ez)
(1.55)
The subscripts have been dropped temporarily for economy of notation.
Integrating the first integral by parts results in
m(z) = D (z,S) d5 g) ) jGzg
oo 00
, D dg az, d( + D S idi(~cz~
oo
 D d i F(i) G(z,i)
O
(1.34)
The first two
gradient of G
integrals cancel, both + and G go to zero as (  c, and the
is zero at 5 = 0. Therefore, the desired integral equation
 D dj F(i) G(z,I)
O
(1.55)
Q(z) = D G(z,g)
.(a + p 1)rxo
sin ~
_iT3~S~a
(a p)nxo
sin 
a
pl= (a p)n t
(1.27)
(2 1)nx,
a+ (2a 1)Jr
In order to simplify Eq. (1.27) further, let us assume that xo is suffi
ciently less than a that the first term of the Maclaurin's series may be
used for the sine terms. Equation (1.27) then becomes
2x
Fpa ao
(1.28)
Sfor all a and p.
Noew inserting Eq. (1.26) and Eq. (1.28) into Eq. (1.19) gives
BI 6$r(z z) 1L (z).
k=1
It may be noted from the above equation that the onedimensional problem
turns out to be not really one dimensional, in the sense that coupling re
mains between the spatial harmonics associated with the x direction.
In order to solve Eq. (1.29), it is convenient to convert it into an
integral equation. Substitute 5 for z and rewrite as
(130)
where
2hx
Fpq() Da akk
L k=1
d~ o
(~2 i~~~ bpa(i)= %$i)
The problem now is to determine the Green's function defined by
Eq. (1.32). The solution to the homogeneous form of Eq. (1.32) is
G(S,z) = Aepz + Cepz
(1.36)
where 5 and z have been switched. The particular integral of Eq. (1.32) may
be found by the technique of the variation of parameters. Let the constants
in Eq. (1.56) depend.on z (they also depend on 5, of course, but this may
be ignored at this point):
G(S,z) = A(z) epz + C(z) epz .
Differentiation of Eq. (1.37) gives
(1.57)
aG(S,z) pz pz dA pz dC pz
= p Ae + p Ce + e + e
(1.48)
Choose A and C such that
dA pz dC pz
e + e = O 3
dz dz
Then differentiate again:
a~(Sz) = pZ Ae~P + Cp e z
az,
(1.59)
dA pz dC pz
p e + p e
dz dz
(1.40)
Inserting Eqs. (1.40) and (1.57) into the original different
for the Green's function results in
A ~p epz gepz + C p epz epz
Apz dC pzp 6(z g)
dzdz D
ial equation
(1.41)
pz dA pz dC 6(z 5)
pe + pe
dz dz D
(1.42)
Equations (1.59) and (1.42) are two equations in dA/dz and dC/dz which may
be solved simultaneously to give
dA 8(z g) epz
dz 2pD
and
dC 8(z 5) epz
dz 2pD
Integrating Eqs. (1.4)) and (1.44) gives
A = e S/2pD + B
and
C = e P/2pD + E
(1.4~)
where B and E are constants to be determined from the boundary conditions.
Using Eqs. (1.45) and (1.46) in Eq. (1.56) gives
G(S~) p pt p(gz) .p(Sz)
G~~)=Bep + Eep e (1.47)
2pD 2pD
or
G(5,z)pz pz sinh p(S z)
G~~)=Be + Ee + (1.48)
pD
The onedimensional Green's function has a discontinuity in slope at the
source point, z = E. Actually, this does not need to be known beforehand,
since straightforward application of the only two boundary conditions now
available leads to the solution G 0 Then, in order to obtain a solution,
before applying the boundary conditions, Eq. (1.48) must be divided into
two parts, one which will apply for z < ( and the other for z > E. Now,
two more boundary conditions are needed, which may be found by integrating
Eq. (1.32) over a small region including the source point:
dzP p (5,z) = 6(zD .5 (1.49)
Assuming that the second derivative is more singular than the function it
self, one obtains, in the limit of small E,
aG(g,z) aG(5,z) 1 (.0
az az D
Therefore the two boundary conditions are continuity of G at z = ( and
Eq. (1.50), which gives the magnitude of the discontinuity in the first
derivative.
The two parts of the solution will be denoted as
L L pz L epz sinh p(E z) (1
G (g~z)= B e E e + pD ,z<((.1
and
GR~eZ =R ,pz R epz sinh p(5 z) (.2
G (g~z)= B e E e + pD ,z>(.(.2
The z = 0 boundary condition is
dGL(S,z)
dz
=pL + pEL cosh pi O
D
(1.55)
z=o
At large z, the boundary condition is
lim cR~g,z) =o lim E epz ep ep ]=
z oz o 2pD
Equation (1.54) immediately determines ER:
2pD
The continuity condition, at z = 5, leads to
BLeP + EL pS = BReP + ERe .
Equation (1.50) becomes
(1.54)
(1.55)
(1.56)
 pELeP + 1
= 
 pBReP + pEReP + pBLe P
(1.57)
Using Eq. (1.55)1 Eqs. (1.53), (1.56), and (1.57) may
be put in matrix form
as
,sh pS
D
2 D
2D
co
p
e" 5
pg
pe
(1.58)
The three solutions to Eq. (1.58) are
p O B
eP e'P EL
pS p5 R
pe pe B
BL sinh pS (1.59)
pD
EL epD (1.60)
pD
BR  (1.61)
2pD
Upon substituting Eqs. (1.55), (159), (1.60), and (1.61) into Eqs. (1.51)
and (1.52) and rearranging, one obtains
3L~gz = e Dosh pz (1.62)
and
pz
cR~g,z) = e Dosh pS z > 5 (1.63)
It may be noticed that Eqs. (1.62) and (1.63) define a function
which is symmetric with respect to interchange of the coordinates of the
source and observer locations (z and 5), thus satisfying one of the nec
essary criteria for a Green's function (9). Equation (1.63) may now be
substituted into Eq. (1.55), giving
.pz d pz
Q(z) = D O~ D F(S) dS
pD dS pD
g=0 0
D e cosh pz Fi) dj (1.6Li)
Reinserting subscripts and using Eq. (1.24), one obtains
pz
Pq p z z~ 5
S e pq
Opq(Z s 2p Dp pq pQ
pq pq
at e (5) (1.65)
pq
Since F q(h) contains 8 functions in j, Eq. (1.j1) may be inserted into
Eq. (1.65) and the integrals in 5 performed, leading to
S e Pq 2hx 
+ z)=ps o eB ohp z4(
mpq( 2p D p Da~ 9 a 6~ k Cs pq k* eq k
pq pq k=1
+ csh Pp, z i, B6a,k ep$ zk p + zk) i (1.66)
k=Mr+1
where M' is determined by zM' < z < zM'+1. If z equals one of the zk's,
it can be taken to be either zM' or zM'+1l since the Green's function is
continuous at the source point.
In order to obtain a solution to this system, the original expansion
of Q(r), Eq. (1.17), must be truncated. Assume that the desired a~pproxi
maztion is
*(r) ~V ZL (z 4I "(xy) (1.67)
=1 9=1
For each value of q, a set of LM simultaneous equations may be obtained
from Eq. (1.66) by writing it for each of the L modes and for z = z ,
z ,...zM.This set determines each of the LM values of 4 z
Writing these values as a vector, where the superscript T indicates the
transpose,
(1.68)
(1.69)
where [I] is the identity matrix, the vector B_ contains the constant
(source) terms, and the matrix [A ] contains the Green's functions. B
Q q
may be written as
Bg [il ~ lqZI l lq 2 Y
(1.70)
whe re
S
aq b2 p D
The [A ] matrix may be partioned into L' submatrices which are each M by
M. The n~p element of the I,J submatrix is given by
(A ~)n~ = Bp (z ) cosh P~lp zeIQn
(1.72)
Sn > p
(A ~)n~ = Bp (z ) cosh pgzn e lq p
(1.73)
, n
where
Og = $1(z,) m1q(z2) '. lg(zM) 2qg(z1' .. mLq(sM)
the set of equations may be put into matrix form. That is,
iAq q + [I] m = B ,
i~ Cq Pgz#
(1.71)
2hx 6C
B, (z )' ,LEoa (1.74)
Solving Eq. (1.69) Q times by inverting the matrix [A ] + [I] gives
the 0 s:
g [A ] + [Il]l B (1.75)
which may be inserted into the summations in Eq. (1.66) The results of
Eq. (1.66) are then inserted into Eq. (1.67) to give the desired approxima
tion to the spatial flux *(r) .
CHAPTER II
USE OF A ONEGROUP TWODIENESIONAL GREEN'S
FUNCTION IN HETEROGENEOUS NEDIA
The problem of the previous chapter will now be extended to two
dimensions by allowing the perturbations to be placed anywhere in the xz
plane contained in the moderating material. Equation (1.1) will still apply,
but now BZa(r) will be taken to be a general function of x and z. Equation
(1.7), with its boundary conditions given in Eqs. (1.11) and (1.13), will
be taken as the starting point, leaving 8L(x,z) general.
In this case the spatial flux Q(r) will be expanded in the single
set of transverse eigenfunctions, defined by
 + B 4 ~(y) = 0 (2.1)
dy2
and the boundary conditions 4 (+b) = 0 and 4m(Y m (y). The eigenfunc
tions are
and the eigenvalues are
[~m (2m 1)T(b(2n
where m = 1,2,3,... Using the following expansions of the flux and the
source,
and
S(r) = Sm(x) em(Y) ,
in Eqs. (1.7) and (1.13) gives
(2.4)
(2.5)
o+ iw
C~ o
m
4 (x, z) 4 (y) = 8L(x, z)
m(x, z) (y)
(2.6)
. D 4m(x, z)
m
1 Z
m
Sm(x) em(y)
(2.7)
b
with day mn(y) gives
b
Operating
O
+~ 
fz n
4 (x, z) = GL(x, z) + (x, z)
(2.8)
whe re
2 ao + iW
n n Do
(2.9
and the corresponding boundary condition
1
2 Sn(
(2.10)
D a
z= o
m,(y) I
z= o
The remaining boundary conditions given in Eq. (1.11) still apply.
It is evident that this system may be treated similarly to the prev
ious problem if a twodimensional Green's function can be found which is
the solution of
+2 aZ 2 n1(x z;xo ,z ) =8 x x) 6(z z ) (2.11)
with boundary conditions
Gnt (+a3c,z;x ,z )= n _li~m G =O .(2.12)
z=o
This Green's function may be found by expanding in the eigenfunctions of
d + k ((x) = 03 #k(+a)= O .(2.13)
These are
(2 1)Jrx anx (.4
Jr(x) = AR cos + BQ sin (.4
Note that k assumes two different values for the two different eigenfunc
tions, the sine and the cosine. The expansion of Gn(x,z;x ,zo) is
[ (2& 1)Jx axx
On(Xjx,zx,z Z) = C (z'zo ) [A (x cos 2a + B (xo) sin a
a (2.15)
Inserting Eq. (2.15) into Eq. (2.11) gives, with arguments dropped for con
venience,
2 2
cos B sin n
I ~2a 2a a
a~e o (2 1)lrx nx
+Z A2 os + BQ sin  ~ p C
a z a
['i O (2 1)J[x inrx] x )hZ 1 (.
cs 2a + BQ sin =6xx)8zz).(.6
a
asi 0^E a gives
a
C.A. 2j 1)n a +1J Aap C.Aja = 8(z z )cos (2 )x
2a8a a 2a
(2.17)
and
 C.B. a + J Bja p C.BJ a = 6(z z ) sin  (2.18)
JJ 8a n a
These two expressions may be rewritten so as to have the lefthand sides
functions of z and zo only and the righthand sides functions of xo only.
They become
a1_ rL 2 f (2j 1)x CO cos no
za2 Ln 2a 2a
= (2.19)
8(z zo) a A.
and
a n j sin 
bz a .(2.20)
6(z zo) a B
Obviously the C.'s in the above two equations are not the same; so the C.
will be replaced by A! in Eq. (2.19) and by B! in Eq. (2.20). Notice also
that both sides of Eqs. (2.19) and (2.20) are equal to constants and no
loss in generality occurs if this constant is taken to be unity in each
equation. With this step, A. and B. may be obtained immediately:
1 (2ij 1)nxo
A.(x ) =cos (2.21)
o a 2a
and
B.x) 1 sin o (2.22)
B. a) a
and with the definitions
42 p (2j 1)n (2.23)
4nd
a p + (2.24)
the two equations for Aj and B:J are
~2 2 A ;(z,zo) = 6(z z ) (2.25)
and
(2.26)
10m
SA! =
zKoo 3
lim
BI = 0
(2.27)
z=o
z=o
plus the function continuities and firstderivative discontinuities at
z = zo. These functions are virtually identical to the onedimensional
Green's function developed in Chapter I and can therefore be immediately
written down, being
e cosh uz
A (z,zo) =
e cosh Ciz
and
e zoshg
Sz < zO
S z > z
(2.28)
>z < zO
e~ cosh gzo
(2.29)
Sz > z
Collecting Eqs. (2.21), (2.22), (2.28), and (2.29) into Eq. (2.15), the
complete Green's function is
a2 B!(i.zoz) = B(z zo) .
The related boundary conditions are
zAI z
pro (26 1)nxo 2 ~t
G x z~oz Zre cosh CLzc o (s 1x
Qnxzxoz L IJa os 2a os 2a
rlzo asrx
e cosh rlz o .rr xx2.0
+ ~~~sin  sin  z
IIa a ao
and
C ~ o (P l o (2 1)rxx
G (~z~ z =cos cos
n o o a 2a 2a
e~ cosh gzo axxo alxl
+ sin  sin  I z > z (2.31)
rla a aJ o
Similarly to the procedure followed previously, this Green's function
may be used to formulate an integral equation for en(x,z). xo and zo are
interchanged with x and z in Eqs. (2.8) and (2.11). The difference between
Gn(xo,z 3x,z) times Eq. (2.8) and @n(xo'zo) times (2.11) is taken. The
oo a
resut i opeate on ith dz d ogvwt ruet n ne
o a
gration limits dropped,
dz~ ~ Sx~ Gn z a o dzoJ dx e _
o o n a2z, o o n n
In the terms containing partial with respect to x, an integration by parts
in x may be done, and similarly with the z partial. All remaining surface
integrals on the lefthand side of Eq. (2.32) cancel and most of the line
integrals are eliminated by the boundary conditions. The remaining terms
give
a as(xo,,zo)
@n(x'z) = dx, Oncx ,ZO;x,z) n zo
a zo=o zo=o
Co a
dzg axo Gn(xo.z 3x.z) BL(x ,z ) "n(x ,z ) (2.53)
o a
Equations (2.10), (2.30), and (2.31) are to be substituted into Eq. (2.53),
along with the desired form of BL(xo~zo), in order to obtain the integral
equation for On(x)z).
If the fons of BL(x,z) of Chapter I is used in Eq. (2.53), it is
straightforward although rather tedious to show that Eq. (2.33) reduces to
Eq. (1.66). To do this, the following expansions must be used:
d (xIz) = e mj(z) O.(x) (2.34)
and
S (x) =Z Sn Q.(x) (2.35)
where
mjx) 1 C (2j 1)rxx
+ x o 2a (2.36)
29
These expansions are justified, since if 6L is symmetric in x [and it is
assumed as before that S(r) is symmetric around the x and y axes ], then
en_ will also be symmetric in x. It is then a matter of performing the xo
1 ( 2p 1) stx
and zo integration in Eq. (2.jj) and operating with Ja cos 2a d
a
to obtain Eq. (1.66) .
CHAPTER III
USE OF AN AGEDIFFUSIGN TWODIMENSIGNAL GREEN'S
FUNCTION IN HETEROGENEOUS MEDIA
The development of this chapter is similar to that of the preceding
chapter except that the theory is to be the continuous slowingdown model
with a thermal group. The geometry is the same semiinfinite moderating
medium, with the perturbations now being small fuel rods which, as before,
run the full length of the y axis. This model is similar to the Feinberg
Galanin treatment of heterogeneous reactors (7) and (8), with two major
exceptions. In this study time dependence is included and the finiteness
of the medium is taken into account. This will result in a finite medium
Fermiage kernel, or Green's function, as well as the corresponding diffu
sion kernel obtained in the previous developments of this work.
For lethargies below thermal, the applicable equation is
D(u)~ o2 (ru,t) + Ca(u) e(r,u,t)
1 ae(r,u, t) aq(r u, t)
=T; + S(r~ut) u '(3.1)
where the nomenclature is similar to that previously used except that now
all quantities are lethargydependent and S and q represent the neutron
source and slowingdown density, respectively. It has been assumed that
Za(u) for u < us, where us is the thermal neutron lethargy, is independent
of position. This is valid as long as the fuel rods are small.
The relation between the flux and slowingdown density applicable
to an infinite medium will be assumed, that is, that
9(r,u,t) = 5 It(u) Q(r,u,t) (3.2)
where g = the average logarithmic energy decrement, and Ct(u) = the total
macroscopic cross section in the moderating media.
The thermal group is described as before, except for the addition
of the source of thermal neutrons from the slowingdown density. This
relation is
[D +E + 1 a ra)= t~ ,t) A(r, t) ,(3.5)
where A(r,t) represents the absorption in the fuel rods and the subscript
s indicates that the quantity is evaluated at the thermal neutron energy.
In order to define the quantities S(r,u,t) and A(r,t), it will be
assumed that (1) all fission neutrons are born at lethargy u = O, (2) all
fuel rods are lines, the kth rod being representable as 8(x xk) 8(z zk)
and (3) the driving source is located at z = O, producing only thermal neu
trons and is of the form Re[S(x,y) 6(z) e in. As before, since the source
is present only at the boundary z = O, it will be included as a boundary
condition rather than being included in Eq. (3.5).
Since this treatment is an extension of the FeinbergGalanin treat
ment, some of Feinberg's notation will be adopted. Define the thermal
constant 7k as the ratio of the total net current of thermal neutrons into
the kth rod to the value of the thermal neutron flux at the surface of the
kfth rod. 4 (r,,'t) will be used for the thermal neutron flux at the
surface of the kth rod. The neutron yield per thermal neutron absorption,
r7, will be assumed to be the same for all rods. With the above assumptions
and definitions, S(r,u,t) and A(ryt) become
s~~~r,~ ~ ut ~)(0)
s~~~r,~ u, t)= e(u k ms(rk t) 6( r ) rk.
and
A(r t) = ksrt) 6(r rk) .(5)
IJ ~k
The ratio of the scattering to total cross sections at the source lethargy
in Eq. (3.4) represents the probability of a neutron starting the slowing
down process.
Due to the separability of space and time for the flux, slowing
down density, and source term, Eqs. ( .1), (3.2), and (3.4) may be combined
to give
[ D(u) + Zalu) + + i LT(u) q (r4u)
2 (0) Z Bu k ( ( 36
= C t(u) "ppy n ~ )7 s rk rr
Similarly, Eqs. (3.5) and (3.5) become
[ D + 2 E e(r) = (r~u U) 7k 6 (r ~l) .(>.7)
The boundary conditions on es(r) will be the same as those used for the
thermal neutron flux previously; the boundary conditions on q(r,u) will
be identical, for u > 0, except that at z = 0, q(r,u) will be taken as
zero on all boundaries. This corresponds to assuming that Q(r~u) = 0 at
all boundaries for u > O, which is consistent with the thermal neutron
flux assumptions.
Defining
iw
Eq. (3.6) may be written as
[ 5t(u) o' + B(u,w) + au ]4(rP)U
M~ultiplying through Eq. (3.9) by the integrating factor
exp B(u',u) du ], Eq. (3.9) may be written as
DZ u) ~ q(r,u) eo
rIB(u ',w)du'
o
I (0)
cosO/ b u) 7k ~s ~k) r gIk)
(3.10)
In order to simplify the notation, define the Fermi age, 1, from
dz D(u) with z(u = 0) = 0
du 5( I~u
(3.11)
and also define a frequencydependent resonance escape probability,
p(u,W), as
p(u,w) =eo
(3.12)
Notice that the above definition is similar to the usual resonance escape
probability except that an additional 1/v absorber has been added by the
oscillating source.
With the two definitions [Eqs. (3.11) and (3.12)] Eq. (3.10) be
comes
I C(0) (r 6 )
1 s~, CtOJ k rk)
(3.15)
Bo 'm~u
+ fru)e
The partial with respect to the Fermi age will be eliminated by
using the Laplace transform. Defining
6(rysa e dr (3.14)
Eq. (3.13) becomes
ye (rs) s 8(rys) +FijT 0 k1 s k r
k (3.15)
The term q(r,0) is zero since q(r,0) was incorporated into the original
equation and is, of course, the negative of the righthand side of Eq.
(3.15). If it had been chosen that this source be introduced as a boundary
condition, it could have been omitted from the original equation and have
been introduced naturally at this point.
The obvious step at this point, consistent with previous methods,
is to expand 9 in the y eigenfunctions given in Eq. (2.2), which were
determined by Eq. (2.1) and the boundary conditions of symmetry and
k(+b) = 0. The required expansions are
e(r,s) = e,(x,z, s) 9,R(Y) (3.16)
and
For future use, it should be noted that Eq. (3.16) implies that
Eq. (3.16) may, in fact, be obtained from Eq. (3.18) by Laplacetransform
ing both sides of Eq. (3.18), where QQ is the inverse transform of 9 .
Inserting Eqs. (3.16) and (3.17) into Eq. (3.15) and operating with
(y)I~ dy gives
Ox
+ 
a2z
a (x'z~S)
Z (0)
= 4 )
7k 6(x xk) 6(z zk
Q(x~z)
(3.19)
whe re
[2(2a 1)n]
S2b
(3.20)
Now 6 (x,z's) is expanded in the x eigenfunctions of Eq. (2.14),
defined in Eq. (2.13):
6 (x'z's); =B 6(z,s) 1 cos
J
6!(z,s) 1 sin j x
Ja a
(2j 1)nx +
2a
(3.21)
Inserting Eq. (3.21) into Eq. (3.19) results in
S(s + B )
in
a
6 j
j
C (2j 1) 1Xax + B
2a j
1
sin j[
a
a8z
a8z
+
cs(2j 1) ax
via 2a
1 sin jC
,A a
1 (2j 1)nx
 0 ) . COS L +
Ja 2a
J
7Ex
a
8'  sin
x (0)
= 0) r
k, 6(x xk) 8(z zk) Q(x'z)
(3.22)
a
S1
(2m 1)71x
cos 2adx gives
Operating on Eq. (3.22) with
a82 z ,s)
7k 8(z zk)
 (s+ B ) 6( z~s)
(2j 1)nxk
2a
1
cos C
Q(xk z)
(3.25)
a
1 .mnx
 sin  dx gives
a
Similarly, operating on Eq. (3.22) with
a8 9(z,s)
+ (s + B ) 9!(z~s)
a2 J
a [
 6!(z,s)
3
Is(0)
= z (0
1 anxk
sin  (kz)
Z k B(z k)
(3.24)
(2j 1 a
2a
 8(z,s) (2j 2al)n r2a
Defining
2 as + B2 + (2j 2al)n (3.25)
and
52 's + B~ + ,a (4.26)
Eqs. (g.2g) and (g.24) may be rewritten as
a26.(z,s)
J 2 B(0~)1 (2 1)x
= I Ok 6(z zk) cos 2a (xz) (3.27)
and
a8 9(z s)
Jr g 9(z,s)
Zs(0) Z kh 8(z zk) J1 sin k +(xk'z) .(3.28)
Ska
The variable s has now been hidden in CI and g; so the above two equa
tions are effectively one dimensional. It may also be noted that Eqs. (3.27)
and (g.28) are Green's function equations. Consider the boundary conditions
which are appropriate for these equations. Aside from the usual conditions
of continuity at z = zk and the discontinuity in the first derivative at
that point, the appropriate conditions must be derived from the original
z acis condition on q(r,T), which are
z=o
Using these relations along with the expansions that have been utilized,
Eqs. (3.16), (3.18), and (3.21), results in the following boundary condi
tions for Eqs. (3.27) and (3.28):
e.(z,s) = lim 9(~)=0(.0
z=o
and
9!(z,s) = lim 9(~)=0.(.1
z=o
Except for the fact that the function itself, rather than the first
derivative, is zero at z = 0, we now have a problem identical to the one
dimensional Green's functions of Chapters I and II. By following procedures
identical to those used in Chapter I, it is easily shown that the solution
to
d29(z,zo)
K (z,zto) = (z zo) (s.52)
dz2
with boundary conditions
6(z,zo) = im 6(z,zo) = 0 (3.SS)
z=o
K2
B~zzo) e sinh K zo
Kzo
e sinh K z
Sz > zo
(g.~4)
Sz < zo
Comparing Eqs. (2.25) and (2.28) with Eqs. (.3.2) and (3.54) shows that the
change in boundary condition has simply resulted in the cosh being replaced
by the sinh. By writing the sinh as the difference of two exponentials in
Eq. (3.54), it is easily seen that the two solutions may be combined into
a form which will prove to be useful later:
K zzo( K(z+zo)
e(z,zo) =e 2K 2K (3.55)
Using this result, the solutions to Eqs. (g.27) and (g.28) may now
be immediately written down, being
Zs(0) 11 p zzk
6 (z, s) = e~ ~
1(2j 1)rxk
7k~ cos 2a e(xk zk)
(~.36)
CI( Es(0)t g zz ['"'
6 (z, s) = e
1 Exk
7k sin  + ( zk)
(~~7)
p";z"zk
S(z+zkl
41
Inserting Eqs. (3.56) and (3.37) into Eq. (3.21) gives
1J(z+zki
e
c (0)
sif
j k
e (x'z's) =
(2j 1)nxx
cos 2a
c (2j 1)~xx
cos 2a
R(xk'zk)
g(z+z )
ax
(3.38)
sin j~Xsin
The next step is to inverseLaplace transform Eq. (3.38), obtaining
the expansion coefficient in Eq. (3.18), which is
Q (x,z,,w = co
6 (x'z's) .
(3.59)
Rewriting the definitions of y and 5 as
I +B
and 52 = s + B' ,
(3.40)
whe re
and B 9 =
(3.41)
it is evident that each term to be inversetransformed is of the form
e~ t (3.42)
lJ zzk
Yk e
 Izz l
m (xk, zk)
B = (2c~ 2bl7n (2j 2al)x 2
22bl ~ a
which has the inverse
 ez24r for zu > O .
Using Ea_. (3.41), the inverse transform of Eq. (3.38), Q (x,z,r), is
(}.43)
Zs(0)
= t
(z+zk)24
2a /r Yk k
(xzw)
cos cos + e
2a 2a
"F
sin j sin
a
4 (xk)zk)
(~.44)
Returning now to the thermalneutron diffusion equation, insert Eqs.
(3.17) and (73.18) [with Ts replacing r in Eq. (3.18)] into Eq. (j3.7)
and
b
b
(jy) dy to obtain
 Ds
a2
a22
Ds B as v
= p(. s (,z,,Ts Ik
k
,(x,z) 6(x xk) 6(z zk)
.(3.45)
Using the expansions
q~~)= z cs(j2alx + z i a
I,~lr ~JQZ o 2 n j,) i x
(g.46)
(zzk) /AT
BE.T
e 3
3m(x,z)
and
j~(zTs) cos (2j 2al)rrx + Q(z,ts sin j~
1
Z
(3.47)
Q (x, z,s 1
a
cos (2a l~ dx gives
in Eq. (3.45) and operating with
D a+
Sdz
+ 
as v
s
(2j 1) xx
Yk cos 2a k
4 (xk~z) S(z 
zk).58
= P( s' j zTs
a
with sin dx gives
aa
Similarly expanding and operating
a 4' (z)
D j
Sdz
[D B~ + Ea
]m' (z)
+ 
v
s
Yk sin a
k
(3.49)
Q(xk~z) 6(z zk)
=p(Ts ,W
Defining
C +
K2 B + s
3 D
(350)
Ds 2aj(2j 1)Jl ]2 I(z)
i j(z)
[D B~ + E
Dsa "j
Z + 
K' ,B' + s
j D
Eqs. (3.48) and (3.49) may be written as
dze
p(r ,w)
 K 4.(z) = s Q(,
0Ds &jts)
Z (2j 1)nxk
7k cos 2a a(xk'z) 8(z zk)
1
+ 
D/
(3.52)
and
dze
 K'2 4! (z) = Q(zT )
03 D~s s
Yk sin 4 ,(x ,z) B(z z )
+D
Comparison of the expansion, Eq. (3.47), with Eq.
(3.44) shows that
(zzk)
Zlk~~km l
(z+zk)2
4r
 e s
s
Z (0)
Q (zTs) c o
(2lj 1)nxk
cos 2a a(xkrzk)
(5.54)
and
(zzk)
s
(z+zk)2
s
e
Z (0) js
Q (z, Ts 2 k k
sin a~ (xkfzk)* (s55)
Since Eq. (3.52) along with Eq. (3.Sk) is so similar to Eqs. (3.SS) and
(555), attention will be restricted to the former set; then the latter
set will be handled by comparison.
Equations (3.52) and (3.Sk) may be combined as
d2 ) K 9 z A Fjr(Lz, z) + ZHjk(z) 6(z zki (.6
k k
whe re
Z (0) 71 p(T ,W) e a
A ss,(357)
aS
(zzk)
Fjki(zzk) lk e [ "s
(z+zk)2
s
e
cos a R(xk~zk '
(3.58)
4 (x z). (3.59)
1 (2j 1)nxk
Hjk~ D~ kgC 2a
Equation (3.56) is virtually identical to Eq. (1.30) and has es
sentially the same boundary conditions, so that the Green's function of
Chapter I may be used to convert Eq. (3.56) to an integral equation. The
appropriate Green's function is determined by
S Ke G.(z,S) = 8(z 5) (3.60)
and the usual homogeneous boundary conditions. The solution to Eq. (3.60),
from previous results, is
Kz
e'K cosh Kz
= K, z<5 (3.61)
As before, the integral equation for mj (z) is found by forming
oo co
Sdl G (g,z) [Eq. (3.56)] di 4 l(S) [Eq. ().60)], after reversing
O O
z and 5 in Eqs. (3.56) and (3.60). This result is
= J di G (Llz) Fjk(ilzk)
o k
+ d i G (i~z) Hjk(i) B(L zk Ojl(z) (3.62)
o k
An integration by parts may be performed on the integrals containing the
derivatives and the integral containing Hjk may be completed, resulting in
md (z) = G (zk'z) Hjk(zk) + Aij di G (i'z) Fjk(iizk)
k o k
a G. (4,z) o
(3.63)
Inserting the boundary conditions into the last term gives
oo
'j(z) = G (zk'z) Hjk(zk) + AQ at G (i,z)
Fjk(S'zk)
gdo [=o
(3.64)
Eq. (3.64) becomes
Inserting A and H.k and assuming that z < z < z 1
M'
k=1
do. (5)
dg
5=o
cosh Kzk 7 cos 2
Kz
e
K
1 z
KD Ja~
'^ (z)
k= M '+1
(2j 1)rxx
7k cos 2a
(x,zk,
4 (xk'zk) + cosh Kz
B ~s
e
o k
Z (0) q p(T ,w)
iKz
Fjk(g~zk)
Lt(0) 2! K Ds a t
(3.65)
+ cosh Kz
d 4. ()
SdgeK Z Fjk(l,zk~
z k
48
It seems appropriate at this point, before continuing into the
mathematical complexities, to simply the notation in Eq. (3.65) and try
to understand its physical significance. To this end call the first term
in Eq. (j.65) H Q, since it represents the contribution which would be
present in the homogeneous case. Also, redefine Eq. (3.58) and (3.59) as
Fjk(z,zk) =ar I(z) W.jk ,xkzk (.6
and
Hjkc(z) = ~Wjk (xk~zk)(.7
where
i(z) e~ [" '' e (3.68)
and
7k (2j 1)nxk
w.k =cos (3.69)
Jk 2a
W~jk is simply a weighting factor depending on the distance f~o~m the center
line of the foil location. N (z) is the slowingdown kernel for the parti
cular geometry chosen, having an obvious resemblance to the infinite:medium
Fermiage kernel. Noticing the factor 1/ /T, it can be seen that the
slowingdown kernel obtained for the present geometry resembles the plane
source infinitemedium slowingdown kernel most closely. This is to be
expeccted, since the problem was reduced to a semiinfinite, onedimensional
case.
Inserting Eqs. (3.66) and (3.67) into Eq. (3.65) and using the
symmetric part of Eq. (3.46) and Eq. (3.61) results in
4 (xz) = W (x H   Wjk G (zk'z) 4 (xk'zk)
j j k
oo
+ Z AJ;J W (x) Wjk Q(xk,zk) S di Mk~(i) G (zk 5) (3.TO)
j k o
where
1 o (2j 1)nx (1
w.(x)= cs(.1
~ 2a
It should now be noted that the Ath mode of the flux distribution consists
of the following parts: (a) the homogeneous term H e, which represents
the propagation of a neutron wave in a homogeneous system, (b) the absorb
ing part, the second term, which accounts for the heterogeneous absorption
of the plates, and (c) the last term, which accounts for the production of
fast neutrons by fissions in the plates, and their subsequent slowing down
into the thermal group.
Returning to the mathematical manipulations, it now remains to
evaluate the two integrals in Eq. (3.65), which are complicated by the
fact that K is a complex quantity. The first integral may be written as
(5zk)2 (S+zk)2
s S
d'i cosh Kg e  e Ejk (j.72)
o k
whe re
jB~ ~ (2 T C ~ 1)nxk) 1 3
Expressing the cosh as the sum of exponentials and performing the squares
in the exponentials, Eq. (3.72) may be reformulated as
Z'2
k
Sd KEb i Kg s+Yr
d( e +
gzk/2Te s
gzk/2T)Ej
Multiplying and combining the various exponentials leads to
zk 6 rs z
S2 Ejk Sd
k o
[K(z /21s )]5
e
[K(zk/21s)](
(3.75)
~lK+(zki/2r1si
By completing the squares in each of the exponentials, Eq. (3.75) may be
rewritten as
(3.74)
i'T~ [[K+(z /21 )]g
s~K 2r k
s2r
s
2k 4s Z
e 2 Ejk S
o
dg (e s
 J~ K
T( Z
1 K
s 2s/
T ~ Zk
1 K
2
SJr~
2
S
 ~r
Zk
2s
I K+
+ / K+ s
(3.76)
Each of the four integrals in Eq. (3.76) is of the form
r a2 z T
which becomes, upon changing the variable of integration to
(377)
v l + Fa ,
2 s
(3.78)
+~ / K
( 2s
t 21 Is
JIS
+ Jftr
e "~ dw
 + J ex
S
O
 S

e dw
O
(3.79)
2 ;I
~2
e dw 
Recognizing the two integrals as error functions, Eq. (3.79) becomes
es
+J; / erf
 + Jr

S
es
(3.80)
where the error function is defined as
x
erf(x) 2_ s
t
(3.81)
dt .
The integrals of Eq. (3.76) now become
(erf [
SfT a
,] + erf [J t a
k Eak 2 s
K +
Zk
s
2
s
2s
 erf 
s
(3.82)
The second integral in Eq. ( .65) is
n
z
di e e ~(g24
(i"*'""Ejk
(3.83)
e
which can be handled in a similar manner to give
2
K z s
s es
2
SEjk e
k
K
(3.84)
s21c
s
t e5 Z
Z
(erf 2z_ F
2 s
erf 2 c + / K+k rfK+
[ s s
s s
(erf z+[
S
K 1 2 erf J. K 
s
1 rf z +
s
e s (1 erf + K +
s
If ~ ~ ~ ~ ~ ~ ~ ~ ~ th qaesi h fcosex in Eqs. (j.82) and
(1.84) are performed, a factor exp [zk s] will result which will cancel
the exp [ zk /4s The two remaining factors are exp [+ K zk] and
exp [Ts iE i which can be taken out of the summation and combined with
the factor ex~ [ Bj s '
Now the expressions for the two integrals, Eqs. (3.82) and (3.84),
may be inserted into Eq. (3.65) to give the complete expression for ". (z).
Thus :
1 Kz M'(2j 1)nxkc
e cohKz7cos
5 K Ds~ kv k Yk 2a
eKz de (5)
oJ,( = K Jtdg
+ (xk'zk) + cosh KE
k= M '+1
(2J 1)n~ exk'*
M Kz zer [ z~~:
k k
k=1 s
s
k F K
binued on following page)
Z (0) 8 p(T ,, ) e j s z
Z ( O) 2 K Ds 2
erf + / Kk ~ + 2 erf k
S S~
(erf k +~r r K erf z+z
+ os K Z jKz z + z2
+ oh zek efk q
k=1 s
Kz ( [ z z (j1nx
~~~ ~ ~ e r K1: Yk cos (2a (xk~zk
Equation (3.85) with Eq. (3.lc6) used for m (xk~zk) is the complete
solution to Eq. (3.52). 4 g(z), the solution to Eq. (;5.53), is similar to
Eq. (3.85) except that K' replaces K, Bj replaces B j, and sin(jrrxk/a)
(2j 1)nxk
replaces cos wherever they appear. Equation (3.85) and the
2a
corresponding equation for O; (z) are inserted into Eq. (3.4~6), which in
turn is inserted into Eq. (3.17) to give the complete expression for 4(r),
the thermal neutron flux at any point in the medium.
It is instructive to see what happens to Eq. (3.85) when the age
diffusion theory is reduced to onegroup diffusion. If no errors have been
made, Eq. (3.85) should reduce to essentially the same form as Eq. (1.66).
To do this, let Tss O and set p(ts,,) = 1. As Ts Ot the arguments of
k k
all the error functions will approach the form or +  .Ob
S S
viously all the arguments will approach +oo on the real axis, depending in
some cases on whether z < zk or z > z The error function for +oo is
simply +1, and so the last tent in Eq. 0.8.5) becomes
iZk 2
z'/4T
e
2
k= 1
C,( O) 2Ks ~
Z2il
k
s "
s
e
M o
z /4T
s
+ cosh Kz e
k1
l' z > z
> k
13 z < zk
7k co 2aS  "a zk ('6
With a little further manipulation Eq. (g.86) can be further reduced to
the form
M' (2j 1) xx
coshKzkk cos 2a a(xk zk)
k= 1
Zs(0) 9
I (0) D l~
e K
L K
M
cosh KzK Z k c ( lrk
+ e Yk cs 2a a(xk zk)
k= M '+1
(3.87)
Insertion of Eq. (3.87) into Eq. (3.85) reduces the latter equation to
2 3o z > zkk
(2j 1)r xk
2a (xk'zk)
2; z > Zk;z
Kz d@. (g [ Z(0)s) Kz]
e 30s1 e
jiK df, D K
F=o Ds ~
~1~ ~CS 2jlnk cosh Kz
coshKzk cos 2a Q(xk'zk K
k=1
Kzk (2j 1)icx
e k cos 2a s(xk'zkj ]388
k M'+1
The resemblance of Eq. (3.88) to Eq. (1.66) is obvious. Except for
the differences in the shapes of the perturbing foils assumed in the two
derivations, by following a procedure similar to that described at the
end of Chapter ,II the equivalence of the two expressions can easily be
shown.
CHAPTER IV
CALCULATIONAL RESULTS
Two computer codes were written to obtain numerical results from
the analytical results of the previous chapters. The simpler of the two
uses the results of Chapter I with no additional analytical simplifications.
The second code uses the agediffusion model developed in Cha~pter III, bu~t
several restrictions had to be made because of computer storage and numeri
cal problems.
The simpler of the two codes performs the procedures given from
Eq. (1.67) to the end of Chapter I. Computer storage limitations resulted
in several limitations on variable sizes. M, the number of foils, was
limited to 6; using the terminology of Eq. (1.67), L was limited to 10, Q
was limited to 5, and the product LM was limited to 30. The output of the
code consisted of 4 (z) for up to 49 uniformly spaced values of z and the
summation of all the flux modes for the same z values. These results were
given by the code for any number of uniformly spaced frequencies of the
driving source.
Since the onegroup, onedimensional computer code almost completely
filled the computer core, it was obvious that the twodimensional age
diffusion results could not be coded in their entirety. Two major simpli
fications were made. One was to restrict the foil positions and flux
calculational points to the z axis (xk = x = O), thereby transforming the
problem to a onedimensional problem. This results in the x expansion,
Eq. (3.46), reducing to
(O,z = ~ (z .jr (4.1)
(2j 1)nxk
Equation (3.79) simplifies somewhat since the factor cos 2a
becomes unity, and since O; (z) is not used in Eq. (4.1), the equation cor
responding to Eq. (3.79) for 4 Q(z) is not needed. Substitution of Eq.
(4.1), with zk replacing z, into Eq. (3.85) along with replacing the cosine
factors with unity gives the basic equation the code must handle as follows:
Kz de.(g
1 e Ja
QOz = K dg
1
Da
s
M'
CrzK Z k cosh Kzk
j k=1
s k eIZ ~ O,zk )
kM'+1l
j s Kz
Zje K 'kj2
cosh Kz
Q(O,zk) K
Zs(0) 2 DsaW
SKzk C" z z
k=1 ' s
efk + K 2 erf 2 + / K
2 s
Kz [ z z r
+ek efk  K erf z k K]
2 sJ 2 s
2 erf zk + / K + cosh Kz
s k=1
Ez + zk ;K~ Kzk
r~f kz Zk 7K] I 1 4(O,zk) (4.2)
The only other simplification that was made was only the inclusion
of the fundamental spatial mode in the source. Formally, this implies that
Eqs. (1.13), (1.22), and (1.24) apply and that S~m 0 for > 1 and m > 1.
In practice, the higher values (x modes) could have been easily included,
but to have included the higher m values (y modes) would have overloaded
the computer core. This is not an insurmountable obstacle, since magnetic
tape storage could have been used but this is a timeconsuring and there
fore expensive approach. Physically, including only the fundamental spatial
mode seems justified since recent experimental work at the University of
Florida indicates that a highly thermalized fundamental mode is experi
mentally feasible (10).
Using Eq. (4.2) rather than Eq. (1.66), the procedure outlined in
Chapter I from Eq. (1.66) to the end of the chapter may be used to obtain
the complete solution to the problem. This will be done later in this
chapter, but, first, attention will be directed to two factors in Eq. (4.2)
which have physical significance that has not been pointed out previously.
If Eq. (3.8) is inserted into Eq. (3.12), with the age to thermal,
Ts, replacing u, the result is
usza(u')du us
p(.ts,w) =e o e o
du'
(4.3)
Noting that the first integral in Eq. (4.5) is the usual resonance escape
probability to thermal, P(Ts) or Plus), Eq. (4.5) may be written as
iwT
(4.4)
where
u
T =
du'
Szt(u') v(u')
(4.5)
or, in tents of the Fermi age,
Ts SD(T) v()d
(4.6)
Obviously, Ts has the units of time, and represents the time required for
the neutron to slow down to thermal energy. The factor exp [i Ts] then
represents the phase shift introduced by the slowingdown time.
The second factor in Eq. (4.2) which can be simplified is
exp[(K~ B j)Ts]. With the use of Eq. (J.50), this immediately becomes
Z + 
as
(K B j)"
e
( *7)
as s
D
= e s
v D
s s 48
The term in the first exponential is the ratio of the age to thenaal to
the square of the diffusion length, .rs/L Remembering that for a point
source L2 is onesixth the mean square distance that a thermal neutron
travels before being absorbed and that T(u) is onesixth the mean square
distance that a neutron travels in attaining lethargy u, the ratio can be
seen to be a measure of the distance traveled while slowing down relative
to the distance traveled at thermal energy. Typical values for this ratio
are 0.13 for graphite and 0.008 for D20.
Now returning to the actual calculational procedure, as before, the
expansion, Eq. (4.1),must be truncated; i.e.,
i (O1,z = 1 (z) (4.9
j=1
Because of the restriction of the source to the fundamental mode only,
Eq. (4.9) simplifies further, since can only assume the value 1.
Restricting j to 1 in Sj does not, of course, eliminate the higher j com
ponents in the flux since the perturbing foils generate higher x modes but
they cannot generate higher y modes since the foils are assumed to extend
the full length of the medium in the y direction. Equation (3.85) with
1 i .( )Ji, inserted for 4 (x z ) and cos ( )' replaced by unity
 01k kk 2a
may now be written JM times for each of the J modes and for each of the M
foils. This set of equations may then be solved simultaneously for the JM
values of mjl(zk~). Writing these values as a vector,
T 11zj 1(z ""11(M 21zj"" J(zM] ,(4.10o)
63
the set of JM equations may be put into matrix form, as before. That is,
where B is the vector of source tents, [I] is the identity matrix, and [A]
is the matrix containing the Green's functions connecting the foils.B
may be written as
(4.12)
where
S1
11 2DKl
(4.13)
it into J2 submatrices,
given by
The [A] matrix is defined, as before, by
each M by M. The n,p element of the I,K
partitioning
submatrix is
(A = p z ) coh K z KIlzn R)
(I,K n,p = Il p) Ch Il p IlI n
= p z )coshK z eKIlz (R,
BIl p) O' KIl n AlI
,(4.14)
S(4.15)
for n > p
for n < p
where
p ( )  
BIl p) KIlDa '
(4.16)
10
I + 
as Yg
D ~s
'1 P(r ,O) e s
s2 D~a K
11 =(0)
t
[A] 0 + [I] 0 = B ,
(4.11)
T 11 1 11 M
B = 1 e ... a11 e ,O..O ,
(4.17)
64
(RIi n~pe~Ki~~i np 2 (Al) (A2) +2(Ai)j
K1(z +z )
7p(,,
(A )
 (A5)
+ cosh KYlzn Ip !e Ilp [(A2)
Z Z
n p
S
z + z
n p
2, 
S
 K ,
(Al) = erf
(A2) = erf
(A3) = erf
(A4) = erf
(A5) = erf
(A6) = erf
(4.19)
(4.20)
(4.21)
(4.22)
(4.23)
(4.24)
Pz JIK1
2 c  + / K,
2, s I
z z
n p
z + z
n p
2#
 K ,
s 1~n,
,~ K1l
z
P
2Jr
S
The vector 0 is obtained from
L = A] [I] B .
(4.25)
 2(A6))
 11 e I (h ]
,(4.18)
The elements of 0 are inserted into Eq. (3.85) with the previously mentioned
simplifications, which is in turn inserted into Eq. (4.9) to give the com
plete answer to 4 (O,z). The entire thermal neutron flux is obtained by
using Eq. (3.17), which gives
s1 1
j=1
Since the representation of the flux is by expansion in eigenfune
tions, an important characteristic of the numerical calculations is the
rate of convergence with increasing numbers of the eigenfunctions, or
spatial modes. Figures 2 and 3 show flux amplitude versus distance along
the z axis with 1, 2, 4, and 10 spatial modes and with the source modulated
at 100 and 1000 eps, respectively. The model is the onegroup representa
tion of Chapter I with the geometry shown in Fig. 1. One cadmium strip,
0.0508 x 1 x 71.12 em, is positioned 8.89 cm from the end of a graphite
assembly 44.45 x 71.12 em, the dimensions of one of the experimental
facilities available at the University of Florida. In addition to the
curves shown in Figs. 2 and 3, calculations were performed with 6 and 8
spatial modes. At the position of slowest convergence, which is in the
neighborhood of the cadmium foil, the ratios of the 6 and 8mode results
to the 10mode results are approximately 15 and 6/o, respectively, for both
100 and lOOOops source frequencies. Figure 4 shows the phase lag versus
z for 100 ops, corresponding to Fig. 2. In this case the change from 6 to
10 modes and from 8 to 10 modes is approximately 5 and 2/o, respectively.
Figure 5 is the phase lag versus z for 1000 ops, corresponding to Fig. 3.
The convergence is faster for this case, approximately 2.5 and 1o for 6 to
r \ 40 MODES
z 4 MODES
o
2 MODES
J (MODE
U.U1
0 5 r0 15 20 25 30 3
DISTANCE ALONG I AXIS (cm)
Fig. 2. OneGrroup Flux Amplitude vs Distance Along the Z Axis for
One Cadmium Foil in Graphite, at 100 ops, with 1, 2, 4, and 10 Sp~atial
Modes.
rO.O
67
0.01
Modes.
10 MODES
4 MODES
2 MODES
4 MODE
O S l0 13 20 25 30 31
DISTANCE ALONG 7 AXIS (cm)
0.7
0.6
0.1
S0.3
S0.
Fig. 4. OneGroup Flux Phase Lag vs Distance Along the Z Axis for
One Cadmium Foil in Graphite, at 100 eps, with 1, 2, 4, and 10 Spatial
Mode s .
69
I~~ MODESD
4 MMODE
a 5 10 15 20 25 30 35
DISTANCE ALONG I AXIS (cm)
Fig. 5. OneGroup Flux Phase Lag vs Distance Along the X Axis for
One Cadmium Foil in Graphite, at 1000 cps, with 1, 2, 4, and 10 Spatial
Modes.
10 and 8 to 10 modes, respectively. Calculations were limited to 10 spatial
modes because of machine storage limitations. Although the convergence is
somewhat slow, it is probably adequate for most purposes, especially since
the problem used is a rather severe test.
Figures 610 show the results of calculations in the same geometry
as described above except that two perturbing foils are present, at
z = 8.89 and 17.78 cm. These foils are either cadmium, as used previously,
or U2ss of the same total absorption cross section as the cadmium strips.
The calculations with cadmium are onegroup and with U2ss are agediffusion.
Figure 6 shows flux amplitude as a function of z for one spatial mode cal
culations using cadmium foils with source frequencies of 200, 500, and 800
eps and using U~ss foils at 200 eps. The unperturbed flux at 200 eps is
also plotted on Fig. 6. Figure 7 shows the same data using six spatial
modes. The unperturbed flux is identical in Figs. 6 and 7 since the
assembly is driven with fundamental mode only in each case. Figure 8
presents the phase angles as a function of z corresponding to the amplitudes
shown in Figs. 6 and 7. The onegroup, sixspatial mode calculations shown
in Figs. 7 and 8 are replotted as a function of source driving frequency in
Figs. 9 and 10. Figure 9 presents flux amplitude versus frequency for z
values of 2, 8, 8.89, 15, and 3O em. Figure 10 presents phase lag for the
same conditions.
A second facility that is available at the University of Florida for
measurements compatible with the geometry restrictions of these calculations
is a water, or D20, tank which is 120.97 x 130.17 cm in the transverse
directions. Figure 11 is a plot of flux magnitude versus z at 200 eps
for an agediffusion calculation in D2O with two 1in.dina natural uranium
r 200 cps, UNPERTURBED
200 cps, AGE DIFFUSION
\ \ 200 cps
'500 cps >ONE GROUP
'800 cps
) (O 20 30 40 50 60 7
DISTANCE ALONG z AXIS (cm)
0.0(
0.00~
Fig. 6. Flux Amplitude vs Distance Along the Z Axis for Two Cadmium
or Uranium Foils in Graphite, at 200, 500, and 800 ops, with One Spatial
Mode .
0.0r
0.00(
O (0 20 30 40
DISTANCE ALONG I AXIS (cm)
Fig. 7. Flux Amplitude vs Distance Along the
or Uranium Foils in Graphite, at 200, 500, and 800
Modes .
50 60 70
Z Axis for Two Cadmium
cps, with Six Spatial
75
6 MODES~MD 0 p
800 cps
6 MODES
6 MODES GE
3 4 MODE
O~~~ MODE 20304 50 cps
DISTANCE ALONG z ALIS (cm)
Fig. 8. Flux Phase Lag vs Distance Along the Z Axis for Two Cadmium
or Uranium Foils in Graphite, at 200, 500, and 800 cps, with One and Six
Spatial Modes.
z= 2cm
8cm
8.89 c
30 cm
0.(
0.05
0.02
0.0(
400 600
SOURCE FREQUENCY (cps)
4000
Fig. 9. OneGroup Flux
Foils in Graphite, at z = 2,
Modes .
Amplitude vs Source Frequency for Two Cadmium
8, 8. 89, 1~ and 3O cm, with Six Spatial
4.0
3.5
3.0
z = 30 cm
2.5
2.0
cc
S20.
O
c 0 0 0080(O
SOREFEUEC cs
Fi.1.OeGopFu hs a sSuc rqec o w amu
Fol nGahta ,8 .89 5 n Ocwt i pta
Moe
(0.0
5.0
2.0 ~ v
J 6 MODES
OO
4 OD
0.5
0.2
O 40 20 30 40 50
DISTANCE ALONG I AXIS (cm)
Fig. 11. AgeDiffusion Flux Amplitude vs Distance Along the Z Axis
for Two Uranium Foils in 'D0, at 200 eps, with One and Six S~patial MIodes.
77
rods 130.17 cm long located at z = 10.998 ad 21.996 em. The obvious
characteristic of the calculation is that the natural uranium rods do not
offer sufficient perturbation to produce dramatic effects. An alternate
and feasible perturbation could be obtained by using single plates from
fuel elements of the MTR type. One of these plates has approximately the
same total absorption cross section as the natural uranium rods. Therefore,
the perturbation would be approximately the same. In order to produce a
significant perturbation with fuel, fairly large U2ss rods would have to
be used.
Since the calculation shown in Fig. 11 has small perturbations, it
presents a good opportunity to show the relative contributions of the
unperturbed, thermalneu~tron diffusion, and slowingdown contributions to
the final answer. These are shown in Fig. 12, where the three curves cor
respond to the three terms in Eq. (3.85). Since it is a sixspatial mode
calculation, the thermalneutron diffusion and slowingdown curves are sums
of the second and third terms in Eq. (3.85).
UNPERTURBE )
THERMAL DIFFUSION
r I SLOWING DOWN
FOIL POSITIONS
O 10 20 30 40 51
DISTANCE ALONG r AXIS (cm)
Fig. 12. Components of AgeDiffusion Flux Amplitude vs Distance
Along the Z Axis for Two Uranium Foils in D20, at 200 eps, with Six
Spatial Modes.
CHAPTER V
SUMMARY AND CONCLUSIONS
The goal of this dissertation has been to predict, as accurately
as possible, the behavior of neutron waves in heterogeneous multiplying
and nonmultiplying media. The primary approach has been similar to the
heterogeneous reactor theory of Feinberg and Galanin in that the continu
ous slowingdown model has been assumed to describe the behavior of the
neutrons above thermal energy and diffusion theory has been assumed to
describe the behavior of the thermal energy neutrons. Rather than assume
that the infinitemedium slowingdown and diffusion kernels are applicable,
a more basic approach was chosen which resulted in slowingdown and diffu
sion kernels for the particular finite geometry that was assumed.
The preceding chapter showed some results of calculations using
this theory. In this chapter, it will be shown how this theory may be
applied and extended.
Application to Critical Reactors
As stated previously, the theoretical developments of this work are,
in several ways, an extension of the FeinbergGalanin heterogeneous reactor
theory. Aside from the obvious feature of predicting the behavior of neu
tron waves in subcritical assemblies, this work takes into account the
finiteness of the assembly. The effect that this has is best shown by
applying the theory to a critical reactor, which can be done by dropping
80
the source terms and changing the boundary conditions slightly. The ab
sorptive Green's function defined by Eq. (3.60) will have the same boundary
conditions as before except that now G.(z,S) itself, rather than the deriv
ative of G.(z,g), will be required to go to zero at z = 0. This change
in boundary condition simply results in G.(z,S) changing to
e sinh KS 51
G (z,S) = K z > E,(5
e'K sinh Kz
= K z < 5 (5.2)
With the new Green's function, Eq. (3.65) applies to the critical assembly
with the added simplification that all the boundary conditions given by
the last term in Eq. (3.65) vanish. Equation (3.65) now becomes
oo
O (z) = G (zk'z) Hjk(z) + Aid di G (i,z) Fjk(II.Izk3
k k o
where G.(Stz) is now given by Eqs. (5.1) and (5.2). Remembering that the
factrs H and F.k both contain
Jkak
it can be seen that, for k = 1,2,...,M, and j = 1,2,...,J, JM homogeneous
equations are available for each value of a. For the system to be critical
the determinant of the coefficients of the 4 Q(zk) must be zero. This of
course gives an equation which can be solved f~or rl (contained in A j) or
whatever other eigenvalue is chosen to determine criticality.
The first difference noticed in comparison with the FeinbergGalanin
theory seems to be simply additional complexity since the FeinbergGalanin
theory only requires the solution of an M by M determinant. The improve
ments are the improved kernels contained in Eq. (5.5). As pointed out in
Chapter III, Aie and Fjk contain the finitemedium slowingdown kernel
rather than the infinitemedium kernel used by FeinbergGalanin. The
Green's function G. also is an improvement, being the finitemedium diffu
sion kernel for the geometry used in this work, rather than the infinite
medium diffusion kernel used by FeinbergGalanin.
Conclusions
The theoretical developments of this dissertation probably represent
a model which is as complex and physically realistic as is practical at this
time. Further improvements, in order to be practical for machine calcula
tions, would certainly be desirable if they simplified the mathematics or
if they eliminated remaining physically undesirable assumptions, or both.
If such improvements are not easily attainable, the primary question remain
ing is how to use the theory best to direct future experimental work in this
area.
Several approaches are possible. One is simply to continue full
calculations of the thermal flux in various configurations, such as those
presented in the preceding chapter. This is a rather brute force technique
but at least has the advantage that calculations are cheaper and much more
easily accomplished than experiments since the machine codes are available.
Some possibly more satisfying approaches can be suggested by reexamination
of some of the results of Chapter III. To this end, write Eq. (3.64) as
1 L k (2 1)x
me (z)= H G (z,zk) cos 2a (xk'zk)
Ee Z 7 (2j 1)TIxk oc4r iq()0G(k~
k o
where
Kz do. (5)
H e (5.6)
je K dS
5=0
and
B .1s
I (0) a p(T ,W) e a
j I (0) 2D D57
t s
atnd Gj and Mk are defined by Eqs. (j.61) and (3.68), respectively. A cor
responding equation exists for Qj (z); thus
1~~ YkJxk
mj (z) = H; zDk sin ak (xk~zk)
+ sin ~ka (xkzk) di Mk(i) G (zk'n ,(58
k o
where G!, H; and E'. are obtained by replacing B with B'~ and K~ with
K'2 where they appear in Eqs. (3.61), (5.6), and (5.7), respectively.
Insertin Eqs. (5.5) and (5.8) into Eq. (3.46) and in turn into Eq. (3.17)
results in
1 (0 1)n 1 (2j 1)nx
+ 1 sin janx H', .zz)C
+ G! (z'zk) Sjk] (xk~izk) + Ea Cjk e(xk'zk)
diM~)G. (zk,8)~ +~ E Sjk e(x'zk S at Mk(i) Gr (zk^>
o k o
(5.9)
whe re
7k c (2j 1)nx (2j 1)nxki
C.Ecscos (5.10)
jk a 2a 2a
and
S. Y sin sin rx (5.11)
jk a a a
Experimental measurements of thermal neutron flux as a function of position
and frequency represent the quantity calculated in Eq. (5.9). These data
may be operated on with
a b
S~ (2j 1)rtx (2 1)n
2a 2t
a o
or with
a b
in sing jx dy,~ cos (2R 2bl)r[y (515)
a b
to yield experimental measurements of the flux components calculated,
respectively, in Eqs. (55) and (5.8). This entire procedure may be per
formed experimentally with the moderator alone, with fuel plates or rods
and with purely absorbing plates or rods of the same macroscopic absorbing
cross section as the fuel assemblies. The technique of replacing fuel as
semblies with comparable pure absorbers has been utilized previously to mea
sure effective delayed neutron fractions (11). The six terms in Eqs. (5.5)
and (5.8) may now be separated. The homogeneous experiment yields Hj and
H' Subtracting the homogeneous term from the pureabsorber data yields
the absorption terms, the second terms in Eqs. (5.5) and (5.8). Similarly,
the difference between the pure absorber and the fuel assembly data yields
the last terms which contain the slowingdown kernels.
Since the above experimental data would presumably be obtained by
cadmiumdifference measurements with a neutron detector which has approxi
mately a 1/v response, it is also of interest to calculate the epicadmium
flux or cadmiumcovered detector response. Equation (g.44) may be rewritten
as
Z (0) Z k" t
e Cjk + e Sjk (xk'zk) (5.14)
The argument r has been added to Mk as a reminder that Mk is defined as in
Eq. (3.68) except that I replaces Ts. Insertion of Eqs. (5.14) into Eq.
(g.18) yields
Z (0)
s(o
if
q(r,, ) = p(T~W)
c (2 1)ny
cos 2b
Sjk
Mk~( z,T)
j k
a(xk)zk)
(5.15)
Cjk
The response of the 1/v detector covered
with cadmium can be taken to be
proportional to
ucd(r, ) 1 du
S ed
= 2E
m
(r, )
epi
u
= edq(rT) eu/2 du ,
(5.16)
(5.15)
where ucd is the lethargy of the
into Eq. (5.16) gives
cadmium cutoff.
Insertion of Eq.
z (0)
s ~0
t
S1 c (20 1)Iy
=' cos 2b
epi(r )
(5.17)
(5.18)
and
e ~j
Z Iejk Cjki jk Sjk ke(i~zk)
j k
whe re
ud B~j .u/~
Iejk ~Cp(u,w) N ~(z,u) e e/ u
ucd B' u
I~k pu~)Mkz~) u./2 du (5.19)
As before, experimental measurements of the quantity represented by Eq.
(5.17) taken as a function of position and frequency may be operated on
with the operators given in Eqs. (5.12) and (5.13). Let the result of
operating with Eq. (5.12) on Eq. (5.17) be called F'Ajk, which is
2 (0) 7 (2j 1) xk
Fejk Zs0) jk Icos 2a (xkr~zk). (5.20)
The nonsyrmmetric part may be called F~ijk and is obtained by operating with
Eq. (5.13) on Eq. (5.17), giving
F'Ajkt Z jk sin k +(xk~zk) (5.21)
Several things should be noted at this point. First, all the inform
ation obtained thus far could have been obtained from experimental data
taken at only one z value. Secondly, each term obtained experimentally
except for the homogeneous term has contained a weighted summation over the
thermal neutron flux at the surface of the foils. Alternatively, by com
bining the thermal constant, yk, with the flux, these summations may be
considered to include the net current into each of the foils.
Assume now that the experimental data and following analysis have
been performed for M values of z, different from the locations of the M
foils, of course. These M values of z should be chosen so as to be fairly
close to the M foil locations in order to maximize the information about
87
the foils in the M sets of measurements. Now assume that this set of
information is available for the synmmetric absorption term, the second
term in Eq. (5.5). A set of M inhomogeneous equations results, with the M
unknowns being the net current into each of the M foils or, if 7k is known,
being the surface flux at each of the M foils.
Now that the net current into each foil is known, this information
may be inserted into some of the other terms to obtain additional informa
tion. For instance, having the net currents and experimental values of
F~jk at M points, the integrals Iljk defined in Eq. (5.18) may be obtained
by solving the set of M equations given by Eq. (5.20). Alternatively, if
the E .j's are assumed to be known, the integrals over the slowingdown
kernel in Eq. (5.5) may be obtained for each foil.
APPENDIX A
THE COMPLEX ARITHMETIC VERSION OF GINV
GINV is a F RTRAN subroutine developed by Burrus et al. (12) for
the purpose of obtaining working inverses to singular matrices or to matrices
which appear singular to a digital computer. It was modified to handle
complex arithmetic and was used in the programs employed for calculations
in this dissertation because, at one point, it appeared that the matrices
in these problems were very ill conditioned. Although this turned out not
to be the case, it was left in because it takes up little more storage than
a standard matrix inverter and gives slightly better answers in practically
all situations.
Theoretically, the inverse of a matrix A, A l, has the property that
A A = I, where I is the identity matrix. This property never is true
precisely, when an inverse is obtained on a digital computer with finite
accuracy. GINV operates so as to find the working inverse A defined as
that matrix for which the root mean square value of (I A A) is a minimum.
The theory of operation will not be given here since it is not
pertinent to this dissertation and should be available soon in the
literature (13).
Pages 9092 contain the F RTRAN listing of the program as modified
for complex arithmetic. Sufficient comment cards are included to enable
a competent programmer to modify the program for ordinary F RTRAN arithmetic.
A few comments may be appropriate on the use of the program. First,
when TAU is set to zero, the routine performs virtually the same operations
as other matrix inverters and should give identical results on well
conditioned matrices. It may also be noted that if the problem to be solved
is to obtain x from Ax = b, where A is an m x n matrix (with n < m) and b is
a given m component vector, GINV automatically obtains the leastsquares
solution to the set of n equations. It is obviously much more versatile
than most matrix routines, although its full versatility was not needed
in this work.
SUBROUTINE GINV(A,MR,NRNGIfAU,UATEMP)
C UPON ENTRY A = A MATRIX WITH NR RUWS ANO NC COLUMNS
L AFTER THE RETURN THE ORIGINAL A IS DESTROYED ANU THE
C ARRAY A CONIAINS THE TRANSPOSE OF T, THE TRANSFURMATlON
C MATRIX (X = T+B) WHICH SELVES THE PAIR OF EQUATIONS
C A+X = 8 ANO TAU leX = 0
C IN THE LEAST SQUARE SENSE.
C MR = IST DIM. NO. OF ARRAY A IN THE CALLING PROGRAMS
C FAU IS A NONNEGATIVE CONSTANT WHICH CONTROLS THE ERROR
C PROPAGAllONN OF THE TRANSFORMATION MATRIX. IF FAU = 0.
C THEN THE RESULTING TRANSFORMATION MATRIX IS THE
C ORDINARY LEAST SQUARES TRANSFORMATION MATRIX.
C (THE INVERSE IF NR = NC)
C ASTEMP AIJo U ARE USED FOR WORKING SPACE 8Y THE ALGORITHM
C AND 00 NOf NECESSARILY CONTAIN ANY RELEVANT NUMBERS AT
i IHE CONCLUSION. ATEMP MUST BE DIMENSIONED AT LEAST NC
C t)Y THE CALLING PROGRAMS. U MUST BE DIMENSIONED AT
C LEASE NC*NC UY THE CALLING PROGRAMS.
C NO. OF MULTIPLICATIONS = NC+*2 (5/2 NR + 2/3 NC)
I DIMENSION A(30,30), U(30,30), ATEMP(30), 0001l),
L DUT2(1), TAUSQ(1), DUM(1)
C THE DIM. ST. AS SENT HY RO\SS 8URRUS FOR FO~RRAN IV WAS
DIMENSION A(MRNC),U(NCINC),ATEMP(NL)
C NOTE THAT IN COMPLEX ARITHMETIC VERSION, IT HAS 8EEN
C ASSUMED THAT MR IS ALSO fHE 2NO DIM. NO. OF A AS WELL
C AS THE 15T.
TAUSQ(1) = CAU**2
TAUSQ(2) = 0.0
C PLACE UNIT MATRIX'IN U
00 5 I = LNC
00) 4 J = 1,N~C
1 4 U(IJ) = (0.0,0.0)
I U( )= .0 0 0
C URTHOGONALILL COMBINED MATRIX (A ABOVE U) BY GRAMM
C SCHMIDTHIL~ftRT METHOD WITH FIRST NR ROWS WEIGHED
C wifH L ANO THE OTHER NC ROWS WEIGHTED WITH 1/TAU. THEN
SREORTHOGONALLZ TO LESSEN ROUNOUFF ERROR.
00 20 1 = 1,NC
IMR = I + MR
II = I 1
IF (II) 2,11,2
2 011 10 LL = 1,2
00 10 J = ,II
JMR = J+MR
I DOT =(0.0,0.0)
1 00T2 =(0.0,0.0)
UO 3 K = 1,NR
DUM(1) = A(KJ)
DUM(2) = A(KJMR)
1 3 DOT = A(KI)*DUM + UOT
00 6 K = 1,J
OUM(1) = U(K,J)
DUM (2) = U(KJMR)
I 6 OOT2 = U(K,1)*0UM + OUT2
I UOT = 00T + 00T2*+TAUSQ
00 8 K = 1,J
1 8 U(KI) = U(K,1) 00T*U(KIJI
00 10 K = IrNR
I 10 A(KI) = A(KI) 00T*A(KJ)
C NORMALI:L THEf COLUMN I OF THE COM81NED MATRIX
1 11 007 =(0.0,0.0)
I DOT2 =(0.0,0.0)
DO 12 K = 1,NR
12 UUf = 00T + A(KI)**2 + A(KIMR)**2
00 14 K = 1,I
14 0012 = 00T2 + U(KII)**2 + U(KIMR)**2
00f  OUT + 00T2*TAUSQ
UUT = SURTF[00T)
00 17 K = 1Ir
1 17 U(K,1) = U(K,I)/UOT
00 19 K = 1,NR
I 19 A(KI) = A(KI)/DOT
20 CONTINUE