MULTIGROUP DIFFUSION METHODS IV

I V

MULTIGROUP

DIFFUSION METHODS

One of the most important uses of numerical methods in nuclear engineering is in the calculation of reactor properties based on the multigroup formulation of diffusion theory. T h e multigroup equations arise as approximations to the well-known age-diffusion theory. It should be noted that the approximations inherent in age-diffusion theory are implicitly contained in the multigroup formulation.

T h e r e are also multigroup formulations based upon the transport theory, and we shall consider such problems in the next chapter. Our present interest will center on diffusion calculations. In this chapter we begin the study with a review of the age-diffusion equations and the adjoint equations. W e then illustrate the construction of the multigroup difference equations and the resultant matrix formulation of the system.

4.1 Age-Diffusion Approximation

For time independent problems the general form of the age-diffusion equations is¹

= ν · [Z)(r, u) V<£(r, u)] - aa(r, u) φ(τ, u) + 5(r, u), (4.1.1) and

<7(r, Kth) = - V · [Z)t h(r) ν φ φ ) ] + af(r) φ φ ) - Sth(r), (4.1.2)

1 See References 1 and 2 for derivation, etc. Note also that the macroscopic cross sections are denoted by σ to avoid confusion with summation operators.

W e shall follow this notation throughout the text.

151

where (r, u) denotes position and lethargy, the subscript ^{e e}t h " denotes the thermal group,

<f>(r,

u) is the flux of neutrons,

q(r,

u) is the slowing down density, Z)(r, u) is the diffusion coefficient, aa( r , u) is the macroscopic absorption coefficient, and S(r, u) is the number of neutrons made per unit volume per unit lethargy. T h e flux is the total distance traveled by neutrons in a unit volume per unit time; the slowing down density is the number of neutrons crossing the lethargy u per unit volume per unit time. T h e left-hand side of Eq. (4.1.1) thus represents the rate of increase of the slowing down density with lethargy. T h e terms of the right-hand side represent, respectively, the rate at which neutrons leak into a unit volume, the rate at which neutrons are absorbed per unit volume, and the rate at which they are created per unit time.

T h e model adopted for the slowing down process is the following:

neutrons are born with lethargies greater than zero, but less than that of thermal energy, i.e., 0 < u < Wth > where wth is the lethargy at thermal cutoff. T h e slowing down equation applies between u = 0 and u = wth · A b o v e wth we assume the slowing down density, q, is zero, and include all neutrons of lethargy greater than #th as thermal neutrons of lethargy wth precisely. ( T h i s model is sufficient for many practical cases.) W e take u = 0 at an energy such that no neutrons are produced by any process at an energy above the reference energy, consequently q(r, 0) = 0.

W e now see that the left-hand member of Eq. (4.1.2) is the rate at which neutrons become thermal per unit volume. T h e terms of the right-hand side represent the rate of loss of thermal neutrons per unit volume, the rate of absorption of neutrons per unit volume, and the rate at which they are created per unit volume. Thus, the two equations merely state that the rate of increase of neutrons is minus the rate of loss plus the rate of creation.

T h e source function, 5 ( r , u), consists of fission sources plus extraneous sources. T h e fission source, say £f(r, # ) , is given by

W e assume the fission spectrum, χ ( « ) , independent of the lethargy, u\

of the incident neutron, and we also assume the number of neutrons per fission, v> to be independent of lethargy. T h e fission spectrum is normalized such that

u) = v

X

{u) [f ^du'

^af(r,

u')<Kr, u

¹

)

+ af t h< M r ) ] . (4.1.3)

ο (4.1.4)

4.1 A G E - D I F F U S I O N A P P R O X I M A T I O N 153

Consistent with the earlier approximations, we have

dux(u) = dux(u). (4.1.5)

J 0 ^ 0

T h e coupling relation between the flux and slowing down density is taken in the form

?(r, u) = £rs(r, u) <£(r, U) , (4.1.6)

where ξ is the average logarithmic energy decrement. Equation (4.1.6) is strictly valid only for an infinite nonabsorbing medium. For simplicity we shall use Eq. (4.1.6); the results to be derived can always be modified for more realistic relationships. Since

? ( r , u ) = 0 , m< 0 , w > w t h , (4.1.7) we could modify the integral in Eq. (4.1.3) to have limits of 0 and «th ·

W e shall use the infinite, range for later purposes. T h e thermal fissions could be included in the integral since

< M r ) = du S(u — wf°° th)<£th(r, w ) . (4.1.8)

J ο W e have assumed

Γ Α χ ( « ) = 0 , (4.1.9)

consequently, the thermal source term consists of extraneous sources only.

By using the coupling relation, Eq. (4.1.6), we may express the basic age-diffusion equations entirely in terms of the flux or slowing down density. I n terms of the flux w e have

- * · [DV4>] + σ&φ - νχ du' σψ - νχσψ <£th = Se , (4.1.10) and

- f a e f l r , « » ) - ^v · [ A h V M r ) ] + °-l^h φφ) = .Stn(r). (4.1.11) T h e first term of Eq. (4.1.11) may be written

-ξσ9φ(τ, nth) = ~ Γ Ai'fciKr, u') 8(u' - uth). (4.1.12)

J ο

T h e set of equations (4.1.10), (4.1.11), and (4.1.12) may be written in the matrix form

•_a_

du

—νχσι th

(£σΒ) - ^V ' t^O l +^{V σ}» - "Χ f ' ^σd* > ^u J ο

/%00

- I du' ias8(u' - ttth); - V · [ Dt hV ] + σξ r se"

>tn(r) J

(4.1.13)

T h e operators are understood to operate on the appropriate argument.

T h e set of equations may be written

[ A - νΒ][ψ] = [ S ] , (4.1.14) where A and Β are formed by factoring the matrix above. T h e functions ψ and S are defined as

Ψ = φ(τ, u)

S = Se(r, u)

Sth(r) (4.1.15)

A similar formulation is possible in terms of the slowing down density and thermal flux. W e shall find such a procedure useful when we consider the multigroup equations.

T h e set of equations (4.1.14) represents a coupled pair of integro- differential equations. T h e objective of the multigroup formulation is to reduce the pair of equations to a coupled set of differential equations, which are, in general, easier to solve.

