Title: Propagation of neutron waves through heterogeneous multiplying and nonmultiplying media
CITATION PDF VIEWER THUMBNAILS PAGE IMAGE ZOOMABLE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00098220/00001
 Material Information
Title: Propagation of neutron waves through heterogeneous multiplying and nonmultiplying media
Physical Description: viii, 106 leaves : illus. ; 28 cm.
Language: English
Creator: Cain, Victor Ralph, 1934-
Publication Date: 1965
Copyright Date: 1965
 Subjects
Subject: Neutrons -- Scattering   ( lcsh )
Nuclear Engineering Sciences thesis Ph. D   ( lcsh )
Dissertations, Academic -- Nuclear Engineering Sciences -- UF   ( lcsh )
Genre: bibliography   ( marcgt )
non-fiction   ( marcgt )
 Notes
Thesis: Thesis--University of Florida, 1965.
Bibliography: Bibliography: leaves 104-105.
General Note: Manuscript copy.
 Record Information
Bibliographic ID: UF00098220
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 000565778
oclc - 13583933
notis - ACZ2197

Downloads

This item has the following downloads:

PDF ( 3 MBs ) ( PDF )


Full Text







PROPAGATION OF NEUTRON WAVES

THROUGH HETEROGENEOUS MULTIPLYING
AND NONMULTIPLYING MEDIA
















By

VICTOR RALPH CAIN


A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMIENUTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY










UNIVERSITY OF FLORIDA
August, 1965








































To Bliz














ACENOWL;EDGMENTS


The author wishes to express his appreciation to his supervisory

committee for their assistance. Special appreciation is due to the

committee chairman, Dr. Rafael Perez, for his unusual patience and

constant encouragement and guidance.

Acknowledgment must also be rendered to the Department of Nuclear

Engineering and to the Computing Center of the University of Florida for

their combined financial support of the computational work for this

dissertation.

Finally, the author wishes to acknowledge his debt to the Oak

Ridge National Laboratory and to the U.S. Atomic Energy Commission for

the assignment which made this work possible.












TABLE OF CONTENTS


PAGE

iii

v

vii

1




4




21




30

58

79



88











101




104


ACKNOWLEDGNENTS . . . . . . . . . . . .

LIST OF FIGURES . . . . . . . . . . . .

ABSTRACT .. . . .. . .. . . . . .. . ..

INTRODUCTION .. . . . . . . .. .. . . .. .

CHAPTER I. USE OF A ONE-GROUP ONE-DIMENSIGNAL GREEN'S

FUNCTION IN HETEROGENEOUS MEDIA .........

II. USE OF A ONE-GROUP TWO-DLKENSIONAL GREEN'S

FUNCTION IN HETEROGENEOUS MEDIA .........

III. USE OF AN AGE-DIFFUSION TWO-DLMENSIONAL GREEN'S

FUNCTION IN HETEROGENEOUS MEDIA .........

IV. CALCULATIONAL RESULTS ...............

V. SUMMARY AND CONCLUSIONS ..............


:ES

A. The Complex Arithmetic Version of GINV .......

B. Digital Computer Techniques for the Calculation

of the Error Function with a Complex

Argument .. . . . . .. . . . ..

C. Galanin's Thermal Constant for a Sinusoidally

Varying Flux .. . .. . . . . . ..

D. Physical Constants Used in Sample Calculations ...

:ES .. . .. . . . . . . .. -- - -


APPENDIC


















REFERENCE













LIST OF FIGURES


FIGURE


PAGE


1. Geometry for the One-Group, One-Dimension Green's Function
Calculation ............................................. 5

2. One-Group Flux Amplitude vs Distance Along the Z Axis for
One Cadmium Foil in Gra~phite, at 100 eps, with 1, 2, 4, and
10 Spatial Modes .................. ................... ....... 66

3. One-Group Flux Amplitude vs Distance Along thle Z Axis for
One Cadmium Foil in Graphite, at 1000 cps, with 1, 2, 4,
and 10 Spatial Modes ................... .................. ... 67

4. One-Group Flux Phase Lag vs_ Distance Along the Z Axis for
One Cadmium Foil in Graphite, at 100 eps, with 1, 2, 4, and
10 Spatial Modes .......................................... 68

5. One-Group Flux Phase Lag vs Distance Along the Z Axis
for One Cadmium Foil in Graphite, atl1000 cps, with 1, 2,

4, and 10 Spatial Modes .................. .................. 69

6. Flux Amplitude vs Distance Along the Z Axis for Two
Cadmian or Uranium Foils in Graphite, at 200, 500, and
800 eps, with One Spatial Mode ................... ...****** 71

7. Flux Amplitude vs Distance Along the Z Axis for Two
Cadmium or Uranium Foils in Graphite, at 200, 500, and
800 eps, with Six Spatial Modes ......................****** 72

8. Flux Phase Lag vs Distance Along the Z Axis for Two
Cadmiun or Uranium Foils in Graphite, at 200, 500,
and 800 eps, with One and Six Spatial Modes ................ 73

9. One-Group Flux Amplitude vs Source Frequency for Two
Cadmium Foils in Graphite, at z = 2, 8, 8.89, 15, and

30 em, with Six Spatial Modes .............................. 74









FIGURE PAGE

10. One-Group Flux Phase Lag vs Source Frequency for Two

Cadnium Foils in Gratphite, at z = 2, 8, 8.89, 15, and

3O cm, with Six Spatial Modes .............................. 75

11. Age-Diffusion Flux Amplitude vs Distance Along the Z Axis

for Two Uranium Foils in D2O, at 200 cps, with One and
Six Spatial Modes ................... ................... .... 76

12. Components of Age-Diffusion Flux Amplitude vs Distance

Along the Z Axis for Two Uranium Foils in D20, at
200 cps, with Six Spatial Modes ............................ 78

13. Regions of Complex Plane Requiring Different Techniques
for Calculating the Error Function ......................... 97












Abstract of Dissertation Presented to the Graduate Council
in Partial Fulfillment of the Requirements for the
Degree of Doctor of Philosophy


PROPAGATION OF NEUTRON WAVES THROUGH HETEROGENEOUS
MULTIPLYING AND NONNUTLTIPLYING MEDIA


By


Victor Ralph Cain


August 1965


Chairman: Dr. Rafael B. Perez
Maj or Department: Nuclear Engineering


Prediction of the behavior of neutron waves in heterogeneous media

is performed using one-group diffusion theory and age-diffusion theory. The

one-group theory is used in two different, finite geometries which allow

reduction to a one-dimensional problem and to a two-dimensional problem.

These cases result in diffusion kernels, or Green's functions, for the two

finite configurations.


The theory is then extended to include Fermi-age, or continuous

slowing-down, theory for higher energy neutrons. The approach is similar

to the Feinberg-Galanin heterogeneous reactor theory except that it is a~p-

plied to a finite geometry and includes time dependence. A finite diffu-

sion kernel is obtained which is similar to the results of the simpler

calculations. In addition, a finite-medium, Fermi-age kernel results which

describes the behavior of the slowing-down neutrons in the finite geometry.









Both the one-group and the age-diffusion developments are used to

calculate numerically sample configurations which are suitable for experi-

mental verification. These include rectangular assemblies of graphite and

heavy water which have poison or fuel rods inserted.


In order to demonstrate better the improvements of this work over

the Feinberg-Galanin theory, an extension to critical assemblies is made.

This offers the possibility of doing criticality calculations for physically

small assemblies which are not amenable to homogenization techniques.


Using the age-diffusion theory results, it is also shown that data

from neutron wave experiments may be processed in such a way as to give

experimental measurements of both the diffusion and the slowing-down

kernel.


v111













INTRODUCTION


In 1960 a program was initiated at the University of Florida in the

field of neutron wave propagation. The initial effort was in the determina-

tion of diffusion and thermalization parameters in moderating media. It was

shown by Perez and Uhrig (1)1 that these parameters were related to the

damping coefficient and phase shift per unit length of the neutron wave.

The work of Hartley (2) and Booth (3) confirmed these predictions, and

lately Perez and Booth (4) were able to perform accurate measurements of

thermalization parameters in graphite using the neutron wave technique.

The above mentioned studies were performed in homogeneous neutron moderat-

ing media. Some work has been done by Booth (5), in heterogeneous media,

and by Denning et al. (6) in two-region moderating media. Denning et al.

were able to show a definite reflected wave from the interface between the

two media which interferes with the incident wave.


The goal of the present work is to predict, as accurately as is fea-

sible, the propagation of neutron waves in heterogeneous multiplying and

nounultiplying media. The incentive to perform this work is twofold. From

the purely academic viewpoint, the propagation of waves through periodic

structures is of sufficient interest in the fields of quantum mechanics,

sound theory, solid-state physics, and others that extension of the





lUnderlined numbers in parentheses refer to corresponding numbers
in the list of references.









previously mentioned work in heterogeneous media becomes de s irab le. That

neutron wave propagation theory is applicable to these other fields may

perhaps seem reasonable by pointing out that the neutron wave propagation

we are concerned with here is a macroscopic phenomenon, qualitatively

similar to sound wave propagation. One simply replaces gas density with

neutron density. Quantitatively, the neutron wave propagation is the more

difficult problem since nuclear reactors are highly dispersive and absorp-

tive media in which disturbances of the neutron field are quickly damped.


Secondly, there is a strong incentive from the practical viewpoint

of gaining understanding about heterogeneous reactors. Traditionally,

homogenization techniques based on the Wigner-Seitz unit-cell theory have

been used to design heterogeneous reactors, and it is only recently that

theories accounting for the interactions between fuel plates are evolving

from the works of Galanin ( ), Feinberg (8), and others.

The study of neutron wave propagation through heterogeneous media

involves analytical techniques and assumptions similar to the developments

used in heterogeneous reactor theory. Hence studies of this type offer the

possibilities of experimental tests of the various theories in a way which

enlarges the realm of application of purely static experiments. Also, any

disturbance of the neutron field in a heterogeneous reactor can be analyzed

as a superposition of neutron waves, thus yielding information on possible

short-range instabilities for a particular reactor design.

Theoretical studies in this field are quite complicated because of

the interplay of the effects of the spatial heterogeneities produced by the

fuel plates and the energy distribution of the neutron population. Before









meaningful but costly experiments can be performed, the theoretical founda-

tions have to be developed, which is the goal set forth in this disserta-

tion.


At first sight, the obvious approach to the theoretical work is

simply to modify the Feinberg-Galanin theory to include time dependence so

that it can be applied to subcritical assemblies. In practice, this is not

satisfactory since most reactors are rather large assemblies. Because of

this the present heterogeneous reactor theory assumes that the moderating

medium in which the fuel is embedded is infinite spatially. Practical

experimental assemblies for neutron wave propagation experiments, however,

tend to be quite small in comparison. Consequently, it was necessary to

start from the basic equations including the finiteness of the feasible

experimental assemblies.














CHAPTER I


USE OF A ONE-GROUP ONE-DIMENSIONAL GREEN'S
FUNCTION IN HETEROGENEOUS MEDIAL


