
Citation 
 Permanent Link:
 https://ufdc.ufl.edu/UF00097627/00001
Material Information
 Title:
 A multigroup eigenfunction expansion technique for selecting poison distributions for optimum reactor power distribution
 Creator:
 Long, Steven Michael, 1946
 Publication Date:
 1972
 Copyright Date:
 1972
 Language:
 English
 Physical Description:
 xiv, 151 leaves. : ; 28 cm.
Subjects
 Subjects / Keywords:
 Adjoints ( jstor )
Diffusion theory ( jstor ) Eigenfunctions ( jstor ) Eigenvalues ( jstor ) Eigenvectors ( jstor ) Error rates ( jstor ) Mathematical variables ( jstor ) Mathematical vectors ( jstor ) Matrices ( jstor ) Poisons ( jstor ) Dissertations, Academic  Nuclear Engineering Sciences  UF ( lcsh ) Mathematical optimization ( lcsh ) Nuclear Engineering Sciences thesis Ph. D ( lcsh ) Nuclear reactors  Problems, exercises, etc ( lcsh )
 Genre:
 bibliography ( marcgt )
nonfiction ( marcgt )
Notes
 Thesis:
 ThesisUniversity of Florida.
 Bibliography:
 Bibliography: leaves 148150.
 Additional Physical Form:
 Also available on World Wide Web
 General Note:
 Typescript.
 General Note:
 Vita.
Record Information
 Source Institution:
 University of Florida
 Holding Location:
 University of Florida
 Rights Management:
 Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for nonprofit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
 Resource Identifier:
 029445173 ( AlephBibNum )
AEG6379 ( NOTIS ) 014267171 ( OCLC )

Downloads 
This item has the following downloads:

Full Text 
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

