Nuclear reactors calculations with the group diffusion equations on digital computers

PROEFSCHRIFT

T ER VERKRIJGING V AN DE GRAAD VAN DOCTOR IN DE TECHNISCHE WETENSCHAP-PEN AAN DE TECHNISCHE HOGESCHOOL TE DELFT OP GEZAG VAN DE RECTOR MAG-NIFICUS IR. H. J. DE WIJS, HOOGLERAAR IN DE AFDELING DER MIJNBOUWKUNDE, VOOR EEN COMMISSIE UIT DE SENAAT TE VERDEDIGEN OP WOENSDAG 22 APRIL

1964 DES NAMIDDAGS TE 2 UUR

DOOR

WILLEM JACOBUS VAN DE LINDT geboren te Rotterdam

Stellingen

1. Het is mogelijk om als standaard van microgolf energie te gebruiken de kracht uitgeoefend op een metalen voor-werpje dat opgehangen is in een trilholte en waarin de microgolf energie wordt toegelaten.

W. J. van de Lindt, "A note concerning forces and torques on spheroidal bodies in cavities" Can. Journalof Physics 33, 1955; 113-117.

2. Bij het opstellen van eindige differentie-schemata, ter benadering van differentiaal-vergelijkingen, verdient het dikwijls aanbeveling de physische toepassing in het oog te houden.

3. In verband met het toepassen van "PERT" netwerken is het zeer wel mogelijk en nuttig,benaderende methodes te geven voor het economisch gebruik van arbeidskrachten.

4. De methodes bedoeld in stelling 3 kunnen zeer goed gebruikt worden voor de totale kostenberekening van projekten als functie van de tijdsduur van het projekt.

5. Het is waarschijnlijk dat de snelheid waarmede een rood bloedlichaampje zuurstof opneemt (of koolzuur afgeeft), vrij sterk afhangt van zijn uitwendige vorm.

6. In overeenstemming met stelling 2 is het mogelijk om bij het numerieke oplossen van de barotropische vorticiteits-vergelijking

J. G. Charney and N. A. Philips, "Numerical integrations of the quasigeostrophic equa-tions for barotropic and simple baroc~inic flows." Journalof Meteorology, April, 1953.

7. Als de vergelijking van Poisson AV - S wordt opgelost met behulp van punt-overrelaxatie langs een horizontale

lijn, dan is er tenminste een orde waarin men door de punten heen kan gaan, waarbij de convergentiesnelheid /2 maal zo groot is.

8. De orde bedoeld in stelling 7 is zeer geschikt indien het integratiedomein een willekeurige vorm bezit.

9. Er bestaat een eenvoudige meetkundige oplossing voor het volgende probleem:

Een stad C moet voorzien worden van aardgas uit twee bronnen A en B. Pijpleidingen van A en B worden ver-bonden in een punt U en gaan dan als één pijpleiding naar C. De drie verschillende stukken pijp hebben in het algemeen verschillende doorsneden en daarom ver-schillende prijzen per lengte-eenheid. Kaak de totale pijpkosten minimaal.

10. De toepassing van experimentele, orthogonale polynomen op het machinaal herkennen van cijfers en letters dient nader onderzocht te worden.

11. Bij het benaderen van functies met behulp van rationale functies moet men voorzichtig zijn: andere dan de . gewenste oplossing kunnen optreden.

DIT PROEFSCHRWT IS GOEDGEKEURD DOOR DE PROMOTOR PROF. DR. E. VAN SPIEGEL

Acknowledgement Summary Chapter 1 Chapter 2 Chapter 3 Chapter 4 Chapter 5 Chapter 6 Chapter 7 Chapter 8 Introduction

Derivation of the group equations

NumericaI method of solution

The finite difference operator

The group diffusion equations as an eigenvalue problem

The inner iteration

Design of a digital computer program

A method for alternating direction iterations Sunimary in Dutch References

•

1 5 27 31 51 75 91 123 137 139

ACKNOWLEDGEMENT

This dissertation is the result of research carried out by the author as part of his task in the Professional

Services Divisional Group of Burroughs Corporation in Pasadena. He is indebted to them for making available the general means to implement this dissertation and would like to thank Dr E. L. Eichhorn, manager of the divisional group, for his continuous encouragement and interest.

The author also gratefully acknowledges the advice and personal intercession of Professor W. A. J. Luxemburg of the California Institute of Technology.

Finally he would like to thank Mrs. J. L. Lundstedt, Mrs. E. A. Potter, and Mrs. S. A. Groves who devoted con-siderable effort to the typing of the drafts.

rectly from the Boltzmann equation, allowing downscattering into an energy level from all other higher levels. These group equations are then solved in two dimensions, with the help of a finite difference method. To approximate the resulting differential equations by finite difference equa-tions, a nine point star of the Varga type is used. With this method a better neutron balance is obtained and a higher accuracy of the criticality results.

In order to be able to show that a real single dominant eigenvalue exists, certain restrictions have to be imposed on the mesh size.

For heterogeneous reactors a difficulty arises as at points on an interface the source function defining the eigenvalue problem, is not defined. The source is redefined in such a way that the total source in an interface is cor-rect though the source for each group may still be incorcor-rect.

A digital computer program is described, with which it is possible to solve the equations in a domain of any shape and having interior regions of arbitrary shape, provided boundaries and interfaces can be approximated by line seg-ments, coinciding with the mesh lines. An example of a two group problem is given. In the last chapter a computational method is described which makes it possible to use alter-nating direct ion iterations with tape limited problems. An accelerat10n method for the outer 1teration is described.

CHAPTER 1

INTRODUCTION

The numerical calculation of nuclear reactor fluxes, even wlth present-day computers, ls still a very lengthy problem. The fundamental equatlon to be solved is the Boltzmann equatlon glven in Chapter 2 (Equation 1). An in-tegratlon of this equation has been done for several simple configuratlons. For lnstance, in the simplest case there is one space dimenslon and the neutrons all have the same energy. The scatterlng is symmetrie around the space axis but is dependent on the direct ion of scattering. This

prob-lem is then discretlzed, giving rlse to a two-dimensional problem. This problem can be generalized by assuming more than one energy and more than one space dimension. This is the so-called transport model. It is used mainly for theo-retical investigatlons and, wlth the advent of faster computers, might ln the future be used for practical re-actor calculations.

The diffuslon model is obtalned from the Boltzmann equation by assuming that scattering is almost isotropie in the laboratory coordinate system. An elastic scattering model gives an explicit expression for the angular depend-ence of the scattering. Introducing this and retaining only the first few terms of a spherical harmonies development of the independent variabIe gives rise to a set of partial dif-ferential equations. Each equation in this set belongs to a

2

different energy. This is shown in Chapter 2. The solution of these equations in three dimensions, is too large a prob-lem for the average size computer and it is therefore usual to consider only a two-dimensional cross section. Escape of neutrons in a direct ion out of the cross section is taken into account by modifying the material constants appearing in the equation. Instead of solving the time-dependent problem, the lat ter is changed into an eigenvalue problem. This is done by multiplying the fission source by a c.onstant, the reactivity. This constant is then determined with the help of the condition that the equations allow a solution which is stabIe and different from zero. How this is done,

is shown in Chapter 3.

The differential equations are then replaced by suit-able finite difference equations (Chapter 4). By using a nine point star and considering the neutron balance in an elementary cell rather than at a grid point, a higher accuracy for the eigenvalue can be obtained.

The linear algebraic system that arises from replacing the differential equations by a set of finite difference equations is too large to solve by matrix invers10n and an 1terat1on method 1s used 1nstead. A choice has to be made from among the several different methods available. All of the alternating d1rection methods are excluded s1nce the use of code words to describe the shape of the reactor cross sect10n and its const1tut1on would be prohib1t1vely

complicated. (It would probably be necessary to give two sets of code words, one for horizontal, the other for verti-cal sweeps.) The use of a line relaxation method was also considered. Here the maximum gain in speed of convergence per sweep is a factor /2. However, the effect on the total running time of a problem is much less. This is due to the following factors:

a) An implicit method requires more computation. b) An appreciable part of the time per scan is used

for calculating the operator. This time is not reduced here.

c) The inner iteration is only One part of the com-plete calculation. (Especially in the case of one inner iteration per outer iteration is th is significant.)

Consequently, the gain, if any, would be small. The method which was finally adopted was the point

over-relaxation method with over-over-relaxation. The neutron fluxes for each energy group are computed and the next approxima-tion to the source funcapproxima-tion is calculated. An inner itera-tion together with a source calculaitera-tion constitutes an outer iteration. The outer iteration will be done with the Power Method. This will be shown in Chapter 3. No acceleration methods have been used, though running times might be re-duced by such techniques.

4

A disadvantage of these acceleration methods is that it requires the storage of the previous souree. Also it

assumes the existance of a complete set of eigenfunctions which have not yet been proven to exist. While each souree problem, solved by an inner iteration, is a self-adjoint problem, the operator applied to the souree is not self-adjoint. The operator M in Chapter 3, for example, is not symmetrie, and the existance of a complete set of eigen-functions belonging to M cannot be shown.

CHAPTER 2 DERIVATION OF THE GROUP EQUATIONS

The behaviour of neutrons in a medium can best be de-scribed with the help of Boltzmann's equation in a phase space.

For the purpose at hand a neutron is conveniently de-scribed by its radius vector r, its speed or scalar velocity v and the direct ion vector

0

of the velocity. Also used is the energy E of the neutron or the lethargy u. The latter is defined by the equation:

Here Eo is some suitable energy level, usually the highest energy that can occur in a problem so that the leth-argy is always positive. If an elementary volume dr du dO of phase space is considered then the angular neutron den-sity N(r,u,O)dr du dO is defined as the number of neutrons having a radius vector lying between rand r + dr, having a velocity vector lying in the solid angle dO around the di-rection vector

0,

and having a lethargy between u and

u + duo A Neutron moving in a medium will interact with the atoms in this medium. This interaction can be anyone of four types.

First the neutron can approach an atom so closely that it collides and is consequently scattered in a direct ion generally different from the one it originally had. In this

6

type of scattering the laws of conservation of energy and momentum apply.

The other three types of interaction are nuclear re-actions and the following can happen. A neutron is captured by the nucleus af ter which i t is expelled again. If the kinetic energy of the neutron is high enough this may be accompanied by the emission of gamma radiation. This will give rise then to inelastic scattering. If no gamma radia-tion is emitted the scattering is elastic again. In the fourth and last case, a neutron is captured by the nucleus which subsequently becomes unstable and disintegrates. This process is called fission. For a self-sustaining chain re-act ion to be possible it is necessary that some of the disintegration products be neutrons again.

In order to describe these processes, the concept of cross section has been introduced. The microscopic cross

I

section is defined as follows: Suppose a be am of neutrons 2

having I neutrons per cm falls on a layer of material of 2

one atom thickness and with NI atoms per cm. Let C be the number of processes that takes place. The nuclear or micro-scopic cross section ~ for this process is then defined as the number of processes that takes place per neutron and per atom. Thus,

C

~

-

The macroscopie cross section

a

is defined as the number of 3

processes not per atom but per cm of the material, so:

a

N₂ IJ.

in which N

2 is the number of atoms per cm 3

of the material. The Boltzmann equation states the continuity of neutrons in phase space.

To derive this equation, the change in angular neutron density in an elementary volume dr du dTI is considered. The following phenomena will cause a change in th is density:

Neutrons will be moving in x-y-z space and thus be moving into and out of the volume element without undergoing collision. The number of neutrons moving in minus the number moving out, per unit time, is:

- vO grad N dr du dO

Some neutrons having a velocity v in direct ion TI will col-lide with atoms. They may either cause fission or change their velocity and direction. In either casa they will be removed from the phase space element. Other neutrons may collide and be permanently taken up by the nucleus without producing fissions. They are absorbed and lost for the purpose of maintaining a chain reaction. All these proc-esses together give rise to the total removal cross section

at.

The total number of neutrons, removed per unit of time, becomes:

AN 0= a" v N (r, u, TI) dr du dIT t

8

in which:

Similarly, neutrons originally having lethargies and direc-tions different from those under consideration, will be scattered into this last group.

Let f(o',u',n,u) dO' du' dO du be the probability that a neutron initially having a lethargy u' and a direct ion

0'

af ter scattering has a lethargy u and a direct ion

0.

More correctly, "having a lethargy u" should read "having a

lethargy lying between u and u + du". The same is true for direction.

The number of neutrons scattered into the volume ele-ment dr du dO by interaction with atoms in the material

becomes:

<Ir du do

J

du'

f

dO' as(r,u') v' N(r,u'

,0')

f(O' ,u' ,O,u) If a source is present, the Boltzmann equation becomes:

aN

-at - v 0 grad N - at v N

+

J

du'

f

dU' as(r,u') v' N(r,u'

,0')

f(TI' ,u' ,O,u) + S (2.1) in which S(r,u,O) is the source function. If the source is due to fission, i t is for all purposes isotropic. If V is the average number of neutrons produced per fission, g(u)du the probability that a fission neutron has a lethargy lying between u and u + du, and af is the fission cross section, the source function can be written as:

S(r,u,o) 41 .. 11 g(u)

I

dUf

J

dIT' aF(r,u') N{r,u',O') v' 1

The 4v comes in because of the isotropy of the fission neutrons. g(u) is also called the spectrum of the fission neutrons. For diffusion approximation calculations only equilibrium states are considered. In the Boltzmann equation

:~

can then be set equal to zero. The boundary conditions of the resulting integro-differential equation pose some problems. Neutrons are not confined to any region. In the case that the reactor region is enclosed by a convex curve neutrons leaving the reactor will not return at some ot her point and in this case the boundary condition is:

N(r,v,TI) = 0

where

TI

is any direct ion into the reactor. This condition will also hold for concave parts of the boundary provided the radius of curvature is sufficiently large. However, for such boundaries with a small radius of curvature, neutrons leaving the reactor at one point can enter the reactor again at another point. In those cases i t will be necessary to enclose the reactor with a convex curve and consider the region between that curve and the actual reactor boundary inside the domain of integration. Suitable equations hold-ing in this region will have to be derived. However, such regions will not be considered here.

The neutron density function depends on six variables. To solve Equation (2.1) on a digital computer an iteration

10

method would almost certainly have to be used. Even for modern high speed computers this would be an almost impos-sible task, and if it were posimpos-sible, the output would be overwhelming and difficult to use in an efficient way. Possibly averaged values, for instanee over v, could be used. However, in that case it would be advantageous to do the averaging before doing the computations. This is what is done in the multi-group diffusion computations. The lethargy range is divided into a few intervals. Neutrons lying in lethargy range u. - u. 1 are considered to have

1 1+

some suitably averaged lethargy

ü

.

.

Cross sections are

1

averaged over the particular lethargy range. Neutrons losing energy due to collisions with atoms are assumed to have the average lethargy of the range they are in until they have lost so much energy that they fall into the next lethargy group.

Going back to Equation (2.1) it can be seen that in each term the angular density function occurs together with the scalar velocity v. This suggests the introduction of

the neu~ron flux v N ~. Equation (2.1) also shows that

in order to proceed, the function f(O',u' ,O,u) needs to be examined more closely.

f(O' ,u' ,O,u) is the probability that a scattered neutron having a velocity v' in the direction

0'

before scattering, has a velocity v in the direct ion

°

af ter

scattering. The usual theory ignores any inelastic

scattering. It also assumes that the atoms with which the neutrons collide are at rest. Neither is true exactly; however, this is areasonabIe approximation. Neutrons from fission have energies up to 10 MEV. In thermal reactors, fission is caused mainly by thermal neutrons (.025 EV.). Thus, high energy fission neutrons have to be slowed down. This is done by letting them collide with the atoms of a "moderator," i.e. a material suitable for this purpose. The absorption cross section of this material should be as small as possible as absorbed neutrons are lost to the chain re-act ion process. It will be shown in this chapter that some of the lighter elements are quite suitable as moderators. The first excited states of these light elements lie in the region of 1 to 4 MEV. Only the fastest fission neutrons will be scattered inelastically. Furthermore, af ter a few

inelastic collisions they will have insufficient energy left for any further inelastic collisions. For neutrons which are close to thermal, the atoms cannot be considered to be at rest. However, over the remainder of the energy spectrum i.e. from about .1 EV. to 1 MEV. the scat tering can be con-sidered elastic with the scat tering atoms at rest. In the case of collision with heavier elements, their first excited states are at energies of about 100 kEV. More accurate and

2 complete data are given by Weinberg and Wigner

12

To derive the transference function fen' ,u' ,n,u) two coordinate systems are used.l The first is the so-called laboratory system or L system in which the scattering atoms are considered to be at rest. The second one is the center of mass system or C system in which the origin is at the center-of-mass of the colliding neutron and atom. In this last system the equations take a particularly simple form. These equations are then transformed back to the L system. In the C system the scattering is assumed to be iso-tropie. This is true up to the high end of the fission neutron spectrum. For high energy neutrons there is a slight bias in the forward direction.2

Transference Function

Suppose a neutron moves with speed vI towards a stationary nucleus of mass number A. The speed of the center of mass in the L system is:

while the speed of the neutron in the C system before col-lision is equal to its original speed minus the speed of the center of mass:

AV_I V

c = vI - v m A+l

Let va and v

b be the speeds of neutron and nucleus respectively, measured in the C system and af ter collision.

From the principle of conservation of energy and momen-turn i t follows that

and:

To obtain the speed of the neutron v

2 in the L system af ter

collision, the speed of the center-of-mass must be added. The lat ter has remained the same because of the conservation of momentum. In Figure 1

a'

is the angle through which the neutron is scattered. This is the angle in the C system. Let

a

be the scattering angle in the L system.

From the figure it can be seen that: 2 2 A2 + 2 A cos

a'

+ 1 v_{2 vI} (A+l)2 (2.2) cos

a

A cos

a'

+ 1 VA2 + 2 A cos

a'

+ 1 (2.3) FIGURE 1

14

From this last formula it can be seen that if scattering in the C system is isotropie, the scattering in the L system in general is not. Only for large A is cos a ~ cos a'.

Equation (2.2) shows that for cos a' = -1 the energy loss suffered by the neutron is maximal. The ratio of the energy af ter to that before collision is:

a ₍A-I) 2 A+l

For hydrogen A=l and this shows that for a neutron colliding with a hydrogen nucleus, the former can lose all its energy

in one collision. A quantity that will be useful later on is the value of cos a, averaged over all solid angles. From Figure 2 it can be seen that:

cos a _1 J'IT

4'IT 0 cos a 2'IT sin a' da'

-1

~ J

+1

~A2

A cos a' + 1 _{d(cos a')}

+ 2 A cos a' + 1 z

FIGURE 2 y

2 3A

Putting cos a = ~ and cos a'

=

~' it follows from Equations (2.2) and (2.3) that E ET from which: (A+l)2 2 + fJo] (2.4) Equation (2.4) states that there is arelation between the angle of

scatt~ring

a and the energy

ratio~.

Now

'0'.'0

=

~

IE

t i - U

while

~~

=

e 2 The transference function f(IT',u',IT,u) must then contain the factor

ö

[TI'. IT

+ A sinh

~

- cosh u-

2 ']

in which ö is the Dirac delta function.

(2.5)

As was mentioned before, scat tering in the C system is usually considered to be isotropie, or in other words, the probability for a neutron to be scattered into the solid

dw'

angle dw' is equal to

lGr'

All neutrons being scattered through an angle between a' and a' + dal pass through the shaded area on the unit sphere in Figure 2. The area of this ring is: dw ~ 2~ sin 8' d8'. The probability of being scattered through an angle a' with respect to the original direction is then:

~ ~

i

sin a'

dB'

~

-

i

d(cos a')

Apparently, equal intervals of cos a* have equal proba-bilities. Now Equation (2.2) shows that equal intervals of

16

COS

e'

correspond to equal intervals of energy. So equal

intervals of energy have equal probability. The allowable

energy range goes from

E

=

E'

2 (cos

e'

= +1)

to

E =

(A_l~2

E'

1 A+l)

The allowable energy range is then: E

2 - E ₁

=

_(A+l)2 4A E'

(cos

e'

-1)

The probability for a neutron to be scattered from an energy level E' to a level between E and E + dE becomes:

(A+l)2 dE

P(E' ,E) dE = 4A

ET

dE u '-u

so that

ET

=

e duo The

transference function expressed in terms of the lethargy u then becomes:

f (IT'

,u' ,O,u)

=

(:~l)

2 eU ' -u • _u

r;:-O" -0

+ A . h u-u' h u-u' ]

L

S1n --2-- - cos --2-- ; (2.6)

where the factor 2v comes in because

0'· 0

has not yet been

integ'rat·ed over the shacled ring of Figure 2

Diffusion Approximation

The flux defined by ~ = N v is a function of the direct ion

O.

If scattering in the L system were isotropie, ~ would

to try to develop ~ in a series of spherica1 harmonies and, as a first approximation, on1y keep the first few terms of this series. 2

In the general anisotropic case:

~(r,E,O)

in which the ~n are the coefficients be10nging to the .L, m

spherica1 harmonies defined by: Pt,m@ ==

eim0 (sin

e

)-m [(t+m)!J

i

dt-m(cos2

e

_l)t

t! 2] O-m)

!

(dcos

e )

l-m

(2.7)

(2.8)

With the help of this equation the first four terms can be written out. They are:

Po 0

_,

1

e i0 sin

e

~ 2-~ e- i0 sin ö v-I, -1 == 1:7

The fo11owing orthogona1ity re1ation exists:

f

dÜ P_{t ,}

,m' Pl,m

=

(2.9)

(2.10) from which the fol1owing coefficients for the first four harmonies can be obtained:

18 411' _{<Po, 0}

f

<p(r,o,E) dO 41T <Pl,O

f

<p(r,

0,

E) cos 8 dO

T

41T <Pl,l

-f

2-~ <p(r,O,E) -Hl

T

e

The total flux density is given by: <PO =

f

<t>(r,o,E) dO

(2.11)

sin 8 do

With the help of Figure 3 the flux in the x, y, and z direct ion are given by:

jx

jy

jz

FIGURE 3

x

_ f

<p(r,O,E) sin 8 cos

0

dO

-f

<t> (r,

0,

E) sin 8 sin

0

dO (2.12)

_ f

<t> (r,

° ,

E) cos 8 dO Q) rn o C) z sin

8

sin ~---~~--~~~~~-y

Equations (2.11) ean be written as: 41T ct>l,O

=

j

3"

Z 41T ct>1,1

=

-1 [jx

-

i jj] "3 /2 41T ct>1,_1

=

+1

~x

+ i jy]

3"

/2

The serial development (2.7) beeomes:

- - 1 ~

ct>(r,o,E)

=

41T Lct>O + 3jz cos

8

3 ( ' _ i

+ '2 J_x 8 + i sin ~ sin 8)

+

~2

(j _x + i j ) (cos

0

sin

e -

i sin

0

sin

e)

y 1 [ct>o + 3j

e

+ 3j sin 8 cos 0 +

..

.

_{= 41T} cos Z X + 3' J_y sin

e

sin 0 +

...

]

4 1 1T [CPo + 3Ö·J<r,E) +

.. J

(2.13) (2.14)

In these equations the lethargy ean be substituted for the energy. Going back to Equation (2.1) and putting :~ equal to zero, Equation (2.14) ean be substituted for v N and integrated over all solid angles.

The following result is obtained:

...

div j + C7_{t ct>o ""}

20

In this equation:

A second equation can be derived by multiplying Equation (2.1) by the direct ion vector

0

and then integrating over all solid angles. The result is:

i

grad

~o

+

at

J

-J

as f

1 (u'-u) ](r,u') du' (2.16) in which:

f₁(u'u)

-and an isotropic source is assumed: SI -=

J

S 0

dO -

0

To obtain the multigroup diffusion equations it is assumed tha

tas

~o (r, u ') does not vary much as a funct ion of u'. Developing this expression in a Taylor series and re-taining only the first term.

Cfs(r,u') 'J(r,u') Rl as(r,u) 'J(r,u)

Substituting this in the right hand side of Equation (2.16), an expression for j in terms of ~O is found:

....

j - - D grad ~O

in which:

Substitution in Equation (2.15) gives: - div D grad ~o + at ~o

=

U'-U

e du I + So (2.17)

The flux ~ in this equation is still a function of the lethargy u.

To derive the group equations, the lethargy range is divided into finite intervals. The flux belonging to the ith interval is defined by:

~

i

=

fU

i N(u) v du u

i _l

Equation (2.17) is then integrated over the lethargy inter-val. If the material constants are fairly constant over this interval no problems arise. However, if they vary

1 2 appreciably some suitable averaging process must be used. ' The resulting equation is:

i-I Di div grad ~i + f

i ~i

=

_j ~ ₌₁ f _1J . . ~_J . + So (2.18) in which the new constants Di' f

i , fij are respectively the diffusion constant, total removal cross section and scatter-ing cross section averaged in a suitable manner over the lethargy interval.

Boundary Conditions

The boundary conditions pose a special problem. The general condition that no flux enters the reactor from

22

without from any direct ion can not be satisfied by the

diffusion equation. It is usually replaced by the condition that the net flow of neutrons entering the reactor is zero. In other words, the flux integrated over all directions into the reactor must be zero. To obtain an expression for this condition Figure 4 is used. The x-y plane is made to form the interface between the two media

A

and

B.

There are four fluxes to be considered. In the positive z direct ion there are the neutron fluxes jA+ and jB+' in the negative direct ion there are jA- and jB-'

The boundary condition between the two media is:

(2.19)

z

FIGURE 4

To obtain jA+ etc. Equation (2.14> is used, which says:

The total flux, going through a unit area of the interface in the positive z direct ion and coming from medium A is given by: . 1 J21T S1l'/2 J_A+ = 411' d0 [<1>0 + 3 TI.}] cos 9

o

0 sin 9 d9 As j = -D

A grad <1>~

o·T

can be written in its components form and the result is:

+ sin 9 sin 0 :;0 + cos 9

::>]

cos 9 sin 9 d9 Performing the integration gives:

The same can be done for the other neutron currents and gives: 1 jB+ =

'4

<1>0 1 <1>0 j _A- =

'4

₊₂ 1 jB- = -₄ 1 <1>0 + 2 1 êl<1>A DA

az

êl<1>B DB öz

24

App1ying the boundary conditions, Equation (2.19) gives: 1 1 31t>A 1 1 ê3It>B

4"

1t>0 - - D _äZ

4"

1t>0

"2

_{DB iJ.z} A 2 A B 1 1 31t>A 1 1 31t>B 4" 1t>0 +

"2

_{DA äZ} 4" 1t>0 +

2

_{DB az} A B

Adding and subtracting these two equations one obtains:

and

1t>0

B

(2.20)

(2.21) If B is a vacuum jB- is zero and according to Equation

(2.19), jA- is a1so zero which resu1ts in the fo11owing con-dit ion on the outside of the reactor:

It> + 2D alt> -= 0

o

A az (2.22) This expression leads to the so-ca11ed extrapo1ation dis-tance. It can be seen from the condition, Equation (2.22), that if the flux is 1inear1y extrapo1ated into the vacuum, the flux becomes zero at a distance

d -= 2D

A (2.23)

from the boundary. More accurate formu1as than Equation (2.23) can be obtained by using more spherica1 harmonics.2 Usua11y d is taken to be:

The source function S(r,u,m -=

i

11 g (u)

f

du'

I

dO' af

(r,

u ') N (r' , u ' , TI') v ' can a1so be written in a discrete form. This gives:

N

S.

=

X.

~ f. ~.

1. 1. j =1 J J

in which f_i is the fission cross section for group i and X is the spectrum of the fission neutrons. The latter does not depend on the energy of the neutrons causing the fission. In reactor ca1cu1ations it is usual to mu1tiply the source S by the constant factor

~

and then determine this constant with the help of the condition that the solution for t ~ ~

is finite. The fo1lowing termino1ogy is then used: k

<

1 the reactor is subcritica1;

k 1 the reactor is just critical;

CHAPTER 3

NUMERICAL METHOD OF SOLUTION

The differential equations for the ~i'S to be solved for each energy group are usually written as:

o

(3.1)

In these equations, (one for each energy group), i and j refer to the energy group. The constants D, r _i , rij' Xi' f

i are both energy and position dependent. However, they are assumed to be constant so long as the material does not change. They are defined as follows:

Di diffusion constant, r

i total removal cross-section,

rij scattering cross-section from energy j to energy i,

Xi product ion spectrum, f

i fission cross-section,

N total number of energy groups, 11k multiplication factor.

The following boundary conditions must be satisfied: Between two regions, the neutron flux must be con-t inuous • Thus :

28

(cp) I (cp) 11

At an outside boundary, or at a plane of symmetry, the general form of the boundary condition is:

A ocp + Bcp = 0

dn

where n is the normal to the bounda~y, and A and Bare (3.2)

(3.3)

constants. At a plane of symmetry, B = O. The problem as stated 1s a homogeneous problem and we can 1mpose an added cond1t10n on the flux that will normalize i t in a suitable manner. Chosen is the condition that on the average one neutron per unit volume of the reactor will be produced; thus:

1

J

N N

V

( I ; X i I; f. CP.) dv = 1

V i=l j=l J J

(3.4)

where V is the total volume of the reactor.

In the next section a description is given of the method used to solve the above stated problem.

Solution

As scattering from one energy group to another is assumed to be from higher to lower energy, the third term on the left in Equation (3.1) is zero for the equation for the highest energy group. First an initial dlstribution for the CPi is assumed. The zero-th approximation to k is

then found by using the normalizing equation for ~.

1

J

N N

V

(~Xi ~ f. ~!) dv = k

V i=l j=l J J

(3.4a)

The following differential equation for the highest energy group is then solved:

(3.5)

For the next lower energy group the third term on the left of Equation (301) can now be approximated by using the

flux just calculated. For the fourth term on the left of Equation (3.1) the zero-th approximation is used. The equation to be solved is then:

div DN_l grad ~N-l +

r

_{N- l} ~N-l

(3.6)

For the third highest energy group the third term on the left is approximated by using the fluxes already cal-culated, while for the fourth term the previous approxima-tion is used.

Af ter the fluxes for all energies have been calculated in this way, k is found from the normalizing conditiony

One full outer-iteration has now been performed. This itera-tive process is now repeated until the value of k no longer changes significantly. The problem is then solved. The reason for calling the iteration an outer iteration is that

30

I

within each iteration a set of differential equations has to be solved. This is usually done by iteration also. Those iterations are referred to as innner iterations.

It is not immediately obvious that k attains a value that does not change appreciably anymore af ter a finite number of outer iterations. In Chapter 5 it will be shown that the problem as stated is an eigenvalue problem and that the method described converges to a solution. However, first it will be shown how the differential operator'is changed into a finite difference operator.

CHAPTER 4 THE FINITE DIFFERENCE OPERATOR

One technique to solve a differentia1 equation on a digita1 computer is to rep1ace the differentia1 operator by a finite difference operator. This is usua11y accomp1ished by covering the domain over which the differentia1 equation is to be integrated with a net of horizontal and vertical lines. The points where those 1ines interseet are numbered and the differentia1 operator replaced by an algebraic expression, relating the function value at point i to those at the surrounding points. There are in fact an infinite number of ways to do this. Depending on accuracy desired and on the function of the dependent variabIe for which this accuracy is required, a choice has to be made. This will be clarified later on in this chapter.

The use of points other than the eight immediate neigh-boring ones gives rise to difficulties in evaluating the function near boundaries, and for this reason, points other than these eight are excluded. For second order differential equations this is always possible. The desired accuracy can always be obtained by choosing the distance between the horizontal and vertical lines sma11 enough. It wi11 be convenient to number the horizontal and vertical lines sepa-rate1y and indicate a column with the index i, a horizontal row with the index j. A point lying on column i as weIl as

32

•

on row j will be indicated by (i,j). Function va lues will be given the subscript i , j , if they belong to the point

(i,j), (e.g. U . . ). The subscript x or y will indicate 1,J

partial differentiation with respect to x or y. For

in-oU

stance U

x will mean

dX .

In Figure 5 the eight neighboring points of point i , j are shown together with their proper indices. To derive the correct algebraic operator, each neighboring function value is expressed in terms of the function at point i , j and the derivitives at this point. This can be done by expanding the function U in a Taylor series about the point i,j. These expansions are then added together, applying some weighting factor to each of the~. This sum must

represent the differential operator with as high an accuracy as possible. Or in other cases, the accuracy of some func-tion of the dependent variabIe has to be made as high as possible.

i-I j+l ili+l i+ j+l

i-I _.i i,j i+ j

i-I i-I i li-l i+l j - l

In homogeneous regions of the reactor the operator to be approximated is

~A+r in which A is the Laplace operator.

First it will be shown that 9-point finite difference approximations to the Laplace operator form a one parameter

family. Then some special cases will be considered.

As the operator actually used will be derived in a different way, some generality will be sacrificed here, and only the case of equidistant grid lines will be considered. It can easily be shown that the A operator is invariant

under a rotation of the coordinate system. From this it

follows immediately that the points i,j+l; i+l,j; i,j-l;

i-l,j; must have the same weighting factor. The same is

true for the remaining four points which also must have a

similar common weighting factor. The Taylor expansions arè

34 Table I <Pi-l,j+l <PiJ" - <P h + <P h +

!

<P h 2 +

!

<P h2 - <P h2 x Y xx yy xy _ 1/6 <P h 3 + 1 <P h3 1 <P h3 + 1/6 <P h3

xxx ~ xxy - ~ xyy yyy

+ 0 (h 4)

<Pi+l, j+l = <Pij + <Px h + _{<Py h} +

!

<pXXh2 +

!

<pyyh 2

+ <pXy h 2

1/6 h 3 1 h 3 1 h 3 1/6 m h3

+ <Pxxx + ~ <Pxxy + ~ <PXyy + ""yyy

+<ph-<p h +

x

y

+ 1 /6 _""'"' _xxx h 3 _ 1 _{~ <Pxxy} h 3 ₊ 1 _{~ <Pxyy} h3 _- 1/6 _<Pyyy h3

_ 1/6 <P _xxx h3 - _~ 1 <P _xxy h3 - 1 _~ <P _xyy h3 - 1/6 <P _yyy h3

+ 0 (h 4) <Pi,j+l <Pij + <P y h +

!

<P h 2 + yy 1/6 h 3 <Pyyy +

o

(h 4 ) <Pi+l,j <Pij + <P h +

!

<P h 2 + 1/6 h3 O(h4) x xx <Pxxx + <Pi , j - l <Pij - <P h +

!

<P h 2 - 1/6 <P h3 + O(h4) y yy yyy <Pi-l,j <Pij - <P x h +

!

<Pxx h 2 - 1/6 _<Pxxx h3 ₊

o

(h4)

The expression for the Lap1ace operator becomes:

(2+

~) Am = 1

r

m. . + m + m + m

v D " , --rT _h~

L

'"

_{1- ,J+} 1 1 "'_{1+ ,J+} . 1 . 1 "'_{1+ ,J-}. 1 . 1 "'i 1 . 1 _-

,J-- 4 (1 + cr ) ct>..

J '

(

4 • 1 )

1J

in which cr is a parameter still to be determined. The case

cr

=

4 gives the weIl known weighting factors

-1 - 4 -1

-4 +20 -4

-1 - 4 -1

which gives an accuracy

o

(h6) in the finite difference

3

approximation of the equation !:J.U

=

O. However this is not

the equation under consideration here and cr

=

4. does not

have any special preference. It wi11 be shown that the

source prob1em is a se1f-adjoint prob1em and thus can be

solved by minimizing aquadratic function. Using this

principle, a star cou1d be found with cr

=

1:

-1 -1 -1 -1 +8 -1 -1 -1 -1

Of particu1ar interest is the case cr

=

2, giving rise to

the star

-1 - 2 -1

-2 +12 -2 -1 - 2 -1

36

This star gives the ~, integrated over an elementary mesh, with the highest accuracy. This is especially important as the eigenvalue is determined by integrating the flux ~ over the volume of the reactor. The general formulas will be derived in the 'next paragraph. However there are two more important cases. The case cr = ~ gives the well-known 5-point formula

-1

-1 +4 -1 -1

while the case cr

o

gives another 5-point star namely

-1 -1

+4

-1 -1

Though not used here, this case deserves special attention. Solving the problem with an extrapolated Liebmann iteration, this difference scheme gives the fastest convergence. It will be shown later that the speed of convergence is the same as that of a line invers ion scheme.

Finite Difference Scheme

The method used here to derive the finite difference 4

scheme is due to Varga. The idea behind the method is that the group equations are actually balance equations. There-fore instead of mak.ing up the balance of neutrons .in a point, as is done to derive a d.ifferential equation, the neutron balance over a f.inite area is made up by integrating the different.ial equation over this finite area. These integrals

are then appróximated by finite difference methods. By doing this the highest accuracy of the neutron balance is obtained, rather than the highest accuracy of the dependent variabIe itself. This will also give the best value of k.

The four quadrants around point ij are indica ted by Rl'. R

2, R3, and R4 as shown in Figure 6. The material constants D and r may be different and are indicated by Dl' D₂, D

3, D

4, etc. The mesh size may be different in different direc-tions and are indicated by hl' h 2 , kl , k2 · The differential equation is first integrated over the rectangle b-d-g-p:

- I

_Is

N

(r1

j

)r <Pj (Di)r t.<Pi ds + (ri)r <Pi ds =

j=~+l

Is

ds S

(4.2)

i-I ij+l ij+l i+l jrl

k l (x,kl) R4

_Rs

à

k_l [b q P (0, y) --_ .. _-...,{x,y) I S4 S3 I

i

h_{2 c} I 1i.h_l :i,-lö _{P i • .1}

!

_{h i+]~} FIGURE 6 Ih 2 (x,O) hl SI _{f tk2S2} d g i-I j - l Rl i j-l

~

i+l ~-1 1K2

38

Each integral can be split into four separate integrals, namely over the smaller rectangles 8

1, 82, 83, 84• As interfaces are required to coincide with horizontal or vertical grid lines the regions Rl' R₂, R

3, R4, are each

homogeneous and the first integral on the left-hand side becomes:

- J

D b.el> ds

S i

Applying Green's theorem and remembering the boundary condition at an interface: + D3

J

h-p-q del> de + D

J

do

4 (4.3) q-b-c

as the line integral over the common bourtdaries q-P, c-P, f-P, and h-P add up to zero. To obtain the finite difference approximation, these integrals are expressed in the values of the function at the nine grid points. For instance, i t is assumed tha t

~

varies linearly over b-c, namely

Integrating this over y between y

o

and y

(4.4)

i

k

D ₄

J

c _Ön dq, d n

k

b

- 3

q"

1 -

1

,J

,

]

(4.5)

The other parts of the line inte~rals are evaluated in a similar way. The integral

r.t

J

q, ds is evaluated by the

S

tra.pezoidal rule. E.g. for the integral over region 4, the following expression for q, is used:

x +

[(

q"

_1,J+

, 1

(4.6)

This expression can be obtained by (see Figure 6) first interpolating linearly between the points (i,j) and (i+l,j) to find q,(x,O) and then between the points (i,j+l) and

(i+l,j+l) to find q,(x,k

l ). Finally q,(x,y) is found by linearly interpolating between the points (x,k

l ) and (x,O) to find q,(x,y).

Integrating this expression for x from 0 to

i

h

2 and y from 0 to

i

k

l giVéS the ~esult:

(4.7)

The operator on the left of Equation (4.2) can then be written

40

Kt

=

a <Pi-l,j+l + b <Pi,j+l + C <Pi+l,j+l + d <Pi_l,j

+ e <Pi j + f <P _{1+ ,J} . 1 . + g <P . _{1- , -}1 j 1 + h <P _1,J-. . 1

+ k <Pi+l,j-l (4.8)

by collecting the proper terms. In this expression D4 (h2 kl) h 2 k l a

- g

k

+

h

_{+ r 4} 64" 1 2 b D3 Chl _ kl) _ D 4 (3h2 _ kl) 3k l (r3 h l + r 4 h 2)

-S-

_{kI hl 8 k l h 2} + 64 f D 2 Ck2 _ hl) D3 (3k l hl) 3hl

-S-

_{hl k 2}

- 8

hl - k l +

s.r

Dl (h2 k2) (h2k2) - 8 k + h _{+ r l}

-s:r.

2 2 g Dl (3h2 _ k2) _ D2 Chl _ k2) 3k2

-S-

_{k 2 h 2 8 k 2 hl + 64} h D 2 (hl k2) (hlk2)

-S-

k + h _{+ r 2} 64 2 1 k (4.9)

Equation (4.8) can be written in matrix notation K = H<t>

K is the row vector, consisting of elements

Kt,

<t> is the flux in each point K of the grid and is a column vector. H is an M x M matrix, where M is the total number of points in the grid. All elements except nine on each horizontal line are zero. The value of the elements not equal to zero are given by expressions (4.9). These equations are valid in homogeneous areas as weIl as on interfaces.

The integrals on the right of Equation (4.2) are treat-ed in the same way as was the second integral on the left.

c' _{r 3}

_6'4

hlk l d' 3h 2

(

r₄k_l)

=6'4

,r_{lk 2} + e' 9k 2 (rlh2 + r₂h_l) + 9k l (r3 h l + _{r 4h2)} (4.10) 64 64 f' 3hl ( _{64 r 2k 2 +} r 3kl) g' h 2k2 r l

6'4

h' 3k2 ₆₄

~

_{r l h2 +} r₂h_l) k' h lk2 r 2

6'4

42

whence,

(4.11) The Jk'S can be thought of as a row vector, the ~'s as a column vector. The matrix, relating these two vectors will be called Hr. In matrix notation Equation (4.11) becomes

The matrix H is an M X M matrix if M is the total number of points in the grid. All coefficients on a horizontal line except nine are zero. These nine coefficients different from zero are given by Equation (4.10).

Later on use will be made of the matrix ~. This will be the same matrix as Hr' except that the r ' s in Equation

(4.10) are replaced by the

X's.

The same will be true for the matrix Hr •

i j

The difference operators in points lying on axis of symmetry and on outside boundaries need special attention. First there is the case of an axis of symmetry. At this axis the boundary ·condition

~

0 holds, in which n is the normal to the axis of symmetry. As a specific case, let Figure 7 represent the situation in which there is a hori-zontal axis of symmetry. If the contributions of R4 and R3

are set to zero by setting D3 and D4 equal to zero, the first term on the left of Equation (4.2) gives the contri-but ion of Rl and R₂ except for the contribution of the line integral over r-R-q. However, as

~

is zero on this line, this contribution is zero. By integrating the other terms in Equation (4.2) only over the regions Rl and R

2, that is, by setting

r;

and

r!,

r;j and r!j,

X~

and

X!

equal to zero, the correct finite difference operator for an axis of

symmetry is obtained. i-l,j on-O q l' P i , j r i+l,j I 8 1 82 i - l , j - l i , j - l i+l,j-l FIGURE 7

For a mixed (logarithmic) boundary condition, the same constants are set equal to zero. However, this time there is a contribution from the boundary r-P-q. _As

_dn

d<1> _{= -}

_A

B _<1>, this has to be integrated over this boundary. The result is:

44

- D2

S

deP di - D

J

deP di + D 2

B 3eP ij+eP i+l, j

hl =

A-dil

1

dil

8 r-P P-q B 3ePi .+eP. 1 . + Dl

_A-

J ₈ 1 - zJ h₂

.

(4.12) It follows then that the operator for a point on sueh a boundary is obtained by setting the material constants of the outside equal to zero and adding (4.12) to (4.9).

The matrix H, the elements of whieh are given by ex-pressions (4.9) is a symmetrie matrix. This ean be seen as

follows. Eaeh point in the grid is coupled only to its eight nearest neighbors. The eoeffieient h

ki (see Figure 8) in the matrix H is determined by the material constants and the dimensions of reetangle a only. This ean be se en from expression (4.9) for e.

FIGURE 8

Similarly hik is determined by the material constants and dimensions of reetangle ~.

However, these are the same as those of rectangle a

and therefore

The matrix H also had diagonal dominancy, that is, the sum of the absolute values of the off-diagonal elements is smaller than or equal to the absolute value of the diagonal element:

However, essentially non-negativity is not always assured. Even though convergence of the Liebmann iteration does not require this last property, existing proofs for the existence of a real dominant eigenvalue k are based on this property.5 Even though it has not been shown necessary but only sufficient for the existence of a real simple

dominant eigenvalue, it is worthwhile to establish condi-tions for H to be essentially non-negative. Birkhoff and

5

Varga define this as follows:

A matrix P is called essentially non-negative if and only if

46 h 2 k1

_>

r 4 h2k1 8 D4 ( - + _k - ) 1 h2 8 D4 3h 2 k1

_{> 3r}

4 h2k1 (4.13 ) ( -

-

- ) k 1 h2 8 D4 3k l h2

_>

3r4 h 2k l ( -

-

- ) h₂ k_l

at any point in the grid, the matrix H is certainly essen-tially non-negative, except perhaps at logarithmic boundary points. Equations (4.13) show that the grid sizes must be chosen carefully. lf the last two equations are satisfied, the first one is automatically satisfied. From the last two equations i t follows immediately that the ratio

~

should never be larger than ~3. lf a finer mesh is neces-sary this should be obtained by gradually changing the grid size. However, it should be pointed out again that viola-tion does not necessarily mean that the outer iteraviola-tion will not converge.

Accuracy Of The Finite Difference Approximation

It can be shown that any nine-point operator belonging

4

to the one parameter family is accurate to h in a homogene-our region if the mesh is square, and accurate to h3 other-wise. It will now be shown that for the operator derived in the preceeding paragraph, the same accuracy holds for ~~ integrated over a mesh. Using Figure 6,

~

can be developed in a Taylor series around point P. Integrating over b-c gives:

(4.14) and a similar expression holds for

~

integrated over p-h, taking into account the change in sign to obtain the inward o<t>

dil'

Let

At an interface the following identities are valid:

(4.15)

Adding 14 and 13 and using these identities, the following expression is obtained:

2 D4h2kl (4) D3hlkl <t>(3) + D4h2kl <t>(4) - _{(1 3}

+

₁₄₎

=

4

<t>xx

+

~74=-~ xx 16 xxy 2 2 + D3h lkl <t>(3) _{16 xxy} + D4h2kl <t>(4) _{16 xxx} 2 D3h lkl _{<t>(3) +}

o

(h 4) (4.16) 16 xXx

48

It was shown in a previous section that

(4.17) and a simiIar expression for

(4.18)

Expressing the different ~'s in Taylor series around P, and taking into account the three identities Equation (4.15), gives for

13

+ IJ

- (I' + I') 3 4

(4.19) Comparing the two integrals, it can be seen that they agree up to the last two terms. In the case that D3 = D4 and the mesh is a square one, these two terms cancel and

4

the truncation error is of order h . Similar expressions hold for the other line integrals.

The five poin~ star also has a truncation error of 4

order h . However, this is the truncation error in the finite difference approximat1on to the differential equa-tions in each mesh point and not the truncation error in the neutron balance taken over a fin1te area.

The integration over a cell is performed by applying the finite operator in Equation (4.10) with

r

l =

r

2

=

r

3

=

r

₄

=

1. In the case that hl

=

h₂

=

k

l

=

k2, this operator gives

Using a Taylor expansion for ~(x,y) and integrating

over an elementary cell gives

J

~

ds

~

~

ij

h2 + 1/24

~

xx

h4 + 1/24

~

yy

h4

The truncation error in the integration is therefore of 4

order h .

For the matrix operator

Hr

can therefore be taken a diagonal matrix without losing accuracy. It will be shown in Chapter 5 that this makes it possible to proof the existence of real eigenvalue for the nine point operator.

In conclusion, the nine point operator makes it possible to obtain a more accurate neutron balance, or alternatively the mesh size can be chosen larger than with a five-point operator if no higher accuracy is required.

CHAPTER 5 CONVERGENCE OF THE OUTER ITERATION

To show the convergence of the method of successive application of the outer iteration it will first be shown that the problem that has to be solved, is an eigenvalue problem.

The differential operator

(5.1) can be changed into a finite difference operator by means of

the method described in Chapter 4. If this finite difference

operator is denoted by H

i , the equations to be solved are of the form:

(5.2)

In this equation ~ is a vector representing the total fission

source and, Xi is the matrix representing the fraction of

neutrons that has energies falling in the ith energy group.

The second term on the right gives the number of neutrons scattered from a higher to a lower energy group.

The source ~ can be written as:

(5.3)

where fj is the matrix of the fission cross sections or that fraction of neutrons in the jth energy group causing fission.

Suppose i t is possible to find a matrix operator Li such that the equation:

52

(5.4) is satisfied. Substituting this in Equation (5.2) gives:

(5.5)

from which it' follows that

(5.6)

-1

under the condition that Hi exists.

Since there is no scat tering from lower energy groups to the higher'ones

(5.7) Equation (5.7) defines ~ while Equation (5.6) expresses Li' for N-l

>

i

>

1, in the L's already found. Thus, these equations define the Li recursively.

From Equation (5.3) and (5.4) then:

Writing:

N

11 - t; f j Lj j-l the eigenvalue problem

is defined.

It should be pOinted out that Equations (5.6) and

(5.8)

(5.9)

(5.10)

(5.7) are N equations in the unknown Li' T~ey define the Li' s whether there is up-scat tering or not., Even though the

matrix M is not calculated explicitly, it will be shown later that its' definition by Equations (5.6), (5.7), and

(5.9) can be used to prove the convergence of the outer iter-ation. A property of the operator H is that it is symmetrie,

-1

as will be shown later, hence the operator H is also symmetrie.

Now X and the f's are diagonal matrices, while r is an upper diagonal matrix. Since the product of a diagonal matrix and a symmetrie matrix is in general not symmetrie, M is not generally symmetrie. Only in the case of a homo-geneous reactor, that is a reactor consisting of only one material, will the matrix M be symmetrie and in 'this case is

it known that the eigenvalues are real and that the eigen-vectors form a complete system.

In practice, M is not explicitly calculated. Rather a souree is

~

N

X + E

i j"i+l

assumed which is normalized. Then rij Lj) is applièd, to obtain the

the operator souree for the ith group. The Li's are then obtained (see formula (5.6»

-1

with the help of the inner-iteration which effects Hi . Finally K is calculated with the help of (5.9). The re-sulting new source function is averaged over the volume of the reactor, giving an improved estimate of the eigenvalue. This method of calculating eigenvalues is known as the power methode It can be shown that the method cOnverges if K has a single dominant real eigenvalue.6 The existence of such an eigenvalue will be considered now.

54

Group Diffusion Equations As An Eigenvalue Problem

It was shown that if a source function

1 N

'1' ""

k

I; fJ. cpJ.

j - l

(5.11)

is defined, the neutron diffusion problem can be defined as an eigenvalue problem of the form:

(5.12) Several questions involved in this were not touched upon. For instance, at points lying on an interface, the source

function is not defined. In genera 1 the fission cross

section will differ for different materials and consequently

'1' is discontinuous at such a point. Integration of the

dif-fusion equations over an elementary cell will not solve all problems and in this chapter a more precise definition of the eigenvalue problem will be given.

When a mathematical model is made of a physical problem, certain characteristics of the lat ter are not necessarily

part of the former. The constant k, as defined for a

physi-cal reactor, must be positive. The eigenvalue k in the

mathematical model mayor may not be positive. The proof of

it being positive is mathematically of interest while

pro-viding further evidence of the sufficiency of the model. In

this chapter a proof will be given of the positiveness of the eigenvalue for some cases that may occur, and the approx-imations made, will be indicated.

In the second part of this chapter it will be shown that each source problem:

- div Di grad ~i +

r

_i ~i

=

S (5.13) is a self-adjoint problem. lts solution therefore minimizes aquadratic function. In the next chapter use is made of this to prove the convergence of the inner iteration. Definition Of The Eigenvalue Problem

As was shown in the previous chapter, the finite differ-ence operator of each group equation was obtained by

integrating it over an elementary cello This results in the following equation:

N N

Hi

~

i " t

J

r

ij

~

j ds +

~

J

X i t f j

~

j ds (5 . 14)

j-i+l cell cell j-l

in which Hi is the finite difference operator defined by Equation

(4.9).

In the homogeneous case, the quantity

N

'l' = t fj ~j j-=l

can be defined at each point and 'l' is continuous.

Equation (5.14) can then be written in the matrix form:

(5.15)

in which Hr and H are matrices obtained from Equation

ij Xi

(4.10) by replac!ng rl' r_{2, r 3 ,}

r

₄ by rij and Xi respec-tively.