The :problem to be considered is the calculation of the transmission

of a thermal neutron wave, according to one-group diffusion theory, in a

simple heterogeneous geometry. The geometry chosen is such that it permits

a practical experimental setup and also results in an analytical form reduc-

ible virtually to a one-dimensional problem.


Figure 1 shows the geometry of the problem. The semi-infinite

parallelepiped of moderating material extends from x = -a to x = a, y = -b

to y = b, and z = 0 to z = oo. The transverse dimensions a and b are assumed

to include the diffusion theory extrapolation distance. The perturbation

in the medium is taken to be in the fonn of very thin sheets, the kth sheet

extending from x = -xo to x = xo and y = -b to y = b and of thickness h,

centered around z = zk. The system is driven by a sinusoidally modulated

source of thermal neutrons placed at the face z = 0.

The time-dependent, one-group diffusion equation describing this

system is


1 at(ryt) = U [Z + 2 (r)3 ] (r, t) (1.1)
v at a a





1The solution obtained in this chapter is an extension of work
originally performed by R. S. Booth.




















































Fig. 1. Geometry for the One-Group, One-Dimension Green's Function
Calculation.









where

v = the neutron speed,

Q(rjt) = the spatial and time-dependent neutron flux,

r = a position vector,

t = time,

D = the diffusion coefficient (units of length),

ye = the La~placian operator,

La = the macroscopic cross section of the moderating media,

6Ca = the space-dependent macroscopic cross section of the

absorbing sheets.

The space-dependent perturbation may be written as




Ba(r) = BlaR( u x xo -u(x x )
k=1


u z -z k -3 ,(1.2)



where Blak the macroscopic cross section of the kth foil, u(y) = the

unit step function, and M = the total number of foils. If h is taken to be

sufficiently small, the first two terms of the Taylor's series expansion

for the unit step function in z may be used in Eq. (1.2), reducing it to




Bla(r) = h lu(x + xg) (x x)] a~,k 8(z z) (1.s)


where 6(y) is the Dirac delta function.








It should be noted that if the sheets contain fissionable material,

8a,k is replaced by 82a,k Ff~k, where v is the neutron multiplicity
and 82f~k is the macroscopic fission cross section of the _kth foil.

The sinusoidally modulated neutron source is assumed to be isotropic

with a Maxwellian energy distribution and can be represented by


