Title Page
 Table of Contents
 List of Illustrations
 List of Tables
 Development of the homogenous medium...
 Angular integration
 Speed dependence
 Difference density method
 Plane slab geometry
 Comparison of n and d methods for...
 The first-flight green's function...
 Solution of the transport...
 Results and conclusions

Title: Spatially dependent integral neutron transport theory for heterogeneous media using homogeneous Green's functions.
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00091574/00001
 Material Information
Title: Spatially dependent integral neutron transport theory for heterogeneous media using homogeneous Green's functions.
Series Title: Spatially dependent integral neutron transport theory for heterogeneous media using homogeneous Green's functions.
Physical Description: Book
Creator: Church, John Phillips,
 Record Information
Bibliographic ID: UF00091574
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 - 000551294
oclc - 13311489

Table of Contents
    Title Page
        Page i
        Page ii
        Page iii
    Table of Contents
        Page iv
        Page v
    List of Illustrations
        Page vi
        Page vii
    List of Tables
        Page viii
        Page 1
        Page 2
        Page 3
    Development of the homogenous medium green's function method
        Page 4
        Page 5
        Page 6
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
    Angular integration
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
    Speed dependence
        Page 28
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
    Difference density method
        Page 36
        Page 37
        Page 38
        Page 39
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
    Plane slab geometry
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
    Comparison of n and d methods for slab geometry
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
    The first-flight green's function for a homogenous medium
        Page 69
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
    Solution of the transport equation
        Page 102
        Page 103
        Page 104
        Page 105
        Page 106
        Page 107
        Page 108
        Page 109
        Page 110
        Page 111
        Page 112
        Page 113
        Page 114
        Page 115
        Page 116
        Page 117
        Page 118
    Results and conclusions
        Page 119
        Page 120
        Page 121
        Page 122
        Page 123
        Page 124
        Page 125
        Page 126
        Page 127
        Page 128
        Page 129
        Page 130
        Page 131
        Page 132
        Page 133
        Page 134
        Page 135
        Page 136
        Page 137
        Page 138
        Page 139
        Page 140
        Page 141
        Page 142
        Page 143
        Page 144
        Page 145
        Page 146
        Page 147
        Page 148
        Page 149
        Page 150
        Page 151
        Page 152
        Page 153
        Page 154
        Page 155
        Page 156
        Page 157
        Page 158
        Page 159
        Page 160
        Page 161
        Page 162
        Page 163
        Page 164
        Page 165
        Page 166
        Page 167
        Page 168
        Page 169
        Page 170
        Page 171
        Page 172
        Page 173
Full Text







April, 1963

To my former teacher

Dr. Dudley E. South


The author wishes to express his appreciation

to his chairman, Dr. G. R. Dalton, for his guidance and

assistance throughout the course of this investigation.

In addition, the author is indebted to the

staff of the Computing Center of the University of Florida

for their assistance, and to the Computing Center and the

Department of Nuclear Engineering for their combined

financial support of the computational work for this


The author has been ably assisted in the prepa-

ration of this dissertation by his wife, Nancy, who typed

the numerous drafts and the final manuscript; its final

form owes much to her concern. It is difficult to express

in words the author's appreciation of her faith in him.

She has made the path to his goal less difficult.




LIST OF TABLES . . . . . . . .

.N1 TR .W I A
















S. . iii

* . vi

S. . viii

. . . . . . . # .



Expansion of Angular Dependence in
Spherical Harmonics and Zeroth
Harmonic Approximation ..

Advantages and Disadvantages of
Homogeneous Green's Function Method


General Speed Dependent Equation .

The Source Term . .. . .

One-Speed Equation . . .


Development of Difference Density
Equations . . . . . . .

Equivalence of Equations for N and D


One-Dimensional Problem . . .

Symmetry Considerations . . .

* J.

S 4

S 15

S 15

S 26


S 28

S 30

S 34

S 36






VII. 1








IX. 2









Differential Equation and Boundary
Conditions . . . . .

