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

Material Information

A multigroup eigenfunction expansion technique for selecting poison distributions for optimum reactor power distribution
Long, Steven Michael, 1946-
Publication Date:
Copyright Date:
Physical Description:
xiv, 151 leaves. : ; 28 cm.


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 )
bibliography ( marcgt )
non-fiction ( marcgt )


Thesis--University of Florida.
Bibliography: leaves 148-150.
Additional Physical Form:
Also available on World Wide Web
General Note:
General Note:

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:
000869354 ( AlephBibNum )
AEG6379 ( NOTIS )
014267171 ( OCLC )


This item has the following downloads:

Full Text

A Multigroup Eigenfunction Expansion Technique
For Selecting Poison Distributions
For Optimum Reactor Power Distribution





This dissertation is dedicated to my lovely wife,


who ungrudgingly shared the effort of completing

the beast during our first months of marriage.


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.








ABSTRACT. . . . .






SET . . . . . . . . .




STUDY . . . . . . . . .


PROGRAM "CODE2" . . . . . .

* . 6

* .12

* .18

* .35

* .47

. .65

. 104

. 113

. iii

. vii

. .viii

S. xxii

. . 1


* .

: : :

. . . . . .


"CODE2" . . . . . . .. . .119

MIZATION CODE. . . . . . . .130

POISON DISTRIBUTION . . . . .... .135

LIST OF REFERENCES . . . . . . . . .. .148

BIOGRAPHICAL SKETCH. . . . . . . . . .151


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


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







D (r)




F(m,i, [N])






energy per fission (in appropriate units)

expansion coefficient associated with function

vector representing a set of expansion

matrix on left-hand side of various eigen-
value formulations of the criticality equa-

axial geometric buckling

matrix on the right-hand side of the general
eigenvalue formulation of the criticality

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

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


[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

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

[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


Cy M


a (k,j)


f (n-k,j)
Ek (r)

ZR (k, j)

Xk or X(k)


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



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


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


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.


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.




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


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

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.



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


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

*(r) (3)
AZf (r)

and that the desired distribution of control poison follows

directly from equation (1) as

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


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

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,
*(r) = ai i(r), (8)

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
After invoking the eigenvalue relationships of equation (5),

this becomes

aiyi i (r)

p (r) i=(9)


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)

The optimum real flux is constructed as in equation (8),
R(r) = aii (),

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
E (r) a (1-()

Ei (F) = (13)

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.



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 Vn n
-V'D k(r)v kk(r)+k (r)k (r)-X (I, ),n(r)
-V*D () ( +Ea f
k- k n- k- k
+ (r)E k (r)- E(n-k,r)n ( )+C(r) (r) = 0. (1

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
-V*D k(r)Vk + (r)+E (r) ()- f) n X"y ( )
k k n( k k
+ER(r) (r)- Es(k-n,r)) (r)+Z (r)k) (r) = 0 (15)
R p

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
P(r) = A E(r) (r). (16)

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

k(F) = K (18)

A2 n s rn (F)

(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.
These solutions would have the form
n (r)
S (r) P* (r) S (r) z (r)
k n=
Z (r) = VD (r)V K
n -n k
S n(V) (r)n rS k (r)P*(r)
n= 1
+xk S n(r),,n()/Sk (r)-f ka(r)-ER)

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).

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)

k n"v n k k k
+X k sn (F) vE (r)/S (r) -a R Ek (r)

+ s (nk,r) Sn (r)/S (r). (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 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
(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.



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)

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)

Esz(n 3k,j) (n,j)+Zp(k,j) (k,j)
-X(k) I vZf(n,j) (n,j)+ZR(k,j) (k,j) = 0, (26)

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


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, i, will have (J.K)


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


(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
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


(l/a-1l/a.)w.[B]j. = 0. (36)

Substituting .i for w., this gives the orthogonality condi-


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)

-X(k) Z (nJ)(n,j) (n,j)+ER(k,j)(i(k,j)
Es(n-k,j) i(n,j) = Yi4i(k,j), (39)

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)

-1 1s (n-k,j)i(n,j) +R(k,j)i(k,j)

= i X(k) IVf (n,j)f(n,j), (42)

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)

-X(k) vE f(n,j) (n,j)+ER(k,j)( (k,j)
-1 Es(n-k,j){i(n,j) = -8io (k,j) i(k,j), (44)

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)




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-


ai = i[Op] (47)

The flux at the control points, at least, is represented

= aii (48)
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

[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-

-[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

[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


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-


[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.



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

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


{[L]- K [S]i = pi [Q]i (59)

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)

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

V = a i (61)

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).

[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)
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'


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,


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..

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.



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)


s(n-kj)0 i(n,j)

SK1 -X(k) vzE(n,j)(n,j), (73)
eff n

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-


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)

- -X(k) vZf (n,j)4 (n,j)+E (k,j)( (k,j)
eff 7 f

Es(n k,j). (n,j) = -cE (k,j) (k,j). (74)

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..

S= aii' (76)

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.

[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

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

-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 .

-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

am m = E aiF(m,i,[N]), for m = 1 to M, (82)


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


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,
'1, gives a scalar equation which may be solved for E to


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-


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

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


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

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)

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
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.



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

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.



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

S-.---A----A ---._..._, .

0 Fast Flux
SThermal Flux

0 1 2 3
0 10 20 30 40

X (cm)



0 10 20 30 40
X (cm)

0 0

0 10 20 30 40

X (cm)



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


0 0
o A

o Fast Flux

& Thermal Flux

S 10 20 30 40
0 I0 20 30 40




X (cm)

I 1 II
) 10 20 30 4C


0 /


0 10 20 30 40

X (cm)


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


The Eigenvalues Associated with the Poison
Eigenvectors for Example 1

























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-


+ + +
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


Optimization Results for Example 1

Results of General Poison Distribution Search Method:

2.75802E 01
2.75777E 01
2.75518E 01
2.73029E 01
2.49391E 01

2.44744E 00
2.44756E 00
2.44878E 00
2.46055E 00
2.57227E 00


Results of Power Vector Expansion Coefficient
Real Solution:

2.75802E 01
2.75778E 01
2.75518E 01
2.73029E 01
2.49391E 01

2.47444E 00
2.47555E 00
2.48788E 00
2.40555E 00
2.52266E 00



Search Method:


2 .7200 1E- 1
-8.81539E 01
4.88696E 01
-2.71797E 01
8.85565E 01



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

4.48286E 02
4.21704E 02
4.33649E 02
4.28244E 02
4.31338E 02



-1.16233E 00
-1.14970E 00
-1.17165E 00
-1.12248E 00
-1.25350E 00

-1.54869E 03
-9.86252E 02
-1.22898E 03
-1.13463E 03
-1.19011E 03






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,


= ia+ +a ~ i (102)

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

P*= b.P. (104)

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


(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
degree that the products m [p]~i, for m ? i, have values
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


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

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



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


Three-Region Poison Optimization Results for Example 2
























































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)


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


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


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)

Region II:

Group (k)

Group (k)

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)



pure water
D(k) E (k)
1.136999 .0005892
0.149380 .0019240




Zf (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.


i 'I

Region I Region II

0 50 140 190
I i
1/ I


Radius (cm)

I I/
Region I I Region II I Reflector

I A Thermal Flux


0-A-he al-Fl I
0 50 140 190
Radius (cm)

Figure 6 Poison Eigenvectors +1 and $i for Example 3