Full Text 
PAGE 1
A Multigroup Eigenfunction Expansion Technique For Selecting Poison Distributions For Optimum Reactor Power Distribution By STEVEN MICHAEL LONG A DISSERTATION PRESENTED TO THE GR/iDUATE COUNCIL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 197 2
PAGE 2
This dissertation is dedicated to my lovely v/ife, Diane, who ungrudgingly shared the effort of completing the beast during our first months of marriage.
PAGE 3
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. CD. 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 appreciation 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.vS. 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. Ill
PAGE 4
TABLE OF CONTENTS Page ACKNOWLEDGEMENTS iii LIST OF TABLES vi LIST OF FIGURES vii NOMENCLATURE viii ABSTRACT xii CHAPTERS: I. INTRODUCTION 1 II. THE ONEGROUP EIGENFUNCTION EXPANSION PROBLEM 6 III. THE POISON DISTRIBUTION PROBLEM IN MULTIGROUP FORMULISM 12 IV. SELECTION OF A MULTIGROUP EIGENFUNCTION SET 18 V. GENERATION OF THE POISON EIGENVECTORS ... .35 VI. THE MULTIGROUP OPTIMIZATION PROBLEM 47 VII. NUrlERICAL CALCULATIONS AND THE PROSPECT FOR PARTIAL EXPANSIONS 65 VIII. CONCLUSIONS AND SUGGESTIONS FOR FURTHER STUDY 104 APPENDICES: I. A DESCRIPTION OF THE TECHNIQUES EMPLOYED BY THE EIGENVECTOR GENERATING PROGRAM "C0DE2" 113 IV
PAGE 5
Page II. LISTING OF THE DEMONSTRATION POISON EIGENVECTOR GENERATION PROGRAM "C0DE2" 119 III. A DESCRIPTION OF THE "GAUSS" OPTIMIZATION CODE 130 IV. LISTING OF THE USER SUPPLIED SUBROUTINES FOR THE USE OF THE "GAUSS" PROGRAM TO FIND THE OPTIMUM CRITICAL POISON DISTRIBUTION 135 LIST OF REFERENCES 148 BIOGRAPHICAL SKETCH 151
PAGE 6
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 ()t and 4)7 for Example 4 . 101 VI
PAGE 7
LIST OF FIGURES Figure Page 1 The [B] Matrix and Flux Vector for the "Fission Eigenfunctions" 30 2 Poison Eigenvectors (f), , ^, 2 ^^'^ 'i'j for Example 1 6 9 3 Adjoint Vectors ^,, tp7, ^j and ip, for Example 1 i x z z 7q 4 Comparison of the Poison Eigenvector, ({;, , for Example 2 As Calculated by C0DE2, by CORA and Analytically 79 5 Flux Distribution for the Optimum ThreeRegion Poison Distribution in Example 2 82 6 Poison Eigenvectors (^t and i^l for Example 3 86 7 Adjoint Vectors ii, and ip, for Example 3 87 8 Power Distribution for the Optimum TwoRegion Poison Distribution In Example 3 92 9 Schematic Representation of Reactor Model for Example 4 100 vii
PAGE 8
NOMENCLATURE A energy per fission (in appropriate units) 3.^ expansion coefficient associated with function i a vector representing a set of expansion coefficients [A] matrix on lefthand side of various eigenvalue formulations of the criticality equation 2 B axial geometric buckling [B] matrix on the righthand side of the general eigenvalue formulation of the criticality equation C(k,i,j) 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 ]^ _ D (r) diffusion coefficient function in energy group k D(k,j) diffusion coefficient at point j in group k E weighted scalar squared error optimization parameter E vector of error components F(m,i,[N]) the scalar product of the poison eigenvector i and the adjoint vector m, weighted by the operator [N] [a ] [F] matrix with elements F(m,i,tN]) G number of control points H number of fueled points I number of eigenvectors i preferred index for designating a particular eigenvector Vlll
PAGE 9
[I] the identity matrix J number of spatial points J (r) zeroth order Bessel function of the first kind [J] the Jacobian matrix 3 M
PAGE 10
S (r) spectrum function relating the flux in group k to the flux in group K S. fast/thermal flux ratio for eigenf unction i 1 1 JL (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 velocity terms of the Laplace transformed kinetic equation [W] diagonal matrix with elements V(j) X distance from centerplane in slab geometry a. eigenvalue associated with the eigenf unction ^ i to a general eigenvalue problem 3. eigenvalue associated with poison eigenfunction i [3] diagonal matrix with elements g. y. eigenvalue associated with the "natural" ^ eigenfunction i 6. Wielandt transformed eigenvalue associated with the mode i 6(i,j) Kronecker delta e scalar coefficient for the adjustment of the magnitude of a poison distribution e. poison distribution magnitude eigenvalue associated with mode i n. zero number i of the function J_(r) X o X. eigenvalue associated with the fission eigen^ function i y. eigenvalue associated with the thermal poison eigenfunction i number of neutrons per fission
PAGE 11
Â„k poison cross section for a unit concentration P , ^'^(j) equals 0^ if point j is a control point, P zero othirwise k [^ ] diagonal matrix with elements o^Cj) j;k(7) absorption cross section function for group k j; (k,j) absorption cross section for group k at ^ point j j;k() fission cross section function for group k ^ (k^j) fission cross section for group k at point j 2k() controllable poison cross section function P for group k 2k() removal cross section function for group k R 2 (k^j) removal cross section for group k at point j Z (n>k,j) cross section at point j for scattering from s group n to group k <^^{t) flux function for group k (j,*(F) ideal flux function ^(k^j) flux at point j in group k A flux vector with elements <)(k,j) ^ . flux eigenvector i X^ or X(k) fission spectrum fraction in group k m u,. adjoint flux at point j in group k adjoint flux vector m eigenvalue associated with kinetic mode i XX
PAGE 12
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 TECPINIQUE 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 distributions as they are selected for investigation by a numerical optimization routine. This method is based in multigroup diffusion theory, and employs the finite difference approximations to the spatial dimensions. First, the eigenfunction expansion problem is reviewed in onegroup theory. Then, the multigroup problem is outlined 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 Xll
PAGE 13
problem is discussed for the case where these unknowns are treated directly as variables. The orthogonality and completeness properties of various 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 twogroup, 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 criticality 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 eigenvalues 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
PAGE 14
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 eigenvectors for expansion. The results of partial expansions in both underspecif ied and completely specified optimizations are discussed. Considerable accuracy of partial expansions is shown in the limited class of examples considered (thermal reactors) . Finally, the author presents his thoughts on the potential 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
PAGE 15
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 multigroup diffusion theory. The spatial dimensions are represented by the finite difference approximation, and the functions 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 representation of the physical problem and lend computational power as well. The problem of determining the control poison distribution for the optimization of some reactor performance parameter has been treated by many authors, using various methods and mathematical formulisms. Harker used a pseudo twogroup diffusion theory that made use of a "fast fission factor", 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
PAGE 16
which would result in a physically realizable flux shape. (2) (3) Haling and Crowther used a more sophisticated twogroup 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 "Pov/er 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 tv/oregion, bare core with control rods ganged in each region. Terney and French used dynamic programming to minimize power peaking over core life, based on a direct flux synthesis method (using two fluxes) , for approximating a tv;ogroup, twodimensional reactor with two gangs of control rods. Wade and Terney applied optimal control theory to a generalized set of design objectives, using a onegroup spatially nodalized reactor model. Jack(7) son used an eigenf unction expansion technique, based on onegroup diffusion theory, to determine the optimum physically realizable power and flux distributions and the poison distribution 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, v/hich reviev;s the onegroup eigenf unction expansion problem.
PAGE 17
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 v;ork and in some kinetics work, little use has been made of the higher modes to these problems. Even for perturbation calculations, only the fundamental eigenf unction (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 solution of group diffusion equations which involved expansion of the problem in the eigenfunctions of the Helmholtz equation for the appropriate geometry. An article published in { 9 ) 1962 by Garabedian and Thomas^ extended the approach to two dimensions. While these methods have the advantages of using wellknown sets of tabulated functions and analytically handling the leakage terms within homogeneous regions, they possess a distinct disadvantage. Because these functions 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 Shibata ^''^' combined this approach with variational calculus in a study to optimize total pov;er and
PAGE 18
conversion ratio. Their independent variables were regionwise 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 in 1961. He illustrated techniques 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 particular problem the "natural modes" for that problem. In (12) 1966, Kaplan and Yasinsky used an analysis of the eigenvalues of the natural modes to predict the behavior of xenon (12) oscillations. Henry, ' in 1963, extended modal approximations 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. (1A) In 19 63, the MULE computer code ^ became available for computing three particular types of eigenvector sets to twogroup diffusion theory reactor models. These eigenfunctions sets included: 1) the "X modes" or"fission eigenfunctions"
PAGE 19
where A is the inverse reactivity for the mode, 2) the "y modes", where y is the concentration multiplier for a given normalized distribution of thermal poison, and 3) the "w modes" or eigenf unctions to the Laplace transformed, timespace separable kinetics problem. As will be discussed in Chapter IV, with the exception of the u modes, these modes are degenerate in the energy dimension, and are limited, for expansion purposes, to expansions of the spatial dimensions, 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 displayed in the determination of the poison distributions and concentrations necessary for criticality in the optimum realizable fluxpower shapes for several example reactors.
PAGE 20
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 treatment, 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 distribution 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 complicated 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)V4)(r)+vZ^(r)(J)(r)Z^(r)(j)(r)Zp(r)())(r) =0 (1) plus the boundary conditions on <}) (r) and/or its first derivative. Standard notation is used except that Z (r) is the macroscopic cross section of the controllable poison.
PAGE 21
The powerflux relationship is P(F) = AE^(F)
PAGE 22
VÂ«D(r)V<^^(r)+vZf (r)(})^(r)E^(r)4)^(r) = Yi<^i(r). (5) These eigenfunctions, the *. After invoking the eigenvalue relationships of equation (5) , this becomes
PAGE 23
(r) Zp(r) = iji . (9) i=l As pointed out by Jackson, this approach still has a failing. Because the optimum pov/er shape, P* (r) , is defined only for the reactor core, where E^(r)>0, but these eigenf unctions 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 Z^(r) = and the necessity of integration 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 param.eter, the adjustment of v;hich might allow the eigenf unction approximation to the ideal incore flux to be close in that core region in which the nearness of the approximation is considered (7 Â• 24) to be most important." Neither suggestion seems particulary 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 necessary for the property of finality, it also introduces an added complexity which seems to destroy much of the power of the method.
PAGE 24
10 Recognizing a better alternative himself, Jackson recommended the exploration of other sets of eigenfunctions which are complete and orthogonal over just the fueled portions of the reactor. He suggested a set he called the "fission eigenfunctions", given by the equation VD(?)Vc^i(F)E^(F)(})^(F) = Aj^vE^(F)i(r: E (F) = ^ . (13) 2l^^*^^^^ i=l
PAGE 25
11 This method at last seems to handle the onegroup problem 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 control rods, but these bring with them the requirement for optimization searches v/hich 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.
PAGE 26
CHAPTER III THE POISON DISTRIBUTION PROBLEM IN MULTIGROUP FORMULISM In contrast to some of the previous authors, 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 v/ill be generally considered to have a cross section in each group. Because of this general treatment, 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 equations has the general form K V . D^ (F) V(})^ (F) +z^ (F) c^^ (F) x^ 2 ^^f ^^^ '^'^ ^^^ n=l K +E^(F)<{)'^(F)2 Sg(n>k,F)(^"(F)+Ep(F) (r) , are required to satisfy homogeneous boundary conditions and to be continuous, along with their k Â— k Â— net currents, D (r)V({) (r) , at internal material interfaces. 12
PAGE 27
13 The multigroup equations are not self adjoint . The adjoint fluxes to the critical system satisfy a set of K equations of the form K vD^(?)yi^^(?)+E^(F)4'^(F)vZÂ£(F) 2^'^'^'^^^^ n=l K +E^(F)i;;^(F)yz^(k>n,F)^''(F)+zJ^(F)i];^(F) = (15) 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 pov/erflux relationship is given by K p(F) = A 5'z^(F)({)^(F) . (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 pov;er shape does not allow the solution for the multigroup fluxes. This occurs because there is no a priori 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
PAGE 28
14 the other fluxes. These relationships have the form 4)^(F) = S^(r)(})^(F), for k = 1 to (K1) . (17) For convenience in notation, the identity relation for the 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 k Â— for the optimum k,r)s''(r)/s'^(r) n=l (19)
PAGE 29
15 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 relationship 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 mixture, 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 Â— absorption cross section is known, and the 2_(i^) of equations (19) must be constrained to E^(F) = N(F)o^. (20) Thus, the real variables of the problem must be considered to be the S (r) unknown spectrum functions, of which there are (K1) , and the N(F) poison distribution function. This gives a total of K unknown functions.
PAGE 30
16 The K relationships among these unknowns are derived by making the substitutions from equations (20) into equations (19) . N(r)a'^ = VD^(r)V S^(r)P*(r) as^F)z? (r) n=l K A2s'^(r)E5(r) n=l S^(r)P*(r) K +x''2 s"" {r)vl'^ (r ) /S^ (r ) Z^ (r ) Z^ (r ) 'R n=l K n=l (21) 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 shov/n there will be K solutions, all but one of which will involve some group fluxes with negative magnitudes. Although the physical realities of the problem allow the selection of the proper solution (the one with all positive group fluxes) , these constraints may only be expressed mathematically as inequalities. S (r) > for all k. (22) Unfortunately, such inequalities are not suitable for
PAGE 31
17 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) v;ill not be valid for arbitrary P(r) and ^^(r) functions, it is best to adopt a viev/point 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 solution to the coupled, nonlinear requirements of determining the optimum realizable power distribution and the associated distribution of a previously selected control material. However, it will be shown that proper choice of an eigenfunction expansion technique can lend considerable power to the search for a solution.
PAGE 32
CHAPTER IV SELECTION OF A MULTIGROUP EIGENFUNCTION SET In order to find a useful set of eigenf unctions , it is first necessary to find some linear portion of the problem 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 ^ ''' ^ 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, cj) (r) . It is nonk Â— k Â— linear overall only due to the control terms, I {r)(p {r) , ir 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 eigenfunctions, some are m.ore 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 eigenfunctions" of Jackson, ' on the other hand, contain a 18
PAGE 33
19 degeneracy in the energy variable, and are suitable only for expansions of the spatial dimension. Other eigenf unctions will be introduced which overcome both of these objections. Before this is done, however, the author will introduce the multigroup, finite difference formulism that must generally be employed to calculate whatever eigenf unctions are to be used. This will allow the comparison of the various eigenfunction sets within the actual context 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 continuousfunctions 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 v/ritten together as a vector. Applying this to the group fluxes of equations (14) , the flux in group k would be represented as 4)^(F) = {4;^(?=F^);4)'^(F=F2);...;
PAGE 34
20 Equations (14) may be rewritten for each spatial position, 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. For example, in onedimensional, slab geometry, the gradient term becomes V.D'^(F)V.(f)^(F) = ^{D^(x)^''(x)} (24) which, for constant diffusion coefficient and evenly spaced x(j), would be approximated as _pk d ^ (x j_ ~ D(k) [(+. (k,j+l)24) (k,j)+(J) (k,jl)]/Ax^. (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 v/ill consist of J*K equations of the general form C (k, j , j +1 ) ({) (k , j +1 ) +C (k , j , j ) (j) (k, j ) +C(k,j,jl),^(k,jl)+E^ (k,j)4)(k,j) K 2 l^{nyk,j)
PAGE 35
21 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 i:p(k,j) = N(j)ap(k), (27) where the a (k) is a set of knovm values. Thus, the finite difference, raultigroup approximation to the criticality condition, as represented by equations (26), has J'(Kl) unknowns: the JÂ«K fluxes, ({)(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 v/ith 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 freedom 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 elements. ' Therefore, in order to expand the multigroup flux vector with a set of eigenvectors, it is necessary 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
PAGE 36
22 vector space may be spanned by only H K 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 arbitrary 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 eigenvalue problems, the resulting J'K equations will be linear in the J*K unknown fluxes. The sets of solutions to these (18 19 ) problems may be found by well established methods, ' v/hichwillbe discussed in the next chapter. These problems may all be v/ritten in matrix notation in the form [A]({.^ = a^[B]({.^, (28) where a. is the generalized eigenvalue. The actual form of 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, ^., will have (J'K) elements . Before considering particular eigenvalue problems, some observations about the general formulation are in order.
PAGE 37
23 The adjoint to equation (28) is given by [A]"^!^^ = a^lBl'^.j;^. (29) Following Hansen and Clark/ ' the orthogonality relationship may be derived as follows. Assuming that [A] has 'an inverse, then the "solution" to equation (28) may be written as (f^/a^ = [A]"'[B]()^; and the "solution" for the adjoint problem is ,i;^/a^ = {[A]'^}~^[B]'^ii;^ = { [A] "^ l"^ [B] ^ij;^ . (30) (31) The adjoint to equation (30) is e^/a^ = {(A]"^[B]}'^6j_ = [B]'^{[A]"^}^6j_. (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, ^ ' 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 desirable orthogonality conditions. The e's satisfy the usual, unweighted condition Q.^^= 6(i,j), and are thus orthogonal only over the whole JÂ«K vector space. A more useful orthogonality condition may be derived for the ijj. adjoint set. If the operator [B] is applied to both sides of equation (30) , the result is
PAGE 38
24 [B](J)^/a^ = IB] [A]"'[B]<})^. (33) Now, if [B](^^ is considered as a vector, its adjoint, denoted by w. , must satisfy the equation w^/a^ = {[A]"^}'^[B]'^w^. (34) By comparison of equations (31) and (34) , it can be seen that V7^ is identical to ip^. By taking the inner product of the transpose of equation (34) with the vector [B] (p . , and subtracting the inner product of w. with equation (33) , the result is 
PAGE 39
25 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 a^ 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 restriction is placed on [B] . The rank of a matrix is equal to the number of rows (or columns) in the matrix that are linearly 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, therefore, have zero determinates. The coefficients of all terms a", v/here n>m, will therefore be zero. Equation (38) will thus be an ra 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 "oj 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, however, such roots do not seem to occur for formulations where the eigenvalue may be considered as the multiplier of
PAGE 40
26 some real physical parameter. Double or higher order repeated 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".^ ' 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+1) .J,^ (k, j+1) +C (k, j , j ) ({)^ (k, j ) +C(k,j,jl)(^^(k, jl)+E^(k,j)cf)^(k,j) K X(k)2 vZ^(n,j)(t)^(n,j)+Zj^(k,j)({)^(k,j) n=l K "2^s^^'^'^'^^*i^'^'^^ = Yi4>i(k,j), (39) n=l where y is the eigenvalue for these modes. In matrix notation, this is [A]<^^ = y^[l]
PAGE 41
27 an inverse unless the reactor is just critical. (That is, det[A] = is just the criticality condition.) The identity 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 orthogonality conditions 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 JK vectors may contain more information 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. Hov/ever , these eigenvectors, with J'K elements, satisfy the orthogonality conditions with their ad joints only over the entire J'K 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
PAGE 42
28 C(k, j,j+l)c})^(k,j+l)+c(k,j,j)(^j_(k,j) +C(k,j, j1) ii)^(k, jl) + E^(k, j) 4,^(k,j) K "2 ^s^^^^'^^*i^^'^^'^^R^^'^^*i^^'^^ n=l K = A^X(k)^ vEÂ£(n,j)(})(n,j) , (42) n=l where A. , the fission mode eigenvalue, is equal to the inverse effective reactivity of the mode i. The orthogonality condition leads to an expression for the expansion coefficients for a given flux vector, 4), given by a^ = i(;T[B]({), (43) that weights any fluxes in unfueled portions of the reactor 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
PAGE 43
29 hand sides of equations (42) are being "driven" by a neutron 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 eigenf unctions that is complete over the desired space can be generated by equations of the form C(k, j,j+l)<})^(k,j+l)+C(k,j,j)(j)^(k,j) .+C(k, j, jl)(j,^(k, jl)+z^(k,j)<^^(k, j) K X(k) 2vZg(n,j)(j)^(n, j)+5:j^(k, j)(()^(k, j) n=l K "2 ^s^^"^^'^^*i^"'^^ = Biap(k,j)
PAGE 44
30 eeCN (N eeo o CM
PAGE 45
31 is to be used. Physically, this is equivalent to specifying that the concentration of poison is the same at each control point. The eigenvalue, 3^, 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 GK, 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]i = 3i[a ]4>^, (45) where [a ] is the [B] matrix for the poison eigenvectors. These eigenvectors still have J'K elements, and so represent 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 straightforward. Consider the finite difference criticality equation for the controlled reactor, [A]c{. = [S ]4). (46) The matrix [Z ] is diagonal and represents any critical
PAGE 46
32 poison distribution, generally with different concentrations at each control point. Like [o ], this critical control matrix, [ I ] , 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 relations. ^i " ^if^pl* (47) The flux at the control points, at least, is represented by * = Z^i^i(48) i=l Now, if (})* is different from (j) , it must differ only at noncontrol points, since the poison eigenvectors are complete over the control points. Let any difference be denoted as e = (j)(j)*. (49) It is easily seen that the constructed flux satisfies the criticality equation (46) by using the eigenvalue relations (45) . I I (50) X JL J. i=l i=l Since the matrices [a ] and [Z ] have zeros on the diagonal for all noncontrol points, the substitution of (^* in the criticality equation yields a relationship dependent only J. J.
PAGE 47
33 upon the fluxes at control points, v/here ^* and (p are identical . I I i=l 1=1 Thus, (+)* also satisfies the criticality relation. Now, by substituting (})*(}) for ({> in equation (46) , [A] (4)*e) = [Ep] ((j)*9) . (52) Since (})* already satisfies equation (46) , this may be reduced to [A] 6 = [Zp] 9. (53) However, 6 must be zero at all control points, and the operator [Z ] multiplies all noncontrol points by zero, so IT [Z ] Q = U. (54) Thus, if the difference vector, 6, exists, it must satisfy [A] e = . (55) So, for a nontrivial solution for 6 to exist, it is necessary that det[A] E 0. (56) If this condition is not satisfied, then 6 is the zero vector 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.
PAGE 48
34 a poison eigenvector will exist with a zero poison eigenvalue. [A]<^^ = 0[ap]({)^ E 0" (57) Thus, if it existed, G must be identical to 4> / ^rid 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 distributing poison at the designated control points. These poison eigenfunctions at last seem to be able to adequately represent the physical situation. For a reactor with G control points (e.g., rods, regions, etc.), there are only G degrees of freedom for control purposes. Another G'(Kl) variables are needed for expansion, however, because the spectrum is not known a priori . This set of eigenvectors provides the minimum number of necessary variables and insures that they are complete over only those spatial points that may actually be directly controlled.
PAGE 49
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 themselves 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 eigenf unctions . 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 eigenf unction. (14) An exception is the MULE code/ which is capable of finding any desired mode for three types of twogroup problems, the fission eigenf unctions, the poison eigenfunctions for a thermal poison only , and the time or kinetic eigenfunctions, with or without delayed neutrons. The fission modes satisfy the equation [L]<})^ = X^[S]*^, (58) where [L] contains all of the leakage and absorption loss 35
PAGE 50
36 terms and [S] contains the fission source terms. As discussed 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 eigenfunctions have a similar degeneracy. They satisfy the relation {[L]^^[S] }<^. = y. (Q]<})., (59) ^eff 111 V7here [Q] is a diagonal matrix with zeros on the diagonal for all the fast group equations and Z 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]^[S]}'4,. = oj.[V]4)., (60) eff 111 where [V] is a diagonal matrix with the inverse of the neutron velocity for the appropriate group on each diagonal location. The rank of [V] is thus 2J, and the eigenfunctions of this type will be complete over energy and space. This problem is very similar to the twogroup poison eigenfunctions 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 eigenf unctions 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
PAGE 51
37 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 v;ith the program while calculations are in progress. The desirability of such an arrangement will become apparent shortly. The actual calculation of an eigenvector is typically ( 18 ) accomplished by the pov;er iteration method. For the general eigenvalue problem of equation (28) , [A]<)). = a^[B](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 i=l If [A] has an inverse, then the eigenvalue problem may be written in the form of equation (30) . /a^ = [A]"^[B]()j_ (30) The effect of operating on vector V by the matrix [A] [B] can be found by substitution for each eigenvector from equation (30) .
PAGE 52
38 [A]~^[B] I Â£^ i i dm* a^c^i/a. (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/a^. If a^^ 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]~"^[B], the eigenvector component v/ith 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 dominant eigenvector will remain. As can be seen from equation (30) , when the guess vector has converged to the eigenvector 4)^, 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 a.(},./a. a . d) . 1^1 a .(J) ./a Â• a . t Â• a JL a . 1 (63) The greater the dominance ratio, the more rapidly vector j will die away compared to vector i.
PAGE 53
39 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.4>^ = [B]"^[A](},j^, (64) and operating repeatedly with the matrix [B] [A]. Still, the vectors with eigenvalues of interm.ediate modulus must be found by some other method. There are tV70 possible CIO) methods. The first is successive elimination 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 sum in equation (61) no longer contains these components, or by modifying the original matrix 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. (18 19) A better method, proposed by Wielandt, ' 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]ag[B]}*^ = 6^[B]^^, (65)
PAGE 54
40 where a is the value of the guess, and the eigenvalues of the transformed equation (64) are related to the original eigenvalues by ^i = '^i"^Â°'g(^^) The dominance ratio for the new equation is calculated using the 6's. Obviously, the smallest 6 will be the one associated with the a. 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 clor jr 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 VJielandt transformation itself may be used as an acceleration technique if the eigenvalue guess is updated using partially converged T . J, (21,22) values of . To aid in guessing close to the eigenvalues, it is helpful to consider the physical significance of the eigenvalue problem chosen. The poison eigenf unctions satisfy equation (45) . [A]^^ = Bi[ap](j)^ (45) The [A] matrix contains all of the terms for the uncontrolled reactor, and matrix [a ] contains the terms for absorption
PAGE 55
41 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 fundamental mode (the most reactive one with flux that is everywhere positive) will require some control poison insertion to reach the critical state. Thus, the associated g, 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 behave, it is instructive to examine a simple example. For a twogi*oup representation of a bare, homogeneous cylindrical reactor of infinite axial extent, the equations for the twogroup fluxes in the poison eigenvalue formulation are given by D(F)V^(})^(F,r)+{E^(F)+Ej^(F)vEf (F)}c},^(F,r) vZÂ£(T)(},^{T,r) = Biap(F)ff,^(F,r) (67a) for the fast group, and D(T)V^())^(T,r)+Z^(T)(^^(T,r)Z^(F)(},^(F,r) = Bj^a (T)(j)^(T,r) ' (67b) for the thermal group. Because the energy and space variables
PAGE 56
42 are separable for this example, the solutions to equations (67) may be written as fluxes of the form ())^(F,r) = S^J^{n^r/R) (68a) '})j, (T,r) = J^(u^r/R) , (68b) where R is the radius of the reactor, J is the zeroth order Bessel function, >i. is the ith zero of J , and Sis the 1 o 1 fast to thermal flux ratio, which is a constant over the entire reactor volume for each mode. Substitution of equations (68) into equation (67) allows the solution for the eigenvalue, 6, and the spectral constant, S., for each mode. These solutions are 6^ = {S^E^(F)D(T) [u^/R]2v^(T) }/ap(T) , (70) where a = aÂ„(F) EÂ„ (F) P ^ b = ap(T){D(F) [n^/R]2 + j^^(F)+5:^(F) } aÂ„(F){D(T) ['f./R]^ + Z^(T)} c = vZ^(T)a (T) . The quadratic nature of the solution is due to the twogroup formulism. In general, the degree of the equation for S. 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. .
PAGE 57
43 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 particular choice of the index, i, the fast and thermal fluxes will have the same spatial shape, differing only in magnitude by the factor S.. Because there are tv70 groups, there are two values of 6 and tv/o 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 (7 0) , the negative S^ will always be associated v/ith a negative B., while the positive S. may have either a positive or negative value of Bj_. In fact, as discussed previously, when the uncontrolled reactor is supercritical, the fundamental mode must have a positive value of 6As the index, i, and thus n, is increased in value, the leakage terms in equations (67), given by D(k)V^J^(Hj_r/R) = D(k) [n^r/R]^jQ(n^r/R) , (71) will become more positive. In order to satisfy equations (67), then, the values of 3 must become more negative. The same conclusion may also be reached by physical reasoning. The innermost region of any higher flux mode, where 0^rÂ£n.R/n., is physically equivalent to the fundamental mode of a smaller reactor. In order to make this smaller reactor critical, less control poison v;ill be needed because of the greater neutron leakage. Thus, again, as n^ increases, Bmust decrease.
PAGE 58
44 For this case, then, the 3's will occur in two sequences, denoted by 6and 6for the values associated with the positive and negative values of S., respectively. These sequences satisfy the relations S^+i < &l (72a) ^i+1 ^ ^i < 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 complicated 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, hov/ever . Experience has showi that, if the flux values associated with noncontrol points are ignored, the remaining values will shov/ the proper number of sign changes in both the energy and space variables. This is of great utility to the person v;ishing to generate a set of eigenvectors using the Wielandt transformation technique. Any guess at an eigenvalue will result in convergence to some eigenvalue and vector. By using the above criteria, it may easily be determined v;hich 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.
PAGE 59
45 for values of 3 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 6 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 3 sequences, variability of the control (completeness) regions betv/een specific problems, and the possibility of finding the same root more than once, greatly complicate the logic and bookkeeping necessary for a computer program that would be capable of finding the complete set of eigenf unctions 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 1300 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 3. 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.
PAGE 60
46 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 C0DE2 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 twogroup, onedimensional calculations in cylindrical geometry. Extentions to more groups, more dimensions and other geometries are straightf orv;ard, follov/ing 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.
PAGE 61
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 equations. 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 distribution gives the "error" in power resulting from the use of this poison distribution. By systematically varying the individual poison concentrations and recalculating the power 47
PAGE 62
48 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 location of, or direction to, the poison distribution yielding the optimum pov;er distribution. This results in a new poison distribution 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 pov/er distributions in the same reactor, one for each poison distribution to be investigated. 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 formulation 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
PAGE 63
49 written to solve neutron diffusion problems. A good example of such codes is CORA/ ' which was written in 1970 for the solution of onedimensional, fewgroup problems. Typically, 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 + l).Mk,j+l)+C(k,j, j)(^(k, j)+C(k,j,jl)cf)(k, j1) + Z^ (k , j ) (^ (k , j ) + Ej^ (k , j ) 4) (k , j ) +e Zp (k , j ) <, (k , j ) K ^ ZgCnvk, j)(j)^(n,j) n^l ^ = rr^X(k) 7 vZ.(n,j)())(n,j) , (73) ^eff ^ ^ 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 accelerating 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
PAGE 64
50 a single search parameter. This parameter, denoted as e in equations (73) , is varied between successive flux solutions (23) according to a hyperbolic extrapolation algorithm 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 magnitude. This author would like to point out that, by selecting the poison eigenvalue formulation, instead, this same problem can be handled with only one solution of the flux equations. In order to include the option of a user specified reactivity other than unity, equations (44) for the poison eigenfluxes may be written as C(k,j,j + l)(;)(k,j+l)+C(k,j,j)<)(k, j) +C(k, j, jl)(})(k,jl)+E^ (k,j)^(k, j) K Â— Â— y Ci<\ K Â— X(k) ^^'^f ('^':3)'? (n,j)+Zj^(k,j)(},(k,j) eff , n=l Z^(nvk,j)4)(n,j) = eZÂ„(k,j)^(k,j) . (74) n=l With K ^^ and the shape of the poison distribution, Z (k,j), specified, equations (74) may be solved by the pov/er iteration method, with acceleration techniques if desired. The solution will be the flux distribution for the reactor configuration with the desired reactivity, and the eigenvalue,
PAGE 65
51 e, will give the necessary magnitude of the poison distribution. 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 criticality, with a given shape of the poison distribution, is only part of the overall problem of finding the optimum poison 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 distributions v/hich result in a critical reactor. Several numerical 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 distribution 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 pov/er distribution must be computed and the appropriate error value (s) must be calculated. The utility of the poison eigenvalue formulation has already been discussed. Using this approach, the solution to a J'K
PAGE 66
52 order matrix eigenvalue equation must be found for each poison distribution investigated for a reactor represented by J space points and K energy groups. This matrix equation has the form [A]* = e[Z^], (75) where [A] contains all of the terms on the left hand side of equations (74) , and [Z ] 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.. M i=l where M is the order of expansion. If M = G'k, where G is the number of control points, then, as discussed in Chapter IV, the M poison eigenvectors will be sufficient to represent 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] 2Vi = ^^VH^i*! ^^'^ i=l i=l
PAGE 67
53 Since the matrices [A] and [E ] are not functions of the index of summation, they may be included inside the sums. Substitution according to the poison eigenvalue relation, equation (45), allov/s the left hand side of equation (77) to be written as M M i=l 1=1 The [Z ] matrix m.ay be separated into the matrix, [N] , giving p just the normalized spatial distribution of the poison concentrations, and the matrix [a ] , 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 ^f^pl2^i^i = ^[^^HCpl^a.^.. (79) X i 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 i=l 1=1 The orthogonality relations for the poison eigenvectors and their ad joints are found from equation (37) by substituting [a ] for [B] .
PAGE 68
54 ^mf'^p^*! " '5(i,j) (81) Using these relations, equations (80) may be simplified to M ^m^m " e^a^F(m,i, [N] ) , for m = 1 to M, (82) i=l where F(m,i,[N]) = ^^^N] [a ]<})^. (83) 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. [e]a = e[F]a , (84) where [3] is a diagonal matrix, with all of the 6 on the m 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 substitution of a G'K order matrix equation for a J'K order equation appears to offer a substantial savings in computational work for evaluating the power distribution produced by each poison distribution considered. In addition, [3] is diagonal and is more easily inverted than is the [A] matrix of the unexpanded problem. However, the eigenvalue nature of equation (84) presents a problem. Generally, there are as many
PAGE 69
55 as M solutions to the equation; just as there are to the eigenvalue equation (45) used to generate the M poison eigenfunctions. As pointed out in Chapter V, a simple power iteration solution of the equation a/Â£ = [3]~'^[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, v/ith the smallest neutron losses due to leakage, requires the greatest amount of poison to attain criticality, and therefore the most positive value of e. Obviously, this problem also occurs for the unexpanded poison eigenvalue formulation, as given by equation (75) . There are tv/o 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 V/ielandt 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
PAGE 70
56 that are close enough, that this method may be expected to work quite well. In practice, no difficulty has been experienced 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 v/ith a uniform distribution. In this case, [Z ] = io^] , and the unexpanded poison search equation (75) becomes identical to the poison eigenfunction equation (45). Thus, the desired value of e is identical to 3,, the eigenvalue of the first poison eigenfunction. This method for predicting e, v/hile accurate, does not allow any prior knowledge or intuition of the optimum to be used in specifying the starting distribution for the poison search. A third method, v/hich the author has found satisfactory, is to use first order perturbation theory to predict the value of e for the fundamental mode associated v/ith 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]<})^ = e[N] [o ]
PAGE 71
57 Using eigenvalue relationship for 4), , this becomes 6^()^ = Â£[N] [Oplcfi^. (87) Multiplying both sides by the transpose of the adjoint flux, T ^^, gives a scalar equation which may be solved for e to yield e = i>1[m[o^]^^/^^ = F(l,l, [N])/3j_. (88) Since F(1,1,[N]) has already been found for use in the expansion 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 Â£, this author has experienced no difficulty in practice. The disadvantage of the V7ielandt transformation approach is that the transformed equation has the form { [6]eg[F] }a = 6[â€¢]a. (89) Thus, a full M order matrix must be inverted, instead of the simple diagonal matrix [6] , for the power iteration with the transformed operator {[6]g [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 hov^r many spatial
PAGE 72
58 dimensions and energy groups have been used and upon the details of scattering among the energy groups. In particular, if only one spatial dimension is considered, the multidiagonal nature of [A] (and therefore of [A] e [Z ]) lends y ir itself to the use of "inner iteration schemes" that take the (23 24) place of the matrix inversion. ' ' 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 particular case. If more than one spatial dimension is considered, solution of the unexpanded equation becomes much m.ore difficult. 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. In such cases, hov;ever, inversion of the G'K matrix would certainly 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 e is designed to preserve the [B] matrix in a simply invertable form. This may be accomplished by observing the following phenomenon. Any eigenvector, 4)., that is a solution to equation (45) will also be a solution to the equation {[Aj+BgtapDcj,^ = (6^63) [OplcJ). = &[io^]^^. (90) That is, reducing all of the 3j^ by the amount 3g is equivalent
PAGE 73
59 to assuming that the concentration of control poison is already 6 at all control points. If B is made large enough that s s the fundamental poison eigenvalue, 3j_, is negative, then physically, the reactor is shut down. The problem of determining the optimum poison distribution then becomes a matter of removing poison from v/hat is already there. The distributions investigated by the search program will be the complementary distributions necessary to make the actual distributions equal to B^ at every control point. That is, for normalized actual distribution [N] , and complementary distribution [N'], the relationship between the magnitudes necessary for criticality of the fundamental mode is e[N][a^] = e' [N'] [a^]+B^[a ] . (91) p y "^ c The flux matrix equation associated with the evaluation of such complementary distributions, [N']/ is {[A]B,[a^]}$ = Â£' [N'] [a 14.. (92) Using the same set of poison eigenvectors and adjoints and the new eigenvalue relations (90) , equation (92) may be expanded by a method analogous to that used in obtaining the expanded equation (84). [B']a = e' [F']a (93) v/here the elements of [F'] are given by analogy to equation
PAGE 74
60 (83) as F'(m,i,[N']) = iJ^^tN'] [ap](})^. (94) Thus, the value of e' which corresponds to the fundamental mode may be varied by varying 3 . As long as poison is only removed at each point, there is no ambiguity. However, if the base level of control poison, 6 , 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 v/ill be positive at some points and negative at others. Such a distribution may be interpreted two v;ays. VJith a negative e', the trial distribution is to be subtracted from the uniform distribution, but with a positive e', it is to be added. Since the trial distribution has both positive and negative elements, for either case, poison v;ill be added at some points in the reactor and removed at others. If the complementary distribution chosen for investigation is such that both positive 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 situation does not arise, it is prudent to chose the value of Q
PAGE 75
61 so that, for all total poison distributions that may reasonably 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 [W ] itself is of no significance because only the product, e'[N], must be negative to designate removal of poison.) The fact that such a selection of 6 will result in s all possible root values of e' being of the same sign is physically obvious. Since the reactor is shutdown by the Â• uniform poison distribution, g , and since poison may only be removed everywhere or added everywhere by the distribution [N']/ obviously, 3' 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 m.ore buckled flux shapes, and require even more poison to be removed for criticality of their modes. Thus, the other values of e' must be greater in absolute value and of the same sign as e\ Â• This insures that the power iteration solution of equation (91) v;ill always be the fundamental mode and that the distribution [N'] will always be considered as subtractive. Unfortunately, equation (91) does not aid in the selection of a value for 6 that will insure behavior for e' as described above. Since the "worst case" trial distribution, e [N] , will not be known before the search program is
PAGE 76
62 run, some other means must be used to select 6 . Certainly, 6 must be greater than 6, , or the reactor will not be initially shut down. If Q 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 3 also serve to increase the absolute magnitudes of all the g!. This tends to decrease the dominance ratio for Bj , causing more iterations to be required for convergence of the flux for each [N'] investigated. For the general case, this may be hard to grasp, but for the special case [N' ] = [I], it is easily illustrated. In that case, e. = B . , and the subtractive distribution [N] = [I] . This allows equation (91) to be simplified to B Â• = e'+B . (95) ^1 X ^s Thus, the dominance ratio for the el, as a function of B^' is I (B2B3)/ (3]_Bg) I Â• Since B2 < 33_ < 65/ this dominance ratio for e, must decrease as B increases (i.e., as the reactor is further shut down) . In order to use this method for insuring that the proper 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 B Â• If the value chosen is too small, the error may be detected v;hen a value of el v/ith an improper sign is computed, or v/hen [N'] contains negative elements. In such an event, it is necessary that the value
PAGE 77
63 of B 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 3 , that corresponds to the best distribution found so far using the old value of B Â• In general, [N] and [a ] will be diagonal with some zero diagonal elements in positions corresponding to noncontrol points. If these zero rows and columns are eliminated, however, then the resulting matrices will possess inverses. With the understanding that this has been done, the solution for the total poison distribution associated with the best complementary, distribution, e ' (1) [N ' (1) ] , found so far using 8 (1), is Â£:[N] = Â£' (1) [N(1)] + Bg(l) [I] , (96) where the fact that the reduced [a ] is diagonal, and thus commutative under matrix multiplication, has been used to simplify the relation. The new complementary distribution, e'(]) [N'(2)], must represent this same total distribution. Equating the two expressions for e[N] and solving yields e'(2)[N'(2)] = Â£' (1) [N(l)] + {3g(l)Bg(2)}[I] . (97) This new distribution, normalized to [N' (2)], may be used as the nev; starting point when the optimization search is begun again using 3 (2) . This process of retransf orming the 3 set and restarting the optimization search could be done manually or automatically by the search program. If an unnecessarily large value of B is picked, however.
PAGE 78
64 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 3 to increase the dominance ratio. Perhaps a routine could be programmed that would automatically decrease 3 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 tv/o methods described above for the solution to the poison distribution evaluation problem, the first, using the Wielandt transformation, is the easier to program and understand. 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 econoraica 1 , however .
PAGE 79
CHAPTER VII NUMERICAL CALCULATIONS AND THE PROSPECT FOR PARTIAL EXPANSIONS In order to illustrate the use of some of the techniques 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, the principles of which are described in Appendix III. The eigenvector expansion technique, which is used for synthesizing the flux and power distributions for each trial poison *C.D. Kylstra, personal communication. 65
PAGE 80
66 distribution generated by GAUSS, is incorporated in the subroutine 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 poison concentration is calculated by the algorithm ERROR = / {lP(j) }^V(j) , (98) J 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 sura 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 techniques, requiring different formulations of the error function, 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 requirements, the Weilandt transformation approach has been used to ensure that the proper Â£ 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 3 transformation method, it is perfectly adequate for
PAGE 81
67 the purpose of demonstrating the eigenvalue expansion technique. In order to gain a slight simplification in the necessary matrix algebra, the equation solved by the UFUNT subroutine has been rewritten from equation (84) . [F]a = (l/Â£) (6]a (99) Using the Wielandt transformation with a guess at the inverse of the desired eigenvalue, (l/e) , this equation may be written as {[F](l/e)^[3] }a = 6[QU, (100) where 6 = (l/Â£)(l/Â£)^. The matrix operator used for iteration purposes is thus { [F] (l/c) ^ [ B] } [0]. In order to provide a check of the optimization program, a simple case, amenable to tv/o solution techniques, was chosen as the first example. A bare, homogeneous infinite half slab 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 constant) 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
PAGE 82
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 condition at centerplane (point 1) , zero flux condition at extrapolated boundary (point 6) . Materials Parameters (in standard metric units) : Group (k) D(k) Ea(k) ZR(k) vZf (k) X(k) 1 1.46 .0122 .0175 .0126 1 2 0.37 .1841 .2666 lattice spacing = 8 cm Control Poison Cross Sections: a (1) 0.1 cm_ per grara/cc a (2) = 0.2 cm per gram/cc
PAGE 83
69 10 20 30 X (cm) 40 20 30 40 X (cm) Figure 2 Poison Eigenvectors ((),, (J),, 4)^ and (})_ for Example 1
PAGE 84
70 'It A, O A A 'o Â° Fast Flux ^ Thermal Flux 10 20 30 X (cm) OA 10 20 X (cm) I 40 a 30 40 10 20 30 40 X (cm) Figure 3 Adjoint Vectors \p. , ^, ^^ and ip^ for Example 1
PAGE 85
71 TABLE 2 The Eigenvalues Associated with the Poison Eigenvectors for Example 1 Mode
PAGE 86
72 The behavior discussed in Chapter V for the poison eigenvalue series may easily be seen in this example. Because the fundamental mode in the unpoisoned reactor is supercritical, g, is positive. All the other 8 grow more negative with increasing index. The 3 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 3< 6, Â• 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 techniques 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 pov/er distribution vector was associated v/ith 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 equation P^(j) = A{Z^(l)*^(l,j)+ZÂ£(2)4)]^(2,j) } = A{Z^{1)S'[T.^{2)}^1{2,J) . (101)
PAGE 87
73 TABLE 3 Optimization Results for Example 1 Results of General Poison Distribution Search Method Point No. 1 2 3 4 5 Fast Flux 75802E 75777E 75518E 73029E 49391E 01 01 01 01 01 Thermal Flux 44744E 44756E 44878E 46055E 57227E 00 00 00 00 00 Power Density 1, 1, 1, 1, 1 Poison Concentration 6.5S350E02 6.54151E02 6.41074E02 5.16011E02 1.02413E01 Results of Power Vector Expansion Coefficient Search Method Real Solution: Point No. 1 2 3 4 Fast Flux 75802E 75778E 75518E 73029E 01 01 01 01 Thermal F lux 47444E 47555E 48788E 40555E 00 00 00 00 2.49391E 01 2.52266E 00 Power Density 1.0 1.0 1.0 1.0 1.0 Poison Concentration 6.55372E02 6.54153E02 6.41076E02 5.16005E02 1.02412E02 Mode Index 1+ 2 + 3 + 4 + 5 + Expansion Coefficient ~2.7 2 00rE 4 2 8 81539E 88696E 71797E 85565E ^1 01 01 01 01 Mode Expansion Index Coefficient T^ 3.82306E~0"2 21.49576E03 31.65447E01 42.21372E01 59.48818E02 Second Mathematical Solution (Not Physically Point No. 1 2 3 4 5 1, 1, 1 9. 1. Fast Flux 19011E 12463E 23898E 86252E 54869E 01 01 01 00 01 Thermal Flux ' 31338E 28244E 33649E 21704E 48286E 00 00 00 00 00 Power Density 1.0 1.0 1.0 1.0 1.0 Possible) : Poison Concentration 1, 1, 1, 1, 1, 16233E 14970E 17165E 12248E 25350E 00 00 00 00 00 Mode Index 1+ 2+ 3 + 4 + 5 + Expansion Coefficient 4, 4 4, 4 4 48286E 21704E 33649E 28244E 31338E 02 02 02 02 02 Mode Index 12345Expansion Coefficient 1 9 1 1 1 54869E 86252E 22898E 13463E 19011E 03 02 03 03 03
PAGE 88
74 A similar series of power vectors, P~, is defined for the (f)j^ series of modes. Since ({>. and <^. have the same spatial shape, and differ only by their fast/thermal flux ratios, S. and S., their associated power vectors, P. and pT, differ only in magnitude. That is, P. = C.pT. J J 'ill The "desired power shape" is everyv;here constant, and, when normalized to a unit pov;er 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 tha poison eigenvectors. 1=1 then the desired pov;er vector, P*, is given by 5 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* = 2^i^i' ^^^^^ i=l The vector equation (104) may be solved by simple matrix inversion to obtain the b.. By equating coefficients for P. in equations (103) and (104), the following relations are derived.
PAGE 89
75 (a"!'+aT/c.) = b. (105) 1 1^ 1 1 If the five new variables, d., are introduced such that 4 = b^/d^, (106) then, from equation (105) , aT = c^b^(d^l)/d^. (107) Thus, the choice of the five variables, d. , allows the computation of a set of ten flux expansion coefficients, a., v/hich allow the construction of a flux satisfying the flat power condition. With these coefficients and the poison eigenvectors, the poison concentrations necessary for criticality 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 eN(k,j) = 2a^6j_<'i(k,j)/2]^i
PAGE 90
76 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 ij; [o ]
PAGE 91
77 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 multigroup diffusion theory program in which a uniform poison concentration, Bi, 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 C0DE2 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 v/ithin 10 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 optimization v;as 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
PAGE 92
78 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 extrapolated boundary (point 11) Materials Parameters (in standard metric units) : Group (k) D(k) Z^(k) Ej^(k) vZ^ (k) X (k) 1 1.461860 0.0103369 0.0169434 0.005872 1 2 0.366341 0.0979940 0.135760 lattice spacing = 15 cm Control Poison Cross Sections: ,_ a (1) = 2.86715 cm1 per gm of B /cc a (2) = 130.693 cm~^ per gm of B''"Â°/cc To achieve regionv/ise 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.
PAGE 93
79 w
PAGE 94
80 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 regions to either side, weighted by the fraction of the reactor 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 reactivity for this materials configuration is 1.0000114, v;hich 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 variation may lead to more than one "optimum realizable pov/er 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 underspecif led , 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
PAGE 95
81 TABLE 5 ThreeRegion Poison Optimization Results for Example 2 Poison
PAGE 96
82 9>n
PAGE 97
83 one of these may be eliminated by physical considerations, however.) If the problem becomes underspecif ied, the hyperspace has fewer dimensions and the error surface may develope extra local minima. This is a phenomenon that is caused by the physically underspecif ied 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,
PAGE 98
84 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 (}>. and 4). 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 poison eigenvectors. The c^>, and t^, 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. However, 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
PAGE 99
85 TABLE 6 Reactor Model for Example 3 Core Configuration: infinite cylinder, two homogeneous fueled regions from radii 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) I 2 V7estinghouse PV7R fuel assemblies, 3% enrichment, 18,000 MV?D/T D(k) E (k) 1.466676 0.364366 .01083982 .09871243 .01648968 vZ^(k) 00542511 1305982 X(k) R egion II ; Group (k) 1 2 Westinghouse PWR fuel assemblies, 10,000 MWD/T D(k) 1.46T8'62 0.366341 .01033693 .09799408 R (k) .01694346 3% enrichment, vE^(k) X(k) .00587230 .1357601 1 Reflector Group (k) I 2 pure water D(k) Z (k) 1.136999 0.149380 0005892 0019240 .04830 vEÂ£(k) ~o7o 0.0 X(k) lattice spacing = 10 cm for all regions Control Po Gp(l) a;(2)
PAGE 100
86 Region I / .o 0Â— oRegion II ^^0Â— Q^O 'Q Â©^ I Reflector ./ / ^ Â— AAA ^ ._AÂ— AÂ— AÂ— AAÂ— AA I 50 140 Radius (cm) 190 Region I AÂ— AÂ— AÂ— < ^Â— AÂ— AÂ— / ^Â—^Â— AÂ— AÂ— AÂ— AÂ— A AÂ— vV I Region II I Reflector ./i 'V/^ 0^^ / / / / / Â° Fast Flux A Thermal Flux / / _o0Â®' Â© 60 50 140 190 Radius (cm) Figure 6 Poison Eigenvectors 4), and (J), for Example 3
PAGE 101
87 + ^1 Region I A/ % Â— oO ^Â— AÂ— AA, ,/ / A / Region II l\\ Reflector / 50 140 <>Â—^ 190 Radius (cm) AÂ— AÂ— AÂ— AÂ—, AÂ— AAÂ— ^__ A^ Region I Region II / Â® Fast Flux A Thermal Flux ,0 / y i_oo ^oo I 50 140 Radius (cm) 190 Figure 7 Adjoint Vectors ^, and t' for Example 3
PAGE 102
88 modes are expected to be small, they may still be present. In order to give the best opportunity for successful solution with partial expansions, the problem was designed so that the interface between the two poison regions concides with the interface between the two fueled regions. This should allow the use to maximum advantage of any representation of an interface at that point which may already be built into the poison eigenvectors. Except for the success of the partial expansions, however, no evidence is offered that there is indeed any advantage to such a coincidence of interfaces. Optimum tworegion poison concentrations were calculated for Example 3 using complete and several partial expansions of the problem. The resulting "optima" are listed in Table 7 as a function of the number of vectors used. As may be seen from these results, the tworegion optimum poison concentrations may be calculated v/ith an accuracy of better than 5 percent using as few as five of the thirty eigenvectors. Table 8 lists the expansion coefficients found for the optimum of each expansion. As can be seen, there is little variation in the modes which give the most significant contributions to the flux shape. This is in large part due to the grossly underspecif ied nature of Example 3. Considerable compensation by the remaining modes for those that have been deleted v/ill be seen when partial expansions of the completely specified Example 1 are considered.
PAGE 103
89 TABLE 7 Comparison of Reactor Configurations Resulting from the TwoRegion Poison Optimization of Example 3 Using Full and Partial Expansions Expansion Designations: A full expansion using complete set of 30 vectors, 1+ through 15+ and 1through 15B partial expansion using 20 vectors, 1+ through 10+ and 1through 10C partial expansion using 10 vectors, 1+ through 5+ and 1through 5D partial expansion using 5 vectors, 1+ through 5+ Solutions for Optimum Two Region Poison Concentrations: Region Exp. A Exp. B Exp. C Exp. D I 1.1740E05 1.1680E05 1.2306E05 1.2310E05 II 6.1978E05 6.1956E05 6.1935E05 6.1936E05 Power Density as a Function of Position Point Exp. A Exp. B Exp . C Exp. D 1.28071 1.35300 1.35318 1.26381 1.33234 1.33249 1.23678 1.27517 1.27524 1.22184 1.19734 1.19736 1.17008 1.12347 1.12348 1.09322 1.09448 1.09452 1.15118 1.14727 1.14730 1.07691 1.10928 1.10936 1.02184 1.03187 1.03196 0.97586 0.96187 0.96192 0.92327 0.91062 0.91061 0.87441 0.87134 0.87130 0.82418 0.83118 0.83114 0.77357 0.78294 0.78291 0.76262 0.76405 0.76403 1.94757 1.93433 1.93423 Worst Error Relative To Full Expansion Solution: Exp. A Exp. B Exp . C Exp. D 0.0% 0.8% 6.5% 6.5% 1
PAGE 104
90 TABLE 8 Comparison of the Expansion Coefficients for the Full and Partial Expansion Optimums for Example 3 Sector
PAGE 105
91 Before leaving Example 3, it is instructive to examine the "optimum realizable power distribution" as computed by the GAUSS program. The shape found for the full expansion solution is plotted in Figure 8. This shape is "optimum" only in the weighted least squared errors sense. That is, it gives the minimum value of the error as computed by equation (86) . This is the form required by the Gaussian optimization technique in order to predict an optimum, but it obviously is not representative of the actual limitations placed on the pov/er distribution, or of the actual penalty associated v;ith a power peak. In practice, it is the maximum power density which must be limited to some safe value, not the average. Any other peaks which come nearly to the same value only serve to further increase total power output. Therefore, one might suspect that the true physically optimum tworegion poison distribution would result when the concentration in region I is decreased to the point that the central peak at point 1, or the interface peak at point 6, coincides with the "twogroup peak" caused by the reflector at point 15. Thus, an optimization technique utilizing a more realistic error function is certainly to be desired. The purpose of this study, hov/ever, has not been to select an optimization technique. It has been, rather, to introduce a powerful technique for the evaluation of the power distribution produced by each poison distribution specified by whatever optimization technique an investigator chooses
PAGE 106
92 o H uoTq.ncftaq.STa :isMOd pazTx^uiJON o Â•H +Â» .Q Â•H V4 +J 0) Â•H Q O (0 Â•H O (U O Â•H tn Â« O 15 d g Â•H P 014 o p o m o Â•H +Â» Â•H 4J w n Â•H Q Q) gj g O X Pa W 00 0) Â•H
PAGE 107
93 to use. In that context, the GAUSS program provides as good a demonstration tool as would a technique providing a more realistic optimum, in that the fluxpower construction methods employed are independent of the error function used. In order to investigate the possibilities for partial expansions in completely specified problems, tests were made using the m.athematical model of Example 1. These tests are not the most general in nature, due to the separability of the energy and space variables and to the small number of points used to represent the space dimension. Hov/ever, tests on a model utilizing a large number of control points would require prohibitively expensive amounts of computer time for the GAUSS program to independently manipulate the poison concentration at each control point. Also, inclusion of an interface between tv;o fueled regions creates an inherently underspecif ied problem (assuming a mesh point is located on the interface.) This is due to the discontinuity of the power density at such interfaces. Two values for the power density must be calculated for the point, while, at most, only one poison concentration may be specified to control them. Example 1 seems sufficient at least to represent two important properties for partial expansions of completely specified problems. Because the series of eigenvectors with negative fast/ thermal flux ratios played such a small part in the solution of the fully expanded problem, the necessity of their inclusion
PAGE 108
94 in the expansion may seem in question. A solution was found for the problem as expanded using only the positive flux ratio eigenvectors, and the results are compared with the full solution in Table 9. Because the five vectors used are still sufficient to expand any power vector, it was still possible for the GAUSS program to find a set of poison concentrations v/hich reduced the error function to zero. However, the resulting poison concentrations are significantly different, especially near the outer boundary, where the spectral shift is greatest from the undisturbed system. Still, these values are accurate within about 6 percent. Even for completely specified problems, then, the elimination of the series of negative flux ratio eigenvectors leads to a reasonable approximation of the optimum poison concentrations. Hov/ever, the examples considered in this study are special cases in a way that has not been mentioned as yet. They are all twogroup representations of thermal reactors, and, as such, allow little variation in the energy spectrum of a critical reactor. All source neutrons appear in the fast group, and, due to its much larger cross section, most fissions are produced in the thermal group. In a fast reactor, the spectrum is not dominated by these phenomena, and introductions of poison material may cause gross spectral shifts. In such cases, elimination of the negative flux ratio eigenvectors may lead to extremely erroneous results. Another phenomenon which may show up in partial expansions
PAGE 109
95 TABLE 9 Results of Partial Expansions for a Completely Specified System (Example 1) Expansion Designations: A partial expansion using 5 vectors, (^t^ through (p't B partial expansion using 4 vectors, (p't through
PAGE 110
96 is compensation by the remaining modal components for those that have been deleted. For instance, in a case where only the positive flux ratio eigenvectors have been used, except that for one vector, (p. has been substituted for 4)., the coefficient for 4). would be greatly increased over Its value m a full solution. This is because the shape, c}) . , is necessary to flatten the power distribution, and it will be used by the GAUSS program, if available, no matter what happens to the spectrum of the flux. Obviously, an intelligent approach must be taken to providing a set of eigenvectors to be used in partial expansions. The vectors which are expected to give the smallest contributions to the solution should be eliminated first. In addition to the change in the predicted optimum values, there is also a decrease in the ability of the GAUSS program to find this optimum as the number of vectors used in the expansion is decreased. This is due to the fact that the power calculated in the expanded equation no longer completely represents the power in the unexpanded formulation. The "error hypersurf ace" in "concentration space" is made increasingly more rippled, until many local minima result. It becomes increasingly likely that the GAUSS program will not be able to distinguish these false minima from the one associated with the physical optimum. Finally, the rippling of the error surface becomes so exaggerated that the optimization program can no longer calculate
PAGE 111
97 meaningful derivatives for the prediction of the optimum point. In general, this phenomenon may be expected to become important sooner as the problem becomes more nearly completely specified. Partial expansions of Example 1 were carried out for several sets of eigenvectors which are not capable of completely representing the "desired power vector". The results are also given in Table 9. These solutions were all terminated by the GAUSS program due to an inability to find meaningful derivatives, rather than due to the minimization criteria for the error value. As may be seen in the table, the expansions using the first three and four positive flux ratio eigenvectors provided poison distributions in the general vicinity of the true values. However, the expansion using the first four negative modes, as well as the first four positive modes, provided a much less satisfactory solution. Although it provided the same four shapes as the expansion using only four positive ratio eigenvectors, it provided twice as many coefficients to be found in solution of the expanded equations. This apparantly gave more occasion for the coefficient compensation phenomenon, producing less valid representations of the flux and power than did the expansion using only positive flux ratio modes. The lesson to be drawn from this experience thus seems to be to eliminate the modes in the order that they contribute to the desired solution, from smallest contribution to largest. Since the solution is not known beforehand, this
PAGE 112
98 method requires some experience and insight, as does any approximation technique. The negative poison concentrations occassionally found in the solutions to these examples represent a problem encountered in any automated poison distribution search. If there is already some poison in the reactor, then the negative concentration may not exceed this magnitude. If there is no initial poison, then the optimum poison concentrations must be constrained to be positive. Because this study is directed at the evaluation of the effects of a trial poison distribution, once selected, rather than the actual selection of the distributions for trial, this problem will not be considered further. There is one more phenomenon of the poison eigenvectors which the author yet wishes to demonstrate. Example 4 is included to illustrate the concepts advanced in Chapter V for the identification of the eigenvectors. The reactor model is a twogroup, fortyfourpoint representation of a cylindrical reactor with three fueled regions and a water reflector. The reactor model is described in Table 10, and shown schematically in Figure 9. The poison eigenvectors (f), and (})Â„ are given in Table 11. All of the even numbered points within fueled regions are control points, giving rise to 32 poison eigenvectors. The fundamental mode, ({)Â•,, is easily identified in that it is everywhere positive. However, without the concept of the region of completeness, it may
PAGE 113
99 TABLE 10 Reactor Model for Example 4 Core Configuration: finite cylinder, three homogeneous fueled regions from radii to 90 cm, 90 to 130 cm, and 130 to 160 cm, pure water reflector from radius 14 to 22 cm, 366 cm core height Mathematical Model: diffusion theory, 2 energy groups, 44 space points, 16 control points, zero current condition at centerline (point 1) , zero flux condition at extrapolated boundary (point 45) , continuous neutron current conditions at internal interfaces (points 19, 27, and 33) Materials Parameters (in standard metric units) : Region I : Westinghouse PWR fuel assemblies, 3% enrichraent, 18,000 MWD/T Group (k) D(k) Z^(k) Z^{k) vE^(k) X(k) a (1) = 2.86715 cm"! per gm of B,q/cc aM2) = 130.693 cm"! per gm of B /cc 1 1.466676 .01083982 .01648968 .00542411 1 2 0.364366 .09871243 .1305982 Region II : V>Jestinghouse PV7R fuel assemblies, 3% enrichment, 10,000 MWD/T Group(k) D(k) E^(k) l^{k) vZ^ (k) X(k) 1 1.461862 .0T03T693 .01694346 .00587230 I" 2 0.366341 .09799408 .1357601 Region III ; Westinghouse PWR fuel assemblies, 3% enrichment, MWD/T Group(k) D(k) Z^(k) Zj^(k) ^^^(k) X(k) 1 1.455021 .00949965 .01768530 .00648183 1 2 0.374245 .08487368 .1307898 Reflector: pure water Group(k) D(k) E^(k) Ej^(k) vE^ (k) X (k) 1 1.136999 .0005892 .04830 0.0 =~ 2 0.149380 .0019240 0.0 lattice spacing = 5 cm for all regions Control Poison Cross Sections ,p(2) = 130.693 cm~Jper gm Control points are those with the even indices 2 through 32 B^ = 7.380E05 cm~^ z
PAGE 114
100 r 3 o H H fn P Â•H O T3 d) O O Â•H +J H 'Ho o 0) o 0) u
PAGE 115
101
PAGE 116
102 g U Xi I in oiH (M (N ro o o o o I I I I w w w w in "^ iH 00 rH PI en O fN rH 00 a\ in n ^ o vO ^ [^ rH .HCNOOOOOrH oooooooo I I I wwwwwwwu >^u3vDOcNoocricrÂ» ocrrofNuocTiooin CTiOCNOJUJin'^rH crv'^cNrHin'^i'inm rH iH rH (N OJ (N o o o o o o W PJ w W W w iH in vD ^ as in in (N in fN 00 ^ CO in vD CO cN Â•Â•^ in "^ 00 CN rH CN ^ iH rCTt CN CN o) rH rin IÂ— icNrH^rofNiHcninroiHcrNincN I I I I i I I I I I I I I P M X (13 3 hi I in eOOrHrHiHOOiHrHfN OOOOOOOOOO III II wwwwwwwwww oon OOU3iHCOrOOVOoOrHr~ 04 nH
PAGE 117
103 be difficult to identify ~. If only the flux values at the control points are considered, however, it may be easily verified from Table 11 that these components exhibit five sign changes in each group and a negative fast/thermal flux ratio everywhere. The vector is thus readily identified as (})j. Finally, some mention should be made of the accuracy possible in using the eigenvector expansion technique. The eigenvectors and adjoint vectors generated for use in Examples 1 through 4 have generally satisfied the orthogonality condition very well. The cross products ii"^ [o U. result in values of 10~^ to 10~^Â° v/hen m ^^ i . The use of M eigenvectors for a complete expansion may, the author estimates, introduce an additional fractional error on the order of /M x lO"^, which is usually rather small compared to the errors associated with the materials parameters and the problem formulation in descrete point, descrete energy group diffusion theory. In using partial expansions, the error will obviously be significantly greater. The accuracy xn such cases may only be gauged by approaching the observed minimum from several directions, to insure that it is indeed the optimum, and by increasing the number of vectors in the expansion to see how this affects the location of the minimum.
PAGE 118
CHAPTER VIII CONCLUSIONS AND SUGGESTIONS FOR FURTHER STUDY In this study, it has been shown that the specification of a physically realizable power distribution gives no a priori insight as to the spectrum of the multigroup fluxes which must produce it. The spectrum must therefore be considered as a set of unknown functions. The spatial distribution of some control material (with known energy dependence of its cross section) must be found v;hich maintains the desired power shape in the critical state. The problem is completely specified, leaving no freedom for arbitrary selection of some of the unknown functions, and a solution is knov/n to exist from physical considerations. Unfortunately, the multigroup criticality equations for the controlled reactor are nonlinear. The control terms contain the products of the unknov/n poison distribution and the unknown spectrum functions. This nonlinearity has prevented the author from obtaining an analytic solution to the general problem, and forced him to consider the approach of using numerical optimization techniques to find solutions to particular problems. The author has introduced a method which he believes will be of interest and utility to those who seek to optimize 104
PAGE 119
105 reactor power distributions through poison management. In comparison to the conventional method, this method considerably reduces the amount of work necessary for adjusting the magnitude of a poison distribution to achieve criticality and for evaluating the resulting power distribution. It is this type of calculation which must be done repeatedly, for different poison distributions in the same reactor, during numerical searches for the optimum poison distribution. This technique may be considered in two parts. First is the poison eigenvalue formulation of the criticality equation and a suitable transformation technique to ensure convergence to the proper eigenvalue. This formulation allov/s for the simultaneous adjustment of the magnitude of the poison distribution to the critical value, and the calculation of the flux shape in the critical reactor. Only one solution of the criticality equations is necessary for the effects of each poison distribution to be investigated. This is a considerable simplification over the conventional, reactivity eigenvalue formulation, where the criticality equations must be solved many times in order to adjust the magnitude of a poison distribution to achieve criticality. An additional advantage of the poison eigenvalue formulism is that it separates the control terms from the rest of the terms in the criticality equation. In this formulation, it is the matrix containing the noncontrol terms which must
PAGE 120
106 be inverted in order to perform the power iteration for the flux. Since this matrix remains unchanged no matter how the control terms are varied, only one matrix inversion is necessary no matter how many poison distributions are investigated. (Improper choice of the method for ensuring convergence to the proper eigenvalue may destroy this property, however.) In the conventional, reactivity eigenvalue formulation, the control and some noncontrol terms are mixed on the same side of the criticality equation, and a matrix inversion (or equivalent) is necessary each time a new flux vector must be found. Two transformation methods have been discussed which ensure convergence to the proper poison eigenvalue. One of these, the V7ielandt transformation, has been demonstrated successfully by num.erical examples. It was chosen for demonstration purposes because it is easily programmed. However, it has the undesirable feature of mixing control and noncontrol terms on the side of the criticality equation which must be inverted. A matrix inversion is thus required each time a new flux vector is found. The other method, the 6 transformation, has been advanced as offering similar benefits without the penalty of additional matrix inversions. Because of the complexity of the programming required to overcome its inherently slov/ convergence rate, this latter technique was not demonstrated. The second part of the technique involves the expansion
PAGE 121
107 of the criticality equation in a set of eigenvectors which individually satisfy the boundary and interface conditions required of the flux in a particular reactor. This set of eigenvectors must be generated for each reactor for which the optimum poison distribution is to be found. They are then used to repeatedly expand the criticality equation, as different poison distributions are added to the uncontrolled materials configuration. The orthogonality relations and regions of completeness were investigated for various set of eigenvectors, which satisfy various eigenvalue equations associated with the criticality equation. The fission eigenvectors v;here shown to be incomplete in the energy dimension, and thus unsuitable for the expansion of a general multigroup flux vector. The poison eigenvectors v;ere introduced and shown to be complete in energy and space for the flux elements at control points, only. It was then proven that this vector set is capable of representing the flux at all points in the reactor for any critical configuration which may be produced by the poison material at the control points. The poison eigenvectors were chosen as the set to be used for the expansion, due to the great simplifications that result from their use. The first of these simplifications is the reduction of the order of the matrices in the criticality equation. The unexpanded problem has an order given by the product of the number of spatial points and the
PAGE 122
108 number of energy groups used to model the reactor. The poison eigenvector expansion reduces this order to the product of the number of control points and the number of energy groups. In addition, the matrix which contains all but the control terms of the criticality equation is reduced to diagonal form. This is true for all reactor models, whether they employ onedimensional or multidimensional representations of the spatial coordinates. Because it is this matrix which must be inverted in order to perform the power iteration solution for the flux, its reduction to diagonal form has the advantage of making this inversion trivial, In essence, the generation of the poison eigenvectors and their adjoints has also provided for the inversion of this matrix. The use of the B transformation, rather than the s Wielandt transformation technique, to ensure convergence to the proper poison eigenvalue, will preserve the diagonal nature of this matrix in the expanded equation. These simplifications are possible because the poison eigenvectors are a complete set which satisfy the criticality equation, reflecting the geometry, the various regions of different materials, the energy dependence of the poison cross sections and the positions at which poison may be introduced. The technique, as a v;hole has been demonstrated by onedimensional, twogroup numerical examples, and the Gaussian least squares optimization method has been found capable of handling the numerical search. These examples have shown the feasibility of the method and allowed comparison of
PAGE 123
109 flux shapes and reactivities with calculations made by a standard diffusion theory code. Agreement with the standard code was excellent for these parameters. In addition, these examples verified the necessity of using a set of eigenvectors which are capable of expanding both the energy and space dimensions of the flux values at the control points. In cases where the most general variations of poison concentrations are to be allowed, it v;as shown that all of the eigenvectors gave significant contributions to the optimum flux shape. In cases where the concentrations are to be somewhat constrained, for instance to regionv;ise variation, the use of relatively few eigenvectors in the expansion has led to quite reasonable approximations to the optimum poison distributions. The method, as presented in this study, offers considerable power for the solution of the general problem, v;hile at the same time allov/ing considerable simplification for approximating the less troublesome cases. Before this technique may be considered a practical tool for the reactor engineer, hov;ever, it will be necessary to embody these concepts in a generalized, wellprogrammed and welldocumented computer code. Toward that end, the author would like to propose the follov/ing areas for investigation. First, the 6 transformation technique should be programmed and tested for reliability and speed. If it can be made to perform as expected, it should be substituted for the
PAGE 124
110 Wielandt transformation technique in ensuring convergence to the fundamental m.ode. In addition, a study should be made of the applications of various numerical optimization methods to this particular approach for choosing poison concentrations in a reactor. The goal should be to choose the method exhibiting good compatability in handling the variables, a form for the error function which most closely approximates the physical situation, and, of course, an approach sufficient to handle the nonlinearity so as to always obtain a solution. Included in this study must be the selection of some means for limiting the search to physically realizable values for the poison concentrations. Certainly, no optimization method will be perfect in all respects. Hov/ever, the proper optimization method, incorporated in a code designed specifically for applying the techniques discussed in this study, should be a quite pov;erful tool in the hands of the reactor engineer. The utility of such a code v/ill, of course, depend upon the availability of easily handled, welldocumented codes for the generation of multigroup, one or more dimensional poison eigenvectors for various reactor geometries. The principles for writing such codes are v/ell in hand, today, and the growing abundance of time sharing, rem.ote computer terminals should offer ample opportunity for using the online operator interaction techniques advocated in this study, Another facet of the overall technique may be of interest v/hen considering partiall; expanded problems with
PAGE 125
Ill regionv/ise poison variation. This involves using an approximation to the set of optimum concentrations, resulting from the solution of a partial expansion, to generate a new set of eigenvectors. By replacing the constant distribution v/ith this "approximately optimum" distribution in generating the vectors, it may be possible to obtain a new set v/hich is better able to represent the true optimum v/ith the same order of partial expansion. In this manner, two sequential partial expansions may prove to offer better accuracy with less effort than would a single partial expansion using tv/ice as m>any vectors. The total number of eigenvectors found would be identical in both cases. However, the order of the matrix equations to be solved in evaluating each poison distribution would be only half as great for the sequential expansion case. This approach also allov/s for an intuitive approximation of the optimum distribution to be used, in place of results from a previous partial expansion. Sufficiently accurate intuition may thus save even more v;ork. It must be noted, however, that the parameters of the numerical search would no longer be true concentrations. They would be a set of pointwise multipliers for changing the concentrations used in generating the new set of eigenvectors. This and other specializations of the general technique may be v/orthy of investigation for developing the full power of the method for use in various special classes of the poison distribution problem.
PAGE 126
APPENDICES
PAGE 127
APPENDIX I A DESCRIPTION OF THE TECHNIQUES EMPLOYED BY THE EIGENVECTOR GENERATING PROGRAM "C0DE2" This computer code was written purely for demonstrational purposes, and is therefore much less general and sophisticated than would be desirable for a standardized code that is designed to be an engineering tool. The program is limited to the calculation of twogroup, onedimensional poison eigenvectors in axisymmetric cylindrical geometry using diffusion theory. Mesh points are placed on the axis, at materials interfaces, at the extrapolated boundary, and at any number of points inbetween. The interval between mesh points must be constant v/ithin each material region. Each region is considered to be of homogeneous composition. Control points, as defined in Chapter IV, may be specified at any or all mesh points, v/ith the exception of the extrapolated boundary point, for which flux values are never explicity calculated. The zero flux condition is used for the extrapolated boundary point, and the zero net current condition (cylindrical axial symmetry) is used at the axis point. Finite axial height of the reactor may be approximated by the specification of a constant axial buckling, Z The derivation of the finite difference approximation 113
PAGE 128
114 from the multigroup analytic equations is well documented (24) m the literature. The derivative therefore will not be repeated here. The appropriate equations for the general problem outlined above, written for point j located at the interface between materials regions I and II, are C ( 1 , j , j + 1 ) 4> (1 , J + 1 ) + C ( 1 , j , j ) ({) ( 1 , j ) C(l,j,jl)(})(l,jl)F(l,j)(J)(2,j) = eG(l,j)(}>(l,j) (llOa) for the fast group, and, for the thermal group, C ( 2 , j , j + 1 ) (}) ( 2 , j + 1 ) + C ( 2 , j , j ) (J) ( 2 , j ) C ( 2 , j , j 1 ) (1) ( 1 , j 1 ) F(2,j)({)(l,j) = BG(2,j)<})(2,j), (110b) where C(k,j,j1) = 7TD^[r(jl)+r(j)]/Cr(j)r(jl)], C(k,j,j+1) = TTDjj[r(j+l)+r (j)]/[r (j+l)r(j)] , C(k,j,j) = C(k,j,j1)+C(k,j,j+1) +Vj(j) (Z^+B^d'^X^vzJ^)j +v^^(j)(z^bVx^zJ)^^, Vj(j) = TT[3r(j)^2r(j)r(jl)r(jl)^]/4 Vjj(j) = 7r[r(j+l)^ + 2r(j)r(j+l)3r(j)^]/4,
PAGE 129
115 F{l,j) F(2,j) G{k,j) V^(j)[vl2l^+V^^(j)[vz2j^^, = V,(j) [ZJ.^V, fj)[zi],,. R'l II ap(j)[V^^V^,]. R' II' In this formulation, the reaction rates are equivalent to what would occur if the flux at point j was the value for the flux everywhere in the cell about mesh point j . (21) Wachspress introduced a technique for solving group diffusion equations which involved writing the equations (110) in "supermatrix" form, where the elements of the supermatrix are themselves matrices. Because this supermatrix is tridiagonal for onedimensional problems, the Choleski algorithm , extended to matrices with matrix elements, may be used in the power series iterations in place of direct matrix inversion. The supermatrix notation of Wachspress has the form [Ul]
PAGE 130
116 [Wj] = [Yj] [Gj] *j = C(l,j,j1) C(2,j,j1) C(l,j,j+1) C(2,j,j+1) G(l,j) G{2,j) 4>(l/j) A Wielandt transformation applied to this formulation does not change the structure. It only replaces the [U j ] matrix with the matrix tUj]6 [Gj ] . Applying the Choleski algorithm to the tridiagonal supermatrix, the supermatrix is factored into upper and lower triangular forms. [XI] [ ] [ ] ... [Wl] [X2] [ ] . . . [ ] [W2] [X3] .. . [I] [Zl] [ ] . [0] [ I ] [Z2] . [0] [ ] [ ] . (^1 *3' v/here = 3 [Gl] [ ] [ ] .. [ ] [G2] [ ] .. [ ] [ ] [G3] .. <^1 02 [Zj] = {[Uj][Wj] [Zj1]}~^[YJ], [Xj] = [Yj]{[Zj]}"^ = [Uj][Wj] [Zj1] (112) (113) (114)
PAGE 131
117 The index j in the above two series may be used from 1 to J (the maximum number of space points) if the definitions are made that [Al] = [0] and [CJ] = [I] . When equation (112) is used for iteration purposes, the known starting vector is substituted for the flux vector on the right hand side to begin the iteration. Thus, the left hand side becomes a known supervecter with vector elements, vj. The flux vector on the right hand side of equation (112) becomes the once iterated supervector with elements vj' . If another (unknown) supervector, v;ith vector elements Tj , is defined by the equation [I] [Zl] [ ] ... CO] [I] [Z2]... [0] [ ] [ I ] ... Â• Â•
PAGE 132
118 now known, the Vj ' may be found from equation (115) by the algorithm VJ' = TJ and Vj = Tj+[Zj]Tjl. (118) ThuS/ the new vector may be found by applying first algorithm (117) and then algorithm (118) . The only matrix inversion necessary is the inversion of the two by tv;o matrices { [Uj][Wj ] [Zj1] } when finding the [Zj] matrices from equation (113) . In general, for a Kenergygroup formulation, only K by K matrices must be inverted, and this need be done only once per set of iterations to find one mode. Once the nev; vector is found, it is normalized to unit magnitude and the iteration is performed repeatedly until only the dominant mode remains, as discussed in Chapter V. The C0DE2 computer program, which uses these techniques to find the poison eigenvectors, is listed in the next appendix. It is written in the Fortran II language for use on the IBM 18 00 computer. This machine v/as chosen because of its availability for online human interaction with the calculations while they are in progress. The desire for such capability is explained in Chapter V.
PAGE 133
APPENDIX II LISTING OF THE DEMONSTPJ\TION POISON EIGENVECTOR GENERATION PROGRAM "C0DE2" Partial Dictionary to the Computer Code Variables: C0DE2 Variables Name or Variable in Text A1(J) fast elements of the [Wj ] matrices A2(j) thermal elements of the [Wj ] matrices B(M/N,J) [Uj] matrices BZS b CI (J) fast elements of the [Yj] matrices C2 (J) thermal elements of the [Yj] matrices EV e t'l(J) fast flux vector F2(J) thermal flux vector FAl ( J) fast adjoint vector FA2 (J) thermal adjoint vector Gl (J) fast elements of the [G j ] matrices <^2 (J) thermal elements of the [Gj ] matrices GAM(M>N,J) [Zj] matrices INTER(I) list of interface indices NC(I) list of control point indices NCP number of control points NP number of flux points NR number of regions 01 (J) fast elements of the Tj vectors 02 (J) thermal elements of the Tj vectors RAD (J) list of radii at flux points SPl al P SP2 Op W(M,N) working 2x2 matrix ^l(J) working vector for fast components ^2 (J) working vector for thermal components 119
PAGE 134
120
PAGE 135
121 a. ro X IÂ•z. o U' K^ to C a. to < o CD CNl LiJ < d X Io U3 II In < O Ll. .00 . CNl 00 o LU 2 o LTl a. o cj c r: o rH II N Â— II CsJ JÂ— ::i w < O O LU o 12 a: Cl. CD ^'^ Â• iH O C il Q. LU ^^ O Q 2 O Â• X S 00 q: ^~y r^ T^ xO Â— Clto :^ o Cl. < ^> LU c^ d: o LA ^ P ^ o s< Il II Gf Â— D ^ o Â— o n:: 00 LU LU Io 2: < LU 00 f~\ X tC LU o cr: z < N iH 00 II Irl UJ II o Â— CL O Â—I LU C3 Â— 1 < 00 O IC_) LU 00 CNl 00 Cl. 00 00 o rÂ— I Q. 2: Â• oo o 00 /"^ 1 ^N ^ *Â— N #^>. ^ ^N rn 2 JO o Â— o o LTl LTl LO X LO O O < CD CL LTl Ln X *,
PAGE 136
122
PAGE 137
123 o UC c Ic:
PAGE 138
124 + + X X CO CO CNI CM CM CNI rH CJ O O .^ v^ '^ ^^ < < o c; >. >. + + Â•>. V rl lH iH CN + + Â— ^^<< ^^^ ^^ Q C;cr rHr< .V KÂ— Â— ++ XX Â— o IÂ— ICNJrHrH ** v^ LU *xuuu. en X oo '^''KH O^' XÂ«' ^^^^ ilnHO Ca Â—llrHI rHCNÂ«Â— '^' nnLU < U.^^rH'^ iÂ— I ULl.rHiH rICNI> + _^_rH CM IIOO >^^v_^ ^s LU'^^tÂ—l Â— 'IÂ— I V ^^'^"^^. rZnX ^^ r^CDCO OiÂ— liÂ— icsls lH rHrH'^'^ <. <. ZD Â— lHilrH LLiILLO CNI Â— Â— NN ++1JL IÂ—) LOLnuV Z>3C_)^^rH'^x CM v^^'iHrH'^'N *s^^ X rlilrH C0C02 Â— il Â— O lH iHCMn>.Q_CL Â— Â— 3 "Â— ' NSX ^^~~^UJ Â— 'rH'^'rH OO''CMnrr Â— Â— UJ 00 COOOCn^ '^Â•^O rl CSI ^N ^^^Â» * * Â— v^v^v^^ ^^ Â— ^ to iHiÂ— Iiir3 CrTQ.X'^X^^ X a. '"^^^ri!"'Â— icNic iÂ—icNi c<Â£ CO c^'^'Luz:^ Â— 'iÂ— i Â— iÂ— t c^nÂ— '^' <<00;OOUU 7Z + ^,N> iiTiÂ— ICsI>>.C')UJCOLlI ZIxrICM N^^^OO X Â— Z:On2I ^v,nn<:3 vXXlZrlCQIcai lÂ— iHXX rHCM llllrlllllhÂ»^I3 .LJiH OII<Â— II II II O Â— CO CO 3 tl I Â— II II '^'^ Â— ir^^^^LU Â— +iHil^ ZZ Â— IILUIILU^ tÂ— iiii^^dCi. Mil no Â— LI.XXII Â— '^^^cro fÂ— 1 o Â— rn^^'Â» rzi3r>.ii CO > '^^ hcT) O'HtÂ— Â— ^1Â— Â— \Â— yÂ— Â— lH tÂ— icNi~ Â— Â— ' Â— ' Â— ' v^ Â— ' i;^; IZilii di^ ^^ Â— UCO co z^ v^v^ CTtO"Â— ICNJrHCNiO Â— >ICM
PAGE 139
to O 125 < ce. UJ
PAGE 140
126 ti
PAGE 141
127 < 3 CO X ir.
PAGE 142
128
PAGE 143
129 X >o LlJ Â— i.; _i O . UCO ti* o n: < rr uj LU < < Â»> 2 u. 1 Â— N < . X 2 UJ Â— Â— CO e: s c:: Hci UJ . < UJ > Â• X s UJ q: 'r: o c; UJ CM o IÂ— 2 X < UJ rH "21 N O o ^^ bo CNl :: zr _ 3 S 2 < < < I UJ UJ "^ '^ SIClQ DUJOT ~3co CNl O^^^"^ Â— cn>c') Â— IÂ— N 2'>'>::; X Â— UJ vQ: CNl rH UJCMiIO EXCSID SLU v_.^N^ QxxIZ OÂ— C3 ^^>'~^ SiTcN ~~,<icNUj 2:cÂ£ro >Â»_CN XCn ^>.^v^q UJIÂ— IÂ— UJ ^Â— V '^zcNi csi;^::^!^ Q< r) Â— CNl 1Â— luj ^Niiuj :i:Â— i_j C^UJ^^ >Q'^iÂ— (CNIC30> '^ Â— 1 UJ Â— mv^ooujiÂ— i2^oc3co csz:^!:: ^1Â— zi3caiv^ ,_,^_^/)Lu Â— rjouj r^rHwO_ IDCicOiÂ— I i<;i^^,.N^.^N ^lÂ— OUJCDZ oca:^: :2! co Â— icnjiicniouj >rlrHCNICJ CCTOO^ hQ ^ Â— IÂ— IUJciU_ Â— OCiO Â— UJ3UJ2 cooooujQco Â— oossr^^or^uo; Â— e:uj iH CNl Â— IÂ— I CNl o ro H~ CD Â» o
PAGE 144
APPENDIX III A DESCRIPTION OF THE "GAUSS" OPTIMIZATION CODE The GAUSS computer code was written by Kylstra* at the University of Florida as a general program for finding the "least squared error" optimum for user specified problems. The user supplies a subroutine, UFUNT, which calculates the error functions for his particular problem. He also supplies the auxiliary input and output subroutines, USEGI and USEGO, which communicate any additional values to be used in UFUNT. The main program uses the methods of Gauss and Legendre for finding the optimum set of variables. Since these methods are derived in reference 25, only the results will be quoted here as they apply to the problem of selecting the optimum poison concentrations. The normalized power density vector, as calculated in UFUNT by the poison eigenvector expansion technique, is subtracted from the flat "ideal power" vector to give an error vector E. The weighted sum of the squared elements of this error vector is the scalar error value to be minimized. *C.D. Kylstra, personal communication. 130
PAGE 145
131 ERROR = E = E^ [W] E, (119) where [W] is the matrix of fractional volumes used to weight the power density. The Gaussian optimization technique uses both the scalar error and the error vector. By calling the UFUNT subroutine repeatedly, varying only one poison concentration from a given set each time, an m by n matrix is established giving the variation of each element of E with each poison concentration parameter, N(j). Here, m is the nun;ber of power densities computed, and n is the number of variable poison concentrations. This matrix is just a finite difference approximation of the Jacobian matrix, [J] , as evaluated at the given set of (25) poison concentrations. ^ ' Using this matrix, the set of parameters, N, for which the scalar error will vanish is predicted by the equation AFf = {[J] [W] [J]'^}"'[J] [W]E. (120) However, it may be that the predicted optimum, N+AN, lies beyond the region in which the assumptions made in replacing [J] with its finite difference approximation are valid. This may be determined in one of two ways. Either the magnitude of the vector AN is greater than some allowed limit, or, when the GAUSS program calls UFUNT to evaluate the scalar error for the parameters N+AN, it finds the error has increased, rather than decreased. In either case, the parameter step size is limited to some smaller value, r, and
PAGE 146
132 the program minimizes the scalar error, E, on the surface of a "hypersphere" in "N space". That is, the problem is constrained such that Ian'I = r^. (121) By writing the Lagrangian for this constrained problem, making a virtual displacement and requiring that the Lagrangian be stationary, the solution may be found as AN' = {[J] [W] [J]'^A[I] }"^[J] [W] E, (122) where A is the Lagrange multiplier of the constraint relation (121) . The value of X may be computed for a given step size from the relation r^ = e'^[J]'^{[J] [W] [J]'^A[I] }~^[J]E. (123) However, this relation is not used due to the added calculations that it would require. Instead, A is varied directly. As A ranges between infinity and zero, r varies between zero and AN, the predicted optimum step size. Thus, the choice of a large value for A restricts the search for an optimum to a small hypersphere near the point at which the [J] matrix was evaluated. The Lagrange multiplier. A, when used in this context, is known as the Levenberg parameter. Using a large initial value for A, AN' is computed from equation (122) . The main GAUSS program then calls UFUNT to evaluate the scalar error for the new parameter set.
PAGE 147
133 N+AN'. If it is indeed less than the error at N, the program continues to increase the step size, r, by decreasing A, until either the scalar error begins to increase again, or the limiting value on AN is reached. The best value of N+AN' is then designated the new optimum, about which another search iteration commences with the construction of a new Jacobian. If no point is found for which the scalar error is smaller than the beginning point for the present iteration, the program constructs a new Jacobian about the same point, using a smaller step size to compute the derivatives. This whole process continues until it is determined that the minimum point lies within some specified convergence distance of the last point for which a Jacobian was constructed. The program then outputs the optimum parameter set. This basic GAUSS program had to be modified in two ways for use with the poison eigenvalue search technique. First, the ability to use a weighting matrix, [W] , was incorporated into the program as reflected by the above equations. Also, the program was modified to retain the value of e, the poison eigenvalue, that was calculated for the point, N, about which the search for AN is taking place, so that it may be used as the starting guess at the value of e for all trial poison distributions, N+AN. The value of Â£ must also be saved for the best point found so far. Without this latter modification, the author experienced
PAGE 148
134 some difficulty with the UFUNT subroutine converging to the wrong root of e. This occurred when the main GAUSS program, returning from an extended Lagrange search to the previously com^puted optimum point, brought a value of e, from the last calculation in UFUNT that was drastically different from the value of e. associated v/ith the optimum point. The user supplied subroutines are listed in Appendix IV to complete the description of the optimization program. This modified version of GAUSS works quite well as a demonstration device. However, because this method is limited to the inherently unrealistic squared errors representation, other optimization methods should perhaps be considered if the techniques advocated in this study are to be incorporated into a standardized code to be used as an engineering tool. (The unrealistic aspects of the least squared error optimization are discussed in Chapter VII.)
PAGE 149
APPENDIX IV LISTING OF THE USER SUPPLIED SUBROUTINES FOR THE USE OF THE "GAUSS" PROGRAM TO FIND THE OPTIMUM CRITICAL POISON DISTRIBUTION Partial Dictionary to the Computer Code Variables: GAUSS Variables A(I) AV(I,J) CIP (J) ERR EPINV EV ( I , J ) F(M,N) FL(L) FUMB(J) INTER (J) ITMAX MNIP NE NIP NP NR NTP NV P(J) SFl(M) SF2 (M) SP(J) TIP (J) WTF(J) Name or Variable in Text expansion coefficients, a^. (J) poison *i(j) Z(I) elements of adjoint vectors, yj_ increments for variation of the concentrations, AN(j) sum of the squared errors l/Â£ elements of eigenvectors, [F] matrix flux vector (fast and thermal elements) errors in local pov;er densities list of interface indices maximum number of iterations for coefficients allov/ed for each poison distribution number of control points number of elements in flux vector number of independently variable poison concentrations number of flux evaluation points number of regions number of power density evaluations number of eigenvectors used in the expansion pov;er densities list of values for Z^ in each region 2 list of values for Zf in each region diagonal elements of [o^] matrix poison concentrations at control points,N(j) weighting factors for each control point in error calculations (diagonal elements of [W]) eigenvalues, e., for the eigenvectors 135
PAGE 150
136 c: uLU Jc; ^ "^^ '^ ^ ClLU S U. to __ X hlH o^ c: N >tr rr :^ cr '^ J c < Â— UJ ir. c O. Â—I CO '^ LL N < Z> lO OO ^^ X c; tH X to CM c 21 C C. un ^ Â— Q. Ll. LL Â— ^^ d. ^^ ^ O O <i X in _ z: uHin LU O . N Q rH H f''^ z: o ^ toÂ— N >> Â— xc; 3 '^ < r: Â— Lu Â—I to N Â— siÂ— Ci'^ Lxj<; IÂ— ( ^^Q'>x2; 2 c; xs ^'XoLuici N z: IÂ— o DlÂ— 1Â— 00< Nil ca Â— vt^ ujxn; Lu Â— z:o 'nqoÂ— >. ^ tÂ— in sO ^ZZ cn 2 Â— IÂ— IÂ— ^^x'^C3'^ Â— 2 X Di Â— ID caOOiHtO X X LU ^(X ir Â— r<^v^iÂ— lxC:^<: (Â— '^ Â— LI. MOOOOZrCtÂ— Â— ^^ Â— to xxDxm Â— v^M _ O Â— O '^^^S> Â— 'Â» XÂ— '^'^^'' LU 1Â— oo Â— ^ .^ <^ '^ ^^to^^ ^^^^,x^f^^Â« iHrHiÂ— Ith ^^ <^ LU hAO xxO KA^^iH.iÂ— I iÂ— IrHrHx^QCslrArArO tO ^ 20 XfH xUJrAO'^O OtOOOO OOOCSIOCr^CDOOO'^* Â— 1U C3^' IÂ— Z:OiÂ— ILntH>rH.iÂ— ICNIit iHrHrH IÂ— l<;CNICMCslCNOLr\ LiJ IÂ— 00 X K\ Â— I Q. X Q_ X Â— X :: X to xUJ n x x x ""^ xO x n x x iÂ— I iÂ— I ^ Z)r3 Â• CL Â— U. Â— CLSOtOQ xCSiHQJQ OQQrAarHQQQQ Â— LU Â— O LU hU. X 7Z TZ il Â— rH Â— LU Â— ^ Â— nj OOOOOOCOCO ID COCQÂ— 21 Zr02s2I2^ ^^ I Â— Â•Â—' Â— ^' I Â— Â— ' (Â— ^^ S^'^'^Â— 'lE'^'l Â— v^^'v^v^l f^ O 3 hOnOltiOcniOO < < < Â— Â— LU3 2:ci_r:^^r:a:zr:onc!iioiic3SQv:iioo'iÂ— "^^^Â— IÂ— ns: CO tov> tÂ— Â— CO O hOiÂ— IrH CNJIÂ— IrH O CD O O csl ha il <Â— 1 1Â— I iH O O rsi CNj CJ O CJ o o
PAGE 151
137 s r: LU LU V CI. N Z^ ZZ 'N 'N n rH ^ s Â— LU V II fH iH ^^ Z UJ Â— II It C^l X r: N 33 LiiH N'N N N 00 11 > _ ^ ^ >_ X LJ ^^ Â— > <;i Â— ^^ Â— '..^ f^ ^^ 7^ ZZ ^ Â— >> UCL Â— I v.'Uj< ooooQ_ ^ Â— Â• fJ ^^ v/ ^' ^ Â— ' UJ 00 rv\tÂ— IrH IÂ— liÂ— IrH C30 CDOO 000 llil>lltlll rHrHrl OQ xQQQ QQQ OOtHOOOLUOOOLU i^z: II zz::^zzz2rz:zzi:D LU'JJ LULlJUJ Â— UJUJLU Â— d; HHJHlIÂ— IÂ— viÂ— Kl^r) CrtCrtOdciCiOi^ Â— Â— OLUI^ i::;202ii:i;o:r:32oc:iui
PAGE 152
138 o o LlI 00 liJ o CO < Uc: o o o r. DT O D h00 UU to I < ^ o o CJ < 2 c; c UJ "" o UJ Â— =: th< 2; uc; o o a. LU Â— rj J < UJ :r = o I f>Â— CO Q Q :3 < II00 >3 00 Â— GLU Ci fO 2 3 O o n; 00 a. Â— o 12 OJ o :r Q Â— Â— LU Hh00 < r) Z5 M o Â— cr UJ s ca 20 Â— =) I00 >2. < o 00 n: Â— UJ o ir. lO cr
PAGE 153
139 X r) c cc c: I< u. X Â— ^2 JQ 2 O N nO. '^ ll Â— II 
PAGE 154
140 ZD Li. =D LU hZD O =3 Z> 00 <: X . cr o "^ c: < < o ij_ c: a: o QCi o o c: 00 UDco ~ CO ^ < ^ o c: o Â— Â— HÂ— 1 h< < O M c: 2: Â— uj r? z: :si uÂ— uho c; c o o r: c; CO hLU CO ULU < o X o hcr 2 LJ Li. Â— 1/7 O < r> zr 00 LU LiJ LU ~ =3 r: IÂ— I I< >> >CO CD LU LU IÂ— LU IIIÂ— LU O Â— IÂ— LU a:
PAGE 155
141
PAGE 156
142 a. Â•K c/) + LiÂ— C ^ Q. ^ Â— a. hu. ^ '^ ^^ hzi fÂ— ^ UJ > 2 V in N 2 IS s IÂ— I * >t V N Il II '^ II ,H rH II Â— Â— Â— II II Â— X q: ^ 00 CO o :d lH * I '> ^^ Â— ^N v_. c3 in CL ^^ Â— I.: r: rH Â— J ^x o ra nHLj_ < a. u. tc ^^ v^ >^ >^ v_* LU UJ CO + r 1 1 ci e; ti i: + I Â— en Â— CO CL ^^ _ + II II UII Â— a. ^ r: o CL O 2 Â— CL Â— O CL + CO r) 13 Â— CO Â— CO rl 25 II <; Â— II ro O 3 O CO o o o CO o o II o c; Q. LU m w CO Ci^ 2 Ic r) ^ Q_ LL V + rl II c: II ci Â— ^^ LU C^J Â«^ II IÂ— I CO o 5 rr Q U. LU CM L'~1CDOOC3C30 O >Â— I CT O O O O Â• O LH Ll LO L^ Ll O La :1OC3QQQlU Q GOOOOOloluO ^' ^*' ^*^ ^^ V_^ N.^ ^ ^ s^_^ ^ Â— LULULULULU< Â— LuS ^ttihHnii3 cii: Â— hQ LL ci c: ci: c: r; o o ci Lu z Â— r^r:^33l::Li.or:arLu ca O '3 O rH LA
PAGE 157
143 tl lU "Z. to Â— CL X <
PAGE 158
144
PAGE 159
145 c^ CNi r'' ^ LT Lc r~cc c c IÂ— . oj r^ iLT, ic r^ o; c" c r' CN N" .^ LT u. r~~ o: c" c: IÂ— c^ r^ .2^ tr IX rH iÂ—i 1Â— I .Â— I iH ci cvi csi CN c^; cv cn: cv C; c^' f^ r^ i^ r^ N~ r<~ K" ^ > > > ::; > > > A^. 1^ o CM N o CSJ in
PAGE 160
146 cri o IÂ— c : N" .^ ir N~ ^ ^ =C r.:i1^ cc cr o t.:s^ atn Ui f< ^~ fÂ— J LT LC r~C, C C. U^. LO C Ln ITi LT LT ir LT C Â•Cvr^ ^ W. IS. t^ tx C C r^; > ^ ^ ^~ t^ ^ ^ ^ ^^^^^~^_^^^^^ UJ Io IX. O UJ < > Jo > 00 CO JJ< o o CL 1^ Â— ^ < < + I II CI. II ^ 3 :3 i>i II 1 Â•TJ Â— o ^ 3 rr < >CO o o o _J O LU r: c II Â— > Â— a o < ^ c^ o Â— CM CO 1 UJ LU < < CO H^ ^ C/) O Â£0 O O < . CD LI. II Â— Q Lr> to JJC3 LTl N LO LO N . o I in < O Â— X CO Â— I a: 3 i^ + II UJ LU O =3 =3 II LU O Â— 2: i^ 1Â— ^ II Â— :r U. i^ ^o < o Q LU LH + < ::^ 1^ II I ui i:: o Â— to II 1 II o 1^ o :> o Q Â— r: Â— Q to to LA La to to N S O CNI to LO i^ + + i^ i^ I 3 I I "O II w ^ II O U. LJL 3 < + < Q O > >CO o o < LA CO r~ Â— O 3 t/) hO > LJL O CO JCD La La La O CM 10 to O UJ II r: 3 hÂ•^ o < o La to > ::z Â— I II o II La O Q i^ II 3 3 u_ Â— Â— < o LU 3 < 3 Q (_3 ::: o Â— Â— ct: CO HO* 2: Q O II >CO > UJ o < _J QLU o La o o o o o C_J
PAGE 161
147 r^r~r^r~r^r~.f^r^ccccocccocccococococcrcc'cr.
PAGE 162
LIST OF REFERENCES 1. Harker, W.H., "Designing Reactor Cores to Achieve Predetermined Power Shapes," Transactions of the American Nuclear Society; 1 , 142144 (1958) . 2. Haling, R.K., "Operating Strategy for Maintaining an Optimum Power Distribution Throughout Life," Proceedings of the Topical Meeting on Nuclear Performance of Power Reactor Cores , 205, San Francisco, American Nuclear Society (September, 1963) . 3. Crowther, R.L., "The PowerControl Method for Solution of Diffusion Equations Describing Large NonLinear Inhomogeneous Reactors," Transactions of the American Nuclear Society; 6 , 276277 (1963) . 4. Motoda, H. and Kawai, T., "A Theory of ControlRod Programming Optimization in TwoRegion Reactors," Nuclear Science and Engineering; 39 , 114118 (1970) 5. Terney, W.B. and French, H. , "Control Rod Programming Optimization Using Dynamic Programming," Nuclear Science and Engineering; 39 , 109114 (1970) . 6. Wade, D.C. and Terney, W.B., "Optimal Control of Nuclear Reactor Depletion," Nuclear Science and Engineering; 45 , 199217 (1971). 7. Jackson, W.R., "An Eigenfunction Approach to the Problem of Reactor Poison Distribution," Master's Thesis, University of Florida, Gainesville (1969) . 8. Foderaro, A. and Garabedian, H.L., "A New Method for the Solution of Group Diffusion Equations," Nuclear Science and Engineering; 8 , 4452 (1960). 9. Garabedian, H.L. and Thomas, D.H., "An Analytic Approach to Two Dim.ensional Reactor Theory," Nuclear Science and Engineering; 14 , 266271 (1962) . 10. Kara, F. and Shibata, H. , "A Study of Core Optimization from the View Point of Neutron Flux Distribution , " Journal of Nuclear Science and Technology; 5 , 323332 (July, 1968) . 148
PAGE 163
149 11. Kaplan, S., "The Property of Finality and the Analysis of Problems in Reactor SpaceTime Kinetics by Various Modal Expansions," Nuclear Scienc e and Engineering; 9 , 357361 (1961X7 ~' 12. Kaplan, S. and Yasinsky, J.B., "Natural Modes of the Xenon Problem with Flow Feedback Â— An Example," Nuclear Science and Engineering; 25 , 430438 (1966) . 13. Henry, A.F., ''The Application of Inhour Modes to the Description of Nonseparable Reactor Transients," Nuclear Science and Engineering; 20 , 338351 (1964) . 14. Ewen, R.L., "MULE~A Fortran Program for the Calculation of Three Types of Overtone Modes," WAPDTM369 (October, 1963) . 15. Sommerfeld, A., Partial Differential Equations in P hysics, Academic Press, New York (1949) . 16. Clark, M. , Jr. and Hansen, K.F., Numerical Methods of Reactor Analysis , Academic Press, Nev/ York (1964) . 17. Shilov, G.E., An Introduction to the Theory of Linear Spaces , PrenticeHall, Inc., Englewood Cliffs, N.J. (1961) . 18. White, P. A., "The Computation of Eigenvalues and Eigenvectors of a Matrix," J. Soc. In d. Aopl . Math: 6, 393437 (1958). Â— 19. Wielandt, H., "Bestimmung hoherer Eigenwerte durch gebrochene Iteration," Bericht der aerodynam ischen Versuchsanstalt Gottingen, Report 44/J/37 (1944) . 20. Hansen, K.F. and Clark, M. , Jr., "Adjoint Functions and Orthogonality Relations," Nuclear Scie nce and E ngineering; 15 , 139141 (1963)"^ 21. Wachspress, E.L., "A Numerical Technique for Solving Group Diffusion Equations," Nuclear Sci ence and Engineering; 8 , 164170 (196077 22. Varga, R.S., "Numerical Methods for Solving MultiDimensional Multigroup Diffusion Equations in Nuclear Reactor Theory," Proceedings of Symposia m Applied Mathematics, Vol. XI , G. Birkhoff, E. P. Wigner, Eds., 164189, American Mathematical Society, Providence (1961).
PAGE 164
150 23. Anderson, E.G. and Putnam, G.E.,"CORA Â— A Few Group Diffusion Theory Code for OneDimensional Reactor Analysis," IN1416 (August, 1970). 24. Meal, D.E. et aA . , "The ARC System OneDimensional Diffusion Theory Capability, DARCID," Argonne National Laboratory, Argonne, Illinois, ANL7715 (May, 1971) . 25. Wilde, D.J., and Beightler, C.S., Foundations of Optimization, PrenticeHall, Inc., Englewood Cliffs, N.J. (1967) .
PAGE 165
BIOGRAPHICAL SKETCH Steven Michael Long v;as born in Bethesda, Maryland, on May 22, 194 6. He was graduated from Walter Johnson High School in June, 1964. In June, 1968, he graduated with high honors from Princeton University, receiving his bachelor's degree in aeronautical engineering. Upon graduation, he was inducted into the Phi Beta Kappa honorary scholastic fraternity. The award of an AEC Special Fellowship provided the financial support necessary for him to enter the University of Florida Graduate School in September, 1968. In December, 1969, he received a Master of Science degree in the Department of Nuclear Engineering Sciences. In continuing his studies at the University of Florida, pursuant to the degree of Doctor of Philosophy, he has held a College of Engineering Fellowship and an assistantship with the nuclear engineering department. On June 17, 1972, he married the former 'Diana Dombey of St. Petersburg, Florida. While at the University of Florida, he has been a student member of the American Nuclear Society. 151
PAGE 166
I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. G . R . Dalton , Chairman Professor of Nuclear Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Professor of Nuclear Engineering I certify that I have read this study and that in rr.y opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, a< a dissertation for the degree of Doctor of Philosophy. N.J. Dia Assistant of Hue 1 ear Engineering I certify that I have read this study and tr.at in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. AA V\, VsvN V i.\,>Koy U.K. Kurzv/eg Associate Professor c Sciences :r Engineering and Mechanics
PAGE 167
I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. /5TtT/K^Ji*<^ D.T. Williams Professor or Aerospace Engineering This dissertation was submitted to the Dean of the College of Engineering and to the Graduate Council, and v;as accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December, 19 72 Dean, College Dean, Graduate School
PAGE 168
v^