Analytic Form of G(v;x x') . ...

Monte Carlo Generation of G(v;xlx')


Spatial and Speed Dependence . .

Iteration Techniques . .. . .


Comparison of the Homogeneous Green's
Function Method with High Order P
and Sn Methods . . . . .

Monte Carlo Generation of the
Green's Function . . . . .

Iteration Convergence . . .

Generation of the Green's Function
for the Multispeed Problem . . .

. . . . . . . . *

DERIVATION OF Gf(v;xIx') . .


4 0 * 0 0 0 . . . .

S 69

* 77

S 90

. 102

. 102

. 115

. 119









* 0

* 0


Figure Title Page

1 The Spatial Dependence of the h Functions
of Equations (2), (7), and (8) . . . 13

2 Plane Geometry Transformation . . . 50

3 Symmetry of Neutron Flow'at the Boundaries
of Unit Cells in An Infinitely Repeating
Lattice . . . . . . .. . . 75

4 Slab Unit Cell with Mirror
Boundaries at x-+b . . . . . .. 85

5 Mirror Symmetry of First-Flight Green's
Function About Plane at x-0 . . . . 97

6 Average Values of the Whole-Cell and
Half-Cell Green's Functions for
Problem 2, Appendix B . . . . . 128

7 Average Value of the Whole-Cell and Half-
Cell Green's Functions for Problem 4,
Appendix B . . . . . . . . 129

8 Effect of Convergence Acceleration
Techniques Applied to Problem 2 of
Appendix B. Solution Iterated Until
Residuals Were Equal to 10-5 of the
Neutron Density at Each Spatial Point . 133

9 Spatial Dependence of Scalar Neutron
Density for Gold-Graphite Cell of Problem 1. 148

10 Spatial Dependence of Scalar Neutron
Density for Fuel-Water Cell of Problem 2 .155

11 Advantage Factor for Fuel-Water Cell of
Problem 2 by Various Approximation Methods 156

12 Advantage Factor Calculated by HGI
Method vs. Number of Iterations for
Fuel-Water Cell of Problem 2 . . . . 158

13 Advantage Factors for Problem 3 by
Various Calculations . . . . . . 161

14 Advantage Factor for Problem 4 by
Various Calculations . . . . . .. 165

15 Spatial Dependence of Scalar Neutron
Density for Cases 1 and 2 of Problem 4 . 167

16 Spatial Dependence of Scalar Neutron
Density for Cases 3 and 4 of Problem 4 . 168






Number Title Page

1 Relative Effort for Methods of Solution
of End-Point and Midpoint Equations for
a Total of M Spatial Intervals . . . 109

2 Breakdown of Computation Time for
Problem 2, Appendix B . .. . .. 134

3 Parameters for Problem 1 . . . . 147

4 Parameters for Problem 2 . . . . 153

5 Parameters for Problem 3 .. . . 159

6 Parameters for Problem 4 . . . . 163



Probably the most fundamental information needed

to describe static behavior of a nuclear reactor is the

speed and spatially dependent neutron collision rate,

that is, Z(r,v)N(r,v)v. Such information is needed in

order to predict long-term kinetic behavior of the reactor,

burn-up, power level, etc. Thus.one needs to know the

neutron density, N(r,v), over the entire speed spectrum

and space of the reactor.

Because of chemical binding and crystalline

effects, the scattering kernels are considerably different

for neutron speeds near thermal than for the higher speeds

implied when one speaks of epithermal and fast speed

ranges. In general, the absorption cross-sections of

reactor fuels also have a different functional dependence

on neutron speed at thermal speeds than at higher speeds.

Finally, because of the tremendous range of speeds in

reactor spectra, and because a critical reactor with no

independent sources is an eigenvalue problem if the entire

speed spectrum is considered, the problem of predicting

the neutron density is usually divided into several speed


The present paper deals exclusively with the

prediction of the spatial and speed dependence of the

steady-state thermal neutron density in heterogeneous

media, although the basic method to be presented is

equally valid for any speed range.

There have, of course, been numerous approaches

to various aspects of the problem which can be very

loosely divided into two classes according to whether

the integro-differential or the integral form of the

transport equation was solved. The methods of solving

the integro-differential equation differ to a large extent

only in the choice of functions used to expand the angular

dependence of the terms in the equation. The expansion

functions include spherical harmonics (1-3), Tchebycheff

polynomials (4-7), and first order polynomials, that is,

trapezoids (8, 9).

The methods of solving the integral equation

have been a bit more varied, due perhaps to a more

intuitive derivation of the equation to be solved. As

Osborn (20) has pointed out, some approaches, particularly

as applied to thin absorbers (10, 11), have separated the

problem into two parts--first, a calculation of the

absorption rate in the absorbing media in terms of the

density at the surface of the absorber and second, a

density depression problem of evaluating the density

at the surface of the absorbing media in terms of the

density in an all-moderating media. Other integral

approaches (12-14), based on transport theory, depend

on the fact that the absorption rate depends only upon

the zeroth angular moment of the neutron density, the

solution for which is independent of higher angular

moments if scattering is isotropic.

The present paper describes the latter type

of approach. In brief, the method to be presented uses

a kernel for a homogeneous finite medium to formulate an

integral transport equation which predicts the speed and

spatially dependent neutron density in a heterogeneous,

finite unit cell.



Assume that the neutron density in a general
three-dimensional medium satisfies the time independent

linear Boltzmann equation, and let N(r,vn)drdvd. denote

the steady-state angular neutron distribution, that is,

the number of neutrons in the volume element dr about r

having speeds in dv about v, and whose directions of

motion lie within the solid angle dn about n; N(r,vn)

will subsequently be called the angular density.

Define V(r,v) to be the probability per unit

time that a neutron of speed v will suffer a scattering

collision in dr about r with a nucleus which is in

motion with respect to the laboratory system. Further,

define y(r,v) to be the probability per unit time that

a neutron of speed v will be absorbed in dr about r by a

nucleus which is in motion with respect to the laboratory

system. Finally, define P(r,v-,v4*-nn) dv dgto be the

probability per unit time that a neutron of velocity v'r.'

will suffer a scattering collision in 4r about r with a

nucleus which is in motion with respect to the laboratory

system such that the neutron emerges from the collision

with speed in dv about v and direction of motion in dn

about n. Then let S(r,vr)drdvdn denote the source term,

that is, the number of neutrons per unit time being added

to the volume element dr about r, having speeds in dv

about v, and whose directions of motion lie within the

solid angle dn about n.

The transport equation which takes into

account the relative velocity between the neutrons

and collision nuclei which are in thermal motion can

then be written as

vn 7-N(r,vn + [V(r,v)+ A(r,v)] N(r,vr

S(r,vrn) + fdv*d' P(r,vN(v,-. N(r,v'rD. (1)

Equation (1) is an integro-differential equation; to

make the problem of its solution determinate it is

necessary to specify boundary conditions. The boundary

conditions follow from the physical interpretation of

N(r,vn). For example, at an interface between two media,

any packet of neutrons characterized by the vectors

r + Rn and vn will contain exactly the same number of

neutrons when it enters one medium as when it left the

other; in other words, N(r + Rn,vn) is a continuous

function of R for r + Rn at the interface. For a

complete discussion of this and other boundary conditions

applied in neutron transport theory, the reader should

consult Davison (2).

Consider the space of interest to be divided

into R regions, the medium of each region being assumed

homogeneous (the derivation of Eq. (1) assumes that the

medium is isotropic) so that the spatially dependent

probabilities V, ), and P in Eq. (1) are constant with

respect to r throughout each individual region.

Define a unit function, h, such that

h(p-s) 1 if p > s

= 0 ,if p s. (2)

For each interface between regions define an index k and

a position vector T such that -1k_ e JA. Then

h(rk-r) h(Ek-l-r) = 1, if Ek-r1 r Ek
= 0, if r e k-, or if r = rk

where r is on the line between rk and :k-l. One can

then use Eq. (3) to separate the r dependence of the

probabilities V, Y, and P in Eq. (1) in a special way.

In general, for any operator 0, one can write that

S(r,)- kl [h(Ek-r)-h(k-l-r)] (k,) (4)

Use Eq. (4) to rewrite Eq. (1) in the form

yvt* N(r,vn) + [ [h(rkE-r)-h(rk_,-r)]

S[v(k,v)+ (k,v)] N(r,vn)
8(r,vr) + 1 [h(rk-r)-h(!k-1-)]

dv'/ df' P(k,vI-v,n'-,O)N(r,v'n. ') (5)

Choose a particular region, say k', and remove
the term for the k' region from the summation on the LHS
of Eq. (5), meanwhile placing all other terms of that
summation on the RHS. Thus Eq. (5) becomes

vrL N(r,v*) + [h(Ek,-r)-h(rk'_-r1)

[V(k',v)+ /(k,v)] N(r,vn)

S(r,v) [(-r)-h( -r]

[V(k,v)+ y(k,v)] N(r,vn)

+ 1 [h(Ek-r)-h(k_-1r)]

dv' dn P(k,vi-v,n'~-N(r,v'n') (6)

Now using the identity

h(p-s) I h(s-p) (7)

one can write

h(!k,-r) h(Ek, -r) [1 h(r-k,)] h(k,-_l-r)

-1 [h(r-r,.)+h(r,_-r)] (8)

Substituting Eq. (8) into Eq. (6) and putting the
[h(r-r,) + h(Ek'~,i-)] term on the RHS gives

va*VN(r,vm) + [V(k',v)+'(k',v)] N(r,vn)

SS(rvn) + [h(r-,)+h(r .,-r)

[V(k',v)+ '(k',v)] N(r,vn)

[ [h(rkr)-h( ~.-"r) [V(k,v)+ '(k,v)] N(r,vn)

+ E [h(rki-r)-_h(!k-1L]
dv' dn' P(kv-'-,v,nt.) N(r,v'n')

Define a function gk' (r,vLr' ,vtn.') by requir-
ing that it satisfy the equation obtained by replacing
the RHS of Eq. (9) by the Dirac delta function product
6 (r-r') d (v-v') & (n-n'), subject to the same boundary
conditions as Eq. (9). Thus

TVgk T (rn r'Iv' 1)

+ [v(k',v)+ W/(k',v)] g' (r,vnlr' ,v'i')

g(r-r')&(v-v')s(n-n') (10)

The function gk (r,vnCr',v'.n') has a physical
significance; gk (r,vnr r',v'n') drdvdn is the number of
neutrons in the volume element dr about r with speeds
in dv about v and with directions of motion lying within
the solid angle dn about n due to a unit point source
at r' emitting one neutron per unit time in direction
n' with speed v' in a first-flight medium which has
the same properties as the medium in region k'.
It is important to understand fully the term
"first-flight medium" as used here. To be more explicit,
a first-flight medium is a medium in which all collisions

with nuclei result in removal of the neutron from the
population. Thus the removal probability for the Green's
function problem is equal to the total collision proba-
bility in medium k' of the real problem, and the neutrons
denoted by gkt(rvn l r*,v') are first-flight neutrons
in the sense that they have suffered no collisions
enroute from the source point to the field point. The
term first-flight neutrons will be used throughout the
following work.
The subscript k' will be dropped from

gk (r,vr'iE*,v't) in the following work, it being
understood that the function depends on the properties
of the particular region, k', chosen.
The function g(rvnjr',v'n') can be used as
a Green's function to convert Eq. (9) to the following
integral equation (see Chapter 7 of reference (16)).

N(r,vn) dr' dv' dn' g(r,vn r v'n) ,v'f)
rt v' n'

+ [h(r'-Ek)+hc(Ek.,-l-') [V(k',v')+ 2'(k',v')]N(r',v'n')

[h(!k-r' )-h(rkl-r')] [V(k' vt)+ Y'(k',v')] N(r' ,v')

+ [h(Ek-r)-h(k-r1]

dv"d P(kvyrv',n"-)N(r'vlg) .

Since Eq. (10) describes the neutron density
in a medium in which all collisions result in removal
of the neutron from the population, the only way a neutron
which is emitted at r' may arrive at r is by direct
flights, that is, with no collisions enroute. It will
be assumed and later justified that events at the bounda-
ries do not change the speed of the neutron. Hence if a
neutron of speed v' is emitted at r' and arrives at r,
its speed must still be v'.
From this discussion it can be concluded that
the speed dependence of g(r,vfir',v'n') includes a delta
function in speed so that

g(r,vn r' ,v'') g(v;r,~ r' ,n.' )(v-v'). (12)

The product solution, Eq. (12), is a result of the fact
that since there can be no scattering in the medium,
then there is no interdependence of nr*n' and v-v'. Note
that g(v;r,- lr',') ) still has an implicit dependence on
v due to the v dependence of V(k',v) and ((k',v) in
Eq. (10).

Substituting Eq. (12) into Eq. (11) and carry-
ing out the integration over v' yields

N(rvn) dr d g(v;r,nl r,f' )S(r,vn')

+[h(r'-rk )+h(r, l-rt)] [V(k' v)+ Z(k',v) N(r',vn)

Sh(Ek-r)-h(-r'i )][V(k,v)+ (k,v)] N(r',vn')

+ [h(rk-r')-h(k1-r')] (13)
k-1l --k -- --

dv" dn" P(k,v"--V ,n"-nr )N(r',v"n ) .

The integration over r' may now be considered.
The meaning of the functions [h(r'-rk,)+h(Ek1-r' ) and
[h(rk,-r')-h(rk_'-r') may be easily deduced from Eqs.
(2), (7) and (8) and is illustrated in Figure 1. The
function [h(r'-!k,)+h( k,l-r')] has unit value every-
where except in region k' where it has zero value, while
[h(k-r')-h(Ekl-r')] has zero value everywhere except
in region k where it has unit value. It can easily be
seen that the r' integration of the second term on the
RHS of Eq. (13) is an integration over all regions except


0 -

h(r-'k +-r' )

0 ----


h(r'-r) + h(rk--r')
0 -


h(r~kt-re )-h(rk' _1-r')
0 ------

Figure 1 The Spatial Dependence of the h Functions
of Equations (2), (7), and (8)

the k'th one while the r' integration of the remaining
step-function terms is a sum over k of r' integration
over each individual region k. Thus one can write
Eq. (13) as

N(r,,vn) dr' dn' g(v;r,, r',n')S(r',vn')
rt n

+ ki [v(k',v)+ (k',v)] V(k,v)+ -k,v)]

dr' dn' g(v;r,nrr',nt)N(r',vn')
r'in k r'

+ Z dr' dnY g(v;r,n jr,nt)
k-1 ~ -
r'in k n' (14)

dv" d"P(k,v"- v,n'L.' )N(r',v"r") .

Note that the k-k' term in the first summation
is automatically zero, so that no special notation is
needed on the summation symbol.



11.1 Expansion of Angular Dependence in Spherical
Harmonics and Zeroth Harmonic Approximation
The angular integration of Eq. (14) may be

performed as follows. Expand all the angular dependent

terms of Eq. (14) in the complete orthonormal set of

spherical harmonics as defined in Chapter 7 of reference

(15). Thus

oo n
N(r,vn)- Nm(r v)y-(n) (15)
n-0O m--n

co n
S(r,vn) E Sm(r,v)Ym(n) (16)
n-O m--n n n

oo a
P(k,v'"Uv,n'!n) Pa(k,vv)yb()b*(). (17)
a-0 b--a a- a

The expansion for the Green's function is more

complicated. Because of possible events at the bounda-

ries enclosing the space of interest, g(v;r,nI r' ,')

depends explicitly on both n and nr' and not just on r*n'

as P(k,vn.v,n'n ) does. For the moment only the n.' de-

pendence of g(v;r,n r',nh) is represented by an expansion

in the complete orthonormal set of complex conjugate
spherical harmonics. Thus

g(v;r,n r'l,n') gs;rn r')Y (). (18)
s-0 t--s

Substitute Eqs. (15) through (18) into Eq. (14) to obtain

N(r,v ndr' dn' gt (v;rnJf)Yt* (A) : m(r',v)Ym(n!)
Ss,t nm
r' n'

+ R v(k',v)+ (kt,v)0 [V(kv)+ Y(k,v)

Sfdr' dn E gg(v;r, lr')Y t*(n) N(r (n'
S s,t s n,-
r'in k n'

+ rl drt dn' g (v;r,nr')Y (')
r'in k n1'

ab a a n -

Equation (19) may be simplified considerably
by interchanging the order of integration and summations
and by using the orthogonality relation of the spherical
harmonics, namely,

n)Yb n, a Smb (20)

where Sn,a is the Kronecker delta function.
After carrying out the angular integration,
Eq. (19) becomes

N(r,vn) dr' E ge(v;r,nlr')Sm(r',v)

+ [V(k ,v)+ ](k ,v) [V(k,v)+ ?'(k,v)1

Sdr g gm(v;r,nr')Nm(r',v)
r'in k

+ R dr' dv" gn(v;rn l')Pn(k,vZ-v)o(r',v")
kl 1 f n,m
r'in k

Equation (21) is an expression for the total angular
neutron density and is exact; no approximation has been
made so far.
Now, the relative isotropy or anisotropy of
the n' dependence of g(v;r r',n') depends on the
boundary conditions applied to Eq. (10) and also on the

properties of region k'. Recalling the physical meaning
of g(v;r,nlr ,n.'), one observes' that n' denotes the

angular dependence of the unit point source in a medium

in which all collisions result in removal of the neutron

from the population. If the space of interest is of
infinite extent, then g(v;r,n r',n') has a double delta

function behavior in direction, namely S(-n*') -(-r')] ,

which merely describes the fact that since scatters

are not considered in the calculation of the Green's
function the direction of motion at r not only must be

the same as the direction of motion when emitted at r',

but also must be pointed along the vector connecting the

field and source points. In this case, the complete set

of coefficients g (v;r,,nir) would be required in Eq. (18)

to satisfactorily describe g(v;r,n r',n').

The problem is entirely different, however, if

one imposes perfectly reflecting boundaries around the

space of interest including the source as would be done
if one were calculating a unit cell in a repeating array.

Collisions or events at the boundaries will reflect

neutrons back into the space of interest with some

direction nL. If the neutron then arrives at r it will
have direction n, such that n + n' and n (r-r'),
although there llbe a definite relation between
although there will be a definite relation between

, n'_, and (r-r*). Further, if the medium characterized
by [V(k',v)+ ,(k',v)] has a mean free path that is long
relative to the distance between boundaries, then the
neutron may be reflected at more than one boundary before
arriving in the element of volume dr about r.
Now the integral of the angular density due to
a unit point source per solid angle, integrated over all
directions of emission n', is simply the angular density
due to an isotropic point source of one neutron per unit
solid angle or a total of 41T neutrons and is denoted by

g(v;r,nar'). Thus

g(v;r,nlr') g(v;r,nlr',n')dn'. (22)

From the discussion above regarding reflections at the
boundaries, it is obvious that in a finite space enclosed
by reflecting boundaries an isotropic point source will
produce a very nearly isotropic angular density

g(v;r, i r').
In order to take advantage of the relative
isotropy of g(v;r,n r'), one proceeds as follows. If
one asks only for information about the number density

N(r,v), the dependence of Eq. (21) may be integrated
out to obtain

N(r,v) dr' Z g(v;rlr )Sm rv)

k V (k ,v)+ k',v) [V(k,v)+ /(k,v)

dr' g(v;r r)lE (r',v)
r'in k
R m
+ dr' dv" n g(Vr'r)P (kv"-v)Nm(r',v")
k-1 n,mn n -
r'in k


N(r,v) -/N(r,v)dr (24)

Expanding the angular dependence of
g(v;rl r') in spherical harmonics gives

g(v;rn ') g(v;r r')Ym(n) (25)

The coefficients gn(v;rlr') in Eq. (23) are recognized
as being equal to the coefficients in the angular
expansion of g(v;r, nr'), a function which has been
discussed above and which has been concluded to be

nearly isotropic for a finite space enclosed by perfectly

reflecting boundaries and containing a medium having a

mean free path long relative to the distance between


Thus the coefficients gn(v;r r') form a rapidly

decreasing sequence for increasing n. Further, since

Eq. (1) includes a scattering term, the angular density

N(r,vn) should be rather isotropic, even for rather

strong absorbers. A unique advantage of the integral

equation method over the integro-differential equation

method thus becomes apparent; namely, the zeroth harmonic

approximation disregards only products of small terms,

that is, one assumes

-1 -1
1 1
0 0 0
g0 N0 g1 N *

g1 N1
1 1

Due to the orthogonality of the spherical harmonics,

there are no cross products to consider such as g0 No

or go N1 which might be of the same order of magnitude
0 0
as g0 NO'

The situation with regard to the scattering
term in Eq. (23) is even more favorable. One need not
assume completely isotropic scattering, but merely make
the less stringent assumption that the product of the
zero order coefficients of the Green's function, the
angular density, and the scattering kernel, is much
greater than the product of higher order coefficients,
that is,

g1 PI Ni1

g0 P NO 9 10 P NO

gl p N1
1 11

On the basis of the above discussion the
zeroth harmonic approximation is made in Eq. (23) and
only the n-0 terms are kept. Thus Eq. (23) becomes

N(r,v) dr' g0(v;r r')So(r',v)


+ [v,(k.v)+/(kv) [V~f,v)+(kv

dr' g00V;r Ir)N (r',v)

r'in k

R O l o t
+ Z dr' dv" go(v;r r')P.(k,v.v)NO(r,v"). (26)
+-1 kl d -
r'in k

The zero order coefficients may be put in a
more recognizable form as follows. Substitute Eq. (15)
into Eq. (24) to obtain

N(r,v) Na(r,v) Ya(n)dn

Nm(r,v) Ym(c) -/f YO*()

-/41 N(r,v)

where use has been made of the fact that Y* (n)
- YO n) 1/--1 Thus

N0(rv) N(r,v) (27)

Similarly, defining

S(r,v) -s (r,vn)dc (28)

and using Eq. (16), one obtains

o 1
S(r,v) S(r,v) (29)
0 --


g(v;rlr) g(v;r,lr')dn (30)

to be the number density of first-flight neutrons at r
having speed v due to a point source at r' emitting 4TT
neutrons isotropically with speed v in a homogeneous
finite medium. Substitution of the angular expansion
Eq. (25) into Eq. (30) then gives the result that

g0(v;r ') 1 g(v;r r') (31)

For the scattering kernel one has from Eq. (17)

P(k,v-.v, -)dn .- Z Pa(k,vi.v)Y b Y b d)
b 0

abb* /yb ,d

P(kv"-.v) ,/4f YO* ()

PO(k,v!',v) (32)

which is just the probability per unit time per unit
speed v that a neutron having speed v" will suffer a
scattering collision in region k and emerge from the
collision with speed v. The zero subscript will be
dropped and the function will subsequently be denoted
simply by P(k,v!3.v).
Finally, substituting Eqs. (27), (29), (31),
and (32) into Eq. (26) gives

N(r,v) dr' g(v;r r')S(r',v)

+ 4 1 [(V(ktv)+ (k',v)] [V(k,v)+ Y(kv)]

Sf dr/ g(v;rlr)N(r',v)
r'in k

k+ kl f dr' g(v;r r') dv"P(k,v'4v)N(r ,v").
r'in k (33)

11.2 Advantages and Disadvantages of Homogeneous
Green's Function Method
The preceding method of converting the integro-

differential transport equation to an integral equation

with a kernel defined by Eq. (10) has at least two

advantages over the more standard methods of solution.

First, as discussed above, the zeroth harmonic approxi-

mation involves discarding only products of small terms

as in Eqs. (22) and (23). For example, rather than neglect-

ing the higher moments of the neutron density, as must

be done when the integro-differential form of the transport

equation is used, one need only neglect the products of

the higher moments of the Green's function times the higher

moments of the neutron density.

A second advantage and a unique feature of the

Green's function method described in this work is that the

Green's function defined by Eq. (10) depends only on the

dimensions and geometry of the outer boundaries of the real

problem for which the neutron density is to be calculated,

the total cross section of the medium in the k'th region

of that problem, and the boundary conditions imposed on

Eq. (10). Once the Green's function has been determined,

it may be used to determine the neutron density in any

heterogeneous medium, the only restrictions being that

the outer boundaries of the neutron density problem must

coincide with those of the Green's function problem,

the boundary conditions must be the same for both g

and N, and the total cross section of the k'-th region

of the neutron density problem must be the same as

that for the homogeneous medium of the Green's function


A disadvantage of the method is implied above;

namely, if the outer boundaries and/or boundary conditions

of the neutron density problem do not agree with those for

which the Green's function has been determined, then a

new Green's function must be determined. Another dis-

advantage is that, as yet, the method has not been

extended to higher order approximations such as keeping

gm N P, g Pn N terms for n 0, although such an

extension to higher moments would, in general, follow

the work of Wilf (29), and would not alter the above

mentioned advantages.



III.1 General Speed Dependent Equation

Equation (33) is a general speed-dependent

equation and is valid over the entire speed spectrum,

the limits of the speed integration being from v" 0

to v" -oo. As indicated in the introduction, the

present investigation will be confined to speeds in

the thermal region. In practice, there exist no true

sources for the thermal energy region other than slowing-

down sources, and hence S(r',v) 0 in Eq. (33). The

speed integration in Eq. (33) are then separated into

two parts, one for 0 < v i v* and one for v =>v*. The

speed v* represents the upper limit of the speed region

of interest, and is characterized by the following

assumptions (17):

1. For v = v*, the neutron density speed distri-
bution is proportional to 1/v2 (that is, the
neutron flux is proportional to 1/E).

2. The probability of a neutron being scattered
from a speed below v* to a speed above v* is

3. The motion of the collision atoms is ignored
with respect to that of the neutrons with'
speeds above v*.

4. The Maxwellian distribution for thermal
neutrons or collision atoms has negligible
magnitude for speeds above v*.
5. For neutron energies above v*, the energy
transfer is computed assuming that the col-
lision atoms are free from chemical bonds.
This implies that v* should correspond to a
neutron speed well above that required to
excite vibrations in the collision molecules,
so that the collision atoms would recoil
freely at the instant of impact.

6. For neutron energies below v*, the energy
transfer kernel P(k,vv-v) will depend on the
particular scattering model chosen. A
discussion of particular energy transfer
kernels and their scattering models is given
in references (17, 31).

The S(r',v) term in Eq. (33) is then modified
to include all neutrons scattered into a speed interval

dv about v from speeds above v*. Thus Eq. (33) may be

rewritten as

N(r,v) drI g(v;r r')S(r',v)

r'in k

+ 4 ( [V(kv)+ /'(k',v)] [V(kv)+ Y(kv)]

dr' g(v;r r)N(r',v)

r'in k


+ k dr' g(v;r r') dv" P(k,v"!.v)N(r',v")

r'in k v"-0 (34)


S(r',v) dv" P(k,v!.v)N(r',v") (35)


III.2 The Source Term

Assumption (5) of III.1 implies that ordinary

slowing down theory may be used to derive the scattering

kernel for v v* (18, 19). Define

--- (36)

Then ov is the minimum speed that a neutron of initial

speed v can have after collision with a target nucleus

of free atom mass m; the energy transfer probability,

that is, the probability that a neutron of energy E'

will after scattering have energy in the interval dE

about E is (see (18))

T"(r,E..E)dE (dE for a2E' E E'

-0 ,for E < a E'

S0 for E E' (37)

or, changing variables,

7l(r,vi.v)dv 7T(rE'-,E)dE dv

so that

7rf(r,v.v) 7r(r,EL-.E) dE

1 dE mv
SE' [la2(r)] dv l-o2(r)] (mv'
["I, vv for wav' v S v'

S0 for v Oav'

0 for v >v' (38)

Now, the differential scattering cross section,
Lg(r,vi-av) dv, that is, the probability that a neutron
of speed v' will be scattered at r into dv about v is
equal to the product of the total scattering cross
section at r for neutrons of speed v' times the proba-
bility that the speed of the neutron after collision
will be in dv about v. Thus

Zg(r,v4..v)dv s(r,v')7r(r,v'.-v)dv. (39)

Then the probability per unit time that a neutron of
speed v* will suffer a scattering collision at r such
that it emerges from the collision with speed in dv
about v is simply

P(r,v'.v)dv v'1%(r,v'Lv)dv (40)

so that for v' > v* one has from Eqs. (38), (39) and
(40) that

P(r,v-.v) 2vv's(r,v)
[1- (r ) v'2

for Ov' s v s v'
-[l-a(r)] v'

0, for v ctv'

0, for v => v

Assumption (1) of III.1 implies

N(r,v') A for v' > v* (42)

where A(r) is an arbitrary constant with respect to
speed and is normalized to give a stationary distribution
for N(r,v'). Substituting Eqs. (41) and (42) into
Eq. (35) yields

S(r',v) dv" 2vls(k,v") A(r')
S1-a2(k) v" v"2
v"'v* [I I

2vs(k)A(r,) v/(k)
-n [ldv" (43)
[1-a2(k)] j v"3 d

where Zs has been assumed to be independent of v in the
speed range v* 5 v" v v/a(k).

Carrying out the integration in Eq. (43) gives

8(r',v) V A(r') v2 (k) if ca(k)v* v s v.
[1-o(k)] V*2 v2

0 if v < a(k)v* (44)

where the constant A(r) determines the source strength;
A(r) will be assumed to equal a spatially independent
constant A, that is, one assumes a spatially flat
slowing-down source distribution. As Honeck (22) has
pointed out, A(r) could, if necessary, be approximated
by the shape of the resonance group neutron density as
predicted by diffusion theory and the use of suitable
resonance group cross sections given in Weinberg and
Wigner (19).

II.3 One-Speed Equation

To obtain the one-speed equation, one must
assume separability of the spatial and speed dependence
of the neutron density. Then integrating Eq. (33)
over all speeds and using Maxwellian averaged cross
sections, one has

N(r) /dr'g(r r')S(r')


+ 1- [v(k')+ $(k')] -[V(k)+ Y(k)) dr' g(r r')N(r')

r'in k

+ 1 V(-k) dr g(r r)N(r (45)
r'in k

where use has been made of the fact that

P(k,vu.v)dv V(k,v") (46)

which is just the total probability per unit time that
a neutron in region k and having speed v" will suffer a
scattering collision.

In the one-speed case motion of the collision
atoms may be ignored, so that in conventional notation
Eq. (45) may be written as

N(r) 1 dr' g(r r')S(r')

+ 1 (k')- (k)] dr' g(r|r')N(r')
r'in k

r'in k
kj g(rfrl')N(r ) (47)

Finally, the last two terms in Eq. (47) may be combined
to obtain

N(r) dr' g(r r')S(r')

+ k t(k')-a(k) dr' g(rr')N(r'). (48)

r'in k



IV.1 Development of Difference Density Equations

An alternative to the preceding method of

describing the real neutron density, N, in a hetero-

geneous medium is to define a difference density, D,

as the deviation from a known density, N', for a similar

or reference space. Although it may appear at the outset

that such a method is analytically identical to that

developed above for the real density N, nevertheless,

a number of authors have found such a procedure to be

advantageous (13) and, in fact, necessary (12) in order

to obtain results for specific infinite geometry problems.

It will be of considerable interest, therefore,

to develop equations for the difference density, D,

parallel to those above for the real density N. A

direct comparison of the two methods will then be made

in order to determine their relative advantages when

applied to the specific problem to be considered in this


Thus the difference density D(r,vn) is defined

such that

D(r,vn) N(r,vri) N'(r,vn) (49)

N(r,vn) = D(r,vn) + N"(r,vn) (50)

Using primed symbols to denote the reference
space, one can write the transport equation for the
assumed known angular density as

v ~-IN'(r,v) + V'(r,v) + (rv) N'(r,Cvn)
.L (51)

S'(r,vn) + dv' dn P'(r,v-.v,n4 n)N' (r,v'n').

Substituting Eq. (50) into Eq. (1), subtracting Eq. (51)
from the resulting equation, and grouping all of the
N'(r,vn) terms on the RHS gives

vnVD(r,vn) + [V(r,v)+ (r,v)] D(rvi)

= S(r,vn) S'(r,vn)

v(r_,v)+ (r,v)-[,v) (v)+ /(r,v)] N'(r,vr)

+ dv dn' P(r,v--v,n'-.*-0D(r,v'n') (52)

+ dv' d.n P(r,vv,nI)-P' (r_,vv,_)] ,n N(rv,',).

As before, Eqs. (3), (4), and (8) are used to separate
the spatial dependence of the terms in Eq. (52). Thus


vn-VD(r,vn) + [v(k',v)+ (k',v)] D(r,vn)

- S(r,vn) S'(r,vn)

+ [h(r-rk,)+h(k'l-'r)] [V(k',v)+ K(k',v)] D(r,vn)

rh(k-r)-h(rkl-r)] [V(k,v)+ K(k,v)] D(r,vn)

- k 1 h(Ek-r)(k-1-^r)] [V(k.v)+ Y(k,v)]

[V'(kv)+ 't(k,v)] N'(r,vr)

+ r h(rk-r)-h(kl,-r)] fdv dn' P(k,vi-*v ,nn)D(r,v'nW)
k=1 f

+ h(k-r)-h('k-1-r) dv' di' [P(k,vI+v n!*n)
-P' (k,vv,L).)] N' (r,v'n').

The boundary conditions imposed on Eqs. (51) and (53)
are the same as those for Eqs. (1) and (9). Thus the
Green's function which satisfies Eq. (10) may be used
to convert Eq. (53) into an integral equation as before.
Then, using Eq. (12), one may perform the integration
over v' to obtain

D(r,vn.) dr' d' g(v;r,g r' n') S(r' ,v')-S' (r' ,vn )

+ [h(r'- ,)+h(k, -rl') [v(k,v)+ ) k',v)] D(r',vn')

- [hr ktr')-h(h-1-r')1 [V(k,v)+ Y(k,v)]D(r',vn' )

[h(Ek-r')-h(ik-lr')] [V(k,v)+ Y(kv)

yv' (k,v)+ 1 (k,v)] N' (r' ,v')

+ kl[h(Ek-r')-h(rek_-r') dv" dr' P(k,v" vn'n*')D(r',v'tn")

+ k E [h(Ek-r' )-h(Ek-i-r') dvt" d [P(kv' mt-.')
i (54)
P' (k,v-vn'".)] N' (r' vn") .

The r' integration for Eq. (54) is identical
to that for Eq. (13) and the result is

D(rv_) dr' dn' g(v;r,n r' ,n)[S(r',v f)-S'(r_',vn)]

+ V(k',v)+ (k',v)] [V(k,v)+ Yk,v)

dr' da g(v;rn ir'L l )D(r' ,v')
r'1in k 1'

+ dr' dn' g(v;r,nlr'l,.')
r'in k n

dv" a" P(k,v" vn'1,-c n')D(r',v"n')

- [V(k,v)+ Y(k,v)] [V'(k,v)+ (k,v)

S dr' dn' g(v;rnl r',n')N'(r',vn!')
rrin k '

+ dr, dn' g(v;r,nlr',n') dv"fdrl" P(kv4I.vI.nIn')
r'in k T.'
P' (k,v"-.v,n "-.')] N' (r' ,vI'") .

The angular integration presents no new
problems and may be performed as in Chapter II. Making
the zeroth harmonic approximation and asking only for
the scalar neutron difference density, one obtains

D(rv) d 1 dri g(v;rjr9)[ (r',v)-S'(r',v)]

+ ([V(k',v)+ (k',v)]
47r k-i
[V(k,v)+ Y(k,v)]) dr' g(v;r r')D(r' ,v)
r'in k

k dr' g(v;rlr') dv" P(k,vU*v)D(r',v")
k-lr -
r'in k


+ 1



Sdr' g(v;rr')N'(r',v)
r'in k

dr' g(v; rr')
r'in k


*dv" P(k,vt v)-P' (k,v"-v) N' (r' ,v").

Writing Eq. (56) specifically for the thermal
spectrum yields

D(rv) dr' g(v;rr')[S(r'v)-S'(r',v)]
r'in k


[V(k',v)+ (k',v)]

[V(k,v)+ /(k,v)] dr' g(v;r r')D(r',v)
r'in k

S[vk,+ (k,v) )]

- [V'(k,v)+ Y'I(kv))



v f

1 R
4+ F i dr' g(v;r r') dv" P(k,v"!v)D(r',v")
r'in k v"'-0

-1 R i (k,v)+ (kv)]

[V'(k,v)+ '(kv dr' g(v;r r')N'(r,v)
r'in k

+ k~ dr' g(v;rlr')
r'in k


dv" [P(k,v,-.u.v)-P' (k,v"--.v)] N' (r' ,v")


where, in region k,

S(r',v) S'(r',v) dv" P(k,v-!.v)D(r',v")

00 (58)
+ dv" [P(k,~,!v)-P'(k,v".v)] N'(r',v")

and the scattering kernels for v > v* are determined
from ordinary slowing down theory as in 111.2.
Finally, the one speed version of Eq. (56) is

D(r) dr' g(r r')S(')- (

1 R

+ -- k

+ Rl

1 R

+1 R
47( N-1

v[zt(k')- It(k)] dr' g(r r')D(r')
r'in k

v Zg(k) dr' g(r r')D(r')
r'in k

v [zt(k)- Z((k) drI g(r rf)N'(r')
r'in k

v[(k)- (k) dr' g(r r')N'(r')
r'in k

or if various terms are combined,

D(r) /dr' g(r r')[S(r')-S'(r')]


+ v k')- Zak)] dr' g(r r)(r
r'in k

S1 v[Za(k)- (k)] dr' g(rr')N'(r') (60)
47r T l k f i
r'in k

IV.2 Equivalence of Equations for N and D
In order to establish the equivalence between
Eqs. (34) and (57) and between Eqs. (48) and (60),
the one-speed equations will be considered first.
Subtracting Eq. (60) from Eq. (48) and using
Eq. (50) gives

N'(r) 1 dr' g(r r')S'(r')

+ k V[Zt(k*)- a(k)] dr' g(r )N(r')
r'in k

+ 1 v[Za(k)- Z(k)] f drv g(rlr')N'(rt)
k 47j
r'in k

_- fdr' g(r r')S'(r')
47r f

+ RF v[ l(k)- .(k)] drt g(rr)N(r'). (61)

r'in k

Comparing Eqs. (61) and (48), one observes that the
requirement that Eqs. (60) and (48) be equivalent is
simply that It(k') -= (k'), that is, the total cross
section for the k'th region of the reference space be
identical to the total cross section of the k'th region
of the real space of interest. This simply says that
the Green's functions in Eqs. (48) and (60) must be
identical, that is, the Green's functions in Eqs, (48)
and (60) must both satisfy Eq. (10) and its boundary
conditions. But this is just the condition imposed in
the original derivation, and hence Zt(k') is equal to

Z (k'); equations (48) and (60) are therefore proved
to be equivalent formulations of the same problem.

Similarly, subtracting Eq. (57) from Eq. (34)
and using Eq. (50) gives

N'(r,v) = dr' g(v;rr')' (r ,)
*E k-lf ~ # ~
r'in k

1 R
47T kl

[V(k,v)+ a(k' ,v)]

Iv(k,v)+ Y(k,v)] dr' g(v;r r')Nt(r',v)
r'in k

+ k dr' g(v;r r') dv" P(k,v"..v)N (r' ,v")
r'in k vIO


k [V(k,v)+ ~(kv)

[V'(k,v)+ f'(k,v)]

r'in g(v;rlr)N'(r',v)
r'in k

Sdr' g(v;rlr')
r'in k

* dv" [P(k,v-t..v)-P' (k,v"!.v) N' (r' ,v")
v" 0





r'in k

dr' g(v;r r')8'(rt,v)

+ 1 [V(k'iv)+ K(k',v)]

['(k,v)+ (k.v)) dr' g(v;rjr')N'(r',v)
r'in k

+ 1 dr' g(v;r r') dv" P'(k,,v'la)N'(r'v,v")
r'in k V1=0

Comparing Eqs. (34) and (62) one again observes
that, just as in the one-speed equation, the requirement
that Eqs. (34) and (57) be equivalent is simply that
[v(k',v)+ (k' ,v)] be equal to [V'(k',v)+ ('(k',v)],
that is, the medium of the k'-th region of the reference
space must be the same medium that is in the k'-th
region of the real space of interest. As has been
mentioned above, this condition was imposed in the
original derivation so that the equivalence of Eqs. (34)
and (57) is established. A discussion of the relative
usefulness of the two formulations will be deferred
until the general equations have been specialized to

slab geometry.



V.1 One-Dimensional Problem
The equations of the previous chapters were
written for a general three-dimensional geometry, and
the Green's function g(v;r r') was interpreted as the

number of first-flight neutrons per unit element of
volume at r having speeds v due to an isotropic point
source at r' emitting 47r neutrons per unit time having
speeds v in a homogeneous medium. Before specializing
Eqs. (34) and (57) to an infinite slab geometry it will
be convenient to define a new Green's function G(v;r r'),
such that

G(v;r g(v;rr') (63)

Thus G(v;r r') is the number density of first-flight
neutrons at r having speed v due to a unit isotropic
point source at r' emitting neutrons of speed v in a
homogeneous medium. Rewriting Eq. (34) in terms of

G(v;rlr') gives

r'in k

+ 1 fV(k',v)+ /(k: ,v) -[V(kv)+ /(k,v)]

*/ dr' G(v;rir')N(r' ,v)
r'in k

v* (64)
+ 17 dr' G(v;r r') dv" P(k,vlvm)N(r',v")
k-.f -I- f
r'in k v"0O

and similarly for Eq. (57).
The transformation to plane geometry may be
performed as follows. Referring to Figure 2, the element
of volume dr' may be expressed as

dr' 2 7Ta da dx'

R r r'

R r r'I (65)

From Figure 2, and Eq. (65) one obtains

R ( x-x'| 2 + a2)1/2 .x-_
In I

r fixed field point
r' variable source point
Y cos G

z /

//'r O
/ X

r x
/ L._



Figure 2 Plane Geometry Transformation

so that for constant Ix-x'l

dR = a da -Ix-x'
a_2 dr .
R 2

dr' 2 7 R dR dx'

_Hx-x'lX 2 d_
S2[Ix-x'I -dx' (66)

where the i integration is from -1l to y-0, and the R
integration is from R jx-x'I to R -oo .
Observing that in plane geometry the r de-
pendence of the functions N, N', D, S, and P is simply
a spatial dependence on x, one may write Eq. (64) for
plane geometry as

N(x,v) 1 dx' S(x',v) G(v;r r') 2 7r R dR
x'in k R- x-x'I

+ 1 [v(k v)+ '(k',v)] (k,v)+ (kv)]

Sdx' N(x',v) G(v;r r') 277R dR
x'in k R-|x-x'l

v* 00
+ 1 dxt dv" P(k,vl-.v)N(x',v") G(v;r r') 27TR dR
xin k v"-O R-Ix-x'i (67)


G(v;xx') G(v;rr ') 2rrR dR (68)
R- x-x'l

The function G(v;x x') is the density of first-flight
neutrons at x having speeds v due to a unit isotropic
plane sourceat x' emitting one neutron per unit time
with speed v in a homogeneous medium. This function is
well known for an infinite medium (see Chapter 2 of
Case, et al. (21)) and will be discussed extensively
for a finite unit cell in Chapter VII.
Using Eq. (68) to rewrite Eq. (67) gives

N(x,v) 4 /dx' G(v;x ,x)S(x,,v)
x'in k

+ k1 [V(k',v)+ /(k',v)]-[V(kv)+ Y(k,v)]

dx, G(v;x x')N(x',v)
x'in k


+ dx' G(v;x x')
x'in k v"0

Similarly, Eq. (57) for plane geometry becomes

D(xv) / dx' G(v;x x')B(x'.v)-S'(x'.v)]
x'in k

+ k1 [(k',v)+ (k',v)]-[V(kv)+ (k,v)]

/dx' G(v;x x')D(x',v)
x'in k

+ f dx' G(v;x x')
x'in k v"-0


3v" P(k,v-L.v)D(x',v")

[V(k,v)+ Y(k,v)] [V (k,v)+ YO (k,v)]

dx' G(v;x x')N'(x',v)
x'in k


+ dx' G(v;x x') dv" P(k,v'-wv)-P'(k,vtlv) N' (x',v").
x'in k v-0

dv" P(k,v"L.v)N(x',v") (69)

V.2 Symmetry Considerations
For subsequent work the space of interest
will be assumed to be a plane geometry symmetrical
about the plane x-0. In the summations over k as
previously discussed, k-l denoted the left-most region
in space. The regions will now be counted differently,
letting k-l denote the region adjacent to, and to the
right of, the x-0 plane. Thus a summation over k from
k-l to k-R is replaced by a summation over k from
k -R/2 to k R/2. Let

R2 "- (71)

denote the number of regions in the half-space. For a
general function f(x) which is symmetric about x-0,
that is, f(x) f(-x), one can then write that

E f(k)J dx' G(v;x x')N(x',v)
x'in k

R f(k) dx' G(v;xlx')N(x',v)
k-R f(k)
x'in k

+ ( f) dx' G(v;x x')N(x',v)
k--R2 k- f(k)
x'in k

0 b

= dx' + dx f(x')G(v;xlx')N(x',v) (72)
x'--b x'-O

where b and -b denote the right and left boundaries of
the space of interest. The notation used is defined

such that

cm c2 7C
+f ... + + k (73)
kl kn k-c1 k-en

C2 cm c2 cm
dx' + ... + dx' f(x') dx'f(x') + ...+ dx'f(x').
x'-c1 X' cn x'c 1 X'-cn

In the first integral in Eq. (72) let x'--y,

dx'=-dy. Observing that f(x)-f(-x) and N(x,v)-N(-x,v),
one can write

0 0

dx'f(x')G(v;x x')N(x',v) dy f(-y)G(v;x -y)N(-y,v)
x' -b y-b


/dy f(y)G(v;xl-y)N(y,v) (75)

Then let y-x' so that Eq. (75) becomes

0 b
Sdx'f(x')G(v;xx')N(x',v) dx'f(x')G(v'x -x')N(x',v).
x-b -0 (76)

Substituting Eq. (76) into Eq. (72) yields

Z f(k) dx'G(v;x x')N(x',v)
x'in k

dx'f(x') [G(v;x x')+G(v;x -x')]N(x',V)
f (k) dx' [G(v;x x')+G(v;x -x')]N(x',v)
x'in k

As noted above, G(v;x x') is the scalar
density of first-flight neutrons at x due to an isotropic
unit plane source at x'. If one defines

H(v;x x') G(v;x x') + G(v;x -x') (78)

to be the half-space Green's function, that is, the
scalar density of first-flight neutrons at x due to
isotropic unit plane sources at both x' and -x', then

equation (69) may be rewritten for the half-space in
the form

N(x,v) dx' H(v;x x')S(x',v)
x'in k

+ kV (V(k',v)+ Y(k',v) -fV(k,v)+ 7(k,v)1

f dx' H(v;x x')N(x',v)
x in k

v* (79)
+ dx H(v;xlx') dv" P(k,v'.v)N(x',v")
x'in k v"-0

while Eq. (70) becomes

D(x,v) R 2dx' H(v;x x') S(x',v)-S'(x',v)]
x'in k

+ 2 ([V(k'.v)+ (k',v)]-[V(kv)+ (kv)]

dx' H(v}x x')D(x',v)
x'in k

+ T dx' H(v;xlx') dv" P(k,v!".v)D(x',v")
x'in k v"-0

- ([V(k,v)+ /(k,v)]-[V'(k,v)+ Y'(k,v)]

dx' H(v;xIx')N'(x',v)
x'in k
+ f dx' H(v;x x') dv" [P(k,v!.Lv)-P'(k,v!..v) N'(x'f,v")
x'in k v"10

The one-speed equations may also be written
for the one-dimensional half-space. Thus Eq. (48)

N(x) dx' H(x x,)S(x')
x >0

+ :2 v[Zt(k')- Ia(k)] fdx' H(x x')N(x') (81)
x'in k

while Eq. (60) becomes

D(x) .dx H(x xI) [(x')-BS(xv)
x') O


+ E2 v[t(k')- Za(k)] dx' H(xx')D(x')
x'in k

i2 v (k)-a( j) dx' H(xJx')N'(x') (82)
x'in k



For the direct comparison of the formulation
in terms of real density,N, with the formulation in

terms of difference density, D, only the one-speed

equations will be considered, namely Eqs. (81) and (82).

Similar arguments would apply to the multispeed equations.

Suppose that all R regions of the reference
space contain the same medium and that the reference

space itself constitutes one spatial period of an

infinitely repeating lattice, that is, the reference

space is a unit cell. Then the neutron density N'(x)
in the reference space is that for an infinite homo-

geneous medium and may be assumed spatially flat.

Since the neutron density and medium are both
spatially uniform in the reference cell, one may assume
that the statement of neutron conservation holds not

only over the entire reference space, but also over any

arbitrarily small element of volume in that space. Thus

at any point

B'(x) S' v IN'. (83)

Substituting Eq. (83) into Eq. (82) and
cancelling terms where possible, one obtains

D(x) -dx H(xlx')S(x')

+ a vt(k')- (k)] dx' H(xjx')D(x')
x'in k

v a(k)N' dx' H(x x') (84)
x'in k

Then using Eq. (50), one may combine terms to write

D(x) dx, H(x x')S(x')

+ k v 2(k') f dx' H(x x')D(x')

v Z.a(k) dx' H(x lx)N(x') (85)
k-1 f
x'in k

Adding and subtracting the quantity

v t(k') fdx' H(x x')N(x')
k-1in k
x'in k

to Eq. (85) gives

D(x) dx' H(x x)S(x')
x 'O

+ 1k2 v[Zt(kt)- Za(k) f dx' H(x x')N(x')
x'in k

I v lt(k') dx H(x x')N'
x'in k

dxt H(x x')S(x')

+ 2 v[t(k')- Za(k)] dx' H(x ')N(x')
x'in k

v I(k')N' dx' H(x x') (86)
x '0

But from Eq. (10)

v Zt(k')H(x x')dx' 1 (87)
x 50

so that using Eqs. (50) and (87), one can reduce Eq. (86)

directly to Eq. (81).

The direct reduction of Eq. (82) to Eq. (81)
having been established, the question remains whether
the use of Eq. (82) has any advantages over the use of
Eq. (81). It will be of particular interest to consider
whether or not the complexity of Eq. (82) can be reduced
by relating the real source S to the source 8' in the
reference space. A few such special cases are considered
below. In each case the assumption is made, as it was
above, that all R regions of the reference space contain
the same medium and that 8' and N' are spatially flat.
Case 1 S(x) 8' Over All Space
Equation (82) then becomes

D(x) 2 v I(k')- a(k)] dx' H(xJx')D(x')
x'in k

va(k)- i(k)]N f dx' H(xxt)

R- 2 tk'- a(k)] dxz H(x x')D(x')
x'in k

v Za(k)N' fdx H(x x') + vE(k'
x'in k

where Eqs. (83) and (87) have been used to obtain the
last term of Eq. (88).

Since S(x) = S' implies that S(x) is spatially
independent, Eq. (87) may be used to show that the first
term in Eq. (81) becomes

dx' H(x xt')(xt) = 8 dx' H(x x')
x'=O x '=O

v St(k') (89)

Thus Eq. (88) can be reduced directly to Eq. (81) in
the same manner that Eq. (82) was reduced to Eq. (81),
so that there is no advantage of the difference density

formulation over that for the real density N. The
reason for this is that, although it is given that

S(x) S S', one must still calculate N' based on 8',
and the N' terms then act as source terms in the differ-
ence density equations.
Case 2 S(x') 0 For Region k-1, S(x') 8'
For Regions k 1

Equation (82) becomes

D(x) dx' H(xIx')S'(x')
x'in 1

+ 2

- 2

v[ZCk)- (k)]

dx' H(xIx')D(x')
x'in k

Sdx' H(x x')N'(x)
:'in k

To specialize further, let R2 2 and k' 2.

Eq. (90) becomes

D(x) - S' dx'
x'in 1

+ v[Zt(2)-~a(1)

H(x x')

dx' H(x x')D(x')
x'in 1

+ v Z~(2) dx' H(x x')D(x')
x'in 2

v[a(1)- a(1)] N' dx' H(xl x').
x'in 1

Using Eq. (83) one may write Eq. (91) as

D(x) v[Zt(2), a(l)] idx' H(xI xtD(x')
x'in 1




+ v zS(2) dx' H(x Ix)D(x')
x'in 2

a(l) N' dxv H(xlx') (92)
x'in 1

or, again substituting Eq. (83) into the last term of
Eq. (92), one finally obtains

D(x) v[Zt(2)- Za()] f dx' H(xx')D(x')
x'in 1

+ v ZS(2) dx' H(x x')D(x')
x'in 2

(1a(') 8'
7 a(2l) dx' H(x x') (93)
x'in 1

For this special case, Eq. (81) becomes

N(x) SI dx' H(xlx')
x'in 2

+ v[t(2)-Za(1)] /dx' H(xlx)N(xt)
x'in 1

+ v Zs(2) dx' H(xlx')N(x') (94)
x'in 2

The obvious difference between Eqs. (93) and
(94) is the source term. In a typical application of
calculating the scalar neutron density in a space
containing a gold foil in a graphite medium, the ratio

Za(1)/Za(2) would become

(1) Au 5.79
1 _7 1.94 x 104
a a
a(2) ^C .000299

It can be seen that in a finite medium containing a
highly absorbing region, even though it be of small
spatial extent, the neutron density will not necessarily
be a small deviation from a known density in a reference
space with no absorber. In fact, the deviation in
magnitude of the densities will be quite large, even
though the spatial shape of the neutron density N may
not be vastly different from the spatial shape of the
known density N'.
It seems obvious now that for problems involv-
ing finite geometries the analytic method of solving for
a difference density, D, offers no advantages over


the direct solution for the unknown density N; the

difference density method may, indeed, be limited in

usefulness to the solution of infinite geometry problems.



VII.1 Differential Equation and Boundary Conditions
The first-flight Green's function
6k' (r,v r',v'n') for a general three-dimensional homo-
geneous medium is defined such that it satisfies Eq. (10),

va* Vgk' (r, vn r', v'n' )

+ [V(k',v)+ (k' ,v)] gk'(rvnr' ,v'n')

&(r-r )S(v-v )S(cn-n'); (95)

further, gk (r,vnr',v'_1') is subject to the same
boundary conditions that are applied to the angular
density N(r,vn).
Considering the RHS of Eq. (9) as a fictitious
source Q(r,vn), one may rewrite Eq. (9) as

vn.7VN(r,v) + [V(k',v)+ ((k',v)]N(r,v)) Q(r,vn) (96)

Thus N(r,vn) represents the first-flight angular density
in a medium as discussed in Chapter I, the removal pro-
bability being equal to the total probability of col-
lision in region k', with a source of neutrons Q(r,vn).
If the homogeneous medium k' is of infinite
extent, then the solution to the one-speed version of
Eq. (95) is known to be (see Chapter 2, of Case, et al.

t(k') -d
gi (rlr S 2 n-r!' )'(n f-r') (97)

where the superscript oo on g signifies that Eq. (97) is
valid only for an infinite medium or for a finite medium
enclosed by a non-re-entrant surface with vacuum boundary
conditions, and r-r' denotes the unit vector (r-r')/jr-r'
The extension to the speed dependent problem is obvious
from the discussion of Chapter I and Eq. (12). Thus

Zt(k',v) r-r'
S0, (r ', ,v'n') e

o(v-v')&(n-n')An ;Cc') .

Integrating Eq. (98) over v' and n! gives

gg?(v;r,n|r') e t-----(k ,v ); (99)
v r-r (99)

integrating over n then yields

k (v;r r') e -Ztk v) (-0oo
v r-r 2 (100)

where use has been made of Eq. (30).
The simple form of Eq.' (100) is due to the
fact that the medium k' is homogeneous so that the
collision probability is not a function of the flight
path of the neutron. As pointed out by Case, et al.
(21), the case of a spatially variable collision proba-
bility is a considerable complication because absolute
coordinates instead of only relative coordinates between
the field and source points enter the problem. It might
be well to point out again that a unique advantage of
the homogeneous Green's function method being described
is that it completely avoids the complication of a
spatially variable collision probability in the evalu-
ation of the Green's function.

From Eqs. (63) and (100) one can write the
expression for G(v;r r') for an infinite medium. Thus

(k';v _r-r'1
G- (v;rr') e(101)
47 v r-Ir' 2

From Eq. (68) one can then evaluate G(v;x x') for an
infinite medium. Thus, using the transformations of
V.1, one obtains

,3- ,,v Ir,-rIt(k'
Gco (v;x x') 27r R dR
47T v r1-r' 2

S e- t(k',v)R -1 dR (102)
R-I x-x'

The exponential integrals or En functions (21, 24) are
defined such that

En(x) e-xu u-n du (103a)

-i pn-2 e-X/P dp (103b)


Sxn-1 e-u u-n du (103c)


so that Eq. (102) may be written

Goo (v;x x') L El [Zt(kv)Xx-x] (104)

The Green's functions discussed above were,
in general, for an infinite medium. The problem
eventually to be considered, however, is that of a unit
cell; that is, the space of interest represents one
unit volume of an infinitely repeating set of unit volumes.
The application of the unit cell approach is
based on the fact that most heterogeneous nuclear reactors
have fuel and moderator arranged in a repeating lattice

structure. For calculational purposes a lattice is sub-

divided into a number of identical unit cells (18, 19).
The spatial and speed dependent neutron density, usually
just the scalar density, is calculated for the unit cell

and is subject to particular boundary conditions. Quanti-

ties such as the thermal utilization are then calculated

for the unit cell and are assumed to be valid throughout

the lattice. The reactor is then treated as a homo-

geneous reactor having the same thermal utilization as

the unit cell, and the scalar neutron density is calcu-

lated for the entire reactor.

In most reactors the lattice is very large

relative to the dimensions of the unit cell, so that

it is not unreasonable to treat the unit cell as if it

were in an infinite lattice. For a complete discussion

of this assumption, see Chapters 7 and 18 of Weinberg

and Wigner (19). If the assumption is then made that

the lattice is infinitely repeating, then at steady-

state the flow of neutrons in dn about direction ni

across a point on the boundary between adjacent unit

cells is equal to the flow in d_ about direction i

where L2 is the mirror image of nl as illustrated in

Figure 3. This is physically equivalent to the

assumption that the cell boundaries are perfectly

reflecting. Hence in the analysis that follows,

neutrons which actually arrive at r in the unit cell

by crossing the cell boundaries on first-flights can

be considered to have originated within the unit cell

at r' and suffered reflections at the cell boundaries

before arriving at r.

S2 180- - - - -

T2 T 1

Figure 3 Symmetry of Neutron Flow at the
Boundaries of Unit Cells in An

problem in the following manner. Due to reflections at

the cell boundaries, first-flight neutrons from an
isotropic point source at r' can arrive at r by various

paths, so that the total density at r is the sum over
all possible paths i, of the densities of neutrons
arriving at r by path i and having traveled a total

path length Li.
Alternatively, the neutrons at r can be con-
sidered to have arrived by first-flights from i image
sources at r located a straight-line distance L from
soure i at located a straht-ine distance L from

the field point r. The total angular density at r is
then the sum over the i image sources of the contri-
butions to the density at r from each image source at
rt. (See Chapter 7 of Morse and Feshbach (16) for a
discussion of the method of images.)
Thus for a general homogeneous medium, one
has from Eqs. (63) and (99) that

G(v;r,n r') e (t( rv)i)
i 4 TTv2 -i (105)

where Li is the total path length traveled by a neutron
arriving at r along a first-flight path from an image
source at ri. From the alternate point of view, Li is
the path length traveled by a neutron along path i as
it suffered reflections at the cell boundaries before

arriving at r.
Integration of Eq. (105) over n yields

-" t(k',v)Li
G(v;r r') 4 7tvvJLi (106)
-O vLi2 (106)
4 7r vLi1

As one would expect, for an infinite medium or for a
finite medium bounded bya non-re-entrant surface with

vacuum boundary conditions Eq. (106) reduces to Eq. (101);
for a finite medium bounded by a re-entrant surface with

vacuum boundary conditions, Eq. (106) reduces to the

infinite medium Green's function with spatially variable

collision probability as given in Chapter 2 of Case, et al.

(21), namely

R (r sR/R, v) ds
G(v;r r') e 2(107)
4 7 vR2

where Zt(s,v) 0 outside the bounding surface.

VII.2 Analytic Form of G(v;xlx')

The physical effect of introducing reflecting
boundaries is to increase the density of neutrons at a

point due to reflections at the boundaries which, in

effect, send neutrons back into the finite unit cell.

The total neutron density at r may be separated into
two parts; namely, that due to first-flights from the
source at r' within the unit cell and that due to

neutrons that may be considered either as arriving

from i image sources located at r' outside the unit
cell or as having suffered one or more reflections at

the boundaries after emission at r'.

The first-flight angular density at r due to
a unit isotropic point source at r' in a homogeneous
medium is simply the infinite medium Green's function
obtained from Eqs. (63) and (99), namely,

0,,nv) rr-,e
oG (v;r,|r,) e- (a ^ ').
47T v r-r - (108)

Let Gfc(v;r,n r') denote the angular density at
r of neutrons which do not arrive by first-flights
from the unit isotropic point source at r' but instead
arrive at r only after having suffered reflections at
the mirror boundaries of the cell, or alternatively,
arrive by first-flights from the i image sources
located at ri. Then let Li,n denote the total path
length traveled by a neutron along a path having a
direction that is characterized by the unit vector
ri-' pointing from the image source at r' to the
j -ji
field point r, the neutron arriving at the field point
r after having suffered n reflections at the mirror
boundaries of the cell. One can then write from Eq. (105)

Gfc(v;rnE I) z e- lt(k',v)Li,n _
n-l i 4 VLln 2 --
n (109)

The double summation over n and i is, of course, equiva-

lent to a single summation over all the image sources.

The total Green's function for the unit cell is then

G(v;r,n Pr') G (v;r,nlr') + Gfc(v'r,rl r') (110)

Thus it is seen that Gfc represents a finite medium

correction term to be added to GO.

In V.1 the plane geometry equations were

obtained by integrating the point source kernel over

all the sources on the yz-plane at x'. For a unit cell

with plane reflecting boundaries, a neutron which

arrives at r after having suffered n reflections, and

thus having traveled a path length Lin, can be con-

sidered to have originated at a fictitious or image

plane source located at a distance pLi,n from the field

plane. The cell boundary at which the neutron made its

last or n-th reflection lies between the field plane and

the image source plane.

The transformation from the point source

Green's function to the plane source Green's function,

that is Eq. (68), must be modified using the above

method of images to locate image plane sources at

distances pLi,n from the field plane. Referring to

Figure 2, G(v;x,plx') is the number of neutrons with

direction cosines between p and + du crossing a ring
element of area 27 Tper unit time, the ring element
being located in the plane at xt. If f is the azimuthal
angle indicated in Figure 2, then

dr du dP (111)

Hence, using the transformations of V.1,

Goo (v;x, Ix')dx dp dx' ,

2 0
x d d G I
fma0 f-2l1 (112)

Since the sources are assumed to be symmetri-
cal about the x-axis, G(v;r,flr') is a function only of
r,r', and u. Thus the integration over c may be
carried out in Eq. (112):

GC (v;x,p x') 47r2/ G(v;rl|r') X- 2 (113)

P O =0

For slab geometry G(v;r,nlr') has a delta function
behavior with respect to p' so that Eq. (108) becomes

Gc (v;rn.r') -

e- t(k',v) (x-x)/p| e r- -r .
4----x- ---------(114)
47yv(Ix-x'l /p)2 27r

The p'-integration may then be carried out in Eq. (113).
The result is

Gc(v;x,pIx') -

- (k',v) r-r


or using Eq. (2),

G C(v;x,p x') -

-Z'(k'.V) I (X-X'f)^l
2v II

-Zt(k,v) (x-x')/P

h(x-x') for f > 0


h(x'-x) for < 0

2v I

where the step function has been introduced to account
for the fact that in an infinite first-flight medium,
the flow of neutrons in the positive P direction is
zero unless x is greater than x', while the flow of
neutrons in the negative p direction is zero unless x
is less than x'. For a more complete discussion see
Chapters 2 and 4 of Case, et al. (21).

In the same manner, one obtains the correspond-
ing expression for the finite medium correction term.

Zt(k',v) (x-xi)/p
GfC(v;x,pl x) e h(x-xi) for p 0
1 2v ||
g(k v) I(x-x) I
Ee -,- / h(xi-x) for p 0
i 2v i

where the single summation represents a sum over all i
image source planes located respectively at xi. The
first two image source planes are illustrated in
Figure 4.
Substituting Eqs. (115) and (116) into
Eq. (110), one can verify that each term in G(v;x,p x')
satisfies the one-dimensional form of the source-free
transport equation corresponding to Eq. (95), namely

oG(v;x,p x')
VU xy + v Zt(k',v) G(v;x,Ix') 0 (117)

where v Zt(k',v) is equivalent to [V(k',v)+ (k' ,v) .
Further, G(v;x,p x') satisfies the reciprocity theorem
(21), that is,

G(v;x,p x') G(v;x' x,-) (118)

Finally, each term of Eqs. (115) and (116) satisfies

the symmetry condition at the boundary x +b that the

flow of neutrons in dp about p is equal to the flow in

dp about -', that is,

G(v;x,u x') x-l b G(v;x,- ix') .x-b (119)
x-+b x-+b

Since each term of the Green's function G(v;x,ujx')

satisfies the original boundary conditions, it can be

concluded that the function itself satisfies the boundary

conditions. One can verify that, in agreement with

sections 4 and 5 of Case, et al. (21), a source plane
of strength qs at x' is equivalent to a discontinuity
of the normal component of the angular current, where

for an isotropic plane source, qs 1/2. The angular

current is simply the number of neutrons crossing
unit area perpendicular to the direction of flow per
unit time and unit solid angle. In particular, for
the infinite medium terms one can write

vpG(v;xuIx') -x-x'+- vuG(v;xGiux') x -x'- q5s() "

In an infinite medium the second term on the LHS

vanishes for p >0 while from Eq. (115) the first

term is equal to 1/2, that is, the source plane value.

In the solution for the scalar neutron density
in Eqs. (69) and (77) one needs the corresponding
scalar Green's function G(v;x x') obtained by integrat-
ing over all p the function G(v;x,p x'). Thus,


G(v;x x') G(v;x,p x')dp (120)

Recall that G(v;x x') is the first-flight
density of neutrons at x due to a unit isotropic plane
source at x' in a homogeneous medium. Then, referring
to Figure 4, one can see that for image sources to the
right of x, that is, x => x, the neutrons reaching x
must have a negative direction cosine u, while for
image sources to the left of x, that is, xi <= x, the
neutrons reaching x must have a positive direction
cosine. Thus the contribution of the image sources is

0 1

Gfc(vI;x x) = Gfc(v;x,plx.) dp + Gfe(v;x,,Ix')di (121)
0- 1 0

where the first and second integrals represent the
contributions of image sources placed respectively to
the right and left of x.


Image Source
Plane at

Figure 4






x-b x-2b-x'
Image Source
Plane at

Slab Unit Cell with Mirror
Boundaries at x-+b

Letting Gtc denote one term of Eq. (116),
one may write the first integral in Eq. (121) as

0 0
f fe(v e- t(k'tv) I(x-X)/IU
Gi (v;x, |')du -d
iI2vl I2

e- Zt(k',v)(x-xi)/l d (122)

the last integral having been obtained from the preced-
ing one by a successive substitution of variables, first
letting y--u, and subsequently letting i-y.
In general then, one can sum over all image
sources and obtain the result that

G (v;x x') f e t(k',v) (x-xi) d (123)

Again, for a first-flight medium the range
of integration over ) for the infinite medium contri-
bution is from 0-O to P1. Thus one can write

G(v;xlx') f e- t(k',v) (x-x')/ d

+ -t(ck',v) (x-x') 4)
i 2v (124)

or from Eq. (103),

G(v;x x') El ( t(k'v)x-x'j]

+ E El[Zt(k',v) x-xil] )

To express G(v;xIx') as a function of the
unit cell dimensions it will be instructive to return
to the physical concept of reflection of neutrons in
order to locate the image sources with respect to the
cell boundaries. Using the notation for path length
previously discussed with regard to Eq. (109), one can
write the expression for GfC(v;x x') in Eq. (125) as

Gc(v;x x') E E El[7t(k',v)ILi.nl] (126)

where PLi,n is, as before, the distance separating the
field plane from an image source plane, the image
plane corresponding to a neutron having made n re-
flections at the mirror boundaries.

In Eq. (126), the summation over i extends

over all possible paths by which a neutron having its

direction cosine equal to u can arrive at the field

point after suffering n and only n reflections at the

boundaries. In plane geometry, it is easily seen that

for a first-flight medium there corresponds to any given

p only two such paths for any n. Hence the summation

in Eq. (126) runs from i-1 to i-2. Figure 4 illustrates

the two possible paths, L1,1 and L2,1 corresponding to

S= 0 and u < 0 respectively, by which a neutron can

arrive at x after suffering only one reflection. The

unit cell in Figure 4 has its midplane at x-O and the

reflecting boundaries at x b. The cell width,

that is, the distance between the two boundaries is 2b.

The expressions for pLi,n can be obtained by

simply adding the total path length a neutron travels

while making n reflections. The expressions for pi,n

depend on the orientation of the source plane and field

plane x relative to each other and to the x-0 plane.

It has been shown in V.2, however, that one need only

consider field planes for positive x, because the

solution of the half-cell is symmetrical with respect

to x.

The expressions for aLin are derived in
Appendix A where it is shown that Eq. (125) may be
written, for x > O, as

G(v;xx') 1. E (k',v)Ix-x']

+ 1 EI[Zt(k',v) 2nb-x'+(-1)nIxl
2v n-1
+ E[Zt(,v)12nb+x'-(-1)nx]) for x = 0.

Note that Eq. (127) can be shown to be term by term
equivalent to a more cumbersome series representation
of G(v;x x') derived by Aswad (25).
Finally, the half-cell Green's function
defined by Eq. (78) becomes

H(v;xlx') = E[l t(kO,v)|x-x'I] + E1[t(k',v)x+x'I

1 co
+ I E 1I[Zt(kv) 12"b'x-'l)n]

+ EL[Zt(k',v) 2nb.+x-(-1)nxi]

+ E,[Zt(k',v) 2nb+x'+(-l)nx]

+ El[Zt(k',v)12nb-x'-(-l)nxi] )(128)

for x =0 and x' O 0

VII.3 Monte Carlo Generation of G(v;xlx')

For complex finite geometries, the complete

analytic form of the first-flight Green's function may

be exceedingly difficult, if not impossible, to derive.

In such cases, it may be advantageous and even necessary

to represent the physical or mathematical system by a

sampling operation satisfying the same probability laws.

Such a process has come to be called a Monte Carlo

method (26).

Essentially, Monte Carlo methods are "paper

experiments." The experiment consists of performing

specific sampling operations many times, the sampling

operations satisfying the criteria mentioned above. As

in many other experimental procedures, the geometry

of the problem does not greatly affect the difficulty

of performing the Monte Carlo experiment. The con-

struction of a function, say neutron density, by

Monte Carlo is conceptually little more difficult for

a square or hexagonal unit cell than it is for a slab

or cylindrical unit cell. This is in striking contrast

to the analytic description of the same function;

there the solution for simple geometries does not, in

general, give any indication of how to proceed to more

complex geometries. Thus there is strong motivation

to generate the plane geometry Green's function G(v;x x')
by Monte Carlo. The transition to cylindrical or even
more complex unit cell geometries would not be difficult;
the general procedure, that is, the Monte Carlo sampling
operations would be essentially the same.
The Monte Carlo generation of G(v;xix'), that
is, the first flight Green's function for a homogeneous
finite medium in an infinite plane geometry with
perfectly reflecting boundaries has been discussed
extensively by Aswad (25) and will be outlined below.
In plane geometry, a cluster of neutrons is
considered to be emitted with unit statistical weight
at the source plane x'. Its direction cosine, P, is
chosen from a random distribution, the distribution
being such as to represent an isotropic source, that is,
a source which has equal probability for emission of
neutrons in an element of solid angle dn about any
direction n; thus in plane geometry the direction
cosine P occurs with equal probability for values
between -1 and +1. The random selection of P is
done by Monte Carlo sampling.
The projection of the path of the neutron
cluster on the x axis is divided into closed intervals

[Xm-l, xmJ [xm xm+1 ... and the statistical

weight of the neutron cluster is computed at points

Xm-l, xm, xm+l, *. The weight of the neutron cluster
at some point Xm, denoted by W(xmjx'), is subject only
to exponential attenuation along the path of travel, the
mean free path between removal collisions being equal
to the reciprocal of the total collision probability
Zt(k',v). Thus

W(xmIx') e- t ( ')/p (129)

If xm+l > xm x', and if W(xmlx') and W(xim+xl')
represent respectively the statistical weight of the
neutron cluster at x, and xm+1, then the quantity
[w(xmJx') W(xm1lx')] must represent the number of
removal collisions per unit time that occurred in the
closed interval [xm, xm+ .
The neutron cluster continues along its
original path until it arrives at a cell boundary
where it is considered to be reflected and sent back
across the cell with direction cosine p', where u*--u,
that is, Ip' Iu." The number of collisions that
occur in the interval [xm, Xm+1 as the neutron cluster
passes through that interval is again recorded. The
process is repeated as the neutron cluster suffers

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