S(r~t) = So + Re [Slr)(= e iwt (1.4)


where i =9 J-, w is the angular frequency of the modulation and So is the
steady component of the source present at z = 0. If the neutron flux is

assumed to be separable in space and time, the oscillating part of the flux
can be represented as (1)


Q(r t) = Re[r (1.5)


Equation (1.5) may be introduced into Eq. (1.1), resulting in

[ + 6Cz (r) 3 Q, v (r) (1.6)


or, in slightly different notation,


D Q- = iw L(x,z) e(r) (1.7)

where

a a vCa (1.8)

Do vD (1.9)


GL(x,z) a 6Ca(r)/D .


(1.10)





8


The boundary conditions will be taken to be


m(r) )
x=+a


= i (r) = 0
z-CD


(1.11)


In addition, at the source plane,


iwt
e


(1.12)


ae(r t) 1
D = -Re [S(r) I
dD z 1 2
z=o z=o


Using Eq. (1.5), this reduces to


ae(r) 1
- D =-S(r) .
az 2
z=o z=o


(1.13)


The spatial flux, Q(r), may be represented by an expansion in the

complete set of transverse eigenfunctions defined by


(9 + B a) 4 m(x1y) = 0


and by the boundary conditions at x = +a and y = +b plus the conditions of

symmetry around x = 0 and y = 0. Normalized, these eigenfunctions and their

related eigenvalues are


1 (2R 1)nx (2m 1)xJy
4 m(x'y) = --cos cos
Am~ 2a 2b


(1.15)


and


(1.16)


wJhere R and m may take on all integral values.


= (r)
y=+b


(1.14)


(20 i--+ 1 (2m 1)Jrn









The expansion of Q(r) is given by


Q(r) = Qm(z) bQm(x,y) .
,m


(1.17)


Substituting Eq. (1.17) into Eq. (1.7) gives


al + iw
Do i im,(x,y) Q~m(z) = BL(x,z) Q~m(x'y) ~Qm(z) .(.8
a~m a~m


dy (xy) esutsin
-b


a


-a





a,m


4 m(z) ,


(1.19)


where


a, + iw
Do


p B
pq9 1pq


(1.20)


and



BLpqim(z)' S ax dy, *F (x~y) 6L(x~z) Qm(x,y)
-a -b


(1.21)


Before continuing further with Eq. (1.19), note that the two boundary

conditions in z have yet to be applied. In order to make use of the source

boundary condition, S(r) must be expanded in the same transverse eigenfunc-

tions. That is,











Am


S(r)
z=o


(1.22)


whe re


a

Sem; =ax
-a


"m(x'y).


It should be pointed out that this expansion, as well as the flux expansion,

assumes symmetry around the x and y axes. If one uses a source containing

asymmnetric components, oQm(x~y) must also include sine tents. Inserting the

expansions, Eqs. (1.17) and (1.22), into the source boundary condition,

Eq. (1.13), gives

de (z) S
am am
- D (1.24)
dz 2
z=o


The boundary condition at z = OD becomes


lim
z- me(z) = O (1.25)


reduction of Eq. (1.19), substitute Eqs. (1.3)

The y integration may be performed immediately,

delta. Also performing the x integration gives


Now returning to the

anid (1.10) into Eq. (1.21).

yielding 8qm, the Kronecker
the result


M
6L (z) =h 8 8Fp 82~~ 8(z z )
pqim D qm p ~
k=1


where 6Fp is given by


b

Sdy S(r) Z=
-b








Now define the Green's function G q(z,S), satisfying


ac; 2~ pqpq(ZS D_( z


(1.32)


and the homogeneous forms of the boundary conditions satisfied by 4 q(z).

Multiply Eq. (1.50) by G q(z,S) and Eq. (1.92) by Op (5). Subtract the
co
two results and operate with d6 to obtain


d i
o


(G(z,S) ddg(5) 8(5 G(2cc ) d"i F(i) G(zi) = ez)


(1.55)


The subscripts have been dropped temporarily for economy of notation.

Integrating the first integral by parts results in


m(z) = D (z,S) d5 g) ) jGzg


oo 00
-, D dg az, d( + D S idi(~cz~


oo
- D d i F(i) G(z,i)
O


(1.34)


The first two

gradient of G


integrals cancel, both + and G go to zero as ( - c, and the
is zero at 5 = 0. Therefore, the desired integral equation


- D dj F(i) G(z,I)
O


(1.55)


Q(z) = D G(z,g)








.(a + p 1)rxo
sin --~---
_iT--3~S~a


(a p)nxo
sin ------
a
pl= (a p)n t


(1.27)


(2 1)nx,

a+ (2a 1)Jr


In order to simplify Eq. (1.27) further, let us assume that xo is suffi-

ciently less than a that the first term of the Maclaurin's series may be

used for the sine terms. Equation (1.27) then becomes


2x
Fpa ao


(1.28)


Sfor all a and p.


Noew inserting Eq. (1.26) and Eq. (1.28) into Eq. (1.19) gives


BI 6$r(z z) 1L (z).
k=1


It may be noted from the above equation that the one-dimensional problem

turns out to be not really one dimensional, in the sense that coupling re-

mains between the spatial harmonics associated with the x direction.

In order to solve Eq. (1.29), it is convenient to convert it into an

integral equation. Substitute 5 for z and rewrite as


(130)


where
2hx
Fpq() Da akk
L k=1


d~ o


(~2 i~~~ bpa(i)= %$i)










The problem now is to determine the Green's function defined by

Eq. (1.32). The solution to the homogeneous form of Eq. (1.32) is


G(S,z) = Ae-pz + Cepz


(1.36)


where 5 and z have been switched. The particular integral of Eq. (1.32) may

be found by the technique of the variation of parameters. Let the constants

in Eq. (1.56) depend.on z (they also depend on 5, of course, but this may

be ignored at this point):


G(S,z) = A(z) e-pz + C(z) epz .


Differentiation of Eq. (1.37) gives


(1.57)


aG(S,z) -pz pz dA -pz dC pz
= p Ae + p Ce + e + -e


(1.48)


Choose A and C such that

dA -pz dC pz
--e + --e = O 3
dz dz


Then differentiate again:


a~(Sz) = pZ Ae~P + Cp e z
az,


(1.59)


dA -pz dC pz
-p e + p --e
dz dz


(1.40)


Inserting Eqs. (1.40) and (1.57) into the original different

for the Green's function results in


A ~p e-pz ge-pz + C p epz epz


A-pz dC pzp 6(z g)
dzdz D


ial equation


(1.41)











-pz dA pz dC 6(z 5)
-pe + pe
dz dz D


(1.42)


Equations (1.59) and (1.42) are two equations in dA/dz and dC/dz which may

be solved simultaneously to give


dA 8(z -g) epz
dz 2pD


and

dC 8(z 5) e-pz
dz 2pD


Integrating Eqs. (1.4)) and (1.44) gives


A = e S/2pD + B


and


C = e P/2pD + E


(1.4~)


where B and E are constants to be determined from the boundary conditions.

Using Eqs. (1.45) and (1.46) in Eq. (1.56) gives

G(S~) -p pt p(g-z) .-p(S-z)
G~~)=Be-p + Eep e (1.47)
2pD 2pD


or

G(5,z)-pz pz sinh p(S z)
G~~)=Be + Ee + (1.48)
pD


The one-dimensional Green's function has a discontinuity in slope at the









source point, z = E. Actually, this does not need to be known beforehand,

since straightforward application of the only two boundary conditions now

available leads to the solution G 0 Then, in order to obtain a solution,

before applying the boundary conditions, Eq. (1.48) must be divided into

two parts, one which will apply for z < ( and the other for z > E. Now,

two more boundary conditions are needed, which may be found by integrating

Eq. (1.32) over a small region including the source point:



dzP p (5,z) =- 6(zD .5 (1.49)



Assuming that the second derivative is more singular than the function it-

self, one obtains, in the limit of small E,


aG(g,z) aG(5,z) 1 (.0
az az D



Therefore the two boundary conditions are continuity of G at z = ( and

Eq. (1.50), which gives the magnitude of the discontinuity in the first

derivative.

The two parts of the solution will be denoted as


L L -pz L epz sinh p(E z) (1
G (g~z)= B e E e + pD ,z<((.1


and


GR~eZ =R ,-pz R epz sinh p(5 z) (.2
G (g~z)= B e E e + pD ,z>(.(.2


The z = 0 boundary condition is









dGL(S,z)
dz


=-pL + pEL cosh pi O
D


(1.55)


z=o


At large z, the boundary condition is


lim cR~g,z) =o lim E epz e-p ep ]=
z oz o 2pD


Equation (1.54) immediately determines ER:



2pD


The continuity condition, at z = 5, leads to


BLe-P + EL pS = BRe-P + ERe .


Equation (1.50) becomes


(1.54)


(1.55)






(1.56)


- pELeP + 1


= -


- pBRe-P + pEReP + pBLe P


(1.57)


Using Eq. (1.55)1 Eqs. (1.53), (1.56), and (1.57) may


be put in matrix form


as


,sh pS
D

2 D

2D


co


-p

e" 5

-pg
-pe


(1.58)


The three solutions to Eq. (1.58) are


p O B

eP e'P EL

pS -p5 R
pe pe B










BL sinh pS (1.59)
pD


EL epD (1.60)
-pD

BR -- (1.61)
2pD


Upon substituting Eqs. (1.55), (1-59), (1.60), and (1.61) into Eqs. (1.51)

and (1.52) and rearranging, one obtains


3L~gz = e Dosh pz (1.62)


and

-pz
cR~g,z) = e Dosh pS z > 5 (1.63)


It may be noticed that Eqs. (1.62) and (1.63) define a function

which is symmetric with respect to interchange of the coordinates of the

source and observer locations (z and 5), thus satisfying one of the nec-

essary criteria for a Green's function (9). Equation (1.63) may now be

substituted into Eq. (1.55), giving

.-pz d -pz
Q(z) = D -O~ D F(S) dS
pD dS pD
g=0 0


D e cosh pz Fi) dj (1.6Li)


Reinserting subscripts and using Eq. (1.24), one obtains







-pz
Pq -p z z~ 5
S e pq
Opq(Z s 2p Dp pq pQ
pq pq



at e (5) (1.65)
pq


Since F q(h) contains 8 functions in j, Eq. (1.j1) may be inserted into

Eq. (1.65) and the integrals in 5 performed, leading to


S e Pq 2hx -
+ z)=ps o eB ohp z4(
mpq( 2p D p Da~ 9 a 6~ k Cs pq k* eq k
pq pq k=1



+ csh Pp, z i, B6a,k e-p$ zk p + zk) i (1.66)
k=Mr+1


where M' is determined by zM' < z < zM'+1. If z equals one of the zk's,

it can be taken to be either zM' or zM'+1l since the Green's function is

continuous at the source point.

In order to obtain a solution to this system, the original expansion

of Q(r), Eq. (1.17), must be truncated. Assume that the desired a~pproxi-

maztion is



*(r) ~V ZL (z 4I "(xy) (1.67)
=1 9=1


For each value of q, a set of LM simultaneous equations may be obtained

from Eq. (1.66) by writing it for each of the L modes and for z = z ,

z ,...zM.This set determines each of the LM values of 4 z








Writing these values as a vector, where the superscript T indicates the

transpose,


(1.68)




(1.69)


where [I] is the identity matrix, the vector B_ contains the constant

(source) terms, and the matrix [A ] contains the Green's functions. B
Q q
may be written as


Bg [il ~ lqZI l lq 2 Y


(1.70)


whe re


S
aq b2 p D


The [A ] matrix may be partioned into L' submatrices which are each M by
M. The n~p element of the I,J submatrix is given by


(A ~)n~ = Bp (z ) cosh P~lp zeIQn


(1.72)


Sn > p


(A ~)n~ = Bp (z ) cosh pgzn e lq p


(1.73)


, n



where


Og = $1(z,) m1q(z2) '. lg(zM) 2qg(z1' .. mLq(sM)


the set of equations may be put into matrix form. That is,


iAq q + [I] m = B ,


i~- Cq Pgz#


(1.71)







2hx 6C
B, (z )' ,LEoa (1.74)


Solving Eq. (1.69) Q times by inverting the matrix [A ] + [I] gives
the 0 s:


g [A ] + [Il]l B (1.75)


which may be inserted into the summations in Eq. (1.66) The results of

Eq. (1.66) are then inserted into Eq. (1.67) to give the desired approxima-
tion to the spatial flux *(r) .













CHAPTER II


USE OF A ONE-GROUP TWO-DIENESIONAL GREEN'S
FUNCTION IN HETEROGENEOUS NEDIA


The problem of the previous chapter will now be extended to two

dimensions by allowing the perturbations to be placed anywhere in the x-z

plane contained in the moderating material. Equation (1.1) will still apply,

but now BZa(r) will be taken to be a general function of x and z. Equation

(1.7), with its boundary conditions given in Eqs. (1.11) and (1.13), will

be taken as the starting point, leaving 8L(x,z) general.

In this case the spatial flux Q(r) will be expanded in the single

set of transverse eigenfunctions, defined by


--- + B 4 ~(y) = 0 (2.1)
dy2


and the boundary conditions 4 (+b) = 0 and 4m(Y m (-y). The eigenfunc-
tions are






and the eigenvalues are



[~m (2m 1)T(b(2n


where m = 1,2,3,... Using the following expansions of the flux and the

source,













and


S(r) = Sm(x) em(Y) ,



in Eqs. (1.7) and (1.13) gives


(2.4)


(2.5)


o+ iw
C~ o
m


4 (x, z) 4 (y) = 8L(x, z)


m(x, z) (y)


(2.6)


-. D 4m(x, z)
m


1 Z

m


Sm(x) em(y)


(2.7)


b
with day mn(y) gives
-b


Operating


O-


+~ -
fz n


4 (x, z) = GL(x, z) + (x, z)


(2.8)


whe re

2 ao + iW
n n Do


(2.9


and the corresponding boundary condition


1
2 Sn(


(2.10)


-D a
z= o


m,(y) I
z= o








The remaining boundary conditions given in Eq. (1.11) still apply.

It is evident that this system may be treated similarly to the prev-

ious problem if a two-dimensional Green's function can be found which is

the solution of


+2 aZ 2 n1(x z;xo ,z ) =-8 x -x) 6(z z ) (2.11)



with boundary conditions



Gnt (+a3c,z;x ,z )= n _li~m G =O .(2.12)
z=o


This Green's function may be found by expanding in the eigenfunctions of


d + k ((x) = 03 #k(+a)= O .(2.13)



These are


(2 1)Jrx anx (.4
Jr(x) = AR cos + BQ sin --(.4


Note that k assumes two different values for the two different eigenfunc-

tions, the sine and the cosine. The expansion of Gn(x,z;x ,zo) is


[ (2& 1)Jx axx
On(Xjx,zx,z Z) = C (z'zo ) [A (x cos 2a + B (xo) sin -a
a (2.15)


Inserting Eq. (2.15) into Eq. (2.11) gives, with arguments dropped for con-

venience,








2 2
cos B -sin n
I ~2a 2a a


a~e o (2 1)lrx nx
+Z A2 os + BQ sin -- -~ p C
a z a


['i O (2 1)J[x inrx] x )hZ 1 (.
cs 2a + BQ sin --=-6x-x)8z-z).(.6



-a



asi 0^E a gives
-a


C.A. 2j- 1)n a +---1J Aa-p C.Aja =- 8(z -z )cos (2 )x
2a8a a 2a
(2.17)


and


- C.B. a + J- Bja -p C.BJ a =- 6(z z ) sin --- (2.18)
JJ 8a n a


These two expressions may be rewritten so as to have the left-hand sides

functions of z and zo only and the right-hand sides functions of xo only.

They become


a1_ rL 2 f (2j 1)x CO cos -no
za2 Ln 2a 2a
= (2.19)
8(z zo) a A.


and










a n j sin ---
bz a .(2.20)
6(z zo) a B


Obviously the C.'s in the above two equations are not the same; so the C.
will be replaced by A! in Eq. (2.19) and by B! in Eq. (2.20). Notice also

that both sides of Eqs. (2.19) and (2.20) are equal to constants and no

loss in generality occurs if this constant is taken to be unity in each

equation. With this step, A. and B. may be obtained immediately:

1 (2ij 1)nxo
A.(x ) =-cos (2.21)
o a 2a


and

B.x) 1 sin -o (2.22)
B. a) a


and with the definitions


42 -p (2j 1)n (2.23)


4nd


-a p + (2.24)


the two equations for Aj and B:J are


-~2 2 A ;(z,zo) = 6(z z ) (2.25)


and









(2.26)


10m
SA! =
z-Koo 3


lim
BI = 0


(2.27)


z=o


z=o


plus the function continuities and first-derivative discontinuities at

z = zo. These functions are virtually identical to the one-dimensional

Green's function developed in Chapter I and can therefore be immediately

written down, being


e cosh uz
A (z,zo) =



e cosh Ciz




and


e zoshg


Sz < zO




S z > z


(2.28)


>z < zO


e-~ cosh gzo


(2.29)


Sz > z


Collecting Eqs. (2.21), (2.22), (2.28), and (2.29) into Eq. (2.15), the

complete Green's function is


a2 B!(i.zoz) = B(z zo) .


The related boundary conditions are


zAI z









-pro (26 1)nxo 2 ~t
G x z~oz Zre cosh CLzc o (s -1x
Qnxzxoz L IJa os 2a os 2a


-rlzo asrx
e cosh rlz o .rr xx2.0
+ ~~~sin ---- sin -- z IIa a ao



and


-C ~ o (P l o (2 1)rxx
G (~z~ z =cos cos
n o o a 2a 2a



e-~ cosh gzo axxo alxl
+ --------sin --- sin --- I z > z (2.31)
rla a aJ o


Similarly to the procedure followed previously, this Green's function

may be used to formulate an integral equation for en(x,z). xo and zo are

interchanged with x and z in Eqs. (2.8) and (2.11). The difference between

Gn(xo,z 3x,z) times Eq. (2.8) and @n(xo'zo) times (2.11) is taken. The
oo a
resut i opeate on ith dz d ogvwt ruet n ne
o -a

gration limits dropped,



dz~ ~ Sx~ Gn z a o dzoJ dx e _



o o n a2z, o o n n





In the terms containing partial with respect to x, an integration by parts

in x may be done, and similarly with the z partial. All remaining surface








integrals on the left-hand side of Eq. (2.32) cancel and most of the line

integrals are eliminated by the boundary conditions. The remaining terms

give


a as(xo,,zo)
@n(x'z) = dx, Oncx ,ZO;x,z) n zo
-a zo=o zo=o

Co a

dzg axo Gn(xo.z 3x.z) BL(x ,z ) "n(x ,z ) (2.53)
o -a


Equations (2.10), (2.30), and (2.31) are to be substituted into Eq. (2.53),

along with the desired form of BL(xo~zo), in order to obtain the integral

equation for On(x)z).

If the fons of BL(x,z) of Chapter I is used in Eq. (2.53), it is

straightforward although rather tedious to show that Eq. (2.33) reduces to

Eq. (1.66). To do this, the following expansions must be used:


d (xIz) = e mj(z) O.(x) (2.34)



and


S (x) =Z Sn Q.(x) (2.35)



where


mjx) 1 C (2j 1)rxx
+ x o 2a (2.36)





29


These expansions are justified, since if 6L is symmetric in x [and it is

assumed as before that S(r) is symmetric around the x and y axes ], then

en_ will also be symmetric in x. It is then a matter of performing the xo


1 ( 2p 1) stx
and zo integration in Eq. (2.jj) and operating with Ja cos 2a d
-a


to obtain Eq. (1.66) .













CHAPTER III


USE OF AN AGE-DIFFUSIGN TWO-DIMENSIGNAL GREEN'S
FUNCTION IN HETEROGENEOUS MEDIA


The development of this chapter is similar to that of the preceding

chapter except that the theory is to be the continuous slowing-down model

with a thermal group. The geometry is the same semi-infinite moderating

medium, with the perturbations now being small fuel rods which, as before,

run the full length of the y axis. This model is similar to the Feinberg-

Galanin treatment of heterogeneous reactors (7) and (8), with two major

exceptions. In this study time dependence is included and the finiteness

of the medium is taken into account. This will result in a finite medium

Fermi-age kernel, or Green's function, as well as the corresponding diffu-

sion kernel obtained in the previous developments of this work.

For lethargies below thermal, the applicable equation is


-D(u)~ o2 (ru,t) + Ca(u) e(r,u,t)


1 ae(r,u, t) aq(r u, t)
=T; + S(r~ut) -u '(3.1)


where the nomenclature is similar to that previously used except that now

all quantities are lethargy-dependent and S and q represent the neutron

source and slowing-down density, respectively. It has been assumed that

Za(u) for u < us, where us is the thermal neutron lethargy, is independent
of position. This is valid as long as the fuel rods are small.







The relation between the flux and slowing-down density applicable

to an infinite medium will be assumed, that is, that


9(r,u,t) = 5 It(u) Q(r,u,t) (3.2)


where g = the average logarithmic energy decrement, and Ct(u) = the total

macroscopic cross section in the moderating media.

The thermal group is described as before, except for the addition

of the source of thermal neutrons from the slowing-down density. This

relation is


[-D +E + -1 a ra)= t~ ,t) -A(r, t) ,(3.5)



where A(r,t) represents the absorption in the fuel rods and the subscript

s indicates that the quantity is evaluated at the thermal neutron energy.

In order to define the quantities S(r,u,t) and A(r,t), it will be

assumed that (1) all fission neutrons are born at lethargy u = O, (2) all

fuel rods are lines, the kth rod being representable as 8(x xk) 8(z zk)

and (3) the driving source is located at z = O, producing only thermal neu-

trons and is of the form Re[S(x,y) 6(z) e in. As before, since the source

is present only at the boundary z = O, it will be included as a boundary

condition rather than being included in Eq. (3.5).

Since this treatment is an extension of the Feinberg-Galanin treat-

ment, some of Feinberg's notation will be adopted. Define the thermal

constant 7k as the ratio of the total net current of thermal neutrons into

the kth rod to the value of the thermal neutron flux at the surface of the

kfth rod. 4 (r,,'t) will be used for the thermal neutron flux at the








surface of the kth rod. The neutron yield per thermal neutron absorption,

r7, will be assumed to be the same for all rods. With the above assumptions
and definitions, S(r,u,t) and A(ryt) become


s~~~r,~ ~ ut ~)(0)
s~~~r,~ u, t)= e(u k ms(rk t) 6( r -) rk.



and


A(r t) = ksrt) 6(r -rk) .(5)
IJ ~k


The ratio of the scattering to total cross sections at the source lethargy

in Eq. (3.4) represents the probability of a neutron starting the slowing-

down process.

Due to the separability of space and time for the flux, slowing-

down density, and source term, Eqs. ( .1), (3.2), and (3.4) may be combined

to give


[- D(u) + Zalu) + + i LT(u) q (r4u)


2 (0) Z Bu k ( ( 36
= C t(u) -"ppy n ~ )7 s rk r-r



Similarly, Eqs. (3.5) and (3.5) become

[- D + 2 E e(r) = (r~u U) 7k 6 (r ~l) .(>.7)








The boundary conditions on es(r) will be the same as those used for the

thermal neutron flux previously; the boundary conditions on q(r,u) will

be identical, for u > 0, except that at z = 0, q(r,u) will be taken as

zero on all boundaries. This corresponds to assuming that Q(r~u) = 0 at

all boundaries for u > O, which is consistent with the thermal neutron

flux assumptions.

Defining

iw




Eq. (3.6) may be written as



[- 5t(u) o' + B(u,w) + au ]4(rP)U






M~ulti-plying through Eq. (3.9) by the integrating factor



exp B(u',u) du ], Eq. (3.9) may be written as











-DZ u) ~ q(r,u) eo


rIB(u ',w)du'
o


I (0)
cosO/ b u) 7k ~s ~k) r gIk)


(3.10)


In order to simplify the notation, define the Fermi age, 1, from


dz D(u) with z(u = 0) = 0
du 5( I~u


(3.11)


and also define a frequency-dependent resonance escape probability,

p(u,W), as




p(u,w) =-eo


(3.12)


Notice that the above definition is similar to the usual resonance escape

probability except that an additional 1/v absorber has been added by the
oscillating source.

With the two definitions [Eqs. (3.11) and (3.12)] Eq. (3.10) be-

comes


I C(0) (r 6 )
1 s~,- CtOJ k -rk)


(3.15)


Bo 'm~u
+ fru)e








The partial with respect to the Fermi age will be eliminated by

using the Laplace transform. Defining



6(rysa e dr (3.14)



Eq. (3.13) becomes



ye (rs) s 8(rys) +FijT 0 k1 s k r
k (3.15)


The term q(r,0) is zero since q(r,0) was incorporated into the original

equation and is, of course, the negative of the right-hand side of Eq.

(3.15). If it had been chosen that this source be introduced as a boundary

condition, it could have been omitted from the original equation and have

been introduced naturally at this point.

The obvious step at this point, consistent with previous methods,

is to expand 9 in the y eigenfunctions given in Eq. (2.2), which were

determined by Eq. (2.1) and the boundary conditions of symmetry and

k(+b) = 0. The required expansions are



e(r,s) = e,(x,z, s) 9,R(Y) (3.16)



and








For future use, it should be noted that Eq. (3.16) implies that






Eq. (3.16) may, in fact, be obtained from Eq. (3.18) by Laplace-transform-

ing both sides of Eq. (3.18), where QQ is the inverse transform of 9 .

Inserting Eqs. (3.16) and (3.17) into Eq. (3.15) and operating with


(y)I~ dy gives


Ox


+ -
a2z


a (x'z~S)


Z (0)
=- 4 )


7k 6(x xk) 6(z zk


Q(x~z)


(3.19)


whe re



[2(2a 1)n]
S2b


(3.20)


Now 6 (x,z's) is expanded in the x eigenfunctions of Eq. (2.14),
defined in Eq. (2.13):


6 (x'z's); =B 6(z,s) 1- cos
J


6!(z,s) 1- sin j x
Ja a


(2j 1)nx +
2a


(3.21)


Inserting Eq. (3.21) into Eq. (3.19) results in


S(s + B )








in
a


6 j
j


C (2j 1) 1Xax + B
2a j


1


sin j[
a


a8z


a8z


+


-cs(2j 1) ax
via 2a


1- sin jC
,A a


1 (2j 1)nx
- 0 ) .-- COS -L--- +
Ja 2a
J


7Ex
a


8' -- sin


x (0)
=- 0) r


k, 6(x xk) 8(z zk) Q(x'z)


(3.22)


a
S1


(2m 1)71x
cos 2adx gives


Operating on Eq. (3.22) with


a82 z ,s)




7k 8(z zk)


- (s+ B ) 6( z~s)


(2j 1)nxk
2a


1
cos C


Q(xk z)


(3.25)


a
1 .mnx
-- sin -- dx gives
-a


Similarly, operating on Eq. (3.22) with


a8 9(z,s)
+ (s + B ) 9!(z~s)
a2 J


a [


- 6!(z,s)
3


Is(0)
=- z (0


1 anxk
--sin ---- (kz)


Z k B(z k)


(3.24)


(2j- 1 a
2a


- 8(z,s) (2j 2al)n r2a








Defining



2 as + B2 + (2j 2al)n (3.25)


and


52 's + B~ + ,a (4.26)


Eqs. (g.2g) and (g.24) may be rewritten as


a26.(z,s)


J 2 B(0~)1 (2 -1)x


= -I Ok 6(z -zk) cos 2a (xz) (3.27)



and

a8 9(z s)
Jr- g 9(z,s)




Zs(0) Z kh 8(z zk) J1- sin k +(xk'z) .(3.28)
Ska


The variable s has now been hidden in CI and g; so the above two equa-

tions are effectively one dimensional. It may also be noted that Eqs. (3.27)

and (g.28) are Green's function equations. Consider the boundary conditions

which are appropriate for these equations. Aside from the usual conditions

of continuity at z = zk and the discontinuity in the first derivative at








that point, the appropriate conditions must be derived from the original

z acis condition on q(r,T), which are




z=o


Using these relations along with the expansions that have been utilized,

Eqs. (3.16), (3.18), and (3.21), results in the following boundary condi-

tions for Eqs. (3.27) and (3.28):


e.(z,s) = lim 9(~)=0(.0
z=o

and


9!(z,s) = lim 9(~)=0.(.1
z=o


Except for the fact that the function itself, rather than the first

derivative, is zero at z = 0, we now have a problem identical to the one-

dimensional Green's functions of Chapters I and II. By following procedures

identical to those used in Chapter I, it is easily shown that the solution

to


d29(z,zo)
K (z,zto) = (z -zo) (s.52)
dz2


with boundary conditions



6(z,zo) = im 6(z,zo) = 0 (3.SS)
z=o







-K2
B~zzo) e sinh K zo


-Kzo
e sinh K z


Sz > zo


(g.~4)


Sz < zo


Comparing Eqs. (2.25) and (2.28) with Eqs. (.3.2) and (3.54) shows that the

change in boundary condition has simply resulted in the cosh being replaced

by the sinh. By writing the sinh as the difference of two exponentials in

Eq. (3.54), it is easily seen that the two solutions may be combined into

a form which will prove to be useful later:

-K z-zo( -K(z+zo)
e(z,zo) =e 2K 2K (3.55)


Using this result, the solutions to Eqs. (g.27) and (g.28) may now

be immediately written down, being


Zs(0) 11 -p z-zk
6 (z, s) = e~ ~


1(2j 1)rxk
7k~ -cos 2a e(xk zk)


(~.36)


CI( Es(0)t -g z-z ['"'
6 (z, s) = e


1 Exk
7k sin ------ + ( zk)


(~~7)


-p";z"zk


-S(z+zkl





41


Inserting Eqs. (3.56) and (3.37) into Eq. (3.21) gives


-1J(z+zki
-e


c (0)
-sif


j k


e (x'z's) =


(2j -1)nxx
cos 2a


c (2j 1)~xx
cos 2a


R(xk'zk)


-g(z+z )


ax


(3.38)


sin j~Xsin


The next step is to inverse-Laplace -transform Eq. (3.38), obtaining

the expansion coefficient in Eq. (3.18), which is


Q (x,z,,w = co


6 (x'z's) .


(3.59)


Rewriting the definitions of y and 5 as


I +B


and 52 = s + B' ,


(3.40)


whe re


and B 9 =


(3.41)


it is evident that each term to be inverse-transformed is of the form



e~- t (3.42)


-lJ z-zk
Yk e


- Iz-z l


m (xk, zk)


B = (2c~ 2bl7n (2j 2al)x 2


22bl ~ a








which has the inverse


--- ez24r for zu > O .



Using Ea_. (3.41), the inverse transform of Eq. (3.38), Q (x,z,-r), is


(}.43)


Zs(0)
= t


-(z+zk)24


2a /r Yk k


(xzw)


cos cos + e
2a 2a


"F


sin j sin
a


4 (xk)zk)


(~.44)


Returning now to the thermal-neutron diffusion equation, insert Eqs.


(3.17) and (73.18) [with Ts replacing -r in Eq. (3.18)] into Eq. (j3.7)


and


b


-b


(jy) dy to obtain


- Ds


a2
a22


Ds B as v


= p(. s (,z,,Ts- Ik
k


,(x,z) 6(x xk) 6(z zk)


.(3.45)


Using the expansions


q~~)= z cs(j2alx + z i a
I,-~lr ~JQZ o 2 n j,) i x


(g.46)


-(z-zk) /AT


-BE.T
e 3


3m(x,z)








and


j~(zTs) cos (2j 2al)rrx + Q(z,ts sin j~


1
Z


(3.47)


Q (x, z,s 1


a

cos (2a l~ dx gives


in Eq. (3.45) and operating with


-D a+
Sdz


+ -
as v
s


(2j 1) xx
Yk cos 2a k


4 (xk~z) S(z -


zk).58


= P( s' j zTs


a
with sin dx gives
-aa


Similarly expanding and operating


a 4' (z)
-D j
Sdz


[D B~ + Ea


]m' (z)


+ -
v
s


Yk sin a


k


(3.49)


Q(xk~z) 6(z zk)


=p(Ts ,W


Defining

C +-
K2 B + s
3 D


(3-50)


Ds 2aj(2j 1)Jl ]2 I(z)


i j(z)


[D B~ + E


Dsa "j








Z + -
K' ,B' + s
j D


Eqs. (3.48) and (3.49) may be written as


dze


-p(r ,w)
- K 4.(z) = s Q(,
0Ds &jts)


Z (2j 1)nxk
7k cos 2a a(xk'z) 8(z zk)


1
+ ----
D/


(3.52)


and


dze


- K'2 4! (z) = Q(zT )
03 D~s s


Yk sin 4 ,(x ,z) B(z z )


+D


Comparison of the expansion, Eq. (3.47), with Eq.


(3.44) shows that


(z-zk)

Zlk~~km l


(z+zk)2
4r
- e s


s


Z (0)
Q (zTs) c o


(2lj 1)nxk
cos 2a a(xkrzk)


(5.54)


and









(z-zk)
s


(z+zk)2
s
-e


Z (0) js
Q (z, Ts 2 k k


sin -a~ (xkfzk)* (s55)


Since Eq. (3.52) along with Eq. (3.Sk) is so similar to Eqs. (3.SS) and

(5-55), attention will be restricted to the former set; then the latter

set will be handled by comparison.


Equations (3.52) and (3.Sk) may be combined as


d2 -) K 9 z -A Fjr(Lz, z) + ZHjk(z) 6(z zki (.6
k k


whe re



Z (0) 71 p(T ,W) e a
A ss,(3-57)
aS


(z-zk)

Fjki(zzk) lk e [ "s


(z+zk)2
s
-e


cos a R(xk~zk '


(3.58)

4 (x z). (3.59)


1 (2j 1)nxk
Hjk~ D~ kgC 2a


Equation (3.56) is virtually identical to Eq. (1.30) and has es-

sentially the same boundary conditions, so that the Green's function of








Chapter I may be used to convert Eq. (3.56) to an integral equation. The

appropriate Green's function is determined by


S- Ke G.(z,S) = 8(z 5) (3.60)



and the usual homogeneous boundary conditions. The solution to Eq. (3.60),

from previous results, is



-Kz

e'K cosh Kz
= K, z<5 (3.61)


As before, the integral equation for mj (z) is found by forming
oo co

Sdl G (g,z) [Eq. (3.56)] di 4 l(S) [Eq. ().60)], after reversing
O O

z and 5 in Eqs. (3.56) and (3.60). This result is









= J di G (Llz) Fjk(ilzk)
o k



+ d i G (i~z) Hjk(i) B(L zk Ojl(z) (3.62)
o k








An integration by parts may be performed on the integrals containing the

derivatives and the integral containing Hjk may be completed, resulting in



md (z) = G (zk'z) Hjk(zk) + Aij di G (i'z) Fjk(iizk)
k o k


a G. (4,z) o


(3.63)


Inserting the boundary conditions into the last term gives


oo

'j(z) = G (zk'z) Hjk(zk) + A-Q at G (i,z)


Fjk(S'zk)


gdo [=o


(3.64)


Eq. (3.64) becomes


Inserting A and H.k and assuming that z < z < z 1


M'
k=1


do. (5)
dg
5=o


cosh Kzk 7 cos 2


-Kz
e
K


1 -z
KD Ja~


'^ (z)


k= M '+1


(2j 1)rxx
7k cos 2a


(x,zk,


4 (xk'zk) + cosh Kz


-B ~s
e


o k


Z (0) q p(T ,w)


i-Kz


Fjk(g~zk)


Lt(0) 2! K Ds a t


(3.65)


+ cosh Kz


d 4. ()


Sdge-K Z Fjk(l,zk~
z k





48


It seems appropriate at this point, before continuing into the

mathematical complexities, to simply the notation in Eq. (3.65) and try

to understand its physical significance. To this end call the first term

in Eq. (j.65) H Q, since it represents the contribution which would be
present in the homogeneous case. Also, redefine Eq. (3.58) and (3.59) as


Fjk(z,zk) =ar I(z) W.jk ,xkzk (.6

and


Hjkc(z) = ~Wjk (xk~zk)(.7

where



i(z) e~ [" -'' e (3.68)


and


7k (2j 1)nxk
w.k =-cos (3.69)
Jk 2a


W~jk is simply a weighting factor depending on the distance f~o~m the center
line of the foil location. N (z) is the slowing-down kernel for the parti-

cular geometry chosen, having an obvious resemblance to the infinite-:medium

Fermi-age kernel. Noticing the factor 1/ /T, it can be seen that the

slowing-down kernel obtained for the present geometry resembles the plane

source infinite-medium slowing-down kernel most closely. This is to be

expeccted, since the problem was reduced to a semi-infinite, one-dimensional
case.








Inserting Eqs. (3.66) and (3.67) into Eq. (3.65) and using the

symmetric part of Eq. (3.46) and Eq. (3.61) results in



4 (xz) = W (x H -- -- Wjk G (zk'z) 4 (xk'zk)
j j k

oo

+ Z AJ;J W (x) Wjk Q(xk,zk) S di Mk~(i) G (zk 5) (3.TO)
j k o

where

1 o (2j 1)nx (1
w.(x)= -cs(.1
~ 2a


It should now be noted that the Ath mode of the flux distribution consists

of the following parts: (a) the homogeneous term H e, which represents


the propagation of a neutron wave in a homogeneous system, (b) the absorb-

ing part, the second term, which accounts for the heterogeneous absorption

of the plates, and (c) the last term, which accounts for the production of

fast neutrons by fissions in the plates, and their subsequent slowing down

into the thermal group.

Returning to the mathematical manipulations, it now remains to

evaluate the two integrals in Eq. (3.65), which are complicated by the

fact that K is a complex quantity. The first integral may be written as

(5-zk)2 (S+zk)2
s S
d'i cosh Kg e - e Ejk (j.72)
o k







whe re

jB~ ~ (2 T C ---~ 1)nxk) 1 3


Expressing the cosh as the sum of exponentials and performing the squares
in the exponentials, Eq. (3.72) may be reformulated as


Z'2
k


Sd KEb i -Kg s+Yr
d( e +


gzk/2Te s


-gzk/2T)Ej


Multiplying and combining the various exponentials leads to


-zk 6 rs z
S2 Ejk Sd
k o


[K-(z /21s )]5
-e


-[K-(zk/21s)](


(3.75)


~-lK+(zki/2r1si


By completing the squares in each of the exponentials, Eq. (3.75) may be
rewritten as


(3.74)


-i'T~ [[K+(z /21 )]g











s~K 2r k


s2r
s


-2k 4s Z
e 2 Ejk S
o


dg (e s



- J~ K-


T( Z
1 K-
s 2s/


T ~ Zk
1 K-


2
SJr~

2

S


- ~r


Zk
2s


I K+


+ / K+ s


(3.76)


Each of the four integrals in Eq. (3.76) is of the form



r a2 z T




which becomes, upon changing the variable of integration to


(377)


v l- + Fa ,
2 -s


(3.78)


+~ / K-
( 2s










t 21 Is
JIS
+ Jftr


e "~ dw


- + J ex
S

O


- S
-
e dw
O
(3.79)


2 ;I


-~2
e dw -


Recognizing the two integrals as error functions, Eq. (3.79) becomes


es


+J; / erf




-- + Jr


-



S


es


(3.80)


where the error function is defined as


x

erf(x) 2_ s


-t


(3.81)


dt .


The integrals of Eq. (3.76) now become


(erf [


-SfT a


,] + erf [J t a













k Eak 2 s


K +


Zk

s


2
s


2s


- erf -
s


(3.82)


The second integral in Eq. ( .65) is


n

z


di e- e ~(-g24


(i"*'""Ejk


(3.83)


e


which can be handled in a similar manner to give


2
K- z s


s es



2


SEjk e
k


K


(3.84)


s21c
s

t e5 Z


Z
(erf 2z_ F
2 s


-erf 2 c + / K+k rfK+
[ s s
s s


(erf --z-+[
S


K 1 2 erf J. K -


s


1 rf -z- +
s


-e s (1 erf + K +
-s








If ~ ~ ~ ~ ~ ~ ~ ~ ~ th qaesi h fcosex in Eqs. (j.82) and


(1.84) are performed, a factor exp [zk s] will result which will cancel

the exp [- zk /4s The two remaining factors are exp [+ K zk] and

exp [Ts iE i which can be taken out of the summation and combined with

the factor ex~ [- Bj -s '

Now the expressions for the two integrals, Eqs. (3.82) and (3.84),

may be inserted into Eq. (3.65) to give the complete expression for ". (z).
Thus :


1 -Kz M'(2j 1)nxkc
e cohKz7cos
5 K Ds~ kv k Yk 2a


e-Kz de (5)
oJ,( = K Jtdg




+ (xk'zk) + cosh KE


k= M '+1


(2J 1)n~ exk'*



M Kz zer [ z~~:
k k
k=1 s





s


-k F K



binued on following page)


Z (0) 8 p(T ,, ) e j s -z
Z ( O) 2 K Ds 2


erf + / Kk ~ + 2 erf --k-
S S~


(erf k +~r r K erf z+z









+ os K Z jKz z + z2
+ oh zek efk q
k=1 s


-Kz ( [ z -z (j-1nx
-~~~ ~ ~ e r K-1: Yk cos (2a (xk~zk




Equation (3.85) with Eq. (3.lc6) used for m (xk~zk) is the complete

solution to Eq. (3.52). 4 g(z), the solution to Eq. (;5.53), is similar to

Eq. (3.85) except that K' replaces K, Bj replaces B j, and sin(jrrxk/a)

(2j 1)nxk
replaces cos wherever they appear. Equation (3.85) and the
2a

corresponding equation for O; (z) are inserted into Eq. (3.4~6), which in
turn is inserted into Eq. (3.17) to give the complete expression for 4(r),

the thermal neutron flux at any point in the medium.

It is instructive to see what happens to Eq. (3.85) when the age-

diffusion theory is reduced to one-group diffusion. If no errors have been

made, Eq. (3.85) should reduce to essentially the same form as Eq. (1.66).

To do this, let Tss O and set p(ts,,) = 1. As Ts Ot the arguments of
-k k
all the error functions will approach the form or + -- .Ob-
S S

viously all the arguments will approach +oo on the real axis, depending in

some cases on whether z < zk or z > z The error function for +oo is

simply +1, and so the last tent in Eq. 0.8.5) becomes









iZk 2


-z'/4T
e
2


k= 1


C,( O) 2Ks ~



Z2-il


k
s "


s
-e


M o
-z /4T
s
+ cosh Kz e
k--1


l' z > z
> k
-13 z < zk


7k co 2aS -- "a zk ('6


With a little further manipulation Eq. (g.86) can be further reduced to

the form


M' (2j 1) xx
coshKzkk cos 2a a(xk zk)
k= 1


Zs(0) 9

I (0) D l~


e -K
L K


M
cosh KzK Z k c ( -lrk
+ e Yk cs 2a a(xk zk)
k= M '+1


(3.87)


Insertion of Eq. (3.87) into Eq. (3.85) reduces the latter equation to


2 3o z > zkk






(2j 1)r xk
2a (xk'zk)


2; z > Zk;z








-Kz d@. (g [ Z(0)s) -Kz]
e 30s1 e
jiK df, D -K
F=o Ds ~


~1~ ~CS 2j-lnk cosh Kz
coshKzk cos 2a Q(xk'zk K
k=1


-Kzk (2j 1)icx
e k cos 2a s(xk'zkj ]388
k M'+1


The resemblance of Eq. (3.88) to Eq. (1.66) is obvious. Except for

the differences in the shapes of the perturbing foils assumed in the two

derivations, by following a procedure similar to that described at the

end of Chapter ,II the equivalence of the two expressions can easily be

shown.












CHAPTER IV


CALCULATIONAL RESULTS


Two computer codes were written to obtain numerical results from

the analytical results of the previous chapters. The simpler of the two

uses the results of Chapter I with no additional analytical simplifications.

The second code uses the age-diffusion model developed in Cha~pter III, bu~t

several restrictions had to be made because of computer storage and numeri-

cal problems.

The simpler of the two codes performs the procedures given from

Eq. (1.67) to the end of Chapter I. Computer storage limitations resulted

in several limitations on variable sizes. M, the number of foils, was

limited to 6; using the terminology of Eq. (1.67), L was limited to 10, Q

was limited to 5, and the product LM was limited to 30. The output of the

code consisted of 4 (z) for up to 49 uniformly spaced values of z and the

summation of all the flux modes for the same z values. These results were

given by the code for any number of uniformly spaced frequencies of the

driving source.

Since the one-group, one-dimensional computer code almost completely

filled the computer core, it was obvious that the two-dimensional age-

diffusion results could not be coded in their entirety. Two major simpli-

fications were made. One was to restrict the foil positions and flux

calculational points to the z axis (xk = x = O), thereby transforming the

problem to a one-dimensional problem. This results in the x expansion,








Eq. (3.46), reducing to



(O,z = ~ (z .jr (4.1)


(2j 1)nxk
Equation (3.79) simplifies somewhat since the factor cos 2a

becomes unity, and since O; (z) is not used in Eq. (4.1), the equation cor-

responding to Eq. (3.79) for 4 Q(z) is not needed. Substitution of Eq.

(4.1), with zk replacing z, into Eq. (3.85) along with replacing the cosine
factors with unity gives the basic equation the code must handle as follows:


-Kz de.(g
1 e Ja
QOz = K dg


1
Da
s


M'
Cr-zK Z k cosh Kzk

j k=1



s k e-IZ ~ O,zk )
k-M'+1l

j s -Kz

Zje K 'kj2


cosh Kz
Q(O,zk) K


Zs(0) 2 DsaW



SKzk C" z z

k=1 -' s


efk + K 2 erf 2 + / K
2 s

-Kz [ z -z r
+ek efk -- K erf z k K]
2 sJ 2 s-









2 erf zk + / K + cosh Kz
s k=1

Ez + zk ;K~ -Kzk





r~f kz Zk 7K] -I 1 4(O,zk) (4.2)



The only other simplification that was made was only the inclusion

of the fundamental spatial mode in the source. Formally, this implies that

Eqs. (1.13), (1.22), and (1.24) apply and that S~m 0 for > 1 and m > 1.

In practice, the higher values (x modes) could have been easily included,

but to have included the higher m values (y modes) would have overloaded

the computer core. This is not an insurmountable obstacle, since magnetic

tape storage could have been used but this is a time-consuring and there-

fore expensive approach. Physically, including only th-e fundamental spatial

mode seems justified since recent experimental work at the University of

Florida indicates that a highly thermalized fundamental mode is experi-

mentally feasible (10).


Using Eq. (4.2) rather than Eq. (1.66), the procedure outlined in

Chapter I from Eq. (1.66) to the end of the chapter may be used to obtain

the complete solution to the problem. This will be done later in this

chapter, but, first, attention will be directed to two factors in Eq. (4.2)

which have physical significance that has not been pointed out previously.

If Eq. (3.8) is inserted into Eq. (3.12), with the age to thermal,

Ts, replacing u, the result is











usza(u')du us

p(.ts,w) =e o e o


du'


(4.3)


Noting that the first integral in Eq. (4.5) is the usual resonance escape

probability to thermal, P(Ts) or Plus), Eq. (4.5) may be written as


-iwT


(4.4)


where


u

T =


du'
Szt(u') v(u')


(4.5)


or, in tents of the Fermi age,



Ts SD(T) v()d


(4.6)


Obviously, Ts has the units of time, and represents the time required for

the neutron to slow down to thermal energy. The factor exp [-i Ts] then

represents the phase shift introduced by the slowing-down time.

The second factor in Eq. (4.2) which can be simplified is

exp[(K~ B j)Ts]. With the use of Eq. (J.50), this immediately becomes


Z + -
as


(K -B j)"
e


( *7)


as s
D
= e s


v D
s s 48








The term in the first exponential is the ratio of the age to thenaal to

the square of the diffusion length, .rs/L Remembering that for a point

source L2 is one-sixth the mean square distance that a thermal neutron

travels before being absorbed and that T(u) is one-sixth the mean square

distance that a neutron travels in attaining lethargy u, the ratio can be

seen to be a measure of the distance traveled while slowing down relative

to the distance traveled at thermal energy. Typical values for this ratio

are 0.13 for graphite and 0.008 for D20.

Now returning to the actual calculational procedure, as before, the

expansion, Eq. (4.1),must be truncated; i.e.,



i (O1,z = 1 (z) (4.9
j=1

Because of the restriction of the source to the fundamental mode only,

Eq. (4.9) simplifies further, since can only assume the value 1.

Restricting j to 1 in Sj does not, of course, eliminate the higher j com-
ponents in the flux since the perturbing foils generate higher x modes but

they cannot generate higher y modes since the foils are assumed to extend

the full length of the medium in the y direction. Equation (3.85) with


1 i .( )Ji, inserted for 4 (x z ) and cos ( -)' replaced by unity
--- 01k kk 2a


may now be written JM times for each of the J modes and for each of the M

foils. This set of equations may then be solved simultaneously for the JM

values of mjl(zk~). Writing these values as a vector,


T 11zj 1(z ""11(M 21zj"" J(zM] ,(4.10o)






63


the set of JM equations may be put into matrix form, as before. That is,


where B is the vector of source tents, [I] is the identity matrix, and [A]

is the matrix containing the Green's functions connecting the foils.B

may be written as


(4.12)


where


S1
11 2DKl


(4.13)


it into J2 submatrices,

given by


The [A] matrix is defined, as before, by

each M by M. The n,p element of the I,K


partitioning

submatrix is


(A = p z ) coh K z -KIlzn- R)
(I,K n,p = Il p) Ch Il p IlI n


= p z )coshK z e-KIlz (R,
BIl p) O' KIl n AlI


,(4.14)



S(4.15)


for n > p



for n < p


where


p ( ) -- --
BIl p) KIlDa '


(4.16)


10
I + -
as Yg
D ~s
'1 P(r ,O) e s
s2 D~a K


11 =(0)
t


[A] 0 + [I] 0 = B ,


(4.11)


T 11 1 11 M
B = 1 e ... a11 e ,O..O ,


(4.17)





64



(RIi n~pe~Ki~~i np 2 (Al) (A2) +2(Ai)j


-K1(z +z )


7p(,,
(A )


- (A5)


+ cosh KYlzn Ip !e Ilp [(A2)


Z Z
n p
S


z + z
n p
2, -
S


-- K ,


(Al) = erf




(A2) = erf




(A3) = erf




(A4) = erf




(A5) = erf




(A6) = erf


(4.19)




(4.20)




(4.21)




(4.22)




(4.23)




(4.24)


Pz JIK1
2 c- - + / K,
2, s I


z z
n p



z + z
n p
2#


-- K ,


s 1~n,




-,~ K1l


z
P
2J-r
S


The vector 0 is obtained from


L = A] [I] B .


(4.25)


- 2(A6))


- 11 -e I (h ]


,(4.18)









The elements of 0 are inserted into Eq. (3.85) with the previously mentioned

simplifications, which is in turn inserted into Eq. (4.9) to give the com-

plete answer to 4 (O,z). The entire thermal neutron flux is obtained by

using Eq. (3.17), which gives



s1 1
j=1


Since the representation of the flux is by expansion in eigenfune-

tions, an important characteristic of the numerical calculations is the

rate of convergence with increasing numbers of the eigenfunctions, or

spatial modes. Figures 2 and 3 show flux amplitude versus distance along

the z axis with 1, 2, 4, and 10 spatial modes and with the source modulated

at 100 and 1000 eps, respectively. The model is the one-group representa-

tion of Chapter I with the geometry shown in Fig. 1. One cadmium strip,

0.0508 x 1 x 71.12 em, is positioned 8.89 cm from the end of a graphite

assembly 44.45 x 71.12 em, the dimensions of one of the experimental

facilities available at the University of Florida. In addition to the

curves shown in Figs. 2 and 3, calculations were performed with 6 and 8

spatial modes. At the position of slowest convergence, which is in the

neighborhood of the cadmium foil, the ratios of the 6- and 8-mode results

to the 10-mode results are approximately 15 and 6/o, respectively, for both

100- and lOOO-ops source frequencies. Figure 4 shows the phase lag versus

z for 100 ops, corresponding to Fig. 2. In this case the change from 6 to

10 modes and from 8 to 10 modes is approximately 5 and 2/o, respectively.

Figure 5 is the phase lag versus z for 1000 ops, corresponding to Fig. 3.

The convergence is faster for this case, approximately 2.5 and 1o for 6 to






































r \ 40 MODES
z 4 MODES
o
2 MODES

J (MODE


U.U1
0 5 r0 15 20 25 30 3
DISTANCE ALONG I AXIS (cm)

Fig. 2. One-Grroup Flux Amplitude vs Distance Along the Z Axis for
One Cadmium Foil in Graphite, at 100 ops, with 1, 2, 4, and 10 Sp~atial
Modes.


rO.O






67






































0.01











Modes.
























10 MODES
4 MODES
2 MODES


4 MODE

































O S l0 13 20 25 30 31
DISTANCE ALONG 7 AXIS (cm)


0.7





0.6








0.1





S0.3












S0.


Fig. 4. One-Group Flux Phase Lag vs Distance Along the Z Axis for
One Cadmium Foil in Graphite, at 100 eps, with 1, 2, 4, and 10 Spatial
Mode s .






69































I~~ MODESD
4 MMODE
















a 5 10 15 20 25 30 35
DISTANCE ALONG I AXIS (cm)

Fig. 5. One-Group Flux Phase Lag vs Distance Along the X Axis for
One Cadmium Foil in Graphite, at 1000 cps, with 1, 2, 4, and 10 Spatial
Modes.









10 and 8 to 10 modes, respectively. Calculations were limited to 10 spatial

modes because of machine storage limitations. Although the convergence is

somewhat slow, it is probably adequate for most purposes, especially since

the problem used is a rather severe test.

Figures 6-10 show the results of calculations in the same geometry

as described above except that two perturbing foils are present, at

z = 8.89 and 17.78 cm. These foils are either cadmium, as used previously,

or U2ss of the same total absorption cross section as the cadmium strips.

The calculations with cadmium are one-group and with U2ss are age-diffusion.

Figure 6 shows flux amplitude as a function of z for one spatial mode cal-

culations using cadmium foils with source frequencies of 200, 500, and 800

eps and using U~ss foils at 200 eps. The unperturbed flux at 200 eps is

also plotted on Fig. 6. Figure 7 shows the same data using six spatial

modes. The unperturbed flux is identical in Figs. 6 and 7 since the

assembly is driven with fundamental mode only in each case. Figure 8

presents the phase angles as a function of z corresponding to the amplitudes

shown in Figs. 6 and 7. The one-group, six-spatial mode calculations shown

in Figs. 7 and 8 are replotted as a function of source driving frequency in

Figs. 9 and 10. Figure 9 presents flux amplitude versus frequency for z

values of 2, 8, 8.89, 15, and 3O em. Figure 10 presents phase lag for the

same conditions.


A second facility that is available at the University of Florida for

measurements compatible with the geometry restrictions of these calculations

is a water, or D20, tank which is 120.97 x 130.17 cm in the transverse

directions. Figure 11 is a plot of flux magnitude versus z at 200 eps

for an age-diffusion calculation in D2O with two 1-in.-dina natural uranium











































r 200 cps, UNPERTURBED



200 cps, AGE DIFFUSION







\ \ 200 cps



'500 cps >ONE GROUP



'800 cps




) (O 20 30 40 50 60 7
DISTANCE ALONG z AXIS (cm)


0.0(














0.00~


Fig. 6. Flux Amplitude vs Distance Along the Z Axis for Two Cadmium
or Uranium Foils in Graphite, at 200, 500, and 800 ops, with One Spatial
Mode .





















































0.0r













0.00(
O (0 20 30 40
DISTANCE ALONG I AXIS (cm)

Fig. 7. Flux Amplitude vs Distance Along the
or Uranium Foils in Graphite, at 200, 500, and 800
Modes .


50 60 70


Z Axis for Two Cadmium
cps, with Six Spatial





75











6 MODES~MD 0 p









800 cps
6 MODES




6 MODES GE


3 4 MODE














O~~~ MODE 20304 50 cps









DISTANCE ALONG z ALIS (cm)

Fig. 8. Flux Phase Lag vs Distance Along the Z Axis for Two Cadmium
or Uranium Foils in Graphite, at 200, 500, and 800 cps, with One and Six
Spatial Modes.



















z= 2cm













8cm







8.89 c

















30 cm


0.(





0.05







0.02





0.0(


400 600
SOURCE FREQUENCY (cps)


4000


Fig. 9. One-Group Flux
Foils in Graphite, at z = 2,
Modes .


Amplitude vs Source Frequency for Two Cadmium
8, 8. 89, 1~ and 3O cm, with Six Spatial












4.0









3.5









3.0




z = 30 cm



2.5





2.0



cc













S20.




O
c 0 0 0080(O
SOREFEUEC cs

Fi.1.OeGopFu hs a sSuc rqec o w amu
Fol nGahta ,8 .89 5 n Ocwt i pta
Moe











(0.0









5.0












2.0 ~- v

-J 6 MODES


OO
4 OD












0.5












0.2










O 40 20 30 40 50
DISTANCE ALONG I AXIS (cm)


Fig. 11. Age-Diffusion Flux Amplitude vs Distance Along the Z Axis
for Two Uranium Foils in 'D0, at 200 eps, with One and Six S~patial MIodes.





77


rods 130.17 cm long located at z = 10.998 ad- 21.996 em. The obvious

characteristic of the calculation is that the natural uranium rods do not

offer sufficient perturbation to produce dramatic effects. An alternate

and feasible perturbation could be obtained by using single plates from

fuel elements of the MTR type. One of these plates has approximately the

same total absorption cross section as the natural uranium rods. Therefore,

the perturbation would be approximately the same. In order to produce a

significant perturbation with fuel, fairly large U2ss rods would have to

be used.

Since the calculation shown in Fig. 11 has small perturbations, it

presents a good opportunity to show the relative contributions of the

unperturbed, thermal-neu~tron diffusion, and slowing-down contributions to

the final answer. These are shown in Fig. 12, where the three curves cor-

respond to the three terms in Eq. (3.85). Since it is a six-spatial mode

calculation, the thermal-neutron diffusion and slowing-down curves are sums

of the second and third terms in Eq. (3.85).





































UNPERTURBE )
























THERMAL DIFFUSION
r I SLOWING DOWN
FOIL POSITIONS

O 10 20 30 40 51
DISTANCE ALONG r AXIS (cm)


Fig. 12. Components of Age-Diffusion Flux Amplitude vs Distance
Along the Z Axis for Two Uranium Foils in D20, at 200 eps, with Six
Spatial Modes.













CHAPTER V


SUMMARY AND CONCLUSIONS


The goal of this dissertation has been to predict, as accurately

as possible, the behavior of neutron waves in heterogeneous multiplying

and nonmultiplying media. The primary approach has been similar to the

heterogeneous reactor theory of Feinberg and Galanin in that the continu-

ous slowing-down model has been assumed to describe the behavior of the

neutrons above thermal energy and diffusion theory has been assumed to

describe the behavior of the thermal energy neutrons. Rather than assume

that the infinite-medium slowing-down and diffusion kernels are applicable,

a more basic approach was chosen which resulted in slowing-down and diffu-

sion kernels for the particular finite geometry that was assumed.


The preceding chapter showed some results of calculations using

this theory. In this chapter, it will be shown how this theory may be

applied and extended.


Application to Critical Reactors


As stated previously, the theoretical developments of this work are,

in several ways, an extension of the Feinberg-Galanin heterogeneous reactor

theory. Aside from the obvious feature of predicting the behavior of neu-

tron waves in subcritical assemblies, this work takes into account the

finiteness of the assembly. The effect that this has is best shown by

applying the theory to a critical reactor, which can be done by dropping





80


the source terms and changing the boundary conditions slightly. The ab-

sorptive Green's function defined by Eq. (3.60) will have the same boundary

conditions as before except that now G.(z,S) itself, rather than the deriv-

ative of G.(z,g), will be required to go to zero at z = 0. This change

in boundary condition simply results in G.(z,S) changing to


e sinh KS 51
G (z,S) = K z > E,(5


e'K sinh Kz
= K z < 5 (5.2)


With the new Green's function, Eq. (3.65) applies to the critical assembly

with the added simplification that all the boundary conditions given by

the last term in Eq. (3.65) vanish. Equation (3.65) now becomes

oo

O (z) = G (zk'z) Hjk(z) + Aid di G (i,z) Fjk(II.Izk3
k k o


where G.(Stz) is now given by Eqs. (5.1) and (5.2). Remembering that the

factrs H and F.k both contain

Jkak




it can be seen that, for k = 1,2,...,M, and j = 1,2,...,J, JM homogeneous

equations are available for each value of a. For the system to be critical

the determinant of the coefficients of the 4 Q(zk) must be zero. This of

course gives an equation which can be solved f~or rl (contained in A j) or
whatever other eigenvalue is chosen to determine criticality.








The first difference noticed in comparison with the Feinberg-Galanin

theory seems to be simply additional complexity since the Feinberg-Galanin

theory only requires the solution of an M by M determinant. The improve-

ments are the improved kernels contained in Eq. (5.5). As pointed out in

Chapter III, Aie and Fjk contain the finite-medium slowing-down kernel

rather than the infinite-medium kernel used by Feinberg-Galanin. The

Green's function G. also is an improvement, being the finite-medium diffu-

sion kernel for the geometry used in this work, rather than the infinite-

medium diffusion kernel used by Feinberg-Galanin.


Conclusions


The theoretical developments of this dissertation probably represent

a model which is as complex and physically realistic as is practical at this

time. Further improvements, in order to be practical for machine calcula-

tions, would certainly be desirable if they simplified the mathematics or

if they eliminated remaining physically undesirable assumptions, or both.

If such improvements are not easily attainable, the primary question remain-

ing is how to use the theory best to direct future experimental work in this

area.


Several approaches are possible. One is simply to continue full

calculations of the thermal flux in various configurations, such as those

presented in the preceding chapter. This is a rather brute force technique

but at least has the advantage that calculations are cheaper and much more

easily accomplished than experiments since the machine codes are available.

Some possibly more satisfying approaches can be suggested by re-examination

of some of the results of Chapter III. To this end, write Eq. (3.64) as









1 L k (2 1)x
me (z)= H G (z,zk) cos 2a (xk'zk)




Ee Z 7 (2j 1)TIxk oc4r iq()0G(k~
k o



where

-Kz do. (5)
H e (5.6)
je K dS
5=0

and
-B .1s
I (0) a p(T ,W) e a
j I (0) 2D D57
t s


atnd Gj and Mk are defined by Eqs. (j.61) and (3.68), respectively. A cor-

responding equation exists for Qj (z); thus


1~~ YkJxk
mj (z) = H; zDk sin -ak (xk~zk)



+ sin ~ka (xkzk) di Mk(i) G (zk'n ,(58
k o


where G!, H; and E'. are obtained by replacing B with B'~ and K~ with
K'2 where they appear in Eqs. (3.61), (5.6), and (5.7), respectively.

Insertin- Eqs. (5.5) and (5.8) into Eq. (3.46) and in turn into Eq. (3.17)
results in








1 (0 -1)n 1 (2j 1)nx




+ 1- sin janx H', .zz)C



+ G! (z'zk) Sjk] (xk~izk) + Ea Cjk e(xk'zk)



diM~)G. (zk,8)~ +~ E Sjk e(x'zk S at Mk(i) Gr (zk^>
o k o
(5.9)


whe re

7k c (2j 1)nx (2j 1)nxki
C.E--cscos (5.10)
jk a 2a 2a


and


S. Y sin sin -rx (5.11)
jk a a a


Experimental measurements of thermal neutron flux as a function of position

and frequency represent the quantity calculated in Eq. (5.9). These data

may be operated on with


a b
S~ (2j 1)rtx (2 1)n
2a 2t
-a -o


or with
a b

in sing jx dy,~ cos (2R 2bl)r[y (515)
-a -b









to yield experimental measurements of the flux components calculated,

respectively, in Eqs. (55) and (5.8). This entire procedure may be per-

formed experimentally with the moderator alone, with fuel plates or rods

and with purely absorbing plates or rods of the same macroscopic absorbing

cross section as the fuel assemblies. The technique of replacing fuel as-

semblies with comparable pure absorbers has been utilized previously to mea-

sure effective delayed neutron fractions (11). The six terms in Eqs. (5.5)

and (5.8) may now be separated. The homogeneous experiment yields Hj and

H' Subtracting the homogeneous term from the pure-absorber data yields
the absorption terms, the second terms in Eqs. (5.5) and (5.8). Similarly,

the difference between the pure absorber and the fuel assembly data yields

the last terms which contain the slowing-down kernels.

Since the above experimental data would presumably be obtained by

cadmium-difference measurements with a neutron detector which has approxi-

mately a 1/v response, it is also of interest to calculate the epicadmium

flux or cadmium-covered detector response. Equation (g.44) may be rewritten

as


Z (0) Z k" t






e Cjk + e Sjk (xk'zk) (5.14)


The argument r has been added to Mk as a reminder that Mk is defined as in

Eq. (3.68) except that I replaces Ts. Insertion of Eqs. (5.14) into Eq.
(g.18) yields








Z (0)
s(o
if


q(r,, ) = p(T~W)


c (2 1)ny
cos 2b


Sjk


Mk~( z,T)
j k


a(xk)zk)


(5.15)


Cjk


The response of the 1/v detector covered


with cadmium can be taken to be


proportional to


ucd(r, ) 1 du
S ed


= 2E
m


(r, )
epi


u

= edq(rT) eu/2 du ,


(5.16)


(5.15)


where ucd is the lethargy of the
into Eq. (5.16) gives


cadmium cutoff.


Insertion of Eq.


z (0)
s ~0
t


S1 c (20 1)Iy
=' cos 2b


epi(r )


(5.17)


(5.18)


and


e ~j


Z Iejk Cjki jk Sjk ke(i~zk)
j k

whe re

ud -B~j .u/~
Iejk ~Cp(u,w) N ~(z,u) e e/ u









ucd -B' u
I~k pu~)Mkz~) u./2 du (5.19)



As before, experimental measurements of the quantity represented by Eq.

(5.17) taken as a function of position and frequency may be operated on

with the operators given in Eqs. (5.12) and (5.13). Let the result of

operating with Eq. (5.12) on Eq. (5.17) be called F'Ajk, which is


2 (0) 7 (2j 1) xk
Fejk Zs0) jk I-cos 2a (xkr~zk). (5.20)



The nonsyrmmetric part may be called F~ijk and is obtained by operating with

Eq. (5.13) on Eq. (5.17), giving



F'Ajkt Z jk sin k +(xk~zk) (5.21)


Several things should be noted at this point. First, all the inform-

ation obtained thus far could have been obtained from experimental data

taken at only one z value. Secondly, each term obtained experimentally

except for the homogeneous term has contained a weighted summation over the

thermal neutron flux at the surface of the foils. Alternatively, by com-

bining the thermal constant, yk, with the flux, these summations may be
considered to include the net current into each of the foils.

Assume now that the experimental data and following analysis have

been performed for M values of z, different from the locations of the M

foils, of course. These M values of z should be chosen so as to be fairly

close to the M foil locations in order to maximize the information about





87


the foils in the M sets of measurements. Now assume that this set of

information is available for the synmmetric absorption term, the second

term in Eq. (5.5). A set of M inhomogeneous equations results, with the M

unknowns being the net current into each of the M foils or, if 7k is known,

being the surface flux at each of the M foils.

Now that the net current into each foil is known, this information

may be inserted into some of the other terms to obtain additional informa-

tion. For instance, having the net currents and experimental values of

F~jk at M points, the integrals Iljk defined in Eq. (5.18) may be obtained
by solving the set of M equations given by Eq. (5.20). Alternatively, if

the E .j's are assumed to be known, the integrals over the slowing-down
kernel in Eq. (5.5) may be obtained for each foil.












APPENDIX A


THE COMPLEX ARITHMETIC VERSION OF GINV


GINV is a F RTRAN subroutine developed by Burrus et al. (12) for

the purpose of obtaining working inverses to singular matrices or to matrices

which appear singular to a digital computer. It was modified to handle

complex arithmetic and was used in the programs employed for calculations

in this dissertation because, at one point, it appeared that the matrices

in these problems were very ill conditioned. Although this turned out not

to be the case, it was left in because it takes up little more storage than

a standard matrix inverter and gives slightly better answers in practically

all situations.

Theoretically, the inverse of a matrix A, A l, has the property that

A- A = I, where I is the identity matrix. This property never is true

precisely, when an inverse is obtained on a digital computer with finite

accuracy. GINV operates so as to find the working inverse A defined as

that matrix for which the root mean square value of (I A A) is a minimum.

The theory of operation will not be given here since it is not

pertinent to this dissertation and should be available soon in the

literature (13).

Pages 90-92 contain the F RTRAN listing of the program as modified

for complex arithmetic. Sufficient comment cards are included to enable

a competent programmer to modify the program for ordinary F RTRAN arithmetic.









A few comments may be appropriate on the use of the program. First,

when TAU is set to zero, the routine performs virtually the same operations

as other matrix inverters and should give identical results on well-

conditioned matrices. It may also be noted that if the problem to be solved

is to obtain x from Ax = b, where A is an m x n matrix (with n < m) and b is

a given m component vector, GINV automatically obtains the least-squares

solution to the set of n equations. It is obviously much more versatile

than most matrix routines, although its full versatility was not needed


in this work.








SUBROUTINE GINV(A,MR,NRNGIfAU,UATEMP)

C UPON ENTRY A = A MATRIX WITH NR RUWS ANO NC COLUMNS
L AFTER THE RETURN THE ORIGINAL A IS DESTROYED ANU THE
C ARRAY A CONIAINS THE TRANSPOSE OF T, THE TRANSFURMATlON
C MATRIX (X = T+B) WHICH SELVES THE PAIR OF EQUATIONS

C A+X = 8 ANO TAU leX = 0

C IN THE LEAST SQUARE SENSE.

C MR = IST DIM. NO. OF ARRAY A IN THE CALLING PROGRAMS
C FAU IS A NON-NEGATIVE CONSTANT WHICH CONTROLS THE ERROR
C PROPAGAllONN OF THE TRANSFORMATION MATRIX. IF FAU = 0.
C THEN THE RESULTING TRANSFORMATION MATRIX IS THE
C ORDINARY LEAST SQUARES TRANSFORMATION MATRIX.
C (THE INVERSE IF NR = NC)
C ASTEMP AIJo U ARE USED FOR WORKING SPACE 8Y THE ALGORITHM
C AND 00 NOf NECESSARILY CONTAIN ANY RELEVANT NUMBERS AT
i IHE CONCLUSION. ATEMP MUST BE DIMENSIONED AT LEAST NC
C t)Y THE CALLING PROGRAMS. U MUST BE DIMENSIONED AT
C LEASE NC*NC UY THE CALLING PROGRAMS.
C NO. OF MULTIPLICATIONS = NC+*2 (5/2 NR + 2/3 NC)

I DIMENSION A(30,30), U(30,30), ATEMP(30), 0001l),
L DUT2(1), TAUSQ(1), DUM(1)

C THE DIM. ST. AS SENT HY RO\SS 8URRUS FOR FO~RRAN IV WAS
DIMENSION A(MRNC),U(NCINC),ATEMP(NL)

C NOTE THAT IN COMPLEX ARITHMETIC VERSION, IT HAS 8EEN
C ASSUMED THAT MR IS ALSO fHE 2NO DIM. NO. OF A AS WELL
C AS THE 15T.
TAUSQ(1) = CAU**2
TAUSQ(2) = 0.0
C PLACE UNIT MATRIX'IN U
00 5 I = LNC
00) 4 J = 1,N~C
1 4 U(IJ) = (0.0,0.0)
I U( )= .0 0 0









C URTHOGONALILL COMBINED MATRIX (A ABOVE U) BY GRAMM-
C SCHMIDT-HIL~ftRT METHOD WITH FIRST NR ROWS WEIGHED
C wifH L ANO THE OTHER NC ROWS WEIGHTED WITH 1/TAU. THEN
SREORTH-OGONALLZ TO LESSEN ROUNOUFF ERROR.

00 20 1 =- 1,NC
IMR = I + MR
II = I 1
IF (II) 2,11,2
2 011 10 LL = 1,2
00 10 J = ,II
JMR = J+MR
I DOT =(0.0,0.0)
1 00T2 =(0.0,0.0)
UO 3 K = 1,NR
DUM(1) = A(KJ)
DUM(2) = -A(KJMR)
1 3 DOT = A(KI)*DUM + UOT


00 6 K = 1,J
OUM(1) = U(K,J)
DUM (2) =- -U(KJMR)
I 6 OOT2 = U(K,1)*0UM + OUT2
I UOT = 00T + 00T2*+TAUSQ
00 8 K = 1,J
1 8 U(KI) = U(K,1) 00T*U(KIJI
00 10 K = IrNR
I 10 A(KI) = A(KI) 00T*A(KJ)

C NORMALI:L THEf COLUMN I OF THE COM81NED MATRIX

1 11 007 =(0.0,0.0)
I DOT2 =(0.0,0.0)
DO 12 K = 1,NR
12 UUf = 00T + A(KI)**2 + A(KIMR)**2
00 14 K = 1,I
14 0012 = 00T2 + U(KII)**2 + U(KIMR)**2
00f -- OUT + 00T2*TAUSQ
UUT = SURTF[00T)
00 17 K = 1Ir
1 17 U(K,1) = U(K,I)/UOT
00 19 K = 1,NR
I 19 A(KI) = A(KI)/DOT
20 CONTINUE




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs