Title Page
 Table of Contents
 List of Tables
 List of Figures
 Volterra series
 Maximum value series
 Nonlinear operations - theory
 Nonlinear operations - applica...
 Examples and results
 Summary and conclusions
 Biographical sketch

Group Title: application of Volterra series and nonlinear operators to nuclear reactor kinetics.
Title: The application of Volterra series and nonlinear operators to nuclear reactor kinetics.
Full Citation
Permanent Link: http://ufdc.ufl.edu/UF00089980/00001
 Material Information
Title: The application of Volterra series and nonlinear operators to nuclear reactor kinetics.
Physical Description: Book
Language: English
Creator: Thierer, Ira
Publisher: Ira Thierer
Publication Date: 1967
 Record Information
Bibliographic ID: UF00089980
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
Resource Identifier: alephbibnum - 000574301
oclc - 13847074

Table of Contents
    Title Page
        Page i
        Page ii
    Table of Contents
        Page iii
        Page iv
    List of Tables
        Page v
    List of Figures
        Page vi
        Page vii
        Page viii
        Page ix
        Page 1
        Page 2
        Page 3
        Page 4
        Page 5
        Page 6
    Volterra series
        Page 7
        Page 8
        Page 9
        Page 10
        Page 11
        Page 12
        Page 13
        Page 14
        Page 15
        Page 16
        Page 17
        Page 18
        Page 19
        Page 20
        Page 21
        Page 22
        Page 23
        Page 24
        Page 25
        Page 26
        Page 27
        Page 28
    Maximum value series
        Page 29
        Page 30
        Page 31
        Page 32
        Page 33
        Page 34
        Page 35
        Page 36
        Page 37
        Page 38
        Page 39
    Nonlinear operations - theory
        Page 40
        Page 41
        Page 42
        Page 43
        Page 44
        Page 45
        Page 46
        Page 47
        Page 48
        Page 49
        Page 50
        Page 51
        Page 52
    Nonlinear operations - application
        Page 53
        Page 54
        Page 55
        Page 56
        Page 57
        Page 58
        Page 59
        Page 60
        Page 61
        Page 62
        Page 63
        Page 64
        Page 65
        Page 66
        Page 67
        Page 68
        Page 69
    Examples and results
        Page 70
        Page 71
        Page 72
        Page 73
        Page 74
        Page 75
        Page 76
        Page 77
        Page 78
        Page 79
        Page 80
        Page 81
        Page 82
        Page 83
    Summary and conclusions
        Page 84
        Page 85
        Page 86
        Page 87
        Page 88
        Page 89
        Page 90
        Page 91
        Page 92
        Page 93
        Page 94
        Page 95
        Page 96
        Page 97
        Page 98
        Page 99
        Page 100
        Page 101
        Page 102
        Page 103
        Page 104
        Page 105
    Biographical sketch
        Page 106
        Page 107
Full Text







April, 1967


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


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.






ABSTRACT . . . . .


















S . . . .. . . .*

INTRODUCTION . . . . . .








. . . . . .

. . . . . .





COMBINATIONS . . . . . . . . .

ASSOCIATION FORMULAS . . . . . . . . .

"TRUE" RESPONSE . . . . . . . . .




















. .

. . .

LIST OF REFERENCES . . . . . . . . . . 104

BIOGRAPHICAL SKETCH. . . . . . . . . . 106


Table Page

1. Constants for the Zero-Power Kernel . . . . . 20

2. Two-Temperature Feedback Model Data . . . . . 24













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 Two-Temperature
Feedback Problem . . . . ....

Linear Response Two-Temperature Feedback

Region of Convergence of the Maximum Value

Regions of Absolute Stability EBWR Model

General Nonlinear System . . ...

11. Combination of Operators .

* .* . .

* .

. .

* .



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 Two-Temperature
Feedback Model Step Function Input .... ...
















. . .



19. Response of Two-Temperature Feedback Model -
Step Input Volterra Series . . . . 82

C-1. Block Diagram of Two-Temperature Feedback Model . . 94

Abstract of Dissertation Presented to the Graduate Council
in Partial Fulfillment of the Requirements for
the Degree of Doctor of Philosophy



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


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


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.



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

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


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 odd-order kernels. Thus results for systems with relatively large

