Hybrid computation of space and energy dependent nuclear reactor kinetics

Material Information

Hybrid computation of space and energy dependent nuclear reactor kinetics
Space and energy dependent nuclear reactor kinetics
Schultz, Kenneth Robert, 1942-
Publication Date:
Physical Description:
xiii, 171 leaves. : illus. ; 28 cm.


Subjects / Keywords:
Analog computers ( jstor )
Coolants ( jstor )
Eigenvalues ( jstor )
Error rates ( jstor )
Hybrid computers ( jstor )
Iterative solutions ( jstor )
Kinetics ( jstor )
Mathematical variables ( jstor )
Neutrons ( jstor )
Time dependence ( jstor )
Dissertations, Academic -- Nuclear Engineering Sciences -- UF ( lcsh )
Dynamics ( lcsh )
Hybrid computers ( lcsh )
Nuclear Engineering Sciences thesis Ph. D ( lcsh )
Nuclear reactors ( lcsh )
bibliography ( marcgt )
non-fiction ( marcgt )


Thesis--University of Florida.
Bibliography: leaves 166-170.
General Note:
Manuscript copy.
General Note:

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
Resource Identifier:
18242399 ( OCLC )


This item has the following downloads:

Full Text



Kenneth R. Schultz

A Dissertation Presented to the Graduate Council of the
University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy



To my Mother and Father


The author would like to express his appreciation

to the members of his supervisory committee. Dr. T. E. Bul-

lock, Dr. M. J. Ohanian, Dr. E. E. Carroll, Dr. R. Gaither

and Dr. Luehr for their assistance during the course of

this research. Dr. M. J. Ohanian and Dr. Thomas Bullock

are due special thanks for the many hours they have devoted

to directing this study. The entire faculty and staff of

the University of Florida Nuclear Engineering Sciences

Department were always helpful beyond the call of duty.

Dr. A. J. Mockel was of special help in the analysis of the

mathematical techniques.

Use of the hybrid computer was provided by the De-

partment of Nuclear Engineering Sciences at the University

of Florida.

The author's studies at the University of Florida

were supported by United States Atomic Energy Commission

Traineeships and Fellowships and a Graduate Assistantship

under National Science Foundation Grant No. GK-2106.

Thanks for this support are gratefully given.

Above all the author would like to thank his

wonderful wife, Mary Lou, for her patience and support

during his six years of graduate study. Without her en-

couragement and help this research could not have been done.



ACKNOWLEDGMENTS ..... ................................

LIST OF TABLES ....... ..............................

LIST OF FIGURES ................................ ...

NOTATION ..........................................

ABSTRACT ...........................................


I. INTRODUCTION ...............................

Space and Energy Dependent Reactor
Kinetics ..............................
Survey of Kinetics Methods ..............
The Use of Hybrid Computers .............

II. DESCRIPTION OF METHOD ......................

Mathematical Model ......................
Description of Hybrid Technique .........
Selection of Iterative Scheme ..........
Description of Hybrid System ............
Implementation ..........................

EXAMPLES ................................

Example 1: A One-group Model ...........
Example 2: A Two-group Model ...........
Example 3: A Two-group Model with
Temperature Feedback ..................

TECHNIQUE ..... .......................

Playback .................... ..........
Accuracy ................................
Run Time ................................


















DISCUSSION .................................

Recommendations for Further Study .......
Conclusions .............................



B. HYBRID KINETICS PROGRAMS ..................

EQUATIONS ..............................

REFERENCES ........................................

BIOGRAPHICAL SKETCH ...............................









Table Page

1. One-group Turkey Point Example ............. 60

2. Two-group Yasinsky and Henry
Example ................................. 64

3. Two Groups with Feedback J. B.
Yasinsky Example ........................ 73

4. Example 3: Reactor Data ................... 146

5. Example 3: Initial Conditions ............. 148

6. Example 3: Pot List ....................... 149


Figure Page

1. Hybrid computer ............................. 12

2. Hybrid diagram ............................. 24

3. Iterative procedure ....................... 25

4. Hybrid flow diagram ......................... 28

5. Example of nonconvergence ................... 43

6. Over-relaxation results .................... 46

7. Analog implementation ...................... 53

8. One-group example-Flux vs. distance ........ 57

9. One-group example-Flux vs. time ............ 58

10. Two-group example-Flux vs. distance ........ 61

11. Two-group example-Flux vs. time ............ 62

12. Example three-Thermal shape function ....... 69

13. Example three-Reactor power vs. time ....... 70

14. Example three-Fuel temperature vs. time .... 71

15. Data playback ............................... 76

16. Sensitivity to offset error ................. 78

17. Computation time vs. number of space
points .................................... 83

18. Computation time vs. integration
interval ................................. 86

19. HYKIN 3 analog neutronics .................. 126



Figure Page

20. HYKIN 3 analog thermal-hydraulics .......... 127

21. Eigenvalue error ........................... 158



Standard reactor physics notation as defined below

is used throughout this dissertation.




Z (E '-E



= Neutron flux

= Delayed neutron precursor density in the
Zth group

= External neutron source

= Total neutron cross section
Neutron cross section for absorption
= Neutron cross section for absorption
= Neutron cross section for fission

) = Neutron cross section for scattering from
group E' to E

D = Neutron diffusion coefficient

v = Neutron fission yield

;) = Energy distribution of neutrons produced
by fission

8 = Total delayed neutron fraction

S= Delayed neutron fraction to th precursor

= Delayed neutron precursor decay constant

Thermal Hydraulics

V = Reactor coolant volume

Volumetric specific heat of coolant

T = Coolant temperature

AH = Heat transfer area per unit coolant volume

U = Conductivity of cladding

h = Heat transfer coefficient

W = Coolant -flow rate

T = Fuel temperature

C = Specific heat of coolant
T. = Coolant inlet temperature
r = Fraction of fission power appearing in

q = Fission power

V = Reactor fuel volume

C = Conversion factor; fissions to energy

Pf = Fuel density

Cf = Specific heat of fuel

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



Kenneth R. Schultz

March, 1971

Chairman: Dr. M. J. Ohanian
Co-Chairman: Dr. Thomas E. Bullock
Major Department: Nuclear Engineering Sciences

The neutron flux in a nuclear reactor may generally

be described as a function of time, space, and energy.

Substantial errors can result if spatial effects are ignored

when studying the kinetics of large reactors or if spectral

effects are ignored in the study of fast reactor dynamics.

Hence, there is substantial current interest in the study of

space and energy dependent nuclear reactor kinetics.

The digital computer computations used in these

studies can be lengthy and hence expensive. An analog com-

puter, which solves differential equations efficiently,

might be used in kinetics studies but solution of realistic

time, space, and energy dependent reactor problems is beyond

the ability of an analog machine alone. By combining a

digital computer with an analog computer, a hybrid machine

is obtained that can possess the best features of both

digital and analog machines. A method of direct numerical

solution of the time, space, and energy dependent neutron

diffusion equation has been developed that utilizes the

unique capabilities of the hybrid computer.

The hybrid method uses an iterative discrete space-

continuous time algorithm. The space derivatives of the

conventional multigroup reactor kinetics equations are

represented by finite difference approximations. The re-

sulting very large set of coupled ordinary differential

equations in time can be solved quite readily on the hybrid

computer by multiplexing the analog hardware and using an

iterative procedure. The flux, precursor, and temperature

equations at one space point are set up on the analog com-

puter. The digital computer supplies the space coupling

terms from adjacent space points; it also observes and stores

the solutions for later playback. Feedback effects are

calculated by the digital machine and are implemented by

means of variable coefficients. Automatic digital rescaling

of the analog equations allows consideration of realistic

reactor dynamics studies involving flux level changes of

many orders of magnitude.

The hybrid method is well suited for solution of

the space and energy dependent reactor kinetics equations

and shows strong potential for use in nuclear reactor


engineering design. It offers reasonable accuracy at

modest cost. Hybrid solutions to typical spatial kinetics

examples and a description of the HYKIN 3 hybrid kinetics

code, which is functionally equivalent to the WIGL 2 digital

kinetics program, are presented.




During the 28.years since Enrico Fermi's group of

physicists brought the first nuclear reactor to critical-

ity on December 2, 1942, nuclear engineering has come a

long way. Reactor power levels have been raised from the

original 500 milliwatts of Fermi's Chicago Pile 1 to

thousands of megawatts. Temperatures in the more exotic

reactors approach 40000 F. This year over 60 per cent of

all the orders placed for electric generating stations

will be for nuclear power plants. In 28 years nuclear

reactors have evolved from primitive physics experiments

into sophisticated pieces of production machinery. Nuclear

engineering has become an established professional field.

Advances in the field of nuclear engineering have brought

about the need for more detailed calculational models.

Space and Energy Dependent Reactor Kinetics

The neutron flux in a reactor may generally be

adequately described as a function of time, space and

energy. The first reactor kinetics studies, however,

made use of the point model equations which ignore the

- 1 -

- 2 -

space and energy dependence of the flux. In the design

of large reactors with high power densities, space de-

pendent effects that had not been of previous concern

began to be noticed. It has been shown[]* that the

reactivitiy worth of a local perturbation can be much

greater than would be calculated by a space independent

analysis. It is necessary, therefore, to consider space

dependent effects in the analysis of nuclear reactor


The need to study the energy variable in detail

is peculiar to fast reactors. In these reactors the

