A Multigroup Eigenfunction Expansion Technique
For Selecting Poison Distributions
For Optimum Reactor Power Distribution
By
STEVEN MICHAEL LONG
A DISSERTATION PRESENTED TO THE GRADUATE
COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTIAL
FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1972
This dissertation is dedicated to my lovely wife,
Diane,
who ungrudgingly shared the effort of completing
the beast during our first months of marriage.
ACKNOWLEDGEMENTS
The author wishes to express his appreciation to the
chairman of his supervisory committee, Dr. G.R. Dalton,
who originally suggested the subject of this study and
whose support and criticisms greatly aided in the completion
of this dissertation. The author thanks Dr. C.D. Kylstra
for his help with the optimization problem and for the use
of his computer optimization code. The author thanks Dr.
N.J. Diaz for his aid and advise in handling many other
computer problems. The author also expresses his apprecia
tion to Dr. D.T. Williams and Dr. U.H. Kurzweg for serving
as members of his supervisory committee.
In addition, the author acknowledges the financial
support provided by the U.S. Atomic Energy Commission in
the form of a Special Fellowship, by the University of
Florida in the form of a College of Engineering Fellowship,
and by the Department of Nuclear Engineering Sciences in
the form of an assistantship.
iii
TABLE OF CONTENTS
Page
ACKNOWLEDGEMENTS. . .
LIST OF TABLES. . . .
LIST OF FIGURES . . .
NOMENCLATURE . . .
ABSTRACT. . . . .
CHAPTERS:
I. INTRODUCTION.
II. THE ONEGROUP
PROBLEM . .
EIGENFUNCTION EXPANSION
III. THE POISON DISTRIBUTION PROBLEM IN
MULTIGROUP FORMULISM. . . . . .
IV. SELECTION OF A MULTIGROUP EIGENFUNCTION
SET . . . . . . . . .
V. GENERATION OF THE POISON EIGENVECTORS
VI. THE MULTIGROUP OPTIMIZATION PROBLEM .
VII. NUMERICAL CALCULATIONS AND THE PROSPECT
FOR PARTIAL EXPANSIONS . . . .
VIII. CONCLUSIONS AND SUGGESTIONS FOR FURTHER
STUDY . . . . . . . . .
DICES:
I. A DESCRIPTION OF THE TECHNIQUES
EMPLOYED BY THE EIGENVECTOR GENERATING
PROGRAM "CODE2" . . . . . .
* . 6
* .12
* .18
* .35
* .47
. .65
. 104
. 113
. iii
. vii
. .viii
S. xxii
. . 1
APPEN
* .
: : :
. . . . . .
Page
II. LISTING OF THE DEMONSTRATION POISON
EIGENVECTOR GENERATION PROGRAM
"CODE2" . . . . . . .. . .119
III. A DESCRIPTION OF THE "GAUSS" OPTI
MIZATION CODE. . . . . . . .130
IV. LISTING OF THE USER SUPPLIED SUBROU
TINES FOR THE USE OF THE "GAUSS"
PROGRAM TO FIND THE OPTIMUM CRITICAL
POISON DISTRIBUTION . . . . .... .135
LIST OF REFERENCES . . . . . . . . .. .148
BIOGRAPHICAL SKETCH. . . . . . . . . .151
LIST OF TABLES
Table Page
1 Reactor Model for Example 1 68
2 The Eigenvalues Associated With the Poison
Eigenvectors for Example 1 71
3 Optimization Results for Example 1 73
4 Reactor Model for Example 2 78
5 ThreeRegion Poison Optimization Results for
Example 2 81
6 Reactor Model for Example 3 85
7 Comparison of Reactor Configurations Resulting
From the TwoRegion Poison Optimization of
Example 3 Using Full and Partial Expansions 89
8 Comparison of the Expansion Coefficients for the
Full and Partial Expansion Optimums for
Example 3 90
9 Results of Partial Expansions for a Completely
Specified System (Example 1) 95
10 Reactor Model for Example 4 99
11 The Poison Eigenvectors 1 and 05 for
Example 4 101
LIST OF FIGURES
Figure Page
1 The [B] Matrix and Flux Vector for the
"Fission Eigenfunctions" 30
+ +
2 Poison Eigenvectors + 1', 2 and <2
for Example 1 69
+ +
3 Adjoint Vectors i,' ~' 2 and ~2 for
Example 1 70
4 CCmparison of the Poison Eigenvector,
,for Example 2 As Calculated by
CODE2, by CORA and Analytically 79
5 Flux Distribution for the Optimum Three
Region Poison Distribution in Example 2 82
6 Poison Eigenvectors P1 and 1i for Example
3 86
7 Adjoint Vectors i1 and 1 for Example 3 87
8 Power Distribution for the Optimum Two
Region Poison Distribution In Example 3 92
9 Schematic Representation of Reactor Model
for Example 4 100
vii
NOMENCLATURE
a
[A]
B2
z
[B]
C(k,i,j)
k
D (r)
D(k,j)
E
E
F(m,i, [N])
[F]
G
H
I
i
energy per fission (in appropriate units)
expansion coefficient associated with function
i
vector representing a set of expansion
coefficients
matrix on lefthand side of various eigen
value formulations of the criticality equa
tion
axial geometric buckling
matrix on the righthand side of the general
eigenvalue formulation of the criticality
equation
coefficient of the flux in group k at point
j in the gradient term of the equation for
the flux in group k at point i
diffusion coefficient function in energy
group k
diffusion coefficient at point j in group k
weighted scalar squared error optimization
parameter
vector of error components
the scalar product of the poison eigenvector
i and the adjoint vector m, weighted by the
operator [N] [o ]
matrix with elements F(m,i,[N])
number of control points
number of fueled points
number of eigenvectors
preferred index for designating a particular
eigenvector
viii
[I] the identity matrix
J number of spatial points
J (r) zeroth order Bessel function of the first kind
[J] the Jacobian matrix
j preferred index for designating a spatial
point
K number of energy groups
Keff reactivity of a reactor
k preferred index for designating an energy group
[L] matrix containing the neutron loss terms of
the criticality equation
M number of vectors used in expansion
m preferred index for designating a particular
adjoint vector
N(r) poison distribution function
N(j) poison concentration at point j
N vector with elements N(j)
[N] diagonal matrix with elements N(j)
n auxiliary index for indicating energy group
[0] zero matrix
P(r) power distribution function
P*(r) ideal power function
P(j) power density at point j
Pi power vector associated with eigenvector 4i
P* ideal power vector
[Q] thermal poison matrix
r radius
r general position vector
sk(r) spectrum function relating the flux in group
k to the flux in group K
S. fast/thermal flux ratio for eigenfunction i
S (spatially constant)
[S] matrix containing the fission source terms
of the criticality equation
V arbitrary vector
V(j) fraction of reactor volume associated with
mesh point j
[V] diagonal matrix containing the inverse velo
city terms of the Laplace transformed kinetic
equation
[W] diagonal matrix with elements V(j)
X distance from centerplane in slab geometry
ai eigenvalue associated with the eigenfunction
i to a general eigenvalue problem
8i eigenvalue associated with poison eigenfunc
tion i
[8] diagonal matrix with elements 8
Yi eigenvalue associated with the "natural"
eigenfunction i
6. Wielandt transformed eigenvalue associated
1 with the mode i
6(i,j) Kr6necker delta
e scalar coefficient for the adjustment of the
magnitude of a poison distribution
e. poison distribution magnitude eigenvalue
1 associated with mode i
n. zero number i of the function Jo(r)
X. eigenvalue associated with the fission eigen
1 function i
Pi eigenvalue associated with the thermal poison
eigenfunction i
v number of neutrons per fission
k
Cy M
op
k(
o(kj)
Za(r)
a (k,j)
k
f (nk,j)
Ek (r)
ZR (k, j)
Xk or X(k)
(k,j)
poison cross section for a unit concentration
equals 0 if point j is a control point,
zero otherwise
diagonal matrix with elements oa(j)
absorption cross section function for group k
absorption cross section for group k at
point j
fission cross section function for group k
fission cross section for group k at point j
controllable poison cross section function
for group k
removal cross section function for group k
removal cross section for group k at point j
cross section at point j for scattering from
group n to group k
flux function for group k
ideal flux function
flux at point j in group k
flux vector with elements f(k,j)
flux eigenvector i
fission spectrum fraction in group k
adjoint flux at point j in group k
adjoint flux vector m
eigenvalue associated with kinetic mode i
Abstract of Dissertation Presented to the Graduate Council
of the University of Florida in Partial Fulfillment of
the Requirements for the Degree of Doctor of Philosophy
A MULTIGROUP EIGENFUNCTION EXPANSION TECHNIQUE
FOR SELECTING POISON DISTRIBUTIONS
FOR OPTIMUM REACTOR POWER DISTRIBUTION
By
Steven Michael Long
December, 1972
Chairman: Dr. G. Ronald Dalton
Major Department: Nuclear Engineering Sciences
The purpose of this study is to explore the capabilities
of an eigenfunction expansion technique for aiding in the
selection of a poison distribution that will produce the
optimum power distribution in a reactor. A method is
introduced which is useful in the repeated evaluation of
power distributions resulting from various poison distri
butions as they are selected for investigation by a numerical
optimization routine. This method is based in multigroup
diffusion theory, and employs the finite difference approxi
mations to the spatial dimensions.
First, the eigenfunction expansion problem is reviewed
in onegroup theory. Then, the multigroup problem is out
lined and the arbitrariness of the fluxpower relationship
in this formulism is discussed. The actual unknowns for the
optimization problem (using a preselected poison material)
are pointed out to be the poison concentrations and the flux
spectral relationships. The inherent nonlinearity of the
xii
problem is discussed for the case where these unknowns are
treated directly as variables.
The orthogonality and completeness properties of vari
ous sets of eigenvectors are discussed regarding the poison
optimization problem. These sets include the natural
modes, the fission modes, the timespace separable kinetic
modes, the thermal poison modes and the multigroup poison
modes.
The multigroup poison eigenvectors are chosen for use
in the expansion technique. A Wielandt transformation method
for generating this set of eigenvectors is explained, and
a computer program, demonstrating the technique, is listed
as an appendix. The poison eigenvectors for several two
group, multiregion reactors are presented.
The necessity for solution of the optimum multigroup
poison distribution problem by numerical search is discussed
in general. The poison eigenvalue formulation of the cri
ticality equation is introduced, and its advantages and
weaknesses, as compared to the usual reactivity eigenvalue
formulation, are discussed. Two methods are suggested for
overcoming the weaknesses. One, a Wielandt transformation
method, is conceptually straightforward, but computationally
clumsy. The second, a transformation of the set of eigen
values of the poison eigenvectors, is potentially a better
computational device, but is conceptually much more complex.
Only the first method is programmed and demonstrated by
numerical example.
xiii
Sample twogroup, onedimensional optimization searches,
using the Gaussian least squared error method, are presented
for several examples in order to demonstrate the techniques
and phenomena that have been discussed. The results are
checked for accuracy with a standard multigroup diffusion
theory code. These calculations are then extended to explore
the possibility of using incomplete sets of poison eigen
vectors for expansion. The results of partial expansions in
both underspecified and completely specified optimizations
are discussed. Considerable accuracy of partial expansions
is shown in the limited class of examples considered (ther
mal reactors).
Finally, the author presents his thoughts on the poten
tial utility of the method and makes suggestions for future
work directed at ultimately incorporating the techniques in
a standardized computer code for use as an engineering tool.
xiv
CHAPTER I
INTRODUCTION
The purpose of this study is to explore the capabilities
of an eigenfunction expansion technique in the problem of
reactor power shaping by neutron poison distribution. The
mathematical model used for the reactor core is that of multi
group diffusion theory. The spatial dimensions are repre
sented by the finite difference approximation, and the func
tions used for expansion are the eigenvectors to the problem
as expressed in this framework. The author believes that
this technique will provide a flexible and realistic repre
sentation of the physical problem and lend computational
power as well.
The problem of determining the control poison distribu
tion for the optimization of some reactor performance para
meter has been treated by many authors, using various methods
and mathematical formulisms. Harker (1) used a pseudo two
group diffusion theory that made use of a "fast fission fac
tor", which allowed him to determine a flux distribution
from a desired power distribution. He could then calculate
the distribution of a thermal poison necessary to achieve
criticality for this desired flux shape. Unfortunately, he
provided no a priori method of guessing a power distribution
which would result in a physically realizable flux shape.
(2) (3)
Haling(2) and Crowther(3) used a more sophisticated two
group theory and a "backwards iteration technique" to find
the flux shape, from which they calculated the necessary
distributions of a thermal poison. Crowther applied this
"Power Control Method" to the problem of minimizing the peak
(4)
to average power ratio over core life. Motoda and Kawai(
applied the Pontryagin Principle of optimal control to the
problem of maximizing core lifetime while limiting the power
peaking ratio. They used onegroup, twomode formulism to
represent a tworegion, bare core with control rods ganged
in each region. Terney and French(5) used dynamic program
ming to minimize power peaking over core life, based on a
direct flux synthesis method (using two fluxes), for ap
proximating a twogroup, twodimensional reactor with two
gangs of control rods. Wade and Terney(6) applied optimal
control theory to a generalized set of design objectives,
using a onegroup spatially nodalized reactor model. Jack
son(7) used an eigenfunction expansion technique, based on
onegroup diffusion theory, to determine the optimum physically
realizable power and flux distributions and the poison dis
tribution necessary to maintain these shapes in the critical
state. Because of the direct relationship to the subject
of this paper, Jackson's work will be described in some
detail in Chapter II, which reviews the onegroup eigenfunc
tion expansion problem.
Eigenvalue problems are well known to the nuclear
engineer in such forms as the reactor criticality equations
and the kinetic equations. However, except in Jackson's
work and in some kinetics work, little use has been made of
the higher modes to these problems. Even for perturbation
calculations, only the fundamental eigenfunction (and its
adjoint) are in general use for the solution of static
problems. Such first moment calculations are only the barest
form of eigenfunction expansion.
Investigative interest in eigenfunction expansion
techniques is not new, even if not widely spread. In 1960,
Foderaro and Garabedian introduced a method for the solu
tion of group diffusion equations which involved expansion
of the problem in the eigenfunctions of the Helmholtz equa
tion for the appropriate geometry. An article published in
1962 by Garabedian and Thomas(9) extended the approach to
two dimensions. While these methods have the advantages of
using wellknown sets of tabulated functions and analyti
cally handling the leakage terms within homogeneous regions,
they possess a distinct disadvantage. Because these func
tions are not actually eigenfunctions of the multiregion
problem being considered, they require individual (but
coupled) expansions for each materials region and require a
large number of terms to adequately represent the solution.
In 1968, Hara and Shibata0) combined this approach with
variational calculus in a study to optimize total power and
conversion ratio. Their independent variables were region
wise fuel concentrations and region boundary locations.
Still these works are best described as expansions in sets
of tabulated, orthonormal functions, rather than as modal
expansions.
True modal expansion techniques for kinetic problems
were explored by Kaplan(11) in 1961. He illustrated tech
niques using three different sets of eigenfunctions, each
generated from a different eigenvalue equation associated
with his problem. He introduced the property of "finality",
that the value of the expansion coefficients are independent
of the order of the expansion, and showed how this property
could be achieved for certain kinetics problems. He termed
the set of eigenfunctions having this property for a parti
cular problem the "natural modes" for that problem. In
1966, Kaplan and Yasinsky(12) used an analysis of the eigen
values of the natural modes to predict the behavior of xenon
oscillations. Henry,13) in 1963, extended modal approxi
mations to the nonseparable spacetime kinetics problems,
defining what he termed the "inhour modes". These occur
in clusters of seven, with similar spatial shapes within
clusters.
(14)
In 1963, the MULE computer code(14) became available
for computing three particular types of eigenvector sets to
twogroup diffusion theory reactor models. These eigenfunc
tions sets included: 1) the "X modes" or"fission eigenfunctions"
where X is the inverse reactivity for the mode, 2) the
"p modes", where p is the concentration multiplier for a
given normalized distribution of thermal poison, and 3)
the "w modes" or eigenfunctions to the Laplace transformed,
timespace separable kinetics problem. As will be discussed
in Chapter IV, with the exception of the w modes, these modes
are degenerate in the energy dimension, and are limited,
for expansion purposes, to expansions of the spatial dimen
sions, alone.
In this study, an expansion technique using a set of
eigenvectors capable of representing both the energy and
space dimensions is developed. A method for generating such
eigenvectors will be presented, and its use demonstrated.
These vectors and the expansion technique will then be dis
played in the determination of the poison distributions and
concentrations necessary for criticality in the optimum
realizable fluxpower shapes for several example reactors.
CHAPTER II
THE ONEGROUP EIGENFUNCTION EXPANSION PROBLEM
Although a treatment of the onegroup problem would
seem to be the simplest case of a general multigroup treat
ment, such generalization is not necessarily possible.
Because the onegroup formulism does not contain the energy
variable, it is amenable to a quite different treatment
than is a multigroup problem. As will be discussed in
Chapter III, this is due to the nonlinear nature in which
the energy variable enters the poison distributionpower dis
tribution problem. Still, it is instructive to review the
onegroup treatment for purposes of later comparison and
for the purpose of introducing some concepts in a less com
plicated context. The onegroup problem and solution was
described by Jackson in reference 7, and his approach will
generally be followed in this chapter.
In the onegroup formulism, the criticality equation
is
V D(r)V (r)+v f(r)((r)ZE (r) (r)z (r)4(r) = 0 (1)
plus the boundary conditions on (r) and/or its first
derivative. Standard notation is used except that E (r) is
the macroscopic cross section of the controllable poison.
The powerflux relationship is
P(r) = AEf(r) (r), (2)
where A is the energy per fission in appropriate units. At
first glance, it may appear that, once a desired power
shape, P*(r), is chosen, the resultant flux could be found
directly from equation (2) as
P*(r)
*(r) (3)
AZf (r)
and that the desired distribution of control poison follows
directly from equation (1) as
V.D(r)V**(r)
Ep(r) =f I(r)Ea () (4)
p f *(a) f a
This is the approach taken by Harken,() Haling(2) and
Crowther.(3) As Jackson(7) pointed out, however, most
reactors of interest will have step variations in Ef(r),
due to the practice of regionwise fuel loading. Thus, the
ideal flux as defined by equation (3) must also have step
variations. Physically, the flux cannot behave in this man
ner. As represented by diffusion theory, the neutron flux,
O(r), and the neutron current, D(r)V7(r), are constrained
to be continuous.
Jackson proposed an eigenfunction expansion technique
to get around this problem. Following Kaplan,() he chose
the set of "natural modes" to his problem as defined by the
equation,
VD(r)V (r)+vEf (r) i(r) ZEa(r)(i(r) = yi i(r). (5)
These eigenfunctions, the iq, form a complete, selfadjoint
set(7) over the entire reactor. They satisfy all of the
boundary and interface conditions required of the real flux,
and, thus, any linear sum of these eigenfunctions will also
be a solution to equation (1).
A physically realizable approximation to the "ideal
flux shape" given by equation (3) can then be constructed
from these eigenfunctions. Using the orthogonality relation
J jdr = 6(i,j), (6)
the expansion coefficients are found by
vR
a .= R i*dr. (7)
These coefficients are found independently of one another,
or the order of expansion, I, and thus the expansion has the
property of finality. The constructed optimum real flux,
I
*(r) = ai i(r), (8)
i=l
approximates the "ideal flux", 4*, in the least squared
errors sense(15) to the degree of expansion, I. Thus, the
optimum poison concentrations may be determined from an
equation analogous to equation (4), but written for p*R
R*
After invoking the eigenvalue relationships of equation (5),
this becomes
I
aiyi i (r)
p (r) i=(9)
i=l
As pointed out by Jackson, this approach still has a
failing. Because the optimum power shape, P*(r), is defined
only for the reactor core, where Ef(r)>0, but these eigen
functions are orthogonal only over the complete reactor,
the expansion coefficients must contain an ambiguity if a
reflector is present. This is due to the arbitrariness of
(*(r) in regions where Zf(r) = 0 and the necessity of inte
gration over those regions. Jackson suggested that the
flux in the reflector might either be assumed unchanged from
the unpoisoned case, or that it be given a constant value,
"treating its magnitude as an adjustable parameter, the
adjustment of which might allow the eigenfunction approxima
tion to the ideal incore flux to be close in that core re
gion in which the nearness of the approximation is considered
to be most important." (724)
Neither suggestion seems particular realistic. The
second has the added disadvantage of somewhat linking the
solutions of the expansion coefficients by evaluating them
together during the optimization search for the reflector
flux magnitude. This not only violates the conditions nec
essary for the property of finality, it also introduces an
added complexity which seems to destroy much of the power of
the method.
Recognizing a better alternative himself, Jackson
recommended the exploration of other sets of eigenfunctions
which are complete and orthogonal over just the fueled por
tions of the reactor. He suggested a set he called the
"fission eigenfunctions", given by the equation
VD(r)V. i(r)E (r)~.(r) = X.vEf(r) .i(r) (10)
1 a 1
and showed that they satisfy the orthogonality relation
I i(r)v f(r)4 (r)dr = 6 (i,j). (11)
Because of the weighting factor, f (r), fluxes in unfueled
regions have no effect on the computation of the expansion
coefficients for this vector set, as they carry zero weight.
The expansion coefficients are thus found from
ai = 1 (r)V f(r)4*(r)dr. (12)
O
The optimum real flux is constructed as in equation (8),
I
R(r) = aii (),
i=l
and the optimum poison concentrations are again determined
by an equation analagous to equation (4). After invoking
the eigenvalue relationships of equation (10), the poison
concentrations are given by
I
E (r) a (1()
Ei (F) = (13)
aili(r)
i=l
This method at last seems to handle the onegroup prob
lem rather well. Although it does not allow for control
poisons in the reflector, that location is not popular for
poisons (other than the soluable type) in the kind of
reactors at which this analysis was aimed. Other sets of
eigenfunctions are possible which allow for reflector con
trol rods, but these bring with them the requirement for
optimization searches which are not otherwise required in
the onegroup formulism. Such a set will be introduced
and used by this author in the following chapters, as he
extends the approach to multigroup formulism, where such
searches apparently cannot be avoided.
CHAPTER III
THE POISON DISTRIBUTION PROBLEM IN MULTIGROUP FORMULISM
In contrast to some of the previous authors,3) who
used modified forms of twogroup formulism in the hope of
thereby obtaining direct solutions, this author has chosen
to use the complete multigroup formulism. In principle, at
least, fissions and births will be considered in all groups.
The control poison will be generally considered to have a
cross section in each group. Because of this general treat
ment, the resulting method should be applicable not only to
the analysis of the current water moderated reactors, but
also tofast reactors.
In the multigroup diffusion theory formulism, for K
groups, the criticality equation becomes a set of K coupled
equations, one for each energy group. Each of these equa
tions has the general form
K
k k k k k Vn n
V'D k(r)v kk(r)+k (r)k (r)X (I, ),n(r)
V*D () ( +Ea f
n=l
K
k k n k k
+ (r)E k (r) E(nk,r)n ( )+C(r) (r) = 0. (1
n=l
The group fluxes, k(r), are required to satisfy homogeneous
boundary conditions and to be continuous, along with their
k k
net currents, D (r)V4 (r), at internal material interfaces.
The multigroup equations are not selfadjoint. The adjoint
fluxes to the critical system satisfy a set of K equations
of the form
K
V*D k(r)Vk + (r)+E (r) () f) n X"y ( )
n=l
K
k k n( k k
+ER(r) (r) Es(kn,r)) (r)+Z (r)k) (r) = 0 (15)
R p
n=l
and the adjoint boundary and interface conditions. Chapter
IV will discuss some of the eigenvalue equations associated
with equations (14) and (15). Orthogonality properties
among the eigenfluxes and their adjoints will be demonstrated
in that chapter for the sets of interest.
The powerflux relationship is given by
K
P(r) = A E(r) (r). (16)
k=l
This equation contains the crux of the problem for the full
treatment of the multigroup formulation. Unlike equation
(2) for the onegroup case, specification of an optimum
power shape does not allow the solution for the multigroup
fluxes. This occurs because there is no a prior optimum
flux spectrum.
To understand the constraint on the energy dependence
of the optimum flux, it is instructive to first proceed as
if it were given. If this were true, each group flux could
be written as the product of a known function and one of
the other fluxes. These relationships have the form
k(F) = k(r)(r), for k = 1 to (K1). (17)
For convenience in notation, the identity relation for the
K
group K flux may be used to define S (r) = 1. Equations
(16) and (17) then form a set of K expressions in K unknown:
flux functions. The equations are linear, and the solutions
for the optimum k(r) would be straightforward, yielding
Sk(r)P*(r)
k(F) = K (18)
A2 n s rn (F)
n=l
(The fact that this could lead to a physically unrealizable
flux, as in equation (3), Chapter II, is ignored here for
the moment.) Then, upon substitution of these functions
into equations (14), all terms would be known except the
Z (r), for which solutions would again be straightforward.
p
These solutions would have the form
K
n (r)
S (r) P* (r) S (r) z (r)
k n=
Z (r) = VD (r)V K
p
n n k
S n(V) (r)n rS k (r)P*(r)
n= 1
K
+xk S n(r),,n()/Sk (r)f ka(r)ER)
n=l
K
n k
+ L Es(nk,r)S (r)/S (r ). (19)
n= 1
The solution, then, for some arbitrary flux spectrum, would
be a different specification of the macroscopic poison
cross section function in each group, with no a priori rela
tionship among them.
To achieve such a distribution of macroscopic poison
cross sections physically, however, would be quite difficult.
It is very unlikely that there exists a set of K poison
materials, each with a finite cross section in only one of
the energy groups, that could be distributed according to
each of the equations (19). More likely, each position, r,
would require a different poison material, probably a mix
ture, with an energy dependence for its cross section which
would match that specified by all of the equations (19) for
point r.
Except for some problems in breeder reactor design,
where it is important to breeding, there is no desire to
specify the neutron energy spectrum. Generally, however, it
is desirable to use only one, preselected material for a
control poison. In that case, the energy dependence of the
k F
absorption cross section is known, and the Z (r) of equa
tions (19) must be constrained to
k k
E (r) = N(r)a (20)
p P
Thus, the real variables of the problem must be considered
to be the S (r) unknown spectrum functions, of which there
are (Kl), and the N(r) poison distribution function. This
gives a total of K unknown functions.
The K relationships among these unknowns are derived
by making the substitutions from equations (20) into equa
tions (19).
K
n n nF)
Sk (r)P*(r) A Sn ( ) rf
N k k n=l
N(r)a = V*D (r)V K
p K
A S (r)E () S (r)P*(r)
n=l
K
k n"v n k k k
+X k sn (F) vE (r)/S (r) a R Ek (r)
n=l
K
+ s (nk,r) Sn (r)/S (r). (21)
n=l
The problem is thus fully specified, leaving no arbitrarily
chosen variables. However, because the K equations (21) are
nonlinear in the variable functions, the set of solutions,
although finite, may contain more than one member. In fact,
experience has shown there will be K solutions, all but one
of which will involve some group fluxes with negative magni
tudes. Although the physical realities of the problem allow
the selection of the proper solution (the one with all posi
tive group fluxes), these constraints may only be expressed
mathematically as inequalities.
Sk(r) > 0 for all k. (22)
Unfortunately, such inequalities are not suitable for
substitution into equations (21) to analytically produce a
linear set of equations. Consequently, equations (21) are
not solved directly.
Considering now the problem that equations (16) and
k
(18) will not be valid for arbitrary P(r) and f(r) func
tions, it is best to adopt a viewpoint somewhat different
than that of Jackson and Chapter I. If P*(r) in equations
(18) is assumed to be chosen such that it is the optimum
function that satisfies equation (16), while allowing the
fluxes to satisfy their continuity conditions, then it is
just the function P* sought by Jackson. Equations (18)
through (21) are all valid if P* is assumed to be this
optimum realizable power distribution.
Thus, the complexity added to the problem, in going
from onegroup to multigroup theory, is that a nonlinearity
has been introduced. In addition to determining P* from
P*, equations (21) must also be solved.
Unfortunately, this author has been unable to find an
eigenfunction set capable of providing an analytical solu
tion to the coupled, nonlinear requirements of determining
the optimum realizable power distribution and the associated
distribution of a previously selected control material. How
ever, it will be shown that proper choice of an eigenfunction
expansion technique can lend considerable power to the search
for a solution.
CHAPTER IV
SELECTION OF A MULTIGROUP EIGENFUNCTION SET
In order to find a useful set of eigenfunctions, it
is first necessary to find some linear portion of the prob
lem suitable for expansion. Because of their nonlinearity,
eigenfunctions associated with equations (21) would be no
easier to find than the direct solution of these equations.
Thus, the "natural modes" or Kaplan(11) for this particular
problem are impractical for expansion purposes because they
also are given by a nonlinear equation.
The flux portion of the problem, as given by equations
(14), is linear in the unknown fluxes, k (r). It is non
k k
linear overall only due to the control terms, E (r)k (r),
because of the unknown distributional function in the poison
cross sections. Eigenfunctions to expressions associated
with equations (14) can thus readily be found if the control
terms are left out of the eigenvalue equations.
Of the several possible sets of multigroup flux eigen
functions, some are more useful than others. As was the
case in the onegroup formulism, the "natural modes" to the
unpoisoned flux problem contain more information than is
needed, or wanted, about the flux in unfueled portions of
the reactor. As will be demonstrated, the "fission eigen
functions" of Jackson,() on the other hand, contain a
degeneracy in the energy variable, and are suitable only
for expansions of the spatial dimension. Other eigenfunc
tions will be introduced which overcome both of these ob
jections. Before this is done, however, the author will
introduce the multigroup, finite difference formulism
that must generally be employed to calculate whatever eigen
functions are to be used. This will allow the comparison
of the various eigenfunction sets within the actual con
text of their computation and use, helping to more clearly
exhibit their various strengths and failings.
The finite difference approximation to the spatial
dimension is quite analogous to the multigroup treatment of
the energy dimension. A finite number of points, r., are
chosen along the ordinates of the reactor, and the con
tinuousfunctions of flux, power, etc., are represented
by their values at these points. Thus, if J points are
chosen, an unknown function becomes a set of J unknown
values, usually written together as a vector. Applying
this to the group fluxes of equations (14), the flux in
group k would be represented as
= kk k
k(r) {k(r=r ); ( r2); ... k(r=rj)}. (23)
In order to simplify notation for the finite difference
equations, the following convention will be used. Any func
tion, f, of energy and space will have its value at point
j in energy group k denoted as f(k,j).
Equations (14) may be rewritten for each spatial posi
tion, r(j), for which a flux value is to be determined.
The second derivatives associated with the gradient terms
may be approximated by the second central difference
method. 16) For example, in onedimensional, slab geometry,
the gradient term becomes
k k d k (X)d k
V'D (r)V.' (r) = {D (x)14 (x)} (24)
which, for constant diffusion coefficient and evenly spaced
x(j), would be approximated as
kd2 k(x) 2
D 2 D(k) [(k,j+l)2 (k,j)+c (k,jl)]/Ax2. (25)
dx
From this it can be seen that in one dimension the set of
J finite difference equations, replacing the kth multigroup
equation (14), each contain three of the flux values in
the kth group and are therefore coupled spatially. Thus
the finite difference representation of the criticality
conditions will consist of J.K equations of the general form
C(k,j,j+l) (k,j+l)+C(k,j,j)c (k,j)
+C(k,j,j1) (kjl)+Za(k,j)) (k,j)
K
Esz(n 3k,j) (n,j)+Zp(k,j) (k,j)
n=l
K
X(k) I vZf(n,j) (n,j)+ZR(k,j) (k,j) = 0, (26)
n=l
where the C's will depend upon choice of geometry, spatial
variation of the diffusion coefficient, and spacing of the
evaluation points, r(j). In analogy to equation (20), the
finite difference representation of the previously chosen
control material is
Sp(k,j) = N(j)op(k), (27)
where the a (k) is a set of known values, Thus, the finite
difference, multigroup approximation to the criticality
condition, as represented by equations (26), has J.(K1)
unknowns: the J.K fluxes, p(k,j), which are coupled in
energy and space, and the J poison concentrations, N(j),
which are coupled to the fluxes of point j at all energies.
The real flux in the reactor is now represented by
only J*K values, or as a vector with that many elements.
An arbitrary flux vector will thus have J*K degrees of
freedom, and, if it is to be represented by a sum of known
vectors, this sum must also have that many degrees of free
dom in the unknown coefficients. Stated in terms of vector
analysis, a set of J.K linearly independent vectors is
necessary to span the space of all vectors with J.K ele
ments.(16,17) Therefore, in order to expand the multi
group flux vector with a set of eigenvectors, it is neces
sary that there be enough of them to span the desired vector
space. If only H of the spatial points are of interest to
the power specification problem, then only H*K flux values
are of interest. This H*K subspace of the general J.K
vector space may be spanned by only HK of the vectors,
provided a set is chosen that is still linearly independent
when only the H*K elements of interest are considered. This
is in contrast to the continuous function formulism, where
an infinite number of functions, forming a complete set
over some space, are required to fully construct an arbi
trary function over that space. Because in the finite
difference case the physical representation of the problem
was truncated at a finite number of points and energies,
a finite number of vectors is sufficient to construct any
vector representing the flux values at each point and energy.
If the control terms of equations (26) are neglected,
and the remaining terms are rearranged into various eigen
value problems, the resulting J*K equations will be linear
in the J*K unknown fluxes. The sets of solutions to these
problems may be found by well established methods,1819
whichwillbe discussed in the next chapter. These problems
may all be written in matrix notation in the form
[A]
where a. is the generalized eigenvalue. The actual form of
1
the matrix operators, [A] and [B], will depend upon the
choice of eigenvalue problem. These matrices will be square,
of order (J'K), and the eigenvectors, i, will have (J.K)
elements.
Before considering particular eigenvalue problems,
some observations about the general formulation are in order.
The adjoint to equation (28) is given by
[A]T = ai[B]T i. (29)
Following Hansen and Clark,(20) the orthogonality relation
ship may be derived as follows. Assuming that [A] has an
inverse, then the "solution" to equation (28) may be written
as
(i/ei = [A]1 B]i, (30)
and the "solution" for the adjoint problem is
i /ai = {[A]T} [B]T4i = {[A]1 } (B]T (31)
The adjoint to equation (30) is
ei/ai = {[A]1 B] }Ti = (B] T(A]1 } .. (32)
Since matrix multiplication is not commutative, in general,
the adjoint to the "solution" need not be the same as the
"solution" to the adjoint equation. Thus, in the terminology
of Hansen and Clark, (20) the operations of adjoining and
"solving" are not generally commutative. Use of this fact
may be made in selecting a set of adjoint vectors with de
sirable orthogonality conditions. The 8's satisfy the usual,
unweighted condition e4 = 6(i,j), and are thus orthogonal
only over the whole J*K vector space. A more useful ortho
gonality condition may be derived for the 4j adjoint set.
If the operator [B] is applied to both sides of equation (30),
the result is
[Bli/ai = [B] [A] [B] i. (33)
Now, if [B](i is considered as a vector, its adjoint, de
noted by wi, must satisfy the equation
T T
wi/ i= {[A] } [B] wi. (34)
By comparison of equations (31) and (34), it can be seen
that wi is identical to i By taking the inner product of
the transpose of equation (34) with the vector [B] ., and
T
subtracting the inner product of wi with equation (33), the
result is
T T 1
= w [B] [A] [B] jw [B] [A] ([B]i = 0. (35)
With the scalar a's removed from the inner products, this
becomes
T
(l/a1l/a.)w.[B]j. = 0. (36)
Substituting .i for w., this gives the orthogonality condi
tion
T(37)
i[B j = 6(i,j). (37)
Another property of interest in selecting a set of
eigenvectors is the number of eigenvectors associated with
a particular formulation. This can be found by considering
the determinate of [A]a[B] as a function of a. In order to
satisfy equation (28), a must satisfy
det{[A]a[B] } = 0. (38)
In general, expanded analytically, equation (38) results in
an equation in a of degree equal to the order of the square
matrices [A] and [B]. The coefficient of the an term will
be equal to the sum of the determinates of all possible
matrices formed by replacing n columns (or rows) of the [A]
matrix with their corresponding columns (or rows) of the
[B] matrix. It has already been specified that [A] must have
an inverse, so det[A] is not zero. However, no such restric
tion is placed on [B]. The rank of a matrix is equal to
the number of rows (or columns) in the matrix that are linear
ly independent. If [B] is of rank m, then all matrices
formed by replacing rows of [A] by more than m corresponding
rows of [B] must contain linearly dependent rows, and, there
fore, have zero determinates. The coefficients of all terms
a, where n>m, will therefore be zero. Equation (38)
will thus be an m degree equation in a and have no more than
m solutions. The number of eigenvalues (and eigenvectors)
to equation(38) will therefore be limited to being less than
or equal to the rank of the matrix [B].
Except for the "t modes" to the Laplace transformed
kinetics problem, where they indicate oscillatory behavior
in time, complex or imaginary solutions for a seem to have
no convenient physical interpretation. In practice, how
ever, such roots do not seem to occur for formulations
where the eigenvalue may be considered as the multiplier of
some real physical parameter. Double or higher order re
peated real roots may occur in special cases of some of
the eigenvalue formulations to be considered herein. In
such cases, it is not always possible to find a complete
set of eigenvectors. If not, methods exist for treating
the problem by supplementing the eigenvectors with a set
of "principal vectors". (16) Such special cases will not
be considered further in this paper. Hereafter, it will
be assumed that, for matrix (B] with rank m, only the cases
for which there are m true eigenvectors are considered.
Considering now the flux "natural modes" associated
with the criticality equations (28) with the control terms
removed, the individual equations are of the form
C(k,j,j+l) ji(k,j+l)+C(k,j,j) i(k,j)
+C(k,j,jl)i (k, j1)+ a (k,j) q (k,j)
K
X(k) Z (nJ)(n,j) (n,j)+ER(k,j)(i(k,j)
n=l
K
Es(nk,j) i(n,j) = Yi4i(k,j), (39)
n=l
where y is the eigenvalue for these modes. In matrix no
tation, this is
[A](i = yi[I] i.. (40)
All of the materials parameters, with the exceptions of the
control poison, are included in matrix [A], which will have
an inverse unless the reactor is just critical. (That is,
det[A] = 0 is just the criticality condition.) The iden
tity matrix, [I], is its own inverse. Thus, for a J point,
K energy group representation of the reactor, there will
be J*K eigenvalues and eigenvectors satisfying the ortho
gonality conditions
[I] i = i = 6(i,j). (41)
These J*K linearly independent eigenvectors (assuming no
degeneracy of roots in [A]) are sufficient to expand any
arbitrary flux vector of the chosen mathematical model.
Still, however, these J*K vectors may contain more infor
mation than is of interest or more than can conveniently
be used. If only H of the J spatial points lie within
fueled regions, where the power is finite, then only HK
flux values need be specified. It would be convenient if
the rest could be neglected. However, these eigenvectors,
with JK elements, satisfy the orthogonality conditions
with their adjoints only over the entire JK range. Thus,
J*K possible variables (the expansion coefficients for all
the eigenvectors) have been created to handle a case where
there are at most H*K degrees of freedom of any interest.
Experience with the onegroup problem suggests the
selection of the "fission modes" as a set of eigenvectors
that would be complete over the fueled region alone. The
individual equations to this formulation are
C(k,j,j+l)pi(k,j+l)+C(k,j,j) i(k,j)
+C(k,j,j1) i(k,jl)+Ea(k,j) i(k,j)
K
1 1s (nk,j)i(n,j) +R(k,j)i(k,j)
n=l
K
= i X(k) IVf (n,j)f(n,j), (42)
n=l
where Xi, the fission mode eigenvalue, is equal to the
inverse effective reactivity of the mode i. The ortho
gonality condition leads to an expression for the expan
sion coefficients for a given flux vector, 0, given by
ai = 9i[B]4, (43)
that weights any fluxes in unfueled portions of the reac
tor by zero. Just as in the onegroup case, only the fluxes
in the fueled regions need to specified to determine the
expansion coefficients. However, the [B] matrix to this
problem will have J sets of linearly dependent rows, K rows
per set. This is due to a physical degeneracy in the energy
dimension of the fission process: the energy spectrum of
the neutrons born in fission is virtually independent of the
energy spectrum of the neutrons causing the fissions.
Thus, as can be seen by the right hand side of equation (42),
the births due to all fissions at each point are summed
and then distributed over the energy groups at the same
point by the function X(k). Physically speaking, the left
hand sides of equations (42) are being "driven" by a neu
tron source with a constant spectrum, X(E), and this
driving source spectrum is independent of the flux spectrum.
The resulting [B] matrix has the form shown in Figure (1).
As can be seen in the figure, the K rows associated with
the K equations for a given point will differ only by
factors of X(k). Thus, all K of the rows for that point
will be linearly dependent. The rows for different space
points are not linearly dependent, however. The rank of
matrix [B] for the "fission modes" is then J, the number of
space points. This set of eigenvectors,although orthogonal
over the proper region of the reactor, is not complete. The
J vectors are insufficient to expand the H.j elements of
interest in the flux vector.
A set of eigenfunctions that is complete over the de
sired space can be generated by equations of the form
C(k,j,j+l) i(k,j+l)+C(k,j,j) i(k,j)
.+C(k,j,j1) i(k,j1)+Ea(k,j ) (k,j)
K
X(k) vE f(n,j) (n,j)+ER(k,j)( (k,j)
n=l
K
1 Es(nk,j){i(n,j) = 8io (k,j) i(k,j), (44)
n=l
where a (k,j) is the poison cross section in the kth group
(for some unit density of poison) for points where control
poison is to be used, and is zero where no control poison
X(2)vCf (1,1)
X(2)vZf (1,1)
X(1) V f (2,1)
X(2)vZf (2,1)
S.. X(l)v f(l,2)
... X(2)v f (1,2)
Figure 1
The [B] Matrix and Flux Vector for the "Fission Eigenfunctions"
The rows of the matrix represent equations for the fluxes in the
first and second energy groups at the first point and the first
and second energy groups at the second point, respectively. The
columns (and thus the fluxes) correspond to the first and second
energy groups at the first point and the first and second energy
groups at the second point, respectively.
c (1,1)
((2,1)
c(1,2)
4(2,2)
X(1) ~f (2,2) ...
X(2)vZf(2,2) ...
is to be used. Physically, this is equivalent to specifying
that the concentration of poison is the same at each con
trol point. The eigenvalue, Bi, varies the magnitude of
all the concentrations to achieve criticality in each mode.
In this case, which the author choses to call the "poison
eigenfunctions", matrix [B] will be diagonal, with zeros
on the diagonal for all the rows corresponding to points
that have no control poison. For a problem with G control
points, matrix [B] will be of rank G*K, and there will be
G*K eigenvectors which are complete in energy and space
only for the regions of the reactor chosen for control. In
matrix form, equations (44) are written as
[A"l i = !i[op]i', (45)
where [a ] is the [B] matrix for the poison eigenvectors.
These eigenvectors still have J*K elements, and so re
present the flux at the noncontrol points as well. The
criticality condition limits the flux at noncontrol points
to values which are representable by sums of the poison
eigenvectors. The proof that the poison eigenvectors span
the "space" of critical configurations is quite straight
forward. Consider the finite difference criticality equa
tion for the controlled reactor,
[A]> = [Z p] (46)
The matrix [E ] is diagonal and represents any critical
poison distribution, generally with different concentra
tions at each control point. Like [ p], this critical con
trol matrix, [Z ], will have zeros on the diagonal for all
noncontrol points. Assuming that a complete set of poison
eigenvectors is available, the expansion coefficients for
the critical flux may be found from the orthogonality rela
tions.
T
ai = i[Op] (47)
The flux at the control points, at least, is represented
by
I
= aii (48)
i=l
Now, if p* is different from p, it must differ only at non
control points, since the poison eigenvectors are complete
over the control points. Let any difference be denoted as
8 = (p *. (49)
It is easily seen that the constructed flux satisfies the
criticality equation (46) by using the eigenvalue relations
(45).
I I
[A]1* = [A] aii =[ap] aijii (50)
i=l i=l
Since the matrices [a ] and [Z 1 have zeros on the diagonal
for all noncontrol points, the substitution of c* in the
criticality equation yields a relationship dependent only
upon the fluxes at control points, where }* and are iden
tical.
I I
[op] aiiP i = Ep aii (51)
i=l i=l
Thus, 4* also satisfies the criticality relation. Now, by
substituting (*0 for in equation (46),
[A] (4*6) = [Zp] (*6). (52)
Since p* already satisfies equation (46), this may be
reduced to
[A] = [Z ] (53)
However, 8 must be zero at all control points, and the
operator [Z ] multiplies all noncontrol points by zero, so
p
[Z p] e = (54)
Thus, if the difference vector, e, exists, it must satisfy
[A]8 0 (55)
So, for a nontrivial solution for 6 to exist, it is necessary
that
det[A] H 0. (56)
If this condition is not satisfied, then 6 is the zero vec
tor and the expanded flux is identical to the actual flux.
If equation (56) is satisfied, the uncontrolled reactor
must be critical in some mode, i. For this case, however,
a poison eigenvector will exist with a zero poison eigen
value.
[A] = 0[op ]( 0 (57)
Thus, if it existed, 6 must be identical to .i, and thus
already included in the expansion of the flux, in which
case it could not be the difference between the expanded and
the true flux. Thus, in all cases where there is a complete
set of poison eigenvectors, they are capable of representing
any critical reactor configuration attainable by distribu
ting poison at the designated control points.
These poison eigenfunctions at last seem to be able
to adequately represent the physical situation. For a reac
tor with G control points (e.g., rods, regions, etc.),
there are only G degrees of freedom for control purposes.
Another G (K1) variables are needed for expansion, however,
because the spectrum is not known a priori. This set of
eigenvectors provides the minimum number of necessary vari
ables and insures that they are complete over only those
spatial points that may actually be directly controlled.
CHAPTER V
GENERATION OF THE POISON ENGENVECTORS
Before the eigenvectors to a particular problem can be
used in the solution of that problem, the eigenvectors them
selves must be obtained. Most computer codes available today
for the solution of the diffusion equations are concerned
only with the fundamental mode of the fission eigenfunctions.
That is, they find the eigenvalue that is the inverse of the
reactivity for the real flux shape in the reactor. Such
codes are designed to find only the one mode for only the
one eigenvalue formulation, which is degenerate in the energy
dimension. Thus, there is little hope for the modification
of this type of code to find the complete set of poison
eigenfunction.
(14)
An exception is the MULE code,(14) which is capable of
finding any desired mode for three types of twogroup prob
lems, the fission eigenfunctions, the poison eigenfunctions
for a thermal poison only, and the time or kinetic eigen
functions, with or without delayed neutrons. The fission
modes satisfy the equation
[L]bi = Xi[S]b i, (58)
where [L] contains all of the leakage and absorption loss
terms and [S] contains the fission source terms. As dis
cussed in Chapter IV, the [S] matrix has a rank equal to
the number of spatial points, J, and thus equation (58)
produces only J eigenvectors. The thermal poison eigen
functions have a similar degeneracy. They satisfy the
relation
{[L] K [S]i = pi [Q]i (59)
eff
where [Q] is a diagonal matrix with zeros on the diagonal
for all the fast group equations and E on the diagonal for
all the thermal group equations. Thus, the rank of [Q] is
also J, and this set of eigenvectors is also insufficient to
fully expand the twogroup problem. The time modes, without
delayed neutron groups, satisfy the equation
{[L] K [S]}i = w l[V]a, (60)
eff
where [V] is a diagonal matrix with the inverse of the neu
tron velocity for the appropriate group on each diagonal
location. The rank of [V] is thus 2J, and the eigenfunc
tions of this type will be complete over energy and space.
This problem is very similar to the twogroup poison eigen
functions desired. If the poison fast and thermal cross
sections were substituted for the inverse velocities in the
fast and thermal groups, respectively, then true twogroup
poison eigenfunctions could be found. However, as written,
the MULE code would impose the unwanted restriction that all
points would be control points. This problem could easily
be overcome by writing an additional input routine to
allow the desired flexibility in control point locations.
Instead of altering MULE, this author has chosen to
write a new eigenvector generator that can be run on a
smaller computer, allowing the programmeroperator to have
more convenient interaction with the program while calcula
tions are in progress. The desirability of such an arrange
ment will become apparent shortly.
The actual calculation of an eigenvector is typically
accomplished by the power iteration method.18) For the
general eigenvalue problem of equation (28),
[A]ji = ai[B]J i, (28)
there will be a set of I solutions, where I is less than or
equal to the order of matrices [A] and [B]. These vectors
can be combined in any arbitrary manner to form a vector
I
V = a i (61)
i=l
If [A] has an inverse, then the eigenvalue problem may be
written in the form of equation (30).
(i/ai = [A]l1B]i (30)
The effect of operating on vector V by the matrix [A]1 [B]
can be found by substitution for each eigenvector from
equation (30).
I I
[A] [B] ai = aii/i (62)
i=l i=l
The vector on the right hand side of equation (62) has thus
had its component of the ith eigenvector changed by the
factor 1/i If ai is greater than one, the contribution
from vector i will be decreased; if it is less than one,
the contribution will be increased. Therefore, by repeatedly
operating on any vector V by the matrix [A]1 B], the eigen
vector component with the smallest absolute value of a will
be the one to grow the fastest or decrease the slowest in
absolute value. If the new vector is renormalized after
each successive operation, then eventually, only this domin
ant eigenvector will remain. As can be seen from equation
(30), when the guess vector has converged to the eigenvector
i, the renormalization factor will be just the inverse of
the associated eigenvalue. How fast the other vectors die
away is governed by the dominance ratio. If vector i has
the smallest value of a, then the ratio of increase in the
contribution of vector i, relative to vector j, between
iterations is given from equation (62) as
aii / a j i (63)
a $ a a
The greater the dominance ratio, the more rapidly vector j
will die away compared to vector i.
Obviously, if matrix [B] also has an inverse, it is
possible to find the largest eigenvalue and associated
eigenvector by writing the problem as
a = [B] [1A]i, (64)
i
and operating repeatedly with the matrix [B] [A]. Still,
the vectors with eigenvalues of intermediate modulus must
be found by some other method. There are two possible
methods. The first is successive elimination(18) of the
roots as they are found. This can be done by choosing a
starting vector that is orthogonal to all previously found
eigenvectors, so that the sumin equation (61) no longer
contains these components, or by modifying the original
matrix(18) to eliminate known roots. This leaves the next
root dominant. There are several disadvantages to such
methods. The dominance ratios for the last modes are usually
not very great, causing many iterations to be necessary for
convergence to the next mode. Accuracy of the later modes
also suffers because the calculations depend upon the results
for all previous modes, and the errors accumulate.
A better method, proposed by Wielandt1819) is to
transform the eigenvalue problem to one in which the desired
root is dominant. If a guess is made for the next value of
a, then the eigenvalue problem may be rewritten as
{[A]a [B]}( i = 6i[B] ,i'
(65)
where ag is the value of the guess, and the eigenvalues of
the transformed equation (64) are related to the original
eigenvalues by
a. = 6L.+ (66)
The dominance ratio for the new equation is calculated using
the 6's. Obviously, the smallest 6 will be the one associated
with the ai closest to the a and the closer the guess,
the faster will be the convergence. In order to find any
eigenvector, it is therefore only necessary to provide a
guess at the eigenvalue which is cloFrr to that eigenvalue
than to any other, and to repeatedly mu tiply some starting
vector by the matrix {[A]a [B]} [B] until convergence is
obtained. Of course, standard acceleration technique may
be used to speed the convergence of the vector. The Wielandt
transformation itself may be used as an acceleration technique
if the eigenvalue guess is updated using partially converged
values of 6.(21,22)
To aid in guessing close to the eigenvalues, it is
helpful to consider the physical significance of the eigen
value problem chosen. The poison eigenfunctions satisfy
equation (45).
[A]bi = i[o p]i (45)
The [A] matrix contains all of the terms for the uncontrolled
reactor, and matrix [o ] contains the terms for absorption
in unit concentrations of poison placed at all control points.
The eigenvalue serves the function of uniformly varying the
amount of poison at all the control points. If the reactor
is supercritical without control poison, then the funda
mental mode (the most reactive one with flux that is every
where positive) will require some control poison insertion
to reach the critical state. Thus, the associated 81 will
be positive. A good estimate of its value can be made from
first order perturbation theory, if the fundamental fission
mode and its adjoint have already been found by one of the
standard diffusion codes. This is usually not necessary,
however.
In order to understand how the rest of the modes be
have, it is instructive to examine a simple example. For
a twogroup representation of a bare, homogeneous cylindri
cal reactor of infinite axial extent, the equations for the
twogroup fluxes in the poison eigenvalue formulation are
given by
D(F)V20i(F,r)+ a ()++ R(F)vZf(F)i(F,r)
v f(T)(i(T,r) = 8iop(F) i(F,r) (67a)
for the fast group, and
D(V2T)V2 aT)i (T a (T,r)ER (F) }.(F,r)
= 8io (T)c (T,r) (67b)
for the thermal group. Because the energy and space variables
are separable for this example, the solutions to equations
(67) may be written as fluxes of the form
i(F,r) = SiJo(nir/R) (68a)
i(T,r) = Jo(nir/R) (68b)
where R is the radius of the reactor, J0 is the zeroth order
Bessel function, ).W is the ith zero of J and Si is the
fast to thermal flux ratio, which is a constant over the
entire reactor volume for each mode. Substitution of equa
tions (68) into equation (67) allows the solution for the
eigenvalue, 6i, and the spectral constant, Si, for each
mode. These solutions are
b b c
S + (69)
i 2a 2a a
i = {SiE (F)D(T)[ ni/R]2_ a (T)}/op(T), (70)
where a = p (F) R(F)
b = 0 (T){D(F) ["./R]2+Z (F)+ZR(F)}
p 1 a +ZR
0 (F) {D(T) [r./R] +a (T)}
}p }
c = vZf(T)o (T).
The quadratic nature of the solution is due to the twogroup
formulism. In general, the degree of the equation for Si
will be equal to the number of energy groups.
The J function is a well known oscillatory function
which crosses the axis at the tabulated values, n..
1
By scaling its argument as (n.r/R), it is made to have zero
value at the extrapolated reactor boundary, R, and to pass
through zero (i1) times within the reactor. For a parti
cular choice of the index, i, the fast and thermal fluxes
will have the same spatial shape, differing only in magnitude
by the factor Si. Because there are two groups, there are
two values of 6 and two values of S associated with each
of these shapes. From equation (69), it may be seen that
one value of S. will always be positive, and the other
always negative. By equation (70), the negative Si will
always be associated with a negative 8B, while the positive
Si may have either a positive or negative value of Bi. In
fact, as discussed previously, when the uncontrolled reactor
is supercritical, the fundamental mode must have a positive
value of 8. As the index, i, and thus ni, is increased in
value, the leakage terms in equations (67), given by
D(k)V J (nir/R) = D(k) [nir/R] 2J(nr/R), (71)
will become more positive. In order to satisfy equations
(67), then, the values of 6 must become more negative. The
same conclusion may also be reached by physical reasoning.
The innermost region of any higher flux mode, where
Orn.R/n., is physically equivalent to the fundamental mode
of a smaller reactor. In order to make this smaller reac
tor critical, less control poison will be needed because
of the greater neutron leakage. Thus, again, as ni increases,
Bi must decrease.
For this case, then, the B's will occur in two sequences,
denoted by B and 8T for the values associated with the posi
tive and negative values of Si, respectively. These sequences
satisfy the relations
+ +
Bi+1 < i (72a)
Bi+l < Bi < 0. (72b)
From the shape of the fluxes, the index, i, may be derived
by counting the number of times the flux passes through zero
and adding one.
In general, the flux shapes of reactors with more com
plicated materials configurations will not be given by the
analytic function, J Over their region of completeness,
they will still behave in a manner similar to that described
in the last paragraph, however. Experience has shown that,
if the flux values associated with noncontrol points are
ignored, the remaining values will show the proper number of
sign changes in both the energy and space variables. This
is of great utility to the person wishing to generate a
set of eigenvectors using the Wielandt transformation
technique. Any guess at an eigenvalue will result in con
vergence to some eigenvalue and vector. By using the above
criteria, it may easily be determined which eigenvalue has
been found, and where the other eigenvalues must be located
in relation to this and previously found eigenvalues. Thus,
a few educated guesses, or, at worst, a few shots in the dark,
for values of 8 will result in enough information about the
magnitude and spacing of the B sequences so that guesses
near the remaining roots may be made quite readily.
While the evaluations and decisions associated with
this initial searching and later prediction for values of
S may easily be handled by the human mind, programming of
the general process would be quite complicated and require
much computer storage. Such complexities as the possible
overlap of the 8 sequences, variability of the control
(completeness) regions between specific problems, and the
possibility of finding the same root more than once, greatly
complicate the logic and bookkeeping necessary for a compu
ter program that would be capable of finding the complete
set of eigenfunctions with no operator interaction.
For these reasons, and for the purpose of illustration,
the author chose to write a sample eigenvector generator
designed to run on the IBM 1800 computer system. Briefly,
the program reads the problem geometry from punched cards,
sets up the iteration matrix, and asks, via the typewriter,
that a guess be made for 8. When this value is entered,
along with the maximum number of iterations to be allowed
on this guess, the program performs the iterations until
convergence of the vector or the iteration limit is reached.
The revised eigenvalue, the vector and, if convergence is
obtained, the adjoint vector are then displayed to the
operator on the typewriter, and the operator is given the
option of having this output punched on cards for later use.
The program then cycles back and asks for another eigenvalue
estimate. Through the use of sense switches, the operator
also has the options of varying the convergence criteria
according to the estimated numerical stability of the various
modes, looking at any of the intermediate vectors in the
iteration process, or specifying that the adjoint be found
regardless of eigenvector convergence.
This program, named CODE2 by the author, is described
more fully in Appendix I, and a listing appears as Appendix
II. Because the program was written only for purposes of
demonstration, it is much less general and sophisticated
than it is allowed to be in theory. It is limited to two
group, onedimensional calculations in cylindrical geometry.
Extentions to more groups, more dimensions and other geome
tries ate straightforward, following established techniques
employed by typical modern diffusion codes. For the purpose
of demonstrating the concepts and techniques introduced
in this paper, however, the program is sufficient as written.
CHAPTER VI
THE MULTIGROUP OPTIMIZATION PROBLEM
As pointed out in Chapter III, the problem of finding
the optimum multigroup flux functions, and the distribution
of a control poison necessary to maintain these flux shapes
in the critical state, results in a set of nonlinear equa
tions. There is no general analytical approach for solving
nonlinear systems, and the author has been unable to find a
specialized technique capable of providing a solution for
this particular problem. Therefore, this author has chosen
the numerical optimization approach to provide individual
solutions for individual reactors.
Numerical search techniques, in general, find solutions
for particular problems based on the results of sampling
calculations performed with various sets of the independent
variables. As applied to the problem of finding the optimum
poison distribution in a particular reactor, these variables
are the poison concentrations at each control point. An
initial set of poison concentrations is specified at the
beginning, and the resulting power distribution is calculated.
Comparison of this power distribution with the optimum dis
tribution gives the "error" in power resulting from the use
of this poison distribution. By systematically varying the
individual poison concentrations andrecalculating the power
distribution, the derivatives of the error are calculated as
functions of the poison concentrations at each control point.
Various algorithms may then be employed to predict the loca
tion of, or direction to, the poison distribution yielding the
optimum power distribution. This results in a new poison dis
tribution that is closer to the optimum, about which the
calculation of the derivatives is then repeated. This process
continues until no better poison distribution may be found.
The actual formulae for varying the poison concentrations
and calculating the error values are dependent upon which
optimization technique is used. One of these will be examined
in more detail, later, when numerical examples are presented.
All of these optimization techniques, however, will
require repeated evaluations of power distributions in the
same reactor, one for each poison distribution to be inves
tigated. In order to evaluate the effects of a given poison
distribution, it is first necessary to adjust its overall
magnitude to obtain criticality. The flux distribution in
the critical reactor must then be found, and, finally, the
power distribution may be calculated. In this chapter, it
will be shown how the proper choice of eigenvalue formula
tion and an eigenvector expansion technique can greatly aid
in this part of the numerical analysis, the evaluation of the
power distribution resulting from a given poison distribution.
For purposes of comparison, it is first instructive to
review the approach typically used by modern computer codes
written to solve neutron diffusion problems. A good example
of such codes is CORA,(23) which was written in 1970 for
the solution of onedimensional, fewgroup problems. Typi
cally, the eigenvalue problem chosen for solution is that
of the "fission eigenfunctions," and the algorithm used is
capable of finding only the fundamental mode and its adjoint.
The equations are similar to equations (42) except that a
separate, specified control term is included. The equation
for the flux in group k at point j has the form
C(k,j,j+1) (k,j+1)+C(k,j,j) (k,j)+C(k,j,j1) (k,j1)
+Ea (k, ) ( ,(k, (kj)())(k,j)+E (k,j) (k,j) (kj) (k,j)
K
s(nkj)0 i(n,j)
n=l
SK1 X(k) vzE(n,j)(n,j), (73)
eff n
n=l
where e is a search parameter to be explained below. The
solution of these equations is obtained by the power iteration
technique, as discussed in Chapter V. A Chebyshev polynomial
extrapolation technique is optionally available for ac
celerating the convergence of the flux vector. Along with
the fluxes, the reactivity is found as the inverse of the
eigenvalue. A search routine is available to automatically
vary the control poison concentrations to achieve criticality.
This is accomplished by repeated solution of the fission
eigenvalue problem with all the control terms modified by
a single search parameter. This parameter, denoted as E in
equations (73), is varied between successive flux solutions
according to a hyperbolic extrapolation algorithm(23) in
order to reach a reactor configuration with the value of
reactivity specified. Thus, the user picks the shape of the
poison distribution, and the program gives the necessary mag
nitude.
This author would like to point out that, by selecting
the poison eigenvalue formulation, instead, this same prob
lem can be handled with only one solution of the flux equa
tions. In order to include the option of a user specified
reactivity other than unity, equations (44) for the poison
eigenfluxes may be written as
j+l)+ ,j+l) (k,j+l)+C(k,j,j))(k,j)
+C(k,j,jl)( (k,jl)+ a (k,j) (k,j)
K
 X(k) vZf (n,j)4 (n,j)+E (k,j)( (k,j)
eff 7 f
n=l
K
Es(n k,j). (n,j) = cE (k,j) (k,j). (74)
n=l
With Keff and the shape of the poison distribution, E (k,j),
specified, equations (74) may be solved by the power itera
tion method, with acceleration techniques if desired. The
solution will be the flux distribution for the reactor con
figuration with the desired reactivity, and the eigenvalue,
c, will give the necessary magnitude of the poison distri
bution. The need for a search for the solution of this
problem may thus be eliminated by the selection of the poison
eigenvalue formulation.
Finding the concentration magnitude necessary for cri
ticality, with a given shape of the poison distribution, is
only part of the overall problem of finding the optimum poi
son distribution to minimize power peaking. As repeatedly
discussed, this total problem is nonlinear when use of a
particular control material is specified. It therefore must
be solved by a numerical search over various poison distri
butions which result in a critical reactor. Several numeri
cal optimization techniques are more or less suitable for
application to this problem. The choice is inmaterial to
the topic of this paper because the methods outlined herein
are suitable for use with all of them. They all require that
an error value, giving the difference between the desired
and actual power shapes, be found for each poison distribu
tion investigated. The specific form of the error depends
upon the method employed and how that method uses previous
results to estimate the desired solution.
For each poison distribution specified by the search
algorithm, the magnitude must be adjusted for criticality,
the fluxes must be found, the power distribution must be
computed and the appropriate error values) must be calculated.
The utility of the poison eigenvalue formulation has already
been discussed. Using this approach, the solution to a J*K
order matrix eigenvalue equation must be found for each poi
son distribution investigated for a reactor represented by
J space points and K energy groups. This matrix equation
has the form
[A]H = E[Z p], (75)
where [A] contains all of the terms on the left hand side
of equations (74), and [E ] contains the cross section terms
for some normalized poison distribution.
The use of the poison eigenvectors to the particular
reactor geometry can aid greatly in the solution of equation
(75). The unknown flux vector may be written in terms of
the known eigenvectors as a function of unknown expansion
coefficients, a..
1
M
S= aii' (76)
i=l
where M is the order of expansion. If M = Gk, where G is
the number of control points, then, as discussed in Chapter
IV, the M poison eigenvectors will be sufficient to repre
sent any flux distribution over the reactor that could result
from an arbitrary poison distribution (with magnitude adjusted
for criticality). Equation (75) may therefore be written
in terms of the poison eigenvectors.
M M
[A] aii = C~ pl] a (77)
i=l i=l
Since the matrices [A] and [Z ] are not functions of the
index of summation, they may be included inside the sums.
Substitution according to the poison eigenvalue relation,
equation (45), allows the left hand side of equation (77)
to be written as
M M
ai[A] = Oaii[apPi. (78)
i=l i=l
The [E ] matrix may be separated into the matrix, [N], giving
just the normalized spatial distribution of the poison con
centrations, and the matrix [ap giving the cross section
terms that would result from unit poison concentrations at
all control points. The right hand side of equation (77)
may then be written as
M M
E[(E] Eaii = E(N][ [p] aii. (79)
i=l i=l
After making these substitutions in equation (77), a set of
M equations may be generated by multiplying both sides of
the equation from the right by the transpose of one of the
M adjoint vectors, m .
M M
i i= E aiT[N] [ap]i (80)
a i m p m
i=l i=l
The orthogonality relations for the poison eigenvectors and
their adjoints are found from equation (37) by substituting
[a ] for [B].
m[ p]i = 6(i,j) (81)
Using these relations, equations (80) may be simplified to
M
am m = E aiF(m,i,[N]), for m = 1 to M, (82)
i=l
where
F(m,i,[N]) = T [N] [cp]i. (83)
m p1
If M = G'K, then solution of the M "component equations" (82)
is completely analogous to the solution of the unexpanded
equations (74). Equations (82) may be written together as a
single matrix equation.
[g]a = :[F]a (84)
where [8] is a diagonal matrix, with all of the Bm on the
diagonal, and [F] is a full M by M matrix with elements given
by equation (83).
The number of control points, G, may be substantially
less than the total number of points, J. Thus, the substi
tution of a G*K order matrix equation for a J'K order equa
tion appears to offer a substantial savings in computational
work for evaluating the power distribution produced by each
poison distribution considered. In addition, [B] is diagonal
and is more easily inverted than is the [A] matrix of the
unexpanded problem. However, the eigenvalue nature of equa
tion (84) presents a problem. Generally, there are as many
as M solutions to the equation, just as there are to the
eigenvalue equation (45) used to generate the M poison eigen
functions. As pointed out in Chapter V, a simple power
iteration solution of the equation
a/ = [ 1[F]a (85)
will converge to whatever mode has the e with absolute
value closest to zero. In general, there is no reason to
expect the fundamental mode to meet this criterion. The
fundamental mode, with the smallest neutron losses due to
leakage, requires the greatest amount of poison to attain
criticality, and therefore the most positive value of e. Ob
viously, this problem also occurs for the unexpanded poison
eigenvalue formulation, as given by equation (75). There
are two approaches to insuring that power iterations converge
to the fundamental mode of equation (84). Both have their
advantages and disadvantages.
The more obvious approach involves the use of the
Wielandt transformation. An initial guess is made for the
proper value of E to produce the fundamental mode for the
first poison distribution to be considered. Then, as the
computer program proceeds in its search from one poison
distribution to the next, the eigenvalue found for the last
solution is used as the guess for the next. If the program
does not drastically change the poison distribution from one
time to the next, the fundamental modes to the two different
reactor configurations should be similar enough, with eigenvalues
that are close enough, that this method may be expected to
work quite well. In practice, no difficulty has been experi
enced by this author from the interference of other modes.
The starting guess for E in the first poison distribution
case could be made in several ways. A very large value of E
could be chosen to insure that it is greater than all the
actual roots. However, the safer this guess is made, the
worse the dominance ratio becomes and the more iterations are
required to reach convergence. Another method is to start
the poison search with a uniform distribution. In this
case, [Z ] = [p ], and the unexpanded poison search equation
(75) becomes identical to the poison eigenfunction equation
(45). Thus, the desired value of E is identical to 81, the
eigenvalue of the first poison eigenfunction. This method
for predicting c, while accurate, does not allow any prior
knowledge or intuition of the optimum to be used in speci
fying the starting distribution for the poison search. A
third method, which the author has found satisfactory, is
to use first order perturbation theory to predict the value
of E for the fundamental mode associated with the desired
starting poison distribution. If the flux for this mode is
assumed to be approximately equal to the flux associated
with a uniform poison distribution, then the first fission
eigenflux may be substituted for the flux in equation (75).
[A] 1 = E[N] (ap]1l
(86)
Using eigenvalue relationship for 1, this becomes
6_ 1 = E[N] [O p]1 (87)
Multiplying both sides by the transpose of the adjoint flux,
T
'1, gives a scalar equation which may be solved for E to
yield
E = T[N] [crp] / = F(l,l,[N])/81. (88)
Since F(1,1,[N]) has already been found for use in the expan
sion equations (84), this method provides a simple means for
generating a reasonably accurate guess at the fundamental
mode e. Although there is no guarantee that a guess generated
by equation (88) will be closest to the maximum positive
value of e, this author has experienced no difficulty in prac
tice.
The disadvantage of the Wielandt transformation approach
is that the transformed equation has the form
{[BIe [F]}a = 6[F]a (89)
Thus, a full M order matrix must be inverted, instead of
the simple diagonal matrix [8], for the power iteration with
the transformed operator {[8]c [F]} [F]. In comparing the
unexpanded and expanded formulisms for complexity of the
power iteration process, note must be made of the special
nature of the [A] matrix in the unexpanded equation (75).
The exact form of [A] will depend upon how many spatial
dimensions and energy groups have been used and upon the de
tails of scattering among the energy groups. In particular,
if only one spatial dimension is considered, the multi
diagonal nature of [A] (and therefore of [A] [E p]) lends
itself to the use of "inner iteration schemes" that take the
place of the matrix inversion.(23,24) The decision of which
is more complicated, the use of such an inner iteration scheme
on a J*K matrix, or the generalized inversion of a G.K(=M)
matrix, must depend upon the relative magnitudes of G and J
and upon the exact nature of the [A] matrix for each parti
cular case. If more than one spatial dimension is considered,
solution of the unexpanded equation becomes much more diffi
cult. The "inner iterations" may be done in each direction
independently, and the directions alternated until convergence
of the entire set of inner iterations is obtained.(16) In
such cases, however, inversion of the GK matrix would cer
tainly seem preferable unless G is very nearly as large as J.
The second method for insuring that the power iteration
method converges to the proper value of c is designed to pre
serve the [B] matrix in a simply invertable form. This may
be accomplished by observing the following phenomenon. Any
eigenvector,
be a solution to the equation
{[A]+Bs[a ]}i = (i8s) a Opi = l[o p]i. (90)
That is, reducing all of the Bi by the amount Bs is equivalent
J> 5
to assuming that the concentration of control poison is already
Bs at all control points. If 0s is made large enough that
the fundamental poison eigenvalue, B1, is negative, then phy
sically, the reactor is shut down. The problem of determin
ing the optimum poison distribution then becomes a matter of
removing poison from what is already there.
The distributions investigated by the search program
will be the complementary distributions necessary to make
the actual distributions equal to 8s at every control point.
That is, for normalized actual distribution [N], and comple
mentary distribution [N'], the relationship between the
magnitudes necessary for criticality of the fundamental mode
is
e[N] [o ] = E' [N'] [p ]+ s[o ] (91)
The flux matrix equation associated with the evaluation of
such complementary distributions, [N'], is
{[A] s [ap] } = c'[N'][o p] (92)
Using the same set of poison eigenvectors and adjoints and
the new eigenvalue relations (90), equation (92) may be ex
panded by a method analogous to that used in obtaining the
expanded equation (84).
[B']a = e' [F']a (93)
where the elements of [F'] are given by analogy to equation
(83) as
F'(m,i, [N']) = m[N'] [o p i. (94)
Thus, the value of e' which corresponds to the fundamental
mode may be varied by varying s .
As long as poison is only removed at each point, there
is no ambiguity. However, if the base level of control
poison, s, is such that some near optimum trial distribution
has values greater at some points, and less at others, then
the complementary distribution to be removed will be posi
tive at some points and negative at others. Such a distri
bution may be interpreted two ways. With a negative e', the
trial distribution is to be subtracted from the uniform dis
tribution, but with a positive e', it is to be added. Since
the trial distribution has both positive and negative ele
ments, for either case, poison will be added at some points
in the reactor and removed at others. If the complementary
distribution chosen for investigation is such that both posi
tive and negative roots for e' exist, then it is impossible
to specify beforehand that the distribution is to be added
or subtracted. The power iteration solution method will
always converge to whichever value of e' is closest to zero.
Computer search programs usually cannot allow for more
than one interpretation of the set of search parameters (in
this case, [N']) chosen for evaluation of the error function.
Therefore, in order to insure that such an ambiguous situa
tion does not arise, it is prudent to chose the value of B
S
so that, for all total poison distributions that may reason
ably be considered, the complementary poison distribution
will be everywhere negative. The search program may then
be expected to restrict itself to trial distributions
that everywhere have the same sign, thus insuring that all
the values of e' are of the opposite sign. (The actual sign
of [N'] itself is of no significance because only the pro
duct, e'[N], must be negative to designate removal of poi
son.) The fact that such a selection of s will result in
all possible root values of e' being of the same sign is
physically obvious. Since the reactor is shutdown by the
uniform poison distribution, Bs, and since poison may only
be removed everywhere or added everywhere by the distribu
tion [N'], obviously, 8', the fundamental mode root must be
of the proper sign to remove poison if that mode is to be
critical. Other roots must correspond to more buckled flux
shapes, and require even more poison to be removed for cri
ticality of their modes. Thus, the other values of e' must
be greater in absolute value and of the same sign as Cl.
This insures that the power iteration solution of equation
(91) will always be the fundamental mode and that the dis
tribution [N'] will always be considered as subtractive.
Unfortunately, equation (91) does not aid in the selec
tion of a value for Bs that will insure behavior for e'
as described above. Since the "worst case" trial distribu
tion, e[N], will not be known before the search program is
run, some other means must be used to select s. Certainly,
Bs must be greater than 81, or the reactor will not be initially
shut down. If 0s is chosen extremely large, then it is also
certain that poison would have to be removed everywhere to
reach the critical state with the optimum power distribution.
However, such large guesses for Bs also serve to increase the
absolute magnitudes of all the 8i. This tends to decrease
the dominance ratio for 81, causing more iterations to be
required for convergence of the flux for each [N'] investi
gated. For the general case, this may be hard to grasp, but
for the special case [N']= [I], it is easily illustrated.
In that case, Ei = Bi, and the subtractive distribution [N] =
[I]. This allows equation (91) to be simplified to
i = Ei+s (95)
Thus, the dominance ratio for the e!, as a function of 6s,
is I (82B)/( 1,s) Since <2 < 1l < s, this dominance
ratio for E1 must decrease as s increases (i.e., as the
reactor is further shut down).
In order to use this method for insuring that the prop
er root is obtained for e in the solution of the expanded
flux equations, it is therefore necessary to make an educated
guess at the best value of Bs. If the value chosen is too
small, the error may be detected when a value of e; with an
improper sign is computed, or when [N'] contains negative
elements. In such an event, it is necessary that the value
of Bs be increased. The progress toward the optimum need
not be lost, however. From equation (91), it is possible
to find the new complementary poison distribution associated
with the new value of Bs, that corresponds to the best dis
tribution found so far using the old value of Bs. In general,
[N] and [ p] will be diagonal with some zero diagonal ele
ments in positions corresponding to noncontrol points. If
these zero rows and columns are eliminated, however, then the
resulting matrices will possess inverses. With the under
standing that this has been done, the solution for the total
poison distribution associated with the best complementary
distribution, (1)[N'(1)], found so far using s (1) is
E[N] = E' (1) [N' (1) ]+ s(l) [I], (96)
1
where the fact that the reduced (p ] is diagonal, and thus
commutative under matrix multiplication, has been used to
simplify the relation. The new complementary distribution,
' (])[N' (2)], must represent this same total distribution.
Equating the two expressions for e[N] and solving yields
E' (2) [N'(2)] = (1) [N'(1)]+{ s(1)B S (2) }[I]. (97)
This new distribution, normalized to [N'(2)], may be used
as the new starting point when the optimization search is
begun again using Bs(2). This process of retransforming
the 8 set and restarting the optimization search could be
done manually or automatically by the search program.
If an unnecessarily large value of Bs is picked, however,
there will be no way of knowing this for certain until the
optimization search is completed. This is too late to
improve convergence speed by decreasing s to increase the
5
dominance ratio. Perhaps a routine could be programmed
that would automatically decrease gs if previously found
values of e' have been consistently large. In lieu of such
a scheme, or perhaps in addition, it would probably be
worthwhile to include a flux convergence acceleration scheme
to help make up for the expected poor dominance ratios.
An example of this type of acceleration scheme is the
Chebyshev polynomial extrapolation method utilized by the
CORA program.
Of the two methods described above for the solution to
the poison distribution evaluation problem, the first, using the
Wielandt transformation, is the easier to program and under
stand. Because of the necessity of inverting the (F] matrix,
however, its utility may be limited. The second method,
involving only the linear modification of the set of poison
eigenvalues, seems to offer much greater computational
power. This method requires complex programming to be made
economical, however.
CHAPTER VII
NUMERICAL CALCULATIONS AND THE PROSPECT
FOR PARTIAL EXPANSIONS
In order to illustrate the use of some of the tech
niques described in previous chapters, their application
to several sample problems will now be presented. This
will serve two purposes. First, some of the principles
previously discussed only in the abstract sense will be
demonstrated by concrete example. These calculations are
then extended to explore the possibility of using less than
the complete set of "component equations" (82) in order
to approximate the solution with less work.
The principles of eigenvector generating programs
were discussed in Chapter V and Appendix I, and a sample
program is listed in Appendix II. The optimization search,
as discussed in Chapter VI, is carried out by the use of a
specially modified version of the GAUSS computer code written
by Kylstra* at the University of Florida. This program is
(25)
based on the least squared error technique of Gauss,(25) the
principles of which are described in Appendix III. The
eigenvector expansion technique, which is used for synthesi
zing the flux and power distributions for each trial poison
*C.D. Kylstra, personal communication.
distribution generated by GAUSS, is incorporated in the sub
routine UFUNT. Listings of this subroutine and the auxiliary
input and output subroutines are given in Appendix IV. The
value of the error function associated with each trial poi
son concentration is calculated by the algorithm
ERROR = {lP(j) }2V(j), (98)
where P(j) is the (normalized) power density at point j and
V(j) is the volume over which this power density applies.
The sum runs over all points j in fueled regions. It is
this sum that is minimized by the GAUSS program. In cases
where the minimum error is not zero, it may be seen from
equation (98) that the solution given is the optimum in
the least squared errors sense. Other optimization tech
niques, requiring different formulations of the error func
tion, could be used if desired. The Gaussian technique
was chosen because of its reputation for providing solutions
in complicated cases where other methods may fail, and
due to its availability to this author in a generalized
computer code.
Because of its relatively simple programming require
ments, the Weilandt transformation approach has been used
to ensure that the proper E root is found in solving the
expanded equation (84). Although, as discussed in Chapter
VI, this method may be more limited in its utility than is
the Bs transformation method, it is perfectly adequate for
5
the purpose of demonstrating the eigenvalue expansion tech
nique. In order to gain a slight simplification in the neces
sary matrix algebra, the equation solved by the UFUNT sub
routine has been rewritten from equation (84).
[F]a = (1/c) ([]a (99)
Using the Wielandt transformation with a guess at the inverse
of the desired eigenvalue, (1/E), this equation may be
written as
([F](1/)G[] }aG = 6[(]a, (100)
where 6 = (1/E)(1/E)G. The matrix operator used for itera
tion purposes is thus {[F](l/E)G[W }1 ( 1.
In order to provide a check of the optimization program,
a simple case, amenable to two solution techniques, was
chosen as the first example. A bare, homogeneous infinite
halfslab reactor was represented by a fivepoint, twoenergy
group model. The materials parameters for Example 1 are
given in Table 1. In this model, there are ten flux values
and, with each point designated a control point, there are
ten poison eigenvectors. Because the energy and space
variables are separable in diffusion theory models of bare,
homogeneous reactors, these vectors may be associated as
pairs, with identical shapes but different (spatially con
stant) fast/thermal flux ratios. The first two pairs of
modes are illustrated in Figure 2 and their adjoints are
shown in Figure 3. The ten eigenvalues are listed in Table. 2.
68
TABLE 1
Reactor Model for Example 1
Core Configuration: bare, homogeneous, infinite slab,
80 cm total thickness.
Mathematical Model: diffusion theory, 2 energy groups,
5 space points representing halfslab, mirror condi
tion at centerplane (point 1), zero flux condition at
extrapolated boundary (point 6).
Materials Parameters (in standard metric units):
Group (k) D(k) Ea(k) ER(k) vZf (k) X(k)
1 1.46 .0122 .0175 .0126 1
2 0.37 .1841 .2666 0
lattice spacing = 8 cm
Control Poison Cross Sections:
(1) = 0. cm per gram/c
o (2) = 0.1 cm per gram/cc
p 1
o (2) = 0.2 cm per gram/cc
p
S.AA ._..._, .
0 Fast Flux
SThermal Flux
0 1 2 3
0 10 20 30 40
X (cm)
2
/
I I I I J
0 10 20 30 40
X (cm)
0
0 0
I I I I
0 10 20 30 40
X (cm)
o
o
0 I I I 3 4
0 10 20 30 40
X (cm)
+ for Example 1
Poison Eigenvectors +I i' 2 and 2 for Example 1
Figure 2
A
0 0
A
o A
o Fast Flux
& Thermal Flux
S 10 20 30 40
0 I0 20 30 40
(
2
0
X (cm)
I 1 II
) 10 20 30 4C
0
0 /
Aa
A A
0 10 20 30 40
X (cm)
2A
I ~ I A
0 10 20
X (cm)
30 40
X (cm)
Figure 3
+ a+
Adjoint Vectors +lI' , 2 and 2 for Example 1
1 1, 2 2
TABLE 2
The Eigenvalues Associated with the Poison
Eigenvectors for Example 1
Eigenvalue
0.047070
0.086924
0.288781
0.467260
0.563996
Mode
Index
1
2
3
4
5
Eigenvalue
1.163727
1.216454
1.316718
1.440360
1.530345
Mode
Index
1+
2+
3+
4+
5+
The behavior discussed in Chapter V for the poison eigen
value series may easily be seen in this example. Because
the fundamental mode in the unpoisoned reactor is super
+ +
critical, 8 is positive. All the other + grow more nega
tive with increasing index. The B series are all negative,
and also become more negative with increasing index. It
is also apparent that, if enough points were added to the
model, the two series may overlap. That is, for some index
i great enough, there would be i < ~ This is one of the
phenomena which complicates any automatic process for finding
the whole set of eigenvalues with no human intervention.
The pointwise optimum poison concentrations for Example
1 were found by the modified GAUSS program, using the tech
niques outlined in Chapter VI. The results are displayed
in Table 3.
In order to provide a check on the functioning of the
optimization coding, this problem was also solved by the
following alternative method. A power distribution vector
was associated with each of these poison eigenvectors, and
+ +
is designated by P for the mode. These power vectors
have only spatial dependence, and thus only five elements
each. The elements of the power vectors, P (j), are related
+
to the elements of the flux vectors, (k,j), by the equa
tion
+ + +
Pi(j) = A{Zf(1) i(l,j)+Zf(2) i(2,j)}
+ +
= A{Zf (1)Si f (2) } (2,j) (101)
f i f I
TABLE 3
Optimization Results for Example 1
Results of General Poison Distribution Search Method:
Fast
Flux
2.75802E 01
2.75777E 01
2.75518E 01
2.73029E 01
2.49391E 01
Thermal
Flux
2.44744E 00
2.44756E 00
2.44878E 00
2.46055E 00
2.57227E 00
Power
Density
1.0
1.0
1.0
1.0
1.0
Results of Power Vector Expansion Coefficient
Real Solution:
Fast
Flux
2.75802E 01
2.75778E 01
2.75518E 01
2.73029E 01
2.49391E 01
Thermal
Flux
2.47444E 00
2.47555E 00
2.48788E 00
2.40555E 00
2.52266E 00
Power
Density
1.0
1.0
1.0
1.0
1.0
Poison
Concentration
6.55390E02
6.54151E02
6.41074E02
5.16011E02
1.02413E01
Search Method:
Poison
Concentration
6.55372E02
6.54153E02
6.41076E02
5.16005E02
1.02412E02
Expansion
Coefficient
2 .7200 1E 1
8.81539E 01
4.88696E 01
2.71797E 01
8.85565E 01
Mode
Index
1
2
3
4
5
Expansion
Coefficient
3.82306E02
1.49576E03
1.66447E01
2.21372E01
9.48818E02
Second Mathematical Solution (Not Physically Possible):
Fast Thermal
Flux Flux
1.19011E 01 4.31338E 00
1.12463E 01 4.28244E 00
1.23898E 01 4.33649E 00
9.86252E 00 4.21704E 00
1.54869E 01 4.48286E 00
Expansion
Coefficient
4.48286E 02
4.21704E 02
4.33649E 02
4.28244E 02
4.31338E 02
Mode
Index
1
2
3
4
5
Power
Density
1.0
1.0
1.0
1.0
1.0
Poison
Concentration
1.16233E 00
1.14970E 00
1.17165E 00
1.12248E 00
1.25350E 00
Expansion
Coefficient
1.54869E 03
9.86252E 02
1.22898E 03
1.13463E 03
1.19011E 03
Point
No.
1
2
3
4
5
Point
No.
1
2
3
4
5
Mode
Index
1+
2+
3+
4+
5+
Point
No.
1
2
3
4
5
Mode
Index
1+
2+
3+
4+
5+
A similar series of power vectors, P., is defined for the
+
Series of modes. Since 4 and 4. have the same spatial
shape, and differ only by their fast/thermal flux ratios,
+ +
S. and S., their associated power vectors, P. and P., dif
fer only in magnitude. That is, P = C.P.
i 1 1
The "desired power shape" is everywhere constant, and,
when normalized to a unit power density, is given for this
example by the vector P* = (1,1,1,1,1). If the flux in the
reactor, p*, that results in this power shape is expanded
as a sum of the poison eigenvectors,
5
= ia+ +a ~ i (102)
i=l
then the desired power vector, P*, is given by
5 5
+ + + +
P* = a.P.+aiP = (a +a./c )P. (103)
1 1 1 1 1 1 i 1*
i=l i=l
+
As may be seen from equation (103) the set of five P. vectors
is alone sufficient to expand the desired power vector.
Such an expansion has the form
5
P*= b.P. (104)
i=l
The vector equation (104) may be solved by simple matrix in
version to obtain the bi. By equating coefficients for P.
in equations (103) and (104), the following relations are
derived.
+4
(ai+ai/ci) = b. (105)
If the five new variables, d., are introduced such that
+
a. = bi/di' (106)
1 1 1
then, from equation (105),
ai = cibi(dil)/d.. (107)
Thus, the choice of the five variables, d., allows the com
putation of a set of ten flux expansion coefficients, a.,
which allow the construction of a flux satisfying the flat
power condition. With these coefficients and the poison
eigenvectors, the poison concentrations necessary for cri
ticality may be found by solving equation (77). After making
the substitutions given by equations (78) and (79), the
solution for each diagonal element of the [N] matrix becomes
10 10
N(k,j) = aiBi (k,j)/ ai i(k,j). (108)
i=l i=l
However, the concentrations for both groups must be the
same at any one point, so there are five constraint equations.
N(l,j) = N(2,j) (109)
Unfortunately, equations (108) are nonlinear in the di when
these variables are substituted for the a.. The problem
was solved by an optimization search on the variables di,
with the optimum point defined by the satisfaction of the
constraint equations (109).
This search was carried out by an unmodified version
of the GAUSS program, and the results are given in Table 3.
As may be seen, both methods result in the same solution,
indicating that the general program does indeed find the
optimum poison distribution. The poison eigenvectors for
this example satisfy their orthogonality relations to the
T
degree that the products m [p]~i, for m ? i, have values
6
on the order of 10 or less. With this fact in mind, it
may be seen from the expansion coefficients for the optimum
flux, given in Table 3, that all of the modes give signifi
cant contributions to the expansion. The modes with negative
fast/thermal flux ratios give only small contributions in
this case. More will be said about the significance of
this series of modes later in the chapter, when partial
expansions are considered.
Before leaving Example 1, it is useful to point out
that another solution exists to equations (108) and (109).
This mathematical solution has negative fast flux components,
and may therefore be excluded from consideration on physi
cal grounds. However, it serves to illustrate the fact that
the nonlinearity of the mathematics does allow for more than
one solution, as discussed in Chapter III. Because a two
group formulation was used in this example, two mathematical
solutions exist. This nonreal solution is also included
in Table 3.
Example 2 was run in order to provide an independent
check of the CODE2 eigenvector generator and the GAUSSUFUNT
eigenvector expanded poison optimization routine. This
example employs a twogroup, tenpoint model for a bare,
homogeneous, infinte cylindrical reactor. The reactor
configuration from Example 2 is given in Table 4.
The fundamental poison eigenvalue and vector were
checked against a calculation performed by the CORA multi
group diffusion theory program in which a uniform poison
concentration, B1, had been added to the materials. The
eigenvector, the CORA generated fluxes, and the analytic
solution are all plotted in Figure 4. As may be seen, the
comparison is excellent. Because the CODE2 and CORA values
occur at the same radial locations, and the two sets of
values are identical to five significant digits, only one
set of points is shown. In addition, the reactivity of the
flux mode calculated by CORA was within 105 of unity, the
value assumed in the generation of the poison eigenvectors.
For the unpoisoned core, CORA and the analytic solution give
a reactivity of 1.0603.
Because the CORA program is not capable of handling
pointwise variations in the poison cross section, the opti
mization was performed for three regions of poison. The
poison concentrations at the control points were constrained
to vary as three groups in the manner described in Table 4.
Because CORA requires that materials interfaces occur on
mesh points, it was necessary to locate the boundaries of
the constant poison regions at the mesh points used for the
eigenvector calculations, in order to preserve identical
TABLE 4
Reactor Model for Example 2
Core Configuration: bare, homogeneous, infinite cylinder,
150 cm radius
Mathematical Model: diffusion theory, 2 energy groups, 10
space points along radius, zero current condition at
centerline (point 1), zero flux condition at extra
polated boundary (point 11)
Materials Parameters (in standard metric units):
Group (k) D(k) Ea(k) ER(k) vf (k) X(k)
1 1.461860 0.0103369 0.0169434 0.005872 1
2 0.366341 0.0979940 0.135760 0
lattice spacing = 15 cm
Control Poison Cross Sections:
a (1) = 2.86715 cm1 per gm of Bl/cc
P 10
a (2) = 130.693 cm per gm of B /cc
p
To achieve regionwise poison variation:
points 1 through 3 are constrained to the same value,
points 5 and 6 are constrained to the same value,
points 8 through 10 are constrained to the same value,
point 4 is constrained to the volume weighted average
of the values in the inner and middle regions,
point 7 is constrained to the volume weighted average
of the values in the middle and outer regions.
0 CODE2 and CORA Values
 Analytic Solution
+
0
0F
Figure 4
Radius (cm)
Comparison of the Poison Eigenvector,
by CODE2, by CORA and Analytically
+ 'o
, for Example 2 as Calculated
mathematical models between the two calculations. In
order to do this, the concentrations at the interface points
were taken as the average of the concentrations in the re
gions to either side, weighted by the fraction of the reac
tor volume that is in each region and associated with that
mesh point. The results of the optimization are given in
Table 5. A comparison calculation was then made for this
optimum configuration using the CORA program. The results
show excellent agreement. The flux values are identical to
at least three, usually four significant digits. These
values are plotted in Figure 5. The CORA calculated reac
tivity for this materials configuration is 1.0000114, which
is also in excellent agreement with the value of unity
assumed in using the eigenvector expansion technique.
It should be pointed out that regionwise poison varia
tion may lead to more than one "optimum realizable power
distribution". This is due to the fact that the error
value, given by equation (98), is then a function of more
points than there are available control parameters. The
problem is thus underspecified, and the minimum value of
the error will not, in general, be zero. If the poison
concentrations available for manipulation by GAUSS are
considered as the ordinates of a "concentration hyperspace",
then the error function may be considered a hypersurface
within that space. As shown already by the two solutions
to Example 1, there are more than one minimum on this error
hypersurface for a completely specified problem. (All but
TABLE 5
ThreeRegion Poison Optimization Results for Example 2
Poison
Region
Point
No.
Poison
Concentration
4.66710E05
4.66710E05
4.66710E05
4.66710E05
9.68959E05
9.68959E05
9.68959E05
9.68959E05
7.47632E06
7.47632E06
7.47632E06
7.47632E06
Power
Density
1.20655
1.18203
1.10942
0.97136
0.97136
0.90423
0.97290
1.19125
1.19125
1.32066
1.06908
0.58144
Volume
Weight
2.77008E03
2.21607E02
4.43213E02
3.04708E02
3.60111E02
8.86427E02
5.30931E02
5.77099E02
1.32964E01
1.55126E01
1.77285E01
1.99446E01
1
2
3
4
4+
5
6
7
7+
8
9
10
01
0
Fast Flux
Poison I Poison IPoison
Region I I Region II Region III
Thermal Flux 
10 20 30 40 50 60 70 80 90 100 110
Radius (cm)
120
Figure 5
Flux Distribution for the Optimum ThreeRegion Poison Distribution in
Example 2
one of these may be eliminated by physical considerations,
however.) If the problem becomes underspecified, the hyper
space has fewer dimensions and the error surface may develop
extra local minima. This is a phenomenon that is caused
by the physically underspecified nature of the problem,
rather than by the techniques used for solving that problem.
Underspecified problems may be created by ganging control
points for regionwise variation or by the specification of
some fueled points that are not control points. In any
event, no matter how these problems are treated, it must be
remembered that an optimization search program may find a
local minimum that is not the lowest minimum. Indeed,
there may be more than one lowest minimum, although this is
unlikely.
Now that it has been established that the eigenvector
expansion technique may be expected to give a reasonable
representation of the unexpanded problem, the possibility
of using expansions with less than the complete number of
vectors may be investigated. Because Examples 1 and 2
both involve cases where the energy and space variables
are separable, partial expansions of these problems are
special cases. If the vectors are eliminated in the pairs,
+
Oi and 4., some of the ability to represent spatial shape
has been lost. However, any shape that may be generated by
the remaining vectors may have its (spatially dependent)
spectrum multiplied by any arbitrary magnitude. The spectrum
of each shape component is still arbitrarily adjustable
between the two vectors with that shape. For a reactor
with materials interfaces, the energy and space variables
are no longer separable, and the vectors i and i will
no longer share the same shape. In such a case, elimination
of some of the vectors may also impare the representation
of the energy spectrum of the flux shapes which may still
be generated.
In order to investigate partial expansions for general
multigroup problems, Example 3 was contrived. The model is
a twogroup, nineteenpoint representation of an infinite
cylindrical reactor with two fueled regions and a water
reflector. The details of this reactor model are given in
Table 6. The fifteen control points result in thirty poi
son eigenvectors. The 21 and I1 modes are displayed in
Figure 6 to illustrate their dissimilarity in shape. In
addition, it may be seen that the fast/thermal flux ratio
is now a function of position in the reactor for each mode.
The adjoint vectors for these modes are shown in Figure 7.
Partial expansions can be expected to be of the most
utility when the poison is to be varied by region. In such
cases, if ever, a few modes might be expected to adequately
represent the somewhat rippled flux shape that results. How
ever, because regionwise variation of poison involves step
changes in poison concentration, all of the eigenvectors
may still be needed to represent such situations with maximum
accuracy. That is, although contributions from the higher
TABLE 6
Reactor Model for Example 3
Core Configuration: infinite cylinder, two homogeneous
fueled regions from radii 0 to 50 cm and 50 to 140 cm,
pure water reflector from radius 140 to 190 cm
Mathematical Model: diffusion theory, 2 energy groups,
19 space points, 15 control points, zero current
condition at centerline (point 1), zero flux condition
at extrapolated boundary (point 20), continuous neutron
current conditions at internal interfaces (points
6 and 15)
Materials Parameters (in standard metric units):
Region I:
Group (k)
1
2
Region II:
Group (k)
1
2
Reflector:
Group (k)
1
2
Westinghouse PWR fuel assemblies, 3% enrichment,
18,000 MWD/T
D(k) a (k) ZR(k) C f (k) X(k)
1.466676 .01083982 .01648968 .00542511 1
0.364366 .09871243 .1305982 0
Westinghouse PWR fuel assemblies, 3% enrichment,
10,000 MWD/T
D(k) E (k) E (k) v, (k) X(k)
1.461862
0.366341
.01033693
.09799408
pure water
D(k) E (k)
1.136999 .0005892
0.149380 .0019240
.01694346
ER(k)
.04830
.00587230
.1357601
Zf (k)
0.0
0.0
1
0
X(k)
lattice spacing = 10 cm for all regions
Control Poison Cross Sections:
Op(l) = 2.86715 cm1 per gm of B /cc
a (2) = 130.693 cm1 per gm of Bl0/cc
Materials interfaces occur at points 6 and 15.
Control points are those numbered 1 through 15.
To achieve regionwise poison concentrations:
points 1 through 5 are constrained to the same value,
points 7 through 14 are constrained to the same value,
point 6 is constrained to the volume weighted average
of the values in the two regions,
point 15 is weighted by the fraction of its associated
volume that is within Region II.
86
i 'I
Region I Region II
Reflector
0 50 140 190
I i
1/ I
I
OO G
Radius (cm)
I I/
Region I I Region II I Reflector
I A Thermal Flux
/
0Ahe alFl I
I I I
0 50 140 190
Radius (cm)
Figure 6 Poison Eigenvectors +1 and $i for Example 3