even-order 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 (1-1)

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

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

time-dependent 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 time-response

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


absolutely necessary. Then linear feedback models have been chosen

because of ease of calculation. For the neutronics models, it appears

that space-time 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.




One of the more easily applied and potentially useful approaches

to the solution of integral or integro-differential 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 point-model 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 point-model. 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

space-time separability, the presence of only monoenergetic neutrons,

and the applicability of diffusion theory, the reduction of the Boltz-

mann equation to the point-model kinetics equations is well established,

These equations, in terms of power rather than neutron population,


dP k(+-S)-1 N
dt 1 i C

dCi kpi
d--- = I---P- XiCi
dt Ii* i

where P. = the reactor power level,

C. = the concentration of the ith delayed neutron group

k = the effective multiplication factor,

i* = the neutron mean lifetime,

Xi the decay constant for the ith group of delayed neutrons,

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 steady-state 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 (2-1)

and appropriate elimination of some of the steady-state components


dSP tko+6k(t)](l-p)-l N
S 15P (t) + -P P8k(t) + X iC(t)
dt 1* 1 iCi(t)
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 steady-state component and the

prefix 5 indicates the transient component of the variable to which

it is attached. The Laplace transforms are taken of Equations (2-2)

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(t-T)fk(T)dT + g(t-T)6k(T)P(T)dT (2-3)
o fo
where g(t) = the inverse Laplace transform of Gl(s), and

1-s Z P

Gl(s) = i= N L
1-k +s s*+k Z
o o s+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 ..,t-T;P(rl)...P(Tn)]
o o (2-4)

where K is some feedback mechanism kernel. In most cases of interest

in this work the very general expression can be reduced to:


kf(t) = K[t-T;P(T)]dT

which can still take into account nonlinearity in the feedback mecha-


Then, if the externally introduced variations in the multipli-

cation factor are denoted by ke(t) the total variation is given by:

5k(t) = ke(t) + K[t-T;P(T)]dT. (2-6)

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 (2-3) results in:

t t T
5P(t) = P g(t-T)ke(T)dT + Po g(t-T)dTf K[T-a;P(a)]da

+ g(t-T)5P(T)ke(T)dT

t 7I
+ f g(t-T)5P(T)dT 7T-oa;P(o)]da. (2-7)
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 (2-7) a model

for the kernel R(t-T;8P(T)] must be chosen. Linear feedback has been

initially selected for the sake of simplicity. The kernel then becomes:

K[t-T;6P(T)] = Kl(t-T)5P(T) (2-8)

where Kl(t) is the linear feedback kernel. Introduction of this feed-

back model into Equation (2-7) yields:
t t T
8p(t) = g(t-T)ke(T)dT + g(t-T)dT Kl(T-r)Po6p(a)da

+ g(t-T)6p(T)ke(T)dT

t T
+ g(t-T)5p(T)dTf Kl(T-a)P6p(a)do (2-9)

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 (2-9) belongs is found, according to Volterra,(1l) by

iteration. An initial solution, 5Po(t), is assumed and substituted

into the right-hand side of the equation to give a first approximation

to the true solution. This approximation is then substituted into the

right-hand 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 (2-9) the initially assumed solution is:

8P((t) = g(t-T)k (,)dT (2-10)

which is equivalent to zero-power reactor kinetics. The first approxi-

mation derived from this is:

pl(t) = g(t-T)k(T)dT + g(t-T)g(Tro)ke(T)ke())doad

+ Pofff g(t-T)g(o-p)K1(T-or)ke(p)dpdad

+ Po f g(t-)g(T-p)g(-g)Kl( )ke(p)k()
dpdadpdT ,

and does not have as simple an interpretation as Equation (2-10)

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 (2-12)

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 ever-increasing 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 (2-9), 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


p(t) <


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

(zero-power kinetics) is a good representation of the total solution.

Substitution of Equations (2-10) and (2-11) into (2-13), combined with

the additional assumption that feedback is negligible, i.e. Kl(t) = 0

for all t, results in:

t T
Itg(t-T)ke (T) g(T-o)k (cr)dadT
---<-- 1 i. (2-14)
g(t-T)ke (T)dT

This yields the following condition on ke(t):


