Citation
A multigroup eigenfunction expansion technique for selecting poison distributions for optimum reactor power distribution

Material Information

Title:
A multigroup eigenfunction expansion technique for selecting poison distributions for optimum reactor power distribution
Creator:
Long, Steven Michael, 1946-
Publication Date:
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 )
non-fiction ( marcgt )

Notes

Thesis:
Thesis--University of Florida.
Bibliography:
Bibliography: leaves 148-150.
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 non-profit 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:

PDF ( 5 MBs ) ( .pdf )

multigroupeigenf00longrich_Page_051.txt

multigroupeigenf00longrich_Page_127.txt

multigroupeigenf00longrich_Page_038.txt

multigroupeigenf00longrich_Page_164.txt

multigroupeigenf00longrich_Page_100.txt

multigroupeigenf00longrich_Page_130.txt

multigroupeigenf00longrich_Page_098.txt

multigroupeigenf00longrich_Page_124.txt

multigroupeigenf00longrich_Page_143.txt

multigroupeigenf00longrich_Page_165.txt

multigroupeigenf00longrich_Page_018.txt

multigroupeigenf00longrich_Page_141.txt

multigroupeigenf00longrich_Page_069.txt

multigroupeigenf00longrich_Page_001.txt

multigroupeigenf00longrich_Page_059.txt

multigroupeigenf00longrich_Page_134.txt

multigroupeigenf00longrich_Page_160.txt

multigroupeigenf00longrich_Page_103.txt

multigroupeigenf00longrich_Page_104.txt

multigroupeigenf00longrich_Page_092.txt

multigroupeigenf00longrich_Page_015.txt

multigroupeigenf00longrich_Page_003.txt

multigroupeigenf00longrich_Page_123.txt

multigroupeigenf00longrich_Page_131.txt

multigroupeigenf00longrich_Page_159.txt

multigroupeigenf00longrich_Page_085.txt

multigroupeigenf00longrich_Page_139.txt

multigroupeigenf00longrich_Page_067.txt

multigroupeigenf00longrich_Page_120.txt

multigroupeigenf00longrich_Page_148.txt

multigroupeigenf00longrich_Page_128.txt

multigroupeigenf00longrich_Page_030.txt

multigroupeigenf00longrich_Page_021.txt

multigroupeigenf00longrich_Page_106.txt

multigroupeigenf00longrich_Page_113.txt

multigroupeigenf00longrich_Page_026.txt

multigroupeigenf00longrich_Page_062.txt

multigroupeigenf00longrich_Page_119.txt

multigroupeigenf00longrich_Page_076.txt

multigroupeigenf00longrich_Page_089.txt

multigroupeigenf00longrich_Page_129.txt

multigroupeigenf00longrich_Page_075.txt

multigroupeigenf00longrich_Page_079.txt

multigroupeigenf00longrich_Page_121.txt

multigroupeigenf00longrich_Page_081.txt

multigroupeigenf00longrich_Page_138.txt

multigroupeigenf00longrich_Page_033.txt

multigroupeigenf00longrich_Page_077.txt

multigroupeigenf00longrich_Page_105.txt

multigroupeigenf00longrich_Page_108.txt

multigroupeigenf00longrich_Page_065.txt

multigroupeigenf00longrich_Page_097.txt

multigroupeigenf00longrich_Page_007.txt

multigroupeigenf00longrich_Page_024.txt

multigroupeigenf00longrich_Page_151.txt

multigroupeigenf00longrich_Page_083.txt

multigroupeigenf00longrich_Page_090.txt

multigroupeigenf00longrich_Page_043.txt

multigroupeigenf00longrich_Page_099.txt

multigroupeigenf00longrich_Page_110.txt

multigroupeigenf00longrich_Page_020.txt

multigroupeigenf00longrich_Page_055.txt

multigroupeigenf00longrich_Page_155.txt

multigroupeigenf00longrich_Page_066.txt

multigroupeigenf00longrich_Page_166.txt

multigroupeigenf00longrich_Page_154.txt

multigroupeigenf00longrich_Page_029.txt

multigroupeigenf00longrich_Page_144.txt

multigroupeigenf00longrich_Page_061.txt

multigroupeigenf00longrich_Page_142.txt

multigroupeigenf00longrich_Page_161.txt

multigroupeigenf00longrich_Page_052.txt

multigroupeigenf00longrich_Page_094.txt

multigroupeigenf00longrich_Page_107.txt

multigroupeigenf00longrich_Page_071.txt

multigroupeigenf00longrich_Page_027.txt

multigroupeigenf00longrich_pdf.txt

multigroupeigenf00longrich_Page_005.txt

multigroupeigenf00longrich_Page_022.txt

multigroupeigenf00longrich_Page_028.txt

multigroupeigenf00longrich_Page_163.txt

multigroupeigenf00longrich_Page_137.txt

multigroupeigenf00longrich_Page_084.txt

multigroupeigenf00longrich_Page_050.txt

multigroupeigenf00longrich_Page_112.txt

multigroupeigenf00longrich_Page_167.txt

multigroupeigenf00longrich_Page_046.txt

multigroupeigenf00longrich_Page_032.txt

multigroupeigenf00longrich_Page_040.txt

multigroupeigenf00longrich_Page_013.txt

multigroupeigenf00longrich_Page_014.txt

multigroupeigenf00longrich_Page_070.txt

multigroupeigenf00longrich_Page_082.txt

multigroupeigenf00longrich_Page_126.txt

multigroupeigenf00longrich_Page_162.txt

multigroupeigenf00longrich_Page_042.txt

multigroupeigenf00longrich_Page_150.txt

multigroupeigenf00longrich_Page_056.txt

multigroupeigenf00longrich_Page_063.txt

multigroupeigenf00longrich_Page_016.txt

multigroupeigenf00longrich_Page_096.txt

multigroupeigenf00longrich_Page_019.txt

multigroupeigenf00longrich_Page_091.txt

multigroupeigenf00longrich_Page_017.txt

multigroupeigenf00longrich_Page_073.txt

multigroupeigenf00longrich_Page_101.txt

multigroupeigenf00longrich_Page_111.txt

multigroupeigenf00longrich_Page_140.txt

multigroupeigenf00longrich_Page_088.txt

multigroupeigenf00longrich_Page_023.txt

multigroupeigenf00longrich_Page_034.txt

multigroupeigenf00longrich_Page_010.txt

multigroupeigenf00longrich_Page_086.txt

multigroupeigenf00longrich_Page_025.txt

multigroupeigenf00longrich_Page_116.txt

multigroupeigenf00longrich_Page_156.txt

multigroupeigenf00longrich_Page_054.txt

multigroupeigenf00longrich_Page_064.txt

multigroupeigenf00longrich_Page_153.txt

multigroupeigenf00longrich_Page_149.txt

multigroupeigenf00longrich_Page_152.txt

multigroupeigenf00longrich_Page_035.txt

multigroupeigenf00longrich_Page_053.txt

multigroupeigenf00longrich_Page_133.txt

multigroupeigenf00longrich_Page_078.txt

multigroupeigenf00longrich_Page_147.txt

multigroupeigenf00longrich_Page_072.txt

multigroupeigenf00longrich_Page_118.txt

multigroupeigenf00longrich_Page_031.txt

multigroupeigenf00longrich_Page_087.txt

multigroupeigenf00longrich_Page_012.txt

multigroupeigenf00longrich_Page_047.txt

multigroupeigenf00longrich_Page_157.txt

multigroupeigenf00longrich_Page_009.txt

multigroupeigenf00longrich_Page_115.txt

multigroupeigenf00longrich_Page_109.txt

multigroupeigenf00longrich_Page_068.txt

multigroupeigenf00longrich_Page_080.txt

multigroupeigenf00longrich_Page_146.txt

multigroupeigenf00longrich_Page_122.txt

multigroupeigenf00longrich_Page_057.txt

multigroupeigenf00longrich_Page_006.txt

multigroupeigenf00longrich_Page_041.txt

multigroupeigenf00longrich_Page_049.txt

multigroupeigenf00longrich_Page_125.txt

multigroupeigenf00longrich_Page_145.txt

multigroupeigenf00longrich_Page_048.txt

multigroupeigenf00longrich_Page_117.txt

multigroupeigenf00longrich_Page_135.txt

multigroupeigenf00longrich_Page_058.txt

multigroupeigenf00longrich_Page_060.txt

multigroupeigenf00longrich_Page_093.txt

multigroupeigenf00longrich_Page_008.txt

multigroupeigenf00longrich_Page_136.txt

multigroupeigenf00longrich_Page_036.txt

multigroupeigenf00longrich_Page_132.txt

multigroupeigenf00longrich_Page_158.txt

multigroupeigenf00longrich_Page_037.txt

multigroupeigenf00longrich_Page_004.txt

multigroupeigenf00longrich_Page_102.txt

multigroupeigenf00longrich_Page_039.txt

multigroupeigenf00longrich_Page_095.txt

multigroupeigenf00longrich_Page_045.txt

multigroupeigenf00longrich_Page_114.txt

multigroupeigenf00longrich_Page_002.txt

multigroupeigenf00longrich_Page_074.txt

multigroupeigenf00longrich_Page_011.txt

multigroupeigenf00longrich_Page_044.txt


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 ONE-GROUP
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 Three-Region Poison Optimization Results for
Example 2 81

6 Reactor Model for Example 3 85

7 Comparison of Reactor Configurations Resulting
From the Two-Region 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 left-hand side of various eigen-
value formulations of the criticality equa-
tion

axial geometric buckling

matrix on the right-hand 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 (n-k,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 one-group theory. Then, the multigroup problem is out-

lined and the arbitrariness of the flux-power 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 time-space 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 two-group, one-dimensional 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 one-group, two-mode formulism to

represent a two-region, 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 two-group, two-dimensional reactor with two

gangs of control rods. Wade and Terney(6) applied optimal

control theory to a generalized set of design objectives,

using a one-group spatially nodalized reactor model. Jack-

son(7) used an eigenfunction expansion technique, based on

one-group 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 one-group 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 well-known 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 space-time 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

two-group 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,

time-space 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 flux-power shapes for several example reactors.












CHAPTER II


THE ONE-GROUP EIGENFUNCTION EXPANSION PROBLEM


Although a treatment of the one-group problem would

seem to be the simplest case of a general multigroup treat-

ment, such generalization is not necessarily possible.

Because the one-group 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 distribution-power dis-

tribution problem. Still, it is instructive to review the

one-group treatment for purposes of later comparison and

for the purpose of introducing some concepts in a less com-

plicated context. The one-group problem and solution was

described by Jackson in reference 7, and his approach will

generally be followed in this chapter.

In the one-group 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 power-flux 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,








V-D(r)V (r)+vEf (r) i(r) -ZEa(r)(i(r) = yi i(r). (5)


These eigenfunctions, the iq, form a complete, self-adjoint

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 in-core 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


V-D(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 one-group 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 one-group 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 two-group 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 to-fast 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(n-k,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 self-adjoint. 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(k-n,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 power-flux 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 one-group 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 (K-1). (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(n-k,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, pre-selected 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 (K-l), 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 one-group 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 one-group 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-

tinuous-functions 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 one-dimensional, 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,j-l)]/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,j-1) (kj-l)+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.(K-1)

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 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 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] j-w [B] [A] ([B]i = 0. (35)


With the scalar a's removed from the inner products, this

becomes

-T
(l/a-1l/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,j-l)i (k, j-1)+ 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(n-k,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 H-K

flux values need be specified. It would be convenient if

the rest could be neglected. However, these eigenvectors,

with J-K elements, satisfy the orthogonality conditions

with their adjoints 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 one-group 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,j-1) i(k,j-l)+Ea(k,j) i(k,j)

K
-1 1s (n-k,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 one-group 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,j-1) i(k,j-1)+Ea(k,j ) (k,j)

K
-X(k) vE f(n,j) (n,j)+ER(k,j)( (k,j)
n=l
K
-1 Es(n-k,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 (K-1) 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 two-group 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 two-group 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 two-group 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 two-group

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 programmer-operator 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 two-group representation of a bare, homogeneous cylindri-

cal reactor of infinite axial extent, the equations for the

two-group 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 two-group

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 (i-1) 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, one-dimensional 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 one-dimensional, few-group 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,j-1) (k,j-1)


+Ea (k, ) ( ,(k, (kj)())(k,j)+E (k,j) (k,j) (kj) (k,j)

K

s(n-kj)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,j-l)( (k,j-l)+ 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 = G-k, 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


{[BI-e [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 G-K 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 = -(i-8s) 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 (82-B)/( 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 = {l-P(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

half-slab reactor was represented by a five-point, two-energy

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 half-slab, 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-.---A----A ---._..._, .








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.55390E-02
6.54151E-02
6.41074E-02
5.16011E-02
-1.02413E-01

Search Method:

Poison
Concentration
6.55372E-02
6.54153E-02
6.41076E-02
5.16005E-02
-1.02412E-02


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.82306E-02
1.49576E-03
-1.66447E-01
2.21372E-01
-9.48818E-02


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(di-l)/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 GAUSS-UFUNT








eigenvector expanded poison optimization routine. This

example employs a two-group, ten-point 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 10-5 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 cm-1 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


Three-Region Poison Optimization Results for Example 2


Poison
Region


Point
No.


Poison
Concentration

4.66710E-05

4.66710E-05

4.66710E-05

4.66710E-05

9.68959E-05

9.68959E-05

9.68959E-05

9.68959E-05

7.47632E-06

7.47632E-06

7.47632E-06

7.47632E-06


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.77008E-03

2.21607E-02

4.43213E-02

3.04708E-02

3.60111E-02

8.86427E-02

5.30931E-02

5.77099E-02

1.32964E-01

1.55126E-01

1.77285E-01

1.99446E-01


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 Three-Region 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 two-group, nineteen-point 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 cm-1 per gm of B /cc
a (2) = 130.693 cm-1 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

/


0-A-he al-Fl 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 ONE-GROUP 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. NUr-lERICAL 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 Three-Region Poison Optimization Results for Example 2 81 6 Reactor Model for Example 3 85 7 Comparison of Reactor Configurations Resulting From the Two-Region 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 left-hand side of various eigenvalue formulations of the criticality equation 2 B axial geometric buckling [B] matrix on the right-hand 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 one-group theory. Then, the multigroup problem is outlined and the arbitrariness of the flux-power 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 time-space 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 two-group, one-dimensional 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 one-group, two-mode formulism to represent a tv/o-region, 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;o-group, two-dimensional reactor with two gangs of control rods. Wade and Terney applied optimal control theory to a generalized set of design objectives, using a one-group spatially nodalized reactor model. Jack(7) son used an eigenf unction expansion technique, based on one-group 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 one-group 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 well-known 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 space-time 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 two-group 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, time-space 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 flux-power shapes for several example reactors.

PAGE 20

CHAPTER II THE ONE-GROUP EIGENFUNCTION EXPANSION PROBLEM Although a treatment of the one-group problem would seem to be the simplest case of a general multigroup treatment, such generalization is not necessarily possible. Because the one-group 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 distribution-power distribution problem. Still, it is instructive to review the one-group treatment for purposes of later comparison and for the purpose of introducing some concepts in a less complicated context. The one-group problem and solution was described by Jackson in reference 7, and his approach will generally be followed in this chapter. In the one-group 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 power-flux 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 in-core 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 V-D(?)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 one-group 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 one-group 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 two-group 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 -v-D^(?)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/er-flux 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 one-group 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 (K-1) . (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, pre-selected 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 (K-1) , 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'^ = V-D^(r)V S^(r)P*(r) a|s^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 one-group 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 one-group 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 continuous-functions 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 one-dimensional, 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,j-l)]/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,j-l),^(k,j-l)+E^ (k,j)4)(k,j) K -2 l^{n-yk,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'(K-l) 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,j-l)(^^(k, j-l)+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 J-K 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 H-K 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 one-group 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, j-1) ii)^(k, j-l) + 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 one-group 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, j-l)(j,^(k, j-l)+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 -e-eCN (N -e-eo 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 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]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'(K-l) 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 two-group 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 two-group 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 two-group 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 two-group 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 programmer-operator 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 two-gi*oup representation of a bare, homogeneous cylindrical reactor of infinite axial extent, the equations for the two-group 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]2-v^(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 two-group 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 (i-1) 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, one-dimensional 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 one-dimensional, few-group 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,j-l)cf)(k, j-1) + Z^ (k , j ) (^ (k , j ) + Ej^ (k , j ) 4) (k , j ) +e Zp (k , j ) <|, (k , j ) K ^ ZgCn-vk, 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, j-l)(})(k,j-l)+E^ (k,j)^(k, j) K — — y Ci<\ K — X(k) ^^'^f ('^':3)'? (n,j)+Zj^(k,j)(},(k,j) eff , n=l Z^(n-vk,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 B-j , 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 (B2-B3)/ (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 = / {l-P(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 five-point, two-energy 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 half-slab, 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.5S350E-02 6.54151E-02 6.41074E-02 5.16011E-02 -1.02413E-01 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.55372E-02 6.54153E-02 6.41076E-02 5.16005E-02 -1.02412E-02 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.49576E-03 3-1.65447E-01 42.21372E-01 5-9.48818E-02 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 two-group, ten-point 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, B-i, 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 cm-1 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 Three-Region 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 two-group, nineteen-point 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 ./ / ^ — A-A-A ^ ._A— A— A— A-A— A-A 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 / / _o-0-®' -© 6-0 50 140 190 Radius (cm) Figure 6 Poison Eigenvectors 4), and (J), for Example 3

PAGE 101

87 + ^1 Region I A/ % — o-O ^— A— A-A, ,/ / A / Region II l\\ Reflector / 50 140 <>—^ 190 Radius (cm) A— A— A— A—, A— A-A— ^__ A^ Region I Region II / ® Fast Flux A Thermal Flux ,0 / y i_o-o ^o-o 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 two-region 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 two-region 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 Two-Region 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.1740E-05 -1.1680E-05 -1.2306E-05 -1.2310E-05 II 6.1978E-05 6.1956E-05 6.1935E-05 6.1936E-05 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 two-region 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 "two-group 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 flux-power 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 two-group 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 two-group, forty-four-point 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.380E-05 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» o-crrofNuocTiooin 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 r-CTt 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 one-dimensional 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, two-group 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 well-documented 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, well-documented 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 two-group, one-dimensional 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,j-l)(})(l,j-l)-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,j-1) = 7TD^[r(j-l)+r(j)]/Cr(j)-r(j-l)], C(k,j,j+1) = TTDjj[r(j+l)+r (j)]/[r (j+l)-r(j)] , C(k,j,j) = C(k,j,j-1)+C(k,j,j+1) +Vj(j) (Z^+B^d'^-X^vzJ^)j +v^^(j)(z^bV-x^zJ)^^, Vj(j) = TT[3r(j)^-2r(j)r(j-l)-r(j-l)^]/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 one-dimensional 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,j-1) C(2,j,j-1) 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] [Zj-1]}~^[YJ], [Xj] = [Yj]{[Zj]}"^ = [Uj]-[Wj] [Zj-1] (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]Tj-l. (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 ] [Zj-1] } when finding the [Zj] matrices from equation (113) . In general, for a K-energy-group 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 ^'^ • i-H O C i-l 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< I-l II Gf — D --^ o — o n:: 00 LU LU Io 2: < LU 00 f~\ X tC LU o cr: z < N i-H 00 II Ir-l 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 C-J O O .^ v^ '-^ ^^ < < o c; >. >. + + •>. V r-l l-H i-H CN + + — ^-^<< ^^--^ ^^ Q C;cr rHr-< .-V K— — ++ XX — o I— ICNJrHrH ** v-^ LU *-xuuu. en X oo '-^'-'K-H O^-' X«' ^^^^ i-lnHO Ca —llrHI rHCN«— '^-' nnLU < U.^^rH'^ i— I U-Ll.rHi-H r-ICNI> + _^_rH CM IIOO >^^v_^ ^-s LU'^^t—l — 'I— I V ^^'^"^^-. rZnX ^^ r^CDCO Oi— li— icsls l-H rHrH'-^'-^ <. <. ZD — l-Hi-lrH LLi-ILLO CNI — — NN ++1JL I—) LOLnu-V Z>3C_)^^rH'-^x CM v^^-'i-HrH'-^'-N *-s^-^ X r-li-lrH C0C02 — i-l — O l-H i-HCMn>.Q_CL — — 3 "— ' NSX ^^~~^UJ — 'r-H'^-'rH OO'-'CMnrr — — UJ 00 COOOCn^ '-^•-^O r-l CSI ^N ^^^-» * * — v^v^v^^ ^^ — -^ to i-Hi— Ii-ir3 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 ZIxr-ICM N^^^OO X — Z:On2I ^-v,-n--n<:3 vXXlZr-lCQIcai l— i-HXX r-HCM llllr-l|llllh-»^I3 .LJi-H OII<|— II II II O — CO CO 3 t-l I — II II '-^'^ — ir^^^^LU — +i-Hi-l^ ZZ — IILUIILU^ t— iiii-^-^d-Ci. Mil no — LI.XXII — '^^-^cro f— 1 o — rn^^'-» rzi3r>.ii CO > '^-^ h-cT) O'-Ht— — ^1— — \— y— — l-H t— icNi~ — — ' — ' — ' v^ — ' i;^; IZilii di^ ^-^ — UCO co z^ v-^v^ CTtO"— ICNJrHCNiO — >-ICM
PAGE 139

to O 125 < ce. UJ

PAGE 140

126 t-i

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 -D-UJOT ~3co CNl O-^^^"^ — cn>c') — I— N 2'->'->::; X — UJ vQ: CNl rH UJCMi-IO 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 — icnji-icniouj >r-lrHCNIC-J CCTOO^ h-Q ^ — I— IUJciU_ — OCiO — UJ-3UJ2 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 hl-H o^ c: N >tr rr :^ cr '^ -J c < — UJ ir. c O. —I CO '-^ LL N < Z> lO OO ^^ X c; t-H 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 D--l— 1— 00< Nil ca — vt^ ujxn; Lu — z:o '-nqo— >. ^ t— in sO ^ZZ cn 2 — I— I— ^^x'-^C3'^ — 2 X Di — ID caOOi-HtO X X LU ^(X ir — r<^v^i— lxC--:^<: (— '-^ — LI. MOOO-OZrC-t— — ^-^ — to xxD-xm — v^M _ O — O '^^^S> — '-» X— '-^'-^---^'-' LU 1— oo — ^ .--^ <-^ '-^ ^^to^^ ^-^^-^,-x^f^^-« i-HrHi— It-h ^^ <^ LU hAO xxO KA^^i-H.i— I i— IrHrHx^QCslrArArO tO ^ 20 Xf-H xUJrAO'^O OtOOOO OOOCSIOCr^CDOOO'-^* — 1-U C3^-' I— Z:Oi— ILnt-H>r-H.i— ICNIi-t i-Hr-HrH 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 xCSi-HQJ-Q OQQrAarHQQQQ — LU — O LU h-U. X 7Z TZ i-l — r-H — LU — -^ — nj OO-OOOOCOCO ID COCQ— 21 Zr02s2I2^ ^-^ I — •—' — ^-' I — — ' (— ^-^ S^-'^-'^— 'lE'^-'l — v^^-'v^v^l f^ O 3 hOnOltiOcniOO < < < — — LU-3 2:ci_r:^^r:a-:zr:onc!i-ioiic3SQv-:iioo--'i— "^-^-^— I— ns: CO tov> t— — CO O hOi— IrH CNJI— Ir-H O CD O O csl ha i-l <— 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 f-H i-H ^^ Z UJ — II It C^-l X r: N -3-3 Lii-H N'-N N N 00 11 > _ ^ ^ >_ X LJ ^-^ — > <;i — ^^ — '-..^ f-^ ^^ 7^ ZZ ^ — >> U-CL — I v.'Uj< ooooQ_ ^ — • fJ ^^ v-/ ^-' ^ — ' UJ 00 rv-\t— Ir-H I— li— IrH C30 CDOO 000 l-li-l>l-lt-ll-l rHr-Hr-l OQ xQQQ QQQ OOt-HOOOLUOOOLU i^z: II zz::^zzz2rz:zzi:D LU'JJ LULlJUJ — UJUJLU — d; H-H-J-Hl-I— I— v-i— Kl^r) CrtCrtOdciCiOi^ — — OLUI^ i::;202ii:i;o:r:32oc:iu-i

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. '-^ l-l — 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 I-l II '-^ II ,-H rH II — — — II II — X q: ^ 00 CO o :d l-H * 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; t-i i: + I — en — CO CL ^^ _ + II II UII — a. ^ r: o CL O 2 — CL — O CL + CO r) 13 — CO — CO r-l 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 + r-l 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 L-l LO L-^ L-l O La -:1-OC3QQQlU Q GOOOOOloluO ^' ^*-' ^*-^ ^^ V_^ N.^ ^ ^ s^_^ -^ — LULULULULU< — LuS ^t-t-i-h-H-ni-i-3 cii: — h-Q 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 t-l 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 r-H i—i 1— I .— I i-H c-i 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 u-i 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 , 142-144 (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 Power-Control Method for Solution of Diffusion Equations Describing Large NonLinear Inhomogeneous Reactors," Transactions of the American Nuclear Society; 6 , 276-277 (1963) . 4. Motoda, H. and Kawai, T., "A Theory of Control-Rod Programming Optimization in Two-Region Reactors," Nuclear Science and Engineering; 39 , 114-118 (1970) 5. Terney, W.B. and French, H. , "Control Rod Programming Optimization Using Dynamic Programming," Nuclear Science and Engineering; 39 , 109-114 (1970) . 6. Wade, D.C. and Terney, W.B., "Optimal Control of Nuclear Reactor Depletion," Nuclear Science and Engineering; 45 , 199-217 (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 , 44-52 (1960). 9. Garabedian, H.L. and Thomas, D.H., "An Analytic Approach to Two Dim.ensional Reactor Theory," Nuclear Science and Engineering; 14 , 266-271 (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 , 323-332 (July, 1968) . 148

PAGE 163

149 11. Kaplan, S., "The Property of Finality and the Analysis of Problems in Reactor Space-Time Kinetics by Various Modal Expansions," Nuclear Scienc e and Engineering; 9 , 357-361 (1961X7 ~' 12. Kaplan, S. and Yasinsky, J.B., "Natural Modes of the Xenon Problem with Flow Feedback — An Example," Nuclear Science and Engineering; 25 , 430-438 (1966) . 13. Henry, A.F., ''The Application of Inhour Modes to the Description of Nonseparable Reactor Transients," Nuclear Science and Engineering; 20 , 338-351 (1964) . 14. Ewen, R.L., "MULE~A Fortran Program for the Calculation of Three Types of Overtone Modes," WAPD-TM369 (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 , Prentice-Hall, 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, 393-437 (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 , 139-141 (1963)"^ 21. Wachspress, E.L., "A Numerical Technique for Solving Group Diffusion Equations," Nuclear Sci ence and Engineering; 8 , 164-170 (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., 164-189, American Mathematical Society, Providence (1961).

PAGE 164

150 23. Anderson, E.G. and Putnam, G.E.,"CORA — A Few Group Diffusion Theory Code for One-Dimensional Reactor Analysis," IN-1416 (August, 1970). 24. Meal, D.E. et aA . , "The ARC System One-Dimensional Diffusion Theory Capability, DARCID," Argonne National Laboratory, Argonne, Illinois, ANL-7715 (May, 1971) . 25. Wilde, D.J., and Beightler, C.S., Foundations of Optimization, Prentice-Hall, 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^