T h e boundary conditions applied to the equations are the usual conditions, namely continuity of current and flux at interfaces, and vanishing of flux at extrapolated boundaries or return current at the actual boundaries. For different lethargies, the extrapolation distances are not all equal. W e shall assume, however, that the extrapolation distance is a constant for all lethargy.

W i t h the given boundary conditions, the equations are now defined.

In the absence of extraneous sources, the equation (4.1.14) becomes

Α ψ = νΒψ , (4.1.16)

which is a generalized eigenvalue problem. For given collections of nuclear data and the geometry of the system, the operators A and Β are completely defined. Solutions of the equation (4.1.16) exist only for certain values of v, i.e., eigenvalues. T h e study of criticality is usually

4.2 A D J O I N T E Q U A T I O N S 155 undertaken by finding a value of ν which is an eigenvalue of Eq. (4.1.16) and comparing with the experimentally determined value of v, say ve . For ν > j /e , the system would be subcritical, etc. Before we outline the computational steps, we first consider the system of adjoint equations.

4.2 Adjoint Equations

In Chapter I we discussed biorthogonal vectors and gave an illustration of the utility of the dual set of vectors. T h e adjoint function (importance function) plays a similar role in reactor calculations. T h e importance function is very useful in multigroup calculations and also in studying perturbation effects, etc.

In order to derive an adjoint set of functions, we review the formation of adjoint vectors in matrix algebra. For the matrix eigenvalue problem

A x , = XfXi , (4.2.1) the eigenvectors of

A^Ty , - = A,y, (4.2.2)

form a biorthogonal set. T h e proof of biorthogonality was found in Section 1.12 from

y J A x . - x J A^ry . = 0 . (4.2.3)

A similar procedure is possible for operators in continuous spaces.

Consider an operator Θ with the corresponding eigenfunction equation

β ψ , = X&f. (4.2.4)

W e shall assume the are all distinct and the complete.² W e define the adjoint operator Θ * from the relation

f d(wo\) [ ψ * Θ ψ . - ψ θ * ψ * ] = 0 , (4.2.5)

J v o l

where thetyf are the eigenfunctions of Θ * with the associated eigenvalues Yj, i-e.,

Θ * ψ * = 7. ψ * . (4.2.6)

2 In many problems of physics, it is not known whether or not the eigenfunctions of an operator are complete. Completeness is frequently assumed nevertheless.

T h e integration in Eq. (4.2.5) is over the entire volume in which the functions are defined. Equation (4.2.5) is merely a generalization of Eq. (4.2.3).

A s a consequence of the definition of the adjoint operator, it is possible to show that the functions and ψ? form a biorthogonal set and further, that the operators Θ and Θ * have the same eigenvalues. T o prove the biorthogonality, we multiply Eq. (4.2.4) by ψ? and integrate over the volume.

(4.2.7) Similarly, multiply Eq. (4.2.6), for the ^th eigenfunction, by and integrate over the volume to obtain

(4.2.8) By subtracting the resultant equations and use of Eq. (4.2.5), we have

(4.2.9) For distinct eigenvalues we therefore conclude the eigenfunctions of the adjoint operator form a biorthogonal set of functions to the original set. Notice that if Θ = Θ * , then the eigenfunctions are orthogonal amongst themselves; such operators are called self-adjoint.

W e can now show that the eigenvalues X{ and yi are equal. T o this end, multiply Eq. (4.2.4) by ψ ? and integrate over all volume. W e find that

(4.2.10) by Eqs. (4.2.5) and (4.2.7). W e assume, without proof, that the ψ * form a complete set and expand Θ * ψ ? in this set

( 4 . 2 . 1 1 )

(4.2.12) Substituting this relation into Eq. (4.2.10), we find that

4.2 A D J O I N T E Q U A T I O N S 157

by Eq. (4.2.10). I f the expansion coefficients be inserted into Eq. (4.2.11), we learn that

Θ*ψ.* = λ.ψ.* . (4.2.13)

Obviously from Eqs. (4.2.6) and (4.2.13) it follows that

y , = A , , (4.2.14) as was to be shown.

W e now consider the construction of the adjoint age-diffusion equa­

tions, using the concepts and formalism introduced in this section. T o this end, we shall have to generalize the previous matrix formalism slightly, in particular, Eqs. (4.1.13). Since we shall find it necessary to integrate over the variables r and w, the explicit dependence of the func­

tions and operators upon r and u must be given. T h e thermal flux and thermal adjoint flux are represented by delta functions in keeping with the model previously introduced.

<£th(r, u) = <£th(r) 8(w - uth),

<t>U*>

^u

)

= ^ * h W § ( « - wth ) .

T h e thermal diffusion equation, in terms of <£th(r, w), is written

- Γ du' ξσ^η' - H th) f l r , u') - Γ du' [V · ( Dt hV ) - a ? ] ^ b ( r , u')

J ο ^J ο

= Sth. (4.2.15)

Similarly, the fission source in the slowing down equation is

vx(u) Γ du' aft h^th(r, u') = vx(u) afth<£th(r) .

J ο

W e define the vectors ψ and ψ * as

φ