oIk ( )|dT <<
e I-max I

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 (2-3). 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:


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 (2-13) 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 (2-15), 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



Two examples have been chosen to demonstrate the effectiveness

of finding time-dependent 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 e-2.63t + 0.2685 e"3.12t

+4.62XO1-4 e-4t 0.0584 e-572t

+0.0514 e"9'09t 8.62X10"3 e"33"33t (2-16)

where Kl(t) is given in units of cents/KW-sec. 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-
cally approaches a constant value of -3.50X10 cents/KW'sec.

The zero-power kinetics kernel to go with this model is the

same one as defined in connection with Equation (2-3). With the con-

stants given in Table I, this kernel becomes:

g(t) = 11.00 + 46.78 e-0O1847t + 9872.5 e-701013t. (217)

For a 8-function 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).


KW.sec -.012




o 2 4 6 8 10
Figure 2. Feedback Function Based on EBWR Data





0 5 10 15 20 25 30 35 40 45 50

Figure 3. First Approximation in the Volterra Series


Constants:for-the Zero-Power Kernel

PI 0.002074

P2 0.004896

X 0.02524 sec"-

X2 0.56409 sec-1

10-4 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( L-1 G1(s)k (S)
lin(t) L PoKl(s)G(s (2-18)

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

(2-16) and (2-17). 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




5 10 15

Figure 4.

k = .001 40

P 104 KI






SI I I I I-10
20 25 30 35 40 45 50

First Iterative Solution EBWR Problem

Scale Scale
for for
k =0.01 k =0.001
e e
8 -\ 0.8
P = 10 Kw




Envelope of Response Curve

-8U -0.8

15 20 25 3
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 two-temperature

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


p is the fraction of neutrons appearing as delayed

neutrons, and


1* is theimean neutron lifetime.

From the model and Appendix C, the feedback kernel is:

Kl(t) = L-i [ (l+T (2-20)

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 6-functions, and with

the data given in Table 2,. Again, for comparison, the linear response

of the system has been found from Equation (2-18). This is shown in

Figure 7 with the same scales as indicated for Figure 5.


Two-Temperature Feedback Model Data

0.125 sec-1
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


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





0.2 0.4 0.6 0.8 1.0 1.2 1.4


Figure 6. First

Iterative Solution Two-Temperature
Feedback Model







k =0.01



Figure 7. Linear Response Two-Temperature Feedback Model

k =0.001





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 two-temperature 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 two-temperature feedback case. Calculation of time

responses for the next iteration, 5p2(t), with twenty-five terms is

almost prohibitive.


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




In the previous chapter it has been noted that the infinite

Volterra series solution for the time-dependent 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


Development of the Maximum Value Series

The maximum value series is obtained from the basic inequality

given in Equation(3-2) below. If:

Ir f f(x)g(x)dx (3-1)


b b
III b f(x)g(x)dxl J If(x)I |g(x)ldx (3-2)
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 (3-3)
a a
because gmax is a constant.

Now the application of Equation (3-3) to the initial approxi-
mation for the Volterra series, Equation (2-10), results in:

IPo(t)I < kM ft |g(t-T)|dT (3-4)

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(t-T)JdT a G
o o
f IK(T)IdT f IJKl(t-T)IdT a K. (3-5)
o o

1The quantities defined in Equation (3-5) 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 right-hand

side of Equation (2-10) becomes:

t 7

1121 = I g(t-T)ke(T) f g(T-C)ke(a)daddTl
0 0

2 t t
S, k2 I g(t-T)| f g(T-0)|dadT
o o

: k[F Ig,-lT12 2 2 (3-6)

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

and the second iteration is immediately seen to give:

ISl1 kMG + G+ o2 o 2K + P kG2K. (3-8)

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) (3-9)
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:

1Plinl -P KG (3-11)

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



Equation (3-9) 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. (3-12)
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 (3-9) 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 right-hand side of Equation

(3-12), 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. (3-13)

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 (3-13) 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 e-folding 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 (3-13). 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


Figure 8. Region of Convergence of the Maximum Value Series

Equation (3-9) is not in a form amenable to calculation, it can be

simplified somewhat. In Equations (3-10) 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 (3-13) are non-negative, the con-

dition P KG<1 needed to obtain the simplification will always exist

if Equation (3-13) is satisfied. For terms up to fourth-order, the

simplifying process then gives:

kMG (kMG)2 (1+P KG)
-P KG (1-P KG) + (1-PoKG)5 (kMG)3

1+3P KG+(P KG)2
+ +3P0KG+(P KG (kMG)4 + ... (3-14)
(1-P KG)7

which is a much better form for calculation of the maximum J1pl.


To demonstrate the conditions derived above, the EBWR model

used in Chapter II has been chosen. Equations (2-16) and (2-17) 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 found-which will keep the system absolutely

T 5 sec.

T = 10 sec.

T = 20 sec.

T = 60 sec.

T 300 sec.


Figure 9. Regions of Absolute Stability EBWR Model







stable for any particular kM. This is the reason for the selection of

the axes used in Figure 9.

From Equation (2-17) it can be seen that the value of G, as

defined in Equation (3-5), 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 (3-14) yields:

18pl < 0.305 + 0.1162 + 0.0533 + 0.0277

+ 0.0157 + 0.00955 ..... (3-15)

where each quantity corresponds to a particular order term in Equation

(3-14) 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 (3-13) 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 (3-15) 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 (1-1) is definitely over-restrictive.

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 all-positive 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, non-maximum

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.




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


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

so-called "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

sets-of 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)

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:

y(t) j hl(i)x(t-T)dT + hT1l, 2)x(t-T1)x(t-T.)dTdT2
-CO -00 -CO

+ ..... + .... ( h ,. )xt(tT,) ...x(t-T

n' n

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


where H f hl(r) -- dT,

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) + .... (4-3)

which can be rewritten as:

y(t) (i + 2H + 3 + ..... +--H + .... [(t)] (4-4)

or, upon combining all the operators into a single operator for the


y(t) = H[x(t)] (4-5)

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

operator, the output due to two simultaneous inputs is:

Y2(t) = H[x(t)'+ z(t)] (4-6)

which, with the use of the definition of H2 given in Equation (4-2),


Y2(t) = f h2 (1 '2)[x(t-T1) + z(t-T1)]rx(t-T2) + z(t-rT2)]

dT ldr2 (4-7)

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) (4-8)

where H2(xz) ff h2(.T2)x(t-Td)z(t-T2)d1dT2., A second-order
-o? -o)
operator acting on a sum of two inputs leads to an output made up of

three second-order 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

X(t) C y (t


L() C y(t)


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 (4-9)


C = A.B (4-10)

respectively. The addition and multiplication of operators are both

commutative and associative, thereby giving the following identities:

A + B B + A (4-lla)

A + (B(C) (A+B) + C (4-11b)


A.B = B.A (4-12a)

A*(BSC) = (A*B)*C (4-12b)

The multiplication operation is also distributive:

A.(-B+C) = (A-B) + (A.C) = (B+C)*A (4-13)

where the last equality results from Equations (4-12a) and (4-9).

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


c A*_B (4-14)

and is associative:

(A*B)*c = A*(B*C) (4-15a)

but is not in general commutative:

A*B 0 B*A. (4-15b)

The most important exception here is the case where A and B are both

linear operators. Then, these operators are always commutative in


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)

(A'B)*C = (A*C).(B*CC)


C*(A+B) ( (*A) + (LCB)

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 second-order 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 (4-7) is still valid. Since the outputs

may be expressed as operators acting on the input to the various

systems, Equation (4-7) can be rewritten as:

y2(t) H2((A[w])2) + 2H2(A[w]B[w])

+ H((B[w]) ) (4-18)

where H2((A[w])2) = H*A. Equation (4-18) gives rise to the general

cascade operation denoted by "o":

(Ho(A.B))[w] = H(A[w].[Bw]) (4-19)

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 (4-20)

A*I = I*A = A. (4-21)

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 one-to-one correspondence

between equal members of the two. The remaining operations of sub-

traction and inversion are derived from Equations (4-20) and (4-21) by

the following:

A+ (-A) = 0 (4-22)

A*A1 = (4-23)

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


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 nth-order 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 (4-i). By -using the standard

Laplace transformation, the first term on the right-hand 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 one-dimensional case, by:

Fn(Sl,...,sn) f fn(tl, ..,tn)exp(-ltl-...-sntn)
0 0
dtl....dt (4-24)

a+j 0+jcO

fn (t1 tn ) = ( )n .... f Fn(s1,...,s)
-jo -ja

x exp(sltl+...+s t )ds.... ds. (4-25)

If Equation (4-24) is applied to the kernels of Equation (4-1),

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 (4-1) 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)
dTI ....dT. (4-26)

The result of applying the transform definition to this is:

Yn(sl...s ) = Hn(s ...,sn)X(sl)'..'X(s). (4-27)

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 (4-1). 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 single-variable function can

then be inverted to give a time response directly. For the second-

and third-order operators this association is given formally by:

G3j w*
H(s) 2j H(s-u,-u)du (4-28)

O-jo a-jc
duldu2 (4-29)

and the association technique is easily extendible to higher order

operators. As an example of association, consider the second-order


H2(S s2) +a (4-30)

which inverts to:

h2(tlt2) = ABexp(-at -bt2) t, t2 2 0 (4-31)

and, upon setting tI = t2 = t, becomes

h2(t,t) = ABexp[-(a+b)t] t 0. (4-32)

The transform of this is now

H(s) = s++ (433)

Thus, if H2(s1,s2) is treated formally by association, the result must

H2( I A B AB
(s) --------- du sab (4A-34)
) s 2j j s-u+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 single-dimensional one. Association

formulas based on Equation (4-28) and using the result obtained in

Equation (4-34) 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.




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

Figure 12. Reduction of a General Reactor System Part I

G =" + +G ..... +G + .... (5-1)
S 1 =2 -n

K 1 + K2 + ..... +K + .... (5-2)
1 -n

where, for any particular case, all or possibly only some of the various

order operators are non-zero. 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. (5-3)

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


Upon application of the distribution law for a cascade over a

sum of operators and with the appropriate general cascade operations,

the right-hand side of Equation (5-4) 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 (5-5a)

2 *--2 +2 (5-5b)

H3 = K1G3 + ~2(Gl'2) + K3 1 (5-5c)

and succeeding terms are found in a similar manner.

Next, the feedback loop containing H is converted to a single

feed-through 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 ]. (5-6)

Again, assuming a multi-order expansion for L, this last equation can

be rewritten as:

L1.2+L3+ ... = I + (H1+H2+H3+.... )*(L1+ +...) (5-7)

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 (5-8a)

L -1*2 + -H*L1 (5-8b)

L3 H*L + 2%eo( .L ) + H 3*L (5-80)
-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 (7-Hl)1 (5-9a)

L = (I-H,) -1*H('1-H 1 (5-9b)
-1 1-1 (5-9b)

+ H *(I-l)'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 (5-10)

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*(I-K!ll)'-1 (5-11a)

_R 2 g*(I-a 1)" 1+ Gl*(-KlGl)1* *[_Kl 2*(.1*1)"

where the quantities found in Equations (5-5) and (5-9) 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 (4-1). 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 point-model given in Chapter II has been used again here.

Equation (2-3), with a normalized power variable, is solved by the

iteration technique of Chapter II to yield the following Volterra


t t Tl
p(t) g(t-T)5k(,)dT + f f g(t-Tltg(Tl-T2)6k(T( 6k(T2)dTldT2
o 0 0

+ f ... g(tTg(tr)(Tl').. g(Trn-1 )k(T) .... k(T)
o0 0
dT n.....d .


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) (5-13a)

G2(sl,s2) G1(s+s2)Gl(s2) (5-13b)

G3(s81s2s3) = GI(s +s2+s3)G1( s2+s3)G(s3) (5-13c)

and similarly for the rest of the components.

Now, with the application of multidimensional transforms to
the operators of Equation (5-11) 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 point-model neutronics given in Chapter II, the various order
neutronics kernels which are based on this model and are given in
Equation (5-13) 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:

Rl(s) l () (5-14a)

R2(s 2 1Kl(s,)Gl(s1) + K2(sl's2)Rl( sl +s2)R1(s8)Rl(s2)


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)

+ 2K2(s1 82+s )R1(s1)R2(2,s3)

2R1(S2+s3)R(s3) )
+ 1l-Kl(s)G(sl) K2(s2's 3)Rl(s2)

K1(s2+s 3)G1( (5-14c)2+s)
+ s1),(, (. J j (5-14c)
1-K l(s2)Gl(s2)

from which it can be seen immediately that the linear operator R
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 (4-1). 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(s-u,U)e(S-)ke(u)du

+() R(s'1 -ue(Sle(ul e(2)dul

+ a...... ..


where Ap(s) and ke(s) are the one-dimensional 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 (5-15), 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

(2-8) 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 (2-8) is based on non-normalized power so that

an adjustment of the operator is required. In the transform domain

this readily gives:

K P K (5-16)

where P is the initial power level of the reactor and is considered

to be constant.

With the feedback operator given in Equation (5-16) and the

neutronics operator given in Equation (5-13), the transforms of the

components of the overall operator for the reactor with linear feed-

back become:

R1(s) = (5-17a)

R2(sl,S2) = Kl(,), ) (5-17b)
1-P K1(1)G1(s1)

R1(s +s2+s3)R2(s2,s3)
R (sSp,s ) = [1
R3( 1s2s 3) = 1-PoK(sl)Gl(sl)

+ PoK(s2+S3)Gl(s2+s3)] (5-17c)

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 (5-15) 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 point-model 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 (5-15), 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 point-model neutronics used in obtaining the transforms in

Equation (5-14) 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 (4-1). 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 s-plane, 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 (5-14). However, if only linear feedback is

assumed, the transforms in Equation (5-17), which are considerably

simpler to handle, can be used in connection with Equation (5-15) to

yield some stability information for inputs which are square integrable,


f k(t)dt < c (5-18)

A special class of inputs, those which are bounded but not square

integrable, are examined briefly in the next chapter.

Equation (5-15) now becomes:

SR(u)k (u)k (s-u)
Ap(s) = R1(s) k (s) + -7f P Kl ( eU)1( du
Ie 2J 1-P Kl(s-u)Gl(s-u)

1+PKI (u )G () -
+ ( )2f R(ul) [1 s-u (sI ke (s-Ul)dU l

I P l(ul'2)Gl(ulu)
Xf 1 l(2) (l-u2)G (g) du2 + ....... (5-19)

where Rl(s) is defined in Equation (5-17a).

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 time-domain kernel in an orthogonal series of exponen-

tials, either real or complex or both. For the transforms involved

in Equation (5-19), the polynomial representations are:

G (s) = T((s+c )/7 (s+dk)
jl k1l

Kl(s) = T (s+em)/ T(s+fn)
mEl n=l

and the denominator of Equation (5-17a) is:


= 7T(s+a )/7T(s+b )
i=l pl P

7f(s+b) 7~ (s+dk) 7(s+fn)
p=1 k=l n-l

and the values of all the ai are a function

readily seen that:

of P It is also

I+PoK(s)GI(s) = 2-[1-PoK(s)GC(s)] (5-20e)

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.


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











and in the limiting case of extremely slow exponential decay:

e() L (5-23)

When the quantities in Equation (5-20) and the input of

Equation (5-23) are substituted into Equation (5-19), 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


s+a ++a +a s+a +dk *
j=l k=l ml m j=1 k=1 J k


where the Bi,Cjk, Djkm and Ejk are all constants. The limiting case

of the input function as given in Equation (5-23) 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 (5-21) also implies that the first three terms on

the right-hand side of'Equation (5-24) will not have any poles in the

right half of the s-plane. 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.


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 (5-21) and

(5-25) are adequate to insure stability of system response based on

kernels up to third order.


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 right-hand side of Equation

(5-15) with Rl(s) given by Equation (5-17a). The linear portion of a

time-response 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 (2-11). 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)

where Gl(s) is the Laplace transform of g(t).

Then, if

PoKl(s)G1(s) < 1 (5-27)

this sum can be put into the closed form:

Plin() 1-P 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 (5-27).

In the next chapter some examples using the operator method

are given and some interesting points are noted.




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 (6-1)
P = k

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 (6-2)

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 (2-2) to Equation (2-3) is used here

to obtain kernels of the type appearing in Equation (5-12). From

these kernels the corresponding transforms of Equation (5-13) can be

derived. These transforms are determined to be:

Gl(s) (6-3a)

1 1 1
G2(s1,s2) A (sl+s2) 2 (6-3b)

G (s1,s,s ) 1 1 (6-3c)
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:

K (s) = o (6-4)

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 non-normalized neutron population as

input. With the transforms of the operators G and K, given in

Equations (6-3a) and (6-4) respectively, the transform of the linear

operator becomes:

s +2s+l
Rl(s) = 2 (6-5)
( s3+2s2+s+no)
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 (5-15) 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 (5-25) reduces to

the same criterion as Equation (5-21) 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, Root-locus, etc.) which can be used to analyze

the transform of Equation (6-5), 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

numerical-ly, 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 step-function 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

~x1O I I r

-j :x
8n In I2 I n
o A A 1.9
:1I '
-2 :/ \ I

-4x102 2 21

4 8 12 16 20 24 28 32 36

Figure 15. "True" Response of System with No Delayed Neutrons Exponential Input


1.2x10- / I \

I iI I I-
I I \


I i I I I
I -2

-4x0o2 -=

p, = 0.05

4 8 12 16 20 24 28 32 36

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 step-function input of pe = 0.05.

Since the feedback is linear, Equation (5-19) can be used for the

calculation with Rl(s) given by Equation (6-5), Gl(s) I/A s, and

Kl(s) given by Equation (6-4).

Figure 17 shows the result of taking the inverse Laplace

transform of the linear terms, the second-order terms, and the whole

sum in Equation (5-19). 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 second-order 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 second-order 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 t-18.5 seconds). Some

comparative data for the two time response curves are tabulated in

Appendix G.


4 8 12 16 20 24 28 32 36


Figure 17.

Operator Method Response of System with No Delayed Neutrons
Step Function Input

Two-Temperature Feedback Model

The second reactor system considered is the two-temperature

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 right-hand side of Equations (2-19)

and (2-20) respectively. Here too, the feedback is linear so that

Equat-ion (5-19) 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+ (6-6)

reveals that the dk, as defined in Equation (5-20d), satisfy the


dk 0 for k = 1,2 (6-7)

so that the stability criterion of Equation (5-25) will be true when-

ever Equation (5-21) 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 second-order terms, and finally the whole sum in Equation (5-19)

to the time domain when the input to the system is a step-function of







Figure 18.

Operator Method Response of Two-Temperature
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 second-order 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 (2-11) with a step-function 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

step-function rather than a 5-function. 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.


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.




4 8 12 16 20

Figure 19.

Response of Two-Temperature 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 second-order terms are found to be adequate to

give good agreement with the true response of the systems, with the

second-order 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 third-order term will also contribute a non-trivial

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

terms but it is definitely not adequate with ohly a linear response

present unless the input is greatly reduced. For ke 0.01 the

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




In the work that has preceded, two methods of analysis of

the time-dependent 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 time-response

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 over-restrictive. Consider the classic example:

1 1 + x + x2 + ... (7-1)

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 first-and second-order terms. The small

number of terms which need be calculated helps point up the superi-

ority of this method over the direct, time-domain, 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 pseudo-random 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


the analytic models used in reactor analysis. This is one of the

more important uses of the work developed here.



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


combination with block diagram manipulation should yield some prelim-

inary stability information. Similarly, the method might hold promise

for coupled-core systems.

There are many other applications of the series approach both

in the reactor field and in a large number of other fields.



Region of Applicability of Zero-Power Kinetics

Start with Equation (2-14):


8po(t)-5p1(t) g(t T)ke( ) g(T-)ke(o)dadT
<< 1 (A-1)
-- p (t) M----------- < -
Po(t) f g(t-T)ke(T)dT


but, for a physically realizable g(t):

t T t
Sg(t-T)k (T) g(T-o)ke(a)dadT [f g(t-T)k (r)dv]2. (A-2)
0 0 0


86p(t)-6p(t) g(t-T)k()d2 (A-3)
"Po) t, << 1 (A-3)
i (g(t-T)ke(T)d


Spo(t)-6pl(t) t
5Po (t) f g(t-T)ke(T)dTv 1 (A-4)

Now, applying Equation (3-2)

t t
If g(t-T)k (T)dTj ~ f g(teT)I Ike(T)I d (A-5)

University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs