THE APPLICATION OF VOLTERRA SERIES
AND NONLINEAR OPERATORS TO NUCLEAR
REACTOR KINETICS
By
IRA THIERER
A DISSERTATION PRESENTED TO THE GRADUATE COUNCIL OF
THE UNIVERSITY OF FLORIDA
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE
DEGREE OF DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
April, 1967
ACKNOWLEDGMENTS
The author wishes to express his sincere appreciation to the
members of his supervisory committee and to Dr. R. E. Uhrig for
allowing him the freedom to undertake the particular line of research
presented in this dissertation. He is grateful to Drs. R. B. Perez,
M. J. Ohanian, and A. P. Sage for the many suggestions they have
offered and the continued counsel they have given throughout the
work.
The author also wishes to extend thanks to the Atomic Energy
Commission and the National Aeronautics and Space Administration for
their support during the period of research.
Appreciation is due to the staff of the University of Florida
Computing Center for their invaluable assistance. Finally, a very
special thanks goes to Mrs. Philamena V. Pearl who undertook the
preparation of the final manuscript.
TABLE OF CONTENTS
Page
ACKNOWLEDGMENTS . . . .
LIST OF TABLES . . . .
LIST OF FIGURES . . .
ABSTRACT . . . . .
Chapter
I.
II.
III.
IV.
V.
VI.
VII.
VIII.
APPENDIX
A.
B.
C.
D.
E.
F.
G.
S . . . .. . . .*
INTRODUCTION . . . . . .
VOLTERRA SERIES . . . . .
MAXIMUM VALUE SERIES . . . .
NONLINEAR OPERATORS THEORY . .
NONLINEAR OPERATORS APPLICATION
EXAMPLES AND RESULTS . . . .
SUMMARY AND CONCLUSIONS . . .
RECOMMENDATIONS . . . . .
. . . . . .
. . . . . .
REGION OF APPLICABILITY OF ZEROPOWER KINETICS . .
CONSTANTS FOR EBWR FEEDBACK. . . . . . . .
A TWOTEMPERATURE FEEDBACK MODEL . . . . . .
HYPERGEOMETRIC FUNCTIONS OF TWO VARIABLES . . .
MULTIDIMENSIONAL TRANSFORMS OF KERNEL
COMBINATIONS . . . . . . . . .
ASSOCIATION FORMULAS . . . . . . . . .
COMPARISON OF OPERATOR METHOD CALCULATION AND
"TRUE" RESPONSE . . . . . . . . .
iii
v
vi
viii
1
7
29
40
'53
70
84
88
91
93
94
97
100
102
103
. .
. . .
LIST OF REFERENCES . . . . . . . . . . 104
BIOGRAPHICAL SKETCH. . . . . . . . . . 106
LIST OF TABLES
Table Page
1. Constants for the ZeroPower Kernel . . . . . 20
2. TwoTemperature Feedback Model Data . . . . . 24
LIST OF FIGURES
Figure
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
Feedback Model for Boiling Water Reactors
Feedback Function Based on EBWR Data ..
First Approximation in the Volterra Series
First Iterative Solution EBWR Problem
Linear Response EBWR . . ..
First Iterative Solution TwoTemperature
Feedback Problem . . . . ....
Linear Response TwoTemperature Feedback
Region of Convergence of the Maximum Value
Regions of Absolute Stability EBWR Model
General Nonlinear System . . ...
11. Combination of Operators .
* .* . .
* .
. .
* .
Model
Series
0 0@OQI
* 0 *
Reduction of a General Reactor System Part I
Reduction of a General Reactor System Part II
* *
14. A Reactor System without Delayed Neutrons .
15. "True" Response of System with No Delayed Neutrons 
Exponential Input . . . . . ..
16. "True" Response of System with No Delayed Neutrons 
Step Function Input . . . . ...
17. Operator Method Response of System with No Delayed
Neutrons Step Function Input . . . ....
18. Operator Method Response of TwoTemperature
Feedback Model Step Function Input .... ...
Page
16
18
19
21
22
25
26
34
36
42
45
54
57
73
. . .
Figure
Page
19. Response of TwoTemperature Feedback Model 
Step Input Volterra Series . . . . 82
C1. Block Diagram of TwoTemperature Feedback Model . . 94
Abstract of Dissertation Presented to the Graduate Council
in Partial Fulfillment of the Requirements for
the Degree of Doctor of Philosophy
THE APPLICATION OF VOLTERRA SERIES AND NONLINEAR
OPERATORS TO NUCLEAR REACTOR KINETICS
By
Ira Thierer
April 22, 1967
Chairman: Dr. Rafael B. Perez
Major Department: Nuclear Engineering Sciences
The analysis of nonlinear reactor systems through the use of
solutions to the equations which are descriptive of the systems have
been investigated. Two different methods of solution have been
examined. The first, solution by Volterra series, deals with the
system strictly in the time domain. The other method, nonlinear
operator analysis, yields results in a multidimensional transform
domain.
The Volterra series solution for a reactor system is derived
from the integral equation representation of the system and some of
its properties are presented. For a particular reactor model, it
has been proven from the Volterra series that a reactor with positive
viii
feedback at all times can have a range of power levels for which it
is absolutely stable. In addition, the basic Volterra expansion is
vital to the operator analysis method. However, it is a poor method
for the calculation of system response.
For the method of analysis by operators, the appropriate
algebra of operators is described and applied to a reactor system.
Also, the concept of the multidimensional transform is presented.
The response of a reactor system is found by a combination of the
operator algebra with the transforms. Based on the operator method,
some stability criteria are developed for systems with linear feed
back. It is found that linear stability is necessary but not suffi
cient as there are some additional restrictions on both the system
and the input to the system when nonlinearities are considered. Based
on comparison with simulated system;response, the operator method is
shown to give results of good accuracy.
Some possible applicationsand extensionsof the operator
method are discussed.
CHAPTER I
INTRODUCTION
In an invited paper presented at the American Nuclear Society
meeting of June, 1962, Jack Chernick begins by stating, "The motiva
tion for the present session may be explained by the bare statement
that reactor dynamics is basically nonlinear. Recent studies of the
stability of simple reactor systems have demonstrated serious inade
quacies in the linear theory."(1) This also describes most directly
the motivation for the work that follows.
Unlike linear systems for which there exists a coherent method
of analysis giving both qualitative and quantitative information, non
linear systems and, of particular interest here, reactor systems, have
been examined by a number of methods depending on the type of informa
tion desired. For reactivity inputs of sufficient amplitude to cause
the reactor system to move out of the linear analysis domain into a
quasilinear one, the describing function technique has been applied.
In this case, a sinusoidal input to the reactor produces an output con
taining the same frequency as the input but with several harmonic
frequencies superimposed on it. The describing function, a function of
both the amplitude and frequency of the system input, is the ratio of
the amplitude of the fundamental component of the output to the ampli
tude of the input. It can be used in stability analysis of the system
being examined once it has been obtained. Wasserman has employed this
1
technique in the analysis of data from SPERT tests.(2) One of the
major drawbacks of the technique is the limited range of input ampli
tudes over which results are valid. Also, when the describing func
tion ,is derived by functional analysis it is seen to depend only on
the oddorder kernels. Thus results for systems with relatively large
evenorder kernels are bound to be incorrect.
A second method attempted for nonlinear description of reactor
kinetics is based on the principles of classical mechanics. The "equa
tions of motion" for a reactor have been determined and Welton actually
has established some stability criteria for nuclear systems.(3) The
method should, in theory, also give the time dependence of various
reactor variables. However, for anything but the very simplest cases
the mathematics of the problem rapidly becomes too unwieldy.
Two additional techniques have been used to obtain purely
qualitative information concerning reactor stability. One deals with
the Laplace transformability of solutions to the point kinetics equa
tions and the Laplace transform of the feedback kernel.(4) Since
these are unidimensional transforms the feedback models must neces
sarily be restricted to linear ones. The authors of this method con
clude that one of the conditions for boundedness of the solutions to
the kinetics equations is:
JG(t)dt < 0 (11)
where G(t) is the feedback kernel of the reactor system. Comment on
this condition is made in Chapter III of this work.
The other technique for solely qualitative purposes is geo
metric and is based on work by Poincare' and Liapunov.(5,6,7) First
the equilibrium operating states of the reactor are determined. Then,
each equilibrium state is asymptotically stable if a continuous posi
tive definite function of the independent variables can be found which
also has a negative total derivative in a given region about the equi
librium state. The equilibrium states must be examined separately and
for a reactor there may be many; However, the main problem here is the
lack of definiteness in the method. The inability to find a Liapunov
function for an equilibrium state is not sufficient to preclude sta
bility of that state. Also, for mathematical manageability, the models
used for the reactor and feedback must be kept simple.
Since all of these methods have found only limited success in
their application and are often difficult to implement mathematically
in anything but the simplest cases, it is felt that a completely differ
ent approach to the nonlinear problem might yield more generally
useful results. But, rather than tackle the problem to obtain a direct
method yielding rigorously correct results, the author has decided to
borrow an approach used commonly in algebra: approximation by series.
This approach has two distinct advantages. First, for timedependent
solutions to the kinetics equations, the series can usually be made to
converge as closely to the true solution as desired by increasing the
number of terms included. Secondly, depending on the formulation of
the series, there exists the definite possibility of obtaining stabil
ity criteria for the reactor system under investigation by examining
either the whole series solution or just part of it.
Actually, some work has been done in this area. An attempt
has been made, and is continuing, to use Wiener functionals in connec
tion with experimental methods to analytically represent nuclear sys
tems.(8,9) However, the work of Barrett (10) is a good indication
that the method of series solutions represents a strong approach to
timedependent solutions and stability analysis in nonlinear systems
and is worth investigation.
Barrett considers a particular nonlinear differential equation
and obtains a Volterra series solution for it. Then, by examining the
regions of convergence of the series it is possible to predict directly
the stability properties of the system represented by the differential
equation. Since regions of convergence of many kinds of series have
been established and a systematic method for obtaining them exists for
the cases where they may not yet be known, the series approach can be
thought of as a reduction of the problem to almost solely an algebraic
one, albeit not necessarily a trivial one.
Although interest usually centers around the stability of
systems, an unstable system may have a sufficiently long time constant
so that there is adequate time for external means to be employed to
offset the instability when it occurs. It is then also of value to
know something about the time response of an unstable system as well
as the stability properties of any system.
In the work that follows two methods, both based on series
representations, are presented for analyzing,reactor systems. The
first method, presented in Chapter II, involves the integral equation
description for reactor kinetics and solution of the equation by
iterative techniques to yield a Volterra series expansion. Some prop
erties of this type of solution are given along with timeresponse
calculations for some examples.
In Chapter III, the maximum values of the integrals appearing
in the Volterra series solution found in the preceding chapter are
used to determine absolute stability bounds for reactors. Associated
with this chapter is a derivation appearing in Appendix D.
The second method of analysis is based on the theory of non
linear operators and is presented in Chapters IV and V. In the first
part, Chapter IV, the algebra of operators is described and the multi
dimensional transforms, used to convert the abstract operators to
calculable quantities, are explained. Application of the operator
theory to a general reactor system is made in the second part, Chapter
V, and some stability conditions are examined. In addition to this
chapter, some properties of operators and their manipulation are
found in Appendices E and F.
Chapter VI presents some examples and results of calculations
based on operator theory. These are compared to results obtained by
other means.
At various points in the work the stability of systems is
discussed. Unless otherwise specified, the type of stability under
consideration is Lagrangian stability where any bounded input to the
system results in an output which remains bounded. Thus a system may
have a limit cycle and it will still be considered a stable system.
Overall, the methods have been kept as general as possible.
No assumptions have been made about the nature of the feedback until
6
absolutely necessary. Then linear feedback models have been chosen
because of ease of calculation. For the neutronics models, it appears
that spacetime separability is necessary. There is no loss of gener
ality in the methods, as the material presented is sufficiently complete
for use with nonlinear feedback models as well as linear ones.
CHAPTER II
VOLTERRA SERIES
Introduction
One of the more easily applied and potentially useful approaches
to the solution of integral or integrodifferential equations is the
functional theory developed about thirty years ago by Volterra.(ll)
In this chapter the functional expansion is derived with particular
application to nuclear reactor systems and its usefulness is examined
by means of specific examples.
The differential equations for pointmodel reactor kinetics are
transformed into a single integral equation for which a functional
series solution is then obtained. The integral equation so obtained is
in a form which makes it applicable over a broader range of conditions
than those implied in the pointmodel. However, in those cases where
this particular integral equation is not valid, this chapter can still
serve as an example for solving the equation based on the reactor
model chosen.
The development presented here is basic to a large part of the
work that follows and reference will be made to it often.
Derivation of the Integral Equation for Reactor Kinetics
In terms of macroscopic phenomena the fundamental expression
describing the kinetics of the reactor neutron field is the Boltzmann
neutron transport equation. However, for the methods of interest here,
this equation requires simplification to be of value. Thus by assuming
spacetime separability, the presence of only monoenergetic neutrons,
and the applicability of diffusion theory, the reduction of the Boltz
mann equation to the pointmodel kinetics equations is well established,
These equations, in terms of power rather than neutron population,
are:
dP k(+S)1 N
i=1
dt 1 i C
(21)
dCi kpi
d = IP XiCi
dt Ii* i
where P. = the reactor power level,
C. = the concentration of the ith delayed neutron group
precursor,
k = the effective multiplication factor,
i* = the neutron mean lifetime,
Xi the decay constant for the ith group of delayed neutrons,
and
Pi = the effective fraction of neutrons appearing in the ih
delayed group.
There is no loss of generality by assuming that the reactor is at a
critical or supercritical state so that a source term can be omitted.
The equations are treated in the standard manner by expressing
all variables as having a steadystate component and a transient com
ponent with no assumptions being made about the relative size of the
two components. Substitution of these components into Equations (21)
and appropriate elimination of some of the steadystate components
leaves:
dSP tko+6k(t)](lp)l N
S 15P (t) + P P8k(t) + X iC(t)
dt 1* 1 iCi(t)
i=l
(22)
d6C [k +Sk(t)]p p(
dt 1 O 5P(t) + k(t) Xci(t)
t 1* P(t)
where the subscript o indicates the steadystate component and the
prefix 5 indicates the transient component of the variable to which
it is attached. The Laplace transforms are taken of Equations (22)
and combined into a single equation by the elimination of the delayed
neutron precursor variables. Inverse Laplace transformation of this
equation yields:
t rt
5P(t) = Po f g(tT)fk(T)dT + g(tT)6k(T)P(T)dT (23)
o fo
where g(t) = the inverse Laplace transform of Gl(s), and
N
1s Z P
Gl(s) = i= N L
1k +s s*+k Z
o o s+i
i=
At this point the variation in the effective multiplication
factor, 6k(t), requires some examining. This variation can arise in
two ways: 1) externally induced variations such as control rod movement
and 2) internally produced variations associated with the generation of
energy within the reactor. The former is generally independent of the
phenomena within the reactor coreand analytic representation depends
on the specific nature of the induced change.
latter depends on the reactor power level and
and is considered to be feedback variations.
may be assumed to be analytically represented
possible case by:(12)
On the other hand, the
its entire past history
This can be nonlinear and
in the most general
kf(t) = r... dTn K[t'Tl ..,tT;P(rl)...P(Tn)]
o o (24)
where K is some feedback mechanism kernel. In most cases of interest
in this work the very general expression can be reduced to:
(25)
kf(t) = K[tT;P(T)]dT
o
which can still take into account nonlinearity in the feedback mecha
nisms.
Then, if the externally introduced variations in the multipli
cation factor are denoted by ke(t) the total variation is given by:
t
5k(t) = ke(t) + K[tT;P(T)]dT. (26)
Here the feedback variation in the effective multiplication
factor is added to the external variation. The actual sign of the feed
back is taken into account in the kernel K. Insertion of this expression
for the total variation into Equation (23) results in:
t t T
5P(t) = P g(tT)ke(T)dT + Po g(tT)dTf K[Ta;P(a)]da
t
+ g(tT)5P(T)ke(T)dT
t 7I
+ f g(tT)5P(T)dT 7Toa;P(o)]da. (27)
o0 0o
This is the integral equation which will be solved by the Volterra
series technique.
Volterra Series Solution
Before anything further can be done with Equation (27) a model
for the kernel R(tT;8P(T)] must be chosen. Linear feedback has been
initially selected for the sake of simplicity. The kernel then becomes:
K[tT;6P(T)] = Kl(tT)5P(T) (28)
where Kl(t) is the linear feedback kernel. Introduction of this feed
back model into Equation (27) yields:
t t T
8p(t) = g(tT)ke(T)dT + g(tT)dT Kl(Tr)Po6p(a)da
+ g(tT)6p(T)ke(T)dT
t T
+ g(tT)5p(T)dTf Kl(Ta)P6p(a)do (29)
where 6p(t) = 8P(t)/Po. It is evident that even with a linear feedback
model this equation remains nonlinear because of the presence of pro
duct terms of the form 5p(t)k (t).
The solution of an integral equation of the class to which
Equation (29) belongs is found, according to Volterra,(1l) by
iteration. An initial solution, 5Po(t), is assumed and substituted
into the righthand side of the equation to give a first approximation
to the true solution. This approximation is then substituted into the
righthand side of the equation to obtain the next approximation and so
on until the total solution is obtained. This is generally an infinite
process and the series obtained is not necessarily convergent.
For Equation (29) the initially assumed solution is:
t
8P((t) = g(tT)k (,)dT (210)
0o
which is equivalent to zeropower reactor kinetics. The first approxi
mation derived from this is:
pl(t) = g(tT)k(T)dT + g(tT)g(Tro)ke(T)ke())doad
+ Pofff g(tT)g(op)K1(Tor)ke(p)dpdad
+ Po f g(t)g(Tp)g(g)Kl( )ke(p)k()
0000
dpdadpdT ,
(211)
and does not have as simple an interpretation as Equation (210)
because of the presence of feedback and nonlinear terms. The next
0 1 2 3
iteration will give terms containing P P P and P and even
o o o o
tually all powers of P will be present. For a nuclear system all
kernels must be physically realizable, i.e.:
g(t) 0
for t < 0.
Kl(t) 0
No further approximations will be shown here because the number
of terms appearing in each succeeding iteration increases rapidly.
In fact, the number of terms in any particular approximation is:
Nm = (Nm,+1)2
No = 1 (212)
where Nm is the number of terms in the mth approximation. It is easily
shown that by the third iteration there are 676 terms present. As yet
no way has been found to write these series of everincreasing numbers
of terms in a convenient form.
One of the points to be noted about the sequence of iterations
is that a particular iteration causes only the addition of more terms
to the terms already present from previous iterations. There is no
change in the form of those already existing terms. Thus the set of
iterations po(t), 6p1(t), 8p2(t), ... can be considered a sequence
of functions defined over the interval [O,t]. In order that the
sequence be uniformly convergent to a limit function 5p(t), which is the
true solution to Equation (29), on the interval, it is necessary and
sufficient that to each E > 0 there correspondssome integer N indepen
dent of t such that
bp (t)6pm(t) < E
whenever NSm and Nsn.(13) Since the iteration numbering scheme has
started with n=O, it is permissible to have N=0.
For 5pn(t), where n>l, no new physical quantities are added to
the solution, the only additions being higher order products of the
kernels present in 6Pl(t). From an intuitive standpoint then, it seems
reasonable that, for many cases, N=1.
Let us look at the case where N=0 and consider the following
inequality:
6po(t)6pl(t)
p(t) <
(213)
which is a representation of the statement that the result of the
first iteration is very nearly equal to the initial assumption. This
equation would then give the regions where the initial assumption
(zeropower kinetics) is a good representation of the total solution.
Substitution of Equations (210) and (211) into (213), combined with
the additional assumption that feedback is negligible, i.e. Kl(t) = 0
for all t, results in:
t T
Itg(tT)ke (T) g(To)k (cr)dadT
< 1 i. (214)
tt
g(tT)ke (T)dT
0
This yields the following condition on ke(t):
(215)
t
oIk ( )dT <<
e Imax I
o
The derivation of this condition and the definition of gmax appear
in Appendix A.
It is to be expected that ke(t) will be very small and an
example will show this clearly. Consider a critical reactor with
Gl(s) as defined in connection with Equation (23). Only one group
of delayed neutrons will be used with p = 0.007, X = 0.08, and 1* =
. Then, for this case:
10 Then, for this case:
gmax e g(0) = lim [sG(s)] = 9930
max S8c
so that:
t
Ske(T)IdT << 1.006 X l0"4
which is indeed a very small total allowable k .
If the feedback is assumed to be negative and not necessarily
negligible, the conditions placed on k (t) from Equation (213) will
involve the initial power level of the reactor in addition to the
kernels Kl(t) and g(t). These conditions could be less restrictive
than the condition in Equation (215), but they are not immediately
obtainable from 6po(t) and 6Pl(t). However, some conditions will be
found from a slightly different treatment of the problem in the next
chapter.
Examples
Two examples have been chosen to demonstrate the effectiveness
of finding timedependent responses by series solution of the integral
equation for reactor kinetics. One of the examples uses a theoretical
feedback model for natural circulation boiling water reactors. The
other one is based on work appearing in a book by Schultz (14) and
will also be used in later chapters for comparison of methods.
Consider the boiling water reactor problem first. The theo
retical background for the feedback model is not covered here,as it is
adequately described by Akcasu who developed it but a block diagram
for the feedback loops is given in Figure 1 where the various gains,
Ck, are related to the reactivity coefficients for the reactor and
Figure 1. Feedback Model for Boiling Water Reactors
the T and T. are the time constants for the processes occurring in
the reactor.(15) For the Experimental Boiling Water Reactor (EBWR), a
set of values for the constants is given in Appendix B. With these
values, the feedback kernel is found to be:
Kl(t) = 3.50X106 + 9.75X10"3 e"o.685t 0.2033 e"2.27t
0.0604 e2.63t + 0.2685 e"3.12t
+4.62XO14 e4t 0.0584 e572t
+0.0514 e"9'09t 8.62X10"3 e"33"33t (216)
where Kl(t) is given in units of cents/KWsec. This function is shown
in Figure 2, where it can be seen that the feedback kernel, although
initially negative, becomes positive after about 1.7 seconds and
remains positive for about 13 seconds more. After this it asymptoti
6
cally approaches a constant value of 3.50X10 cents/KW'sec.
The zeropower kinetics kernel to go with this model is the
same one as defined in connection with Equation (23). With the con
stants given in Table I, this kernel becomes:
g(t) = 11.00 + 46.78 e0O1847t + 9872.5 e701013t. (217)
For a 8function input and a given reactor power level, the
time response of the EBWR system can now be calculated based on the
various iterative approximations. For an initial operating power level
of 10 megawatts (MW) and for input multiplication variations of ke
0.01 and ke = 0.001 Figure 3 shows the initial approximation, 6p (t).
K,(t)
cents
KW.sec .012
.016
.020
.024
I I I I I
o 2 4 6 8 10
Time(sec.)
Figure 2. Feedback Function Based on EBWR Data
8po(t)
0.3
0.2
0.1
0 5 10 15 20 25 30 35 40 45 50
Time(sec.)
Figure 3. First Approximation in the Volterra Series
TABLE I
Constants:forthe ZeroPower Kernel
PI 0.002074
P2 0.004896
X 0.02524 sec"
X2 0.56409 sec1
104 sec.
Figure 4 shows the approximation to the solution from the first iter
ation, 6Pl(t), for the same initial power level and input variations.
For purposes of comparison, the time response has also been
calculated using linear theory. This would be a good representation
of the system for small variations in the input. The smaller of the
two inputs used to obtain Figure 4 and Figure 5 would cause the system
to fall into the linear class or come very close to it. The linear
response is given by the standard equation:
()Pli( L1 G1(s)k (S)
lin(t) L PoKl(s)G(s (218)
where e(s) is the Laplace transform of the input variations, L" is
the inverse Laplace transform operator, and the rest of the notation
has already been defined. Figure 5 shows 8plin(t) for the EBWR
model given in Figure 1 and for the associated kernels, Equations
(216) and (217). The scale on the left is for ke 0.01 and the
scale on the right is for ke = 0.001. Since the input is an inpulse,
the only change in the linear response due to a change in the
0.3
s8p(t)
0.2
5 10 15
Figure 4.
k = .001 40
e
P 104 KI
o
30
8p1(t)
20
a.
10
0
10
SI I I I I10
20 25 30 35 40 45 50
Time(sec.)
First Iterative Solution EBWR Problem
Scale Scale
for for
k =0.01 k =0.001
e e
8 \ 0.8
P = 10 Kw
0.4
5Plin(t)
0
Envelope of Response Curve
S0.8
8U 0.8
15 20 25 3
Time(sec.)
Figure 5. Linear Response EBWR
magnitude of the input is in the amplitude of the response function.
The shape of the function is not affected.
The second example involves a reactor with twotemperature
feedback mechanisms.(14) A block diagram for this reactor model and a
derivation of the feedback kernel are given in Appendix C. The zero
power kernel is based on a single effective group of delayed neutrons
and is given by:
s+T 1
g(t) = L i[ (s
where X is the effective decay constant for the delayed neutron
precursors,
p is the fraction of neutrons appearing as delayed
neutrons, and
(219)
1* is theimean neutron lifetime.
From the model and Appendix C, the feedback kernel is:
Kl(t) = Li [ (l+T (220)
where A, T1, and T2 are related to the heat transfer and feedback
parameters of the reactor. The data used in the calculation of the
time response of the reactor are shown in Table 2..
Figure 6 is the result of the first iteration, 85p(t), for
input variations ke 0.01 and k = 0.001, both 6functions, and with
the data given in Table 2,. Again, for comparison, the linear response
of the system has been found from Equation (218). This is shown in
Figure 7 with the same scales as indicated for Figure 5.
TABLE 2
TwoTemperature Feedback Model Data
1
0.125 sec1
1" 4
1* 10 sec.
P5 x 10"3
T1 1.0 sec.
T2 0.1 sec.
A 1.6 X 10'7/BTU
Po 100 MW
Results
It is immediately evident from Figures 3, 4, and 5 that the
initial approximation and first iteration only begin to approach the
true solution for the time response in EBWR. If the case ke 0.001
e
is considered, the linear response should be very close to the true
response. The first iterative solution is just beginning to show some
of the oscillatory nature found linearly,while the initial approxima
tion does not show any. In addition, the linear solution asymptotically
approaches zero which neither of the iterative approximations does.
However, the 5pl(t) solution is generally closer in magnitude to the
linear solution than is 5po(t).
For ke = 0.01, it is unlikely that the linear solution is
P1l(t)
6
8
10
0.2 0.4 0.6 0.8 1.0 1.2 1.4
Time(sec.)
Figure 6. First
Iterative Solution TwoTemperature
Feedback Model
200
4oo00
6pl(t)
800
1000
1.6
Scale
for
k =0.01
e
plin
Time(sec.)
Figure 7. Linear Response TwoTemperature Feedback Model
Scale
for
k =0.001
e
.1
.2
.3
0.4
still a good approximation to the true solution and the comparison of
iteratively obtained solutions to it would have little meaning. From
Figures 3 and 4 it appears that the approximations are diverging,
although the shape of the time response remains the same. However,
without calculation of further approximations there is no way of pre
dicting what the sequence of iterations will do.
In the case of the twotemperature feedback model, Figures 6
and 7 show that the first iteration, 8pl(t), already exhibits somewhat
the same shape as the linear response. The magnitude of the iteratively
obtained response is quite different from the linear response but the
difference is relatively smaller for the smaller variation in input.
One further word needs to be added about the results before any
conclusions are drawn. The amount of calculational effort required to
obtain the first iteration, 8pl(t), which only has four terms, is large
even for the simpler twotemperature feedback case. Calculation of time
responses for the next iteration, 5p2(t), with twentyfive terms is
almost prohibitive.
Conclusions
Although a Volterra series can theoretically be obtained, as
outlined in this chapter, as a solution to the integral form of the
kinetics equations, the length and complexity of the series makes it
nearly impossible to determine its convergence properties. Without
this knowledge, nothing can be said about the quality of the various
iterative approximations to the final solution. If convergence is
assumed, the results show that for anything but very small variations
in the input to a reactor system, a first iteration on an initial
approximation to the true solution is not sufficient to obtain a good
representation of the solution. From a practical standpoint, however,
a calculation of the second iteration is not feasible even with the
aid of highspeed digital computation devices.
In the following chapters use will be made of the Volterra
series solution to find convergence and time responses less directly
but, particularly in the case of time response, more accurately.
CHAPTER III
MAXIMUM VALUE SERIES
Introduction
In the previous chapter it has been noted that the infinite
Volterra series solution for the timedependent response of reactor
power is too complex to yield sufficiently useful information directly.
Nevertheless, it is desirable to obtain some convergence properties
for the series so that the first few iterations in obtaining the solu
tion could be used for calculational purposes, even under limited
circumstances. One of the ways to do this is to find a series which
will always be equal to or greater than the true series solution. If
convergence can be obtained for this larger series, it is then guaran
teed for the true series. In this chapter, the larger series is
determined and is referred to as the maximum value series because it
is based on the maximum values of the integrals of the various kernels
involved.
Development of the Maximum Value Series
The maximum value series is obtained from the basic inequality
given in Equation(32) below. If:
b
Ir f f(x)g(x)dx (31)
a
then:
b b
III b f(x)g(x)dxl J If(x)I g(x)ldx (32)
a a
where If(x)j = magnitude of f(x).
If jg(x)[ is always less than or equal to some maximum value,
gmax' on the interval (a,b], then:
b b(
I 1 f If(x)lgmaxdx'= gmax f if(x)ldx (33)
a a
because gmax is a constant.
Now the application of Equation (33) to the initial approxi
mation for the Volterra series, Equation (210), results in:
IPo(t)I < kM ft g(tT)dT (34)
0
where kM is the maximum value of Ike(T)I on the interval [O,t]. The
following quantities can be defined:1
t t
f g(ir)dT f jg(tT)JdT a G
o o
t
f IK(T)IdT f IJKl(tT)IdT a K. (35)
o o
1The quantities defined in Equation (35) are represented by
symbols having neither subscript nor argument. The time, t, is implied.
Since all the ranges of integration in the various approxima
tions to the complete solution are over the interval [0,t] or some
smaller interval, an expression like the second term on the righthand
side of Equation (210) becomes:
t 7
1121 = I g(tT)ke(T) f g(TC)ke(a)daddTl
0 0
2 t t
S, k2 I g(tT) f g(T0)dadT
o o
t
: k[F Ig,lT12 2 2 (36)
0
Thus, many terms which have the same kernels but differing
orders of integration of these kernels can be shown to have the same
maximum value over the interval (0,t]. This greatly simplifies the
infinite series.
With the definitions above, the first iteration can be
rewritten as:
16pol s kMG (37)
and the second iteration is immediately seen to give:
ISl1 kMG + G+ o2 o 2K + P kG2K. (38)
A similar inequality exists for each approximation until the final
solution is obtained. For the final solution it can be shown, after
some arithmetic manipulation, that the maximum value infinite series
is given by:
PI (klG)1 Lj+)!(i (P KG) (39)
j=0 j!i(J+'). 0 (i+)
Consider the linear terms only. For these, j=0 and the series
obtained is:
Plinl [1 + PoKG + (PoKG)2 + ..... ](k G). (3.10)
When P KG<1, this reduces to:
kMG
1Plinl P KG (311)
which is very much like the familiar form of the linear transfer func
tion of a system. The condition required to be able to reduce the
series to a closed form is related\to a linearly determined stability
condition.
Convergence
Equation (39) can be shown to be similar to a hypergeometric
function of higher order in two variables.(16) In fact, in the nota
tion used for this type of function:
0    0
Ip P F k( i+2 kG, 1. (312)
1 a,2
Appendix D has a definition of this notation. The only difference
between the hypergeometric functions tabulated in the literature and
the.summation in Equation (39) is the appearance, in the latter case,
of one of the running summation indices in the coefficient as well as
the exponential part of each term. This does not make any difference
in obtaining convergence criteria.
In considering convergence of the righthand side of Equation
(312), it is obvious that the constant will not have any effect,so
that convergence need only be obtained for the hypergeometric portion
of the expression. It has been determined that this expression con
verges when:
(G) + (PKG)2 < 1. (313)
A derivation of this condition also appears in Appendix D. The region
of convergence associated with this condition is shown in Figure 8.
In all of the above work a time, T, is implied for the upper
limit of integration of the kernels. Convergence is guaranteed up to
this time if the condition in Equation (313) is satisfied. For con
vergence over all times the integrals must be taken to infinity. How
ever, for reactor systems, convergence for times up to about a minute
is all that is necessary because it can usually be assumed that some
external means can be applied to control divergent power fluctuations
if the divergence has ah efolding time of the order of several seconds
or more.
For any particular initial operating power level, Po, and for a
known input variation in multiplication factor the absolute stability
of the reactor system can be checked by use of Equation (313). Once
it is determined that the system is stable, i.e. the maximum value
series for the system is convergent, it may be desirable to know
what the maximum value of the series is because this is the maxi
mum value of any output power variations. Although the series in
PoKG
Figure 8. Region of Convergence of the Maximum Value Series
Equation (39) is not in a form amenable to calculation, it can be
simplified somewhat. In Equations (310) and (3rll) the linear terms
have been extracted and reduced to a simpler form with a certain added
condition, P KG<1. The second, j=l; third, j=2; and higher order
terms can be handled in a similar manner with the same condition.
Since all the quantities in Equation (313) are nonnegative, the con
dition P KG<1 needed to obtain the simplification will always exist
if Equation (313) is satisfied. For terms up to fourthorder, the
simplifying process then gives:
kMG (kMG)2 (1+P KG)
P KG (1P KG) + (1PoKG)5 (kMG)3
1+3P KG+(P KG)2
+ +3P0KG+(P KG (kMG)4 + ... (314)
(1P KG)7
which is a much better form for calculation of the maximum J1pl.
Results
To demonstrate the conditions derived above, the EBWR model
used in Chapter II has been chosen. Equations (216) and (217) give
the kernels Kl(t) and g(t) respectively. Figure 9 shows the regions
of convergence for this model. The various curves are for different
integration times.
For any particular time, T, for which convergence of the
maximum value series is desired, the values of G and K are constant.
Thus, the allowable values of kM can be found for any particular initial
power level or, in an inverse manner, the range of allowable initial
power levels can be foundwhich will keep the system absolutely
T 5 sec.
T = 10 sec.
T = 20 sec.
T = 60 sec.
T 300 sec.
2.4x10"3
Figure 9. Regions of Absolute Stability EBWR Model
P
(KW)
1200
1000
8oo
6o0
4oo
stable for any particular kM. This is the reason for the selection of
the axes used in Figure 9.
From Equation (217) it can be seen that the value of G, as
defined in Equation (35), approaches infinity when the integration
time goes to infinity as a limit. This is reflected in the decreasing
region of convergence with increasing integration time in Figure 9.
As an example of the calculation for the maximum power vari
ation, an integration time, T=20 seconds, is considered and a point
which falls within the region of, convergence for that particular time
is chosen: kM = 4X10" and Po = 200 KW. For T=20 seconds, the values
of K and G are 1.64XO10 /KW and 610 respectively. With these quanti
ties, Equation (314) yields:
18pl < 0.305 + 0.1162 + 0.0533 + 0.0277
+ 0.0157 + 0.00955 ..... (315)
where each quantity corresponds to a particular order term in Equation
(314) and the sum approaches an approximate value of j6p( 0.535.
This means that, for an initial operating power level of 200 KW and
an input variation of 4X10"4, the reactor system will be stable and,
upon remembering that 6p = 8P/Po, the fluctuation of output power can
never exceed 107 KW around the initial power level. This seems to be
a large fluctuation but it also must be remembered that the feedback
kernel, Kl(t), is treated as if the feedback is positive at all times
in order to obtain K.
Conclusions and Comments
Since the conditions established in this chapter are based on
maximum values of integrals, they can be considered only as bounds
beyond which information becomes ambiguous. For example, it can be
stated that a system is absolutely stable if its kernels are such that
the inequality in Equation (313) is satisfied. On the other hand, if
the inequality is not satisfied, the system may or may not be unstable.
The stability will depend on how much difference exists between the
maximum values of the integrals and their true values. If the maxi
mum value series diverges but the various integrals have true values
which are sufficiently smaller than the maximum values, the true series
could converge. This ambiguity shows up most when the feedback in the
system is always negative.
In addition, the maximum value of the variations in the output
has been calculated from the maximum value series. For the particular
data given previously, intuition suggests that the true range of vari
ations is very much smaller than the maximum value indicated by
Equation (315) because the allowed maximum value of the input varia
tion is very small and the feedback is negative for at least part of
the time.
However, neither of the above conclusions can overshadow the
most significant conclusion that can be drawn from this chapter. It
has been proven here that a system can be stable under a range of input
and operating conditions even with feedback that is positive for all
times of interest. The derivation of the maximum value series has
been accomplished using only the magnitudes of the kernels and the
conditions obtained from this series are, therefore, applicable to
positive as well as negative feedback. This indicates that the condi
tion given in Equation (11) is definitely overrestrictive.
If the feedback kernel shown in Figure 2 was positive for all
times but with the same magnitude the negative portions now have, the
results in Figure 9 would be unchanged. The only difference in the
problem would be in the degree of proximity to the maximum value
quantities the system would achieve. The allpositive feedback case
would give true variations in power that are closer to the maximum
than would partially negative feedback.
The simplification of the maximum value series to a sum of
terms, each being of different order in kM, and the similarity between
the linear term in this series and the standard linear transfer func
tion for a system leads to speculation that the actual, nonmaximum
value series could also be reduced somehow to a sum of terms arranged
in increasing orders of ke(t). Calculation of time responses should
be easier for this type of sumand stability conditions might be more
readily available than with the infinite series of terms, most of which
contain multiple integrals.
In the next chapter a more compact series of the type just
indicated is derived by a completely different method but use is made
of the Volterra series already obtained in Chapter II.
CHAPTER IV
NONLINEAR OPERATORS THEORY
Introduction
The series solutions examined in the previous chapters have
been obtained by iterative techniques and theoretically do not present
total system response characteristics at any level of approximation.
For example, with the iterative series it takes an infinite number of
terms to obtain the linear responseor, for that matter, any other
order response of a system. However, certain information is more
readily extracted from this iteratively obtained series if there is
some knowledge of the form of the final infinite series. This has been
demonstrated in the preceding chapter where absolute stability bounds
have been established.
In order to obtain more of the information present in solutions
to the reactor kinetics problems and to achieve a more convenient form
for calculation of these solutions, a completely different approach to
the problem is taken in this chapter. The reactor system, described
in terms of several operators, is analyzed by using the algebra of
operators and multidimensional Laplace transforms to obtain a time re
sponse for the system. This approach gives a more compact series as
well as one in which response characteristics of any particular order
such as linear response appear in a relatively simple closed form.
The rigorous theory backing up the work with operators has been
examined and describedby Brilliant,(l7) George,(18) and Zames(19) and
will not be presented here. There are also some applications to
electrical engineering problems by Parente.(20) However, a description
of the algebra of operators is included because of its possible appli
cation to almost all types of systems.
Definitions
In order to classify the concept of an operator in the hier
archy of functions, it is first necessary to establish the fundamental
concept of a function. In abstract terms, a function is any determi
nate correspondence between two sets of objects, usually called the
domain and the range. The most commonly encountered functions, the
socalled "real functions," are those whose domain and range are sets
of real numbers. An example of this is the function y = exp(ax); a>O,
with its domain specified as the set of real numbers on the interval
[0,0] and its range specified as the set of real numbers on the
interval [0,1].
The domain and range of a function need not be restricted to
setsof real numbers. It is permissible for either or both to be a
set of functions. Then, a function whose domain is a set of real
functions and whose range is a set of real numbers is called a func
tional. The definite integral falls into this category. If both
the domain and range are sets of real functions, the function obtained
is known as an operator1. A good example of this is Laplace
Operators will be denoted by underlined capitals, e.g. H.
transformation which is the correspondence between a real function of
time in the domain and a real function of frequency in the range.
Now consider the general system shown in Figure 10 where the
x(t) Nonlinear y(t)
System
H
General Nonlinear System
Figure 10.
input and output are members of sets of real functions. In many cases
the mathematical description of an analytic nonlinear system can be
given in terms of a Volterra series of the type obtained in Chapter II:
CO O
y(t) j hl(i)x(tT)dT + hT1l, 2)x(tT1)x(tT.)dTdT2
CO 00 CO
+ ..... + .... ( h ,. )xt(tT,) ...x(tT
n' n
dTl......dTn
(41)
where each term of the series is an operation on the input function
which produces a function that will be a part of the total output. Thus:
Yl(t) Hl[x(t)]
y2(t) H21[x(t)]
yn(t) = .[x(t)]
(42)
where H f hl(r)  dT,
rCO
A2 = Jf h2(Tl1'2) dT1dr2, and in general
00 2C
n n .. ....n
The total output is then:
y(t) = yl(t) + y2(t) + .... + yn(t) + .... (43)
which can be rewritten as:
y(t) (i + 2H + 3 + ..... +H + .... [(t)] (44)
or, upon combining all the operators into a single operator for the
system:
y(t) = H[x(t)] (45)
where H = H + H2 + 3 ...... + H + ....
The order of an operator, indicated by the subscript on the
operator, is found by determining the order a constant multiplier to
the input function has at the output. For example, if the input is
x(t) = Az(t), and the operator is H2[x(t)], the output is y2(t) =
A H2[z(t)], indicating a secondorder operation. It can be seen that
operators of different orders can be added to produce new operators.
Also, the preceding shows how a system can be represented by a general
overall operator.
A problem arises when a system has multiple simultaneous inputs.
The principle of superposition holds for linear systems but certainly
can not be generally applied in other cases. For a secondorder
operator, the output due to two simultaneous inputs is:
Y2(t) = H[x(t)'+ z(t)] (46)
which, with the use of the definition of H2 given in Equation (42),
becomes:
Y2(t) = f h2 (1 '2)[x(tT1) + z(tT1)]rx(tT2) + z(trT2)]
dT ldr2 (47)
Upon expanding and using the operator definition in reverse, this
can easily be shown to be:2
y2(t) H= (x) + H2(xz) + 11(z) (48)
where H2(xz) ff h2(.T2)x(tTd)z(tT2)d1dT2., A secondorder
o? o)
operator acting on a sum of two inputs leads to an output made up of
three secondorder terms. Similar expansions can be derived for higher
order operators and for many simultaneous input functions.
Algebra of Operators
Figure 11 shows pictorially four methods for combination of
systems. Each system is represented by its characteristic operator
W which gives the correspondence between the input functions and the
2The time variable, t, will be omitted where it is clearly
understood.
X(t) C y (t
(a)
Addition
L() C y(t)
(b)
Multiplication
xCt B A
xt C
Simple Cascade
y(t) w(t) y(t)
P. .2 C 
J (d)
General Cascade
Figure 11. Combination of Operators
output functions for that system. The combined system is represented
by a single operator which is related to the individual system oper
ators according to the algebra that follows:
An input can be applied to two systems simultaneously and
their outputs can be either added, as shown in Figure lla, or multi
plied, as in Figure lib. These processes are represented, in the
standard notation, by:
C A + B (49)
and
C = A.B (410)
respectively. The addition and multiplication of operators are both
commutative and associative, thereby giving the following identities:
A + B B + A (4lla)
A + (B(C) (A+B) + C (411b)
and
A.B = B.A (412a)
A*(BSC) = (A*B)*C (412b)
The multiplication operation is also distributive:
A.(B+C) = (AB) + (A.C) = (B+C)*A (413)
where the last equality results from Equations (412a) and (49).
The third combination, shown in Figure llc, is the result of
applying an input to one system (B) and using the output of that
system as input to the second system (A). The output from the second
system is usually the one of interest. This is the "Cascade" operation
and is indicated by the symbol "*". This operation is given by the
expression:
c A*_B (414)
and is associative:
(A*B)*c = A*(B*C) (415a)
but is not in general commutative:
A*B 0 B*A. (415b)
The most important exception here is the case where A and B are both
linear operators. Then, these operators are always commutative in
cascade.
If the cascade operation precedes, in time, an addition or
multiplication, although the notation places the cascade operator
after the algebraic operation, the cascade operation is distributive.
If the cascade comes after either an addition or multiplication it
is generally not distributive Thus:
(A+B)*C = (A*C) + (B*C)
(416)
(A'B)*C = (A*C).(B*CC)
but
C*(A+B) ( (*A) + (LCB)
(417)
C*(A'B) s (c*A).(CCB)
In this latter case, with the cascade following any other
operation, it is important to determine what the combined operations
will yield. Consider a secondorder operator following a sum of two
operators of any orders, Figure lid with only the H2 component of H
present. If x(t) and z(t) are now the outputs of the operators A and
B, respectively, Equation (47) is still valid. Since the outputs
may be expressed as operators acting on the input to the various
systems, Equation (47) can be rewritten as:
y2(t) H2((A[w])2) + 2H2(A[w]B[w])
+ H((B[w]) ) (418)
where H2((A[w])2) = H*A. Equation (418) gives rise to the general
cascade operation denoted by "o":
(Ho(A.B))[w] = H(A[w].[Bw]) (419)
where 2o(A2) = H**A.
Two further operators which must be included to aid in making
the algebra complete are best explained by the following equations:
A + ( + A = A (420)
A*I = I*A = A. (421)
The zero operator, 0, can be thought of as a system for which there is
no output no matter what the input may be (e.g. an open circuit),
0[x] = 0, or, in more formal terminology, the zero operator has a
range which is a null set. The identity operator, I, is representa
tive of a system which gives an output equal to the input (e.g. a
short circuit), 1[x) = x. Formally the identity operator has a range
which is identical to its domain with a onetoone correspondence
between equal members of the two. The remaining operations of sub
traction and inversion are derived from Equations (420) and (421) by
the following:
A+ (A) = 0 (422)
A*A1 = (423)
The question of the existence of an inverse for a given operator is
not considered here. It is assumed that, for all operators of interest
in this work, unique inverses exist but are not necessarily physically
realizable.
It has already been noted that operators of different orders
can be added to produce a new operator. Similarly different order
operators can be multiplied and cascaded. The operator resulting from
the multiplication of an m order operator with an nthorder operator
is of order m+n. The cascade of these same two operators gives an
operator of order mn and their sum is an operator containing both
th th
m and n order terms. The order of other more complex arrangements
can be derived from these few fundamental properties.
Multidimensional Transforms
In order to practically apply the operator concept to systems
analysis there must be some way to relate the operators for systems to
kernels of the kind given in Equation (4i). By using the standard
Laplace transformation, the first term on the righthand side of that
equation can be transformed to yield the transfer function Hi(s),
which is nothing more than a linear operator. For higher order systems
a multidimensional, unilateral Laplace transform pair can be defined,
in a manner similar to the onedimensional case, by:
Fn(Sl,...,sn) f fn(tl, ..,tn)exp(ltl...sntn)
0 0
dtl....dt (424)
a+j 0+jcO
fn (t1 tn ) = ( )n .... f Fn(s1,...,s)
jo ja
x exp(sltl+...+s t )ds.... ds. (425)
If Equation (424) is applied to the kernels of Equation (41),
it is easily shown that, for various methods of combination of the
kernels resulting from addition, multiplication, and cascading of the
associated operators, the transforms of the kernels are added or
multiplied and their arguments depend on the method of combination of
the kernels. The combinations of transforms are given in Appendix E.
Consider the general term in Equation (41) and artificially
introduce tl, t2, ..., tn so that the term can be written:
Yn(tl,...,t) = .... f hn(T1'''"n )x(tl 1)...x(tnTn)
*CO CO
dTI ....dT. (426)
The result of applying the transform definition to this is:
Yn(sl...s ) = Hn(s ...,sn)X(sl)'..'X(s). (427)
Thus knowledge of the system operator and input in the transform domain
leads to an output form which can be inverted to the time domain.
Then, for the time response, all the ti must be set equal to yield
results equivalent to those obtained with Equation (41). It is,
however, more convenient to use George's method of association of
variables in the transform domain (18) to reduce the system function
to a single transform variable. This singlevariable function can
then be inverted to give a time response directly. For the second
and thirdorder operators this association is given formally by:
G3j w*
H(s) 2j H(su,u)du (428)
Ojo ajc
duldu2 (429)
and the association technique is easily extendible to higher order
operators. As an example of association, consider the secondorder
operator:
A B
H2(S s2) +a (430)
which inverts to:
h2(tlt2) = ABexp(at bt2) t, t2 2 0 (431)
and, upon setting tI = t2 = t, becomes
h2(t,t) = ABexp[(a+b)t] t 0. (432)
The transform of this is now
H(s) = s++ (433)
Thus, if H2(s1,s2) is treated formally by association, the result must
be:
H2( I A B AB
(s)  du sab (4A34)
) s 2j j su+a u+b s +a+b (
When several frequency variables are present they can be
reduced to a single variable by judicious selection and association of
two of the variables at a time until all have been associated. This
method replaces the many contour integration that would ordinarily
have to be performed, one for each frequency variable, and thus gives
a fairly simple and rapid way of reducing an unmanageable multidimen
sional function to a useful singledimensional one. Association
formulas based on Equation (428) and using the result obtained in
Equation (434) are given in Appendix F for a number of cases.
In summary, the algebra presented in this chapter helps reduce
the analysis of combined nonlinear systems to a manipulation of block
diagrams. Of course, the kernels for the various components must be
good analytic representations of the individual systems in order to
obtain accurate operators for the combination. Thus, the problem of
finding adequate models for certain systems still exists although
their combination has now been simplified.
CHAPTER V
NONLINEAR OPERATORS APPLICATION
Introduction
In this chapter the algebra of operators just presented is
combined with kernels derived in the manner given in Chapter II in
order to obtain an operator expansion for a nuclear reactor system.
Then the multidimensional transforms are applied to the expansion to
produce a form from which time response can be calculated and also
from which something can be inferred about the stability of the system.
Operator Reduction of Reactor Systems
The general reactor system to be considered is the one shown,
in standard representation, in Figure 12a and consists of a neutronics
operator, G, in the forward part of the loop and a feedback operator,
K. The neutronics operator gives the response of the reactor power as
an output for a timevarying multiplication factor as an input. The
feedback operator relates a changing multiplication factor to the
reactor power. At present the operators can be left as general as the
above statements will allow,but it is assumed that both overall opera
tors, G and K, can be represented by a sum of operators of various
orders up to a theoretically possible infinite order. Thus:
k +
e
Figure 12. Reduction of a General Reactor System Part I
G =" + +G ..... +G + .... (51)
S 1 =2 n
K 1 + K2 + ..... +K + .... (52)
1 n
where, for any particular case, all or possibly only some of the various
order operators are nonzero. For example, only K, is present for
linear feedback and all the other K are zero operators.
In order to facilitate the reduction of the system to a single
operator, the system is rearranged to become the system shown in
Figure 12b. This arrangement is easily seen to be equivalent to the
arrangement in Figure 12a. Now the cascade of K and G can be repre
sented by a single operator H, as in Figure 12c, where H is given by
the relation:
H K*G. (53)
If H is also assumed to be composed of a series of various order
operators, this relation can be rewritten in terms of the expansions
of H, G, and K as:
H + H + H + .... ( K2++K3 + ...)(+G 3+...).
(54)
Upon application of the distribution law for a cascade over a
sum of operators and with the appropriate general cascade operations,
the righthand side of Equation (54) becomes a sum of terms of
various orders involving both G. and K Equating terms of equal
order results in a series of equations giving the components of H as
a function of the components of G and K. For the first few terms
these equations are:
Hl cXI*G1 (55a)
2 *2 +2 (55b)
H3 = K1G3 + ~2(Gl'2) + K3 1 (55c)
and succeeding terms are found in a similar manner.
Next, the feedback loop containing H is converted to a single
feedthrough operation to obtain the system in Figure 13a. If k is
the output of system L when ke is the input then, in words, the output
of L is equal to the identity operation on the input plus the result
of system H acting on the output of L. In operator notation, with k
as the input, this combination becomes:
L(k = (I+HL)[k ]. (56)
Again, assuming a multiorder expansion for L, this last equation can
be rewritten as:
L1.2+L3+ ... = I + (H1+H2+H3+.... )*(L1+ +...) (57)
When the proper operator algebra is performed and it is noted that the
identity operator is always linear, a series of equations is obtained
which must be solved for the components of L as a function of the com
ponents of H. The first few of these equations are:
L h = + Hl*L (58a)
L 1*2 + H*L1 (58b)
L3 H*L + 2%eo( .L ) + H 3*L (580)
3 1 3 =~l2 31()
Figure 13.
Reduction of a General Reactor System Part II
By using the subtraction and inversion operations the compo
nents of L are found explicitly from these equations and are:
11 (7Hl)1 (59a)
L = (IH,) 1*H('1H 1 (59b)
1 11 (59b)
+ H *(Il)'l } (5.9c)
Finally the L and G operators are cascaded to produce the overall
reactor system operator R in Figure 13b:
R = G*L (510)
An expansion of the type already given combined with the proper oper
ator algebra produces the following equations for the first two compo
nents of R:
R1 GI*(IK!ll)'1 (511a)
_R 2 g*(Ia 1)" 1+ Gl*(KlGl)1* *[_Kl 2*(.1*1)"
(5llb)
where the quantities found in Equations (55) and (59) have been
replaced appropriately. The remaining components of R are found in
an identical manner.
Thus, the whole reactor system has been described by an overall
operator made up of operators of individual orders. The next require
ment is the development of a means of obtaining useful information
from the overall operator R. It is assumed that the component opera
tors, R, are related to kernels ri(tl,...,ti) in a Volterra series
expansion of the type given in Equation (41). Then the multidimen
sional transforms defined in Chapter IV and Appendix E can be used.
The feedback operator K can vary greatly from one reactor to
another and can not be further specified without knowledge of the kind
of reactor under consideration. On the other hand, one neutronics
operator can adequately serve for several reactor models. For conven
ience, the pointmodel given in Chapter II has been used again here.
Equation (23), with a normalized power variable, is solved by the
iteration technique of Chapter II to yield the following Volterra
series:
t t Tl
p(t) g(tT)5k(,)dT + f f g(tTltg(TlT2)6k(T( 6k(T2)dTldT2
o 0 0
+ f ... g(tTg(tr)(Tl').. g(Trn1 )k(T) .... k(T)
o0 0
dT n.....d .
(512)
When the multidimensional transforms are applied to the kernels
of this equation, with the reminder that all the kernels must be phys
ically realizable, it can be shown that the transforms of the components
of the overall neutronics operator, G, are given by:
G1(s) = Gl(s) (513a)
G2(sl,s2) G1(s+s2)Gl(s2) (513b)
G3(s81s2s3) = GI(s +s2+s3)G1( s2+s3)G(s3) (513c)
and similarly for the rest of the components.
Now, with the application of multidimensional transforms to
the operators of Equation (511) and the use of the combination rules
of Appendix E, a series of transforms of the components of R is ob
tained. Since the results presented in the next chapter depend on
the pointmodel neutronics given in Chapter II, the various order
neutronics kernels which are based on this model and are given in
Equation (513) are also used.
Higher order kernels can be obtained in terms of lower order
ones and for ease of presentation, although not necessarily for ease
of calculation, they are given this way in the equations below. The
first few transforms are then:
Gl(s)
Rl(s) l () (514a)
1Kl(s)G1(s)
Rl(sl+sg)Rl(s2)
R2(s 2 1Kl(s,)Gl(s1) + K2(sl's2)Rl( sl +s2)R1(s8)Rl(s2)
(514b)
R (s1,s2,s ) Ri(Sl+s2+s 3) K 3(S1,S,s3 )R(S1)R(S2)R1(S)
+ G1(S2 +S3)Rl(s3)
lKl(s1)Gl(s1)][lKl(s2)G(s2)]
+ 2K2(s1 82+s )R1(s1)R2(2,s3)
2R1(S2+s3)R(s3) )
+ 1lKl(s)G(sl) K2(s2's 3)Rl(s2)
K1(s2+s 3)G1( (514c)2+s)
+ s1),(, (. J j (514c)
1K l(s2)Gl(s2)
from which it can be seen immediately that the linear operator R
1
gives a transform identical to the standard transform obtained from
linear system theory.
The output power for a reactor system is given as a sum of the
operator transforms where each is multiplied the proper number of times
by transforms of the input function. This is easily derived by
applying multidimensional transforms to Equation (41). The output so
obtained has a theoretically infinite number of variables but these
can be reduced to a single transform variable by George's association
method already noted. Then the output power can be given by:
Ap(s) = R1(s)ke(s) + R2(su,U)e(S)ke(u)du
+() R(s'1 ue(Sle(ul e(2)dul
+ a...... ..
(515)
where Ap(s) and ke(s) are the onedimensional Laplace transforms of
5p(t) and ke(t) respectively. Substitution of the component trans
forms of R into this equation and association of the variables lead
to a form which can be inverted to the time domain by standard Laplace
transform techniques to obtain the time response of the reactor power.
Linear Feedback
One of the more important reactor cases which can be reduced
in the manner given above is the case of linear feedback. As indicated
earlier, this does not imply that the reactor system becomes linear,
since that quality depends a great deal on the magnitude of the input
function to the reactor. The final solution for the reactor power,
Equation (515), still contains all the higher order components of
the reactor operator but they are simplified considerably by the
absence of components of K which have orders greater than K1. The
overall feedback operator is based on the kernel given in Equation
(28) which has as its transform K1(s).
In connection with linear feedback and the neutronics model
already used, it should be noted that the output function of the
neutronics operator is a normalized power. However, the feedback
operation of Equation (28) is based on nonnormalized power so that
an adjustment of the operator is required. In the transform domain
this readily gives:
K P K (516)
ol
where P is the initial power level of the reactor and is considered
to be constant.
With the feedback operator given in Equation (516) and the
neutronics operator given in Equation (513), the transforms of the
components of the overall operator for the reactor with linear feed
back become:
Gl(s)
R1(s) = (517a)
1PoKl(s)G1(s)
Rl(sl+S2)R1(s2)
R2(sl,S2) = Kl(,), ) (517b)
1P K1(1)G1(s1)
R1(s +s2+s3)R2(s2,s3)
R (sSp,s ) = [1
R3( 1s2s 3) = 1PoK(sl)Gl(sl)
+ PoK(s2+S3)Gl(s2+s3)] (517c)
where, again, the higher order transforms are given in terms of the
lower order ones for ease of presentation only.
When these transforms are used in Equation (515) and the
result is inverted to obtain the time response, the problem becomes
identical to the one examined in Chapter II: the time response of a
reactor with pointmodel neutronics and linear feedback. Thus, a
means of comparison of the two types of series is available.
A Heuristic Argument for Stability Conditions
The transforms of the operators and the final expression for
the power of a reactor system as a function of the change in multipli
cation factor are in forms which make them amenable for use in the
determination of some stability properties of the system. In this
section some comment is made about possible conditions necessary for
stability but no rigorous mathematical proofs are presented to
accompany the statements.
From the nature of the transforms of the components of R, it
can be shown that the sum of terms making up the final transform of
the reactor power, Equation (515), converges or, in other words, the
reactor is stable when each individual term in the sum converges.
This is particularly significant because the first term in the sum is
the standard linear term which can be derived equally well from other
methods. Therefore, linear stability isa necessary condition even
before other terms are examined. This condition is not restricted to
just the pointmodel neutronics used in obtaining the transforms in
Equation (514) but is true for all cases where the neutronics opera
tor is derivable from kernels in a Volterra series of the type indi
cated in Equation (41). The analyses for linear stability are well
established and a conclusion of instability for a system obviates
further examination of the higher order transforms.
For a linearly stable system, one for which there are no poles
of R1(s) in the right half of the splane, a question still arises
regarding the convergence of the higher order transforms. No general
method exists for analyzing the operator transforms of the type pre
sented in Equation (514). However, if only linear feedback is
assumed, the transforms in Equation (517), which are considerably
simpler to handle, can be used in connection with Equation (515) to
yield some stability information for inputs which are square integrable,
i.e.:
f k(t)dt < c (518)
o
A special class of inputs, those which are bounded but not square
integrable, are examined briefly in the next chapter.
Equation (515) now becomes:
SR(u)k (u)k (su)
Ap(s) = R1(s) k (s) + 7f P Kl ( eU)1( du
Ie 2J 1P Kl(su)Gl(su)
1+PKI (u )G () 
+ ( )2f R(ul) [1 su (sI ke (sUl)dU l
I P l(ul'2)Gl(ulu)
Xf 1 l(2) (lu2)G (g) du2 + ....... (519)
where Rl(s) is defined in Equation (517a).
At this point it is assumed that all transforms can be repre
sented as ratios of polynomials or, even simpler, as ratios of pro
ducts of linear factors. Generally, if a transform does not exist in
this ratio form, it can be converted to it by means of an orthogonal
ization process in the frequency domain,(21) This is equivalent to
expanding the timedomain kernel in an orthogonal series of exponen
tials, either real or complex or both. For the transforms involved
in Equation (519), the polynomial representations are:
J K
G (s) = T((s+c )/7 (s+dk)
jl k1l
M N
Kl(s) = T (s+em)/ T(s+fn)
mEl n=l
and the denominator of Equation (517a) is:
1PoKl(s)G1(s)
I P
= 7T(s+a )/7T(s+b )
i=l pl P
P K N
7f(s+b) 7~ (s+dk) 7(s+fn)
p=1 k=l nl
and the values of all the ai are a function
readily seen that:
of P It is also
o
I+PoK(s)GI(s) = 2[1PoK(s)GC(s)] (520e)
In all cases the constants ai, bi, etc.,may be complex. The estab
lished condition for linear stability is:
Re ai 1 0
for all i.
(521)
The input is also assumed to be composed of exponentials and
can be written, in the transform domian, as:
() A
e l qs+z
qPl q
(522)
J
M
(520a)
(520b)
where
I=P
(520c)
(520d)
67
and in the limiting case of extremely slow exponential decay:
q
e() L (523)
q=l
When the quantities in Equation (520) and the input of
Equation (523) are substituted into Equation (519), association of
variables is performed according to the formulas of Appendix F, and
the partial fraction method is used, the resulting equation becomes:
Bi C
Ap(s) s + s+ jka
il j=l k=l s+a+a
I I I I K
s+a ++a +a s+a +dk *
j=l k=l ml m j=1 k=1 J k
(524)
where the Bi,Cjk, Djkm and Ejk are all constants. The limiting case
of the input function as given in Equation (523) will not contribute
anything during the process of association of the variables.
Since linear stability must be guaranteed for the system to
be stable, Equation (521) also implies that the first three terms on
the righthand side of'Equation (524) will not have any poles in the
right half of the splane. However, for the fourth term not to have
any poles in that half plane the following must be true:
Re (aj+dk) s 0 for all j and k.
(525)
This places an added restriction on the neutronics operator but does
not necessarily imply that all the values of dk must have Re dk 0.
Similar consideration of higher order terms could lead to
further restrictions but the conditions given in Equations (521) and
(525) are adequate to insure stability of system response based on
kernels up to third order.
Summary
In the preceding two chapters an examination has taken place
of a number of facets of nonlinear operator theory and its general
application to nuclear reactors. In Chapter IV the algebra of oper
ators has been presented along with multidimensional transforms which
are a means of converting the abstract operators to calculable quanti
ties based on Volterra series expansions. In Chapter V block
diagram manipulation through use of operators is demonstrated with a
reactor system as the example.
With the combined use of the block diagram manipulation
and the techniques of the former chapter, the time response of the
system is readily obtained as a sum of terms in the transform domain.
Since the conditions selected for the reactor example in this chapter
are identical to those presented in Chapter II, the sums obtained by
the two methods, when the operator method sum is converted to the
time domain, should be equivalent. It is possible to show that, at
least linearly, they are identical.
From the operator method the linear response of the reactor
system is given by the first term on the righthand side of Equation
(515) with Rl(s) given by Equation (517a). The linear portion of a
timeresponse produced by the iterative technique is relatively simple
to determine for any order of approximation. The first two terms of
this type appear in Equation (211). For comparison with the operator
method, Laplace transformation of the linear terms present in the
final infinite series which is found by the method in Chapter II
gives the series:
Aplin(s) [Gl(s) + PoGl(s)Kl(s) + P2G,(s)Kl(s) + ...]ke(s)
(526)
where Gl(s) is the Laplace transform of g(t).
Then, if
PoKl(s)G1(s) < 1 (527)
this sum can be put into the closed form:
Gl(s)ke(s)
Plin() 1P oK(s)G(s) (28)
which is identical to the form obtained directly by operators. It
should be noted that the system must be stable before the Volterra
series will converge to the closed form. The operator method gives
the same closed form but the only requirement is the existence of an
inverse operator. This is less restrictive than the inequality in
Equation (527).
In the next chapter some examples using the operator method
are given and some interesting points are noted.
CHAPTER VI
EXAMPLES AND RESULTS
Introduction
In the previous two chapters a method of system analysis
using the theory of nonlinear operators has been developed and applied
to a general reactor system. In addition, some conditions necessary
for reactor stability have been determined for feedbacks which can be
represented by strictly linear operators.
Some examples involving the use of this method to analyze
specific reactor systems are given in this chapter. The results
obtained from operator analysis are compared to results obtained
in other ways.
Reactor with Negligible Delayed Neutron Contribution
The first example is based on the reactor system in Figure 14
where the neutronics equation in the forward loop is given in the
time domain and the feedback is shown in the transform domain. The
input function is "reactivity" which is often defined in simple
reactor models as:
k 1
e (61)
P = k
e
where ke is the effective multiplication factor for the reactor. The
total input to the neutronics operator is assumed to be made up of
an external reactivity, pe, and a reactivity, pf, which is the output
of the feedback operator.
For this particular example it is assumed that there are no
delayed neutrons contributing to the total neutron population. The
elementary reactor kinetics equation is then:
dn pn
dt (62)
where n is the neutron population, a quantity which determines the
power of the reactor, and A is the neutron generation time. The
method used to convert Equation (22) to Equation (23) is used here
to obtain kernels of the type appearing in Equation (512). From
these kernels the corresponding transforms of Equation (513) can be
derived. These transforms are determined to be:
Gl(s) (63a)
As
1 1 1
G2(s1,s2) A (sl+s2) 2 (63b)
G (s1,s,s ) 1 1 (63c)
G3( ) A3 (s1+s2+s3) ( 2+s3) s3
for the operators of the first, second, and third order.
The feedback operation is linear and is given by:
n
K (s) = o (64)
(1+s)2
whereas described previously, the initial neutron population must be
introduced to adjust for the fact that the neutronics operator yields
a normalized neutron population as output but under ordinary conditions
the feedback operator requires a nonnormalized neutron population as
input. With the transforms of the operators G and K, given in
Equations (63a) and (64) respectively, the transform of the linear
operator becomes:
s +2s+l
Rl(s) = 2 (65)
( s3+2s2+s+no)
A
The transforms of the higher order components of R are also easily
obtained once G and K are known. These transforms and the transform
1 1
of the input function are used in Equation (515) to obtain the
response of the reactor system.
For the moment, the stability of this particular reactor
system will be considered. It is immediately obvious from Gl(s), in
this case, that the stability criterion in Equation (525) reduces to
the same criterion as Equation (521) so that only linear stability
need be established to guarantee stability of the reactor output due
to kernels of, at least, up to third order. With any of several well
known methods (Nyquist, Rootlocus, etc.) which can be used to analyze
the transform of Equation (65), it is found that the system of
Figure 14 is stable for no/A on the closed interval Ono/A <2 for
bounded inputs. However, one of the conditions required in order to
establish stability criteria for the nonlinear system in the previous
chapter is that the input to the system be square integrable. This
Figure 14.
A Reactor System without Delayed Neutrons
is slightly more restrictive than the purely linear theory condition
but more accurate in the nonlinear case.
In order to test the stability criterion and also to provide
a means of comparison for the operator method calculations, the same
system has been analyzed by analog simulation. Actually a digital
computer code entitled MIMIC,(20 which does the analog simulation
numerically, has been used to calculate the "true" response of the
system at various normalized initial neutron population levels for a
variety of input functions. This code has been used successfully at
the University of Florida for more complex simulation problems.(23)
Figure 15 shows the "true" response to an exponentially
decaying input which is large enough in magnitude, ke = .05 exp(3t)
to cause the system to act nonlinearly but which is also square
integrable. It can be seen that the system is stable for no/A on
the closed interval ONno/A !2 but definitely unstable for no/A >2.
To be compared to this is Figure 16 which shows the "true" response of
the system to a stepfunction input. The step is a bounded function
but is not square integrable. In this figure it is seen that the
system is stable on the open interval On o/A <2 but unstable at the
Transition point n /A =2.
In addition to demonstrating the accuracy of stability criteria
obtained, even if only obtained heuristically, from operator theory,.
it is desirable to also show the accuracy of the time response of
a system as predicted by the operator method. In order to obtain the
time response for the system shown in Figure 14, a value of no/A 2
4xlo' P1i U.u: expk(.5rj
a2 2.0
2
~x1O I I r
j :x
8n In I2 I n
o A A 1.9
:1:1
:1I '
:1I.
2 :/ \ I
2x10
4x102 2 21
4 8 12 16 20 24 28 32 36
Time(sec.)
Figure 15. "True" Response of System with No Delayed Neutrons Exponential Input
1.6x1o"1
1.2x10 / I \
I iI I I
I I \
I I II I
I i I I I
I 2
4x0o2 =
p, = 0.05
SI I I I I I I
4 8 12 16 20 24 28 32 36
Time(sec.)
Figure 16. "True" Response of System with No Delayed Neutrons Step Function Input
Figure 16. "True" Response of System with No Delayed Neutrons Step Function Input
has been selected along with a stepfunction input of pe = 0.05.
Since the feedback is linear, Equation (519) can be used for the
calculation with Rl(s) given by Equation (65), Gl(s) I/A s, and
Kl(s) given by Equation (64).
Figure 17 shows the result of taking the inverse Laplace
transform of the linear terms, the secondorder terms, and the whole
sum in Equation (519). It is to be noted that the linear response
neither converges nor diverges but is oscillatory as linear theory
predicts for the bounded input and the initial neutron population at
the transition point. However, the secondorder term is definitely
divergent and it is this term and higher order terms which cause the
imposition of the additional restriction of square integrability on
the input. The secondorder term is dependent on the square of the
magnitude of the input so that it can reasonably be expected to play
a less significant part as the magnitude of the input function is
decreased. At some point, the linear term dominates sufficiently to
cause the system to appear to be a purely linear one.
The total response curve in Figure 17 and the "true" response
curve for no/A = 2 in Figure 16 show identical time responses. Based
on the numerical results, the two curves exhibit less than one percent
difference over most of the time scale shown except when the time
response is very small in magnitude (e.g. at t18.5 seconds). Some
comparative data for the two time response curves are tabulated in
Appendix G.
2.0x101
4 8 12 16 20 24 28 32 36
Time(sec.)
Figure 17.
Operator Method Response of System with No Delayed Neutrons
Step Function Input
TwoTemperature Feedback Model
The second reactor system considered is the twotemperature
feedback model and numerical constants presented in Chapter II and
Appendix C. The necessary kernel transforms of G1 and K1 are given
by the quantities in brackets on the righthand side of Equations (219)
and (220) respectively. Here too, the feedback is linear so that
Equation (519) can be used to calculate the sum from which the time
response can be obtained.
A look at Gl(s), which is given by:
Gl(S) S+ (66)
1*s(s+3/1*)
reveals that the dk, as defined in Equation (520d), satisfy the
relation:
dk 0 for k = 1,2 (67)
so that the stability criterion of Equation (525) will be true when
ever Equation (521) is true. For positive feedback from both tempera
ture loops the system is always unstable. Otherwise the system sta
bility depends on the sign of the feedback from each loop and the
relative magnitudes of the temperature coefficients. The values for
the various constants in Chapter II have been selected to insure
stability of the model.
Figure 18 shows the result of transforming the linear terms,
the secondorder terms, and finally the whole sum in Equation (519)
to the time domain when the input to the system is a stepfunction of
10x103
6xl03
2xl03
0
2xl03
12
Time(sec.)
Figure 18.
Operator Method Response of TwoTemperature
Feedback Model Step Function Input
k = 0.05 and the system constants are given in Table 2 of Chapter II.
In this case, it can be seen that the secondorder term is significant
for times up to about nine seconds but the system is essentially
linear for times greater than that. Since the linear response can be
checked by an independent calculation, it is assumed that the repre
sentation of the total response, which is very close to linear for
much of the time shown, is accurate. This is supported strongly by the
accuracy of the)operator method calculation in the preceding example.
For comparative purposes, the time response for the same model
has been calculated from the first few terms of the Volterra series as
given in Equation (211) with a stepfunction input given by ke = 0.05.
The result of the calculation is shown in Figure 19 which is different
from the results shown in Figure 6 because the magnitude of the input
is now larger (k = 0.05 compared to k = 0.01) and the input is a
stepfunction rather than a 5function. Since the system is known to
be stable, it is obvious from a comparison of Figures 18 and 19 that
many more terms are required in the Volterra series than have been
calculated in order that meaningful results may be obtained.
Summary
Based on the results from the two sample reactor models pre
sented in this chapter, several points are particularly to be noted.
First, the stability criteria derived from the operator method for
systems with linear feedback have been found to be valid for the sys
tems examined and the types of inputs to them. The criteria are also
seen to be relatively simple to apply.
1x105
2x105
3xl05
4 8 12 16 20
Time(sec.)
Figure 19.
Response of TwoTemperature Feedback Model
Step Input Volterra Series
Secondly, the magnitudes of the inputs in both cases have been
chosen sufficiently large so that the systems will respond nonlinearly
and linear analysis will definitely be in error. For the input
selected, the first and secondorder terms are found to be adequate to
give good agreement with the true response of the systems, with the
secondorder term contributing 25 percent or more of the total response.
In the system with no delayed neutrons, the first example, it is
probable that the thirdorder term will also contribute a nontrivial
response for the times greater than thirty seconds and will help
increase the accuracy of the calculated response. Generally, the
accuracy of the calculation is good with only the first twoorder
terms but it is definitely not adequate with ohly a linear response
present unless the input is greatly reduced. For ke 0.01 the
secondorder term would only contribute about 5 percent of the total
response as opposed to the 25 percent cited above for ke = 0.05.
Finally, in the comparison of the operator method with direct
Volterra series calculation, it is obvious that the calculations
based on operator theory give more accurate results for less calcula
tional effort. The Volterra series requires the calculation of many
more terms than is physically practical in order to obtain reasonable
results.
CHAPTER VII
SUMMARY AND CONCLUSIONS
In the work that has preceded, two methods of analysis of
the timedependent characteristics of nuclear systems have been
developed. Both methods yield system response in the form of an
infinite series of terms but the type of series obtained is differ
ent for the two cases.
The first method is the Volterra series expansion obtained
by iteration on an integral equation description of the system and
quantities are dealtwith strictly in the time domain. The second
method is based on the theory of nonlinear operators and draws on the
Volterra series to relate the operators to physical quantities. Then
the series resulting from operator manipulation is in the transform
domain. However, techniques are examined for converting this series
into a time response.
From the theory and results presented in this work, a number
of conclusions can be drawn about the characteristics of each method
and about the comparative abilities of the two methods. First, the
Volterra series solution is not practical for calculating timeresponse
characteristics of systems because it requires the calculation of
many more terms than is physically feasible to produce results of
only fair accuracy. Also, the conditions required for reduction of
the Volterra series to a series of terms, each term being only of a
single order, are overrestrictive. Consider the classic example:
1 1 + x + x2 + ... (71)
1x
where the ratio on the left exists and is finite for all x except
x = 1, but the series on the right converges only for Ixl < 1. The
Volterra series exhibits the same type of behavior.
On the other hand, based on the Volterra expansion, it can
be concluded and has been proven that, for a particular neutronics
model, reactor systems with positive feedback definitely can be
stable over some range of power levels. In addition, the Volterra
method provides a theoretical basis for relating the system operators
of the second method to quantities which are physically meaningful.
Next, considering the operator method, the accuracy of time
response obtained with this method is very good for a calculated
response made up of just firstand secondorder terms. The small
number of terms which need be calculated helps point up the superi
ority of this method over the direct, timedomain, Volterra series.
The stability criteria derived from the operator method for linear
feedback systems are slightly more restrictive on both the input to
the system and the forward loop operation than the standard criteria
of linear system theory.
Finally, the operator method can be used to convert more
complex nonlinear systems to a series of calculable kernel trans
forms by block diagram manipulation if the component operations of
the system are each expressible in Volterra series form. This is the
case for a wide variety of systems.
Concluding Comments
From the theory of random signals, a theoretical method for
measuring the various order kernels in a Volterra series of the type
used with the operator analysis has been developed by Wiener.(24)
An extension of this work in order to put it on a more practical
level has been carried out by Hooper and Gyftopoulos (9,25) who base
their proposed measurements on a pseudorandom ternary input to the
system under analysis. The proper higher order crosscorrelations
between the output and input then give the desired kernels in the
time domain. Thus, all the "true" kernels for a reactor system can
be experimentally determined.
This dissertation gives an analytic method.of determining
these same kernels. In Chapter VI it has been shown that, for a
particular reactor model, the operator method can give very accurate
predictions of the time response of a system. This implies that
there must be accurate prediction of the system kernels for the model
chosen. However, the models for the neutronics and feedback opera
tors are only approximate and the calculations of "true" time
response can be only as good as the models chosen to represent the
components of the total system.
A comparison of the kernels calculated by the operator method
with the kernels determined experimentally by some method like the
one mentioned above would give a good indication of the accuracy of
87
the analytic models used in reactor analysis. This is one of the
more important uses of the work developed here.
CHAPTER VIII
RECOMMENDATIONS
It is widely observed in technical work that the attempted
solution of one type of problem always leads to the uncovering of a
myriad of associated problems. So too, the investigation 6f series
solutions for reactor systems undertaken in this work, particularly
the solutions obtained through the use of operators, has led to some
questions which should prove interesting to investigate.
The operator method has been shown to successfully yield
accurate results for systems with linear feedback. However, the
method has been developed for a general feedback,so a natural question
arises as to the kind of results which would come from a system with
nonlinear feedback. Certainly different stability criteria would be
derived and these would probably be more dependent on the feedback
function than in the linear case.
The neutronics operator used in the examples has been based
on a point model. A more complex core neutronics model should give
better comparisons to actual reactor systems.
Some of the other possible reactor applications for which the
operator method holds some promise involve long term effects such as
xenon tilt. Also, in nuclear rockets, the Volterra series in
89
combination with block diagram manipulation should yield some prelim
inary stability information. Similarly, the method might hold promise
for coupledcore systems.
There are many other applications of the series approach both
in the reactor field and in a large number of other fields.
APPENDICES
APPENDIX A
Region of Applicability of ZeroPower Kinetics
Start with Equation (214):
t
8po(t)5p1(t) g(t T)ke( ) g(T)ke(o)dadT
<< 1 (A1)
 p (t) M < 
Po(t) f g(tT)ke(T)dT
0
but, for a physically realizable g(t):
t T t
Sg(tT)k (T) g(To)ke(a)dadT [f g(tT)k (r)dv]2. (A2)
0 0 0
Therefore
t
86p(t)6p(t) g(tT)k()d2 (A3)
"Po) t, << 1 (A3)
i (g(tT)ke(T)d
0
or
Spo(t)6pl(t) t
5Po (t) f g(tT)ke(T)dTv 1 (A4)
Now, applying Equation (32)
t t
If g(tT)k (T)dTj ~ f g(teT)I Ike(T)I d (A5)
0