* _ [**(*. «οι

(4.2.16)

(4.2.17)

(4.2.18)

T h e functional relation (4.2.5) may now be written in matrix form as

A*h(^r' ^W)J (4.2.19) where (#*)^- is the z)'th element of the adjoint operator. T h e matrix elements a{j are given in Eq. (4.1.13), (4.2.15) and (4.2.16). W e desire to find the corresponding elements ( α * )ί ;· .

T o find the adjoint, consider first the elements alx and ( f l * )n . T h e element allL is defined in Eq. (4.1.13). T h e contribution to the functional form above is then

j dr du φ* αηφ = j dr duφ (α*)ηφ* . (4.2.20)

Consider first the leakage term. W e recall the vector identity

V · ax = aV · χ + χ · Va , (4.2.21) where a is a scalar and χ a vector. W e have

j dr f άηφ*ν- [DVf\ = j dr f du[V · (φ*Όνφ) - Όνφ · V < £ * ]. (4.2.22)

T h e first term in the brackets may be transformed to a surface integral, and, if we require the f u n c t i o n ^ * to vanish on the surface, as the bound­

ary conditions that partly determine this function, then the entire first term vanishes. T h e second term is again expanded using Eq. (4.2.21) to yield

- J

^dr

J

^du Όνφ · νφ* =

- j

^dr

j

^du[V

·

(φϋνφ*) - φν · ( D V ^ * ) ] .

° ° (4.2.23) T h e first term of the bracket again vanishes, since φ is zero on the bound­

ary. Consequently we have

J dr du φ*ν · (Ώνφ) = j dr J°° du φν · (ϋνφ*). (4.2.24) T h e operator V · (DV) is thus self-adjoint.

T h e adjoint to the remaining terms in the element αλ1 may be found in a similar manner. By integration by parts, we find the adjoint to the

4.2 A D J O I N T E Q U A T I O N S 159 derivative term, while we merely use the defining equation (4.2.20) to

find the adjoint to the fission source term. Combining the results yields

( « * ) u = - f" e γη - V · (Z)V) + σα - va^u) f du' X(u') . (4.2.25)

T h e salient features to note in Eq. (4.2.25) are the sign reversal in the derivative with respect to lethargy and the interchange of fission cross section and spectrum as the weighting function for the fission source.

W e shall have further comments on the implications later. Physically, the importance is the net number of neutrons ultimately produced in the system per unit time. T h e importance of a test particle is equal to the sum of the importances of its progeny. Various conservation equations can be written down from this basic concept (see Reference 3).

I n a similar manner, it is easily shown that

( * * ) « = f du[-V · ( DT HV ) + aah] . (4.2.26) J ο

N o t e that the off-diagonal elements are transposed, however. F r o m Eq.

(4.2.19) w e have

j dr j άηφ*α12φχγι = j dr J du φι^α*)*^* > (4.2.27) and similarly for #21 and ( # * )12 · By carrying out the detailed steps, we then have

- f* e ³* ^^{φ U)} - V · [ D V^ * ( r , ii)] + aa^ * ( r , u)

= va^u) f" du' γ{η')φ*{ν, u') + 8(u - uih) (a^x^*h(r) , (4.2.28)

J ο and

- V · [ Z Wt*h( r ) ] + aaVt*h(r) = va\^h Γ </M'x(W')^*(r, W). (4.2.29)

J ο

T h e solution of the age-diffusion equations and the adjoint equations can be obtained (in principle) by a method of successive approximations.

W e consider first the original equations in the form

Α ψ - νΒψ . (4.2.30) Denote the initial approximation as ψ0 . T h e homogeneous source terms

in the equations are evaluated as f⁰⁰

vex(u) du' af(tt')<£0(r, w ' ) ,

J ο

and

with vQ the known number of neutrons per fission. Similarly, the source term for the thermal group is merely

ζσβψ0(Γ, wt h) .

Having obtained the sources, we solve the homogeneous form of Eq.

( 4 . 1 . 1 3 ) . Denote this solution as ψχ. Formally we solve Eq. ( 4 . 2 . 3 0 ) in the form

ψχ = νβΑ "¹Β ψ0 = s νβΟ ψ0. ( 4 . 2 . 3 1 ) I f ψ0 were actually a solution of the homogeneous equations, then the

eigenvalue ν would be given as

ι _ Η ν ο ΐ ) £ * Ψ . Ο ψ .

" Η ν ο 1 ) £ Λ ί ψβψβ

W e can compute an approximate value for ν even though'4»0 may not be a solution using ( 4 . 2 . 3 1 ) . T h e functional ( 4 . 2 . 3 2 ) becomes

ve / φ ο Ο ί Γ ^ Ψ , Ψ χ

" Η ν ο 1 ) /0" Λ ι ψ „ Ψ0'

Successive iterations give successive estimates of ν from the functional

Η ν ο Ι Ϊ ^ Α ι φ , φ , ' ( 4 . 2 . 3 4 ) where/) is the iteration index, and is the value of ψ on the/>th iteration.

As shown in Section 3 . 6 , the value of vejv will approach an asymptotic value which is the largest value of \jv corresponding to the smallest eigenvalue v, which we desire.

T h e use of the adjoint function will speed the convergence to the asymptotic value. First, if the adjoint equations are written in the form

Α * ψ * = ν Β * ψ * , ( 4 . 2 . 3 5 )

4.3 F O R M A T I O N O F M U L T I G R O U P E Q U A T I O N S 161

and if w e assume a trial solution ψ£, w e may carry out an iteration and compute ν according to

ve ί Φ ο 1 ) £ Α ψ * ψ ί «

^ ( ν ο 1 ) / > ψ * ψρ* (4.2.36) T h e value of ν obtained by using the function and adjoint should agree after each iteration, but will not in general. T h u s a check on the accuracy of the calculation is provided. M o r e importantly, the adjoint function is biorthogonal to the original function. Consequently, w e can use the functional

and expect a more rapid convergence. For instance, if ψ * were exact, then regardless of the error in ψρ and ψρ +1 the exact value for ν would be obtained, because ψ * would be orthogonal to an error in ψρ and Ψρ+ι (^{s ee} Section 3.6).

T h e adjoint is thus seen to be useful both as a check of the calculation and to speed convergence to the estimate of the multiplication. T h e procedure outlined above is simple in principle, but the actual calcula­

tions could not be performed except in extraordinary cases. T h e object of the multigroup method is to reduce the computational burden to a more tractable form.

4.3 Formation of Multigroup Equations

T h e multigroup equations are formed by dividing the lethargy space into discrete segments. T h e formation may be based upon the equation in terms of the slowing down density or the slowing down flux. W e shall consider the difference in the two forms shortly. T o begin we consider the slowing down equation (4.1.1) where we shall assume the reactor homogeneous in regions, and hence D and aa are not functions of position within a given region. W e divide the lethargy interval 0 ^ u ^ #th into G subintervals, and w e denote the points of division ug , g = 0, 1, 2, G, such that 0 = u0 < ux < ... < uG = wth . T h e interval ug — ug_x is denoted Aug . T h e intervals are not necessarily equally spaced. W e integrate Eq. (4.1.1) from ug_x to ug obtaining

q(*> Ο - fa, u^) = f° <fo[DV²<£(r, u) - aa^ ( r , « ) ] + f° du 5(r, u).

T h e source function is

S{r,u)

⁼

v

^X

(u)Q(r)

⁺

S

e

(r,u)

( 4 . 3 . 2 ) with

Q(r) = Γ du' a

^t

(u') φ(ν,

« ' ) + a*,Vth(r). ( 4 . 3 . 3 ) J 0

W e now assume any external sources are separable in terms of space and lethargy dependence. Thus, we take

S e( r , « ) = S e( r) Z ( n ) , ( 4 . 3 . 4 ) and compute the coefficients

ΧΒ =

ΑΤ„Γ ^

^{X(M) (4}

·

³

·

⁵⁾

and

y—1

T h e integrated source function is then given as

Γ°

^{du S( r , u) =}

v

Xg

Au

g

Q(r)

+

Z

g

Au

9

S

e

(r).

( 4 . 3 . 7 )

T h e remaining terms of Eq. ( 4 . 3 . 1 ) can be treated in different ways;

we consider an analytic procedure first. W e define the averages

J * W w< £ ( r , u)

<PA*) = and

= — - (

⁴

·

³

·

⁸

)

_ Ju rf«^a(M)^(r,u)

<ra9 = , T · ( 4 . 3 . 9 ) Aug$g

A similar definition for DG applies. Equation (4.3.1) may hence be written

qg(r) - q^t) = DgV^g{t) Aug - σ^β{χ)Δηα

+ v

^Xg

Q{r)

^{Aug +}

Z

g

Au

g

S

e

(r);

(j> = 1, 2 , G ). ( 4 . 3 . 1 0 )

4.3 F O R M A T I O N O F M U L T I G R O U P E Q U A T I O N S 163 T h e thermal group can be added to the above set of equations to give G + 1 simultaneous equations. T h e thermal group is given the subscript G + 1. T h e equation is

-DG+1V^G+l{r) + aa G + 1< £G + 1( r ) = gG(r) + S? ( r ). (4.3.11) Before continuing we must obtain some relation between the average flux in a group and the slowing down density at the group boundaries.

In order to construct such a relation, a number of approximations is possible. W e shall discuss several of the most common presently, but first we consider an alternative form for the equations.

Notice that the average values of the coefficients are computed by flux weighting. A s one proceeds in the computation obtaining more and more accurate approximations for the flux, the coefficients must be recomputed. T h i s is a time consuming process, and we look for formula­

tions that might be simpler to apply. However, the averages themselves as defined above are exact. A n y simpler formulation that requires less computation will be less accurate. Consequently, we are trading accuracy for computational ease.

W e recall that the slowing down density is, in general, a much more nearly constant function of the lethargy than the flux. I f the basic equa­

tions are written in terms of the slowing down density, then the averages may be simplified. In terms of q, the basic equation is

dq(r,u) D(u) aa(w) , λ c/ . / yl - v n

= ΐΜΰ)^{v 9 ( r}'^{u )} - l ^ u j^{i ( r}'^{u )} + ^{S ( r}'^{u )} • ^{( 4}·³·^{1 2)} W e integrate Eq. (4.3.12) from ug_x to ug and note that if <?(r, u) is slowly varying with lethargy and if Aug small enough, then we may approximate the average of a product as the product of averages to an accuracy of 0(Au). W e define

*'Ζβτν

⁽⁴

·

³

·

¹³⁾

Z>, = ± p d u ^ , (4.3.14)

and

(4.3.15)

F r o m the above equations we find the set of multigroup equations becomes approximately

+

v

Xg

O(r)

Aug + ZgAugS,{v\ (g = 1, 2 ,G ) . (4.3.16) T o compare approximations, consider the term involving the diffusion coefficients, Dg and D g . By definition

fj du D(u) flr, u) Q du D(u) <?(r, u)tfas(u)

Β = ⁹~^x = ^g~¹ (4 3 17)

f < du${v,u) p duq(riU)lias(u)

Obviously D g = Dg if and only if #(r, u) is constant over the interval.

I f #(r, u) varies slightly, then the difference in accuracy of the two approximations is small, and we may use the simple second form. Hence­

forth, we shall use the simpler form of the equations as given in Eq.

(4.3.16).

In order to relate the average slowing down density to the value at the end points, several standard approximations are possible, and we consider only two here. T h e approximations are

?7(r) « qg{r) (4.3.18)

and

Φ) *> ΗΦ) + fe-iWl- (4.3.19) That is, the function is constant and equal to its end value (4.3.18), or

the function varies linearly in the interval. By considering the Taylor series expansion, it is apparent that the approximation (4.3.18) is accurate to order Au(dqg/du) whereas (4.3.19) is accurate to order Au(dqg/du — dqg-i/du) and may yield greater accuracy.

By use of the approximation in Eq. (4.3.18), the equations can be written in terms of the flux if we take

?9(r) = ( ^s) A ( r ) . (4.3.20)

W e then have

W W -

°Mr) = - +*M -

^{- 1}

Γ'

du

S(r,

u), (4.3.21)

4.3 F O R M A T I O N O F M U L T I G R O U P E Q U A T I O N S 165 with ag = a^g + (l/ap). T h e source integral is given by Eq. (4.3.7). I f we take q(r, u) constant in each interval, the function Q(r) becomes

Q(r) = Χ β Μ ' ) + * h (4.3.22)

9-1

from Eq. (4.3.9) with

A - t f b * / ^ * ^ . (4.3.23)

T h e expansion (4.3.22) is also accurate to 0(Au) and is consistent with earlier approximations. W e define the function fg as

- It M ^r) +

^o4VtJ

- ^e (r) .

^(4.3.24)

Likewise, for the thermal group, w e define

fo+i = - ( i » . ) c * c ( ' ) - S l V ) · (4.3.25)

T h e entire set of multigroup equations (4.3.11) and (4.3.21) can be written DWJLr) - *A(r) = fJLt); (£ = 1,2,.... G + 1 ) . (4.3.26) T h e truncation error associated with the equations is 0(Au). In con­

sidering the numerical solution of the set (4.3.26), w e remark that similar equations are obtained for the second approximation to the average slowing down density, Eq. (4.3.19) (see problem 6).

T h e computational procedure for using the multigroup equations is straightforward. A l l of the coefficients can be computed once and for all when the energy groups are selected. A trial solution, say ψ0 , is assumed and the source terms fg are computed. T h e basic diffusion equation for each group is solved to yield a new approximation, ψχ . N o t e that the group equations must be solved in the order g = 1, g = 2, . . . , £ = G + 1· W e shall consider the problem of the spatial

dependence shortly; first we consider the formation of the adjoint multi- group equations.

4.4 Adjoint Multigroup Equations

T h e adjoint age-diffusion equations were given by Eqs. (4.2.28) and (4.2.29). W e consider the lethargy interval divided into subintervals and integrate the basic equation (4.2.28) from ug_x to ug to obtain

- W O + * ; _l (r ) = [ ^ ^ W , . ) + f ^ ^ * ( r , a )

+ VP ^ ) 0 * (^{Λ ρ} ) ( 4 - 4 . 1 ) where

ρ* ( Γ ) = Γάη'χ(η')φ*(τ,Η'). (4.4.2)

J 0

W e again assume the reactor homogeneous by regions.

In order to perform the integrations, we must concern ourselves with appropriate averages for the integrands. T o this end, assume g Φ G and further that there is no leakage, absorption, or source. T h e n Eq.

(4.4.1) merely states that

W = ^ * - i W - ' (4-4-3) On the other hand, the adjoint slowing down density would obey the relation

(fr.),"¹ ? * W = (4-4-4)

in the absence of sources, etc. Since the cross sections are not usually constant, w e conclude that the adjoint flux is more nearly constant than the adjoint slowing down density. I n order to use the simplest possible approximations to compute coefficients, w e see w e must deal with the adjoint flux instead of adjoint slowing down density in contradistinction to the basic equations.

W e again define the average coefficients by Eqs. (4.3.13, 14, and 15) as previously computed in the last section. Further, w e define the quantity

(4.4.5)

4.5 M U L T I G R O U P D I F F E R E N C E E Q U A T I O N S 167 T h e equations (4.4.1) become

" Φι = ".[DXfi W - * * * * , ( r ) ] + v i A& * ( r ) + <£t*h(r) SgG .

(4.4.6) Notice that the equation above is similar to Eq. (4.3.10), except for the sign reversal of the left-hand side and the appearance of the thermal adjoint flux in the equation. In carrying out computations for the adjoint flux, it is convenient to begin with the thermal group (index G + 1) and then compute in the direction of decreasing g. T h u s the thermal group acts as a source for the distribution, and we follow neutrons moving down the lethargy scale. T h e r e is no limit to the lower end of the lethargy scale for the adjoint neutrons. However, due to leakage and absorption, the function continually decreases with w, and we introduce little error by ignoring adjoint neutrons of lethargy less than 0.

T h e relation between the average adjoint flux and the flux at lethargy interfaces is taken to be either

Φ*=ei

^{- (4·4·7)}

or

-

⁼

£ + _

( 4 A 8 )

W e consider the first form. T h e adjoint group equations are then

W i - ' A * - i = / , (⁴·⁴· 9 )

with

/„=-—- ^±Q*(r) - &

^8gG , (4.4.10)

OLg OLg OLg

where ag = a^g -f- l/ocg . A n expression similar to (4.4.9) is obtained for the approximation (4.4.8). T h e adjoint source terms may likewise be found as sums over g in terms of the

T h e computation of the adjoint flux is similar to the computation of the flux, save for the direction of integration. Notice that the same coefficients are available for both the flux and adjoint flux. T h i s provides an enormous simplification of the algebra.

4.5 Multigroup Difference Equations

In the earlier sections of this chapter, we obtained the multigroup equations by dividing the lethargy variable and obtaining a set of

simultaneous differential equations. T h e next step in the process is to divide the spatial interval into discrete segments. T h e results will provide us with a set of simultaneous algebraic equations for which the methods introduced in other chapters are applicable.

T h e basic diffusion equation for each group is given in Eq. (4.3.26).

T h e coefficients are assumed a constant within any given region, but may vary across boundaries. T h e differential operator for the leakage may be written

(«·»

ρ = 0, plane;

ρ = 1, cylinder;

ρ = 2, sphere.

W e consider the interval 0 ^ r ^ a and divide the interval into Κ subintervals and denote the interpolation points as rk , 0 < k < K. T h e spacings are denoted Ark rk , and are not necessarily equal.

W e do assume, however, that an interpolation point falls at each interface between regions. W e now integrate Eq. (4.3.26) from *v1/2 to rk+1/2 to obtain

dr _{rk + l/2} ^{— rPD„ -}3r _rk-l/2

= J dr(?a„

+f,J)r" + j

Λ<σΛ+/«)»*·

( 4 - 5 . 2 )

Λ - 1 /2

T h e material properties may be discontinuous at rk , and hence the integral is factored into portions whose integrands are continuous.

Consider the integral on the right from rk to rk+1;2 . W e expand the integrand in a Taylor series and assume the interval sufficiently small that we may neglect all but the first term; thus we truncate to terms of order Ar. T h e integral is then

f^k+1'² dr(a^g + fg) rP * [a,(r+) <f>g(rk) + /g(r+)] i*fc ^ , (4.5.3)

J rk ^Z

where the + sign denotes the value of the factors obtained as r -> r^k from the right. A similar result is obtained for the integral from rk_1/2 to

4.5 M U L T I G R O U P D I F F E R E N C E E Q U A T I O N S 169

rk except that arguments rk will be used to denote left-hand limits. For rk not on an interface, the integral becomes

f*^{+ l /2} άτ(σ9φ9 + /.) r ' = [ag(rk) φϋ{τκ) + fg(rk)] rpk ( Ι ^ λ ψ ± ± ) , (4.5.4)

which is merely the arithmetic average.

T o eliminate the derivative terms, we use the simple approximation 3φ

dr

Φη+ι — Φ*:

rk+l/2 ^rk+l (4.5.5)

of accuracy 0(Ar). Using the difference approximation (4.5.5) in the left-hand side of E q . (4.5.2), w e obtain a 3-point difference equation

ag.k<kg,k+l — ΚΛΦΟΛ + ^cg^g,k-\ = ^wg,k > (4.5.6) where the coefficients are-given by

c„,k = ^r*-"*°°-«-v* , (4.5.7b)

ark-\

= y VrJM) + Ar^UrDl (4.5.7c)

L ^rk+l/2 Dg.k+l/2 . ^rk-l/2 Dg.k-1/2 , ^rk r +\ A V . „ / - \ A V η

°y.k = j ^ r 1 r y iPoVk) Ark + a^rA,) J r ^ J . (4.5.7d) M o r e accurate 3-point relations may be derived by considering the next order terms in the expansion of the integrands of Eq. (4.5.2). Even higher order expansions may be considered, but then the basic difference relation is no longer a 3-point formula. T h e derivation of the difference equation is similar to derivations considered earlier. Notice that to retain a 3-point difference approximation to the second derivative, w e achieved accuracy of 0(Ar). Frequently variable coefficients have the effect of increasing the truncation error for an approximation which would otherwise be more accurate. A particular property of the present derivation is the symmetry of the coefficient matrix, that is agtk_Y = cgk , as is readily shown. Properties of real symmetric matrices that can be exploited practically and theoretically are pointed out in Chapter I .

Further, much less time is consumed in computing their elements. I f we take the usual boundary conditions that φ0 = φκ = 0, then an equa­

tion of the form (4.5.6) applies for each group and for 1 ^ k ^ Κ — 1.

Notice that the coefficients can be computed once and for all, and hence the entire set of equations is determined.

Analogous techniques are applicable to two dimensional reactors.

Details are left for a problem. Similar procedures may be used to obtain difference equations for the adjoint flux.

4.6 Matrix Form of Multigroup Equations

In order to simplify the equations, it is convenient to consider the matrix form of the multigroup equations. W e consider first the one dimensional case. T o this end, w e first reduce the spatial dependence to matrix form and then express the group dependence in matrix form. F o r a given group, the relevant difference equation was given by Eq. (4.5.6).

W e define the flux vector in the gth group, and as

and

Φθ,Ι Φθ.2

Φ9.Κ-Ι.

(4.6.1)

"9,1 '9,2

υ9,Κ-ΐΛ

(4.6.2)

where w e take φ9>0 = φ9>κ = 0 (the extrapolated boundary) and similarly

0,0 — ^Wg,K

T h e set of equations for the £th group are then 0.

with

(4.6.3)

4.6 M A T R I X F O R M O F M U L T I G R O U P E Q U A T I O N S 171 T h e vector depends upon the source terms f^g(r) which may be discontinuous across boundaries. F r o m Eq. (4.3.24), the simplest form of fg is

f _ _ ( f r s ) g - i Φ„-ι _ vXyAug \ψ ι _ ZgAugSe(r) (A

,

where j 8G +1 = σψ.

T h e elements of the vector involve the factor fg at point k evaluated at the left- and right-hand limits. Since the flux is assumed continuous everywhere, the different limits affect only the coefficients. Thus, the vector wg can be written

G+ l ίΤ=1

where is the matrix containing the slowing down factors, R^-^ is the fission source matrix, and Fg the extraneous source vector. Physically, we know that every element of , H g, and is nonnegative for all T h e construction of the matrices is easily carried out. For S(J we have

with

'S.Jtl 0 o st t 92

0 0

rj Ark ( f as)f /_l t; H r £ Ark_x ( g a s) , , _ l t fc

(4.6.7)

(4.6.8)

W e have considered only scattering between adjacent groups thus far.

Similarly,

with

(4.6.9)

* ί γ n f e v r ^ <^r*⁺> + ^r>< - f - τ Μ ^ - ^ •^{s ( 4}·⁶·^{1 0)}

T h e fission sources are

Rg'->g.\ 0 R-a'^a.2

with

g'-*o.k

XgAug

0 0

XgAug

Φ9'.2

J>g',K-\_

VgifaUg.k-

(4.6.11)

(4.6.12)

W e remark that scattering beyond adjacent groups, e.g., inelastic scattering, may readily be included in the above matrices. T h e scattering matrix Sg contains the coefficients for scattering from group g — 1 to g.

If w e include transfer between more than adjacent groups, the scattering matrix is summed. T h e scattering contribution to the vector wg may be written

with Sg^g = Sg , where is defined in E q . (4.6.7).

T h e matrices Sg^g would be of the form

So'-)? —

ag'-*g.l

0 σ(

0

g'-*g,2

0 0

Jg '-*g.K-l.

(4.6.13)

where the crg'^gn are transfer cross sections at space point n. N o t e that if down scattering only is permitted, = 0 for g' > g. T h e cross section for elastic transfer between nonadjacent groups may be found by considering the collision mechanics. I f the mass of the scattering nucleus is known, and if the slowing down density is assumed constant within a group, then the transfer probabilities are readily computed (see problem 9). I n the case of inelastic scattering, the cross section must be known from experiment or computed from a nuclear model, e.g., the optical model. Similarly, upscattering cross sections must be com­

puted from a theoretical scattering model.

T h e matrices and R^-^ are diagonal with respect to spatial indices since neutrons change their group at the point of nuclear interaction.

Casting the spatial dependence of the multigroup difference equations

4.6 M A T R I X F O R M O F M U L T I G R O U P E Q U A T I O N S 173 into matrix form has been completed; we turn next to the group depend­

ence.

Define the extended vectors ψ and w as

Ψ =

and

w

'Ψι ψ2

w2

T h e entire set of equations can be written Α ψ = w with

A =

Χ Ο

Ο A2

Ο Ο

T h e vector w can be written in the form w = - β ψ - vR<\> — F ,

(4.6.14)

(4.6.15)

(4.6.16)

(4.6.17)

(4.6.18) where the matrices S, R, and F are obvious. In the case of no extraneous sources, we have

Α ψ = - ( S + vR)4> (4.6.19) or

( A + S)<|> = - Λ ψ . W e define

and hence

Τ = — ( A + S)

Τ ψ

=

^{ι ^ ψ .}

(4.6.20)

(4.6.21)

N o t e that Eq. (4.6.21) is one of the same form as Eq. (4.2.30). Of course, the operators are significantly simpler in the finite difference approxima­

tion.

T h e matrix Τ has the form

Τ =

χ ο ο S2 A2 Ο

Ο Ο Ο ο

ο ο

ο (4.6.22)

I f we include scattering over adjacent groups, then Τ has the form

Γ A, ο Ο

S² A2 Ο ο

Si->3 s³ A3 Ο ο

L s t - S2- * g

(4.6.23)

T h e off-diagonal submatrices are all diagonal, whereas the Ag are tri- diagonal. It is easily shown that for either form T^{_ 1} exists and is non- negative (problem 10). T h e basic equation is then

(4.6.24) where G is nonnegative.

T h e matrix R represents the fission source terms and may have zeros in certain columns. In particular, the zeros occur for any column not in a fuel region. W e may permute the rows and columns of G and the corresponding rows of ψ such that the equations take the form

with

G = Gn G12

Ο G22

(4.6.25)

(4.6.26) T h e submatrices Gn and G22 are square; i.e., G' is reducible. T h e matrix G' is also nonnegative and hence possesses a largest eigenvalue which is positive. Furthermore, the eigenvector corresponding to the largest eigenvalue has nonnegative components. In other words, there is a smallest real ν for which the reactor is critical, and the corresponding

4.7 N U M E R I C A L S O L U T I O N O F T H E M U L T I G R O U P E Q U A T I O N S 175 flux mode is nowhere negative. If none of the columns of R is zero, then G is irreducible. W e then conclude there is a smallest real ν for which the reactor is critical, and the corresponding flux mode is positive every­

where. Analogous results may be derived for higher dimension problems.

T h e adjoint equations may also be differenced and a similar eigenvalue relation obtained. It is important to notice that the formulation of the adjoint equations in Section 4.4 will produce a matrix eigenvalue equation of the form

*1 = ΰ * ψ * , (4.6.27) but G * will not be the transpose of matrix G as defined in Eq. (4.6.24).

T h i s circumstance occurs because we computed the adjoint of the age- diffusion equations and then formed the difference equations in energy and space, instead of vice versa. T h i s phenomenon is another illustration of the noncommutability of the operation of computing the adjoint with other operations. However, if, in the formation of the multigroup adjoint equations, we have assumed q* to be approximately constant from one group to another, instead of then the G and G * would be transposed.

4.7 Numerical Solution of the Multigroup Equations

T h e basic equations given in the last section may be solved by a variety of methods. T h e usual procedure is iterative, rather than direct inversion. As an example of the procedure, we outline the steps for one possible method.

Based upon the reactor materials and the relevant cross sections, the lethargy intervals are selected and the coefficients are computed. T h e coefficients may be computed by hand, graphically, or by numerical integration, depending upon the information and equipment available.

A n initial estimate of the flux in each group, say ψ° , is made.³ W i t h the initial estimates and coefficients, the vectors w|J are computed using

3 One of the many commonly used ways of estimating the flux in each group is to calculate it for an infinite assembly first. This calculation is performed explicitly very easily by starting from the group of highest energy and working down. The size of the thermal source of fission neutrons can be selected arbitrarily.

Only the fission spectrum is important at this stage of the calculation, and this spectrum is known. By knowing the flux in all groups of higher energy than the one in question, the flux in the group of interest can be computed. The thermal flux so calculated will not in general agree with that assumed, because the com­

position of the infinite system assumed is not critical.

Eq. (4.6.6). Recall that the vectors w£ are the sources for the gth group.

T h e flux in group 1, ψ ? , is now corrected, i.e., we solve the diffusion equation

A ^ , = w ? (4.7.1) for the flux at all points of the reactor. For a one-dimensional problem,

Eq. (4.7.1) might be solved by matrix factorization since A1 is then tridiagonal. ( T h e equations could also be solved by inversion or iteration).

L e t the solution be denoted ψ } . In like manner, solve the diffusion equations

Α , ψ , = w» (4.7.2) in the remaining groups. T h e vector ψ¹

ψΐ = Ψ¹* (4.7.3)

| ψισ +ι . onve ( ζ ι Ψ¹, ί ι Ψ¹) = ( Ψ ° , Ψ ° ) ,

represents the first iteration. It is convenient to rescale the vector ψ¹ such that

(4.7.4) where ζ1 is the scale factor. T h e scaling, sometimes called renormaliza- tion, prevents the solution from diverging. W i t h the new scaled vectors we recompute the source vectors in each group. T h u s we find . W i t h the corrected sources, we again solve the group diffusion equations at all points of the reactor to obtain ψ², etc. After a sufficient number of iterations, say p> then has converged, that is

|ψρ+ι _ ψ Ρ I < e , (4.7.5)

where e is a convergence criterion, chosen by the user. In words, ine­

quality (4.7.5) requires the flux on the p + 1st iteration to differ from that on the pth iteration by less than some arbitrarily chosen constant at all points of the reactor and in all speed groups. In testing for convergence one must be very careful in the application of the above inequality. T h e inequality merely states that the difference between two successive itera­

tions be less than some quantity. It does not state that the error in any particular iteration is less than some criterion. Convergence may be slow, or the difference may have a maximum or minimum. T h e iteration pro­

cess just described is sometimes called the "inner iteration."

4.7 N U M E R I C A L S O L U T I O N O F T H E M U L T I G R O U P E Q U A T I O N S 177 Begin

Select space and lethargy

intervals

Compute group cross-sections

and coefficients

Generate initial group fluxes

Inner iteration

Compute group sources

Solve for group fluxes

Rescale flux vector

Flux vector converged?

No Yes

Outer iteration

Reactor critical?

Yes

No

Stop Adjust g and /or ι geometry

naterials

Recompute group cross-sections

and coefficients

F I G . 4.7.1. Flow diagram for computing the critical size or mass of a reactor.

T h e scale factors ζμ will approach an asymptotic value determined by the state of the assembly. I f ζμ < 1, then the flux is growing with each iteration; the assembly is supercritical. Similarly, for ζμ > 1, the assembly is subcritical. If the critical size of a reactor is sought, then the adjustable parameters are altered according to 1 — ζμ . T h e adjustment of the gross properties is called an "outer iteration." After completing the outer iteration, the inner iterations are begun again. T h e process continues until the growth factors ζρ equal 1, which implies the system is critical. T h e scale factors ζμ are in fact the reciprocal of the multi­

plication factor.

It is convenient to outline the logic of the above process in a flow chart.

Figure 4.7.1 is a pictorial representation of the steps.

It is important to recall that the method just outlined is one of many possibilities. Frequently programs are written for which there is only one inner iteration per outer iteration. Other variants exist, and we refer to some operating codes in the references for more detail (see References

10 through 14).

References

The derivations and physical interpretations for the age-diffusion equations are exhaustively considered in / and 2. A very general discussion of adjoint operators is found in 4. T h e use of the importance function in reactor physics is well discussed in I and 2. Further interpretations of the importance function are found in 3. Additional discussion is found in 5. T h e formation of the multi- group equations is considered in many references. For particularly interesting discussions, see 2, 5, 6, and 7. T h e multigroup difference equations are considered in 5 through 8. T h e development in the text is predominantly from 5. References 8 and 9 contain very rigorous discussions of the numerical solution of the multi- group difference equations. A brief list of particular multigroup codes is given in 10 through 13. Reference 14 is an abstract of many different nuclear codes.

1. Weinberg, A . M . , and Wigner, E. P., " T h e Physical Theory of Neutron Chain Reactors." Univ. of Chicago Press, Chicago, 1956.

2. Meghreblian, R. V., and Holmes, D . K., "Reactor Analysis." M c G r a w - Hill, N e w York, 1960.

3. Robkin, Μ . Α., and Clark, M . , Integral reactor theory; orthogonality and importance. Nuclear Set. Eng. 8, 437 (1960).

4. Morse, P. M . , and Feshbach, H., "Methods of Theoretical Physics." M c G r a w - Hill, N e w York, 1953.

5. Marchuk, G . I., "Numerical Methods for Reactor Calculations" (translation).

Consultants Bureau, N e w York, 1959.

6. Ehrlich, R., and Hurwitz, H., Multigroup methods for neutron diffusion problems. Nucleonics 12, N o . 2, 23 (1954).

7. BirkhofT, G., and Wigner, E. P., eds., "Proceedings of the Eleventh Symposium in Applied Mathematics." A m . Math. S o c , Providence, Rhode Island, 1961.

8. Birkhoff, G., and Wigner, E. P., eds., "Proceedings of the Eleventh Symposium

P R O B L E M S 179 in Applied Mathematics," pp. 164-189. A m . Math. S o c , Providence, Rhode Island, 1961.

9. BirkhofF, G . , and Varga, R. S., Reactor criticality and non-negative matrices.

J. Soc. Ind. Appl. Math. 6, 354 (1958).

JO. Wachspress, E. L., C U R E : a generalized two-space dimension multigroup coding of the 704. KAPL-1724 ( M a y , 1957).

II. Cadwell, W . R., Dorsey, J. P., Henderson, H . B., Liska, J. M . , Mandell, J. P., and Suggs, M . C , P D Q - 3 : a program for the solution of the neutron diffusion equation in two dimensions on the I B M 704. W A P D - T M - 1 7 9 . J2. Brinkley, F. W . , and Mills, C. B., A one dimensional intermediate reactor

computing program. LA-2161 ( M a y , 1959).

J3. Stuart, R. N., Canfield, Ε. H., Dougherty, Ε. E., and Stone, S. P., Zoom, a one dimensional, multigroup neutron diffusion theory reactor code for the I B M 204. U C R L - 5 2 9 3 (November 1958).

14. A m . Nuclear S o c , Math. Comp. Div., "Abstract of Nuclear Codes," ( M . Butler, ed.), N o . 1-80 (1962).

1. Derive the critical equation from the age-diffusion equation for an infinite reactor with thermal fissions only. Interpret the results.

2. Generalize the results of problem 1 to an infinite reactor with fissions occurring at all lethargies between 0 ^ u ^ utu .

3. Write the slowing down equation for an infinite medium in terms of the slowing down density, and derive an integral equation for the solution. Find the adjoint of the integral equation and show that solutions of the integral equations and its adjoint are biorthogonal.

4. For problem 3 above, derive the adjoint to the differential slowing down equation and then integrate the adjoint equation. Is the result the same as the adjoint equation in problem 3 ? Explain.

5. Derive the multigroup equations for a two-group approximation and express the criticality condition as a critical determinant.

6. Show that the multigroup equations obtained by using Eq. (4.3.19) are of the form of Eq. (4.3.26) with the following definitions:

Problems

Λ = - ( ^ s ) g.

and

7. Consider the basic group diffusion equation (4.3.26 in text) in a two-region slab. Derive a difference approximation which is accurate to order (Δχ)² in each region and across the interface.

8. Derive a five-point coefficient matrix for a two-dimensional diffusion problem in r, ζ coordinates. Show that the coefficient matrix is symmetric.

9. For elastic scattering, calculate Sitj in terms of the energies et , e i +, , 1

and €,·+ 1, assuming that the slowing down density is constant within a group, for the following cases:

1. €i + 1 > <X€j

2. C{ > OC€j J OC€j > €i + 1 > OC€j + 1 3. €, > OL€j J €i + 1 < OL€j + 1

4. oc€i+1 < €{ < lx€j ; €< +1 < a e i1 +

5. €{ < (X€j + 1

6. €t < a e , ; + 1 > ocej+1

where

oc = (M — \)/(M + l )²

and Μ is the atomic mass. Suppose the groups are of width €f/€i +1 = V\Q.

Which cases apply to Hf D, C, very heavy nuclides ? 10. Show that the matrix T^{- 1} [Eq. (4.6.23)] is nonnegative.

11. Show that the approximation (£σβ),^_1<7, = ( f ^ . ) ^¹^ - ! leads to a matrix eigenvalue problem

ψ * = νΟ*ψ*

such that the G * is the transpose of the corresponding eigenvalue problem for the flux.

The following sequence of problems provides an elementary discussion of the numerical solution of the reactor kinetics equations.

12. T h e infinite medium kinetics equations for G delayed groups may be written dn(t)

_

^{p —}

dt ~~ Λ

Έ£± = h η ( ί) _ λ , Ο ( ί ) , « = 1 , 2 G ,

where η is the neutron density, C^l the zth delayed group precursor, p the reactivity, Λ the generation time, At the zth group decay factor, j3, the ith group fission yield fraction, and β = Σ , j8, the total delayed yield.

(a) Define the vector

η C¹

Ψ

^{I ^}

C^G