NUCLEAR DESIGN OF A SUBCRITICAL ASSEMBLY DRIVEN 
BY ISOTOPIC NEUTRON SOURCES 
BY 
MAAKUU BULMUO TUOR
A DISSERTATION PRESENTED TO THE UNIVERSITY OF GHANA,
LEGON,
FOR THE DEGREE OF MASTER OF PHILOSOPHY 
(M. PHIL) IN PHYSICS
SEPTEMBER 1993
University of Ghana          http://ugspace.ug.edu.gh
*37333
"TLj24e £  d ^ r v \
University of Ghana          http://ugspace.ug.edu.gh
DECLARATION
I declare that, except for references to other people 
this thesis is the result of my own research and that 
neither in part nor in whole been presented elsewhere for 
degree.
A .
(B. T. Maakuu) (Dr. E. H. K. Akaho)
Student Supervisor
s work, 
it has 
another
University of Ghana          http://ugspace.ug.edu.gh
DEDICATION
This work is dedicated to my father and mother 
for their love and care
University of Ghana          http://ugspace.ug.edu.gh
TABLE OF CONTENTS
Abstract (i)
Acknowledgements (ii)
CHAPTER ONE: INTRODUCTION 1
CHAPTER TWO: LITERATURE REVIEW 4
2.1 Introduction 4
2.2 Basic Reactor Physics 5
2.2.1 Reactor Physics Equations 6
2.2.1.1 Neutron Transport Equation 6
2.2.1.2 Neutron Diffusion Equation 8
2.2.2 Methods of Solution of Diffusion 
Equation 11
2.2.2.1 Analytical Method 11
2. 2. 2.2 Numerical Method 19
2.2.2.3 Finite Difference Method 21
2.3 Subcritical Assemblies 24
2.4 Principles of Activation Analysis 29
2.4.1 Neutron Activation Analysis Equipment 31
2.4.2 Isotopic Sources 34
2.5 Concluding Remarks 3 6
CHAPTER THREE: THEORY OF THE CODE SUNDES 3 8
3.1 Introduction 3 8
3.2 Mathematical Model 3 8
3.2.1 Matrix Form of Multigroup Equations 51
3.2.2 Computational Procedure 53
3.2.3 Inner Iteration - 58
University of Ghana          http://ugspace.ug.edu.gh
3.2.4 Normalization of Fluxes 60
3.2.5 Input Description 61
3.2.6 Output Description 63
3.3 Concluding Remarks 64
CHAPTER FOUR: NUCLEAR DESIGN OF THE NEUTRON MULTIPLIER 65
4.1 Introduction 65
4.2 Application of The Code 65
4.2.1 Single and Two- Region Problem 65
4.2.2 Application to Multiregion Problem 69
4.2.3 Description of The Accepted Geometry 88
4.3 Description of The Mechanical Features of
The Neutron Multiplier 89
CHAPTER FIVE: CONCLUSIONS AND RECOMMENDATIONS 93
REFERENCES 97
APPENDIX 103
University of Ghana          http://ugspace.ug.edu.gh
ABSTRACT
A feasibility study of a conceptual nuclear design of a 
facility capable of producing high thermal neutron fluxes 
using isotopic neutron sources in a multiplying medium was carried 
out. The one-dimensional multigroup neutron diffusion equation was 
solved using the finite difference technique. A computer code 
SUNDES was written in FORTRAN 77 programming language and used to 
study the effect of reflectors, shielding materials and strength 
of isotopic neutron sources on the production levels of neutron 
fluxes. The neutronic calculations showed that with a homogeneous 
mixture of 20% enriched U02 and Be as multiplying medium, BeO as
reflector, Al as cladding and concrete as shield, thermal neutron
7 2fluxes as high as 1x10 n/cm -s could be produced. Simultaneous 
irradiation of samples in two different regions at different 
fluxes is possible in the assembly. The analysis also revealed 
that different orders of thermal neutron fluxes could be produced 
depending on the strength of the driving isotopic neutron sources.
i
University of Ghana          http://ugspace.ug.edu.gh
ACKNOWLEDGEMENTS
I wish to express my profound gratitude to Dr. E. H. K. Akaho, 
under whose supervision the work for this thesis was carried out, 
for his unqualified guidance and concern. My gratitude also goes 
to Prof. E. K. Agyei and Dr. K. A. Owusu-Ansah ray co-supervisors, 
for their valuable suggestions and words of encouragement 
throughout this endeavour. Copious thanks are due to Profs. J. K. 
A. Amuzu, G. K. Tetteh, Dr. E. K. Osae and all the lecturers of 
the Department of Physics, University of Ghana, Legon, for their 
concern and moral support.
I gladly acknowledge my colleagues, Messrs Justin Agbotse, 
Anim-Sampong and Guggisberg Amoh for their various and much 
needed assistance.
I extend my heartfelt thanks to the Director of NNRI, GAEC, 
for allowing me the use of their facilities. To Dr. K. A. Danso, 
Dr. H. 0. Boadu, Mr. M. S. Mahama all of the Department of Nuclear 
Engineering and Mr. G. Emi-Reynolds of the Radiation Protection 
Board, GAEC, I say many thanks for your support.
Finally, I thank all friends and well wishers who in diverse 
ways contributed to the successful completion of this work.
il
University of Ghana          http://ugspace.ug.edu.gh
CHAPTER ONE
INTRODUCTION
Neutron activation analysis (NAA) is one of the most powerful 
nuclear techniques for multi-element analysis. Though the first 
activation analysis in history was performed in 1936 [1] , The
technique was not widely developed until the 1960s when Ge(Li) 
detectors with high resolution and efficiency were developed. The 
technique has since become an important means for the super-trace, 
trace, semi-micro and normal analysis. Developing countries with 
small laboratories cannot explore this technique fully due to 
their inability to acquire nuclear reactors and generators to 
produce high neutron fluxes for the purpose. Such laboratories 
most often use isotopic sources such as Am/Be, Pu/Be etc in a non­
multiplying medium. These devices are not capable of producing 
high thermal neutron fluxes for the analysis.
In this study a computer model for the nuclear design of a
subcritical assembly driven by isotopic neutron sources in a multi
-plying medium (neutron multiplier) which is capable of producing
7 2thermal neutron fluxes of the order of 10 n/cm -s is provided. It 
is hoped that this device will be a substitute for high neutron 
flux generating devices for certain neutron activation analysis.
In the neutronic calculation for the design of the 
subcritical assembly, the one-dimensional multi-group neutron 
diffusion equation is solved. It is a conservative equation which 
takes account of both neutrons gained in the assembly through
1
University of Ghana          http://ugspace.ug.edu.gh
fission, isotopic sources, scattering processes and those lost as 
a result of diffusion and structural absorption. Several methods 
of simulations such as the Monte Carlo method, numerical method, 
analytical method etc are available [2]. The Monte Carlo method is 
rather too complex for this analysis. It is difficult to use the 
analytical method to solve practical problems for multi-region 
geometries. The numerical method transforms the analytical 
equations characterizing the system into a set of algorithms and| 
numerical equations. Based on these algorithms and numerical v 
equations, computer programmes can then be developed.
Two well known approaches are available in the literature for 
numerical simulation; the finite element method and the finite 
difference method. The finite element method uses a triangular or 
tetrahedral element to set up a variational formulation of the 
problem which is then solved by optimization technique [3] . The 
finite difference method covers the region under consideration by 
a mesh consisting of horizontal and vertical lines and seeks the 
approximate values of the solution at the intersection [4] . For 
this work, it was found that the finite difference technique is 
more appropriate for the analysis.
In Chapter Two of this work, a review of the basic nuclear 
reactor physics and methods used for calculations will be 
presented. Special reference will be on critical and subcritical 
assemblies. The analytical and numerical methods of solution to 
the neutron diffusion or transport equations are discussed. The 
early works on subcritical assemblies and their importance are
2
University of Ghana          http://ugspace.ug.edu.gh
also presented. The basic principles of neutron activation 
analysis, its applications and description of a typical neutron 
activation analysis experimental set-up located at the National 
Nuclear Research Institute of Ghana Atomic Energy Commission are 
discussed. Last but not the least, the aims and objectives of this 
work are outlined.
In Chapter Three, the mathematical model developed using the 
multi-group diffusion equation is presented. The computational 
flowchart and algorithms based on the finite difference technique 
are presented. A description of the input and output of the code 
are also given.
Chapter Four discusses the results of the preliminary 
investigations of the nuclear design. This is followed by a 
description of the neutron multiplier. Finally, a description of 
the mechanical features of the neutron multiplier and its mode of 
operation are presented.
The conclusions and recommendations on the work are contained 
in Chapter Five.
3
University of Ghana          http://ugspace.ug.edu.gh
CHAPTER TWO
LITERATURE REVIEW
2.1 INTRODUCTION
Since the advent of nuclear reactors, neutron activation
analysis has become a powerful tool, for multi-element analysis.
This nuclear technique has however been an elusive expedition for
most developing countries and small laboratories. They cannot
afford nuclear reactors which produce high neutron fluxes. Other
high neutron flux generating devices such as accelerators are
equally expensive. Alternative sources must therefore be sought
for the purpose. It is envisaged that a form of subcritical
assembly could be designed to achieve thermal fluxes greater than 
7 21 x 10 n/cm -s for neutron activation analysis.
The area of computational reactor physics is extensively 
studied in various nuclear institutes. A wealth of literature 
exists and can be found in reactor physics books [2,6,7,8]. In the 
present review, not all the detailed information in the literature 
will be considered. A review of the nuclear reactor theory and 
methods used for reactor physics calculations for the determination 
of neutron fluxes in nuclear reactors with special reference to 
critical and subcritical assemblies is first presented. The analy­
tical and numerical methods and techniques that are normally used 
to provide solutions to neutron diffusion or transport equations 
are also discussed. This is followed by a discussion of the early 
works on subcritical assemblies and their importance. Next is a
4
University of Ghana          http://ugspace.ug.edu.gh
discussion on the principles underlying neutron activation techni­
que and a description of a typical neutron activation analysis 
experimental equipment using Am/Be source immersed in a pool of 
de-ionized water located at the National Nuclear Research 
Institute, Kwabenya-Accra.
Finally, the aims and objectives of the conceptual nuclear 
design of a subcritical assembly driven by isotopic neutron 
sources are stated.
2.2 Basic Reactor Physics
Basic nuclear reactor physics is a very broad field of 
physics covering cross sections, transport theory, diffusion 
theory, reactor kinetics, multi-group theory etc [6, 7].
For the purpose of this work, emphasis will be on the 
mathematical methods for analyzing the behaviour of neutrons 
namely the transport theory and the diffusion theory. The 
analytical and numei'ical methods of solutions will also be 
discussed.
The detailed behaviour of neutrons in a nuclear reactor may 
be formulated mathematically by considering the production and 
collisions of neutrons of a particular energy moving in a particu­
lar direction and then integrating over all directions and energies 
Neutrons are born at fast energy and slowed down to thermal energy 
as a result of a large number of collisions. Some are absorbed and 
others diffuse out of the reactor. Their distribution pattern can
5
University of Ghana          http://ugspace.ug.edu.gh
thus be predicted by mathematical models based on either the 
transport theory or the diffusion theory.
2.2.1 Reactor Physics Equations
2.2.1.1 Neutron Transport Equation
The neutron transport equation is a linear equation. It is 
fundamental and exact in describing the neutron population in 
nuclear reactors. Though the method is expensive to solve on 
digital computers, it is often used in areas near boundaries of 
reactors or near highly absorbing materials such as fuel rods or 
coolant elements where the diffusion equation is inaccurate.
The energy dependent Boltzmann transport equation may be 
written for a steady state as [2, 9].
fi.V0(r,E,Q) + Z t </>(r,E,Q) = dE'
/
dC2 E  ( r , E'-»E, <p (r,E' , Q ' )
S
E 1
+ S(r,E,Q) (2 .1)
where 0(r,E,Q) is the number of neutrons of energy E per unit 
energy crossing a unit surface at position r per unit time going
in a unit solid angle centered in the direction Q. E (r,E'-»E,s
is the probability per unit path length that a neutron at r and 
going in a direction Q' with energy E' is scattered into a unit 
solid angle centered at Q and a unit energy interval centered at
6
University of Ghana          http://ugspace.ug.edu.gh
E. S(r,E,fi) is the number of neutrons with energy E created per 
unit volume at position r going in a solid angle centered at Q.
The multi-group equations can be cast by dividing the energy 
range into G groups of arbitrary width AE = E - E . For simpli- 
city let the initial and final energies be labelled Eq and E^ 
respectively. Let also the thermal group be indexed by G+l. 
Equation (2.1) can thus be integrated from E to E and simpl-
ified to obtain
Q.7<£g (r,E,Q,) + Zg (r.£2)0g (r,Q) = Qg (r,fi) + Sg (r,Q) g = 1,2,..., (
(2 .2 )
Where
Q (r, E, £2) = dE' dD'E (r, E'-»E, Q'-»Q)s
and the group cross sections averaged over the appropriate flux is:
2g (r,n) = dEZfc(r,E)0 (r,E',Q)
sg-i
„E
dE0 (r,E,Q)
V i
4*i,:
%
and
7
University of Ghana          http://ugspace.ug.edu.gh
Qg (r,fi) dEQ(r,E,n)0(r,E,fi)
3g-i
dE0(r,E,fi)
ig - i
The thermal equation is formulated to read
Q.V0G+1 (r,C2) + ZG+1(r)0 (r,Q) = QG+1(r,fi) + SG+1(r,Q) (2.3)
Using equations (2.2), (2.3) and boundary conditions the desired
multi-group transport equations are approximated as
Q.V<pg ( r,Q) + S^(r)0(r,n) = Qg (r,Q) + Sg (r,fi), g= 1,2 ...,6+1
(2.4)
2.2.1.2 Neutron Diffusion Equation
The neutrons in a reactor have energies spanning the range 
from lOMev down to less than O.Olev. The neutron-nuclear
cross sections are sensitively dependent on the incident neutron
energy. The neutron energies are discretized into energy intervals 
or groups. A multi-group diffusion equation can be successfully 
arrived at by adopting the concept of neutron balance for any
given energy group. An account is taken of the way neutrons can
8
University of Ghana          http://ugspace.ug.edu.gh
enter or leave this group. For a typical group, g the neutron 
balance is represented as:
'
Time rate of change Change due
■
Absorption Sourceneutrons
of neutrons in group g to leakage in group g
+ appearing 
in group g^
f
Neutrons
r
Neutrons
scattering + scattering (2.5)
out of group g into group g
This neutron balance equation is mathematically represented as:
1 dip 
v Stg
G
V.D V<p -  X tp + S + Z 0 + > Z  , <j> , g=l,2,g Yg ag^g g sg^g sg'^g^g'
g' =1
. , G (2.5)
where the source term due to fission is given by the expression 
G
S = x  ) v  ,Y .r ,cp , +  S
g gzL g f g g
ext
g (2.7)
g' =i
To satisfactorily account for the average behaviour of 
neutrons in each group, the energy-dependent diffusion equation
_L_ dA - V.DVcfr + U(r,E,t)
at ^
dE,Ea (E'->E)0(r,E,,t) + Sext(r,E,t)oO
9
University of Ghana          http://ugspace.ug.edu.gh
r“
+X (E) dE'v,(E')Sf (E')0(r,E / , t) (2.8)
J 0
is integrated over a given energy group, E -,< E < E . Other
y-J* y
energy-dependent equations could also be used as a starting point 
for the development of the multi-group diffusion equation [9, 10]. 
Further simplification gives the multi-group diffusion equation as:
G G
f£g = V.D V0 - Z. 0 + V s  , <p , +  X V V ,0 , + S (2.9)v St^  g tg^g sg'->gyg' 'cg(/_ g' fg,vg' g
y g'=l g'=l
g = 1, 2...,G
If it is assumed that neutrons can only scatter to the lower 
groups, then the scatter term may be simplified to read:
G G-l
, 0 , = > 2 , 0 , + Z 0 (2 .10)sg ->g g sg'->grg' sg->gyg
g'=l g'=l
s i  09
and hence the multi-group equation is rewritten as'p' /
A
G-l .. . .
— —  ^ g  = V.D ?0 - E. 0 + > £  , 0 , + E  0 + Sv ^  g *g tg^g sg'-^g' sg^g^g g
g u g'=l
G
+ * g Z v zf g ' V  g= 1,2 G (2-11}
g'=i
10
University of Ghana          http://ugspace.ug.edu.gh
For a steady state the LHS of equation (2.11) becomes zero thus,
G-l G
0 = V.D V0 - Z. + >Z , d> , +  Z <b + S + X > V , Z r , d ,g tg^g sg'->g*g' sg->gvg g g/_ g' fg' g'
g'=l g'=l
(2 .12)
In practice, only a certain number of groups are dealt with 
in reactor physics calculations. The most commonly used is the two 
multi-group calculations. In collapsing the groups for group 
constants calculations, the particular type of reactor being 
analyzed and its operating conditions with respect to fuel loading, 
isotopic composition, temperature and coolant conditions are
considered. Analytical and numerical schemes can be employed to 
provide a solution to the multi-group equation. A brief outline of 
the analytical and numerical methods of solution to the neutron 
diffusion equation is presented. Emphasis will be placed on the 
finite difference method, a numerical technique which is of
interest in this work.
2.2.2 Methods of Solution of Diffusion Equation
2.2.2.1 Analytical Method
Analytical solutions to neutron diffusion equation are
provided in the literature. In this work, the solution as related
to subcritical reactor physics will be discussed.
A feasibility study using a neutron source located at the
11
University of Ghana          http://ugspace.ug.edu.gh
center of a multiplying medium for different fuel to moderator 
ratio was carried out by Akaho and Danso [11]. The Fermi age theory 
[12] was used to study the flux production levels in unreflected 
subcritical assemblies. Solid moderators as Be, BeO and graphite 
were used with different fuel enrichments.
For a homogeneous mixture of fuel and moderator in a 
cylindrical assembly the time dependent thermal equation is 
written as
Sn(r,z,t) = D720T(rjZ(t) Za0(r,z,t) + S(r,z,t) (2.13)
where n(r,z,t) is the number of thermal neutrons and is related to 
the thermal flux by the equation
$(r,z,t) = 2vn(r,z,t) (2.14)
V~n
Substituting equation (2.14) into (2.13) gives the expression
Vn d c p ( r , z  , t )  = DV20(r,z,t) E 0(r,z,t) + S(r,z,t) (2.15)
2v St
Further simplification of equation (2.15) gives:
L2V2$(r,z,t) <p( r,z,t) + |q(r,z,rT,t) = td |^vr,z't5 (2.16)
a
where q(r,z,T , t) is the slowing down density in the absence of
12
University of Ghana          http://ugspace.ug.edu.gh
resonance absorption. The number of neutrons which become 
thermalized is pq(r,z,tt, t) . By the dictates of Fermi age theory 
the equation
a2q(r,z,rT,t) = 3q(r,z,xT,t) (2.17)
Sr2 St
must be satisfied. It will be assumed here that the neutrons from 
fission and isotopic sources are emitted mono-energetically at 
zero lethargy (t = 0) . In the slowing down model, neutrons start 
to slow down soon after being emitted from the sources. The source 
condition satisfied at t = 0 is:
SS (r) S (z)
q(r,z,0,t) = ---------- + TjfeZ 0 (r,z,t) (2.18)
n r l  a
The total number of fission neutrons produced per cm at the point 
(r,z) is accounted for by the second term on the RHS of equation 
(2.18) . Substituting k /p for 7)fc in equation (2.18) gives:
SS(r)5(z) k^ Zl 0 (r, z, t)
q(r,z,0,t) = — --------  + ---------
7rrl p
(2.19)
Assume that in a cylindrical geometry
* (r,z,t) = ^ A mn(t)Jo 
n=odd
' '
X r mrzn cos
R H
(2 .20)
13
University of Ghana          http://ugspace.ug.edu.gh
and
q(r,z,T,t)
= r
C (T, t) J mn o
' r
X r n n zn cos
R HV
(2 . 21 )
— B  T 7where C „ ( z , t )  = T e mn mn mn
For an infinite cylinder of radius R and height H, the infinite 
source at r = 0 is given by the equation:
S (r, z) =
Sfi(r)S(z) ^  Jo (Xnr/R)
7lR2 H - In=l J.(X ) n=odd 1 n
f '
X r n n z
n cos
I R . H
(2 .22)
Substituting equations (2.20), (2.21) and (2.22) into (2.19) gives
the following expression:
(t)e BmnTJ__ mn 
n=odd
r '
X r nn zn cos
R H
2i  V  J o ( x n ^ E > V j J X n r  
vn^ — J 2 , A 1,0 Rn=l J* (Xn) n=odd
( •
X r n n z
n cos
I R . H
+ oo aV^A (t)J 
—  Z_ mn o
p n=odd
' r
X r n n zn cos
R H
(2.23)
From the neutron diffusion equation for a cylindrical geometry
V I  = -B 4 T mn’ (2.24)
14
University of Ghana          http://ugspace.ug.edu.gh
where mn
r r 2iX nnn +
[r ' H'
. R' and H' are the extrapolated radius
and height of the cylinder respectively. Substituting equation 
(2.20), (2.23) and (2.24) into (2.16) and replacing t with x T gives
the expression
-L2 V b 2 A (t)J T / mn mn o
n=odd
r - \
X r nrczn cos
R' H' J
n  _ V"
TT
A (t) J mn o
n=odd
' '
X r n n zn cos
R' H'
anfeddl I k  J1 (xl> P
. k Z A (t)1 oo a mn+ ~B2 TT e mn J
'
X r n7TZn COS
R' H'V.
= t £
dA J 
dtmn °
•
X r nTTZn cos R'R' (2.25)
Simplifying equation (2.25) further gives the expression:
mn
v -B2 X\ k e mn(X)
i-l2b2T mn
I
2
00 J (X r/R') k e mn77 o m oo2PS
£ k v/ ^2 ,
a ” ^ 0  J1{V
dAmn
(i+l2b2 )dtT mn
(2.26)
-B2 tk e mn T t,a
L e t  kmn= ---------^ ^ —  and  tjnn = -------2— 2 t h e n  e q u a t i o n  ( 2 . 2 6 )
15
University of Ghana          http://ugspace.ug.edu.gh
becomes
dAmn
mn dt - <knrn-1)Am„ +
mn (2.27)
Integrating equation (2.27) subject to the boundary conditions 
0 (0) = 0  and JQ = 1 yields:
A 2PS
k
mn S k v a to
mn V  1
- 1 / t2mn L— - «T i
0 mn tmn 1
|1— ■ J, (X ) n=l 1 n
(2 .28)
The expression for the thermal flux is thus obtained by 
substituting equation (2.28) into (2.19) to give:
-(r,z,t) 2 PSZ k v a oo
00
[ kmn
h ■ kmn-1
(k -1)te mn mn - 1 x
' t i
X r m cos n n z
R' H'
(2.29)
The time behaviour of the flux following the insertion of the
source depends on the magnitudes of kmn. The eigenvalue of
increases monotonically while kmn decreases monotonically ie
k . >k >k r. .... In subcritical situations, k < 1 thus all the ml m3 m5 mn
time dependent exponents in equation (2.29) will be negative and 
the flux then approaches a steady state as:
16
University of Ghana          http://ugspace.ug.edu.gh
0 (r,z) = 2 PS2 k v a oo
00 CO f
' I  X  i s -m=l,2 n=l' mn J2i(xn) °
' 1 5
X r n n zm cos
R' H'
(2.30)
For m = n = 1
1 (r, t) = 2 PS£ k v a a>
'ef f
1 kef f' (2.405)
2 .405
R' cos
n z
H'
(2.31)
k e ~ B \
where keff 2 2
1-B L_
and B
f 2 f: *2 .405 n
R'
+ H'
For equation (2.31) to be solved, the resonance escape probability 
p, and infinite multiplication factor k etc must be known for the 
fuel-moderator system. These parameters are evaluated from the 
following expressions [13-19].
p = exp N238Ieff' (2.32)
where -^e f f  is the effective resonance integral and is given by the 
expression:
^   X*
Teff “ 3 ’9
N cr m s
TTK 0 .415
.■-/'-J /
I ’M 1
M
(2.33:
The infinite multiplication factor is given by the expression:
k = eiipf, (2.34)
17
University of Ghana          http://ugspace.ug.edu.gh
where the symbols have their usual meanings. 
The fast fission factor f is given by
235
cr + ter238
f = 235
cr
. 238 , / , mter + (r/a) crd CL
(2.35)
N
r = mN and t =
N‘238
u N'235
1 - e , where e is the enrichment.
Table 1 [11, 20] contains the values £ , L and t for various
3. 1 X
moderators required to obtain the flux from equation (2.31)
Table 1: Thermal Neutron Diffusion Parameters of Moderators.
Moderator E (cms X)
T 2 / 2, LT (cm )
t t
Be 1.04 x 10“3 480 102
BeO 6.00 x io~4 790 100
Graphite 2 .40 x
1O
 
«—1 3500 368
Different moderators, fuel enrichment, moderator to fuel ratio and 
source strengths were used on the basis of the above analysis to 
conclude that Be-moderator for 20% enriched uranium has the 
highest infinite multiplication factor. The variation of the 
infinite multiplication factor with different moderator to fuel
18
University of Ghana          http://ugspace.ug.edu.gh
ratios for 20% enriched uranium is shown in Figure 2.1.
This technique of analysis, however, is applicable to 
source(s) placed at the center of the assembly not at various 
regions. The analytical solution is expected to be difficult in 
analyzing multiple source problems in the multi-regions because as 
given in equation (2.28) the source problem is not easily defined.
A review of the numerical approach to the solution of 
diffusion equation which is practical is presented in the next 
section.
2.2.2.2 Numerical Method
Both analytical and numerical methods are used in solving the 
neutron diffusion equation. Basically, the analytical solution is 
suitable for one-speed ie neutrons diffusion in homogeneous 
reactors. In real reactors, however, calculations are based on the 
heterogeneous nature of the core. Not only must one consider 
non uniformities corresponding to fuel pellets, cladding materials, 
moderator, coolant element, but spatial variation as well [6] . 
These complex situations cannot be catered for by the analytical 
method and hence must be discarded in favour of a direct numerical 
solution of the diffusion equation.
In the use of the numerical method , the differential 
diffusion equation is rewritten in finite difference form and the 
resulting difference equations are solved using digital computers. 
Difference equations can be formulated for plane, cylindrical and 
spherical geometries. There are other numerical methods such as
19
University of Ghana          http://ugspace.ug.edu.gh
k
-
i
n
f
i
n
i
t
y
Figure 2'-1: Variation of the factor kco with Different Nm/Nu 
values for 20^ Enrichment of Uranium, [11-1 •
20.
University of Ghana          http://ugspace.ug.edu.gh
the finite element method, Monte Carlo method etc for solving 
diffusion equations. Details of these methods can be found in 
references [2, 3] . Only the finite difference method shall be
discussed.
2.2.2.3 Finite difference method
For simplicity, assume the differential diffusion equation
Dd20
dx
2 + £a0 (x) = S(x) (2.36)
is to be solved numerically subject to the boundary condition for 
a finite slab of width a, 0(0) = 0(a) = 0. But for real reactors, 
it is assumed that group constants D(x) and Z (x) for the various3.
subregions are constant. Suppose the spatial variable x, is 
divided into a set of N+l discrete points as shown in Figure 2.2
Xi-Ai/2 X.+ (Ai+l)/2
-n - X i , X
XN- X.N
Figure 2.2 Mesh Points with intervals
There exists a number of schemes that can be used to generate a
difference eqiiation representing equation (2.36) on this mesh.
Taylor series expansion can be used to derive a central difference 
2 2formula for d 0/dx . However, the commonly used scheme is to
21
University of Ghana          http://ugspace.ug.edu.gh
integrate the original difference equation over any arbitrary mesh 
interval X^ Ai/2 and X.+(Ai + l)/2 say, where Ai is the interval 
between mesh points. Thus integrating equation (2.36) term by term 
yields
,Xi+ (Ai + 1) /2 
dxZ (x) <p (x)
d.
X±-Ai/2
E  . 0 .airi
Ai Ai+1 
2 + 2 (2.37)
„Xi+ (Ai+1)/2
dx d D(x)d# _ D (x) d<j> 
3x dx ~ dx
Xi~Ai/2
X. + (Ai+1)/2
X^-Ai/2
(2.38)
and ,Xi+ (Ai + 1) /2 
dxS(x) 
X^-Ai/2
S ,i
Ai , (Ai+1)
2 + 2 (2.39)
Equation (2.38) is simplified further by applying a simple two 
point difference formula on ^  to give
d0
dx ^i+1 ^i (2.40)
X±+(Ai+1)/2 Ai+1
d <p 
dx
Xi-Ai/2
V  (
Ai
i-l
(2.41)
22
University of Ghana          http://ugspace.ug.edu.gh
and using the centered average for D;
D(X± +(Ai+l)/2) = § Di+1 + Di i, i + 1 (2.42)
D(Xi-Ai/2) = ~ Di-1 + Di (2.43)
Substituting equations (2.37) (2.43) into (2.36) a set of diffe­
rence equations are obtained as:
ai,i-l^i-1 + ai,i+l^i+1 Si i = 1,2 ...,N-l (2.44)
where
ai,i-l
D ,+D. . 1 i-l
Ai Ai+Ai+1
a . . = S +1,1 a
Di _ + D . D. _ + D . l+l 1 + l-l l
Ai+1 Ai Ai+Ai+1
and
ai,i+1
D. . + D. l+l l
Ai Ai+Ai+i
Similarly, equations of the type as in Eq.(2.44) with different 
coefficients can be derived for two and three dimensional problems 
[6]. The essential task involved in solving the 3-point finite 
difference equation (2.44) is to select a suitable algorithm for
23
University of Ghana          http://ugspace.ug.edu.gh
its solution on digital computers. It can be solved either by the 
Gaussian elimination or by any other iterative method as the
Subcritical assemblies were designed for the investigation of
definite problems dealing with neutron physics and core technology
[21-23]. In these devices a chain reaction cannot be sustained
without the introduction of external neutron sources. The
introduction of these sources offset losses due to diffusion
through the sides of the assembly and structural absorption. Such
an assembly can be obtained by simply surrounding a neutron source
with a moderator and enriched uranium. Subcritical assemblies are
smaller than critical assemblies and consequently requires less
material. Other advantages are that no complicated and expensive
control mechanisms and instrumentation systems are required. The
low flux levels at which these facilities operate present no
problems with respect to radiation protection. Subcritical
facilities are much cheaper to design and construct as compared to
nuclear reactors and are particularly attractive for training
reactor physics students. A typical schematic drawing of a subcri
tical assembly based on this principle is shown in Figure 2.3.
Some critical assemblies also use isotopic sources to drive them
in the presence of a multiplying medium. A typical type of such a
reactor using 20% enriched compressed mixture of U-,0o inJ o
Gauss-Seidal method.
2.3 Subcritical assemblies
%
24
University of Ghana          http://ugspace.ug.edu.gh
Fig.2.3 Schematic drawing of subcritical uranium assembly 
RS-155[26 j
(1) Active zone, containing granular U02 (U enriched to 36% in 
U235) mixed with polyethylene (660g U235);
(2) lead shield;
d
(3) shield of paraffin loaded with 5% boron carbide;
(4) water shield.
25.
University of Ghana          http://ugspace.ug.edu.gh
polyethylene powder is shown in Figure 2.4 and the neutron flux 
distribution within it is presented in Figure 2.5.
The degree of subcriticality of an assembly is measured by 
the value of the effective multiplication factor, [24] . The
effective multiplication factor is defined as the ratio of the 
number of neutrons at a given time to that of the preceding time 
[14-16]. This is mathematically expressed as
number of neutrons in one generation n(t)
K = =
number of neutrons m  the preceding generation n(t-l)
(2.45)
This parameter can also be used to define criticality and super- 
criticality. The assembly is said to be subcritical if the 
effective multiplication factor is less than unity and critical 
when it is unity. If the value is greater than unity then the 
assembly is supercritical.
The usefulness of subcritical assemblies cannot be over 
emphasized. They serve as research as well as teaching tools and 
can be used to predict the optimal arrangement, fuel proportion 
and spacing and for obtaining other important data relating to the 
design and construction of a proposed full-size operating reactor 
[24, 25] . Works by Kerten and Went [25] on subcritical aqueous
suspension reactors have shown that the concentration of fuel as 
well as the enrichment of the fissile components can be adjusted 
to the experimentally desired values. The PUK subcritical facility 
designed to have Keff equal to 0.95 [20, 25] is capable of
26
University of Ghana          http://ugspace.ug.edu.gh
1 Core 9 Reflector of Neutron Detector
2 Reflector (Graphite) 10 Experimentation Channel
3 Inner Reactor Vessel 11 Thermal Column
L Lead Shielding 12 Neutron Source
5 Water Shielding 13 Drive fo r Neutron Source
6 Control Plate - H Hoist fo r Reactor Core Half
7 Drive for Control Plate 15 Motor fo r Hoist
8 Neutron Detector
Figure 2.4 Reactor SUR [23]
27.
University of Ghana          http://ugspace.ug.edu.gh
Figure 2.5: Relative Neutron Flux D is tr ib u t io n  [23].
28.
University of Ghana          http://ugspace.ug.edu.gh
measuring the effect of water gaps in a lattice. It was used to 
study methods for preventing flux peaking and to supply 
experimental data for testing calculational techniques. 
Researchers utilized subcritical assemblies and other facilities 
to measure the effective multiplication factor of damaged cores 
such as TMI-2 and also to estimate its infinite multiplication 
factor distribution [25]. Further work [21, 23, 25] have also been 
carried out in connection with the use of subcritical assembly for 
specific reactor physics problems.
As stated above subcritical and critical assemblies played 
complimentary roles for nuclear reactors. They can be designed to 
find application in neutron activation analysis technique. Before 
presenting work in this area in Chapter Three, the principles 
underlying neutron activation analysis and the equipment presently 
in use in most scientific research laboratories are very briefly 
discussed. *
Neutron activation analysis was proposed since the advent of 
nuclear reactors . G. von Hevesy and H. Levy used the technique 
for the analysis of dysprosium in rare earths using an isotopic 
neutron source [2 6] . It has since developed and become a powerful 
nuclear analytical technique for multi-element analysis. The 
technique is very sensitive and more advantageous when the analysis 
involve trace elements and also when neutron flux of the order
29
University of Ghana          http://ugspace.ug.edu.gh
1012n/cm2-s is used. Moreover, the technique is non-destructive and 
hence the original sample may be recovered after analysis. Small 
quantities of the samples are required for analysis and the 
results are not influenced by the chemical state of the elements 
being investigated [26-28]. Additionally, the samples are not 
contaminated after irradiation which is of prime concern in trace 
elements analysis.
The high neutron fluxes produced by nuclear reactors are 
particularly useful for the determination of extremely small 
quantities of a large number of elements using neutron activation 
analysis. The technique is highly sensitive for trace elements in 
a variety of matrices in many fields of study such as:
(i) high purity materials
(ii) biological
(iii) environmental
(iv) geological
(v) forensic science etc [23, 26-29] .
When the sample is exposed to high flux of thermal neutrons, 
some elements in the sample absorb the neutrons, become unstable 
and disintegrate giving off characteristic rays. Each element 
gives off a unique characteristic rays. Thus the element can be 
identified and quantified if possible. This is done by using 
Si(Li), Ge(Li) or Nal detectors to detect the emitted rays. The 
signals from the detector are amplified in a spectroscopy 
amplifier and passed on to a multichannel analyzer which converts 
the analog signals to digital for analysis.
30
University of Ghana          http://ugspace.ug.edu.gh
2.4.1 Neutron Activation Analysis Equipment
Various experimental equipment using radioisotopes for neutron 
activation analysis are being used in many laboratories. A 
schematic drawing of one such facility is shown in Figure 2.6. The 
experimental equipment of the Ghana Atomic Energy Commission 
(G.A.E.C) is discussed here.
The G.A.E.C radioisotope experimental facility is a typical 
example of such facilities for NAA and is shown in Figure 2.7. It 
consists of a 2 0Ci Am/Be radioisotope neutron source, a pneumatic 
transfer system, a High Purity Germanium detector (HPGE), a 
Canberra spectroscopy amplifier, an oscilloscope, a Canberra 
series 3 0 Multichannel Analyzer (MCA), an APPLE II PLUS Micro­
computer and a printer.
The 2 0Ci Am/Be radioisotope source is cylindrical in shape 
and has a source strength of 105n/s. The source is centrally 
located in a fibre glass tank filled with de-ionized water (non­
multiplying medium) which serves as moderator. As an additional 
measure to prevent any radiation hazards through leakage from the 
source the fibre tank is surrounded with concrete blocks.
Samples are transferred to and from the source in rabbits by 
means of the pneumatic transfer system. It is operated at a pre­
ssure of 15psi and has a sample travel time of 3 seconds [30,31] . 
The HPGE detector is maintained at liquid nitrogen temperature 
(77°K) in a dewar flask to keep down dark currents. It is operated 
at a bias voltage of 3,000V and has a resolution of 2.25KeV for
31
University of Ghana          http://ugspace.ug.edu.gh
32.
University of Ghana          http://ugspace.ug.edu.gh
IR
R
A
D
IA
T
IO
N
 
E
N
D
S
Figure 2.7: PNEUMATIC TRANSFER SYSTEM
SHOWING TUBE CONNECTIONS[30,31]
University of Ghana          http://ugspace.ug.edu.gh
the 1332KeV photo peak of Co-SO. Characteristic rays from the 
irradiated sample are detected as electrical signals, amplified by 
the Canberra spectroscopy amplifier and processed in the Canberra 
series 30 multichannel analyzer. Since these electrical signals 
are analog, they are often converted to digital for easy analysis. 
An APPLE II PLUS micro-computer coupled to the MCA is used for 
analyzing the data.
2.4.2 Isotopic Sources
The major source of obtaining high thermal neutron flux for
activation analysis is from nuclear reactors. There are however
other neutron sources such as accelerators, generators, and
isotopic sources. The most commonly available are the isotopic
sources. The choice of a particular source is determined by
factors such as price, half-life, neutron yield, required
shielding from gamma rays, toxicity, small size and of course
availability [26, 28, 31] .
There are three main types of isotopic neutron sources namely:
alpha-emitters, spontaneous fission sources and gamma-emitters.
Alpha-emitters produce neutrons through an (a,n) reaction. These
radioactive alpha emitters can introduce (a,n) reactions on some
target nuclides. A typical reaction using beryllium as target
9 2.2
material is Be(a,n) C. Table 2 [2 6] shows a number of target
« o-i o
nuclides with their Q-values and neutron yields using Po
alpha-emitters as neutron source. Beryllium gives the highest Q-
34
University of Ghana          http://ugspace.ug.edu.gh
and neutron yield values and is generally preferred for neutron 
activation analysis purposes.
Target Nuclide Q-value Neutron Yield per 106 alphas
7Li -2 . 79 2 .6
9Be 5.70 80
ioB 1.06 13
lie 0.16 26
13C 2 .22 10
18o -0 . 70 29
19„p -1.95 12
In the case of spontaneous fission, some transuranium elements
disintegrate not only by the a-decay, but also by spontaneous
fission. This process releases several neutrons. The most
252frequently used of these transuraniums is Cf. It is small in
size, has similar neutron spectrum with fission uranium but has a
relatively short half-life of 2.65 years and also very expensive.
For gamma-emitters, i n , n ) reactions occur only when the incident
r-energy exceeds the binding energy of the neutron in the nucleus.
Most nuclei have very high binding energies, thus in practice only 
124Sb is used as a gamma source because of its low binding energy. 
The neutron flux energy obtained is very high but serious handling 
problems arise as the sr-dose rate is also quite high.
35
University of Ghana          http://ugspace.ug.edu.gh
2.5 Concluding Remarks
The reviewed research work as related to subcritical assem­
blies has revealed that most developing countries and small 
laboratories cannot fully explore the activation technique because 
they cannot afford nuclear reactors. They therefore depend on 
facilities which use radioisotopes in non-multiplying medium for 
the technique. The thermal neutron flux emanating from such 
facilities is low. With the low flux very small fractions of 
samples are normally activated. The result is that photo peaks are 
not so much different from the background. This presents a problem 
when trace elements are to be analyzed. Poor detection limits are 
therefore associated with these devices. Samples requiring high 
fluxes take longer time to irradiate. Another major problem 
associated with these facilities is that only one irradiation 
channel is available. Samples cannot therefore be irradiated simul­
taneously. This, coupled with the low flux makes the technique 
time consuming.
It will be of interest therefore to provide a conceptual
nuclear design of a multi-region subcritical assembly driven by
isotopic neutron sources in a multiplying medium (neutron
multiplier) for neutron activation analysis. The envisaged
multi-region assembly will be capable of providing thermal neutron
7 2flux as high as 1x10 n/cm -s. It is to be realized by surrounding 
radioisotope neutron source(s) with a multiplying medium. In this 
work the multiplying medium to be applied is U02 (20% enriched) .
36
University of Ghana          http://ugspace.ug.edu.gh
Multiple isotopic neutron sources will also be introduced to drive 
the multi-region assembly. This will increase the neutron flux and 
also provide for the possibility of creating more than one irradia­
tion site in different regions. This will permit simultaneous 
irradiation of samples. Attention will therefore be focused on the 
selection of a homogeneous fuel-moderator system, reflecting 
material, shielding materials and isotopic source strengths to
determine optimal conditions for the achievement of neutron flux
7 2of the order of 10 n/cm -s for neutron activation analysis and 
other reactor physics experiments.
Neutronic calculations for the design of the subcritical 
assembly will involve solving the one-dimensional (1-D) steady 
state multi-group neutron diffusion equation. The multi-groups will 
be collapsed to two energy groups, fast and thermal. A numerical 
scheme based on finite difference method is to be used to develop 
a computer code. The code will solve the 1-D multi-group neutron 
diffusion equation with the possibility of sources placed at 
different positions. It will be written in FORTRAN 77 programming 
language for IBM PC and will have the capability of investigating 
the following:
(i) choice of reflector
(ii) effect of multiple sources
(iii) choice of shielding materials.
37
University of Ghana          http://ugspace.ug.edu.gh
CHAPTER THREE
THEORY OF THE CODE SUNDES 
3.1 INTRODUCTION
In the concluding remarks of Chapter Two, it was stated that 
there is the need for the design of a multi-region subcritical 
assembly using multiple sources to achieve high fluxes for NAA. 
The neutronic design of nuclear reactors, subcritical assemblies 
and other nuclear devices must necessarily go with the ability 
to predict how neutrons are distributed throughout the system. 
It is necessary for one to solve a modified multi-group diffusion 
equation. Detailed calculations based on transport theory are 
very complicated and difficult to solve. The present analysis has 
therefore been restricted to the solution of two group, one­
dimensional diffusion equation.
In this chapter, the mathematical model, numerical methods 
and computational procedures used in the development of the SUNDES 
code which solves the neutron diffusion equation in the presence 
of isotopic sources is discussed.
3.2 MATHEMATICAL MODEL
For a steady state to be maintained in a subcritical assembly, 
external neutron sources must be introduced. The neutrons from the 
external sources, fission and scattering processes offset losses 
from the system through diffusion and structural absorption. A 
spectrum of neutrons, T)(r,E) is emitted by the sources and are
38
University of Ghana         http://ugspace.ug.edu.gh
elastically scattered in the assembly and slowed down, part of 
these are absorbed by the multiplying medium to cause fission. As 
a result of the fission process a new generation of neutron 
spectrum, ^(r,E) is produced. In the presence of the external 
(isotopic) sources, equation (2.6) is thus modified and is 
presented mathematically as
Since there is a wide range of energy groups in the assembly, 
equation (3.1) may be written in the form of multi-group diffusion 
equation for any energy group, g and space point, k in the form;
where <5 (r->r ) is the Kronecker delta function which is defined as
-D (r, E) V20 (r, E) + Xt (r,E)0(r,E) Es (r,E'->E)«(r,E)dE'
0
+ *(r,E')i>gEfg(r,E')0(r,E')dE'
J Q
(3.1)
G G
I * I Xg',kvg£fg' , k ^ g ' , k
g'=i
g'*g
g=l
N G
(3.2)
n=l g'=l
39
University of Ghana          http://ugspace.ug.edu.gh
,1 for r = r.
6 (r->r ) = n
n
0 for r * r
The differential operator for the leakage term may be written as
-D , V20 , = -D i -a fg/k rg,k g,k r dr
a  d d> , 
r drg (3.3)
where
a
I 0: for plane geometry
1: for cylindrical geometry
2: for spherical geometry
and r a variable, is defined for the respective geometries as;
slab thickness (a=o)
0 < r s Rc : radius of cylinder (a=l)
0 s r s R ; radius of sphere (a=2) s
Let the local variable r be divided into K+l intervals (not 
necessarily equally spaced) as depicted in Figure.3.1. Whence the 
following terms are defined;
rk-1/2 = 2 (rk+ rk-l> 
rk+l/2 = |(rk+ rk+1)
Ak- (rk -rk-l}
’fCf /|fet j k,j. .
■i-r
40
University of Ghana          http://ugspace.ug.edu.gh
Ak+ = (rk + r rk}
0 = 0 (a=0) 
v<p = o (a=1 ,2 )
k-1/2 k+1/2
I I
1 1 1 i |
0 1 2
i
k-1
1
k k+l
0=  0
K K+l
Figure 3.1: Mesh points with intervals
Assuming that the current is continuous at r^ and then integrating 
the differential equation (3.2) between the limits an<^
rk+l/2 t*ie leakage term is obtained as follows;
-Dg,k
'k+1/2 ,  
d 
dr
"k-1/2
r d<t> , (r, )
drg,k k
dr = D , ra, . d g,k+ k-1/2 ^g,k+(rk-l/2]
n a d 
4 g,k+rk+l/2 dr ^g,k+(rk+l/2]
(3.4)
where the quantities on the positive and negative sides of r^ are 
denoted by + and - respectively. From the terms defined earlier;
*g,k(rk ’ *g,k(rk-l/2 >
g g , k - lrk-l/2> ' -------------------------
Ak-
41
University of Ghana          http://ugspace.ug.edu.gh
and
g g , k + (rk+l/2)
^g,k(rk+l} - ^g,k{rk)
Ak+
(3.S)
Substituting equations (3.5) and (3.6) into (3.4) gives
rk+:■ /  2
g, k
d
dr r<X-^a k (rk} drg,K K
dr = D , rf _ /„ g,k+ k-1/2 ^g,k{rk) - ^g,k(rk-l}
'k-1/2
- D t
g,k+ k+l/2
Ak+
Ak-
V k ^ k + l 3 - ^g, k {rk) (3.7)
In the next integration the last term on the LHS and the first 
term on the RHS of equation (3.2) are combined, i.e.
‘k+l/2 G
'k-1/2
^tg,k<^ g,k(rk) ^  Esg'^g,k <j> (r )
g' =1 g '*■ K
g'*g
a .r dr
£tg,k-^g,k- ^ k'
G
I ^sg' -»g, k^g' ,k^rk^
'k-1/2
g'=i 
g' *g
r dr
42
University of Ghana          http://ugspace.ug.edu.gh
_rk+l/2
^tg,k+^g,k+ k [_ , ^sg'->g,k^g',
g'=l
g'*g
k+ (rk)
a , r dr
(3.8)
Simplifying equation (3.8) further it becomes
cx+l a+1
rk rk-1/2
G
a+1 ^tg,k-^g,k-(rk} " X !  ^sg'-»g,k-^g',k-(rk}
g'=i
g'*g
a+1 a+1
rk+l/2 rk
a+1 ^tg,k+^g,k+(rk} s^g'-*g, k*V , k+(rk}
g' =1  
g'-g
(3.9)
If the fission neutrons term and isotopic neutrons term in Eq.(3.2) 
are combined and integrated an expression of the form below is 
obtained;
r,
Sg',k(rk) + Qg',k(rk)
’k+1/2
a , r ar =
rk-l/2
Sg',k + Qg,k
a+1 a+1
rk+l/2 rk-l/2
a+1
where G
Sg',k{rk} ~ Y . X g',k(^ f }g',k^g'
(3.10)
(3.11)
g'=i
and N G
V . k  ‘ Z I " g ' . * q9 V
n=l g'=l
(3.12)
43
University of Ghana          http://ugspace.ug.edu.gh
Using equations (3.7) (3.10), the simplified integrated equation
becomes:
Pg,k-rk-l/2
A k -
^g,k(rk)- *g,k{rk-l> D T g,k+ k+l/2
Ak+
^g, k (rk+l} ” ^ g , ^
a+1 a+1
' k - 1 / 2
a+1
Z t g , k - * g , k - { r k ) z s^g'-»g, k-^g, k- (rk}
g'=i
g'*g
„a+l
rk+l/2
a+1
a
^ t g , k - <t>g , k + { r k ) I
g'=i
g'*g
s^g'-»g, k+^g, k+ (rk}
a+1 a+1
■k+l/2" k-1/2
a+1 Sg',k + Qg',k
(3.13]
To simplify work, denote the following by;
a
k " rk-l/2 
Crk‘rk-1)
a
pk " rk+l/2 
‘(rk+i-rk)
University of Ghana          http://ugspace.ug.edu.gh
, a+1 
^k = (rk
a+1 % 
rk-l/2
(a+1)
« - r-ra+1 
k _ k+l/2
a+1, 
rk )
(a+1)
? - (ra+1 ?k “ k+l/2
a+1 
rk-1/2
(a+1)
and
G
U, - I ^sg'->g,k- ^ g ^ k - ^ k 5
g'=l
g'*g
W, I
g'=i
g'*g
Esg'^g,k+0g',k+(rk}
thus substituting the above definitions into equation (3.13) gives
D , x, g,k- k ^g/k (rk} - ^g,k{rk-l} Dg,k+Pk ^g,k(rk+l) - ^g,k(rk}
'
+ Mk ^tg,k-^g,k-(rk} Uk + ^tg,k+^g,k+(rk5 Wk
= Sg',k + Qg' ,k (3.14)
Simplifying equation (3.14) further gives the expression
45
University of Ghana          http://ugspace.ug.edu.gh
-D
g,k-pk^g,k^rk+l) D 1 Mg,k- k Dg,k+Pk+ Mk^tg,k-+ Pk^tg,k+ ‘g,k(rk>
-Dg , k - V g , k <rk-1> - ?k(Sg',k + v . k ’ + “k ^  ^sg'^g, k-^g', k <rk>
g'=l
g'*g
G
+ Pk ^  Esg ' ^ g ,k + ^g '  , k (rk) (3.15)
g'=i 
g'*g
If it is still assumed that the material properties remained 
the same within the delta neighbourhood of r^ then the following 
approximations are made;
y = y = y d = d = dsg' g,k- ^sg' g,k+ ^sg' g,k ' g,k- g,k+ g,k
and ^tg,k- ~ ^tg,k+ " ^tg,k
Using the above approximations, equation (3.15) is rewritten as 
a finite difference equation in the form
■Dg,kpk^g,k(rk+1) +
"Dg,kTk^g,k ^rk-l^ ~ ^k
Dg,kTk + Dg,kPk + ^k^tg,k
G
g',k g' , k ^sg'->g,k^g',k
g'=i 
g'*g
* g»k(rk^
(3.16)
46
University of Ghana          http://ugspace.ug.edu.gh
For brevity and without loss of generality let us designate
^g,k(rk+l) = g, k+l ^g,k(rk^ g,k and ^g,k(rk- ) = g, k-1
then the finite difference equation (3.16) for energy group g, and 
position k, on the radius becomes a nonhomogeneous simultaneous 
linear equation of the form
-ak,k+i^k+i + ak,k^k - ak,k- A = i  = ^k (3-17)
The subscript g has been omitted for simplicity and the 
coefficients are given by
ak,k+l ng,k^k '
ak,k-l Dg,kTk
ak,k ak,k+l + ak,k-1 + ^k^tg,k
°g / k " ?k Sg',k + Qg',k + ^sg'->g, k^g',k
g' =1 
g'  *g
Collapsing the multi-group analysis to two (fast and thermal) 
groups and further assuming that neutrons can only scatter from 
higher energies to lower energies and that all neutrons are born 
at fast energy then
47
University of Ghana          http://ugspace.ug.edu.gh
and
crg,k
?k (Slk + Qlk ) for faSt ^rouP (9=1)
?kl 2,k ^ lk for thermal group (g=2)
(3 .18)
where sl]£ - £  Xg. >k lrZf >g . ,k*g. ,,
g=i
N G
Qg',k= X! Z  ^g''k^g',] 
n=l g'=l
To evaluate Q , , , the expression for q , , must be obtained, g , k g , k
q , , is dependent on the type of source and its location in the g / -K
assembly. For the purpose of this work, three types of sources 
are considered:
(a) a cylindrical source located at the center of the assembly
(b) concentration of point sources at the center and at any other 
position
(c) the combination of a cylindrical source at the center and 
placement of point sources at any other part.
For case (a), if the source strength is q then the 
expression for qlk everywhere in the assembly is of the form
48
University of Ghana          http://ugspace.ug.edu.gh
where = 2.405rk/rkl, while rkl and h are the radius and height 
of the assembly respectively. In the evaluation of qlk for the mesh 
points k, Bessel functions [6, 33, 34] for small and Y are used:
2 4 6
J (X, ) = 1  k + k k +
° 4 64 2304
(3 .2 0 )
Y, Y,3 Y,5T /v \ k k k
l l k' “ 2 16 + 384 (3.21)
For the case of point sources concentrated at the center of the 
assembly qlk is expressed simply as
dlk
0 rk > rx
(3.22)
and any other place apart from the center of the assembly as
, q r, = r ,^o, n k -j
£lk (3.23)
rj < rk rkz
where j and kz are the first and last mesh points of any 
particular zone.
Apart from the source conditions, boundary conditions are
University of Ghana          http://ugspace.ug.edu.gh
often applied in arriving at reasonable numerical solutions to 
differential equations. Equations of the type (3.17) have infinite 
number of solutions but every physical case has a single function 
representing the flux. In order to obtain such a function, some 
restrictions or boundary conditions must be imposed upon the 
general solution. These conditions are dictated by the nature of 
the problem. In the derivation of equation (3.17) it was assumed 
that (i) the current is continuous between rk j_/2 and rk+l/2' 
and (ii) the neutron flux must be finite and non-negative in the 
region r^ and zero at the outermost boundary. Apply the boundary 
condition that the flux is maximum at the center ie 70 = 0, the 
first equation of the series is given by
3g,lr3/2
^g,2 ^g,1
r - t  2 1
+  ^ a g , l ^ g , 1
a+1 a+l'
r3/2 - rl
a + 1
a+1
r2
a+l’
a + 1
G
I
g'=l
g '*g
Jsg'->g, 1 g,0 , , + S .,+Q ,,l Y ’ , 1  g,1 vg,l (3.24)
where the subscript 1 denotes the mesh point at center of the 
assembly and the coefficients a and a are respectively-L f ± ± f Z
defined as:
al,l Dg,lr3/2 + ^ag,1
a+1 a+l' 
r3/2 rl
a + 1
50
University of Ghana          http://ugspace.ug.edu.gh
and
t-, a+1
al,2 “ ~Dg,lr3/2
Also imposing condition (ii), the last equation of the series for 
the last mesh point is thus obtained as
Dg,klTkl + ^ag,kl^kl &g,kl “ TklDg,kl^g,kl-1
G
= Ukl ^  ESg'-»g,kl*g' ,kl + Sg,kl + Qg,kl
g '= i
g'*g
(3.25)
where the coefficients are fixed by virtue of the condition as
and
akl,kl-l “ Dg,klTkl
akl,kl ~ Dg,klTkl + ^ag,klMkl
Several other conditions can be imposed but are however ignored at 
this particular stage of development of the code SUNDES.
3.2 1 Matrix form of multi-group equations.
In order to simplify the equations, it is convenient to 
consider the matrix form of the multi-group equations. The equation
51
University of Ghana          http://ugspace.ug.edu.gh
(3.17) obtained for any group g and position k in the assembly can 
be transformed into a matrix equation of the form
G
Ag,k^g,k B , +g,k
g'=l
g ' *g
(3.26)
where
ail ai2
'
^g/i
a21 a22 a23
■o- 5? X*
ii
^g, 2
akl-l,kl-2 akl-l,kl-l akl-l,kl ^g# ki-i
i—1 1
i—1
j*i—1
M
akl,kl
c , ,g'->g, k
'kl-l
'kl
b.
g,k
kl-l
\l
The matrix C , , represents the scatter term whilst B , is forg'->g,k ^ g,k
the sum of fission and source terms. The solution to this matrix
equation provides the solution for the multi-group diffusion 
equation. The numerical method and techniques applied in its 
solution are presented below.
52
University of Ghana          http://ugspace.ug.edu.gh
3.2.2 Computational Procedure
A computer code named SUNDES; an acronym for Subcritical 
Nuclear DESign was developed to solve the matrix equation (3.26). 
The steps providing a suitable solution to the matrix equation as 
solved using the code are as follows;
(a) Guess the initial values of <p^  ^ and <p„ ^ Calculate the fission 
density, S and source density, Q.
(b) determine the coefficients of the matrix A ,1 , K
and °ik = Sl,k+ Ql,k
(c) Solve for new values of <j>^  ^ using Gauss-Seidal method
(d) Introduce ^obtained in (b) and determine the coefficients 
for the second group and solve for <p^  ^
(e) Calculate the approximate value for fission density, fd
(f) Introduce the values of ^ and ^ and repeat steps (b) 
through (d) until convergence is achieved in <p2 ^ and fd.
The flow chart of iterative procedure of the code is presented in 
Figure 3.2. First, the input data consisting of geometry
53
University of Ghana          http://ugspace.ug.edu.gh
University of Ghana          http://ugspace.ug.edu.gh
University of Ghana          http://ugspace.ug.edu.gh
Increase iteration number
Figure 3.2 Computational Flowchart for SUNDES
University of Ghana          http://ugspace.ug.edu.gh
parameters, macroscopic group constants etc. are read. There are 
three types of sources and are discussed in section 3.1. The 
required source is selected and the number of mesh points deter­
mined based on neutron free path. Mesh points along the radius of 
assembly are computed followed by setting the counter for the 
fluxes.
Subroutine PARA calculates the parameters /3^ , and while 
COEFF calculates the coefficients of the matrix A. The 
convergence criterion is applied to test for convergence and if 
that is achieved the normalized fluxes are written out and 
programme terminated. For non-convergence and depending on the 
type of source, there are three subroutines SUMP, DUMP and TUMP to 
calculate matrix B to satisfy source conditions. Subroutine SIMPR, 
determines the quantities of neutrons from fission. Another sub­
routine GSITER solves the matrix equation (3.26) and the flux is 
tested for convergence. If convergence is still not achieved, the 
energy groups are increased. The multiplication of matrices when 
the energy groups are increased is performed by a subroutine 
PMATRX and recycled through GSITER until convergence is achieved.
For an outer iteration loop to begin, the previous outer 
iteration fluxes are examined and the fission density computed 
from the equation
The computation of effective multiplication factor will not be of
G
(3.27).
g=i
57
University of Ghana          http://ugspace.ug.edu.gh
concern here since a fixed source problem is being considered. The 
group fluxes are rather simply allowed to find their own levels 
until convergence is achieved, i.e.
fd(l) fd(l 1}
fd (i)
(3.28)
where the limit of e has been chosen to be as small as 10" .
3.2.3 Inner iteration
The neutron diffusion equation represented by equation (3.26) 
is solved in the inner loop for any group g. Gauss-Seidal method
[6] is used to solve the nonhomogeneous linear equation. A general 
expression for the solution of the simultaneous equations is given
♦ i
(m)
a. .li
i-1 i-1 n
’i - i
i=i
a . . <b . 
iD 3
(m) a . ,<j). (m-1)
j=i+l j=i+l
(3.29)
where the subscripts denote the iteration numbers.
The test for convergence for the iteration method is by 
comparing the changes in the <p^ values between successive 
iterations ie
58
University of Ghana          http://ugspace.ug.edu.gh
where Stp. = A - <* .t"m and e„ is a tolerance.*1 ri l 2
There are several acceleration schemes to increase the 
convergence, such as coarse mesh rebalance, Chebyshev, synthetic 
method and overrelaxation [35]. For this work, the rate of conver­
gence is improved by applying the technique of overrelaxation 
which is easily programmable. Equation (3.29) can be used to
( TT\ \
estimate new values of < p / which, provided the process is
convergent, are closer to the required solutions than the previous
<p 1 ^ . Estimation of the new values of <p ■ ^  becomes i ri
and w is the overrelaxation factor. Varga [36] has exhaustively 
discussed the theory of predicting its optimum value. It is found 
to be within a range 1 < u < 2 . A n  optimum value of 1.5 was
(3.31)
f n
where (3.32)
_ g
chosen to provide convergence to a tolerance limit of 10 for 20 
iterations in this work.
University of Ghana          http://ugspace.ug.edu.gh
3.2.4 Normalization of Fluxes
The code SUNDES provides for the normalization of neutron 
fluxes to desired power levels in the assembly. The power equation 
is given by
r G
P X  = M  )  ^g,k^fg,dVk 
J V„g=l
(3 .33)
where V^ is the volume of the assembly, £ = 3.2 x 10 11watts/sec 
per fission and
dV, =
1
2Trr
4rrr
a=Q (slab) 
od=l (cylinder) 
a =2 (sphere)
If the fission cross section is considered a constant in each of 
the multi-regions 1 of total number of regions, nl, then the power 
equation (3.33) may be rewritten as
G nl
V  £  I  £fg,l?g,l4v? 
g=ll=l
(3.34)
where AV-, =
(R1+1-El)HL
n(Rl+l'RlIH
l lRl+l-Rl>
for a=0 (slab) 
for a=l (cylinder) 
for a =2 (sphere)
For a  = 0, H and L are the transverse dimensions of the slab. For 
a =1, H is the height of the cylinder.
60
University of Ghana          http://ugspace.ug.edu.gh
Having solved the matrix equation (3.26), the fluxes should
be normalized to the total power of the subcritical assembly. This
requires a normalizing constant f  such that <p , = tp<p , andg / k y , k
relates the power through the equation
P
ip = G
V. Z Tr: 1 (I) , dv“fg,kyg,k k
g=i
X
(3.35)
3.2 5 Input Description
The input data for the code SUNDES consists of two sections
(i) the geometry and convergence criterion
(ii) nuclear properties of macroscopic group constants, buckling 
and fractions of fission and isotopic neutrons in energy 
group g.
SUNDES. INP file is defined on UNIT 1 through which all the 
input data enter the main program SUNDES.FOR (see Appendix A) . 
SUNDES.INP is written according to a specified format and in the 
order listed below.
61
University of Ghana          http://ugspace.ug.edu.gh
Section (i)
Card Format Parameter Mathematical
symbol
Description
i4
i4
i4
i4
ig
ng
nl
nmax
f8.5 erl
G
nl
nmax
f 8.5
geometrical factor 
(o£=0,1, 2) 
total number of 
energy groups 
total number of 
homogeneous zones 
maximum number of 
iterations 
convergence crit­
erion for fission 
density
input power,watts
Section (ii)
The macroscopic group constants for the various zones starting 
from the center of the assembly for energy group, g denoted kg is 
as follows
62
University of Ghana          http://ugspace.ug.edu.gh
Card Format Parameter Mathematical Description 
Symbol______ ___________
nl (f8.3) sigf(kg,l) Vg ^ f g (r>E) macroscopic fission
cross section
nl(f8.3) siga(kg,l) V (r,E) macroscopic absorp-ag
tion cross section
10 nl(f8.3) d(kg,l) D^(r,E) diffusion constant
11 nl(f8.3) sigs(kg,l) V (r,E'-»E) macroscopic scatter-sg
ing cross section
12 nl(f8.3) b(kg,l) B (r,E) transverse buckling
y
13 nl(f8.3) q(kg,l) q , quantity of neutrons
g
from isotopic sources
14 nl(f8.3) fq(kg,l) x (r,E) fraction of fission
neutrons
15 nl(f8.3) fq(kg,l) tj (r,E) fraction of isotopic
neutrons
3.2.6 Output Description
Output from the code SUNDES.FOR can be obtained on screen or 
on printed paper. The output file, SUNDES.OUT is written in UNIT 6 
The information displayed by the code is in the following order: 
first the name of the code then the geometrical characteristics of 
the assembly and other important input parameters. Group constants 
of all the zones are also printed in this order: fission
63
University of Ghana          http://ugspace.ug.edu.gh
cross section, i^ Tf absorption cross section, £a scattering cross 
section, Es1_j,2 diffusion constant, D transverse buckling B, fission 
spectrum, x isotopic spectrum, ij and the quantity of neutrons from 
the isotopic sources, q.
The other portion of the output consists of the iteration 
number, i and the fission density, fd. Fast and thermal flux 
distributions at the various mesh points starting from the center 
of the assembly are also printed.
3.3 Concluding remarks.
The computer code SUNDES, which has been developed based on 
the mathematical model discussed above will be used for the design 
of the neutron multiplier. Preliminary test cases will be conducted 
using the code. These test cases ultimately enable a complete 
realization of the design of the neutron multiplier.
64
University of Ghana          http://ugspace.ug.edu.gh
CHAPTER FOUR 
NUCLEAR DESIGN OF THE NEUTRON MULTIPLIER
4.1 INTRODUCTION
Feasibility studies are normally carried out before embarking 
on detailed construction. Computer codes are employed to perform 
neutronic, thermal hydraulics, structural calculations etc. 
Additionally, various processes associated with the functioning of 
the nuclear equipment are thoroughly investigated. The control rod 
worth, fuel burn up, dose rate, radiation shielding etc are 
studied by the use of computer codes. The trend of results and 
values of neutron fluxes are normally used to select the geometry 
of the nuclear device.
In this chapter, the SUNDES code is tested by using nuclear 
data generated from WIMSPC for single and two region problems. 
This will be followed by its application to multi-region problems 
to obtain maximum thermal neutron fluxes generated by the neutron 
multiplier. A description of the mechanical aspects and operating 
characteristics of the device is also presented.
4.2 APPLICATION OF THE CODE.
4.2.1 Single and Two-Region Problem.
Geometrical features and nuclear properties of macroscopic
65
University of Ghana          http://ugspace.ug.edu.gh
group constants are needed in neutronic modelling. For the 
verification of the code SUNDES, macroscopic group constants
(Dl'D2'^al'^a2'^fl'^f2 and ^sl 2} were generated using WIMSPC 
[37], a PC version of Winfrith Improved Multi-group Scheme (WIMS).
[38-40]. WIMS is a general lattice code based on transport theory.
It contains a library of elements with specific identification
numbers from which elements of the assembly are selected. Data
cards are available through which geometry specifications and
material compositions enter the code. The programme then computes
macroscopic cross sections for each homogeneous region.
Macroscopic cross sections were generated for a single
homogeneous region of 2 0% enriched U02 and Be moderator. With an
isotopic source located at the center and/or at different
positions of the region the neutron flux distribution was studied.
Another case study was a two region problem which consists
basically of the initial single homogeneous region and then a
water region as reflector.
For the one homogeneous region with the presence of a source 
9of strength 8.45 x 10 n/s at the center, the flux distributions of 
both fast and thermal neutron is as shown in Figure 4.1. As the 
position of the source is shifted to 6cm and 12cm away from the 
center, the trend of the neutron flux is observed to change and is 
presented in Figure 4.2. The effect of maintaining point sources 
of equal strengths at the center, followed by one at 9cm and then
66
University of Ghana          http://ugspace.ug.edu.gh
Ne
ut
ro
n 
Fl
ux
 
(n
/c
Figure 4.1: Variation of Neutron Fluxes in a
Homogeneous Mixture of 20% U0o + Be.
67 .
University of Ghana          http://ugspace.ug.edu.gh
W  &C>CH>t>&Ot>t>ft'Q0PD&C>C>l>0
O'
■60''
O'90S'
66-
University of Ghana          http://ugspace.ug.edu.gh
another at 15cm is presented in Figure 4.3. Figure 4.4 represents 
the behaviour of the neutron fluxes when equal strength of sources 
are kept at 0, 6, 12, and 18cm from the center of the homogeneous 
mixture. It can be seen from the plots that the neutron fluxes are 
influenced by the position of the source. In the two region 
problem where water (H.-,0) is used as a reflector there is 
thermalization of neutrons in this region. In the thermalization 
region the collision takes place effectively with the H„,0 molecules 
and some neutrons are scattered back into the fuel region. This 
results in a peaking in the thermal flux. This observation as seen 
in Figure 4.5 is consistent with reactor physics predictions for a 
fuel region surrounded by a reflector. For practical purposes, a 
nuclear device will have multi-regions of fuel, moderator, 
coolant, reflector and shield. In the next section the code will 
be applied to various multi-regions with the ultimate aim of 
achieving high flux of thermal neutrons in the neutron multiplier.
4.2.2 Application to Multi-regions.
As stated, for the complete realization of the aim of this 
work, multi-region problems were considered. A multi-region 
consisting of nine (9) distinct homogeneous regions was adopted 
for the conceptual nuclear design of the neutron multiplier. The 
nine concentric radial rings are hereby denoted as A,B,C,D,E,F,G,H
69
University of Ghana          http://ugspace.ug.edu.gh
"\ V H - W  W  \ VVMr^HHHHHH-^
b fe & (VS-fe-D  & O & C O D O  &-& -ft~6-t> ft fe fe & 0 & O fe-fe-fe-p
University of Ghana          http://ugspace.ug.edu.gh
Th
er
ma
l 
Fl
ux
 
(
n
/
c
m
2 
-s
)
1Q51
H
1 0
_ 1 | | | | - |
0 5 10 15  2 0  2 5  .30
Radius (cm)
Figure 4.4: Effect of Source Located at 1, 6, 12
and 18cm on Thermal Flux.
71'.
University of Ghana          http://ugspace.ug.edu.gh
Ne
ut
ro
n 
Fl
ux
 
(n
/c
Radius (cm)
Figure 4 .5 : V a r i a t i o n  o f  N e u t r o n  F l u x e s  in a Two-Reg ion  
H o m o a e n e o u s ' :  M i x t u r e  o f  20% uo 2 +  Be.
72 .
University of Ghana          http://ugspace.ug.edu.gh
and I starting from the center. A horizontal cross section of the 
geometry is shown in Figure 4,6.
The first concentric ring A is occupied with a homogeneous 
mixture of fuel and moderator (20% enriched UC>2 + Be) . Regions B,D 
and F consist of AI material as cladding which separate regions A,
C and E from each other. Region C is occupied by a moderator. 
Again in region E there is a homogeneous fuel-moderator mixture to 
cause more fission in the system. Regions G and H constitute the 
reflector and shielding materials respectively. Finally, region I 
is ordinary concrete to act as biological shield.
A study was carried out to determine the choice of a 
reflector, shield and source strength to optimize neutron flux 
yield. The effect of different solid reflectors such as Be, BeO, C 
and liquid reflector (H„0) in region G on neutron fluxes was 
thoroughly investigated. Macroscopic cross sections are first 
generated for the various homogeneous zones of the multi-region 
problem. A typical input data for WIMS for this geometry is listed 
in Table 4.1. These two-group constants computed by the lattice 
code form part of the input data for SUNDES as listed in Table 4.2. 
A comparison of the variation of thermal flux with radius of 
assembly for the four different reflectors is presented in Figure 
4.7. It is observed that Be as reflector produced the highest flux 
though not very much different in value from that of BeO. The 
closeness in their values could be attributed to the fact that they
73
University of Ghana          http://ugspace.ug.edu.gh
A -----  UC>2 + Moderator (Be)
B -----  A1 Cladding
C -----  C Reflector
D -----  D A1 Cladding
E -----  UO^ + Moderator (Be)
F ---- —  A1 Cladding
G -----  Reflector
H -----  A1 Shield
I -----  Concrete
IC ----  Irradiation Channels
Figure 4.6: Horizontal Cross Section of Neutron Multiplier.
74 .
University of Ghana          http://ugspace.ug.edu.gh
CELL 7
SEQUENCE 1 i,
M E S H  21 
NGROUP 28 
NMATERIAL 5 
NREGION 9 2 
PREOUT 
INITIATE
* SUPERCELL CALN FOR SOURCE + BeO (SURO) *
FEWGROUP 3 5 9 14 21 25 26 27 28 29 30 32 33 35 $
38 40 41 43 45 47 49 52 55 58 60 63 66 69
MESH 2 1 3 1 4 1 4 1 4
MATERIAL 1 -1 303.0 1 235.4 1.39889E-05 2238.4 2.657891E-4 $ 
16 5.315782E-04 9 1.231023E-01 
MATERIAL 2 2.7 303.0 2 S
27 0.9947 29 0.16E-2 55 0.01E-2 11 0.0001E-2 56 0.32E-2 $
58 0.03E-2 7 6.0E-6 63 0.C12E-2 112 1.0E-6 
MATERIAL 3 1 303.0 3 20C1 0 111111 16 0.888889 
MATERIAL 4 2.95 303.0 2 9 0.3600 16 0.64
MATERIAL 5 2.3 303.0 2 56 .014 2001 .01 16 .531 29 .35 S
23 .016 27 0.034 12 .001
ANNULUS 1 3 . 0 1
ANNULUS 2 4 . 0 2
ANNULUS 3 9.0 4
ANNULUS 4 10.0 2
ANNULUS 5 32.0 1
ANNULUS 6 33. 0 2
ANNULUS 7 40.0 4
ANNULUS 8 41.0 2 •!
ANNULUS 9 43.0 5 ,
ARRAY 1 <1 4 6.0 0.0)
ARRAY 2 (1 6 36.5 0.0)
RODSUB 1 1 1.30 3 
RODSUB 1 2 1.85 2 
RODSUB 2 1 1.30 3 
RODSUB 2 2 1.85 2 
*BELL 1.15581 
POWERC 0 0 
BEG INC 
THERMAL 10 
REGION 11 11 
LEAKAGE 7
BUCKLING 1.33445E-03 
BEGINC
Table 4.1: WIMS Input Data for (20% enriched U0?_+Be)
75 .
University of Ghana          http://ugspace.ug.edu.gh
2
1
22 . 00 
0 . 2 
2 
9 
14 
68 
0 . 20 
20
. 0000 
3.00000 
.01549 
.01089 
.00869 
.00592 
.00238 
.43375 
.53685 
.00134 
.00134 
.00000 
.00000 
0.826+14 
.77e+14 
8. 45e+09 
. 00e+00 ■!
1. 0(?i000 : 
1. 00000 
1 .00000 
0 . 00000
1
1 .00000  
.00000 
.00000 
.01160 
.00816 
.00018 
3.48285 
3.36493 
0.00134 
0.00134 
0.00000 
'0,00000 
.82e+00 
.77e+00 
.00e+00 
00e +00 
1.00000 
1 .00000 
1 .00000 
0.00000
5.00000 
0.00000 
0.00000 
0.01825 
0.01340 
0.16572 
0.15976 
0.35474 
0.00134 
0.00134 
0.00000 
0.00000 
.82e+00
77e+00 
,00e+00 
.00e+00
1.00000 
1.00000 
1.00000 
0.00000
0000  2 
000? 
0000 
0112 
0055 
0005 
4989 
3572 
00134 
00134 
0.00000 
0.00000 
.82e+00 
.77e+00 
.00e+00 
.00e+00 
1 .00000 
1.00000 
1 .00000 
0.00000
00000 
0.01413 
0.00619 
0.00795 
0.00381 
0.37762 
0.43862 
0.56610 
0.00134 
0 00134 
0.00000 
0.00000 
.82e+00 
.77e+00 
.00e+00 
.00e+00 
'1. 00000 
1.00000 
1.00000 
0.00000
1 .
0.00000 
0,00000 
0.01139 
0.00746 
0.00057 
3.49181 
3.36958 
0,00134 
0.00134 
0.00000 
0 . 00000 
.'82e+00 
.77e+00 
.00e+00 
. 00e+i00 
1.00000 
1. 0 0 0 0 0  
1.00000 
0.00000
00000 
00000 
01846 
01469 
37776 
15732 
31216 
0.00134 
0.00134 
0.00000 
0.00000 
.82e+00 
.77e+00 
.00e+00 
.00e+00 
1.00000 
1.00000 
1.00000 
0.00000
1,00000 
0.00000 
0.00000 
0.01179 
0.00986 
0.00014 
3.47456 
3.42266 
0.00134 
0.00134 
0.00000 
0.00000 
.82e+00 
.77e+00 
.00e+00 
.00e+00 
1.00000  
1.00000 
1.00000  
0.00000
2.00000 
0 . ;-i00O 
0 .00003 
0.00 796 
0. 005*: 3 
0.01320 
0.48911
0.64401 
0.00134 
0.00134 
0.00000
.82e+00 
77e+®0 
.00e+00 
00e+00 
1 .00000 
1.00000 
1.00000 
0.00000
Table 4 . 2 :  SUNDES.INP Input Data for Multigion Problem
76 .
University of Ghana          http://ugspace.ug.edu.gh
N
eu
tr
on
 
Fl
ux
 
(n
/
Radius (cm)
Figure 4.7: Comparison of Thermal .Neutron Fluxes
For Four Different Reflectors.
77 .
University of Ghana          http://ugspace.ug.edu.gh
both have the same microscopic absorption cross sections but Be 
has a slightly higher microscopic scattering cross section as in 
Table 4.3. From the plot, it can be observed that HgO has the 
lowest value of flux in the central region A, but rather tends to 
have the highest values from region C. This is not surprising 
because H^O as reflector has the highest microscopic scattering 
cross section as compared with the rest of the tested reflectors. 
By scattering collision, more neutrons are returned into the fuel 
region thus producing additional fissions. Although water is the 
cheapest among the tested materials and it produced the highest 
flux, it was not selected as reflector for region G for the present 
analysis in order to avoid leakage and possible contamination. 
Graphite produced the least flux and could therefore not be a 
possible choice as reflector. The choice of either Be or BeO as 
reflector was not too obvious as a result of their closeness in 
flux. An advantage of using either Be or BeO is that reactivity 
will increase with the presence of fast neutrons available from 
(n, 2n) and photo neutrons from (y, n) reactions. However, BeO was 
preferred to Be because it is much less expensive and easier to 
fabricate.
Investigations were also conducted to select shielding 
material among aluminium (AI), lead (Pb) and stainless steel (SS) 
for region H. It was observed that Pb was the most efficient, 
followed by SS and then Al as seen in Figure 4.8. However, the
78
University of Ghana          http://ugspace.ug.edu.gh
Element or 
Molecule
Microscopic
Absorption
Crossection (barns)
Microscopic 
Scattering 
Crossection (barns)
Be 0 .0 0 9 5 7.0
BeO 0 .0 095 6.8
C 0 .0 0 3 4 4.8
H 0  
2
0 .6 64 103
Table 4.3: Nuclear Properties of Be, BeO, C and H 0
University of Ghana          http://ugspace.ug.edu.gh
Ne
ut
ro
n 
Fl
ux
 
(n
/c
m
 
- 
S)
Radius (cm)
Figure 4.8: Effect of AI, SS, and Pb Shields in
Region H on Thermal Flux.
80.
University of Ghana          http://ugspace.ug.edu.gh
differences in neutron attenuation of these materials were not 
significant. A1 was therefore selected because it is light in 
weight and less expensive. Its choice is also consistent with the 
low power of two milliwatts produced in the multiplier since it 
has found wide application in the construction of low power 
research reactors with less effect of radiation damage.
Multiple source effect was investigated by placing sources of 
different strengths in various zones of the neutron multiplier. 
The effect of the presence of source in the fuel-moderator region 
E can be seen in Figure 4.9. Its effect in region G is illustrated 
in Figure 4.10. The effect of placing sources in both regions E and 
G was also studied and the trend is shown in Figure 4.11. It was 
observed that the fluxes increased generally in the assembly 
and became more pronounced in the regions where the sources were 
located. Different source strengths have different effects on the 
values of thermal fluxes produced in the regions E and G. This can 
be seen in Figures 4.12 and 4.13 respectively. The two plots 
reveal that the values of the thermal fluxes achieved in the 
neutron multiplier depended on the source strength.
Based on the results obtained, it was decided that the 
configuration with BeO in region G and A1 as shield in region H 
could be the accepted geometry of the assembly. Figure 4.14 
presents the variation of neutron flux (fast and thermal) with 
radius of the assembly.
81
University of Ghana          http://ugspace.ug.edu.gh
M
eu
tr
nn
 
Fl
ux
 
(n
/c
rn
Radius ( cm )
Figure 4.9: Effect of Source Strength on Thermal 
Neutron Flux in Region E.
82.
University of Ghana          http://ugspace.ug.edu.gh
o 
\
X3
o
-t-JD
10 3020
Radius ( c m )
Figure 4.10: Effect of Source Strength on Thermal
Neutron Flux in Region G.
40
83.
University of Ghana          http://ugspace.ug.edu.gh
Ne
ut
ro
n 
Fl
ux
(J
c
] 0 5 =
1 0 4 =
%
\
>
v
Source Strength =  8 ,45  E- i -09n /
i
3010 20
Radius (2m )
Figure 4.11: Effect of Source Strengths on Thermal
Neutron Flux.in Regions E and G.
University of Ghana          http://ugspace.ug.edu.gh
Ne
ut
ro
n 
Fl
ux
 
(n
/
c
m
Figure 4.12: Effect of Source Strengths on Thermal
Flux in Region E.
85.
University of Ghana          http://ugspace.ug.edu.gh
N
eu
tr
on
 
Fl
ux
 
(n
/
Radius (cm)
Figure 4.1-3: Effect of Source Strengths on Thermal
Flux in Region G.
86.
University of Ghana          http://ugspace.ug.edu.gh
Radius (cm 1 
Figure 4.14: Variation of Neutron Fluxes with Radius.
87 .
University of Ghana          http://ugspace.ug.edu.gh
4.2.3 Description of Accepted Geometry.
The selected assembly which is cylindrical in shape has
height 86cm and radius 43cm; the height to diameter ratio being
9unity. A radioisotope source of strength 8.45 x 10 n/s is embedded 
in the homogeneous mixture of fuel and moderator. Region B is 
constructed of AI cladding of thickness 1cm which separates the 
fuel-moderator mixture in A from a 5cm water region, C. Region C 
contains six inner irradiation channels each of diameter 26mm. A 
1cm AI cladding in region D also separates the water in C from 
another fuel-moderator mixture of thickness 22cm. Region G is BeO 
reflector of thickness 7cm. This is separated from the fuel- 
moderator mixture in E by yet another 1cm thick AI cladding in 
region F. The BeO reflector has six outer irradiation channels 
each of radius 35mm for samples irradiation at different flux. The 
last two regions, H and I are 1cm Al cladding and 2cm thick 
concrete respectively which serve as shield for preventing leakage 
of neutrons from the assembly.
7 2Thermal neutron fluxes greater than 1 x 10 n/cm -s in regions 
of irradiation which is the main objective of the study were 
achieved in the assembly.
88
University of Ghana          http://ugspace.ug.edu.gh
4.3 DESCRIPTION OF MECHANICAL FEATURES OF THE NEUTRON MULTIPLIER.
The mechanical aspect of the multiplier is basically a 
lifting device for introducing or withdrawing neutron sources in 
and out of the multiplier. The vertical cross section is shown in 
Figure 4.15 . The mechanical engineering design features are
presented in reference [42] . For this work, only its technical 
features and operational characteristics are presented.
Four supports for the multiplier (1) bear it from underneath. 
The supports are bolted onto a concrete tank (2) which contains 
water into which the neutron sources (3) are submerged when 
outside the multiplier. The tank with water also serves as a shield 
against the radiation from the neutron sources and r-emitting 
sources from fission fragments when out of the multiplier.
The sources are mounted on thin metal rods (4) which are 
fitted vertically onto a metal plate (5) . The rods and plate 
materials are made of A1. At such low powers of operation it is 
expected that Al will be resistant to heat, corrosion and 
radiation while ensuring adequate strength for reliability.
The plate is bolted rigidly to an arm (6) which is also 
connected rigidly to an extension of a movable nut (7) . A screw 
(8) runs through the nut and is supported in thrust bearing (9) at 
the upper end. The screw is keyed to a spur gear (10) at its upper 
end. The spur gear is meshed with a spur gear pinion (11) which
89
University of Ghana          http://ugspace.ug.edu.gh
1. Supports 8. Screw
2. Concrete tank 9. Thrust bearing
3. Neutron source(s) 10. Spur gear
4. Thin Metal rods 11. Spur gear pinion
5. Metal plate 12. Electric motor
6. Arm 13. Steel pillar
7. Movable nut
Figure 4.15= Vertical Cross Section of Neutron Multiplier
University of Ghana          http://ugspace.ug.edu.gh
is keyed to an electric motor (12) mounted against a vertically- 
supported cylindrical steel pillar (13). The pillar also serves as 
a guide for the nut and provides a counter moment against that set 
up by the weight of the plate and arm, tending to bend the screw 
inwards towards the pillar.
A projection on the movable nut actuates an electrical switch 
at the two extreme ends of its movements up and down the screw. 
This switch is one of two alternative switches for switching the 
motor on or off. The second switch is located in the control room. 
When one switch is on the other must be off in order to close the 
circuit and vice versa.
The switch in the control room also has a provision to change 
the polarity of the motor. At the same time this switch is kept 
'on'or 'off', it reverses the polarity of the motor. The reversal 
of the polarity of the motor is necessary to change the direction 
of the rotation of the motor and hence to raise or lower the nut 
for the purpose of introducing the sources into or out of the 
multiplier.
When the sources are in their lowest position in the water, 
the first switch would have been in the 'off' position. At the same 
time the second switch would be off. To raise the sources into the 
the multiplier, the second switch is kept in the 'on/up' position 
with the polarity of the motor automatically reversed. The motor 
is thus started and the direction of its rotation is such as would
University of Ghana          http://ugspace.ug.edu.gh
raise the nut and consequently the sources upwards. Near its 
highest position 'on' its travel upwards, the projection on the nut 
would actuate the first switch to put it in the on position. This 
will open the circuit and thus put off the motor. The nut will 
remain at its highest position and will not overhaul downwards by 
virtue of its design.
In lowering the sources into the water, the second switch in 
the control room is put off and the polarity of the motor would be 
automatically reversed. The motor is thus started and the 
direction of its rotation is such as would lower the nut together 
with the sources. Near its lowest position, when the sources would 
have been totally immersed in the water, the projection on the nut 
will actuate the first switch putting it off. This will open the 
circuit and cut off power supply to the motor. The nut together 
with the sources would then remain in their lowest positions.
In effect, the movement of the sources is simply controlled 
in the control room by the second switch by putting it in the 
'on/up' or 'off/down' position as the situation demands.
92
University of Ghana          http://ugspace.ug.edu.gh
CHAPTER FIVE.
CONCLUSIONS AND RECOMMENDATIONS.
With reference to the primary objectives of this work as 
outlined in Chapter Two, a finite difference scheme was used to 
obtain a numerical solution to the one-dimensional neutron 
diffusion equation. A computer programme dubbed 'SUNDES' written 
in FORTRAN 77 programming language for an IBM.PC was developed 
based on the numerical scheme. With the aid of the computer 
package an extensive study was carried out to present a nuclear 
design of the neutron multiplier.
The trend of neutron flux distribution was found to conform
to a large extent to reactor physics predictions. It can be said
that with an isotopic neutron source placed at the center of the
assembly and with the appropriate materials used in the various
regions as presented in the design description, thermal neutron
7 2fluxes higher than 1 x 10 n/cm -s were obtained in the assembly for 
NAA and other reactor physics experiments. This is the main
objective of the study and this clearly shows that it is techni­
cally feasible to achieve such levels of fluxes in a neutron
multiplier.
The results revealed that different fluxes are produced
depending on the source strength. It is therefore very satisfying 
to mention here that since different geological and biological
93
University of Ghana          http://ugspace.ug.edu.gh
samples require different fluxes, this design would be appropriate 
The desired fluxes were obtained by introducing the appropriate 
source strength. With such high fluxes obtained in the neutron 
multiplier, the efficiency of NAA can be improved. Samples will be 
sufficiently exposed to the neutron flux thus increasing the 
probability of excitation and hence good detection limits. The 
design has provided for simultaneous irradiation of samples at the 
same flux and also at different fluxes in different regions. This 
will ultimately mitigate, the time consuming nature of the medieval 
radioisotope source in a non-multiplying medium used for NAA.
The mechanical design of the neutron multiplier is simple, 
affording the possibility of manufacturing most of the components 
locally. It will be easy to install and operate. Exposure of 
personnel and/or the general public to radiation is minimized by 
the shielding and cement materials. Another advantage of the 
design is the inherent ability of the screw not to overhaul. This 
makes it possible to maintain the sources at rest anywhere between 
and at the extreme ends of its travel.
From the foregoing results and discussions, a conceptual 
nuclear design of a subcritical assembly driven by isotopic 
neutron sources (neutron multiplier) has been achieved with great 
success. It may be recommended however that, for actual 
realization of the design, construction and installation of the 
subcritical assembly the following should be considered;
94
University of Ghana          http://ugspace.ug.edu.gh
(1) The neutron multiplier is compact and the possibility exists 
for the leakage of neutrons in the axial direction. A detailed 
nuclear design in two-dimensions (r, z) will be required to 
correctly account for the leakage. In case there is much 
leakage, then reflectors would be needed for the bottom and 
top of the assembly. The author is not aware of the presence 
of any computer package available in literature developed for 
subcritical reactor physics analysis employing multiple sources 
to drive such assemblies. There is therefore the need to 
develop another version of SUNDES for this aspect of neutronic 
design.
(2) A comprehensive thermal and striictural analysis be carried out 
for the assembly.
(3) In addition to a remote control device to be incorporated in 
the mechanical design, a design of irradiation systems for 
transfer of rabbit capsules for activation analysis must be 
provided.
(4) Instrumentation to detect parameters such as temperature, 
water level in the tank and neutron flux levels within the 
regions of interest such as near irradiation sites are needed. 
Protection limits for flux and temperatures shoiild be incorpo­
rated in the design so that the system could be shutdown if 
limits are exceeded.
(5) Radiological consequences must be determined to ascertain the
95
University of Ghana          http://ugspace.ug.edu.gh
safety of the neutron multiplier after it has aged. Procedures 
for decommissioning and waste disposal must be studied and 
subjected to approval by the competent national authority.
(6) Finally, an exercise to determine an economic feasibility 
with respect to cost of fuel, moderator, reflector, structural 
materials and operating cost must be carried out and compared 
with the cost of low power research reactors.
96
University of Ghana          http://ugspace.ug.edu.gh
REFERENCES
1 Wang, K. , Neutron Activation Technique with MNSR, China.
Institute of Atomic Energy (1993).
2 Melville Clark, Jr. and Kent, F. H. , Numerical Methods of
Reactor Analysis, Academic Press Inc. London (1964) .
3 Nakamura, S., Computational Methods in Engineering and Science
with Applications to Fluid Dynamics and Nuclear Systems, John 
Wiley and Son, Inc. (1977)
4 Bathe, K. and Wilson, E. L., Numerical Methods in Finite
Element Analysis, Prenctice - Hall (1976) .
5 Forsythe, G. E. and Wasow, W. R., Finite Difference Method for
Partial Differential Equations, Wiley (1960) .
6 Duderstadt, J. J. and Hamilton, L. O., Nuclear Reactor Analysis, 
John Wiley and Sons, Inc., New York (1976).
7 Levine, S. H., Workshop on Reactor Physics Calculations for
Applications in Nuclear Technology, ITCP Trieste, Italy (12
F e b . -16 Mar. 1990) .
8 Reactor Physics for Developing Countries and Nuclear Spectros-
9 7
University of Ghana          http://ugspace.ug.edu.gh
copy Research, World Scientific Publishing Co. Pte. Ltd.(1986)
9 Weinderg, A. M. and Wigner, E. P., The Physical Theory of
Neutron Chain Reactors, University of Chicago Press (1958) .
10 Bell, G. I. and Glasstone, Nuclear Reactor Theory, Van
Nostrand, Princeton, J. N. (1970).
11 Akaho, E. H. K. and Danso, K. A., Analysis of a Homogeneous 
Unreflected Subcritical Assembly Driven by Isotopic Sources, 
GAEC-NNRI (Nov. 1990) .
12 Meem, J. L., Two Group Reactor Theory, Gordon and Breach
Science Publishers, New York (1964).
13 Williams, M. M. R., Lecture Notes on Nuclear Reactor Theory 
Parti, Faculty of Engineering, Queen Mary College, University 
of London. (1970).
14 Levin, V. E., Nuclear Physics and Nuclear Reactors, MIR
Publications, Moscow (1981) .
15 Bennet, D. J. , The Elements of Nuclear Power, Longman Group 
Ltd. (1972).
16 Larmash. J. R., Introduction to Nuclear Engineering, Addision- 
Weslay Publishing Company Inc. (1975).
98
University of Ghana          http://ugspace.ug.edu.gh
17 Pollard, E. C., and Davidson, W. L., Applied Nuclear Physics,
John Wiley and Sons Inc., (1951).
18 Klimov, A, Nuclear Physics and Nuclear, Reactors, MIR
Publishers Moscow (1982).
19 Manfield, W. K., Lecture Notes on Elementary Nuclear Physics
and Reactor Physics, Faculty of Engineering, Queen Mary
College, University of London. (1970).
20 Reactor Physics Constants, US Atomic Energy Commission 
Report, ANL 5800 (1963).
21 Martelly, J., A Subcritical Experimental Setup, The Neutrostat, 
Proceedings of the 2nd. United Nations international Conference 
on the Peacefull Uses of Atomic Energy,
Vol.12, Geneva (1 13 Sept. 1958).
22 Compbell, C. G. and Grant, I. T., Critical and Subcritical
Experiments with Two-group Correlation of Result, Proceedings 
of the 2nd. United Nations International Conference on the 
Peacefull Uses of Atomic Energy, Vol.12, Geneva (1 - 13 Sept. 
1958) .
23 Technology and Uses of Low Power Research Reactors,
99
University of Ghana          http://ugspace.ug.edu.gh
IAEA, Vienna, IAEA-TECDOC 384. (1986)
24 Osae, E. K., Optimization Design and Construction of Polythene 
Moderated and Reflected Natural Uranium Subcritical Assembly, 
M .Sc., The Polytechnique of the South Bank. (1977).
25 Exponential and Critical Experiments, Proceedings of a Sympo­
sium, Amsterdam Vol.l, ( 2 - 6  Sept. 1963).
26 Isotopic Neutron Sources for Neutron Activation Analysis, 
IAEA, Vienna, IAEA - TECDOC 465.(1988).
27 Kock, R. C., Activation Analysis Handbook, Academic Press New 
York (1960) .
28 From Ideas to Applications (Proc. Meeting. Costa Rica Jose). 
IAEA, Vienna. (1978) .
29 Hunt, S. E., Fission, Fusion and Energy Crisis, Pergamon 
Press(1980) .
30 Osae, E. K. , The Principles of Neutron Activation Analysis, 
GAEC Technical Report. (1988).
31 Maakuu, B. T., Determination of Mineral Elements in Mutants 
of Mung Beans Using Neutron Activation Analysis, Project 
Report, Department of Physics, University of Ghana Legon (1988)
100
University of Ghana          http://ugspace.ug.edu.gh
32 Prince. W. J. , Nuclear Radiation Detection, McGraw Hill Inc. 
New York (1958)
33 Stephensen, G. , Mathematical Methods for Science Students, 
2nd. Edition, Longman Group Lt. (1973).
34 Watson, G. N., A Treatise on the Theory of Bessel Functions, 
Cambridge At The University Press (1944).
35 Maiorino, J. R., Workshop on Reactor Physics Calculations for 
Applications in Nuclear Technology, ITCP, Trieste, Italy (12 
Feb.- IS Mar. 1990).
36 Vorga, R. S., Matrix Iteration, Prentice Hall, London (1963)
37 Huynh, D. P., WIMSPC A Version of WIMS D/4 Adapted to AT
286/386 with Math Coprocessor 80287/80387 (1990).
38 Halsall, M. J., A Summary of WIMSD4 Input Option, AEEW -M1327, 
Atomic Energy Establishment, Winfrith, UK, (July, 1990) .
39 Askew, J. R. Fayers, F. J. and Kemshell, F. B., A General 
Description of the Lattice Code WIMS, Journal of the British 
Nuclear Energy Society 5, 4, 564 (1966).
40 Fayers, F. J. , Davison, W., George, C. H. and Halsall, M. J. ,
101
University of Ghana          http://ugspace.ug.edu.gh
LWR WIMS, A moderator Computer Code for the Evaluation of 
Light Water Reactor Lattices, Part 1, Description of Methods, 
UKAEA Report AEEW R 785, Winfrith, England (1972).
41 Akaho, E. H. K., Anim-Sampong, S. and Maakuu, B. T., Calcula­
tions for the Core Configurations of the Miniature Neutron 
Source Reactor, GAEC NNRI -RT - 13 (Aug.1992).
42 Mahama, M. S., and Akaho, E. H. K. , Mechanical Design 
Features of A Neutron Multiplier. GAEC -NNRI-RT-15 (Dec.1992).
102
University of Ghana          http://ugspace.ug.edu.gh
APPENDIX A: PROGRAM LISTING FOR SUNDES.FOR
c --------------------------------------------------------------
c The code SUNDES solves one-dimensional two-group neutron
c diffusion equation for a sub-critical assembly driven by
c isotopic neutron sources,Multiplying medium is U-235.
c ---------------------------------------------------------------
c For enquiries; Reactor Simulation Group
c National Nuclear Research Institute
c Ghana Atomic Energy■Commissi on
c P.O.Box 80
c Legon,Accra
c Ghana,Vest Africa
C — ---------------------------------------------------
c Options:
c ----------------
c <a) Choice of three geometries;Slab,Cylinder and Sphere
c
c <b> Determination of keff«
u
c <c> Variation of neutron flux with radius of assembly.
c ----- ----------------------------------------------------------
c List of parameters used in the program,
c '
c i radius along assembly
c 1— homogeneuos zone for material
c nl-total number of zones
c fn-fraction of fission spectrum in neutron energy groups
c fq-fraction of source spectrum for for neutron energy groups
c b— transverse buckling
c kg-energy group
c ng-total number of energy groups
c h— interval
c !k— position on radius of assembly
c kl-number of intervals
c ia-number of intervals for multiplying region
c rtmax-maxi mum number of allowed iterations
c ig-type of geometry<® for slab, 1 for cylinder and 2 for sphere,'
c sigf-macroscopic fission cross-section
c siga-macroscopic absorption cross-section
c sigs-macroscopic scattering cross-section
c scf— ratio of netrons released per fission to the thermal
c  energy released per fission(n/kWs)
c sigt-total neutron cross-section (cm-1)
c d--diffusion coefficient(cm-1>
c tau-parameter reqd.for calculating coefficient for matrix
c rf— radius of assembly(cm)
c p— power within the assembly(kV)
c erl-minimum allowed value for convergence criterion
103..
University of Ghana          http://ugspace.ug.edu.gh
c __________________________ _ __ _ _______________________
c Main Program
c -- ---------------------------------------------------------
dimension phi<2 ,100),siga(2 ,10),d<2,10),b <2,10) 
dimension scf <2,10),sigt(2 , 10),sigs(2,2, 10) 
dimens ion ml(10),dv(100),f (200),sigma <2 , 100) 
dimension scm<50),dd(50), gam(50),fi(50) 
dimension c (100,100),phi i(100),sx(100),f1(2 , 100) 
dimension g (100,100),xn(100),gx<100),ax(100) , efc(50) 
dimension a (2 ,100,100),bb(2 ,100),xx(2,100) 
dimension delr(10),noptsr(100),sp(2,10)
dimension sff <2,30),sa(2,30),sgs(2,2,30),of <2,30).buk<2,30; 
dimension fc(2,30),st(2,30)
c
common /mk/ hr(10,60)
common /ff/ fn(2,10),fq(2 ,10)
common /fq/ sigf(2,10),q(2 ,10)
common /xy/ x(100);y <100),w(100)
common /sq/ s(100),se(100),qq(100),qj (100
common /aa/ a (2,100,100)
common / TV / r (100)
common /tt/ tan(100) , beta(100),rho(100)
c onunon / u e / -u (100)
common /s t / eta (100)
common /gt/ sr(2,100),sm<2 ,100),sd(2 ,100)
common /dw/ ak(50 >
common /fr/ sdx <21 100),sdm<2,100)
common /bb/ bb(2,100),x x (2,100)
open <unit=l,file=’sundes,inp',status= ’ old
open <unit=6 , f i le= * sundes. out' , stat-us= ’ new’ )
c
c -------------------------------------------------------------------------------
c input for geometrical data for subassembly
C      ■ ■
read <1,19) j n
read (1 , 20) ig
read <1,21) rf
read (1 , 22) P
read <1,23) ng
read (1 , 24) nl
read <1,25) jl
read <1,25) kl
read (1 , 26) ka
read (1,27) h
read (1 , 28) nmax
read Cl,29) er 1
read < 1, 34) (delr (1)
104.
University of Ghana          http://ugspace.ug.edu.gh
input for nuclear data for the assembly
read <1,30) ((sigf(kg,1),1=1,nl),kg=l,ng)
read(1,31) (Csiga(kg,1),1=1,nl),kg=l,ng)
read(l,32) ((sigs(ng,kg,1),1 = 1,nl),kg=1, (ng-1))
read(l,33) ((d(kg,1),1=1,nl),kg=l,ng)
read(1,34) ((b(kg,l),l=l,nl),kg=l,ng)
read(1,34) ((sp(kg,1),1=1,nl),kg=l,ng)
read(1,35) ((scf(kg,1),1=1,nl),kg=l,ng)
read(1,36) ((q(kg,l),l=l,nl),kg=l,ng)
read(l,37) ((fn(kg,l),l=l,nl),kg=l,ng)
read(1,38) ((fq(kg,1),1=1,nl),kg=l,ng)
determination of number of mesh points along the radial dire 
sl=0 .
do 44 kg=l,ng 
do 45 1=1,nl
if(kg.gt.l) then 
sl=sl+sigs(ng,kg,1)
sigt(kg, 1)=siga(kg,l)+sl+(d(kg,1)*(b(kg,1)*#2)) 
else
sigt(kg,1)=siga(kg,1) +Cd(kg,l)*(b(kg,1)**2))
end if
cont i nue 
continue
do 2 1=1,nl 
es=0 .0
do 1 kg=l,ng
es=es+siga(kg,l)+sp(kg,1) 
noptsr(1)=delr(1)*es+l.0
number of points to be less than 20 per region 
if(noptsr(1).It.20)noptsr(1)=20 
continue
sumr=0. 0 
do 3 1=1,nl 
sumr=sumr+noptsr(1)
if(sumr.It.50)go to 5 
flt2r=sumr 
do 4 1=1,nl 
fltlr=noptsr(1) 
flt3r=fltlr/flt2r*50 
noptsr Cl)=flt3r
continue
University of Ghana          http://ugspace.ug.edu.gh
computation of mesh spacings for radial direction
c
c
j=®
do 6 i=l,nl
j =j +noptsr(i)
if (i,eq.1)jfl=j
if <i.eq. (nl-2))jf2=j
continue
nr=j +1
j sr=ng*nr
sumr=0 ,0
do 7 i=l,nl
noptr=noptsr < i)
if <1.gt.1) sum=sum+noptsr<i-1)
do 7 j=l,noptr
do 7 kg=l,ng
der=delr <i)
l=sumr+j
sff <kg,l)=sigf(kg,i) 
sa(kg,l>=siga(kg,i) 
sgsCng,kg,1)=sigs(ng,kg,i) 
df(kg, 1)=d(kg,i) 
buk(kg,1)=b(kg,i > 
fc(kg,l)=fn(kg,i) 
st(kg,l)=sigt(kg,i)
hr(i,1)=der/noptsr(i >
write(*,*)kg,1,noptsr(i)
7 continue
j r=l
hr(nl,nr)=hr(nl,nr-1)
sumd=0 . 
k=l
; do 227 k=l,nr
do 337 1=1,nl 
do 447 m=l, noptr
; wr i te (>K,-K) m, 1, k
surod=sumd+hr(1, m) 
r (k) =suxnd
: write(*,*>k,r(k>
k=k+l 
447 continue
337 continue
kl=k-l
v/rite (#, ^ ) k 1
ka=kl-2
106.
University of Ghana          http://ugspace.ug.edu.gh
ore la,x-1, 00 
nc=20 
c kl=nl*j1
c -----------------------------------------------------------
c subroutine PARA to calculate radii at positions k ,tau.rho,
c beta
c
c
c
c
write(%,*} kl
call para(noptr, nl, kl, ig)
pi=3.142 
term=l.
do 51 1=1j nl 
ml (i > =1
51 continue
; summation of isotopic sources in the system
smq=0 ,
do 404 kg’-1, ng 
do 405 1=1,nl 
smq=smq+q (kg, ,1)
405 continue
404 continue
c guess of initial fluxes for g=l ng
do 52 kg=1,ng 
do 53 k=l,kl 
phi(kg,k)=term 
53 continue
52 continue
c
c coefficients for zone 1
c
do 333 kg=1,ng 
c kg=l
c
1=1
k=l
c r(k)=hr(1,1)
r (k)=.070
r32=0. 5*(r(k+1)+r(k> )
dl= (( (r32X* (ig+1) >~ <r (k> ** (ig+1 > > > / <ig+15 > 
bx= (<r(k+l)-r(k)>/<(r32**ig> *d(kg, 1>)> 
tx=r (k) *bx* ( (sigt (kg, 1 > *d (kg, 1) ) **0 • 5) 
t, tx=0 .
a(kg, k, 1)-(1,+(sigt(kg,1)*dl*bx)+tx>
10(7.
University of Ghana          http://ugspace.ug.edu.gh
a (kg,k,2)=~1. 
c write (*, -*) j 1, kl
do 334 m=3,kl 
a (kg, k, m) =0 ,
334 continue
c coefficients for points k=2, , . . , kl
c
do 336 1=1,nl
if (1,eq.1)then 
J=2
kz=noptr
c
else 
j =kz+l 
kz=kz+noptr 
c if <1. gt. 1) write <*, *) 1, j , kz
end if
cal1 coef f(j,kz,kg,1,kl,d,sigt,a)
c
336 continue
c
l=nl
kll=kl-l
kl2=kl-2
kl3=kl-3
uk= < (r(kl)** (ig+1) ) - (rke** (ig+1) ) )/(ig+l) 
tk=(rke*>Kig)/(r(kl)-r(kll)) 
vk= ( Csiga (kg, l)Xd(kg, 1>)**0.5) 
a (kg, kl, kll)==- (tk#d (kg, 1))
c
do 338 m=l,kl2 
a (kg,kl, m)=0 .
338 continue
c
a (kg, kl, kl) = <d (kg, 1) Xtk) + (siga (kg, 1) *uk> + (vk* <r (kl) **ig; )
c
333 continue
c set iteration i = l
c
i = l
c
c summation of fission sources and isotopic sources
778 do 88 1-1,nl
i f (1.e q . 1) t h e n
j  = l
kz=noptr
108.
University of Ghana          http://ugspace.ug.edu.gh
cc
c
c
c
c
c
305
306
307
c
88
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
8 7 6
else 
j =kz+l 
kz=kz+noptr 
end If
ak <1) = 1. 
ak(2)=1.
caln. of values using subroutine SUMP,DUMP and TUMP
test for the kind of source used
if(jn.eq.l) go to 305 
if<jn.eq.2> go to 306 
if(jn.eq.3) go to 307
cal1 sump <i,j,kz,1,ng,kl, ig,pi,scf,phi,dv) 
go to 88
call dump <j,kz,1,ng,ig,pi,scf,phi,dv)
go to 88
call tump(j,kz,l,ng,kl,ig,pi,scf,phi, dv)
cont inue
m=kl-l
for fission sources in the entire volume 
of the assembly
call simpr(x, kl, m, h, smf> 
fi < i)=smf
for isotopic sources in the assembly
call simpr (w, kl, m, h, smw) ^
calculation of subcritical multiplication, m,i“ ! 
scm < i ) = (smq+smf ) /smq " .
calculation of factors for estimation ofkeff 
tt=smf/smq 
gam(i)=1. /scm(i) 
dd <i)=smf/smq
spq=®,0 
srp=0 .0
do 876 k=1, <kl) 
spq=spq+x(k) 
srp=sr p+se(k)
cont inue
108.
University of Ghana          http://ugspace.ug.edu.gh
px=3. 2e-ll*smf/2, 5
ax(i)=p/px
if(i.eq.1) go to 779
3 check for convergence
gt=fi Ci-l)/fi (i>
if(i.gt. 2) ak (i)=ak(i-1>*<fi<i-l)/fi<i))
: calculation of matrix B for g=l and C for g=2,
779 kg=l
c
do 102 k=l,kl 
sigma <kg,k)=sr <kg,k) 
102 continue
c
c
if(kg.eq.l) go to 500
777 kkg=kg-l
do 55 kg=l,kkg
c
do 888 1=1,nl 
if (1.eq. 1) then 
j = l
kz=noptr
c
else 
j =kz + l 
kz=kz+noptr
c
end if
c computation of coefficients for matrix C
do 222 k=j,kz 
pliii (k)=phi (kg, k) 
do 223 m=l, kl 
mk=m-k
c
if(mk.eq.0) then 
c (k, m) =sigs (ng, kg, 1) 
e Ise
c <k, m)=0•
c
end if
c
223 continue
222 continue
c
888 continue
110':
University of Ghana          http://ugspace.ug.edu.gh
product of matrix C and phii
55
809
500
7009
718
667
244
234
236
235
407
19
20 
21 
22
23
24
call pmatrx(c,phi 1,k l ,sx)
cont inue
do 809 k= l ,kl
sigma(kg, k)=eta(k)*sx(k)
continue
neqn=kl
do 7009 k=l,kl 
phi(kg,k)=0. 
cont inue 
do 718 k= l ,kl
sigma(kg,k)=sigma(kg,k)/eta(k) 
a (kg, k, m) =a (kg, k, ra)/eta (k) 
continue
cal1 gsiter(a,sigma,phi,kg,kl,orelax, nc) 
kg=kg+l
if(kg.le.ng) go to 777
if(kg.gt,ng) i=i+l
if(i,eq.nmax) go to 667
go to 778
do 234 kg=l,ng
do 244 k = l ,kl
f 1 (kg,k)=ax(i-l)*phi(kg, k)
continue
cont inue
str=0.
do 235 kg=1,ng 
do 236 k='l,kl
str=str+sigf(kg,1)*f1 (kg,k)
cont inue
continue
srq=0.
do 407 k=l,kl 
srq=srq+sd(l,k)/eta(k) 
continue
sdq=(srq+str)/srq 
f ormat(i4) 
format(14) 
format(f 8.3 > 
format(f8.3) 
format(14) 
format(i4)
University of Ghana          http://ugspace.ug.edu.gh
: 25 format(i4)
26 format(i4)
27 formatCfS.3)
28 format <i4)
29 format(f8,3)
30 format(9f8,5)
31 format(9f8,5)
32 format(9f8 .5)
33 format(9f8 .5)
34 format(9f8 .5)
35 format(9e8.2)
36 format(9e8.2)
37 format(9f8.5)
38 format<9f8.5)
output of nuclear and geometrical data of the assembly
write(6,551) 
write <6,552) 
write(6,553) 
write(6,554) 
write(6,555) 
write(6,556) 
write <6,557) 
write(6,558) 
write(6,559) 
write(6,660) 
write <6,661) 
write(6,662) 
write <6,663) 
wr ite(6,664) 
write(6,665)
if(ig.eq.0) wrlte(6,771) 
if(ig.eq.l) writ.e(6,772) 
if(ig.eq.2) write(6,773)
write(6,774)
write(6,775) rf 
write(6,776) ng 
write(6,7771) nl 
wr ite(6,7781) kl 
write(6,7791) erl 
write <6,800)
write (6,801) 
write(6,802) 
write(6,803) 
write(6,804) 
wr ite(6,805)
(ml(i),1=1,nl)
( (sigf (kg, 1) , 1 = 1, nl) , kg=l, ng) 
((siga(kg,1),1=1,nl),kg=l,ng)
((sigs(ng,kg,1),1=1,n l ),kg=l, (ng-1))
112.
University of Ghana          http://ugspace.ug.edu.gh
cc
c
c
c
c
551
552
553
554
555
556
557
558
559
660
661
662
663
664
665
771
772
773
c
774
775
776
777
778:
c779
c
800
Q
801
802
803
804
805
806
807
808
80S
810
write(6,806) ((d(kg, 1> , 1 = 1,nl>,kg=l,ng; 
write (6,807) (<b(kg, 1), 1 = 1, nl),kg=l,ng) 
write(6,808) ((scf(kg, 1) , 1=1, nl) , kg=l,ng)
write(6,8099) ((q(kg, 1),1 = 1,nl),kg=l,ng)
write(6,810) ((f n(kg,1),1 = 1,nl),kg=l,ng) 
write (6,811) ((fq(kg,1), 1=1,nl),kg=l,ng) 
write(6,812)
results of calulation of keff and neutron fluxes
write(6,813) 
write(6,814)
format(’ 
format(’ 
format(’ 
format(’ 
format(’ 
format(’ 
format(’ 
format(’ 
format(’ 
f ormat(' 
f ormat(’ 
format(’ 
format(’ 
format(1 
format(’
1 1 i i i ■ i i ! — \ j ----- ------  ;
■ ■ i I ; \  : : i i i
■ i i 1 ! \ ! ! / ; ;
**********************************************'
*   ,
*
%
*
si/ Nl' sL/ si/ -V  'J' "4/  -Ji- '&■ si/ \1/  -J/ sj/ sV si/ -4 /  sV sV sL* sL- .V SL> .1/  si/ sL- -J/ Nl/ si/ A' -U/ SL- ~i/ si/ ^  sL- Jy si/ -Js si/ -U "J/ 1// f \  /ys /js / | \  /y\ /[s /p. ,/f\ / | \  /js .'*y\ .'js /|S ./JS /|s .'JN -'js /|s 'J\ /] \ /f \ /^\ ,*j\ /f \  "[N ./ys .''fs /|s ,/js /Js yjs ,-jS /Js /f% /|*. /
E.H.K. Akaho & B.T. Maakuu 
Reactor Simulation Group 
Dept, of Reactor Technology 
National Nuclear Research Institute 
Ghana Atomic Energy Commission 
Kwabenya,Accra.
format (//’ 
format(//’ 
format(//’
format(’
Slab
Cylindrical 
Spherical
Geometry
Geometry
Geometry
//.> 
//> 
1 //)
format(5x, ’Radius of assembly(cm) ’ ,f8.3)
format(5x, ’Total number of energy groups ’ ,i4 >
format(5x, ’Total number of material zones ’ , i4 )
format (5x, ’ Number of finite difference points ' , i4 •'
format(5x,’Convergence criterion allowed for keff. ’,fl0.5)
format(1
format(//,5x, Nuclear Data //)
format(5x, ’1:‘,5il0)
format (5x, ’sigf:’,9f 8.5)
format(5x, ’siga;’,9f8.5)
format(5x,'sigs:*,9f8.5)
format(5x,’ d: ,9f8 .5)
format(5x,’ b: ,9f 8.5) ^  ’
format(5k ,'scf: ’,9e8,2)
m f?S iflf 3 3format(5x, ’ q; ,9e8.2)
f ormat(5x,’ f n: , 9f8 .5)
113..
University of Ghana          http://ugspace.ug.edu.gh
c 811 
812
c
813
814
818
817
820
821
822
823
825
824 
819
format(5x, ’ fq: ’ ,9f8,5)
format (5x, ' --------------------------------------------- ’ )
format(/,5x,’i d x m keff',
format (5x, ’ —  —  —  —  -----’
calculated values of subcrit, mult.,m ,keff & neutron fluxes
do 817 i=2, (ninax-1)
write (6,818) i,dd <i),gam<i),scm(i),ak (i ) 
format(6x, 13, 8x, 5el0. 3) 
cont inue
do 819 kg=1,ng 
if(kg.eq.l) write <6,820)
format</,15x,’Mat. Zone Radius Fast Neutron Fluxes’,/)
write <6,821)
format (5x, ’      ’ )
if <kg.eq,2) write(6, 822)
format(/,5x,’Mat. Zone Radius Thermal Neutron Fluxes
write <6,823)
f ormat (6x, ’      *
do 824 k=l,kl 
1= <k+noptr-1)/noptr
write(6,825) 1,r(k),phi <kg,k),f1(kg,k) 
format(10x,i3, 10x,f7.3,15x,2el5.5) 
continue 
continue
close (1) 
close(6)
end
subroutine Para determines parameters for the 
evaluation of coefficients for the matrix A
subroutine para(noptr,nl,kl,ig)
common /mk/ hr(10,60) 
common /rr/ r(100)
common /tt/ tau(100),beta(100),rho(100) 
common /ue/ u(100) 
common /et/ eta(100)
c
kll=kl+l
ii:4.
University of Ghana          http://ugspace.ug.edu.gh
k— 1
r32=0 . 5* (r (k+1) +r (k) )
eta (k)= (r32*X (ig+1))-(r(k)XX (ig+1) }/(ig+1) 
k=2
do 14 1=1,nl 
do 15 m=l,noptr
tau(k) - ( (0 . 5* (r (k) +r (k-1) ) ) **ig)/r;r Cl, m) 
rho(k)=((0.5X(r(k+l)+r(k)))XXig)/hr(1,m)
u(k) = ( (r <k) XX(ig+1) )- (0. 5* (r (k)+r (k-1) (ig+1) ) )/(ig+l)
beta(k)= ( (0 .5X(r(k+1)+r(k))X*(ig+1))-((r(k))X*(ig+1)))
X /(ig+1)
eta(k)=u(k)+beta(k)
k=k+l
m=m+l
continue
continue
rkl=0.5*(r(kl)+r<kl-l>)
eta(kl)=((r(kl)XX(ig+1))-(rklXX(ig+1)))/(ig+1) 
eta(kl)=eta(kl-l)
return
end
subroutine Coeff calculates the coefficients for tridiagonal 
matrix A for point k=2....  (kl-1)
subroutine coeff(j,kz,kg, l,kl,d,sigt,a)
dimension a(2 ,100,100) 
dimension sigt(2, 10),d (2 , 10)
common /aa/ a(2,100,100)
common /tt/ tau(100),beta(100),rho(100) 
common /ue/ u(l®0) 
common /et/ eta(100)
do 1 k=j,kz 
do 2 m=1, k1
rak=m-k
lower coefficients,m=k-l
if(ink,eq.-1) then 
a(kg, k, m)=-d(kg, 1)Xtau(k)
115.
University of Ghana          http://ugspace.ug.edu.gh
upper coefficients m=k+l
else if (mk.eq.l) then 
a(kg,k,m)=-d(kg,1)*rho(k)
diagonal coefficients
else if(mk,eq,0)then
a (kg, k, m) =+d (kg, 1) *tau (k) + (d (kg, l)*rho Ck)) 
+ (eta(k)*sigt(kg,1))
for any other coefficients
else
a(kg, k, m)=0.
end if
continue
noptr=7
if(1.eq,1) a(kg, k, m)=a(kg, noptr, m) 
continue
return
end
subroutine Sump determines the summation of fission sources 
and isotopes.It also calculates points for integration by 
Simpson,s Rule
subroutine sump(i,j,k z ,1,ng,kl,ig,pi,scf,phi,dv)
dimension d v (100),p h i (2,100),s c f (2,10)
common / t t / r(100)
common /fq/ s i g f (2,10),q (2,10)
common /if/ f n (2,10),f q (2,10)
common /xy/ x (100),y (100),w(100)
common /sq/ s (100),s e (100),q q (100),q j (100)
common /et/ e t a (100)
common /gt/ s r (2, 100),sm(2,100),s d (2 , 100) 
common /fr/ sdx<2,100),sdm(2,100) 
common /dw/ ak(50)
do 22 k=j,kz 
if(lg,eq.0) dv(k)=l. 
if(ig.eq.l) dv(k)=2.Xpi*r(k) 
if (ig. eq. 2) dv(k)=4. *pi* (r (k) **2) 
cont inue
initialise values to zero 
do 222 k=j,kz
116.
University of Ghana          http://ugspace.ug.edu.gh
cc
c
222
c
sb=0 . 
ss=®. 
sq=0 . 
sw=0. 
sz=0 .
do 66 kg=l,ng 
if <i.gt. l)then
sb=sb+(f n(kg,1)Xsigf(kg,1> Xphi(kg,k)) 
e lse
sb=sb+(fn(kg,1)Xsigf(kg,1)Xphi(kg,k)) 
end if
ss=ss+(fq(kg,1)Xq(kg, 1)) 
write (6, X)kg, 1, k, phi (kg, k)
sq=sq+(((sigf(kg,l)Xphi(kg,k))/scf(kg,1) ) ) 
if(i . gt, 1)then
sw=sw+ ( (f n (kg, 1) Xsigf (kg, 1) Xphi (kg, k))) 
e lse
sw=sw+((fn(kg,1)Xsigf(kg,1)Xphi(kg, k))) 
end if
sz=ss+(fq(kg,1)*q(kg, 1>)
c
66 continue
x(k)=sbXdv(k) 
y(k> =ssXdv(k) 
w(k)=sqXdv(k) 
s(k)=swXeta(k) 
se(k)=sw 
qq(k)=sz
c
si(1,k)=s(k) 
c the case of infinite line source(ii>0)
yy=2,405
aj1=(yy/2.)-((yyXX3)/16.)+(yyXX5)/384. 
bb= (pi* (r (kl)XX2)X2. *r (kl)X (aj 1XX2))
xx=((2.405Xr(k))/r(kl))
aj 0=(1.-((xxXX2)/4. ) + ((xxXX4>/64.)-((xxXX6)/2304. ))
sd(1,k)=eta(k)X((qq(k)*aj0)/bb)
if(s d (1,k).It.0 .)sd(l,k)=0.0
sr(l,k)=sd(l,k)+sm(l, k)
continue
return
end
University of Ghana          http://ugspace.ug.edu.gh
c subroutine DUMP calculates only point source problem<jj>0>
c -----------------------------------------------------------
subroutine dumpCj,kz,1,ng,ig,pi,scf,phi,dv)
c
dimension dv<100),phi(2,100),scf(2,10)
common /rr/ r<100)
common /fq/ sigf<2 ,10),q (2,10)
common / i f /  fn(2,10),fq (2, 10)
common /xy/ x(100),y (100),w(100)
common /sq/ s (100 ) , se (.100) , qq (100) , qj (100)
common /et/ eta(100)
common /gt/ sr(2, 100),sm(2 , 100),sd(2, 100)
c
do 22 k=j,kz
c
if(ig.eq,0) dv(k)=l, 
if(ig.eq.l) dv(k)=2 .*pi*r(k) 
if(ig.eq. 2) dv(k)=4,*pi*(r(k)**2)
22 continue
do 222 k=j,kz
c
sb=0 . 
ss=0. 
sq=0 . 
sw=0. 
sz=0.
do 66 kg=l,ng
sb=sb+(fn(kg,1)*sigf(kg,1)t phi(kg,k)) 
ss=ss+(fq(kg,1)*q(kg,1))
sq=sq+(((sigf(kg,l)Xphi(kg,k))/scf(kg,1))) 
sw=sw+((fn(kg,1)Xsigf(kg,1)*phi(kg,k))) 
sz=sz+(fq(kg,1)*q(kg,1))
66 continue
x(k)=sb*dv(k) 
y(k)=ss*dv(k) 
w(k)=sq*dv(k) 
s (It) =swteta (k) 
se <k)=sw 
qq(k)=eta Ck)*sz 
sm(l, k)=s(k)
if Ck.eq.j) then 
sd(1,k)=qq(j) 
else
sd(1,k)=0 . 
end if
wr ite f 6 , ) 1, k , sd (. 1, k )
University of Ghana          http://ugspace.ug.edu.gh
s j i (1 ,  k) = s  (k)
222
o
c
c
2
c
c
c
c
c
sz=2,405
aj 2= <zz/2 . ) - ( <zz**3>/16 . ) + <zz**5)/384. 
cc= <pi* <r (kl) **2) * (a j 2 *.*2) ) 
vv=((2.405*r(k)>/r(kl))
a j  3 =  < 1 .  -  ( ( v v X X 2 ) / 4 . ) + < < v v * * 4 ) / 6 4 .  > - a v v * * 6 ) / 2 3 0 4 . > ;
if <r<k). It,r(kl)) then 
sdxCl,k) = <(qj (k)Xaj 3)/cc) 
else
sdxd, k) =0 , 
end if
if<Cl.gt,1),and, <k.eq.j)) then 
sdm <1,k)=qq <j) 
e lse
sdirKl, k)=0, 
end if
sr(1,k)=sdx(1,k)+ sdm(1, k)+sm(l, k) 
continue
return
end
subroutine Simpr subprogram for Simpson's Eule
subrout ine simpr <f,ndim,nstep, h, sint)
dimension fCndim+l)
apply Simpson's rule
sint=0 .
do 2 j = 2 , nstep,2
si nt=sint + f (j-l)+4. *f <j )+f ( j +1) 
sint=sint*h/3.
return
end
subroutine Pmatrx calculates the product of matrices
subroutine pmatrx (g, sen, ni , gx)
dimension g (100, 100),xn(100),gx<100) 
multiplication of matrices
do 2 i=l,ni
c 119.
University of Ghana          http://ugspace.ug.edu.gh
c do 1 j=l,ni
gg=gg+8<i> J >*xn(j>
1 continue 
gx ( i >=gg
2 continue
return
c
end
c ------------------------------------------------------------
c subroutine Gsiter for solving simultaneous linear equations
c using Gauss-Seidel Method
c
c
subroutine gsiter(a, bb, xx, kg,meqn,orelax, nc)
c
dimension bb(2,100),xx(2,100)
c
dimension a(2,100,100) 
neqn=meqn
c
c set up iteration loop
do 5 iter=l,nc
sumx=0. 
sumdx=0.
c. obtain new estimate for each unknown in turn
do 3 i=1,neqn 
deltax=bb(kg,i) 
do 2 j = 1,neqn
2 deltax=deltax-a(kg,i,j)Xxx(kg,j)
c wr ite (>', *) kg, i , a (kg, i, i >
deltax=deltax/a(kg, i, i)
sumdx=sumdx+abs(deltax)
xx(kg,i)=xx(kg,i)+deltax*orelax
3 suinx=sumx+abs (xx (kg, i ) ) 
return
5 continue
return
end
120.
University of Ghana          http://ugspace.ug.edu.gh