system multiplication is a strong function of the neutron

energy spectrum. The predominant feedback effects, the

coolant density effect and the Doppler effect, are also

very spectrum dependent. It is, therefore, important to

have knowledge of the neutron energy spectrum during a

transient if the dynamic behavior of a fast reactor is to

be calculated accurately. Thus it can be seen that a

definite need exists for a method of nuclear reactor

analysis that can consider time, space and energy depend-

ence and that includes realistic feedback representation.

Superior numbers in square brackets refer to the
numbered references in the Bibliography.

- 3 -

Survey of Reactor Kinetics Methods

Before embarking on any program to develop a new

method of reactor kinetics analysis, we must first examine

what has already been done.

The basic equations that describe the behavior of

a population of neutrons is the Boltzmann neutron transport

equation. The full seven dimensions of the Boltzmann phase

space are too much to cope with for practical kinetic cal-

culations. In a large reactor the flux is essentially

isotropic. By ignoring the angular dependence of the

neutron flux, the time, space and energy dependent diffusion

model is obtained. 3]

Neutron Flux Equation

STET (r,E,t) = V D(r,E,t)VQ(r,E,t) Zt(r,E,t)

x @(r,E,t) + J dE' Z (r,E' + Et) (r,E',t)

+ x(E)X C (r,t) + x(E) dE' (1 (E')

x vCf(r,E',t)D(rE',t) + Sext(r,E,t)

Delayed Neutron Precursor Equation

8C (r,t)
t dE'vZ (r,E',t)O(r,E',t) X C (r,t)

- 4 -

These are the basic equations from which most

reactor kinetics analyses begin. The approximation resul-

ting from completely ignoring the angular dependence of

the neutron flux may seem crude, but it has been shown[2]

that the results of reactor kinetics studies using the dif-

fusion approximation compare quite favorably with the results

of much more rigorous treatments.

The treatment of this equation may be divided into

two general categories:

Point model studies.--In these studies the time

dependence of the neutron flux is assumed to be separable

from the other independent variables, space and energy.

4(f,E,t) = f(l,E) T(t)

This may be made rigorous by integrating over space

and energy rather than simply assuming separability.[48]

The rigorous integral parameters of the equations however

are then functions of time. This means that the "rigorous

point model" equations are still as much work to solve as

the time, space and energy dependent diffusion equation


Space and energy dependent studies.--These studies

attempt to account for the fact that the reactor kinetics

equation is nonseparable.

A great deal of work has been done using the point

model equations. This includes analytic studies and

- 5 -

numerical methods using digital, analog, and hybrid

computers. Extensive work on analytic stability analyses

has been done using these equations. Despite the attention

the point model equations have enjoyed due to their com-

parative simplicity, the lack of ability to consider

spatial and spectral effects means that point model

studies, by themselves, are not sufficient for reactor

kinetics analyses.

Those studies that do not assume that the kinetics

equations are separable but rather attempt to account for

the interaction of space and energy dependence as well,

may further be subdivided and classified into four categories:

1. Analytic Solutions

2. Numerical Treatments

3. Series Expansion

4. Factorizing Methods

Analytic solutions that include time, space and

energy dependence are possible for special cases. These

cases include the small variation, linearized problems that

can be handled by means of Laplace and Fourier transforms

of the time and space variables. These lead to space and

energy dependent transfer functions.[-12] The limitation

to small perturbation linear problems however restricts

the usefulness of these analytic transformation methods.

The reactor kinetics equations are analytically unsolvable

when realistic feedback effects are considered. Even with

- 6 -

no feedback and very simple, uniform system changes, the

analytic study includes expansion in infinite-series to

handle the space and energy dependence and involves a great

deal of computation.

Analytic methods of stability analysis have been

developed that consider the complete set of space and

energy dependent partial differential equations that de-

scribes the flux.[13,14] These are purely stability studies

and do not calculate the flux behavior.

Complete numerical treatments have been developed

such as the WIGLE and FREAK codes. 11 5-22] The only approxi-

mations inherent in these approaches are those due to dis-

cretization of variables and use of finite differences for

the derivative terms. The drawback of these methods is

cost. Since the required computations are both involved

and numerous, they are time consuming. On the large digital

computers necessary to handle such a problem, time and money

go hand in hand. Since computation time rises at least

linearly and possibly with the square of the number of

energy groups used, detailed spectral dependence is ex-

pensive to include and requires still larger computers.

There is no theoretical reason why the full seven-

dimensional Boltzmann equation, complete with delayed

neutron equations and all feedback effects, cannot be solved

numerically. From a practical point of view, however, this

is impossible anduannecessary. The greatest obstacle to

- 7 -

direct numerical computation is the small size required

for the time step. If the equations are linear, an implicit

time differencing scheme may be used. If the equation co-

efficients change with time, the time steps must be kept

quite short to insure accuracy. If an explicit time dif-

ferencing scheme is used, computational stability restricts

the maximum time step size to on the order of the prompt

neutron lifetime. In a fast reactor this is on the order

of 10-7 seconds. Because of this the large all-digital

codes using direct numerical integration can be slow and


Following the advice of Fermi who said, "When in

doubt, make a series expansion,"[21] much work has been

done along the line of expansion of the flux in a series

of time independent modes or nodes with time dependent


4(F,E,t) = ) mi(F,E) Pi(t)

This technique can be cumbersome, particularly in the case

of reactors with feedback which causes the flux equation

coefficients to vary with time and may cause the modes used

at the beginning of the problem to be completely unsuited

for the system at a later time.[1,24-3s] Error is always

introduced when the series is truncated after some finite

number of terms.

- 8 -

The factorizing approach[3l637] separates most of

the time dependence out of the kinetics equation:

<(r,E,t) = ((r,E,t) (t)

where 4 is a weak function of time. This method is further

subclassified by the approximations now made.

1. Improved quasistatic method uses large time
steps for 2t, small steps for -.

2. Quasistatic method ignores .
3. Adiabatic approximation ignores both
and 3.
4. Linearized adiabatic method assumes separability
and then uses the new value of ( to calculate a
new shape function 4 every so often, which in
turn is used to continue calculation of (.

At one extreme (4.) the factorization approach is the

simple point model method. At the other extreme (1.), it

is the complete numerical treatment. Accordingly, its

accuracy and cost may be adjusted to obtain a compromise.

The most successful approximate method of space

and energy dependent reactor kinetics analysis developed

to date appears to be the quasistatic approximation code

QX 1 developed by Ott and Meneley[36,37] at Argonne

National Laboratory.

Ott's results[37] show that the quasistatic

approximation yields results for fast reactors that are

almost identical to the "exact" results. This would cause

one to question the advisability of attempting to develop

- 9 -

a more detailed method of fast reactor kinetics analysis.

However, it must be noted that the "exact" standard

used for this comparison is the two-group WIGLE code. Just

as Yasinsky and HenryP[] found large discrepancies in

certain cases when comparing the results of point model

and space dependent kinetics analyses, it is quite possible

that detailed analysis of the energy variable will yield

new insight on fast reactor kinetics and control. Recent

work by Hitchcock, Hansen, and Henry[38] indicates that more

detailed spectral analysis may indeed be necessary in reactor

kinetics studies.

Use of an electronic analog computer, which treats

time as a continuous variable, eliminates the time step

size problems inherent in a digital computer approach. The

lack of storage and retrieval capability on an analog com-

puter, however, limits its use to one set of differential

equations that is solved simultaneously. On the largest

analog computers available a maximum of about 100 first-

order equations (50 second-order equations, etc.) can be

handled simultaneously. Due to this limit, a realistic

time, space and energy dependent reactor problem cannot be

solved on an analog computer.

Despite the fact that much work has been done during

recent years in the area of space and energy dependent

kinetics, it can be seen from the above discussions that

- 10 -

a great deal of interest remains in further development

in the field.

The Use of Hybrid Computers

Until quite recently all of the current methods of

space and energy dependent reactor kinetics analysis have

been implemented only on digital computers. A number of

different methods for solution of partial differential equa-

tions using hybrid computers have been developed.[3as The

hybrid computer is naturally suited for many reactor kinetics

problems. Several of the space and energy dependent kinetics

methods currently in use could readily be implemented on

a hybrid system. Recent applications of hybrid computers

to the solution of nuclear engineering problems[[4~1 and

the new hybrid facilities at Westinghouse and Babcock and

Wilcox among others indicate a growing interest in hybrid

computation by the nuclear industry.

Use of hybrid computers in the analysis of spatial

kinetics problems has been reported in three recent papers.

Saphier and Yiftah of the Israel Atomic Energy Commission[41-.

have attempted a one-dimensional kinetics analysis using a

discrete space--continuous time method. Their method is a

simplified version of the technique developed by the present

author. Feedback effects are not included and the amplitude

scale of the problems considered is fixed and thus of quite

limited range. Unfortunately Saphier and Yiftah were unable

- 11 -

to make their hybrid technique work due to hardware problems.

Wozny, Bartells and Bailey at Purdue[42] used a continuous

space--discrete time procedure to study a one-dimensional

Xenon oscillation problem. Spatial integration was done

on the analog portion of the hybrid computer using

Vichnevetsky's decomposition method. [43 Carter et al.

at Argonne National Laboratory[4. used a hybrid computer

to tackle an ambitious two dimensional kinetics problem

including temperature feedback effects. They used a fairly

crude nodal representation of the core equations. These

papers further indicate growing interest in nuclear appli-

cation of hybrid computers.

A hybrid computer consists of an analog computer

and a digital computer connected together with interface

devices, as shown in Figure 1, in such a manner that the

two may operate as one unit. Analog computers are particu-

larly well suited for the solution of large sets of simul-

taneous ordinary differential equations. The equations may

be intricately cross-coupled and highly nonlinear without

affecting the ease of their solution or the high speed at

which the solution may be obtained. By combining the

analog computer's ability to solve simultaneous differential

equations quickly with a digital computer's prodigious

memory and ability to make complex logic decisions rapidly,

we can obtain a hybrid system that combines the best features

of both analog and digital machines.

- 12 -










Figure 1. Hybrid computer.


'OU_ 1~



m ,~YY nc-- -~

~--YI-YI---- -

- ~~-CI~DI~~-D



- 13 -

Hybrid computers have been used extensively in the

past, primarily by the aerospace industry, to-help facilitate

what are basically large analog simulations. In recent

years interest has shifted from hybrid simulation to hybrid

computation. Both the analog and digital machines partici-

pate fully and equally in the calculations. One of the

fields of major interest for hybrid computation is the

solution of partial differential equations.["3 9,. The

methods used may generally be divided into two categories:

1. Continuous Space-Discrete Time (CSDT) Methods

2. Discrete Space--Continuous Time (DSCT) Methods

The names are self-descriptive. In CSDT methods the partial

differential equations are reduced to sets of coupled ordinary

differential equations in space by replacing the time de-

rivatives with finite difference approximations. The space

equations are solved on the analog machine at each point in

time. This method has been used to study Xenon oscilla-

tions in nuclear reactors. [42 The Xenon problem has a very

slow time scale. For normal kinetics studies the need to

use a very small time step makes the CSDT method impractical.

The method is further limited to one space dimension.

In the DSCT method the space derivatives are replaced

by finite difference approximations. This yields a large

set of coupled ordinary differential equations in time.

- 14 -

These are solved on the analog portion of the computer,

treating time as a continuous variable.

The technique developed at the University of Florida

is a discrete space--continuous time method of direct nu-

merical solution to the time, space, and energy dependent

neutron diffusion equation that utilizes the unique capa-

bilities of the hybrid computer. The method appears to be

well suited to the solution of the kinetics equations but

certain limitations are imposed by the characteristics

of nuclear problems. Specifically, the means of data

playback, the computation time, and the method's accuracy

are affected by the type of problems to be solved. Auto-

matic rescaling is imperative due to the large range of the

flux encountered in typical reactor problems.



Mathematical Model

The equations of the diffusion theory model as

presented in Chapter I describe the space, time and energy

dependent behavior of the neutron flux in a reactor and

comprise the basic nuclear model for most reactor kinetics


1) Flux

1 4 V DVQ E + Z (E' E) (E') dE'
v t t oJ s

+ X(E) I v(E')(1 ((E')) Zf(E') (E') dE'

+ I Xt(E) XkC (1)

2) Delayed Neutron Precursor

ac = X (v(E') v(E') Ef(E') D(E')) dE' C


q,D,jfEt and s = f(r,E,t)

Ct = C (r,t) and v,v and X = f(E)

- 15 -

- 16 -

The system feedback effects are contained explicitly

in the time dependent coefficients of the above equations.

The feedback terms that must be included are density changes

in the fuel and coolant due to.temperature and pressure

changes, thermal distortion of the core, and microscopic

cross-section changes due to the Doppler effect. To cal-

culate the effects of these feedback phenomena, the

temperature must be calculated along with the flux level

throughout the reactor. An approximation to the actual

temperature distribution is the two-component represen-


Fuel Temperature

dTf Vc 1 1 -1
pfCf dT ( r)q V AHh+ (T T- ) (3)
ff dt Vf AHU AHho f c

Coolant Temperature

S c t1 1 f Tc + V i
c + (T TT

+ rq V. (4)


q = C I dE Ef(r,E,t) (r,E,t)

- 17 -

From the temperature distributions it is possible

to calculate the time dependent coefficients of the flux

and precursor equations. Each of the parameters may be

represented as

aT c

This completes the feedback loop. Because of this feedback

the reactor kinetics equations are nonlinear. Standard

reactor notation is used throughout.

A standard method of treating the energy dependence

of equation (1) is the multigroup approach. The energy

range (0 to 10 Mev) is divided into a number of segments

or groups. The flux equation, (1), is then integrated

over each segment, one by one, obtaining for each group

one equation. An approximation obtained from a static

calculation is used to represent the energy dependence of

the flux in the energy integrals and suitable energy

averages are defined for the group constants:

E E (r,Et)dE'
A k .k E a,
(k r,t) (r,E',t) dE' ak (r,t) E
E Ek
k-l k
E 'i(r,E',t) dE'


- 18 -

The resulting time dependent multigroup diffusion

equations are the beginning of most space and energy

dependent reactor kinetics analyses.

kth Prompt Neutron Energy Group: 1 < k < G

a4, G
1 kk + G
V- Xsk ,k
Vk -D Vk k t k k =1 '-k k

+ Xk(l 8) v fk'
k '=1 k'

+ Xk XCQ + Sk (5)

th Delayed Neutron Precursor Group: 1 < Z < L

3C G
= 8 ~ k~ X C (6)
k=l k


k' Dk3 tk s ,k'IkSk and C, = f(r,t)

The thermal-hydraulic equations remain unchanged

except for replacing the energy integral in the fission

energy term with a sum over all the energy groups:

q(r,t) = C X Ef (r,t) Dk(r,t)
k=l k

- 19 -

To handle the space dependence by direct numerical

techniques, the space derivatives are replaced by finite

difference notation using the central second difference.

In cartesian coordinates, for example,

d2z j+1 20 + j-1
d x. (Ax.) 2

The central difference is used because it is the most ac-

curate of the three-point finite difference approximations.

Using this approximation the multigroup time dependent

finite difference model of the reactor is obtained. (For

simplicity of notation slab geometry and one space dimen-

sion are used.)

kth energy group:

d, Oj+1 24 + ~j-1 G
1 k D k k k X
k dt k (Axj 2 k'=l Sk'k k?

G j j
+ k tk k + Xkl k=1 fk k'

+ Xk C (7)

kth precursor group:

dC G
k= Vfz k X (8)
k=l k

- 20 -

where the subscripts and superscripts denote:

j = space point number 1 < j < N

k = energy group number 1 < k < G

A = delayed neutron precursor 1 < k < L
group number

Since the thermal-hydraulic equations contain no

space derivatives, they remain unchanged. A set of tem-

perature equations may be defined for each space point:

Thermal-Hydraulic Equations:

j dT v
PfC d (1 r)qj + AH1 -3
f f dt V AHU AHho f c

9p HC dT (i
at Ao T-) + 2 V

fj H H
x T. Tj + rq [] (10)
in c V


q3 = C E (t) k(t) (11)
k=l k k

And feedback

a (t) = G(0) + (T (t)) (T(0)
k k D f 2

p aj
+ "kapJ(t) Pi (0)J (12)
ap c c c I

- 21 -

The standard procedure for solving this large set

of coupled ordinary differential equations on a digital

computer is to replace the time derivatives with finite

difference approximations and use a standard numerical

integration routine. If the equations are linear, an

implicit time differencing scheme may be used. The calcu-

lation involves either spatial iterations or the inversion

of a large but fairly sparse matrix at each time step. If

the equation coefficients change with time, due to feedback

effects or an external perturbation, the time steps must be

kept quite short to insure accuracy. If an explicit time

differencing scheme is used, computational stability requires

that the time step size must be smaller than the inverse of

the largest eigenvalue of the coefficient matrix, which is

related to and on the order of the prompt neutron lifetime.

In a fast reactor this is on the order of 10-7 seconds.

Because of this, the large all-digital codes using direct

numerical integration can be time consuming and costly.

By modeling equations (7) through (11) on the analog

portion of a hybrid computer, time is treated as a contin-

uous variable. Thus time step size limitations are avoided.

This set of equations is too large to be modeled simulta-

neously on even large analog computers. By setting up a

block of only a few space points at a time on the analog

portion of a hybrid computer and using digital control and

- 22 -

memory as the analog integrates each space block

sequentially, the equations can be solved quite readily.

In this study the space block contained one space point.

The algorithm is an iterative discrete space-continuous

time method with variable coefficients.[46]

The iterative discrete space-continuous time method

has been investigated by Howe and Hsul 46 for nonnuclear

application. Their study must be considered a preliminary

investigation because the hybrid procedure was simulated

using a digital computer. They considered very few space

points and simple "model" problems with constant coeffi-

cients. In so doing this eliminated many of the physical

realities that seriously affect the hybrid method. In the

present study realistic reactor problems were tackled using

actual hybrid hardware.

Description of Hybrid Technique

The solution proceeds in the following manner: the

flux, precursor and temperature equations at one space

block are solved on the analog computer for a short time

interval. The coupling terms from neighboring space points

are supplied by the digital computer which also samples and

stores the solution for later playback during the calcu-

lation of the adjacent space blocks. With the analog

computer operating under digital control the equations are

- 23 -

solved at each space block, one by one, using coupling

terms stored in, and then retrieved from the digital memory.

Figure 2 shows a hybrid schematic diagram.

Use of the second central difference and sequential

solutions means that this procedure must involve iteration.

Equation (5) which describes the time behavior of the neu-

tron flux at space point j requires as input the solutions

of 4(t) at the adjacent points j 1 and j + 1. Since the

points are handled sequentially, a solution for j-l (t) is

already available for playback from digital storage. Dj+ (t)

however has not yet been obtained. Hence one must use an

estimate of what Dj+1(t) will be. Since it is only an

estimate, the solution obtained for J (t) will be in error.

When this procedure has been done for all space points, a

first iteration .estimate of all J (t) will have been ob-

tained. The entire procedure must then be repeated using

the previous solutions as the estimate for j+l1(t), yielding

a second iteration estimate of J (t) at all space points j.

The procedure is repeated again and again until subsequent

iterations agree within some convergence criteria. Figure

3 illustrates the iterative procedure. After the solution

has converged at all space points for the first integration

interval, the final values of flux, precursor and temperature

are used as initial conditions for another time interval.

This procedure may be repeated to follow the solution as

-- 24 -




Figure 2. Hybrid diagram.





.-- i--r

- 25 -

(p)th iteration of ( +(t) from
digital storage
(p)th iteration of J (t)
(p+l)th iteration of i (t) being
generated on the analog
(p+l)th iteration of Dj-1(t) from
digital storage

xj-I xj Xj

Figure 3. Iterative procedure.

- 26 -

long as is desired. An initial guess must be used to start

the iteration procedure. For the first time interval the

flux is assumed to remain constant. For subsequent time

intervals linear extrapolation of the previous solution is

used as the initial guess.

External perturbations, time variant parameters and

feedback effects are calculated by the digital computer and

are implemented on the analog machine by means of variable

coefficients in the differential equations.

On an analog computer accuracy is fixed at a

percentage of the machine's full scale range by noise.

Hence a small signal is less accurate than a large one.

Because of this, the useful dynamic range of the dependent

variables is limited to about a factor of 20. If the value

of a variable falls below 1/20 of full scale, accuracy *

becomes very poor. The restriction to a factor of 20

change in dependent variables is a severe limitation for

nuclear reactor kinetics studies. In these problems it is

often of interest to follow the reactor flux levels through

changes of many orders of magnitude. Therefore an auto-

matic rescaling feature was incorporated in the hybrid

method which frees it from the usual restriction of the

analog computer dynamic range. If at any time the solu-

tion at a space point exceeds the machine range of 100

volts, that point may be automatically rescaled and the

- 27 -

solution continued. At each point the scale factor is

adjusted so that the maximum value of the solution during

the time interval being considered is slightly under 100

volts. This allows the best use of the computer's avail-

able range. Although the automatic rescale feature is

quite straightforward, it appears to be a new idea and is

particularly applicable to nuclear reactor problems. A

flow diagram of the hybrid algorithm is presented in Figure


Although the scale factors are automatically

adjusted at each space point as the iterative scheme con-

verges upon the solution, the scale factors are constant

during each integration interval. Thus it is important

that some foresight be used in the selection of the length

of the integration interval. This length must be such that

the changes in the dependent variables do not exceed a

factor of 20. This in no way restricts the total range

that can be handled. Many integration intervals may be

stacked end to end, each with up to a factor of 20 change.

The automatic rescaling is handled by simply

auctioning the scale factor to be used to the highest

bidder among the variables to be used at the space point

being computed. If any variable exceeds 100 volts, its

scale factor is doubled, the auction repeated, and the

calculation performed again. This scaling is relative

- 28 -

Figure 4. Hybrid flow diagram.

- 29 -

scaling. That is, the ratio of any two variables at a

space point is not changed although their magnitudes are.

Space and time dependent absolute scaling may also be

added to the hybrid kinetics method quite easily. In this

way the magnitude of all variables at all space points and

at all times may be kept within the desired 5 to 100 volt

range. Absolute scaling has not been made automatic; the

scale factors at various space points and times must be

selected and fixed before the solution is performed. With

a hybrid computer that had digital pots and digitally con-

trolled remote patching for input and feedback impedances

the absolute scaling could also be performed automatically.

It is questionable if this would be desirable. The level

of complexity of the hybrid program would be significantly

raised and the same effect could be obtained by simply *

using more forethought with the nonautomatic scaling before

the problem was run.

Analog computer scaling can be a particularly nasty

problem for reactor kinetics studies because of the dif-

ferences in behavior of variables. The fluxes and pre-

cursors can vary by several orders of magnitude at different

space points while the temperature distributions are

generally fairly flat. The flux values at any space point

may change by many orders of magnitude during a reactor

transient while the precursors and temperatures will not

- 30 -

be expected to change by more than a factor of 10. Because

of these differences the scaling that is appropriate at one

time and a given space point may be completely wrong at a

different time and/or location. By using the space and time

dependent absolute scaling and automatic relative rescaling

described above the hybrid kinetics method can readily

handle these problems.

Even the most clever techniques cannot solve all

scaling problems. It is always possible that the equation

coefficients dictated by the scaling are not obtainable on

the analog machine. Gains in the range from 10-4 to 102

for summers and 10-2 to 104 for integrators are generally

available with reasonable accuracy on analog computers.

If scaling of the equations yields coefficients outside

of this range, the equations cannot be handled by the analog

machine. It should be noted however that the factor of 106

range in possible coefficient values is quite large. If a

problem cannot be scaled to within these limits, the ma-

chinery is probably trying to indicate what close inspection

of the physical model should have revealed. That is, that

the model is unnecessarily exact, and that a suitable ap-

proximation may be made that will simplify the problem

considerable at virtually no cost in accuracy.

Some of the approximations that are often made in

digital kinetic studies are particularly useful in the

- 31 -

hybrid method. If 1/vk in equation (5) is assumed to be

zero for the higher energy groups, the flux equations for

these groups are reduced to algebraic relationships. This

approximation is useful and valid for studying slow tran-

sients. If the precursor concentrations in equation (5)

are assumed constant, equation (6) may be eliminated. For

rapid transients this approximation is valid. Use of one

or the other of these approximations, whichever is appli-

cable, allows the analog integrator gains to be kept within

the available 106 range quite easily.

Selection of Iteration Scheme

Iteration is a standard technique for the solution

of second-order boundary value problems.['17. Iteration is

necessary because of the coupling terms from adjacent space

points that arise from use of the central difference repre-

sentation of the second derivative in space. Use of the

backwards second difference to represent the space deriva-

tives of the diffusion equation would eliminate the need for

iteration. This scheme, however, is computationally unsta-

ble[47.] and will not work. The unsuitability of this method

has been experimentally verified.

A number of different iterative schemes exist. Three

schemes were investigated analytically and experimentally

for use in the hybrid solution of reactor kinetics problems.

- 32 -

These schemes were: Simultaneous Displacements, Successive

Displacements and Odd/Even Iteration. Methods of iter-

ation acceleration were investigated in an attempt to

improve upon these schemes. In order to analyze these

techniques the model of the reactor was simplified. No

feedback effects were included, and one energy group was

used. The resulting finite difference flux equation is:

1 dO D (j j+l
v a- X- '(

- 24) + ij-1) -_ + (1-_)vE f3

+ XC3 + S3

dCJ- B Xc
dt f


C3 (t) =





2N -
(t) = 2N eYit
i(t) = e i

C (t) = c eYt

Retaining only the one dominant terms of each series

S(t) =3 eYit

CJ(t) c3 eYkt

Note that


- 33

J(t) y C3

k k
c (t) ykc

P = Yi

eYit = Yij (t)

eykt = ykC (t)


A = yk



so that

ACj = 3zEf j ACj

c = VZf j

thus equation (13) becomes

1 r< = D3 (j+l
v Ax21

- 20 + ~j-l) + ((l-8)vZ Ea

+ vy X' + j sj
f X-- + A ext



J (t) A (t)

cJ (t) ACj (t)

- 34 -

D (j+l 1- r 2D
1 j + A+I v a v Ax2

+ Sj =0 (15)
D ext

This may be represented as

j+1 + aj3 + j- +j = 0 (16)


aj -K2 + D V +a (1-(1 )] (17)
D3 a V (+A f


_j A Ax2 Sj
D ext

The aj and B3 are assumed constant in time. The

particular iteration scheme used will define a relation-

ship using equation (16) whereby the (p + 1)th iterate for

P may be obtained from the information from the pth iter-

ation. Equation (16) may be represented as

AT + ly = 0 (18)

where A is the N x N matrix

- 35 -

ca' 1 0 0 1 0

1 a2 1 0

0 1 a3 1

A = 0 0 1 a4 1

1 0
N-1 1

0 01 a N

(t) p2

= and y =

N (t) BN

N is the number of space points used.

To solve equation (18) without directly inverting

the matrix A, iterative techniques may be defined that

allow a (p + l)th estimate of the solution T to be ob-

tained given a (p)th estimate. The iteration equation is

'(p + 1) = BO(p) + Cy (19)

where B is the iteration matrix.

In order for the iteration scheme to be successful,

it must be consistent and convergent. Consistency means

- 36 -

that iteration using the true solution to the original

differential equation must return the true solution:

Dt = BDt + Cy (20)

Convergence means that an error in the initial trial

function will die away as the iterations progress. The

rate of convergence of an iteration scheme is a measure

of its efficiency and hence its relative merit.

Let the true solution be defined as the pth

iterate plus the pth error:

t = (p) + E(p) (21)

Substitute equation (21) into equation (20) and subtract

equation (19)

(p + 1) = BE(p) (22)

Similarly with equation (18),

AE = 0

Equation (22) shows that the error obeys the homogeneous

version of the iteration equation and

E(p) = BP~(0)

The iteration scheme is convergent if the magnitudes of

all the eigenvalues of the matrix B are less than one.

- 37 -

Assuming that B has a complete set of eigenvectors and

eigenvalues given by Be. = IX., one may expand E(0) in

terms of these eigenvectors:

E(0) = ai.
1- 1


s(l) = BE(O)

B a.e

= C ai.X.iei

E(P) = e a.2e.

If the eigenvalues are ordered such that

l11l > 1x21 > *** > IxI


E(p) = a1X1 e

e(p) = aleP n Alei (23)

p Jn A or e- where
That is, the error decays as e n or e where

v = -kn X1 is defined as the rate of convergence. Thus

the eigenvalues X. of the iteration matrix, B, determine

the convergence rate of the iteration scheme. If XA11 < 1,

- 38 -

then the scheme is convergent. In the general case the

result still holds where ei now includes pseudo-eigenvectors

of B.

The three iterative schemes were examined for

convergence and consistency and their rates of convergence

were compared. The details are presented in Appendix A.

The results are summarized below.

A. Simultaneous Displacements

The equation for the method of simultaneous

displacements is

-aJj (p + 1) = j+l(p) + 0j-l(p) + gj (24)

The solutions at the various space points are

calculated sequentially. Hence the (p + 1)th solution

for j-1 is available during the solution of O (p + 1).

This information, however, is not used. The complete

(p + l)th iteration is performed using only pth information.

This scheme is consistent and is convergent as long as

a01 > 2(1 -()) for all j (25)

where N is the number of space points and a. is defined

by equation (17). The rate of convergence is given by

v = -kn [(1 )- )) (26)

- 39 -

B. Successive Displacements

The iterative equation for this method is

-aj (p + 1) = j+1(p) + j-l(p + 1) + 0 (27)

Again the solutions at the various space points are

calculated sequentially. In contrast to the method of

simultaneous displacements the information of the

(p + l)th iteration is used as soon as it is available.

This scheme is also consistent and is convergent if

(a2 >4(1 T ) 2)2 for all j (28)

The rate of convergence is given by

v -n 4(1 ( ) (29)

C. Odd/Even Iteration

In this iterative method the solutions are not

obtained sequentially. Rather, the order of solution is

1,3,5,...N; then 2,4,6,...N-1. For the "odd" solutions

the iterative equation is

-aJj (p + 1) = j+p) + j-l(p) + j (30)

For the "even" solutions,

-aJ (p + 1) = Oj+1(p + 1) + 0j-l(p + 1) + J


- 40 -

is used. This method is equivalent to successive displace-

ments. It is consistent and is convergent if

(aj)2 > 4(1 () 2 (32)

The rate of convergence is given by

v) -n 4(1 (4 1 )2) (33)
a2 N

Since Odd/Even Iteration is equivalent to

Successive Displacements and slightly more difficult to

implement it will not be considered further.

Squaring the convergence criteria for Simultaneous

Displacements, equation (25), we obtain

a > 4(1 )2)2 (34)

This is the same as that of Successive Displacements,

equation (28). Furthermore the rate of convergence of

Successive Displacements is twice that of Simultaneous


v -[n 4(l (I) 2)
V [T 2
suc a
Vsim -n [(l (7

2 n( 2) + Zn( (T)

Zn(_) + n(l l(,) 2

- 41 -

1 2 1]
2 kn( 2) + I- 2)
2 1 2
( l) + T)2
-a 2


Experiment verifies that the method of Successive Dis-

placements does indeed converge twice as fast as Simul-

taneous Displacements. Of the three iterative methods

examined so far, Successive Displacements is the best.

As a consequence it was selected for use in the hybrid

method of reactor kinetics analysis.

From equation (15) it is evident that the

possibility exists that the solution will not converge.

If the nuclear parameters of the system being studied

are such that

2 1 T 2
(a) < 4(1 -() )2 at any point j (35)

then the largest eigenvalue of the iteration matrix will

be greater than unity and the iteration will diverge.

Substitute into equation (25) the definition of a ,

equation (17). If
x1 1 1
D a v A+X 2 N


then the iteration will not converge.

- 42 -

Several important conclusions may be drawn from

equation (36).

1. If the term

[Za + (1- -1 -v
a v A+- fl

is negative, it is possible that the iterations
will not converge.

2. As the space mesh size Ax is made smaller, the
left-hand side of equation (36) is decreased,
making the inequality more easily met. Thus,
the likelihood of nonconvergence is increased.

3. Convergence is improved by increasing Ea and
decreasing vZf, that is, by considering a less
reactive system.

These observations have been experimentally

verified. It was further discovered that a nonconvergent

calculation could be made to converge by time scaling the

equations to run more slowly on the analog computer. The

existence of a regime of nonconvergence depends upon the

specific parameters of the system being studied. If such

a regime exists, it may be avoided by time scaling the

problem down, by considering a less reactive perturbation,

or by using a coarser space mesh. Figure 5 shows the

results of a typical problem which was convergent if

calculated for 26 space points, but was nonconvergent if

32 or more space points were used.

- 43 -

.-- N = 32
.- y. .. N = 26


S 4-

(t) at reactor
center line

0 2 .3 .4 .5


Figure 5. Example of nonconvergence.


- 44 -

To verify that the simultaneous solution of more

than one flux energy group would not affect the iterative

procedure, the analysis was extended to include more

groups. Because of difficulties that later were found

to be due to a scaling error, the effect of multigroup

equations in a reflector region was of particular concern.

The analysis is presented in Appendix A. The results of

this analysis show that the addition of more energy groups

should not affect the iterative scheme adversely.

Iteration Acceleration

A number of schemes are available to use in

accelerating an iterative method.[47] Two schemes,

over-relaxation and asymptotic extrapolation, were ap-

plied to the method of successive displacements without

marked success. Successive over-relaxation is defined

by the equation

O (p + 1) = y (0+(p) + j-l(p + 1) + i3]

+ (1 y) J(p)

where y is the over-relaxation factor. It is theoretically

possible to increase the rate of convergence greatly by

proper choice of y.[47] The maximum rate of convergence

with optimum y is

- 45 -

1 4 2
1 / -v cos (-)
v = -an Ca N
1 + l- 4 -os2
1 -2 cos ( )
aj2 N

In practice the "optimum" over-relaxation factor is

generally determined experimentally. Unfortunately, the

expectation for great improvement in rate of convergence

was not borne out. An all-digital simulation of the

hybrid algorithm was written to check the entire hybrid

procedure initially in the absence of any possible hard-

ware related problems. This simulation was used to

determine the over-relaxation factor. Although this

factor will depend upon the characteristics of the

particular reactor being modeled, some typical reactor

models yielded "optimum" factors in the neighborhood of

1.70. The reduction in number of iterations necessary

was considerable as shown in Figure 6. With over-

relaxation factors greater than 1.70, the solution

oscillated. Factors greater than 1.75 caused nondamped

oscillations; the solution did not converge.

When over-relaxation was used on the true hybrid

model rather than the digital simulation, the solution

would not converge for factors larger than 1.24. As

before, the solution did not diverge until y exceeded

1.75, but for y between 1.24 and 1.75 small, bounded, but

- 46 -

.-.-- 16-bit digital simulation
...>.. 12-bit simulation or hybrid method

Results shown for a
typical reactor example

"a' I


1.50 1.75

Figure 6. Over-relaxation results.





0 1.0

- 47 -

nondamped oscillations about the true solution occurred.

By including some of the physical realities of the hybrid

algorithm in the all-digital simulation, it was determined

that these oscillations are caused by the 12-bit reso-

lution of the analog-to-digital converter. Why the 12-bit

accuracy of the ADC's should cause this oscillation is not

known. These results are summarized in Figure 6. It is

possible that a nonstationary technique that allows the

over-relaxation factor to be reduced as the solution

converges could improve the situation. This possibility

was not investigated.

Use of the 1.24 over-relaxation factor does result

in up to a 30 per cent reduction in number of iterations

necessary for convergence. The maximum benefit is ob-

tained when the initial guess for the solution is very

poor. The number of iterations that are necessary may be

reduced by using short integration intervals and by

extrapolating the solution from the previous interval

for use as the initial guess.

The nonstationary convergence acceleration method

known as asymptotic extrapolation was also tried.[46]

The residual of an iteration is defined as

r(p + 1) = Q(p + 1) #(p)

Using the definition of the iteration error,

- 48 -

E(P) = true p)

r(p + ) = + + 1) E(p)

Thus from equation (23)

r(p + 1) = Xr(p)


True = T(p) + r(p) + r(P + 1) + ... + r(m)

S(p) + (p) + Xr(p) + X2(p) + XrrCp)

S(p) + (p)


true (p) + F(p) r(p)
true F-(p)
r(p 1)

This equation defines the asymptotic extrapolation technique.

The transients introduced by the extrapolation took so long

to die out that any acceleration of convergence was lost.

Description of the Hybrid System

The hybrid technique was developed using the

Hybrid Computation facility at the University of Florida

Department of Nuclear Engineering Sciences. The digital

portion of this facility consists of an IBM-1800 process

- 49 -

control computer with 16K of core storage. The analog

portion is an Applied Dynamics-80 volt analog computer

with 80 amplifiers. The interface between the two

machines is provided by 16 direct digital lines to and

from the analog control patchboard, a 16 channel, 12-bit

50 KHz Analog to Digital Converter, an 8 channel, 15-bit

50 KHz Multiplying Digital to Analog Converter and a

17-bit BCD digital voltmeter with remote addressing. An

extensive library of general purpose hybrid subroutines

is available. These subroutines were written at the

University of Florida in assembly language and allow

FORTRAN control over the hybrid facility.

The relatively modest number of digital to analog

converter channels available on the facility restricts

the scope of the problem that may be attempted. Eight

channels will handle at most four energy groups. By

doubling up on the DACs, four precursor groups may be

treated at the same time. If temperature and other

equations requiring DAC channels are desired, the maximum

number of energy groups is reduced.

The most complicated problem run in the course of

this research was a one-dimensional model using two energy

groups, one delayed neutron precursor group, and one fuel

and one coolant temperature at each space point. Only one

space point at a time was modeled on the analog. Up to

- 50 -

50 space points were used. Fuel and coolant temperature

feedback are included in the model. This version of the

hybrid kinetics program, which is called HYKIN 3, was

modeled after and is equivalent to the WIGL 2 digital

spatial kinetics code..[20] HYKIN 3 is described in

detail in Appendix C.

On a larger hybrid system such as the new facility

at Westinghouse with 32 DAC channels up to 16 energy

groups with 6 delayed neutron groups and 10 temperature

equations could readily be solved. Alternately, larger

space blocks could be used and the solution time reduced

for a few group problem. Since the equations are solved

simultaneously on the analog for each space block the

addition of more equations does not increase the computation



The hybrid reactor kinetics program is implemented

in a very straightforward manner. The differential

equations for the flux, precursors, and temperature at

each space point as given by equations (7) through (11)

are programmed on the analog computer. The feedback

effects, equation (12) are calculated on the digital


- 51 -

In typical reactor kinetics problems with temperature

feedback, the flux behavior is very sensitive to changes in

temperature. Because of this care must be used to properly

scale the temperature equations to insure good accuracy.

To aid in this scaling the slowly varying equations, those

for the precursor and temperatures, were divided into two

parts. The variables were separated into a constant term

and a time dependent change term, e.g.,

Tf (t) = Tf(0) + AT (t)

Only the change term is computed on the analog computer.

The constant term is kept in the digital memory and is

played back to the analog solution of the change term as

required. This allows the change term to be computed

with greater accuracy than could be obtained if the

entire variable were modeled on the analog. The initial

conditions for the analog equations are supplied by the

digital machine via DAC channels. The coupling terms

from space points (j 1) and (j + 1) in the flux equations

are also supplied by the digital machine. The complete

solutions in time of 3j+l and ~j-1 are available in the

digital memory and are fed to the analog as a sequence of

points and slopes between these points using two DAC

channels per equation. First-order interpolation with

these points and slopes is performed on the analog computer

- 52 -

to reconstruct the coupling terms. Those coefficients of

the equations that are affected by temperature feedback,

and any others that vary in time, are supplied in a

similar manner. The digital computer observes the

reactor temperature during the pth iteration and uses

this information to calculate the feedback effects and

new coefficients for the (p + 1)th iteration. These are

output to the analog as points and slopes at the appro-

priate times to coincide with the flux solution. The

coefficients are reconstructed by analog first-order

interpolation from these points and slopes. The coef-

ficients are then multiplied by the appropriate variable

using analog multipliers.

Equation coefficients that are space dependent

but are not functions of time are implemented using

digitally controlled switches and manual pots. A much

better method would be to use digital pots but these are

not available on the University of Florida Nuclear Engi-

neering Hybrid facility. Servo-set pots would not be

adequate as they are generally quite slow, requiring

several seconds to set. Their use would increase the run

time and hence cost of the hybrid method by many orders

of magnitude. A simplified schematic illustrating this

implementation is shown in Figure 7.

- 53 -

= Di j+1 + j-1+
l Ax2

+ (1--P) VE 2

+ (2D + E
Ax2 a

- (1-f)v) }
f 1 ( 1

+ XcJ

i j
)1 = D (x,t)
1=) Use analog interpolation with points and
,3 (t) slopes from DACs
i = C (x,t)
a a

S= (x) Use switches or digital pots
constant Use manual pot
X = constant Use manual pot

D +1 + -1)






Analog implementation.


Figure 7.

- 54 -

Timing is very important in this hybrid method.

It is essential that data be played back at the same time

during an integration interval as it was taken. An

analog clock consisting of two adjustable pulsers wired

together to form a digital oscillator was used as the

master time base. The digital computer was slaved to

this clock and it in turn controlled the analog computer.

This timing mechanism proved adequate although a crystal

clock would be preferable if one were available. By

carefully studying the timing of the hybrid subroutines

and by making use of programmed delays and natural delays

due to switch closure and amplifier reset times, it was

possible to ensure that the data was played back at very

close to the same time during an integration internal as

it was sampled.

All the digital programming was done in FORTRAN

using general purpose assembly language hybrid subroutines.

This allows easy modification of the programs which was

imperative during development of the method.

Data output is fairly slow with the present DAC

control program. Each call to the routine outputs only

one value and requires 274 seconds. The DAC is capable

of outputting data at 50 KHz or 20 seconds each. This

capability ought to be utilized by making the program output

whole arrays at a single call. This would speed up the

- 55 -

hybrid program somewhat. Now that the time and other

adjustments are determined for the particular pair of

analog and digital computers used, it would be possible

to make a special purpose hybrid kinetics subroutine.

This could be written-in assembly language and would

handle all data input, output, and interpolater control

functions. This would offer no advantage beyond elegance,

however, since the analog interpolaters limit the speed

of the solution, and they are going as fast as they can

already. The essential parts of the digital portion of

the hybrid program are shown below:

Do for each space point

Output initial conditions and begin main analog

-Do for each time sample point

Do for each flux and coefficient using analog

Output on DACs the point and slope for

Wait for clock

Output control to: enable DACs, reset analog

Sample solution for all fluxes, precursors, and

Wait 440 seconds for switch closure and reset

Output control to: disable DACs, operate analog

-Reset analog



During the period from June, 1968, to the present,

work has taken place at the University of Florida on devel-

opment and investigation of the hybrid method of space and

energy dependent reactor kinetics analysis. A large number

of one and two energy group calculations have been made.

The problems generally considered are those of a reactor,

initially critical and at steady state, experiencing a

localized perturbation at time zero. One-dimension slab

geometry was used in the problems considered. The hybrid

method is not restricted to any particular geometry nor is

it limited to only one dimension. Howe and Hsu 6] have

shown that extension of the technique to more dimensions

is straightforward. From 5 to 50 space points were used.

Some typical examples are shown below.

Example 1: A One-group Model

Figures 8 and 9 show the results from a hybrid

computation of a one energy group model of a large thermal

reactor. The reactor is the Westinghouse PWR at Turkey

Point, Florida. The example shown used 24 space points

across the core, which is 388 cm. thick. The perturbation

- 56 -

- 57 -


7 S



2 *

T 0.0 SEC-7,

6 12 18

0 .2 .4 .6 .8 1.0


One-group example-Flux vs. distance.

Figure 8.

- 58 -


7 6

3 -

I -!------------------ l2

0 .1 .2 .3 .4

Figure 9. One-group example-Flux vs. time.

- 59 -

is a step change of -1.03 per cent in Ea in the left-most

one-third of the core. The perturbation has a reactivity

worth of about 50 and results in a delayed supercritical

transient which was followed for 400 milliseconds. The

computation was performed with unity time scaling and took

two integration intervals of 200 milliseconds each. The

intervals each contained 20 sample points and required 35

and 16 iterations respectively to converge. The total

solution run time was 6.5 minutes on the AD 80/IBM 1800

hybrid system. At $25 per hour the computation cost'was

$2.70. A fairly efficient all-digital code using direct

Euler integration was used to obtain a "theoretical" solu-

tion to the same problem using the same nuclear model. This

solution is shown by the dotted line in Figure 9. The

hybrid solution agrees quite well with the all-digital

result. The digital solution required 3.5 minutes on the

IBM 1800 and, at $15 per hour, cost $0.90. These figures

are summarized in Table 1.

- 60 -


One-group Turkey Point Example
(24 Space Points)

Hybrid All-digital
Euler Integration

(Minutes) 6.5 3.5

(dollars) 2.70 0.90

System AD 80/IBM 1800 IBM 1800

Example 2: A Two-group Model

The results of a two-group hybrid computation are

shown in Figures 10 and 11. The reactor model is that used

by Yasinsky and Henry in their classic paper on spatial

kinetics.[1] The example is a 240 cm thick, three-region

slab reactor experiencing a perturbation in the left-most

one-quarter of the core. In the example shown eleven space

points were used across the core. This is the number used

by Yasinsky and Henry. The results are virtually identical

to those obtained when the same problem was also solved with

up to 50 space points.

- 61 -





3 6 9
---1-l-- ------J-------j--j---
.2 .4 .6 .8 1.0

Figure 10. Two-group example-Flux vs. distance.



IT I r 1I E





- 62 -


Y. Q H.

10 uA, A
10 ,. POINT?


i02 3

10 I


I 0 ---1---a----+---+--+--+-+
0 2 4 6 8 10

Figure 11. Two-group example-Flux vs. time.

- 63 -

The reactor is initially critical and at steady

state. At time zero the fast and thermal group fission

cross sections in the left one-quarter of the core are

increased by a step of +9.6 per cent. The step increase

is followed by a ramp decrease to -9.6 per cent in 10 mil-

liseconds. The perturbation has a maximum reactivity worth

of about $8. The transient is thus highly prompt super-

critical. The calculation was performed using x 1/100

scaling. Twenty integration intervals of 50 milliseconds

were used. The solution was sampled five times during each

integration interval. Each interval required on the average

nine iterations to converge. Again, an all-digital Euler

integration code was used to obtain a "theoretical" solu-

tion for comparison with the hybrid result. The all-digital

result is shown in the figures by the dotted line. The

WIGLE code was also used for comparison. All three solu-

tions agree very closely.

The hybrid computation with two prompt neutron

groups, two delayed neutron precursor groups and 11 space

points used 20 integration intervals with an average of

nine iterations for each one. The solution required 9.6

minutes on the AD 80/IBM 1800 hybrid system. Time on the

hybrid facility costs $25 per hour so the hybrid computation

cost is $3.80. The same problem solved by an all-digital

explicit Euler integration scheme required 0.65 minutes on

- 64 -

the University of Florida IBM 360/65. At $500 per hour,

the solution cost is $5.40. Using the WIGLE code with one

delayed neutron group and a five-microsecond time step

size, the solution required 1.44 minutes on the IBM 360 and

thus cost $12.00. The results are summarized in Table 2.


Two-group Yasinsky and Henry Example
(11 Space Points)

Hybrid All-digital


(Minutes) 9.6 0.65 1.44

(dollars) 3.80 5.40 12.00

System AD 80/IBM 1800 IBM 360 IBM 360

In their calculations Yasinsky and Henry ignored

the effects of delayed neutrons. This may have been done

as an economy measure since their work was done in 1965 when

only slower second-generation digital computers were avail-

able and since the run time and hence cost rises at least

linearly with the number of equations included at each space

point. Since more equations may be added to the hybrid

simulation at each space point with no additional computation

- 65 -

time cost, hybrid solutions were obtained using two delayed

neutron groups as well as with the delayed neutrons ignored.

There is a marked difference between the two models. In

both cases the hybrid solution agrees well with the digital

solution of the same model.

Example 3: A Two-group Model with Temperature Feedback

The results of hybrid computation of a model of a

control rod ejection accident are shown in Figures 12, 13

and 14. The reactor model contains two prompt neutron

energy groups, one delayed neutron precursor and one fuel

and one coolant temperature equation'at each space point.

Feedback arising from the fuel temperature, Doppler effect,

and coolant density changes is included. The reactor is a

light water breeder reactor with a fairly fast neutron en-

ergy spectrum. The core has seven regions in slab geometry

and is 300 cm. thick. Thirty-two space points were used.

A control rod is removed from the left-hand side of the core

over a period of .5 seconds. This is modeled as a ramp-wise

decrease in a in this region. The total worth of the
control rod is $9.34. The solution was followed for 600

milliseconds. The example is one used by J. B. Yasinsky

in a recent comparison of spatial and point model kinetics.[481

Time scaling of x 1/10 was used together with 48 integration

intervals of 125 milliseconds, each containing five sample

- 66 -

points. As many as 58 and as few as 6 iterations were

required for the solution within an interval to converge

depending upon the severity of the flux changes during that

interval. The example is a physically realistic problem.

It is typical of one that reactor designers must solve as

part of a reactor safety analysis. The problem fully dem-

onstrates the capabilities of the hybrid system.

The hybrid solution differs substantially from the

results obtained by Yasinsky using the WIGL 2 spatial ki-

netics code. The data given in Yasinsky's paper is sparse

but should be sufficient, given some ingenuity, to recon-

struct the complete problem definition and the appropriate

initial conditions. Despite this, the solutions do not

agree. Although Dr. Yasinsky was interested and sympathetic,

he had discarded any further information he once had con-

cerning the example. Unfortunately the WIGL 2 code is not

available at the University of Florida at the present time.

To verify that the problem being studied with the hybrid

technique was indeed the same as that considered by Dr.

Yasinsky, an all-digital spatial kinetics code, ADCAL 3,

was written to solve the problem digitally. This code uses

the same two-group nuclear model and the same thermal-hy-

draulics model as that used by the HYKIN 3 hybrid code and

the WIGL 2 digital code. It uses a simple Euler integration


- 67 -

The ADCAL 3 results are shown in Figures 12, 13

and 14 along with Yasinsky's results and the results of

the hybrid computation. The hybrid and ADCAL 3 results

agree reasonably well. They do not agree with Yasinsky's

data. This indicates that the problem studied with the

hybrid and ADCAL 3 codes is not precisely the same one

that Yasinsky solved. One or more of the parameters of

the nuclear or thermal-hydraulic models differs from those

used by Yasinsky. A typographical error in Yasinsky's

paper or a misinterpretation of his problem definition

could easily cause this difference. All the problem vari-

ables were initially at steady state and the spatial shapes

obtained agree with those presented by Yasinsky. Because

of this it is suspected that the error lies somewhere in

the feedback terms that dictate the system's time response.

The coefficients of the power equation which relates flux

level to the time behavior of the fuel temperature, the

relationship between fuel and coolant temperature and the

feedback coefficients which relate the temperatures back

to the time behavior of the flux are all possible culprits.

Regardless of the disagreement with Yasinsky's

results, the ADCAL 3 code provides a digital solution with

which to compare the hybrid results. This comparison shows

that the hybrid results are not exceptionally good for this

example. The thermal flux shape function is shown in

- 68 -

Figure 12. The function obtained by hybrid computation

with the HYKIN 3 code agrees quite well with the results

from ADCAL 3. Reactor power as a function of time is

shown in Figure 13. The hybrid results underestimate the

peak power by 16 per cent. The hybrid solution also

appears to be about 25 milliseconds early. The average

fuel temperature in the reactor center is shown in

Figure 14. The hybrid results underestimate the peak

fuel temperature by 13 per cent.

The cause of the poor accuracy in this example

is suspected to be the scaling of the flux and temperature

equations. During the course of the problem the flux

changes amplitude by over a factor of 100 while fuel

temperature changes only a factor of two and coolant

temperature changes less than a factor of 1.1. If the

scaling is set such that the flux is ten volts while the

temperatures are 100 volts at time zero then at the flux

peaks the flux is 100 volts and the temperatures are less

than ten volts. This scaling is marginal since the flux

is very sensitive to temperature through the feedback

equations. By separating the temperature equations into a

constant term and a time dependent change term and only

solving for the change term on the analog the accuracy is

improved. Without this scheme the problem is unsolvable.

Unfortunately even this does not appear to be enough.

Because the change equations contain the constant terms

- 69 -




270 300

0 5 10 15 20. 25 30


Figure 12. Example three-Thermal shape function.

0 30



J 1.0



- 70 -

20 f

E- \1

15 II II
o1 I i \

o -
o I I1
o I

o r 1 I /
o 0 I I I

5k I


0 .1 .2 .3 .4 .5 .6

Figure 13. Example three--Reactor power vs. time.

- 71 -







.1 .2 .3 .4 .5 .6 .7


Figure 14. Example three-Fuel temperature vs. time.

- 72 -

the scaling will result in very small voltage values for

the temperatures at the flux peaks and very small voltages

for the fluxes at other times.

Time dependent absolute scaling is needed. The

scale factor between flux and temperature ought to be

changed as the flux rises and falls. This would allow

both flux and temperature voltage values to be kept

between 50 and 100 volts during the entire solution which

would improve accuracy for both flux and temperature.

Unfortunately, implementation of this idea would require

digital pots or switchable amplifier impedances, neither

of which are available on the University of Florida Nuclear

Engineering Sciences hybrid computer.

The total hybrid run time was 4.0 hours. At $25

per hour on the hybrid facility,.this cost $100. The

ADCAL 3 Euler integration code is simple and straight-

forward. It is also slow. Run time was 30 minutes on the

IBM 360. At $500 per hour, the ADCAL 3 solution cost was

$250. The WIGL 2 run time was estimated at 3.65 minutes

on a CDC 6600 from data supplied by Dr. Yasinsky. The CDC

machine is generally considered to be four times faster

than the IBM 360/65 so the estimated WIGL 2 cost is $215.

These numbers are summarized in Table 3.

- 73 -


Two Groups with Feedback J. B. Yasinsky Example
(32 Space Points)





Run time
(Minutes 240 3.65 30

(Dollars) 100 125 250

Machine AD 80/IBM 1800 CDC 6600 IBM 360

As can be seen from the examples above, the hybrid

method does work and appears to be well suited to the so-

lution of reactor kinetics problems. It offers quite

reasonable accuracy at costs that are competitive with

all-digital methods. The method does have limitations,

however. These are discussed in the next chapter.

__ _____i__



The results of the work at the University of

Florida indicate that the hybrid method is well suited to

the solution of the kinetics equations but certain limita-

tions are imposed by the characteristics of nuclear problems.

Specifically, the means of data playback, the computation

time, and the method's accuracy are affected by the type

of problems to be solved. The results are summarized below.


A standard method of data reconstruction used in

hybrid methods is zeroth-order interpolation with the

playback one-half cycle advanced.[46] This form of play-

back is not adequate for the class of nuclear reactor

problems considered in this study. There are two reasons

for this. First, the fluxes at adjacent space points are

very closely coupled in reactors. Hence, the solution at

a point is strongly affected by the shape of the playback

at the neighboring points. The problem becomes more severe

as the space mesh is made smaller. The second reason for

the inadequacy of zeroth-order playback is the fact that a

- 74 -.

- 75 -

fairly slow sample rate was used. In order to maintain

maximum flexibility and ease of modification in the hybrid

program the coding was done in FORTRAN using general purpose

hybrid subroutines. Because of this there is a minimum

"sample-update playback" cycle time. It is about 5 milli-

seconds for one group, 10 milliseconds for 2 group, and 15

milliseconds for 2 group with feedback. The problems were

scaled to allow the analog to run as fastas possible in

order to minimize the total computation time. This means

that the flux level can change appreciably between two

samples. The combination of large flux changes and close

coupling between space points makes zeroth-order playback

inadequate for all but very small numbers of space points.

This point is illustrated in Figure 15.

Fortunately, first-order interpolation between

data points is easily implemented on the analog computer.

At each time sample point during an integration interval

the digital computer outputs the value and also the slope

of the fluxes at the neighboring space points via two

DAC channels. An analog integrator uses these values to

reconstruct ~J+1(t) and Qj-l(t) as series of straight

lines between time sample points. Furthermore, since the

DAC channels that provide flux and precursor initial condi-

tions are not used during the integration cycle when the

interpolation is needed, first-order analog playback may be

- 76 -










r ... _. .I

T 1

Figure 15. Data playback.

- 77

added without using any extra interface hardware. Use

of "open loop" interpolation-that is, of simply supplying

the new slope at each sample point-is not sufficiently

accurate since errors accumulate. It is necessary to supply

both the point and the slope.


The hybrid method is sensitive to offset or gain

differences between the analog computer and the hybrid

interface devices. Random error and the coarse (12-bit)

resolution of the ADC do not affect the accuracy but offset

or gain errors which are not random have a strong effect.

The reason for this sensitivity is the iterative nature

of the solution. A small error repeated many times will

accumulate. Figure 16 illustrates this sensitivity. Some

possible sources of nonrandom offset and linearity errors

are the ADC and the DAC, the analog multipliers and analog


By checking the analog components occasionally for

correct operation and by adjusting the ADC and DAC offsets

and gains frequently or by using the digital computer to

compensate for these problems, errors due to these sources

can be kept within tolerable limits. It is readily possible

to keep the total error in the solution below 5 per cent.

On the University of Florida Nuclear Engineering Hybrid

- 78 -

N r"
N U lork f m'-.




Figure 16. Sensitivity to offset error.





U~ t ___

- 79 -

Computer it was found that the ADC offset would drift

up and down by one or two bits from day to day. Although

this is only .050 volts out of 100 volts, the high sensi-

tivity to offset of the iterative method made even this

small error unacceptable. By measuring and compensating

for the offset before each run, the problem can be elimi-


Aggravating the sensitivity problem is the fact

that the time derivative of the flux in the reactor kinetics

equations is defined as a large multiple of the small dif-

ference between very large terms:

S= v(V DVD Za + (1 )VZ ff + XC)

A small error in the production, absorbtion, or diffusion

terms can make a large difference in the calculated flux

derivative. For less sensitive equations the fact that

the hybrid computer is not ideal would have much less ef-

fect upon the solution accuracy. Even with the very sensi-

tive flux equations, the realities of the hybrid machine

may be accounted for quite easily.

Run time

No definite numbers can be given for the computation

time due to the iterative nature of the hybrid method. The

run time is a direct function of the number of iterations

- 80 -

necessary for convergence. This depends upon a number

of variables: the number of space points, the length of

the integration interval, the degree of coupling between

space points, the severity of the perturbation being studied

and the initial guess used to start the iteration procedure.

As was pointed out in Chapter II the iteration

process may be represented as an exponential decay with

the number of iterations done of the error function E(r,t).

The initial value of this error is the difference between

the true solution and the initial guess used.

g (r,t) = qt(r,t) c1 (r,t)


Sp(r,t) = (r;,t)eVP


v is the rate of convergence and

p is the iteration number

When Ep+l p is smaller than some convergence

criteria, p+l is assumed to be small enough and the itera-

tion is considered to have converged. For a particular

set of It and 0 and a given convergence criterion, the

error reduction necessary for convergence is eVP. Thus

vp is a constant. Computation time is proportional to

the number of space points to be iterated, N, multiplied

- 81 -

by the number of iterations necessary for convergence.

For vp constant the number of iterations is proportional

to 1/v. Thus computing time, T, is proportional to N/v:

T N/v.

From Chapter II,

v = n A = -2 An(- cos( ))


T cc

a N
-2 n(-Z cos( ))


cos() 1 (N)2

N 2N
-2 An( ) + (-n(1- 1)


nn(1 l ) T 2


T cc N
T (2 2 n(
N Jai.

For the method of successive displacements

I1A = 2 1
l = I 2 +x (E + P (1--f)v)
=2+ a v -

- 82 -

and Ax = L for one-dimensional geometry so finally


N2 + (E + (l-r)vE
N2D a v .

It is difficult to pick any general functional

dependence out of this proportionality. Using the one-

group nuclear parameters from Example 1 of Chapter III,

run times were estimated for various N. From this data,

computation time appears to be proportional to N2*6 in

the region of interest. Empirically, run time was found

to vary with N2'3 for this example. These results are

shown in Figure 17. Run time for Example 2 of Chapter

III, which is a completely different nuclear problem,

was also found experimentally to rise as N2"3, For

comparison, the computation time for an all-digital code

using an explicit Euler integration scheme will rise in

proportion to the cube of the number of space points.

Run time for an efficient implicit few-group code such

as WIGLE may rise as slowly as N10. However, as more

energy groups or other equations are added to these codes

the coefficient matrix of the computation equation becomes

less sparse and the dependence of run time upon N will

become stronger. In the worst case it could be proportional

to N2.0

-- 83 -






5 10

20 50













- 84 -

These computation time estimates hold only for

one-dimensional calculations. The rate of convergence v

depends upon Ax, the number of computations to be performed

upon N. Thus the total computation time will be proportional

to NLo' and to ( For two- or three-dimensional sys-

tems the same T N2*3 functional dependence is expected to

hold. As estimate of the computation time for a two- or

three-dimensional problem, however, cannot be made from one-

dimensional run times by simply multiplying by the ratio of

the number of space points to the 2.3 power. Rather, if Ax

and the other nuclear parameters are the same, the two- or

three-dimension computation times may be estimated as simply

N 1.
T = T .(N2D)
2 ID J D N

Unfortunately most nuclear problems are rather

closely coupled in space and hence are slowly convergent;

they require a large number of iterations. As the -number of

iterations goes up, the solution time rises and the solution

accuracy falls. In order to improve upon the hybrid method

the number of iterations necessary must be reduced. This is

discussed further in Chapter V.

The total computation time is reduced as the length

of the integration interval is reduced, which agrees with the

results of Howe and Hsu.[46] This is explained by the fact

that the error in the initial guess made by extrapolation of

- 85 -

the previous solution becomes smaller as the extrapolation

is shortened. The rate of convergence of the iterative

scheme, given by equation (29) in Chapter II, is not

affected by the length of the integration interval. The

total error in the initial guess, however, does increase as

the interval is lengthened. When the initial guess used is

P(t) = GD, constant in time, then the total computation time

remains about constant regardless of the length of the inte-

gration interval. When extrapolation is used a limit is

reached, however, at which run time increases again as the

integration interval is shortened. This is shown in

Figure 18. Shorter integration periods cause fewer iterations

to be required but also decrease the "duty cycle," the frac-

tion of time spent getting useful data from the analog. The

minimum point will depend upon the particular problem being

solved and the particular computer used.

The major benefit of the hybrid method is the fact

that it allows more equations to be added at each space

point, up to the limit of the available hardware, with no

computation time penalty. Thus the hybrid procedure allows

the addition of more energy and precursor groups or more

temperature equations virtually free. In an all-digital

analysis the run time will vary at least linearly with the

number of equations at each point.

- 86 -


400 -

z /

T 300
300 -

0 0/
E -

00 200 300 400

Figure 18. Computation time vs.Data is fontegration interval.
0 problem with 26 space


0 ____
0 100 200 300 400

Figure 18. Computation time vs. integration interval.



The results presented in previous chapters show

that the hybrid computer technique may be successfully

applied to solution of space and energy dependent reactor

kinetics problems. The method works and is easily imple-

mented on a very modest sized hybrid facility. The method's

accuracy is quite acceptable. It has no inherent limita-

tions; it can readily include feedback effects, it is not

limited to any particular geometry or number of dimensions.

The number of energy groups, precursors, temperature equa-

tions and so forth that may be included in the system

being modeled depends only upon the amount of hardware

available. The cost of the method in its present stage of

development is comptative with all digital codes if modest

numbers of space points are used.

The Discrete Space--Continuous Time method is a

known hybrid technique[39,41,45) and the space and energy

dependent reactor kinetics problem is one that has received

much attention.[1,2,5-37, 12,4'44,48] This study is apparently

the first successful application of the hybrid technique

to the kinetics problem. It also appears to be the first

successful application of the iterative DSCT method to

- 87 -

Full Text
xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd