A THEORETICAL AND EXPERIMENTAL STUDY
OF NEUTRON WAVE PROPAGATION IN
MODERATING MEDIA
By
RAY STURGIS BOOTH
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMlENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
April, 1965
Copyrinht by
Ray Sturgis Booth
and
Department of :Iuclear Engineering!
University ~of Floridla
19 65
ACKiNOWJLEDGMENlTS
The author wishes to express his sincere appreciation
to the members of his supervisory committee for their
assistance in the preparation of this dissertation. Special
recognition is due Dr. Rafael B. Perez for his continual
guidance and encouragement. His advice and profound in-
sight have been invaluable throughout all phases of this
work.
The author is grateful to Dr. M. J. Ohanian for his
major contribution to the theoretical aspects of the
dissertation through stimulating discussions and the contri-
bution of useful computer codes. Important in the experi-
mental work has been the generous aid of Robert H. Hartley
and his careful and exacting experimental technique during
the acquisition of mutually valuable data.
Other persons whose contributions to this work are
gratefully recognized include Robert B. Razminas for the
sharing of his DIAGONAL-CORRECTOR code, the staff of the
University of Florida Computing Center for valuable
assistance in connection with many codes, George Fogle for
the contribution of his technical know-how to the experi-
ment, and Barbara Gyles, who took more than a typist's
interest in turning out the final manuscript.
iii
The author is deeply indebted to his wife, Judy,
for her devoted encouragement, understanding, and patience.
Her valuable help in the preparation of the manuscript is
also appreciated.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS ....
LIST OF TABLES ....
LIST OF FIGURES ....
LIST OF SYMBOLS ....
. . . ii
. . . vi
vii
. . . . .
.
ABSTRACT
. . . . xiii
* *
Chapter
I.
II.
INTRODUCTION ..
. . .1
REVIEW OF THE LITERATURE
III.
IV.
THE BOLTZMANN EQUATION IN THE DIFFUSION
APPROXIMATION WITHIN THE FRAMEWORK OF
NEUTRON THIERMALIZATION THEORY .....
THE EXPERIMENT IN GRAPHIITE ......
. 11
. 58
. 99
V. COMPARISON OF THE THEORY WITH THE
EXPERIMENT ..............
VI.
VII.
RELATIONSHIPS BETWEEN THE NEUTRON WAVE
EXPERIMENT AND DIE-AWAY EXPERIMENT IN
THE DETERMINATION OF DIFFUSION AND
THERMALIZATION PARAMETERS .......
CONCLUSIONS AND RECOMMENDATIONS ....
124
136
194
199
178
. .
. .
. .
LIST OF REFERENCES ....... ........
APPENDIX . . . . .
BIOGRAPHICAL SKETCH ...............
. .
LIST OF TABLES
Table
1.
2.
3.
4.
5.
6.
Page
The Energy Points Used in the Discrete
Energy Representation and the Calculated
Values of the Diffusion Coefficient at
These Points . . . .
Parameters Used in the Numerical
Calculations . . . .
Ligenvalues of the Space and Energy
Dependent Boltzmann Operator Using the
SUMMIT Kernel and the MASSMK Kernel ....
The Experimental and Theoretical Values for
the p(n) Coefficients ...........
The Experimental and Theoretical Values
for the Q(n) Coefficients .........
The Static and Dynamic Parameters .....
39
41
45
122
132
133
Parameters as Determined from Various
Die-Away Experiments Compared with the
Results of This Work .........
. .134
LIST OF FIGURES
Figure
1. Flow Diagram of Several of the Computer
Codes Used in the Theory ...... ..
2. Several Eigenfunctions of the Space and
Energy Dependent Boltzmann Operator
Using the SUMMIT Kernel. .........
3. Several Eigenfunctions of the Space and
Energy Dependent Boltzmann Operator
Using the MASSMK Kernel. .........
4. Real Component of the Complex Inverse
Relaxation Length for the Fundamental and
First Four Energy Harmonics as Calculated
Using the MASSMK Kernel. .........
5. Imaginary Component of the Complex Inverse
Relaxation Length for the Fundamental and
First Four Energy Harmonics as Calculated
Using the MASSMK Kernel..........
6. Real Component of the Complex Inverse
Relaxation Length for the Fundamental and
First Three Energy Harmonics as Calculated
Using the SUMMIT Kernel..... ....
7. Imaginary Component of the Complex Inverse
Relaxation Length for the Fundamental and
First Three Energy Hlarmonics as Calculated
Using the SUMMIT Kernel. .........
8. Real Component of the Complex Inverse
Relaxation Length for the Fundamental and
First Two Space Hlarmonics as Calculated
Using the MASSMK Kernel.
9. Imaginary Component of the Complex Inverse
Relaxation Length for the Fundamental and
First Two Space Harmonics as Calculated
Using the MASSMK Kernel. .........
Page
.43
. Q
.50
.54
.56
vii
LIST OF FIGURES (Continued)
Figure Page
10. Arrangement of the Experiment and Details
of the Thermalizing Tank (from Refer-
ence 14). . * * 60
11. A Block Diagram of the Electronic Equip-
ment Used in the Wave Experiment. ...... 66
12. Time Sequence of Electrical Signals Associ-
ated with the Electronic Equipment Used in
the Wave Experiment . 68
13. A Block Diagram of Electronic Equipment
Used in the Pulse Propagation Experiment
(from Reference 14) .. . 71
14. Thermal Wave of the 200 ops Wave Experiment
for Several z Positions .. .... . 72
15. Thermal Pulse for the Pulse Propagation
Experiment for Several z Positions (from
Reference 14) . . * 74
16. Amplitude of the 400 cps Wave as Obtained
from the First Harmonic of the 200 ops Wave
Data and the Fundamental of the 400 ops
Wave Data .. ... .. .. * 82
17. Phase of the 400 ops Wave as Obtained from the
First Harmonic of the 200 ops Wave Data and
Fundamental of the 400 ops Wave Data. .. 83
18. Combined Amplitude Curves for Several
Frequencies . . . 93
19. Combined Phase Curves for Several of the
Frequencies .. .. .. * * 94
20. Comparison of the Theoretically Derived a
to the Experimental Results ......... 100
21. Comparison of the Theoretically Derived 5
to the Experimental Results ......... 101
22. Comparison of the Theoretically Derived 2aS
to the Experimental Results ........ 106
viii
LIST OF FIGURES (Continued)
Figure Page
23. Comparison of the Theoretically Derived
22 52 to the Experimental Results for
Frequencies up to 1,400 cps. .. ... 107
24. Comparison of the Theoretically Derived
a2 52 to the Experimental Results for
Frequencies up to 400 cps. ... .. .. 108
25. The Ptn) as a Function of the Maximum Fre-
quency Included in the Least-Squares Fit 112
26. The Experimentally Determined Real Part of
1 versus Distance from the Source. .. .. 115
27. The Experimentally Determined Imaginary
Part of t versus Distance from the Source. 116
28. The Experimentally Determined a' and 5' as
a Function of Frequency. .. .. .. .. 117
29. The First Frequency Derivatives of 2aS and
a2 2 VerSUS . . . 118
30. A Pictorial Representation of the Die-Away
and Neutron Wave Experiments .. .. .. 126
LIST OF SYMBIOLS
8: Transverse buckling of the fundamental
space mode (cm-2)
B2 Square of the fundamental complex inverse
relaxation length for BJ=0. Also used as
a general eigenvalue in one part of the
theory (cm-2)
Bn Square of the complex inverse relaxation
length for B:=0 associated with the higher
energy modes (cm-2.)
B K) The Kh coefficient of the Maclaurin's
series expansion of 82 (Cm-2sec )
C, Diffusion cooling coefficient (cm sec'l)
D(E) Energy dependent diffusion coefficient
Do Energy a ve ra ged diffusion coefficient
(cm2sec-1)
E rleutron energy (ev)
7, Coefficient of the B* term in the power
series expansion of X in power of B2
(cm sec-1)
Lmn Matrix elements of the energy expansion
in eigenfunctions of the P operator
(cm-2 Sec)
M(E) Maxwellian distribution (number density)
P Energy and space dependent Boltzmann
operator
(n)
P The coefficients of the power series fit
of the experimental data (cm' seen)
0(n) Coefficients of the power series expansion
of p2 (Cm2)n-1
LIST OF SYMOBLS (continued)
S, The nth ratio associated with the
fundamental eigenvalue
S(K) The Kt coefficient in the Maclaurin's
n series of Sn (secK)
t Time (sec)
u(I-J) Heaviside unit step function
vo Reference neutron speed for thermal
neutrons (2.198xl05 cm sec'l)
Wn(q,v) The nth Fourier width of the measured
pulse in the pulse propagation experiment
x,y,z Space variables of the rectangular
coordinate system (cm)
ao v Eac (sec-1)
6mn Kronecker delta
n Variable used in the discrete energy
representation (n = (E/KT)1/2/(1+(E/KT)1/2)
X Decay constant measured in the die-
away experiment or calculated from the
Q(nI with ao0 (sec-l)
1W Decay constant calculated from the Q(n)
with so0 (sec'l)
v,f Frequency (cps)
02 Squared fundamental complex inverse
relaxation length. Also used in one
section of the theory as the general
eigenvalue (cm'3)
noThe nth eigenvalue of the P operator
(cm-2>
Zao Energy averaged macroscopic absorption
cross section (cm'l)
LIST OF SYMBOLS (continued)
Energy dependent macroscopic scattering
cross section (cm'l)
Energy dependent macroscopic absorption
cross section (cm-l)
Energy dependent second Legendre moment
of the angular dependent scattering cross
section
Energy dependent differential scattering
cross section (cm-l/ev)
Energy dependent macroscopic total
collision cross section (cm-1)
Neutron flux (n cm-2sec' eV-1)
Spatial eigenfunction used in the
space expansion
Summetrized neutron flux
The nth eigenfunction of the P Operator
The measured pulse of the pulse propa-
gation experiment
The nth Fourier moment of ~(q,t)
Frequency (rad/sec)
Laplacian operator
Laplacian operator for the x and y
variables
C,(E)
Cs(E',E)
ct(E)
e, ,(x,y)
X(E)
Qn(E)
d(q,t)
lt (q,v)
V2
X11
Abstract of Dissertation Presented to the Graduate Council
in Partial Fulfillment of the Requirements for
the Degree of Doctor of Philosophy
A THEORETICAL AND EXPERIMENTAL STUDY
OF NEUTRON WAVE PROPAGATION IN
MODERATING MEDIA
By
Ray Sturg~is Booth
April, 1965
Chairman: Dr. Rafael B. Perez
Major Department: Nuclear Engineering
The space, energy, and time or frequency dependent
Boltzmann equation in the diffusion approximation has been
used within the framework of neutron thermalization theory
to predict thermal neutron wave propagation in a moderat-
ing medium. This theory was then compared with the
experiment conducted on a graphite pile.
In the theoretical section the energy dependence
was expanded in the space-energy eigenfunctions of the
Boltzmann operator, and the transverse-space dependence
was expanded in the orthogonal solutions to the Helmholtz
equation with homogeneous boundary conditions. It was
then shown that each of these space and energy modes thus
formed propagate with characteristic complex inverse
xiii
relaxation lengths. The perturbation method developed by
Perez was used to expand these complex eigenvalues in a
Maclaurin's series of the frequency variable. In cases
where this series did not rapidly converge a direct
solution of the complex secular determinant was obtained.
Two scattering kernels were selected for implementation
into the theory, the MASSMK kernel and the SUMMIT kernel
based primarily upon derivations by Parks. The numeri-
cal calculations of the theory were performed on an
IBM-709 digital computer and founded upon the previous
work of Ohanian.
In the experimental section an accurate technique
for conducting neutron wave experiments in multiplying or
diffusing media was presented along with the results of
the graphite experiment. The excitation of neutron waves,
both by modulated and pulsed sources, and the electronic
systems associated with their detection were discussed.
Special attention was given to the analysis of data for
which several computer codes were developed.
In the comparison of theory and experiment it was
shown that the neutron wave technique allows for the
separation of diffusion and thermalization effects, thereby
yielding a more exacting comparison of theoretical and
experimental results. Polynomial fits of functions formed
X1V
from the experimental results yielded coefficients to be
compared with the theoretically derived terms of the
perturbation expansion. The results of the thermal pulse
propagation experiment were also analyzed by a method
suggested by Moore, and confirmation of his predictions
was established.
The diffusion and thermalization parameters
determined from the data analysis were then compared with
results of several die-away experiments and with the
theory. The complementary relationships between the
neutron wave and die-away experiments were also presented.
The diffusion and thermalization parameters
determined in this work for graphite of density 1.60 gm/cm3
were vCa = (91.'1.)sec-1, vD = (2.16'.01)xl05 cm2sec-1,
C = (39._'2.)x105 cmksec-1, F = (12 .'2.)x107 cm6sec-1,
CHAPTER I
INTRODUCTION
Placement of a modulated or pulsed source of thermal
neutrons at one boundary of a nuclear assembly creates a
disturbance which can be analyzed in terms of wave components
propagating from this source with frequency dependent phase
shift and attenuation factors. The importance of the
neutron wave experiment is based upon the fact that these
two independently measurable quantities are strongly dependent
upon the diffusion and thermalization properties of the
assembly. Thus accurate determinations of physical con-
stants and severe tests of theoretical models are possible
through an experimental determination of the frequency
dependent dispersive properties (dispersion law) of the
system.
Several experimental methods for investigation of
moderating and multiplicative assemblies are now in exten-
sive use. They include, among others, exponential experiments,
"poisoning" experiments, and the die-away technique. It is
interesting to note that the neutron wave technique, which
is the subject of this dissertation, is related to each of
the three methods mentioned above. Basically, the neutron
wave experiment is an extension of the exponential experiment
in that the constant source of the latter is replaced by a
modulated or pulsed source. One then measures the complex
inverse relaxation length, p, of the neutron waves as they
travel through the assembly. An extrapolation to zero fre-
quency of the real part of p yields the same value, Bm a
obtained by the exponential experiment. Using a modulated
source, the material buckling is obtained through an extrap-
olation of a curve, which is nearly constant near the zero
frequency intercept, rather than by measuring a single
point on the curve. Also the imaginary part of p is still
available for determining another parameter or comparing
with the theory. Still another possibility with the neutron
wave experiment is the appearance of resonance phenomena
arising from the splitting of each spatial mode into various
energy modes which propagate with different complex inverse
relaxation lengths.
Increasing the frequency of modulation in a neutron
wave experiment can be considered equivalent to "poisoning"
the medium with a 1/v poison. Therefore, determining the
bucklings and relaxation lengths at several poison concen-
trations in the "poisoning" experiment and finding the complex
relaxation length as a function of frequency in the wave
experiment are complementary measurements. Moreover, the
frequency-dependent, energy spectrum of the neutron1 wave
experiment bears a close resemblance to the energy spectrum
associated with the "poisoning" experiment as a function of
the poison concentration. In the wave experiment, however,
solids as well as liquids can be homogeneously "poisoned."
The neutron wave technique and the die-away experi-
ment also bear a close resemblance. The modulus of the
complex inverse relaxation length can be thought of as a
negative buckling. Its relationship to w/v, the "poisoning"
term, is quite similar to the relationship between the geo-
metric buckling and the decay constant as measured by the
die-away technique. In the die-away experiment one determines
X, the lowest time-energy eigenvalue of the Boltzmann operator,
as a function of the geometric buckling. This method, though
quite effective for the determination of diffusion parameters,
has severe limitations when applied to the determination of
thermalization parameters unless the time-dependent, energy
spectrum is measured. At high bucklings, thermalization
effects cause the curve of X versus buckling to deviate
slightly from the linear relationship predicted by diffusion
theory. This second-order effect occurs precisely at the
range of buckling where the least confidence can be given
to the experimental results; for it is in this range of
buckling that the extrapolation distance is most important,
transport corrections become significant, and the energy
dependence takes the longest time to stabilize. Also an
asymptotic decay is never attained in some crystalline modera-
tors because of the so-called "trapped neutrons" which decay
in time much slower than those associated with the
fundamental eigenvalue. Thermalization effects are of
first-order significance in the neutron wave experiment.
Diffusion theory predicts that the real part of p2 should
be independent of frequency, whereas it is in fact nearly
linearly related to the square of the frequency. The
measurement of thermalization parameters via the neutron
wave technique can be accomplished in large systems where
angular dependence and extrapolation distances are of
comparatively minor importance.
The purpose of this dissertation is to present both
theoretical and experimental investigations of the propa-
gation of neutron waves in a moderating medium. The
Boltzmann equation in the diffusion approximation is used
as the theoretical model, and its energy dependence is
expanded in the space-energy eigenfunctions of the Boltzmann
operator. A perturbation method is used to expand the
square of the complex inverse relaxation length in a
Maclaurin's series of the jw variable and in the x variable
which is related to the frequency of excitation by the
expression w = xex. In the experimental section an accurate
technique for conducting neutron wave experiments in multi-
plying or diffusing media is presented along with the
experiment in graphite. Special attention is then given
to the analysis of the data thereby obtained and to its
comparison with the theory.
* Underlined numbers in parentheses refer to the List
of References.
CHAPTER II
REVIEW OF THE LITERATURE
The propagation of neutron waves caused by
oscillating a thermal-neutron absorber in a chain-reacting
pile was investigated by Weinberg and Schweinler (1)."
Weinberg and Wigner (2) derived equations for the atten-
uation distance, wave length, and propagation velocity as
a function of frequency for waves excited by a localized
source in an infinite diffusing or multiplying medium.
It was shown that the diffusion coefficient and the prompt
material buckling could be determined by measuring the
propagation characteristics of these neutron waves.
Experimental measurements of diffusion parameters in
graphite and heavy water were performed by Raievski and
Horowitz (3) by using a mechanical source and a frequency
range which varied from 50 ops to 500 ops. Their results
were in good agreement with theoretical predictions, but
no attempt was made to extend the theory to include therma-
lization effects.
Preliminary experiments in graphite and in water-
moderated subcritical assemblies were done by Uhrig (4,5).
He used one-group diffusion theory to derive equations
relating the prompt material buckling and diffusion co-
efficient to the neutron wave's phase shift and attenuation
and compared his results to experiments conducted on the
graphite, subcritical reactor at lowa State University (6).
A mechanical source was used in these experiments; and the
basic unit in the detection system was a dual, coincidence-
anticoincidence switching circuit which directed the detector
pulses to four slave scalers so that each scaler counted
during a different phase interval of the modulated source
(7). It was determined that the waves do follow the ten-
dency predicted by theory; however, the wave velocity was
in all cases higher than expected, with the deviation becoming
greater as the oscillation frequency increased.
The theoretical work of Uhrig was extended by Booth
(8) to a two-group diffusion theory model. It was shown
that when two energy groups are considered two thermal waves
propagate in the fundamental spatial mode, that the relative
magnitude of the Fermi age compared with the square of the
diffusion length is of predominant importance in determining
the propagating characteristics of the two waves, that K,
is a second-order effect with respect to these same charac-
teristics, and that the material buckling can be measured
from experiments conducted in the central region of the
assembly. Several methods for the detection of neutron
waves were also presented by Booth.
Perez and Ulhrig (9) investigated the propagation
of neutron waves in a moderating media with a one-group
diffusion model and a thermalization model. In the
thermalization model the energy dependence was expanded
in the Associated Laguerre Polynomials of order one. The
complex inverse relaxation length was expressed in a power
series of the ratio w/vD, and it was shown that the co-
efficients of this expansion are functions of the thermal-
ization and diffusion properties of the moderator. The
authors demonstrated that the neutron wave technique affords
a complementary method to the "poisoning" technique and that
the difficulties encountered in the die-away experiment
because of the trap effect should not be present in a
neutron wave experiment.
The diffusion and cooling coefficients of graphite
were determined experimentally by Perez, Booth, and Hartley
(10) by the neutron wave technique. The value of DO agreed
with previous measurements, but the value of CO was a factor
of two larger than that obtained by Starr and Price (11)
using the die-away technique. The propagation of neutron
waves in a subcritical assembly using multigroup diffusion
theory was investigated by Perez, Booth, Denning, and
Hartley (12). It was shown that for "n" groups each spatial
mode splits into "n"' damped waves of the same frequency but
with different phase shift and attenuation factors, and that
the interaction of these waves presents the possibility for
measurable resonance effects which are quite sensitive to
the lattice parameters of the system. An attempt was made
to demonstrate thermalization effects in moderators by the
reflection of neutron waves by Denning, Booth, and Perez
(13). In the experiment a sinusoidally modulated source
of neutrons was placed at one face of a 26.7 cm. long
graphite block with transverse buckling of 4.0 x 10-3 cm-2
which preceded a 62 cm. long heavy-water region of the same
transverse area. Theoretical investigations indicated that
if any reflection occurred at the graphite, heavy-water inter-
face it would have to be caused by energy modes which pro-
pagate only when the energy dependence of the neutrons differs
from a Maxwellian distribution. This was qualitatively con-
firmed by the experiment.
The propagation of thermal pulses was demonstrated by
H-artley (14) and Miley and Doshi (15). Hartley produced
thermal pulses of about one millisecond width by use of a
neutron generator and thermalizing tank, recorded the pulses
at various distances from the source, and compared the
results in the time domain to a one-group diffusion theory
model. An extension of the theory to that of Perez and
Uhrig met with severe difficulties when the inversion to
the time domain of the complicated Laplace-transformed flux
was encountered. The disagreement between the experiment
and the one-group diffusion theory showed the sensitivity
of the experiment to thermalization effects. Miley and
Doshi used the TRIGA Reactor as a source and thermal pulses
were produced of width ranging between 32 and 66 milli-
seconds. The measured asymptotic velocity of the pulse
yielded a value of Do; but because of the large pulse width
(i.e., lack of high frequency content), the more interesting
thermalization effects observed by Hartley were not detected.
Diffusion theory was compared with the experiment, and
corrections were made for the source distribution.
A method of comparing theory and experiment in the
neutron wave experiment with pulsed, thermal sources was
suggested by Moore (16) so that the inversion of the theory
to the time domain (which proved so difficult for Hartley)
would not be necessary. In essence, Moore suggested that
comparisons be made in the frequency domain by numerically
determining Fourier moments and widths of the data measured
in the time domain. Courier widths were defined and shown
to be proportional to the derivatives with respect to
frequency of the dispersion law of the medium. This theory
was shown to be especially applicable when the fundamental
solution of the dispersion law was expressed in the per-
turbation expansion introduced by Perez and Uhrig (9, 17).
Kunaish (18) used time, energy, and space dependent
neutron thermalization theory in the consistent Pl approx-
imation to develop a theory for comparison with neutron
wave and die-away experiments. The non-Maxwellian component
of the energy dependence was expanded in a Taylor series,
a heavy-gas scattering kernel was used, and the results pro-
vided an accurate description of the neutron wave experi-
ment for source frequencies up to 500 ops.
CHAPTER III
THE BOLTZPIANN EQUATION IN THE DIFFUSION
APPROXIMATION WITHIN THE FRAMEWORK
OF NEUTRON THERMALIZATION THEORY
Presentation of the Theoretical Model
and the Method Used for its Solut-ion
Introduction
Placement of a sinusoidally modulated or pulsed
source of thermal neutrons at one boundary of a moderating
medium generates thermal neutron waves which propagate from
the source with frequency and space dependent amplitude
and phase. The positional variation of this frequency
dependent thermal neutron flux is what one measures experi-
mentally, and it will receive the principal emphasis in
the theory which follows. Attention will be given to the
space, energy,and time or frequency dependent thermal
neutron flux in a moderating medium within the framework
of neutron thermalization theory. The main objective will
be to calculate the complex inverse relaxation length as a
function of frequency for the fundamental space and energy
mode.
The energy distribution of thermal neutrons in a
finite moderating medium deviates from a Maxwellian distri-
bution for several reasons. The preferential leakage of
higher energy neutrons from the assembly produces a
"cooling" effect upon the spectrum, and the preferential
absorption of slower neutrons causes a "hardening" effect.
The prediction of this spectrum is further complicated by
the thermal velocity of the scattering nuclei and the
chemical binding and crystalline effects in the moderator.
The probability of a scattering interaction between
neutrons and the moderator nuclei is expressed by the
differential scattering cross section Es(E',E) for energy
transfer from E' to E. It is the essential physical
parameter in the thermal region, and the prediction of
transient and asymptotic phenomena can be only as accurate
as this scattering kernel. For this reason, care has been
taken to select as accurate a scattering kernel as possible.
The more accurate scattering kernels are much too complicated
to be considered analytically; and when using them, one must
employ numerical methods.
The Boltzmann Equation in the Diffusion Approximation
To investigate the effect of the energy spectrum upon
thermal neutron wave propagation, one can use the Boltzmann
equation in the diffusion approximation. For a finite,
isotropic, homogeneous, and non-multiplying medium it is
expressed as
v at -D(E)V24(r E t)) [L (E) + Eg(E)]M(r,E,t) +
3[Et(E) Esi 3 ~
o
Es(E',E)4(r,E',t)dE'
(1)
where
t = the time variable
r = the position vector
E = the energy variable
v = the speed variable
2 = the Laplacian operator
4-(r,E,t)drdE= the total flux of neutrons
at time t in the volume
element d? about point 3
which have energies between
E and E + dE (neutrons/cm2/
sec/ev.)
Ca(E)O(r,E,t)dE = the total number of absorptions
per unit volume about r and per
unit time at time t due to
neutrons with energy between
E and E + dE
Es(E)O(r,E,t)dE = the total number of scattering
collisions per unit volume about
r and per unit time at time t due
to neutrons with energy between E
and E + dE:
CB(E',E)4(rE',t)dE'dE the total number of scattering
collisions per unit volume about
i"and per unit time at time t
due to neutrons of initial energy
between E' and E' + dE' and of
final energy between E and E + dE
D(E) = the energy dependent diffusion
coefficient (cm.)
Within the diffusion approximation, D(E) is given by the
expression
D(E)=
(2)
where
CtEt Cs( a ,E (3)
and
Csl(E)= the second Legendre moment of the
angular dependent scattering cross
section
TEs i II n ~~s(E0*G')d(n*G')
-1
From the definitions of cs(E) and Es(E',E) it is evident
that
Cs(E) = I Zs(E,E')dE' (4)
In the experiment the source is placed near the
exposed, z=0 face of a long rectangular stack of cadmium-
covered moderator. Thus homogeneous boundary conditions can
be applied at extrapolated boundaries of the x and y variables,
the source condition is applied at z=0, and the flux can be
specified to vanish at an infinite value for the z variable
(i.e., end effects are considered negligible since measure-
ments will not be made at the far end of the pile).
Specifically, the boundary conditions that should be applied
to equation (1) for the solution of this problem are
4(Sa,y,z,E,t) = 0 (5)
$(x, Sb,z,E,t) = 0 (6)
Q(x,y,=,E,t) = 0 (7)
and
32
-D(E -- (x,y,z,E,t) =; 1/2 S(x,y,E,t) (8)
where S(x,y,E,t) is the space, energy, and time dependent
source of thermal neutrons.
Under these conditions, it will be advantageous to
expand the flux in an orthogonal set of functions in the
x and y variables. In view of the linear character of the
problem and since our main interest is focused upon the
frequency dependence of the flux, it can also be separated
into a steady part(due to the constant component of the
source) and an oscillatory part. Considering only the
oscillatory component of the solution to equation (1), we
can expand it in the form
O(r,E,t) =Re [ Z tem(x~y)q (zEI) exp(jwt)]
(9)
where 40,m(x,y) satisfies the Helmholtz equation with
homogeneous boundary conditions, which is expressed as
(V2 + B 2 ) 4 (xIy) =
x,y Itn tm (10)
with
B2 = the transverse buckling of the e,m
Iit,m
spatial mode and
Vxy the Laplacian operator for the x and y
variables.
The substitution of equation (9) into equation (1)
with the use of equation (10), multiplication of the result
by 4p,k(x~y), and integration over the x and y variables
results in
E-()Ip,k aE sc(E) + D(E) -2 3 mp,k(Z'E) +
+ I Es(E' ,E)6plk(Z'E')dE' = 0 (11)
In the experiment, one would like to measure only the
fundamental spatial mode so it is of major concern in the
theory. Since the equations for the various spatial modes
differ only by the magnitude of the transverse buckling,
one can, without serious loss in generality, remove the
subscripts. The symbols B~ and g(zE) are now understood
to represent the transverse buckling and the flux associated
with the fundamental spatial mode, respectively. It is also
convenient at this point to express equation (11) in a more
compact, operator formalism by use of the substitution*
* The author wishes to express his thanks to Dr. Frank J.
Munno for helpful discussions in connection with this
substitution.
~VNE: 7/DE~
~JD(E)(12)
where M(E), the Maxwellian distribution, is expressed as
= C(E/kTj12
M(E)
exp(-E/kT)
(13)
By substituting equation~S2i (1 at(2) into equation (11) and
dividing the result by vM(E)D(E), one obtains
C-2 8,E LS(E) + 92
D(E) D(E) az2
vD(E)
I (EE)~zE) vM(E)
w S(E,)(~
+ /
dE' 0
(14)
a2
~-B2 + --
Saz
- ]~X(z,E) + PX(z,E) = 0
vD(L)
(15)
where
E (E',E) v'M(E')'
PX(z,E) = I dE'(dE(L=--~--
o () Z())
(E'-E)[--- + ]} x~z,E')
D(E) D(E)
(16)
Solution of the Eigenvalue Problem Associated with the
space and energy uependent doiramann Up~erator
It is instructive to consider equation (15) for the
case when the transverse leakage and frequency of excitation
are both zero (i.e., a semi-infinite plane with a constant
source). The incentive for this deviation comes from the
work of Ohanian (19) and Ohanian and Daitch (20) in their
theoretical investigation of the die-away experiment. We
shall assume that the solution to the equation thus formed
may be written as
XazE) = a Jl(E) exp[-p, z] (17)
n=1
with the understanding that the above expansion is convergent
only if the qn(E) form a complete set. Substitution of
the above expansion into equation (15) with both B~ and
w equal to zero results in
i apno n(E)exp(-pnOz Pnn(E)exp(-p oz)
n=1 n=1 (18)
If the initial assumption of completeness is correct, the
~n(E) will be linearly independent. Therefore, equation (18)
must hold separately for each value of "n" and reduces to
-no n,(E) =Pn (19)
Thu pn and 9 (E) are the eigenvalues and corresponding
eigenfunctions of the P operator. An inspection of this
operator with reference to the principle of detailed balance,
v'Cs(E',E)M(E')dE'dE = vLs(E,E')M(E)dEdE'
(20)
reveals that it is both real and symmetric. Thus its eigen-
functions must be real, satisfy the orthogonality condition,
I Q (E)9m(E)dE = 0 for n/m, (21)
and be linearly independent (21). This latter characteristic
assures that the expansion expressed by equation (17) and
the subsequent manipulations were valid.
The facility of the substitution expressed by
equation (12) and of the investigation of the semi-infinite,
zero-frequency problem now becomes evident. If one expands
the energy dependence of equation (15) in terms of the
eigenfunctions obtained above, equation (19) can be used
to exchange for the complicated P n,(E) expression its simple
equivalent -p2 9 (E). Thus the solution of equation (15) by
no n
this method is not directly complicated by the consideration
of realistic and therefore involved energy exchange kernels
since more difficult kernelS Simply result in different
eigenfunctions for use in the expansiOn. Or, expressed
differently, the complicated characteristics of the energy
exchange kernel are carried to the solution of equation (15)
by the eigenfunctions and eigenvalues obtained from the P
operator. Moreover, the homogeneous equations and secular
determinant resulting from the above mentioned expansion
will be much simpler than that obtained by using other
orthogonal functions which do not bear as close a relation-
ship to the actual solution.
The eigenfunctions and their corresponding eigen-
values also have physical significance. They are the
thermalization theory solution to the exponential experi-
ment in a semi-infinite plane. The fl(E) function is the
fundamental energy mode which is exponentially attenuated
with distance from the source, and its corresponding eigen-
value is the square of the asymptotic inverse relaxation
length. The one-group diffusion theory equivalent to this
eigenvalue is 1/L2. Thus n,(E) and pno represent, respec-
tively, the eigenfunctions and eigenvalues of the space and
energy dependent Boltzmann operator and are what one en-
counters when performing the exponential experiment in a
very large moderating medium. More important physical
significance can be assigned to these eigenvalues in the
context of this work when one recalls that the exponential
experiment can be considered a neutron wave experiment
?: This argument is based upon the assumption that the
difficulty involved in the solution to the eigenvalue
problem is not dependent upon the complexity of the
scattering kernel. This is, in general, the case since,
as explained in the second part of this chapter, a
discrete energy representation of the eigenvalue problem
is solved on a computer. Thus more involved scattering
kernels simply result in different numbers being read
as input to the code.
performed at zero frequency. From this viewpoint it becomes
apparent that the pno are the zero intercepts of the squared
complex inverse relaxation lengths plotted as a function of
frequency. A knowledge of these zero intercepts is useful
because their closeness to the fundamental eigenvalue is
quite indicative of the amount of higher-energy-mode con-
tamination predicted by the particular scattering kernel
under investigation.
Further insight into the nature of these eigen-
functions and eigenvalues can be obtained by a more detailed
investigation of the eigenvalue problem represented by
equation (19). With the P operator explicitly expressed
and some terms rearranged, equation (19) becomes
C 8 2 S#(E')d
D(E) ~"o nE ME1DE
(22)
For the special case in which Ea(E) = 0, the lowest eigen-
value is zero and its corresponding eigenfunction is
~VM(E)D(E); which is easily verified by direct substitution
and by use of the principle of detailed balance. This is
consistent with the fact that in an infinite medium with
no absorption the flux will attain an equilibrium Maxwellian
distribution in energy and a constant value in space. From
the arrangement of equation (22) it is evident that it will
be satisfied by two types of solutions:
1. P2o C Z (E)/D(E) for any value of E. This
results in discrete eigenvalues and the more common
continuous eigenfunctions as solutions.
2. pno = Zt(E)/D(E). In this case the expression
enclosed in brackets is zero. The right-hand side of
equation (22) is not necessarily zero, however, so that if
the eigenfunctions are to satisfy the equation in this
circumstance they must exhibit in part a singular behavior.
Moreover, condition (2) can be satisfied for any value of
the continuous variable E under the condition p2
no
[E (E)/D(E)]MIN. Thus above this limit point a continuous
eigenvalue spectrum is obtained.
Expansion of the FreQuency Dependent Eauation in
agen unctions of the pace and nery De endent
Boltzmann Operator
The energy dependence of the flux can now be developed
in terms of the eigenfunctions mentioned above. An exponential
dependence upon the z variable will be assumed so that the
solution to equation (15) can be expressed in the form
X~,E = AJl,(E)exp(-pz) (23)
n=1
*A more rigorous expression for this expansion would be
x(z,E) = IA(d) (d)(E +I A A(c)( 2) >(C)(E, 2 )dp2)exp[-pz]
n Jn o 0 0
n=1 CC (E)/D(E)]MIN
where the n(E) satisfy the eigenvalue equation, and they
have been normalized so that
S,( 90 m(E)dE = 6mn
o (24)
Substitution of equation (23) into equation (15) results
18
n 1 [-B+ p2 A (exp[-pz])q(EE) +
+ P(exp[-pz])P9 (E) = (25)
n=1
which, when multiplied by 9m(E) and integrated over all
energy, yields
C(p2 B2 -p2 )6 L jw]A =
I no mn mn n
n=1
m = 1,2,3,*** (26)
where the superecripts "d" and "c" represent, respectively,
the discrete and continuum sets. However, in the actual
solution a discrete energy formalism is used and calcula-
tions are performed only at a finite number of discrete
energy values. Eigenvalues in the continuum, which
actually are a function of energy, and their corresponding
eigenfunctions are obtained only at particular energy
points. Thus equation (23) correctly expresses the actual
expansion used, which is the discrete approximation to the
more rigorous and less convenient representation.
where use has been made of equations (19) and (24) and
Lmn m 9(E) [1/vD(E)] n(E)dE
o (27)
Equation (26) expresses in a compact form a set of linear,
homogeneous, algebraic equations for the modal amplitudes
An and the eigenvalues p2. One is able to solve the set
only for the relative magnitudes of the An, and the number
of eigenvalues which allow a solution will be equal to the
number of equations considered. Also there will be a
different set of ratios for each eigenvalue. The compatibility
condition for the above set of homogeneous equations yields
the secular determinant or dispersion law
S(B2 p2 )6 -jwb 6 (28)
no mn mn
with
B2 = p2 2 (29)
which can be used to determine the complex eigenvalues that
will allow solutions for equation (26). These eigenvalues
are the squared complex inverse relaxation lengths which
are to be compared with the experiment. It is satisfying
to note that for w=0 the secular determinant, equation (28),
reduces to a diagonal determinant yielding the p2 for the
no
B2 eigenValues.
Determination of the Fundamental Complex Inverse
Rel axati on Le~ngth from the Complex Homogenneous Equations
The perturbation method developed by Perez and
Uhrig (9) and systematized by Kunaish (18) will be used
to develop each of the eigenvalues and modal amplitudes
in a power series expansion of the frequency variable.
This power series representation is used because it allows
one to compare very easily various theoretical models
with one another and with the experimental results. Once
the coefficients are calculated, the eigenvalue can be
determined for any frequency within the radius of con-
vergence of the series. To initiate this expansion, the
first equation and the first term of each of the other
equations in the homogeneous set represented by equation
(26) are explicitly expressed. Considering the first "N"
equations of this set, one obtains
(B2 P10) L11hw LlnjwSn = 0
n=2
for m = 1 (30)
-Lmlj C (B2 80 mn mnjw ) Sn = 0
n=2
for m = 2,3,4,*** (31)
where each equation has been divided by Al so that the set
could be solved for the ratios
S, = A /Al n = 2,3,4,*** (32)
and for the eigenvalue B2. In this method, B2 and each
of its corresponding ratios, S are expanded in a Maclaurin's
series. Thus we let
I B(K) jwK
B2 C
K=0 K! (33)
and
Sn = (34)
K=0 K!
where
(K)
B B2
aK ju (35)
and
(K) K
naK je=0.. (36)
The Kth coefficient of B2 and those of each of the S, are
obtained from equations (30) and (31), respectively, after
differentiating them "K" times with respect to jw and
setting jw equal to zero.
It should be emphasized at this point that we are
presently expanding only the fundamental eigenvalue along
with its corresponding ratios. However, it will be shown
that the other N-1 eigenvalues and their ratios are obtained
in exactly the same way after rearranging the initial
equations and redefining some terms.
By letting jw=0 in equation (30) and equation (31),
it can easily be determined that
p2
(37)
and
(0)
S, = 0
n = 2,3,4,***
(38)
Using the identity
A(x)B(x)] =
ax
I
IK) AK(x)B(I-K)(x)
K=0
(39)
where
7 1!
KK!IK !
(40)
one can differentiate equation (30) and equation (31) "I"
times to obtain at jw=0
I I
- ull- ) 1 1
K=0 n=2
B
T ) (I-K)=0
K In K1 n
- 2 6
10 10
-L 6
11 II
(41)
N
L 6 -u(I-2) ( )(B(1)6 -L )S(I1
ml Il 1 an mn n
n=2
N
+ r
I
I j( (B(K) p2 6 )6 LS(I-K)= 0
K ~no KO an an K
-L 6
n=2 K=0
m = 2,3,L4,***
(42)
where u(I-2) is the unit-step function. By evaluating the
summation over "K" in equation (41) and rearranging terms,
one obtains
N
+ u(I-2) C : Ll Sn(I-1)
n=2
B
= 2 4
10 IO
+ L 4
11 Il
(43)
which can be used to determine each 8 if the S(I1
are known. Equation (42) can be simplified by separating
out from the rest of the series in the "K" summation the
terms for which K=0 and K=1. By performing this manipulation,
(I)
extracting the S coefficient from the term for which
(I)
K=0, and rearranging terms, an expression for the Sm
ratios is obtained
(I)
S
m
(P2 -p2 )
mo 10
N I-1
u(I-3) (K) (K)S~(I-K)
n=2 K=2
; m = 2,3,4,***
(P2 p2 )
mo 10
The summation over "K" in the above expression is extended
(0)
only to "I-1" since Sm is zero, and the step functions
are introduced because some terms were extracted from the
summations.
The procedure used for evaluating the B(I and S(I
is initiated by first determining B(1) and S(1) from
equations (43) and (44), respectively. Knowing these terms,
B()can be calculated from equation (43) and placed in
equation (lt4) in order to find the S(2). This process of
alternately using the two expressions to first calculate
BI then when all the B(j and S(J are known for
J*
*
Then using the B(I in equation (33), the frequency dependent
complex inverse relaxation length for the fundamental mode
can be evaluated.
It may happen, however, that the series does not
converge very rapidly. In fact, this is the case for one
of the scattering kernels selected. For this situation, a
more convenient series can be obtained by a change of
variable of the form
jw = jxex
(45)
One can then derive a power series expansion for B2 in
terms of X. Recognizing that
(46)
andusing the procedure as explained with the w expansion,
one can obtain expressions for the B(I and the S(I for
the power series expansions in the x variable. As in the
w expansion
(0) 2
B =P
10 (47)
(0)
S
m
and all of the
are zero. For I = 1, one finds
j11
(48)
( 9)
(1)
Sm =jLml/ pmo 10)
and, in general, for I i 2
m =2,3,4,**
N I-1
L E [K) jKLlnS~(I-K)
n=2 K=1
8 )
= jILll -
(50)
dn
ix-x (jXX) = n
drn|70
and I-1 N I-1
jI~ml KB(K)S'I-K) K)jKLmnS(I
K=1 n=2 K=1
S =
(P2 p2
mo 10
m = 2,3,4,*** (51)
In the determination of p2 as a function of frequency, both
of these power series expansions will be used.
An alternative to the series representation is to
solve the complex secular determinant directly for several
values of w. For reasons previously discussed, this method
is not as appealing as the ones derived above; and it will
be used only as a last resort, that is, for cases in which
both of the series converge slowly or not at all.
Determination of the Higher Order Complex
Inverse Relaxation Lengths
It was stated that the equations obtained for the
(I) (I)
B and the Sm were derived to evaluate only the funda-
mental p2. To calculate the other eigenvalues, one can
rearrange the equations expressed by equation (26) so that
the same expressions can be used to calculate coefficients
associated with all of the eigenvalues. To explain this
procedure, consider the case where m=3 in equation (26).
Allowing n to range from one to three, one obtains
C~2 B2 102 11 jA1 L12j 2 L13 jw3
(52)
- L21jWA1 + ~2 _: p20)0 L22jU1A2 L23jWA3 =0
(53)
-L jWA-L 2jwA2 C(p2 B 2 p2) Ljw]A = 0
Let us specify that the eigenvalue associated with the
second equation is desired (i.e., the eigenvalue whose
zero intercept is p20+ B ), If we exchange equations (53)
and (52) and rearrange terms in each of the equations, the
result is
C(p2 :2 2 :) L2jW]A2 LljWA1 L 3j.A3 = 0
(55)
-L 2jWA + Ep2 2 102 )-L jm)A L 3jwA =
(56)
-L 2j A2 LljWA1 [(pZ B2 2 ) 0, L3jWlA = 0
(57)
An inspection of these last three equations and the ones
from which they were derived indicates that both sets are
of exactly the same form since they differ only in sub-
scripts. Thus the same expressions used for the solution
of the fundamental eigenvalue are also valid for the second
if care is taken to correctly rearrange the L and p2
mn no
matrix elements and to assign subscripts to them correspond-
ing to their new positions before applying the equations.
By analogy, the arrangement of terms for calculating the
nth eigenvalue and its ratios is obtained by exchanging
the first and the nth equations and then exchanging the
first and nth terms in all of the equations.
Numerical Calculations Associated
with the 'Iheory
Introduction
The derivations thus far presented are founded upon
the theoretical formalizations of Perez and 'Ihric (9) and
Ohanian and Daitch (20). The work of the former authors
indicated the method by which, beginning with the Boltzmann
equation in the diffusion approximation, one could formulate
a power series expansion for the complex inverse relaxation
length. They used Associated Laguerre Polynomials in ex-
panding the energy dependence, whereas here the eigenfunctions
of the space and energy dependent Boltzmann operator form
the orthogonal functions for this expansion. This more
elegant method has practical value since the secular deter-
minant used to determine the eigenvalues, B2, has a much
simpler form and one need not calculate the complicated
matrix elements of the scattering operator. As a by-
product of this method, one obtains all of the zero inter-
cepts of the Bn ie, n
For the determination of the B2 and their
corresponding eigenfunctions, use was made of the numerical
techniques developed by Ohanian for his calculation of the
eigenfunctions of the time and energy dependent Boltzmann
operator. In fact, after minor modification, several of
the computer codes developed by Ohanian were used in this
work.
Evaluation of the Eigenvalues and Einenfunctions of
he Space and Energy UDeendent Boltzmann Operator
in_ the bittusiono Approximation
The Discrete Energy Formulation of the Eigenvalue
Problem -- In the discrete energy representation, equation
(15) becomes a set of "J" equations of the form
L, i) s,(E )
02 +(E.) = r K(E ,E.)" (E.)AE. 1 + (E )
noj= n (E ) D(E )
i = 1,2,3,***,J
(58)
where
K(EE i)
iVM(q)1D( )D(ED ) (59)
11=1 (60)
bEj depends upon the numerical integration scheme selected,
"J" is equal to the number of energy points, and the prime
on the summation indicates that the term involving
Es Ei'Ei) has been eliminated from the summation in the
interest of greater numerical accuracy (19). An inspection
of the above equations with all terms included reveals that
for i=j the term E (E.,E.) in the j summation in equation
(58) should be canceled by an identical term in the
expression .It was determined by Ohanian that
D(Ei)
greater accuracy could be achieved by elimination of these
terms altogether before introducing equation (58) into the
computer. This is a very fortunate development since, as
explained later in this section, it was necessary to adjust
the diagonal elements of the scattering kernel calculated
by the SUMMIT code. Since the diagonal elements are not
needed in this numerical solution to the eigenvalue equation,
this modification should not have any direct influence upon
that solution.
The set of equations represented by equation (58)
can be more simply expressed as
no n -n
where
(J9 ) is a column vector composed of the
values of the @n eigenfunctions at
each of the ] energy points,
(bE) is a diagonal matrix containing the
weights for the numerical integration,
and
(P) is a square matrix representing the
P operator.
The elements of the (P) matrix are given by the expression
J'
ijj D(E g) D(E )
(62)
We have now arrived at the equation which is to be
solved by the EIGENVALUE code. Note that the symmetry
of the P operator is preserved in the discrete representa-
tion and that since (AE) is simply a diagonal matrix to
provide weighting for the integration the numerical eigen-
value problem also contains a symmetric operator.
Evaluation of the Soatte~rin~g Kernel and F~ormulation
of the Symmetr~ized P Operator -- All of the theoretical
calculations were performed on an IBM 709 digital computer.
In essence, the implementation of the theory derived in the
preceding section was accomplished through the use of an
echelon of computer codes, beginning with one which calcu-
lated the scattering kernel and ending with the EIGENVALUE
code. Each code made use of the results of one or more of
its predecessors and, in turn, contributed to the codes
farther on in the chain. The format for the discussion in
this section can therefore be logically based upon the
flow diagram shown in Figure 1.~ In the discussion which
follows, pertinent information pertaining to the computa-
tional methods will be presented. For a complete discussion
of the workings and input-output information for the com-
puter codes used the reader is referred to the reports
submitted by the author to his committee supervisor.
One of the more important considerations in the
implementation of this numerical method was the determina-
tion of what variable to use in order to most advantageously
express the energy dependence. This problem was solved by
Ohanian (19) by the use of a variable defined as
n= (E/KT)1/2 / (1 + (E/KT)1/2) (63)
The 33 energy points, corresponding to 32 equally spaced
increments of the n variable, that he used were adopted for
this work and are listed in Table 1. Also shown in the
table are the D(E ) as obtained by the use of SUMMIT,
PINNACLE, and DIAGONAL CORRECTOR codes. For the numerical
integration encountered, except those in SUMMIT and
*This computation was done under the direct supervision of
Dr. M. J. Ohanian, with initial inspiration from Dr. R. B.
Perez. The author is grateful to Robert B. Razminas for
the use of his DIAGON~AL CORRFECTOR CODE and for his con-
tribution to the matching of several of the codes.
c,
Em30
O
O
O C
O o
,O O
v0 9$
popr 4
iJ C Cp a
E aH
GC hO O
H-(O b e
v0 ppy C
0000
HcCCOO
YA WH C)
O CJ
CH
O r
*H
to C
Icr
>O
km,
c (
5:
h
r,
W
r
rl
W
v
LI1
IABLE 1
THEC ENERGY POINTS USEO IN fHE DISCRETE ENERGY REPRESENTATION
AND THE CALCULATED VALUES OF THE DIFFUSION
CUtFFICIENI AT THESE POINTS
ENERGY POINfS
(tV.)
0.19489899E-04
0.82478099E-04
0.19665200E-03
0.37111399E-03
0.61667299E-03
0.94622400E-03
0.13752100E-02
0.19221999E-02
0.26096600E-02
0.34648699E-02
0.45211899E-02
0.58196499E-02
0.74110299E-02
0.93586800E-02
0.11742199t-01
0.14662799E-01
0.18249500E-01
0.22669999E-01
0.28143399E-01
0.34960400E-01
0.43512300t-01
0.54334600E-01
0.68173099E-01
0.86089300E-01
0.10963000E-00
0.L4111499E-00
0.L8413600E-00
0.24448000E-00
0.33191600E-00
0.46393500E-00
0.67426600E 00
0.10345899E 01
0.17191599E 01
DIFFUSION COEF.
(CM.)
0.90540472E
0.16644209E
0.28357907E
0.38656164E
0.50007596E
0.61531631E
0.77092278E
0.51939L03E
0.70397882E
0.94295391E
0.12733748t
0.96080945E
0.95736114E
0.91105106E
0.11272739tr
0.10718941E
0.81394893t
0.79333499E
0.84320298E
0.82992L54E
0.82778580E
0.809777086
0.77186675E
0.78670610t
0.79197136E
0.77184059E
0.78138231t
0.76182067E
0.78947944E
0.91681817E
0.96415816t
0.148592776
0.15770610E
PINNACLE, Simpson's rule was used since most of the functions
to be integrated were generally smooth. Comparatively few
quantities were needed as input parameters, and the more
important ones are listed in Table 2.
Figure 1 shows the sequence of the numerical cal-
culations and the input and output of each code in the chain.
For the first task, that of calculating the scattering kernel,
two different models were used. The first, and more realistic,
kernel was formulated by Parks (22,23) and incorporated by
Bell (24) into a computer code called SUMMIT. The Pl kernel
was calculated from the output of SUMMIT by the PINNACLE
code. The code multiplied the angular dependent scattering
kernel from SUMMIT by the cosine of the scattering angle,
performed a Gaussian integration over the angle, and multi-
plied the result by an appropriate energy dependent factor.
The second scattering kernel considered was a free-gas energy-
exchange kernel obtained by the use of the MASSMK code (25).
Note that the calculations involving the MASSMK kernel
differed from those of the SUMMIT code only by the scattering
kernel Es(EiE~j). All other inputs to the OMEGA-ZERO SYMBOL
via DIAGONAL CORRECTOR (26) were the same as those obtained
from the SUMMIT kernel.
The scattering kernel calculated by SUMMIT does not
include the elastic contribution. For this reason, it was
necessary to correct the diagonal elements of this kernel
so that the scattering cross sections obtained from it were
consistent with values taken from the "barn book" (BNL-325).
Parameter Value Units
TABLE 2
PARAMETERS USED INI THE NUMERICAL
CALCULATIONS
Mass of Graphite
ao
Vo
0I (number density)
kT
12.011
90.04
2.198x105
0.0803
0.02526
amu.
sec'1
cm/sec
n/Cm3
ev.
0.027027
4.6949
barns
as FRE~E
This was accomplished in the DIAGONAL CORRECTOR code.
As shown in Figure 1, several other quantities were also
obtained from this code. The (P) matrix was calculated
using the OMEGA-ZERO SYMBOL code. This code was obtained
by modifying the SYMBOL code written by Ohanian to account
for the differences between the (P) matrix of equation (62)
and the one obtained by him for his eigenvalue problem.
The eigenvalues and eigenfunctions of the space and
energy dependent Boltzmann operator were computed by the
EIGENVALUE code using the Jacobi method of matrix diagonali-
zation (Appendix A of reference (19)), This was not the
same code as used by Ohanian but it was similar to it.~
The code used in this work was modified to accept the
diagonalization when the smallest diagonal element became
more than 1016 greater than the largest off-diagonal element.
Presentation of the Eigenfunctions and Eigenvalues -
Several of the eigenfunctions and all of the eigenvalues
computed by the method outlined above are presented,
respectively, in Figure 2 and Table 3 for the SUMMIT
kernel and Figure 3 and Table 3 for the MASSMK kernel.
Some problems were encountered in drawing the eigenfunctions
* The code was obtained by Dr. Ohanian and converted
for use on the University of Florida IBM 709 computer
by the author.
(2)
(31)
(33)x1/2
(11)
(31
(1)
3
(2)
Log of Energy
(pev.)
(31)
( 11)
( 21)
ig. 2. Several Figenfunctions of the Space and Energy Dependent
Boltzmann Operator Using the SUMMIT Kernel.
g,. 3. Several Eigenfunctions of the Space and Energy Dependent
Boltzmann Operator Using the MASSMK Kernel.
TABLE 3
EIGENVALUES OF THE SPACE AND ENERGY DEPENDEN1
BULFLMANN OPERATOR USING THE SUMMIT KERNEL
AND THE MASSMK KERNcL
SUMMIT
MASSMK
0.42994279E-03
0.30594949E-02
0.143L0712E-01
0.63025152E-01
0.80414867E-01
0.10734916E-00
0.98474810E-01
0.12307265E-00
0.14447561E-00
0.17164934E-00
0.22728613t-00
0.24452724E-00
0.26831973E-00
0.L6281100E-00
0.30214746t-00
0.34459867E-00
0.32484750E-00
0.35953601E-00
0.39137711E-00
0.40926518E000
0.38388654E-00
0.43537266E-00
0.44520015E-00
0.43975324E-00
0.44780591E-00
0.44983647E-00
0.47610G663-00
0.4824065LE-00
0.50220024E 00
0.71235961E 00
0.10231198C OL
0.12931190E 01
0.48972284E 01
0.42844241E-03
0.30875468E-02
0.37078650E-02
0.44850244E-02
0.65955998E-02
0.10662339E-01
0.14949699E-01
0.17229941E-01
0.20076749E-01
0.19529168E-01
0.20320657E-01
0.21378068E-01
0.22301081E-01
0.24511521E-01
0.2775267tE-01
0.29222453E-01
0.35146919E-01
0.40474408E-01
0.366L9057E-01
0.43478405E-01
0.50642820E-01
0.49833021E-01
0.55298758E-01
0.65017583E-01
0.77065729E-01
0.81092229E-01
0.91628064E-01
0.11717452E-00
0.12339253t-00
0.10396006E-00
0.13176729E-00
0.14078926E-00
0.1900513'jE-00
because they can have more than one singularity and in the
discrete representation it is difficult to distinguish
between multiple discontinuities and oscillations. The
possibility for multiple disconitinuities is attributed to
the many fluctuations in Eg(E)/D(E). This function can be
identical to a particular eigenvalue for several energy
values thereby causing the corresponding eigenfunction to
have several discontinuities. From the characteristics of
the eigenfunctions, it appears that only the fundamental
eigenvalue is discrete; and the remaining ones represent
particular eigenvalues in the continuum. As shown in Table
3, beyond p20 the larger eigenvalues associated with the
SUMMIT kernel are closer to the fundamental eigenvalues than
are the ones from the MASSMK kernel.
The Energy Expansion of the Freauency Dependent
boltzmann Lquation inth theDifusion Approximation
The expansion of the energy dependence of equation
(15) in terms of the eigenfunctions of the P Operator re-
quires the evaluation of only one matrix element, Lmns
expressed by equation (27). In the discrete energy repre-
sentation this integration becomes
L, = mi1 nE)g
i=1 v(Ei)D(Ei) (64)
where AEi is the weighting factor for the numerical integration
A L-MATRIX ELEMENTS code was written which evaluated these
matrix elements using both Simpson's rule and Weddle's rule
(27) and which also checked the orthonormality condition.
The numerical integration were evaluated using the n variable
since the eigenfunctions obtained from the EIGENVALUE code
were expressed in this variable and the discrete energy
representation was based upon the selection of equal n
increments for the evaluation of the energy points. It was
felt that since each of the continuum eigenfunctions exhibited
at least one discontinuity, and since D(Ei) is not a smooth
function, it would be necessary to apply very elaborate
integration schemes in order to obtain accurate results. If
any integration other than Simpson's rule was used, it was
found that for some of the eigenfunctions the orthonormality
check was not satisfied (i.e., the integration over energy of
all ~m i) times 9 (E ) was not necessarily much less than
the integration of all [$m(E )]2) ggypp fac isntsrrs
since the eigenfunctions obtained from the EIGENVALL'E code
were obtained with a Simpson's rule weighting function.
For this reason, it was felt that continuity of the numerical
methods, that is, the use of Simpson's rule throughout, would
result in the more consistent set of Lmn elements. Some
justification for this reasoning is found in the work of
Ohanian (19) since a similar integration was required in his
calculation of the cooling coefficient for water and accurate
results were obtained."
A pragmatic consideration of this calculation provides
some satisfaction since its use resulted in theoretical
curves which were in good agreement with the experiment.
Also the diagonal elements of the Lan matrix were, in
general, much larger than the off-diagonal elements since
in the evaluation of the former, the large discontinuity
of both eigenfunctions occurred at the same energy point.
Thus after a careful investigation of the integrands arising
in the evaluation of the Lmn matrix elements, some justifi-
cation can be obtained for their accuracy. But, in the
general case, the numerical integration of such "wild"
functions by Simpson's rule cannot be defended.
Evaluation and Presentation of the Complex
Inverse Relaxarion Lengrhs
With the Lmn matrix elements determined, one can
apply equations (43) and (4R) to evaluate the coefficients
of the expansion of p2 in powers of (jo) and equations (50)
and (51) to evaluate the same expansion in powers of x.
In cases where neither series converges rapidly, the
*Also in calculations: performed by Razmr~inas- mid the author
under the direction of Ohanian, the cooling coefficient
for graphite was calculated by Ohanian's method. It
agreed quite well with the one obtained by Honeck using
the same kernel.
eigenvalues were evaluated directly from the secular
determinant at particular values of frequency. Computer
codes were written to calculate the coefficients for both
the jw expansion and the X expansion of any of the 33 p2
along with their corresponding ratios, S If more than
a few terms were used in the jw expansion, however, it
was necessary to divide the equations by L11 and derive
expressions for the expansion in powers of Llljw. This
was because at high frequencies, wn became larger than
1038, the largest number acceptable to the computer. This
problem did not arise in the X expansion since for a
frequency of about 1000 ops X need be only about 7.
Problems were also encountered in the coefficients. For
the jw and Llljw expansions, the series could not be
carried past the fifteenth term since,in the former case,
underflow occurred in the computer and in the latter case,
overflow occurred. In the x expansion the coefficients did
not change as rapidly and no underflow or overflow problem
occurred; but since factorial 33 is nearly 1038, the series
was cut off at the 33rd term. This, however, was not the
limiting factor in the X expansion for it appeared that
after the 20th term round off errors built up rapidly.
Shown in Figure Ir and Figure 5 are the frequency
dependent real and imaginary components of the complex
inverse relaxation length for the fundamental energy mode
If m I I \ I -I I j
ce cu cy e- e- 0
E*
a,C
co
CC
O <
,OQ
Ca
XC,
00
Hrr
X *Hl
CLE
Ek
OW
COk
arO
O
C *,
Ok
MO r
*
0 *j
0
*H
Q)
o 'c- '
o ri
o OO
hC
o COl
SXO
0t 0
o OCI
O CO
oc
on o
E U
0~ UH
a .sa
E 0*
( O/pm) 3
and the first four energy harmonics as generated by using
the MASSMK kernel. Presented in Figure 6 and Figure 7
are these same eigenvalues up to the third harmonic ob-
tained by use of the SUMMIT code. Figure 8 and Figure 9
show the fundamental and first two space harmonics for
the MASSMK calculation. The transverse bucklings used in
this last set of curves are indicated on the figures, and
for all other curves a value of 5.96 x 104 om2 wsue
since this was the calculated value for the graphite pile
investigated. The omega expansion was used in the genera-
tion of all the MASSMK curves, and the method presented
previously for the generation of higher energy harmonics
was successfully applied with this scattering kernel. Six
terms were sufficient for convergence in this case. For
the SUMMIIT kernel, however, the omega expansion did not
converge when as many as 15 terms were used; and for this
scattering kernel the X expansion and the direct numerical
solution of the complex secular determinant yielded the
eigenvalues*" As shown in Figures 6 and 7, the 18-term
X expansion appears to have converged since agreement exists
between it and the direct solution for the fundamental
" The code used for this calculation was obtained by
Victor Cain from the SHARE Program Library. The code,
called CIGENVALUES 0" COMPLEX YATRICES, has as its
identification number NU EIG4.
r*
tH
Q0
o
to
COv
0 *H
m CF0
o cXc
to HQ
o HO
o *H
Hr XG,
00c
xE
o ,
to a 0
C\C
Le CO
OO
E v
OOc
N *
*
r( rl
(~,U13) D
00+
o
a Ukb
o OC
CD
CO
O cl
o *
o clm
o a Ur
O 4 C
m a
o mC
o FOE
X>
c 40
0 0 EC
rl v O 6
OF
CFO
Ocr
o v C
UN
o do~-
*H O
Hi C/3
(ur~/pe;r>
S,C
t~e
CC
O *H
c u
HQ
X0
oO:
H *H
Ac
O -
Q~O
[2 *H
3
*Hr
0 0
a
0(
v, h
o k C
o OC-
O I
SbC
,CC:
C3
0 J
o C
c CO
CO
O Xa
o Eu
o Oa
rl Ua
,CO
OIJ
AC *r
o OO
aC (
HE~-
a
Co
m rl
H *c
b^C a
*HC
O to
o
Hl 0
(ur~/peJ) t,
eigenvalue. The higher energy mode eigenvalues for the
SUMMIT kernel were calculated by the direct method.
Some concern can be caused by a careful inspection
of the higher energy modes. For higher frequencies, these
modes appear to be less attenuated than the fundamental
mode; and, therefore, modal contamination is predicted.
The less accurate MASSMK kernel does not predict as serious
energy mode contamination. One must recall, however, that
in order for these higher energy modes to propagate, they
must be excited by the source or the medium; and even in
the finite, absorbing medium to be investigated they will
account for only a small percentage of the energy spectrum
compared with the fundamental Maxwellian distribution.
Also Kunaish has shown that for the heavy mass scattering
kernel the total energy spectrum, calculated from as many
as nine of his energy harmonics, will propagate with a
complex inverse relaxation length which is within 1
per cent of that of the fundamental.
From the curves showing the higher spatial mode
eigenvalues, it can be concluded that for this experiment
they will be of secondary importance in comparison with the
higher energy modes. Calculations by Perez (9) and Kunaish
(18) predict that higher spatial contamination should be
reduced to a negligible amount at points farther than 20
cm. from the source. The calculations presented here
confirm these predictions.
CHAPTER IV
THE EXPERIMENT IN GRAPHITE
Introduction
The experimental techniques currently used at the
University of Florida for the production and detection of
neutron pulses and waves were developed over a period of
nearly four years through the efforts of many persons.
As a half-time graduate assistant on the neutron wave
project, the author has been intimately involved in all
phases of this developmental work. The methods for
analyzing the data have also evolved over a period of
years. The principal contribution of the author has been
in this area; for he has, under the guidance of his chair-
man and with suggestions from co-workers, written every
computer code required for the data analysis and analyzed
or co-analyzed all the data taken thus far in the neutron
wave project. Also this work represents by no means the
first, nor the last, experimental results to be derived
from the neutron wave project, as is evident from an in-
spection of the Review of the Literature.
The medium investigated in this experiment was a
rectangular pile of AGOT graphite 44.7 cm. by 71 cm. by
125 cm. covered, except on the source face, with 30 mil.
thick cadmium. The discussion of the experiment which
follows is divided into three general areas. Consideration
is given to the methods used to excite the thermal neutron
waves. Attention is then focused upon the systems employed
in their detection. The third section deals with the
analysis of data and will be given special emphasis.
The Source
The source of thermal neutron waves used in this
experiment was a Cockcroft-Walton type neutron generator
placed in a thermalizing tank as shown in Figure 10. The
neutron generator, manufactured by Texas Nuclear Corporation
of Austin, Texas, produced either a sinusoidally modulated
or a pulsed intensity of 14 Mev neutrons which were then
slowed down to thermal energy and directed to the graphite
pile by the thermalizing tank. The sinusoidally modulated
intensity of neutrons was obtained by varying the extraction
voltage of the ion source bottle, thereby controlling the
intensity of the deuterium beam which impinged upon the
tritium target. The pre-acceleration pulsing unit consists
of deflection plates, Einzel lens, and focus lens. This
unit deflects the deuterium beam to the base of the Einzel
lens during the wait time between pulses and allows it to
Fig,. 10.
Arrangement of the Experiment and Details
of the Thermalizing Tank ( from Re ference 14) .
pass on to the target for the duration of the pulse. For
the data obtained with the modulated source, the frequency
varied between 50 and 1226 cps, the ratio of maximum to
minimum neutron count rate was roughly 9 to 1, and harmonic
and 60 cycle contamination were quite small. For data taken
in the pulsing mode, the maximum to minimum thermal neutron
yield was about 1000 to 1; and the total pulse width was
about 1 millisecond at the a position in the graphite which
was closest to the source.~ More detailed discussions of
the sinusoidal and pulsed mode are given, respectively, in
References (8) and (14).
As previously mentioned, the thermalizing tank was
used to thermalize the 14 Mev neutron waves produced by the
neutron generator and then transport them to the external
port of the tank where the graphite medium was placed. This
was accomplished by placing light water, steel bricks, lead,
and graphite in the tank and facing the tank with a poured-
paraffin section. By systematically varying the dimensions
of the various materials and carving out a 8 cm. by 8 cm.
square hole in the center of the paraffin section, the
*' It should be mentioned that the attainment of such a
high-quality source was due primarily to the efforts
of Robert H. Hartley and George W. Fogle.
uncovered-to-covered count rate of a 1 atm. He3 detector
was optimized at 50 to 1 at the source interface of the
graphite. The initial thinking associated with this
thermalizing tank is presented in Reference (8). The
unwanted, but ever present, fast background emitted from
the tank was eliminated experimentally as discussed below
and in the section dealing with the data analysis.
The Accumulation of Data
It is fitting at this point to indicate precisely
what data are to be presented and how it was recorded. The
main objective of the neutron wave experiment is to measure
the frequency dependent complex inverse relaxation length,
p, of the fundamental space and energy mode. The most
obvious method to accomplish this objective is to modulate
Isinusoidally the neutron current which enters one boundary,
z=0, of the medium under investigation and to measure the
amplitude and phase of this neutron wave as it propagates
from this source. The slope of the semilog plot of amplitude
versus z is the real part of p, or the attenuation factor,
a; while a plot of the phase of the wave versus z is linear
with a slope equal to the imaginary part of p, or the phase
shift factor, 5. The source is modulated for as many fre-
quencies as desired by the experimenter. The process out-
lined above was used to obtain the wave data presented in
this work. Still another way of obtaining the same information
is to record the dispersion of a pulse of thermal neutrons as
it propagates from the source boundary. The amplitude and
phase of any frequency component of this pulse at a particu-
lar z position is obtained by numerically performing a Fourier
transform of the measured, time response at that position.
These methods are not new. The analogy in electrical theory
to the discussion above is that the frequency characteristics
of a "black box" can be determined by either directly measur-
ing its response at various frequencies or by performing a
Fourier transform of its unit impulse response. One method
is not necessarily superior to the other for all applications,
but in this measurement it was found convenient to obtain
most of the data by exciting the neutron waves with thermal
pulses.
The pulse propagation data were collected by Hartley,
and first appeared in his master's thesis (14). The
purpose of his investigation was accomplished by considering
data in the time domain~ only. In this work the analysis of
Hlartley's data is extended beyond that which he completed
since Fourier transforms are performed under the Fguidance
of Reference (16) and the complex inverse relaxation length
is obtained. Actually, in the final comparison with the
theory, the wave data and pulse propagation data are combined
as discussed in the next section.
From the above discussion it is evident that it will
be necessary to discuss two data accumulation systems, that
associated with a sine wave input and the one used by
Hartley in his pulse propagation experiment. Before con-
sidering these systems individually, however, it is con-
venient to present those characteristics which were common
to both experimental set-ups.
In each of the experiments use was made of a movable
detector which could be located at several z positions, and
a stationary, standard detector placed in the thermalizing
tank. To facilitate accurate placement of the movable detector
and cause as little perturbation as possible, length-wise
holes were bored in several of the graphite stringers to
accommodate the detector and its preamplifier. For measure-
ments at a particular z position, a detector-filled stringer
was inserted into the graphite pile so that its center was
aligned with the transverse center of the pile and its length
was in the transverse plane. By using several of these
stringers with holes bored at different distances from their
front face, the detector could be moved from the source in
about 1 cm. increments from 1 cm. to 72 cm. Two measurements
were made at each position of the movable detector; one
where both epicadmium and thermal neutrons entered the
moderating medium, the uncovered run, and a second where
thermal neutron penetration was blocked with a cadmium shutter,
the covered run. The standard detector was used to correct
for source variations and to facilitate subtraction of the
epicadmium background from the movable detector data of
the uncovered run, thereby obtaining the thermal wave. The
same source (i.e., neutron generator and thermalizing tank)
was used in both experiments, and changing from the wave
experiment to the pulse experiment was accomplished by
changing accordingly the mode of operation of the neutron
generator. The same detector systems were also common to
both measurements. They consisted of He3 detectors and Nuclear
Chicago Corporation oreamnplifiers and scalers. The system
dead time was about 3 microseconds,which was found to be
adequate for our purposes. The arrangement of electronics
beyond the scalers was unique for each of the two experi-
ments and they will now be discussed separately.
Electronic Eguipment Associated with the Wave Experiment
A block diagram of the electronics used in the wave
experiment is shown in Figure 11. The main component of this
system is the multi-channel analyzer manufactured by Technical
Measurement Corporation. A modified 212 plug-in-unit is
used in the multi-channel analyzer so that it stores the
input pulses in a channel corresponding to their time of
arrival after an initial trigger. The analysis time of
each of the channels is stepwise variable from 10 to 2,560
microseconds so that the accumulation may cover a wide
range in frequency. The remaining circuitry is arranged so
1
Fig. 11. A Block Diagram of the Electronic Equipment Used
in the Wave Experiment.
Movable
Graphite Pile
Cadmium Sheet
that for two periods of the input sine wave the output
of the standard detector and that of the movable detector
are recorded, respectively, in the first half-memory and
second half-memory of the multi-channel analyzer. This is
accomplished by alternately allowing the pulses from each
detector to reach the 212 plug-in-unit, controlling their
storage in the multi-channel analyzer, and adjusting its
channel width so that the total sweep time is slightly less
than four periods of the input frequency.
A time-sequence diagram for the system is shown in
Figure 12. The HP202A function generator provides a pulse
at the positive peak of its sine wave output to the neutron
generator. This pulse train is fed to a Tektronix type
161 pulse generator which, when triggered by a sync. pulse,
produces another pulse whose width is almost as
long as two periods of the input sine wave. This output
pulse from the Tektronix type 161 pulse generator, which
is triggered by every other sync. pulse, serves two purposes.
It is used as a trigger for a Tektronix type 163 pulse gen-
erator whose output, connected to the trigger input of the
212 plug-in-unit, initiates the sweep of the multi-channel
analyzer. Thus the analyzer sweep always begins at the
same point in the cycle of the input wave. The pulse train
is also connected to the external trigger input of the
Tektronix type 53/54C plug-in-unit, and it controls the
switching of this gate circuit. Other inputs to the Tek-
tronix type 53/54C plug-in-unit are the pulses from the
T
3V53/54C M~emory Sig~nal
T 2T 3T 4T
Time (Sine Wave Periods)
Cycle Time of Elec-
tronic Equipment
Period of Input
Sine Wdave
Neutron Source
Intensity
cync. Pulse from
HP202Ai Function
Generator
2T-Pulse from the
Tektronix 161 Pulse
Generator
10V
12V Trigger Tnput from
Tektronix 163 Pulse
Generator
53/54C Output Sig~ndl
Fig. 12. Time Sequence of Electrical Signals Associated
with the Electronic Equipment Used in the Wave
Experiment.
-- II (20V
Standard Movable
Detector Detector
standard detector and those from the movable detector.
The Tektronix 53/54C plug-in-unit contains two identical
amplifier channels; and it is usually used with a dual
trace oscilloscope so that two signals, one connected to
input A and the other to input B, appear on alternate
sweeps of the oscilloscope. In this application the pulses
from the two detectors are allowed to pass alternately to
the input of the 212 plug-in-unit. An alternate output
from the 53/54C plug-in-unit is connected to the EXT input
of the multi-channel analyzer, and it i9 used to control
the half-memory in which the detector pulses are being
stored. Thus the 53/54C plug-in-unit controls which
detector pulses are being recorded by the multi-channel
analyzer and also controls the half-memory in which the
pulses are being stored. This eliminates the possibility
of the storage locations for the two detectors being inter-
changed during the accumulation of data.
Referring again to Figure 12, the functioning of the
detection system during one analysiS sweep of the multi-
channel analyzer is as follows:
1. A trigger from the Tektronix type 163 pulse
generator starts the multi-channel analyzer on its analysis
sweep. The pulses from the standard detector will be stored
in its first half-memory, and then it will stop and wait
for another trigger.
2. The next trigger will start the multi-channel
analyzer again, but now the Tektronix type 53/54C plug-in-
unit is passing only the movable detector signal and holding
the multi-channel analyzer in the second half of its memory.
3. The next trigger will again start the cycle,
which is repeated until enough statistics have been obtained.
Electronic Equipment Associated with the Thermal
_Neutron Pulse Propagation ELxperiment
A block diagram of the circuitry used by Hartley
is shown in Figure 13. For his experimentzonly the integral
counts of the standard detector were needed. The output of
the movable detector was recorded in the multi-channel analyzer
after shaping in a Tektronix type 163 pulse generator. The
sweep of the multi-channel analyzer was initiated by a signal
synchronized with the start of the neutron burst from the
neutron generator.
The Data Analysis
Introduction
The focal point in the analysis of neutron wave
data is the counts per channel recorded for the movable
detector by the multi-channel analyzer. To indicate the
typical appearance of these data, the 200 ops thermal wave
is shown in Figure 14 for three of the positions used in
the wave experiment; and the thermal pulse at several
Fig.13.A Block Diagram of the Electronic Equipment Used
in the Pulse Propagation Experiment (from Refer-
ence 14).
72
00v
I I I I I
*H
OJ
0. 0
is IIF
Nv N
II II II
,C
EC r:
*~ *
I I I I I E
0( \I I N --
~O~ x u x ~auuP~3 ~ad s~uno3
positions is shown in Figure 15. In the wave experiment
these data are composed of a fundamental sine wave con-
taminated by harmonic distortion, and a Fourier series
expansion yields the amplitude and phase of the fundamental
thermal wave. The thermal pulse, on the other hand, was
considered to be composed of neutron waves in a continuous
spectrum from zero frequency up to a high frequency cutoff.
Thus the amplitude and phase of a particular frequency was
obtained by numerically Fourier transforming the pulse at
the frequency of interest. After the attainment of the
amplitude and phase points versus z position for several
frequencies, the data analysis and comparison with theory
can proceed in an identical manner for both the wave and
pulse propagation experiments. In fact, in this work the
amplitude and phase curves obtained by the pulse propagation
experiment and those from the wave experiment are normalized
together before determining their slope. In this way the
complex inverse relaxation length, p, arises from measure-
ments in both the time domain and the frequency domain. The
fact that these two sets of data can be normalized together
effectively before proceeding with the evaluation of p
is significant since it indicates that the basic assumptions
associated with their analysis were correct and also points
out the fundamental equivalence in these two methods of
excitation of neutron waves.
105 / z = 14.78 cm.
z = 29.67
cm.*
105
a 1
G*.
103
z = 59.23 cm.
102
0 1200 2400 3600 4800 6000
Time (Microseconds)
Fig. 15. Thermal Pulse for the Pulse Propagation Experiment
for Several z Positions (from Reference 14).
In the discussion which follows in this section,
the data evaluation will be carried through the calculation
of the real and imaginary components of the complex inverse
relaxation length. The data acquired from wave experiments
will be considered separately from those of the pulse pro-
pagation experiment in obtaining the amplitude and phase as
a function of position and frequency. Then these data will
be normalized together for the evaluation of p.
The data presented here has been taken over a period
of two years in two separate wave experiments and in Hartley's
pulse propagation experiment. The curves of amplitude and
phase versus z position have been determined for 65 fre-
quencies ranging from zero to 1450 ops. For each of these
curves a minimum of 14 z positions have been taken; and in
many, where all three experiments have been normalized to-
gether, up to 40 z positions are considered. Moreover, to
obtain each of the 14 to 40 points on the 130 amplitude and
phase curves either a Fourier transform was performed over
the 256 channel output of the multi-channel analyzer or a
Fourier series expansion was made with between 50 and 120
channels per cycle.
The complex inverse relaxation lengths, obtained
from least-squares fits of the amplitude and phase curves,
were derived from a careful examination of from 2 to 8 fits
of each of these curves in which points with largest error
were systematically eliminated.- The presentation of all
of these data would be impractical. One can, however,
present some of the data which are indicative of the
results obtained and show sample calculations to indicate
the method of data analysis.
Determination of the Amplitude and Phase of the
Fundamental Thermal Wave of the Wave Data
For the initial analysis of the wave data use was
made of the University of Florida Non-Linear Least Squares
(UF-NLLS) code (28). Several new subroutines were written
for this code, and it was incorporated into another main
program to facilitate a Fourier series expansion up to
ten harmonics of the thermal neutron wave. The bulk of
the numerical calculations were performed in the WAVE-
ESTPAR subroutines which provided estimates of the para-
meters to the UF-NLLS code. The WAVE-ESTPAR subroutine will
now be discussed.
Information was first read into the code which
allowed it to calculate the period and number of channels
per cycle of the fundamental sine wave variation. The code
then read in the counts per channel recorded in the multi-
channel analyzer for the standard detector with no cadmium
placed between the source and the graphite pile (bare
standard data), for the standard detector with the cadmium
sheet in place so that only neutrons of epicadmium energy
entered the graphite (covered standard data), and for the
movable-detector counts per channel for the latter situation
(covered movable data). The data of interest, that arising
only from the thermal content of the source, are contained
in the movable detector counts per channel for the bare
run (bare movable data); but this information is contaminated
by counts arising from the epicadmium content of the source.
The fundamental amplitude, phase, and harmonic content up
to the tenth harmonic of the former three sets of data
were then calculated. From the two sets of standard data
the source variation from the covered to the uncovered
runs for each of the frequencies analyzed could be determined.
This information was then used to generate correction
factors to be applied to each of the frequency components
of the covered movable data. A set of data was generated
which represented the contribution of epicadmium neutrons
from the source to the bare movable data. This was
accomplished by individually multiplying the fundamental
and each of the harmonics of the covered movable data by
the ratio of the standards content of the same frequency for
the covered and uncovered runs and by using a subroutine
NEWCOV which synthesized the generated waves. This generated
epicadmium content of the bare movable data was then sub-
tracted from it to obtain the thermal neutron wave arising
only from the thermal neutron content of the source. Note
that this background subtraction has been made individually
for each frequency; therefore, changes in the harmonic
content of the source from one run to the next should not
affect the result. This elaborate method of background
subtraction was necessary because of the instability of
the neutron generator.
The fundamental component of the thermal wave and
as many of its first ten harmonics as desired were then
determinedsand the results were normalized to the bare
standard data. Corrections were also made for changes in
the constant component of the source as determined from
the total counts of the detectors for both runs. The data
were normalized so that the channel width of the multi-
channel analyzer, its total number of sweeps, counting
time, dead time between channels, and channels per cycle
could be varied without affecting the results. The print-
out of WAVE-ESTPAR included the normalized amplitude and
phase of the frequencies analyzed for the two detectors for
both runs so that any variation in the source was recorded.
The numerical methods used in WAVE-ESTPAR were rather
crude compared to those which could have been applied. The
correction for dead time between channels was made by cal-
culating the average of the counts accumulated in two con-
secutive channels, multiplying this by the dead time between
channels, and adding the result to the counts of the first
of these two channels.b The numerical integration associated
*The 212 plug-in-unit used for this experiment requires
10 microseconds to store the counts accumulated in a
channel before moving to the next channel. During this
storage time no input pulses are counted.
with the calculation of the Fourier coefficients were
approximated by a summation. The method can, however,
be justified for these data. The operation of the multi-
channel analyzer is such that it will store in each channel
all of the counts accumulated during the count-time of that
channel. In effect, the counts recorded in a channel are
the integral of the detector count rate over the time interval
allotted to the channel. If exact adjustment could be made
for the dead time between channels, the summation of the
counts accumulated in all the channels would be identical to
the integral of the count rate over the total sweep of the
multi-channel analyzer. The weighted summation performed by
the code approximates the exact integral over a period of
the fundamental sinusoidal variation by a weighted summation
of integrals over each channel. In each of these "channel
integrals" the sinusoidal function has been effectively moved
outside the integral, and the integrated count rate for the
time interval of the channel is multiplied by the value of
the sinusoidal function at the center of the time interval.
This approximation is mathematically expressed by
T N tj+r
I sin nwt y(t)dt sin nw(tj + t/2) I y(t)dt
oj=1 tj
(65)
where T is the period of the fundamental sinusoidal
variation of the count rate, y(t), r is the channel width,
ti is the time when the jth channel is opened; and it has
been assumed that the dead time was zero and that there were
exactly N channels per cycle. The accuracy of the above
expression is based primarily upon the value of N, which
varied between 50 and 120 for this work.
To prevent inaccuracies in the data, the harmonic
distortion, fast neutron content, and D.C. component of the
source were all held at as low a value as possible; and
every attempt was made to hold the neutron generator in
stable operation during all of the runs of a particular
frequency. These objectives were well realized, especially
for the second wave experiment; and for these data the
elaborate corrections seemed somewhat extreme. The original
purpose of WAVE-ESTPAR was to extract the thermal wave from
the bare movable data, keep track of variations in the source,
and provide estimates to the UF-NLLS code of the amplitude
and phase of the fundamental thermal wave. The UF-NLLS code
was then to be used to modify these initial estimates to
values which yielded the best fit to the experimental results.
Actually, WAVE-ESTPAR did such a good job of estimating the
parameters that changes in them by the UF-NLLS code were
made only in the third and fourth decimal places. In the
interest of saving computer time, the parameters calculated
by WAVE-ESTPAR were accepted for the majority of the data
considered.
To indicate the accuracy achieved by this technique,
Figures 16 and 17 show the amplitude and phase versus z
position of the 400 cycle data obtained from runs made
when the fundamental sinusoidal variation was 200 cps and
then 400 ops. For the 200 ops runs, the magnitude of the
400 ops distortion was about 3 per cent. As shown in the
Figures, the 400 ops results obtained from this small per-
centage of the thermal source agree quite well with that
for the 400 ops data. In fact, plotting the first harmonic
of a particular frequency became a very useful way to
determine if any mistakes had been made in the data input
since this small component of the measured signal was very
sensitive to inaccuracies or inconsistencies.
Analysis of the Pulse Data ~by Moore's Method"
Some enlightening derivations by Moore (16,17) provide
the guidelines for the analysis of data obtained from the
pulse propagation experiment. Moore expressed Fourier
moments, b,(q,v) ,of the measured data as
os -2wivt
In(+q,v) = I dt tn e rC(q,t)
..o (66)
* The author is indebted to Dr. M. N. Moore for several
stimulating discussions concerning the analysis of
pulse propagation data.
I I1I I 1I 1I
5 10 15 20 25 30 35 40
z ( cm. )
Fig. 16. Amplitude of the 400 ops Wave as Obtained from the
First Harmonic of the 200 cos Wave Data and the
Fundamental of the 400 ops Wave Data.
* Fundamental o
+ = First Harmonj
S= Normalization
,f 400 ops Exp
ic of 200 ops Exp
SPoint
2
10
8t
6
*=Fundamental of
+=First Harmonic
of 200 cps Exp
~jP $= Normalization Exp
I I
2.0
1. 5
1.0
0.5
z (cm.)
Fig. 17. Phase of the 400 cps Wave as Obtained from
the First Harmonic of the 200 cps Wave Data
and the Fundamental of the 400 ops Wave Data.
where #($,t) is the experimentally measured pulse at space
point q as a function of time, t, and v is the frequency in
cycls pe secnd. e then defined Fourier widths, W (q,V),
by the expression
-2nivt
Wn0~,v) = E I dt(t-t)n e
4( ,t)]/q( ,v)
(67)
where
t = l9 oqv)
(68)
and proved that, neglecting higher space and energy mode
contamination, these widths are related to the frequency
derivatives of the system dispersion law by the expression
Wn (<,V 2)n I Gn **R2,93,B: ) +
dn
+ d E-l~vBl 1(V,B )] J
n> 1
(69)
where
Gn 9n log A10 23 92'Q3,B )
(70)
and
~o(q,v) = Al0 ~)23 42,43,B )exp (E-al(V,B ) +igl(vB )]z 1
*
* |