THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING THE HARMONICS METHOD FOR NUCLEAR REACTOR ANALYSIS Charles A. Stevens A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan Department of Nuclear Engineering 19o1 December, 1961 IP-540

ACKNOWLEDGMENTS The author would like to express his appreciation to Professor Paul F. Zweifel for his guidance and advice. In addition, sincere gratitude is expressed to the following: To the Atomic Energy Commission and the Babcock and Wilcox Company for supplying financial assistances during the course of this work. To the staff of the Computing Center of the University of Michigan. To Atomic.Power Development Associates for the use of their AIM program. To my wife for her assistance in preparing this manuscript. ii

TABLE OF CONTENTS Page ACKNOWILEGMENTS.......... o.......... o.. ii LIST OF TABLES.................................................. iv LIST OF FIGURES........................................... o.o. v INTRODUCTION...................................................... 1 CHAPTERS I. THE STATIC PROBE o o....o...................... 3 II. KINETICS.,....... o a........ a o o o o. o.. a. o o. o. o. a o o. 29 ITI. FURTHER PROPERTIES OF THE HARMONICS METHOD................. 41 1. The Eigenvalue k as an Upper Bound.................. 41 2. Convergence.o............... o............... o 44 3. Comparison of Harmonics Method with Finite Difference Schemes......,,,,.........,....... o 52 4. Applications, eo.. o o....... o o..*....... o........... < 57 IV, DETAILS OF THE COMPUTER PROGRAM................., 69 _4 Input. a 4e a 0 o o ~ o o o a 0 e e e ee e e e a e e a a 0 a 69. Card Format............,...,...,.............. 6971 2. Card Format. o, oe o e o o o a **. o e e.,...e ee e e o eeae oe.e, 71 3 Output 0 0... 0.... 0..... oe o... o o o....e........ 73 4. Further Details.................................... 74 5. Timing........ o,,.................................. 76 6. Flow Diagram.. o o o o....... a.................... 76 APPENDIX. COMPUTER PROGRAM........ o o. o o............oo.. e... 81 BIBLIOGRAPHY.. o o............................. o. o 89 iii

LIST OF TABLES Table Page I Parameters for Illustrative Reactor 1................. 42 II keff and Prompt Neutron Lifetime for Illustrative Reactor 1, for Various Values of M............... 43 III Iterates for Illustrative Reactor 2 in Eleven Mode Approximation................................. 47 IV Parameters for Illustrative Reactor 2................. 47 V Harmonic Coefficients of Fast Flux for Illustrative Reactor 1, for Different Values of M............00.... 50 VI Harmonic Coefficients of Thermal Flux for Illustrative Reactor 1, for Different Values of M................. 51 VII k f for Illustrative Reactor 2 for Different Numbers of Harmonics and Mesh Spacings......................... 56 VIII Reduced Spatial Dependent Two Group Cross Section Set.. 65 IX Reduced Space Independent Two Group Cross Section Set.. 66 iv

LIST OF FIGURES Figure Page 1 Iterated Thermal Fluxes for Initial Guesses of k = 0~2, 3.0 for Illustrative Reactor 2 in Eleven Mode Approximation.,................... 45 2 Iterated Fast Fluxes for Initial Guesses of k = 0,2, 3 0 for Illustrative Reactor 2 in Eleven Mode Approximation..................................... 46 3 Thermal Flux in Illustrative Reactor 2 for 7, 11, and 20 Mode Approximations....... o 9....... o......... 54 4 Fast Flux in Illustrative Reactor 2 for 7, 11, and 20 Mode Approximations........... o.................. o. 55 5 Thermal Flux for Goertzel-Wilkins Minimum Mass Problem.. 59 6 Comparison of Two Group Harmonic Calculation (Using Space Dependent Coefficients) with 18 Group AIM Code Having Space Independent Cross Sections................ 63 7 Comparison of Two Group Harmonic Calculation (Using Space Independent Coefficients) with 18 Group AIM Code Having Space Independent Cross Sections..........<....... 67 v

INTRODUCTION The so-called harmonics method was first used for the steadystate nuclear reactor problem by Goertzel and Garabedian (1). Further work using the method was later carried out by Edlund and Noderer (2). In 1959, Garabedian and Leffert extended the method in order to deal with the non-steady state problem (3). Also, at that time Holway (4) added a perturbation scheme to the harmonics method, and was the first to use it for problems having more than one energy group. Furthermore, he was the first to allow discontinuous diffusion coefficients. Ih 1960, Foderaro and Garabedian (5) attacked the two energy-group problem without recourse to the perturbation method of Holway. The same authors also solved the two-group kinetics equation by means of the harmonics method (6), in effect a continuation of the work of Garabedian and Leffert mentioned previously. This thesis is a further development of the harmonics methodo In the first section, only the static case is considered. The present development is felt to be superior to the previous ones for the following reasons: (a) It is applicable to any number of energy groups, Previously a maximum of two has been considered. The computer program which is included with this thesis is, however, limited to a maximum of four groups. The coide is not a production code; it has been written to check the methods developed in this reporto -1

-2(b) Group constants may vary arbitrarily with positiono Previous developments allowed only constant cross sections within a subregion of the reactor. (c) Group skipping is permitted; that is, neutrons may be scattered from one group to any or all groups below ito Also, neutrons may cause fissions in any group, with resulting neutrons emitted into any group. These features are not included in most fewgroup codes in use today. Up-scattering presents no difficulty to the solution of the multigroup equations by the harmonics method, but it has not been included in this work. (d) An iteration procedure, based on a variational principle, is used to find the eigenvalue of the static diffusion equationo This procedure is faster and more accurate than other procedures previously used in conjunction with the harmonics method. (e) The adjoint fluxes are calculated as well as the fluxes. The second section deals with some aspects of the kinetics of nuclear reactors. These are also attacked by the harmonics method, and some of the kinetics calculations are also incorporated into the codeo The third section presents some further properties of the harmonics method with other methods in use todayo Some calculations for which the methods described in this work are particularly well suited are also in. cludedo The fourth section describes the computer program to an extent that the reader may use it for his own problems, if he so desires,

CHAPTER I THE STATIC PROBLEM In the multigroup diffusion approximation to the Boltzmann equation the set of equations to be solved, for one dimensional slab geometry, is K = dr =^1; ()=0 For i. For i = 1, the third term does not appear. The symbols in Equation (1) are defined as 0i = neutron flux in the ith group, Di = diffusion constant for the ith group, aai = total macroscopic removal cross section from the ith group, Di(r + (B)i = transverse buckling for the ith group, ci = macroscopic absorption cross section for the ith group, abs a.. = macroscopic transfer cross section from group j to group i, Xi = fraction of all fission neutrons which are emitted into group i, vafc = rg n the average number of neutrons produced by a fission in the th group multiplied by the macroscopic fission cross section in the ith group, -3

-4N = number of groups, x = distance measured from the origin, k = the eigenvalue, or the static effective multiplication factor. The harmonics method Will be used to find a solution to Equation (1) such that (a) Xi(x) > 0 everywhere. (b) Xi(x) is bounded, (c) 0i(x) vanishes on energy-independent outer boundaries. (d) Xi(x) = 0i(-x), or di(O) = 0 (symmetry about x = 0.) dx (e) Di(x) d/i(x) is continuous everywhere. ddx (f) di is continuous within each sub-region of the reactor. dx Requirement (d) is imposed for conveniences it is not necessary. The eigenvalue k corresponding to the solution satisfying conditions (a) to (e) is the effective multiplication factor of the reactor. This is well known, see, for example, Reference (21). It is also known (29) that this eigenvalue, which will frequently be denoted by keff, is positive, real, and larger in magnitude than all other eigenvalues. If Equation (1) is multiplied by a continuous, bounded function G(x) which vanishes at the outer boundaries of the reactor, and is then integrated over the reactor volume the result is 6() JG ) 0 (X) I- ( ()X) l e (2) + Gi W fG (X) 1(t)0^ 0. j =m

-5The reactor is assumed to consist of R sub-regions, each having Di(x), a.i(X), etc., which are continuous. Points of discontinuity of any of these space-dependent coefficients define the region boundaries. We may carry out the integrations required in Equation (2) by summing the contributions of all the sub-regions. In particular, let us now examine the first term. Carrying out the integration of the first term by parts, there results r= XI'rR (3) Xr In the above notation r represents the region number, xr is the position of the interface between the rth and (r+l)th region. D'(xr) is the diffusion coefficient at the left of xr; D(Xrl) is the diffusion coefficient to the right of Xr,,l The other notation is obvious. The integration by parts performed above is valid only if, within each r, Di (xX)d(x) is continuous and its derivative is bounded, These condx ditions are fulfilled by virtue of Equation (1) and conditions (b) and (e) below Equation (1). The order of summation in Equation ('3) may be rearranged to yield

-6~(* rdr ),W d~ +) 5 5 urica c e.5 R-I ~~~~ - > -(xr)[D;(Y (r+. Ds (Xr)'(Xr)] r=i I e ) (4) D, (x) 6 (X) $ - "7 v L But it has been specified that G(x) vanishes on the surface of the reactor. Hence the first term on the right vanishes. Also, the second term on the right vanishes by virtue of condition (e) under Equation (1). Hence, R = / /G~x) c rYO (X O,) Q R C/ /d " ^ ^ 3r 3 5 ~X. Equation (5) is actually a consequence of Green's theorem and applies equally to cylindrical or spherical geometries. Equation (2) thus becomes WX dx X()u -, ra i N ift C - fR 5W 1 X) X) (a xx ~~ o + -( +^~ ^W^W^W^ d:(6) K ^?)^ <

-7Because of conditions (b), (c), (d), and (f), the fluxes may be represented by an infinite sum of eigenfunctions of the scalar Helmholtz equation, i-BV -O (7) dXZ V which vanish at the outer boundaries and are symmetric about x = 0. The reader should not confuse the physical quantity v which appears as vafj with the subscript v which appears in Equation (7). These eigenfunctions form a complete set for functions with similar properties so there is no question of the validity of such an expansion. Because of the symmetry, it will suffice to concentrate on the interval [0,L], unless otherwise specified. The orthonormal eigenfunctions of Equation (7) are'j W-?CoS(X (8) where B Ti-4 )J[ (9) The continuity requirements on the fluxes are sufficient to allow us to represent their derivatives by a term. by term differentiation of their Fourier series expansions (7). Accordingly, M (*W= a;^ I~(x) (10) CL (x) M V= I

-8In principle, M should be infinite. But in practice, it is necessary to keep only a finite number of eigenfunctions. Until now, G(x) has not been specified. Any iV(x) has the properties required of G(x). Thus let 6(X) = 1 (X)/ = I1,2, M (12) If Equations (10), (11), and (12) are substituted into Equation (6), there results *^.~L ~L' TJ 1. (X Vd( (-X)x (X)c X} hXdX (X) (X) C(X) dx ]o =-' gj A / dl =/, l | / M At this point, it is best to introduce more concise notation. Let L L (!?) E, -~^ S [ D a 6 (X) 10 ^r A ( where 5ij is the Kronecker delta, J-Jt

-9and let L _ 7 = fJ h (x) V$ (x) T wx X (15) o0 Equation (13) now becomes ZXEE2- sv7-=o (16) =-t V- l 1'= 122, ^, N] Mx = 1 M. Equation (16) represents a set of MN homogeneous linear equations with unknowns ai, The coefficient matrix may conveniently be arranged in either of two equivalent forms. First let Eg } v a. v (17) K In what is denoted as the (Q v)ij representation, Equation (16) is written in matrix form as QA= 0 (18) where 11 11 11 12 lN 11 ll ~ Q11'11 Q^..........Q(19) 21 9= | O a 0 0 9 s a a |1 Q 1............. (zg

-10and 1 al A1~~~~~~~~ = (2o) a2 1 * "1 N EM In the (Q ) representation, Equation (18) applies with 11 12 iN 11 12 1N 11ll ill ~" 1 Qll1 2 Q12.. Q1M 21 Qll QN1 Q11 21 Nl NN and 1 a2 1 a2 _a_

-11In the (Qv)i) representation the matrix Q consists of N submatrices, each having degree M X M, so that Q has degree MN x MN. In the (Q j) representation, Q consists of M submatrices, each having degree N x M; again Q has degree MN x MN, In this work, the (Q~v)ij representation will be used exclusively. It is well known from the theory of linear algebraic equations that Equation (18) can be solved for the unknowns at, within an arbitrary multiplicative constant, only if the determinant of the matrix Q vanishes, Accordingly, we must set et ((Q)=c (E-~S) =o. (23) If k - 1, Equation (23) is the critical condition of the reactor. If det(E-S) 0, the value of k which satisfies Equation (23), and yields a, corresponding to Xi(x) which are everywhere > 0, is the effective multiplication factor. Of course, the a t are obtained from Equation (16), once Equation (23) is satisfied. At first glance, it appears that there are MN values of k which will satisfy Equation (23). However, because of the high singularity of the matrix S, only M of these are not equal to zero. This observation will now be proved, Some of the properties of the matrices dealt with will become clear during the proof, which is the primary motivation for carrying it out. First note that S may be written as a product of an MN x M matrix with a M x MN matrix. Remember also that, in principle, M > 0oo

-12xi 1 r 1 1 1 2 N 0~ V^afll Tf 12'' VCaflM VCfll4 Crfl M xL 1 1 X2' | | Vdf2l V2 f22 2 afML S = X2 O 0 ~ ~ ~ v 1 XN 0 oo XN 0 XN XN The rank of the X matrix is obviously equal to,M because it consists of N diagonal matrices, each of rank IM It will now be shown that the vaf matrix also has rank M; in fact it will be shown that each of the N suibmatrices of which it is composed are non-singu3aro An arbitrary submatrix is made up of elements of the form (X0) W (x) ( x) dx. If the (vaf)j submatrices were singular, it would be possible to express the elements of one row, say the nth as a linear combination of the other rows, Thus V | Y D M M

-13But since f(x) is quite arbitrary the above equality holds only if the integrands are equal, which implies M /=l However, this violates the condition that the 4n are linearly independent. Thus, it must be concluded that the submatrices are nonsingular. The matrix S will then also have rank equal to M. Satisfying Equation (23) is equivalent to obtaining the eigenvalues k of the matrix E -S, providing that E exists. E is a lower triangular block matrix. For a moment consider the value of the determinant of E. For the square matrix E, whose degree is n = MN, one method for evaluating det E is led E =L "(v.../ 4)e,,.. eden where the sum is over all permutations vl,.., vn of 1,.., n and C(Ovl.0o, vn) is +1 or -1 according as the sequence vl,'.., Vn has an even number or an odd number ofinversions (8). But it is clear that all permutations which include an element from a submatrix below the diagonal blocks must also include at least one element from a submatrix above the diagonal blocks. But the latter are all null matrices.

14 — Accordingly, the off-diagonal blocks do not contribute to the value of the determinant. Therefore, for the purpose of calculating det E, E may be regarded as a block diagonal matrix. But it is known that the determinant of such a matrix is equal to the product of the determinants of the diagonal blocks (9). Each of these submatrices is non-singular, as can be shown in the same way as was done for the (vaf)j submatriceso Accordingly, their determinants must differ from zero, the determinant of E differs from zero, E-1 exists If S has rank M and E1 has rank MN, E-1S must have rank M (8). The diagonal matrix consisting of the eigenvalues k must be equivalent to E'IS, and it must therefore have the same rank M. Therefore there can only be M non-zero eigenvalues k. This completes the proof. Of the M non-zero eigenvalues, only one corresponds to fluxes which are everywhere > 0, This is the root which is all important in this work, In principle, the static reactor problem can be solved by following the above procedure. In practice, if one recalls that the matrix Q has dimension MN (which may be quite large), the numerical difficulties in finding the value of k which satisfies Equation (23) may be quite formidable, Attention will now be devoted to finding an efficient method of obtaining k, Equation (1) may be written in matrix operator form ~L U (24)

-15where L is a N x N matrix whose coefficients are La (x) - F D(x)1 -< ( )1 -- i ) (25) 0 and X are column vectors consisting of the N components Xi(x) and Xi, respectively. Note that X is not the same as the X matrix of page 12 U is a row vector whose components are vafj(x). An equation which is closely related to Equation (24) is IiL " 1 v K U e(" )(26) Where L*, U*, and X* are respectively the transposed matrices of L, U, and X as defined above, The solutions 0* of Equation (26) are the adjoint fluxes. Equation (26) may be written in the expanded form d D i(X), Q -(X) a x) x ( N T. K(+') h =(27) +-1,j2 /^1 ag, ix) f^or 0 f For i = N, the third term does not appear. Note how Equation (27) compares to Equation (1). If Equation (27) is solved by the same method used to solve Equation (1) the equation which results is 21 - ^S(28) 1/ —,,,,=a

-16j where the bV are defined by ^ i W= XM (29) v~ I Equations (28) and (29) are analogous to Equations (16) and (10), respectively. The i and j indices are interchanged on the E and S matrices. Actually, 1i and v may also be interchanged, for it is obvious from Equations (14) and (15) that the E and S matrices are not affected by an interchange of the lower indices. The interchange of both upper and lower indices corresponds to taking the transpose of the matrix Q, defined by Equation (19). Thus to obtain the solution for the b., it is necessary to solve the matrix equation 9 B= (30) where Q* is the transpose of Q and B is defined by b1 b1 2 bM B3 1 2 b1,N bM For a non-trivial solution to Equation (30), det (Q*) must be made to vanish by adjusting k, Once k is determined, one may solve for B by means of Equation (30). It is well-known that the same value of k which makes det Q vanish, will also make det (Q*-) vanish (l0). Of course,

-17it is just as difficult to obtain keff from the adjoint calculation as it was from the flux calculation, and for the same reason. If Equation (24) is pre-multiplied by 0*, and the resulting equation is integrated over the volume of the reactor, there results an expression for k which is "stationary" (11). 4 (X) ( V(X) (X) dX K = -.. (32) jf(x) L (xYQ(x d If the matrix multiplications required in Equation (32) are carried out, there results ~F ^^ )W,.(X) V N(X) WX r= \ —----- - (33) N z;_ C (x) L,, (x) d J x If trial functions are chosen for 0i(x) and 0i(x), the difference between the exact value of k and the k calculated from Equation (33) will be proportional to second order differences between the exact fluxes and adjoints and their respective trail functions. Let the trial functions for the fluxes and adjoints be represented by Equations (10) and (29). Accordingly,,f(x;)= a)- (rx) (10) V=I

-18and x )= = b (X) (29) Insertion of the above expressions into Equation (33) yields hf 1 M M e e ___ ___ __= j i" l =. - (34) ( 1\/ M M_, i i Of course, one must now determine the av and bv which, when inserted into Equations (10) and (29), yield the best trial functions. The usual procedure is to set K - Sag3(0^~~~~~~ ~(35) and IK, _ (36) jj j In this way, the aj and bj will be adjusted so that k is V V stationary; i. e., either a maximum or a minimum. Performing the differentiations indicated above yields the result that k will be stationary if Equations (10) and (29) are satisfied, Of course, this is not much help because these are the same equations we have tried to avoid having to solve, for reasons given previously. Nevertheless, it is informative to note that the coefficients of the ffelmholtz eigenfunctions adjust themselves so that k is stationary. This contrasts with the usual expansion of a function in a Fourier series. In the latter case,s the coefficients adjust ves so that series best approximates the function in a least squares sense0

-19The usual methods of calculating the best parameters for use in the trial functions lead to Equations (10) and (29). As pointed out above, this leads to an interesting conclusion but it sheds little light on the numerical difficulties of obtaining the parameters. Another method of obtaining the parameters will now be discussed. Consider a set of n linear homogeneous algebraic equations having n unknownss adlx-lt'LWa X + o + Ct,,.X, - 0 (37)'LX,.XI += + ct - 0,~ Suppose that this set of equations has one and only one solution and that xl # 0. These conditions can be fulfilled only if det A = 0, and if any matrix obtained by deleting the first column and any row is non-singular, Accordingly, the set of equations aCk XL + a 32 X3 a p OL s c, ak (38) 6a32 + X' which is obtained by deleting the first equation, settting x = 1 (the normalization is arbitrary), and moving the first column to the right hand side, must have a solution.

-20In an engineering problem, the accuracy with wich the coefficients of the A matrix is known may be limited to very few significant figures, say 3, Furthermore, if the manipulations are carried out with, say, 8 place accuracy, it is in general not possible to have det A = 0, exactly. Nevertheless one formally writes the set of Equations (38) and solves for x2, x3,... xno The reliability of the answers, however, may be questionable, This depends on whether or not the matrix of the coefficients in Equation (38) is ill-conditioned (12). If it is ill-conditioned, the answers may be of little value. If, on the other hand, the matrix is not ill-conditioned, the answers are as reliable as the coefficients in the matrix A. The problem of determining whether a matrix is ill-conditioned is a difficult one. The easiest method for determining this is to solve the problem formally and to observe whether or not the answers with one's intuition, experimental results, or a completely different method of solving the physical problem which led to Equation (37), if these are feasible, Now assume that the coefficients of the matrix A contain a parameter which is to be adjusted so that det A = 0, in order that a solution to Equation (37) existo If this parameter is not known, a nontrivial solution for the xi's cannot be obtained. Nevertheless, if the parameter is known "approximately," a formal solution for Equation (38), may yield a solution for the xi's which are "approximately" correct. That is, the solution obtained in this manner is approximately the same as the solution of Equations (37) with the correct parameter. If the matrix of coefficients on the left side of Equation (38) is not ill-conditioned, this would be expected as a consequence of the previous discussion,

-21Now, once again consider Equation (16) N M E E a5[E - J - o (16) -1 1= l 1 First, note that al cannot be zero. This can easily be shown. Since M (X) =X^ A((39) the orthonormality of the ~r (x) functions yields a, = J'f, (40) in particular, a f, (x) i, () dx (41) Since 0l(x) and fl(x) are everywhere > O, it follows that al > Oo Also, there can only be one independent solution for the a3 for otherwise the 0i(x) would be multi-valued, a physical impossibilityo Of course, we are concerned only with 0i(x) > 0. There are other solutions to Equation (16), corresponding to different values of k, and to 0i(x) which are negative over regions of the reactor, but these are of no concern in this work. Thus the set of Equations (16) satisfies the requirements put on Equation (37). Accordingly, take an "intelligent" guess, based on physical knowledge, for the effective multiplication factor; i. e the value of k which will make det Q = O. Furthermore, set al = 1, and solve for 1,,. a,., 2 by discarding the first row of the matrix for ~a2, a3, a.. o1

-22Q defined by Equation (19) and moving the first column to the right hand side of Equation (18). These values are accepted as approximate values - still assuming that the matrix obtained by deleting the first row and column of Q is not ill-conditioned. By using the same guess for k, setting b = 1, and proceeding 1 as above, Equation (29) is solved approximately, for b, b, oo, bM 2' 3' M9 b2 b Now Equation (34) is used to calculate the next guess fx k, presumably more nearly accurate than the initial guess. Once again the al and bi are obtained; this time the approximation is closer to their correct values, and k is recalculated by Equation (34). This procedure is repeated until two successive values of k differ by less than a specified small c. At this time, the process has convergedo The final values of the ai and bi are the correct solutions and k is the multiplication factor of the reactoro The fluxes and adjoints may then be calculated from Equations (10) and (28), A check on the validity of the final answers is provided as N M follows. The quantity Z Z aj Qli represents the product of the first row of Q and A, and must equal zero if the a, are correct. Because the first row of Q does not enter the reduced simultaneous equations solution, this procedure provides a simple and valid checko A similar check may be i N M j jl performed for the b by evaluating l Zi b 1l After any iteration, even if the problem has not converged, the two quantities calculated above must be equal. This simple observation provides a useful check after each iteration, and is particularly informative when one wonders about round-off

-23errors. That these two quantities are indeed equal is easily shown, Consider a partitioned matrix of the form CL c is a square, non-singular matrix; a and b are row and column matrices, respectively, The two checks provided above correspond to evaluating ax T and y b, where x is the solution of ex = -b and y is the solution of T T T c y = Pa, or y c = -a Hence it is clear that ax = -aclb = yTb The superscript T denotes the transpose. This is to be differentiated from the * which has been used previously to denote the adjointo It frequently happens, however, that the adjoint and transpose of a matrix are equal The reader should have no difficulty with this aspect of the notation. The iterative procedure described above for the solution of homogeneous linear equations is actually an extension of the method of W, Kohn (13). Kohn's method is more limited in that it is restricted to equations of the form Ax m Xx where A is a symmetric matrix. In this work the method has been applied to Eix X Sx where neither E nor S are symmetric. It is pointed out in Reference (14) that Kohnts method is

-24basically the same as Wielandt's method of "broken iteration." The latter, however, requires a good initial guess for the eigenvalue but it does not work if the guess corresponds to the exact value. The method given in this work is not subject to this limitation. For the purpose of carrying out numerical calculations, it is desirable to have available formulas for the evaluation of the matrix elements Eij and Sij. [Iv [IV By inserting Equation (8) into Equations (14) and (15), and performing some trivial trigonometric manipulations, one obtains A V rV V L [+i[Di(X)(i )( \! L+(Xa J - i (X) COS (-V)TD D(K)(V (42) *( — ) (L ) + t^ x)]ixM e^(x ^os (AA+V- Tr)7 dxx and L 5t, = W Jv)x) [cAr-v) + 0 Cos( x-1 dx (43) A L ) fL L L o In the code which accompanies this work, the necessary integrations are carried out by Simpson's rule, unless a cross section is constant within a subregion of the reactor. In the latter case, the integrals are very easily carried out without recourse to Simpsonss rule.

-25The method of harmonics is also applicable to multi-dimensional problems. In particular the finite cylinder, with arbitrary space dependent cross-sections in the r and z directions, will now be considered in some detail Of course the method is applicable to other geometries, but the finite cylinder serves as a good example of the two dimensional problem. The appropriate multigroup diffusion equations are r B r W (r, IJ Bj {ro ) e; Da t<* e);A ( _ i ( r/:.):~,, (,,;, ~.) (44) i-1 u +, ) (r ) (r, ) = ~ As in the one dimensional problem, multiply by a function G(x,z), integrate over the volume of the reactor, and apply Green's theorem to get R H 0o 6- I~ (45) + E~ $yc () r, (r Z) 4r ( r) KT: j',=

-26It is assumed that the problem has axial symmetry. The radius of the cylinder is R and its height is 2H. These distances, as before, include extrapolated lengths. As before, expand the fluxes in a series of eigenfunctions of the Helmholtz equation which vanish at the outer boundaries. The appropriate expansion is M^^ f^^zz. Z rr^ it7J. (46) V:I\ 1'1=r( where an are the zeros of the Bessel function Jo(x). Note that Mz and Mr need not necessarily be equal; that is, the cut-off approximation need not be the same in the axial and radial directions. Because of the piecewise continuity of the derivatives of the fluxes, they may be represented by and J M, aaVcsl O -r J - ^/~ /r: d r~~ ~ (48) R ^(r^^ <T7(L4)i a, r (49) /< jt1 I) #] (V-I5 63 7~~ D~(gt i " /442~ ~

-27Inserting Equations (46) through (49) into Equation (45) yields NMe Mr A I XA, vmK w v - M 5,v' (<o),,j 9 Me j~ "-IMJ~ /;f1 (50) where mLr^ = - -4',1^'~ v-L) TJ-[f BCt0 0 o i ~ d~~ c, 3,)Rf(, -, )dr + (- - )(_ ( 0 R o o R R and H 0 R (52) 0 rv+g(za)J *J)J

-28Equation (50) can be solved only if det Q vanishes; Q, as before, is defined by E-S/k. There are several ways in which one may conveniently arrange the matrix elements in the matrix Q. One such arrangement can be formed by placing M matrices of the form given in Equation (19) in a square array. The degree of the resulting matrix is obviously NMzMr. As before, one runs into the same problem of having to adjust k so that det Q vanishes. However, the variational principle which was used for the one-dimensional problem applies equally well to the two-dimensional case. Thus, in principle, one may obtain the fluxes, adjoints, and k for multi-dimensional problems by using the harmonics method. As in finite difference methods of solution (17), however, one has to deal with larger matrices than for the one-dimensional caseo The two-dimensional problem has not been programmed for a computing machine. Therefore it is not possible, at this time, to analyze the merits of the harmonics approach to multi-dimensional problems.

CHAPTER II KINETICS In this section, the harmonics method will be used in the analysis of some aspects of the kinetic behavior of reactors. The time-dependent multigroup equations, with no feedback, may be written as IL b J )L D (x) co 1.i ( J Ac < -9()</)+-/vf y^,(rJ (53)? i is the fraction of delayed neutrons from the l-th precursor which are emitted into the i-th energy group. These are normalized so that J_^ — l - r /={/y..-I ziwhere I is the number of delayed neutron precursors. The delayed neutron fractions Pi are defined so that -29-29

-30It is assumed that for t < 0, the reactor was in the steady state (k=l) with fluxes Oi(x,O). At t=O, some material is inserted (or removed), resulting in the coefficients of Equation (42). For t > 0, it is assumed that the coefficients are independent of time (no feedback). Furthermore, there is no external source. First, take Laplace transforms in Equation (42). Thus let (%a)= p (xe, It) (54) Equation (42) becomes r H (X) d (X, w)- < ^<) ) (X/ )- jX ((X w) If —- iw i {+ / - fX/) j dI' 11ri1 N (55) Impose the following requirements on T(x,t) (a) ji(-L,t) = Oi(L,t) = 0 for all t. (b) Oi(x,t) = Oi(-x,t) for all x and t. (c) Oi(x,t) is bounded and has continuous second-order space derivatives inside each subregion for t < o. (d) 0i(x,t) and Di(x)dii(t ) are continuous everywhere. dx

-31Also, let G(x) be a function, independent of time, which also satisfies the above conditions. If Equation (55) is multiplied by G(x) and integrated over the reactor volume there results, after applying Greents theorem and conditions (a) through (b) - D in (x)d i (x, (xdX - d(x)( (X) X, W) cX 1~~d.,~~ / (X) ( )''. _ _ I - -(X).L = )(); (s')X) (58) equations, as was done for the static case. Accordingly, and The fluxes, oi(x,O), in the initial reactor which existed at i t < 0, and the corresponding values of a are obtained by the methods given previously. These are, therefore, known quantities in the present discussion.

-32If G(x) is chosen to be 4 (x), and Equations (57) and (58) are inserted into Equation (45), there results the set of algebraic equations V - Z^)[- [ IB; X(,-^ ^J^ /-;,/ / -, d (NI (59) where E1i is defined by Equation (14), r-v.v'i~r~i3Y = Jig (X) V (Xw (X)(60) and Shiv - id J /v; Ix)its,+2 ox)61) In principle, the set of Equations (59) may be solved for each cio() by the use of Cramer's Rule. That is, each ci() is equal to the ratio of two determinants, each having the form of Q in Equation (l9)" The denominator is d -t (- E - + - while the numerator is the same, except that the [(i-l)M + v]-th column is replaced by the right hand side of Equation (59).

-33It is possible to write C (w)=M = r (62) where the numerators for each c1 as well as the common denominator have been factored. Ki is merely a constant, which generally differs for each p and i. It is obvious from the form of Equation (59) that the expansions of the determinants involved consist only of polynomials in W. Therefore Equation (62) is a valid representation of the solutions for c(W). The denominator in Equation (62) indicates that there are M(N+I) roots. This can be proved quite easily. The cOs which appear with the source term —due to delayed neutrons —contribute MI roots. This can be shown in precisely the same way as was done previously for the k eigen values. The (i-) elements form a non-singular matrix. Therefore they vi p-v contribute MN roots. Thus, as indicated above, the total number of roots is M(N+I). Physically, this indicates that each mode contributes N+I roots, and that the finite neutron velocities, or neutron lifetimes, contribute roots in exactly the same manner as the finite precursor delay times do. From the manner in which the numerators of Equation (62) are formed, it is obvious that they will contain one root less than the denominator has. If there are no multiple roots in the denominator, the Laplace transform ci(L) may be inverted to yield (Reference 15) Ma ( N /4 1 (65)

-34where ~p i [(^ 1)( -^' -( M, += (64) From Equation (46), the complete solution for the fluxes is -( v' e (65) or, by changing the order of summation, tvM(N-L) M Lt xtze t -A v 2(x) (66) From physical considerations, it is expected that the j j s are not only distinct, but real. This statement is not easy to prove, but is has been shown to be true for a variety of simple cases.(3) The assumptions made in arriving at Equation (65) are (1) The fluxes are well represented by the first M terms of a Fourier Series expansion. (2) No feedback. (3) No external source. It has not been assumed that the space-time problem, Equation (51), is separable. Thus, subject to the above restrictions, Equation (66) is the complete solution to Equation (51). Now recall that the c)j's were the roots of deS (3ZJ- i t}- ~}+SI-R)4tL t F _ ~(67)

-35Equation (67) may be regarded as a generalized Inhour Equation. It is more general than the conventional Inhour Equation because it does not assume that the space and time dependence of the fluxes in Equation (53) are separable. Thus, the manner in which the flux shape changes with time may be obtained from Equation (67). The conventional Inhour Equation has (I+1) roots, while Equation (67) has M(I+N) roots. The problem of finding all the roots of Equation (67) is thus considerably more difficult than finding the roots of the conventional Inhour Equation. If one waits long enough, the fluxes will reach their asymptotic behavior, (X t) Zk) IV v (68) where 13, the algebraically largest root of Equation (67), is the reciprocal of the stable period..d(x), defined by 1 v-I is often referred to as the dynamic flux. Thus far, a formal procedure has been given for the calculation of the time-dependent fluxes, subject to the restrictions given previously. In principle, one may obtain the solutions to any desired accuracy by choosing a large enough value of M. In practice, however, one must find the roots of Equation (67). This is the only serious disadvantage of this method because there is no fast and accurate method of obtaining the roots.

-36In order to circumvent the difficulty mentioned above, a more conventional approach to the time-dependent problem will now be taken. If, at t=O, the magnitude of the change of the coefficients of Equation (53) is small (small reactivity change), it is safe to assume that the space-time problem is separable; that is, the space dependence of the flux will not appreciably change with time. Accordingly, let Lt i(n, )==.( ) e (70) Upon insertion of Equation (70) into Equation (553) there results Xx) _ DC ( x) - < W ( x 4 ) t, ix) 1 + 0(~-f51t: G ~ gi(x),1 (X) (71) In matrix operator form, Equation (71) may be written as (L)t - V (72) where f is a column vector having components ). L G )- ( it A -/ - is a diagonal matrix whose components are -, v Vi

-37L and U are the matrices defined below Equation (25). The column vector Od of Equation (61) is not equal to the 0 of Equation (24) —except in the steady state (=O0, k=l). Od(X) is frequently referred to as the dynamic flux [see Equation (69)]; the ~ of Equation (24) is known as the static flux. Equation (72) may be multiplied by some arbitrary weighting function G(x) and the resulting equation integrated to yield x)C[Lt Jx d6 - U/g rX (73) At this point, the dynamic fluxes will be approximated by the static fluxes. It should be emphasized that this approximation is not rigorous, but in most applications it results in no difficulty. For a thorough discussion of the difference between the static and dynamic fluxes, see Reference 16. Accordingly, Equation (73) becomes /^"r(L+i= +fL+ (74) where G(x) has been chosen to be the static adjoint flux matrix defined by Equation (26). Previously, when the static case was analyzed, no mention was made of delayed neutrons. This omission introduces only a neglible error in the calculation of the static fluxes, adjoints, and keffp which were the only quantities of interest in section 1. If the time behavior is to be calculated, the presence of delayed neutrons cannot be ignored. Some results of section 1 will presently be used, and these should now be modified so as to include delayed neutronso This can be

-38done simply by replacing Xi by Accordingly, Equation (32) becomes K=;f(x)L (x) (x x(75) where X / represents a column matrix whose N components are + Equations (74) and (75) may be rearranged to yield the conventional Inhour equation, _ _ ~.g t ^ 41-(rH (76) -* where the prompt neutron lifetime O is Pr _ t [{i'^ ^) /+>;] }^^ (77) the effective delayed neutron fractions are iA w ~{^ -^ = / (78)

-39the reactivity p is defined by, K-0I (79) and the period T is simply T- -- (80) If, as in Chapter I, is approximated by X, Equations (68) and (69) become lo* ~P+1 O t - |/i y d" (81) and'(,f) t~J#+ V'= u (> (82) If each component of the flux and adjoint matrices is expanded according to Equations (10) and (28), respectively, Equation (81) becomes N M M( ~t -P N N M AM ^'= 1 I /4s4- VI (83) where ( V-<(7t (X) "v (X)J2c (84)

-40and Sij is defined by Equation (15). Similarly Equation (82) becomes ~I-v N N M M t=1 /-t $A= V-I (85) where 14 = i,t () x) (86) The results derived above enable one to compute effective delayed neutron fractions, by Equation (85), and the prompt neutron lifetime, by Equation (83), of a reactor very simply. After having solved the static problem, all that is required to compute the above kinetic parameters are a few additions and multiplications. These have been incorporated into the computer program which accompanies this thesis. Together with k, which is calculated by the methods given in Chapter I, the calculated parameters may be inserted into the conventional Inhour equation, Equation (76), to obtain the time behavior of the reactor.

CHAPTER III FURTHER PROPERTIES OF THE HARMONICS METHOD Many interesting conclusions may be obtained by solving the multigroup diffusion equations by the method described earlier. These will be discussed in this section. 1. The eigenvalue k as an upper bound. Although the set of multigroup equations is not self-adjoint, the set behaves in many ways as if it were self-adjoint. For instance, there is always a real eigenvalue k. Furthermore^ the value obtained for keff by the harmonics method is always an upper bound. This observation has not been proved mathematically, but it has not been violated in a very large number of trials. As the number of harmonics M is increased, keff decreases monotonically toward its correct value. This is because an increase of M corresponds to a better trial function for the fluxes and adjoints. This observation can be quite useful. For instance if for M = 1, a value for keff which is less than unity is obtained, one knows immediately that the reactor is sub-scritical. This is so because it is known that the correct value of keff is less than the one obtained for M = 1. Thus it is often possible to determine if a reactor is sub-critical by means of a trial function which is quite poor. Such calculations are often so trivial that they can easily be carried out by hand in a few minutes. The dependence of the prompt neutron lifetime, by contrast, can be either above or below its correct value. This behavior of keff and e* is shown in Table II for illustrative reactor 1, whose properties are given in Table I. -41

-42TABLE I PARAMETERS FOR ILLUSTRATIVE REACTOR 1 Core Reflector Width (cm) 20 40 Group 1,a 0.02 0.04 D 1.15 1.10 Vcf 0.0 0.0'1- 2 0.02 0.04 X 1.0 1.0 v x 106 52.0 52.0 Group 2 aa -0.0005 x + 0.045 0.02 D 0.20 0o15 vaf 0.05 0.0 X 0.0 0.0 v x 106 0.22 0.22

-43TABLE II keff AND PROMPT NEUTRON LIFETIME FOR ILLUSTRATIVE REACTOR 1, FOR VARIOUS VALUES OF M M keff x 105 eff P 1 1.02321 3.8347 1.00oo543 3.8150 ~5 1.00327 3.7980 7 1.00157 3.7970 9 1.00141 3.7972 11 1.00121 537970 13 1.00110 3.7971 15 1.00109 3.7971 20 1.00101 357971

-442. Convergence. There are actually two calculational procedures which must converge before the problem can be considered solved. The first procedure is the convergence, for a given M, of the iterative solution of the simultaneous homogeneous algebraic equations. The second convergence question arises when one wishes to know how large M must be; that is, how many Fourier coefficients must be retained in order to arrive at satisfactory solutions for keff, the fluxes, and adjoints. These convergence problems will now be discussed in some detail. The question which arises in connection with the first convergence problem is: "How good should the initial guess for keff be in order to assure good convergence?" It has been found, by actual computation, that, no great care is required for the initial guess. Naturally, the better the guess, the faster the convergence. However, even for very poor guesses, the number of iterations required is surprisingly small. Table III shows, for illustrative reactor 2 whose properties are given in Table IV, the number of iterations required for convergence for different initial guesses of keffo In this case, M has been chosen to be eleven. It can be seen that, although in some cases the initial guess is very poor, the necessary number of iterations is never very large. Figures 1 and 2 illustrate the rate of convergence of the fluxes for a poor initial guess of keff. It may be noted that the fluxes approach their converged values quite rapidly. Usually, in reactor criticality calculations, one has a fair idea of the value of keff before attacking the problem. Using that estimate should never lead to complications concerning convergence. When referring to Figure 2 it should be noted that the term "thermal flux"

1.2 X FIRST ITERATE FOR K =0.2 u I~ 0 FIRST ITERATE FOR K =3.0 - SECOND ITERATE FOR K =0.2 AND K = 3.0, CONVERGED FLUX 1.0 -- 1.9 X 0.8 — X —-- -J Q: \ _ _ _ ____ _ _ X-). w x 0.6 X X o~~~~~~~~~~ ~ \>11 Y w X X X L N X 0 0.0 0O.2 -0 —- -- -- -- -- -- -- -- -- -- -- - - - ^, — z 0.3 X_ — _ _ ____ _ ____ 0.2 - _ _ _ _ _ _ X 0.1 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 DISTANCE FROM CENTERLINE Figure 1. Iterated Thermal Fluxes for Initial Guesses of k w 0,2, 3,0 For Illustrative Reactor 2 in Eleven Mode Approximation,

1)~~~.1I~ ~ 0 X 0.9 I__ _ K X_ X x X X X Q8 Q7 x x X 0 0.44 0. i__ —___________ ___ X FIRST ITERATE FOR K=0.2 02 0 FIRST ITERATE FOR K =3.0 - SECOND ITERATE FOR K=02 AND K=3.0 I CONVERGED FLUX E Ol X~~,C~' -0.1 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 DISTANCE FROM CENTERLINE Figure 2. Iterated Fast Fluxes for Initial Guesses of k - 0.2, 3.0 for Illustrative Reactor 2 in Eleven Mode Approximation.

-47TABLE III ITERATES FOR ILLUSTRATIVE REACTOR 2 IN ELEVEN MODE APPROXIMATION Initial guess First iterate Second iterate Third iterate of keff 0.2 0.95109 1.00150 1.00151 1.2 1.00147 1.00151 3.0 1.00137 1.00151 4.0 1.00138 1.00151 5.0 o001o4o 1.00151 TABLE IV PARAMETERS FOR ILLUSTRATIVE REACTOR 2 Region 1 Region 2 Region 3 Region 4 Width (cm) 10 2 10 8 Group 1 a o0.02 0.04 0.02 0.04 D 1.15 1.10 1.15 1.10 Vaf 0.0 0.0 0.0 0.0 aL2 0,02 0.04 002 0o04 X 1.0 lo0 1.0 1.0 Group 2 ra 0.04 0.01 0.4 0.00.04 002 D 0.20 0.15 0.20 0.15 va.0 0. 0.05 0.0 00 0 X 0.0 0o0 0.0 0,0

-48is used somewhat loosely in this work. For our purposes the thermal flux is used synonimously with the flux in the lowest energy groupo Physically, of course, the lowest energy group may include neutrons whose energies are higher than thermal. Calculations show that the number of iterations required is independent of the number of energy groups. The second convergence problem is more difficult to deal with. However, it is possible to make some general statements concerning it. The number of energy groups used does not affect the number of modes to be used. An appropriate choice for M depends primarily on the spatial variation of the group flux which shows the greatest deviation from a pure cosine, usually the thermal flux. The presence of local disturbances of the flux, due to water gap peaking for instance. necessitates a larger value of M than would otherwise be needed. In such a case one may estimate an appropriate choice of M quite easily. The value of M chosen should be such that cos(M - 2) L has a half-period approximately equal to the range of the local disturbance, or MjL/d, where d is the range of the disturbance. It is apparent that if the total dimension of the reactor far exceeds the dimensions of local disturbances, the value of M will have to be large. Other considerations enter the choice of M. One can easily see, from the form of Equation (42), that as M is increased there comes a time when the Ei elements, for high > and v will completely dominate the matrix. p-v 1 1 it2f (ltx This is due to the presence of the term (v-_) (p-) (( )2 Di(x) sin ( -2)12 2 L x2L sin (v- 1). This will obviously occur sooner if L is small. Furthermore 2of the diagonal sub-matrices will the diagonal elements Eu for high [, of the diagonal sub-matrices will

-49dominate the sub-matrices because the above integrals are small for i / v. If Di(x) is constant over the entire reactor, all off-diagonal elements Eii ziV t ~ v, are zero. Thus the Eii elements are large compared to other elements, and the ai as obtained from the iterative cal ti and the a, as obtained from the iterative calculation, may be obtained directly by dividing the right hand side of the equations by E1i. Hence for large p, the ai are determined independently and they are all small. At this stage it may be concluded that convergence has been achieved. If Di(x) is not constant, off-diagonal terms will generally not vanish, but they will also eventually become small compared to the diagonal elements. In this case, however, it will be necessary to go to a higher value of M. In general, therefore, the harmonics method works best for those reactors which are small, which have few local disturbances, and which have constant diffusion coefficients. These conditions are of course not necessary, but they insure good results with low values of M. The disadvantage of having to use a high value of M is merely one of computation time; a larger value of M means that the degree of the matrices involved becomes larger. There have been no roundoff errors in any calculation carried out thus far, so it appears that the method is generally applicable. Note that if the fluxes were known exactly, the first M coefficients of a Fourier expansion of these fluxes would not be the same as those obtained by an M mode approximation by use of the harmonics method. This is so because the two methods converge differently; the Fourier series converges in a least-square sense while the harmonics method makes keff an extremum. This has been mentioned earlier. Of course, as M becomes large the two sets of coefficients approach each other. In fact, even for small M, the difference is not usually significant.

TABLE V HARMONIC COEFFICIENTS OF FAST FLUX FOR ILLUSTRATIVE REACTOR 1, FOR DIFFERENT VALUES OF M M. M 1 2 3 4 5 6 7 8 9 10 1 a2.04672.06868.06844.0o6818.06974.07086.07111.07111.07143 a3 -.05710 -.05876 -.05844 -.05848 -.05864 -.05864 -.05866 -.05870 1 a4.01970.01942.01880.01865.01865.01858.01853 1 a5 -.00072 -.000oo86.00020.00020.00026.00031 a6 -.00589 -.oo616 -.00616.00616 -.00619 6 a7.00412.00412.00406.00405 a8 -.00044 -.00035 -.00030 al -.00165 -.00171 10.00145

TABLE VI HARMONIC COEFFICIENTS OF THERMAL FLUX FOR ILLUSTRATIVE REACTOR 1, FOR DIFFERENT VALUESOF M M a 1 2 54 5 6 7 8 9 10 a2.5260.5570.5527.5316.5517.5517.5516.5516.5516.5516 1 a2 -.1121 -.1010 -.og09605 -.09625 -.09551 -.09483 -.09486 -.09472 -.09455 a2.05285.05206.05190.05145.05151.05150.05118.05115 a2. -.02590 -.02996 -.02441 -.02451 -.02429 -.02415 -.02412 a2 -.00oo67 -.00824 -,.00857 -.008562 -.008666 -.008744 5~~~~~~~~~~~~~~~~~~~~~~ a2.01537.01629.01627.01620.01629 ~~~~~~~a2 ~~~-.oio0 6 -.01002 -.009761 -.009745 a2 -.oo0001956 -.0006751 -.000859 88 2 005001.005396 a2 -.o004154

-52Further indication of the convergence of the fluxes for increasing M is provided by Tables V and VI, in which the harmonic coefficients of the fluxes for illustrative reactor 1 are given. The fluxes are normalized so that the first coefficient of the fast flux is 1.0. Note that as M is increased, the dominant coefficients stabilize to their final values. The data given is for values of M from 1 to 10. M = 10 is sufficient for convergence for illustrative reactor 1. 3. Comparison of harmonics method with finite difference schemes. The details of the finite difference method of solving the multigroup diffusion equations will not be repeated here. They have been widely published; see for instance Reference (17). A comparison between the harmonics method and the finite difference method seems in order. In particular, it is of interest to know the relationship between the number of mesh points (for the finite difference approach) and the number of harmonics (for the harmonics approach) which have to be used to obtain a solution with the same accuracy. An advantage of the finite difference approach is that one can choose a different mesh spacing for each region. Thus a small region which significantly affects the flux shape can be covered by many closely spaced mesh points, while other regions are covered with more widely separated mesh points. In the use of the harmonics method it is not possible to require greater accuracy in one region than in another; one must increase the number of harmonics throughout the entire reactor. Another advantage of the finite difference approach is that, for one-dimensional problems, the matrices involved are of the tri-diagonal type. These are easier to work with than the

-53completely filled matrices which one deals with in using the harmonics method(l7). One advantage of the harmonics approach is that the calculated value for keff, for a given M, is always an upper bound. The usefulness of this fact has been mentioned earlier. Another advantage is that an approximate Fourier series representing the fluxes is automatically provided. Such a series is frequently desired in reactor analysis. For instance, Gyorey (18) utilizes such an expansion in his study of the xenon oscillation problem, although he uses a sine series as well as a cosine series. Still another advantage of the method used in this work is that the flux and adjoint flux problems are solved simultaneously, frequently in less time than is required to solve either problem by finite differences. In the finite difference approach one begins by guessing a source distribution. In order to minimize the number of iterations required, it is desirable to make an intelligent guess for the source. This frequently requires some insight concerning the problem at hand. The method described in this paper begins with a guess for keff~ This quantity is usually known to a fair degree of accuracy before the solution is even attempted. Furthermore, it has been demonstrated that even for a poor initial guess the number of iterations is small. This is a distinct advantage of the harmonics approach as described in this paper. Figures 3 and 4 illustrate the fluxes in the 7, 11, and 20 mode approximations for illustrative reactor 2. Results obtained from the WANDA code (25), by using a large number of mesh points to assure a correct solution agree exactly with the curves drawn for the 20 mode approximation. Also note that a seven mode approximation yields a better representation of the

_ _GI ON _ I I I RGREGION I REGION 3 REGION 4 1.0 -I -- I —--- 0.7 x 0.6. f f-II MODES-_. 0.7 2E y~l~y 7 MODES w 0.6 NI 0Q4 0.3 E O. Z I I. 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 DISTANCE FROM CENTERLINE Figure 3. Thermal Flux in Illustrative Reactor 2 for 7, 11 and 20 Mode Approximations.

1.2 --- I.I 1.0 0.9 0.8 x 1iIAND 20 0.7 MODES LL., 0.6 w 0.5 iE 0.4 0 z REGION I z REGION 3 REGION 4 Q3 -o_3 WI 0.2 0.1 - 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 DISTANCE FROM CENTERLINE Figure 4. Fast Flux in Illustrative Reactor 2 for 7, 1, and 20 Mode Approximations.

-56TABLE VII keff FOR ILLUSTRATIVE REACTOR 2 FOR DIFFERENT NUMBERS OF HARMONICS AND MESH SPACINGS Modes keff 7 1. 0054 9 1.0027 11 1.0015 13 1.0013 15 1.0005 19 0.9999 20 0.9997 Mesh spacing* keff 2, 2, 2, 2 1.0072 4, 2, 4, 4 1.0029 10, 8, 10, 8 0.9999 20, b, 20, 16 0.9995 * 4, 2, 4, 4 implies that regions 1, 3, and 4 are subdivided into 4 mesh spacings, region 2 into 2 mesh spacings.

-57fast flux than of the thermal flux. Thus the relatively poorer representation of the latter does not affect that of the former, although the two are not obtained independently. Further insight concerning the comparison of the harmonics approach with the finite difference method may be gained by examining Table VII. Values of keff, for different mesh spacings and different M, are given for illustrative reactor 2. 4. Applications. It should be emphasized that the applications of the harmonics method to reactor calculations are the same as those of the finite difference methods; the harmonics method merely represents a different approach to the same problems which have been solved successfully by the finite difference methods for a number of years. It is of course possible to solve, by means of finite difference methods, those reactor problems which have arbitrary space-dependent coefficients in the diffusion equations. However, the finite difference codes which exist today allow only constant coefficients within a region of the reactor. The harmonics code, written in connection with this work, does allow arbitrary space dependent coefficients. Accordingly the applications of the harmonics method which are described below are based on this additional feature. The first application is to the well known Goertzel minimum critical mass problem (19). The details of the Goertzel method of analysis will not be repeated here, but results obtained from it will be checked by the harmonics method. The code written in connection with this paper is particularly well

-58suited for such problems because it allows the diffusion equation coefficients to be arbitrary functions of position, and it accepts these in analytic form. The illustrative problem is one which has been solved by Wilkins (20) in connection with his study of the minimum total mass problem, which is closely related to the Goertzel problem mentioned above. Consider a slab reactor, consisting only of water and pure U235 whose properties are: thickness of core = 14.072 cm. thickness of reactor, including reflector = 18.832 cm. Group Core Reflector 1 D =1.00 D = 1.00 aa= 0.018 a = 0.03185 vaf =0.0 vaf = 0.0'al 2= 0.03185 al- 2 == 003185 2 D = 0.1588 D 0.1588 a =? ca = 0.01955 vaf =? vaf =0.0 Wilkins shows that for the material properties given above, the critical reactor will have minimum total mass if the fuel is distributed so that Ca2(re) = 0.03755 + 0.272 cos x a2(c~re) = 5.4

-591.2 - - 0.9 — X1 ~~_________ I' 1ORE! REFE 0.8 0.7.1-I. — -I I =r 0.6 ---- --- _ —- I - 0.5 O.3 Q2 _ __I 0.I 0 I 2 3 4 5 6 7 8 9 10 DISTANCE FROM CENTERLINE Figure 5. Thermal Flux for Goertzel-Wilkins Minimum Mass Problem.

-60and va( = 0.0376 + 0.565 cos x f2(core) 5.4 According to the Goertzel theory, such a fuel distribution should yield a perfectly flat thermal flux in the core. The thermal flux for this problem has been obtained by using the harmonics code with M = 12, and plotted in Figure 5. This value of M is more than sufficient for convergence. It can be seen that the flat thermal flux in the core is very well reproduced by the harmonics method. The reason for solving the above problem is to illustrate the usefulness of a computer program which permits the added freedom in specifying cross sections mentioned above. The method could also be used to solve the flat flux problem, by a trial procedure, for reactors which do not have the many restrictions of the Goertzel problem. For instance, one could apply the method to a nonthermal reactor. Of course, any reactor in which the material composition is distributed arbitrarily can be handled in the same manner. Another problem to which the code might profitably be applied is one in which the neutron spectrum varies appreciably in different positions of the reactor. The usual approach to such a problem is to calculate the diffusion equation coefficients by appropriately averaging the microscopic cross sections over an energy spectrum corresponding to each region of the reactor. It is thus assumed that within each region the neutron spectrum is the same everywhere. Such an assumption is frequently not justified. For instance, temperature gradients or interfaces between two regions whose material compositions are very different could very well result in neutron spectra which

-61vary considerably within a region of the reactor. In order to take into account these effects it would be preferable to use a different spectrum at each point within a region. Such a procedure would result in diffusion equation coefficients which are not constant within a region, and the results obtained would be more indicative of the true physical situation. The problem of obtaining the appropriate space-dependent neutron energy spectra for the averaging process is by no means a trivial one. In fact, to this day there exists no satisfactory method for obtaining these. A possible attack to the latter problem might be made by comparing keff and fluxes obtained by experiment with those obtained by calculations using different space-dependent spectra. In this way some light might be shed on this important space-energy problem. A related application which is well suited to our code is that of replacing a many-group calculation which uses constant coefficients in each region with a few-group calculation whose coefficients vary continuously with position. The latter is a much faster computation, but that is compensated by the fact that it is more difficult to obtain the parameters for the calculation. This objection can be removed, however, as workers in the field obtain a better "feel" for working with space-dependent cross sections. To illustrate the above, consider a reactor, in which the materials are distributed symmetrically about the origin. On each side of the origin are three regions of thicknesses 7.97 cm., 2 cm., and 10 cm., in which materials are distributed as follows:

-62Element Atomic Densities (atoms/cc x 10-24) Region 1 Region 2 Region 3 U255.00231 --.0000667 u238.00911 H1 ---.0667.0667 016 --.0333.0333 12.0564 The transverse buckling was chosen to be zero. Note that Region 1 is composed of atomic densities such that the neutron spectrum is fast, region 2 is a water gap; and region 3 has a material composition typical of a thermal reactor. Primarily because of region 1, a large number of energy groups is ordinarily required for a meaningful analysis. Accordingly, the AIM program (26), which is based on a finite difference calculation, was used to calculate the fluxes by using 18 energy groups. The cross sections used for the calculation will not be reported here, but they are based on the 18 group cross section set given in Reference (27). The 18 group cross sections were then reduced to an equivalent space-dependent 2 group cross section set by combining the first 8 groups into one larger group, and the lower 10 groups into the second group. This was done by flux weighting the 18 group constants set. For instance, the absorption cross section for group 2 of the reduced set was calculated by fr(X) m n onts (xT )s(X for many points in the reactor. The result is a abs(2) (x) which varies

-631.0 1.0 -- REDUCED GROUP HARMONIC 0.9 RESULTS X 18 GROUP AIM RESULTS 0.8 0.7 w ~- -- -- -- FLUX. -— 0.5 -j ~- 0.4 k FLUX DISTANCE FROM CENTERLINE Figure 6. Comparison of Two Group Harmonic Calculation (Using Space Dependent Coefficients) With 18 Group AIM Code Having Space Independent Cross Sections.

-64continuously with position. The other reduced group constants were obtained in a similar fashion. The reduced group constants are given in Table VIII. Also the group 1 and group 2 fluxes are shown in Figure 6, where they are 8 18 compared with their equivalents, Z Xi (x) and Z Xi (x), as obtained from i-l~ i=g the AIM calculation. The agreement is seen to be quite good. The computed values of keff are.997 and 1.07 as obtained from the 18 group calculation and the equivalent two group calculation, respectively. The two group calculation, done with the harmonics code, was done for M = 20. Having noted the good agreement of the two calculations, particularly insofar as the fluxes are concerned, much further analysis for this reactor could be done with the faster two group calculation. The two methods are almost equivalent because the flux weighting of the coefficients conserves the total absorption and neutron production rates at each point. The leakage rate is not conserved by flux weighting however. In order to do so, it would be more appropriate to define DJ, for instance, by YX / ix'dX D2 would be defined in a similar way. Flux weighting of the diffusion coefficients is obviously a far simpler procedure. Furthermore it is seen to yield fluxes which agree quite well with the 18 group results. The agreement in the values keff is not very good however. This is to be attributed to the method of reducing the cross sections; not as an indication of poor performance of the harmonics method. In addition to the calculation described above, a comparison was made between the results obtained from the 18 group calculation and a 2 group

TABLE VIII REDUCED SPATIAL DEPENDENT TWO GROUP CROSS SECTION SET Group 1 Region 1 Region 2 Region 3 D 1.50 1.05 1.05 a (.0519x2 -.103x + 9.94) 10-3.00585x -.0087.04962.o14o 0 < x < 3.33 v f.0003x +.0130 3.33 < x < 6.65 0.0004.000592x +.0111 6.65 < x < 7.97 a12, (.ol8x2 -.036x +.80) 10-3.00585x -.0087.04962 Group 2 D.643 0 < x < 5 -.02x +.455.223.163 0 < x < 5.76.ooo8 +.0277 9.96 < x < 12.97 5a |.1225x2 - 1.48x + 4.6o 5.76<x< 7.97.00254x -.0134.00042x +.0325 12.97 < x < 16.l.0393 x > 16.1.135 o < x < 5.45 Vcf.174x -.815 5.45 < x < 7.20 0 -.00089x2.0 x +.00 x +.268.724x - 4.76 7.20 < x < 7.97

-66TABLE IX REDUCED SPACE INDEPENDENT TWO GROUP CROSS SECTION SET Region 1 Region 2 Region 3 Group 1 D 1.50 1.05 1.05 aa.010.044.o4962 vaf.014 0.0.0004 l< 2.0008.045.0496 Group 2 D.643.275.223 a.163.010.039 vaf.135 0.0.050 f

-671.0 —' —. ____ ___ --- REDUCED GROUP HARMONIC RESULTS ^ x **^ I Ix 18 GROUP AIM RESULTS 0Q9 0 0.8 GROUP I )n J7 FLUX ---- w L \ N Q5 - GROUP 2 0i 0.1 0 2 4 6 8 10 12 14 16 18 20 DISTANCE FROM CENTERLINE Figure 7, Comparison of Two Group Harmonic Calculation (Using Space Independent Coefficients) with 18 Group AIM Code Having Space Independent Cross Sections.

-68calculation using constant cross sections. The constant cross sections were taken from the space dependent set by choosing values far from region boundaries. These were felt to be most representative for the material compositions of the regions. The cross sections are given in Table IX. The value of keff obtained by using the constant two group set was.888 which is considerably worse than that obtained by using the space dependent set. Figure 7 compares the fluxes obtained by the constant two group set with those obtained from the 18 group AIM calculation. The agreement is not so good as that obtained with the space dependent cross section calculation.

CHAPTER IV DETAILS OF THE COMPUTER PROGRAM A computer program, written in "MAD" (23), which is based on the theory given previously has been written for the IBM 704 with 8192 word capacity. The details of the code will be given now. 1. Input A list of required input is given below. The manner in which these are to be punched on cards will be given later. M -- number of modes desired. N -- number of energy groups. REG -- number of regions to the right of center-line. Y(1), Y(2),..., Y(REG) -- These are the positions, in cm., of the region boundaries, measured from the centerline of the reactor. Y(n), for n < REG, is the position of the interface between regions n-l and n. Y(REG) corresponds to the half-thickness of the reactor. K1 -- This is an estimate of keff. It is not necessary that it be a good guess, but it is desirable as it will sometimes reduce the computation time. EPSK - This quantity represents the convergence criterion. EPSK is the accuracy to which keff is desired, for a given M. XI(1), XI(2),..., Xl(N) -- The fraction of neutrons, resulting from fission, which are emitted into group i is given by XI(i). POINTS -- The number of points in the reactor at which the fluxes, or adjoints if desired, are to be computed and printed. These points are -69

-70automatically spaced at equal intervals of length Ax = L/POINTS. The quantity POINTS does not enter the computation. It merely determines the amount of information to be printed. The fluxes and / or adjoints at the edge of the reactor, whose values are zero, are not printed. P -- This quantity is equal to 2 if a printout of the adjoint fluxes is desired; otherwise P = 1. KIN -- This quantity is equal to 1 if a calculation and printout of prompt neutron lifetime and effective delayed neutron fractions are desired. If these are not desired, KIN = 2. Q(1),...,Q(REG) -- These are the coefficients of the diffusion equation. These will consist of several cards, one for each coefficient and each energy group, as will be described later. If a coefficient is constant within a region, its value is printed on a card. However, if the coefficient is not constant, the value assigned to is -1.0, and an analytic function describing the space dependence of the coefficient must be inserted in the form of an external function (or subroutine). This will be described in greater detail later. Cards for the following input quantities are to be punched only if KIN = 1. ND -- the number of delayed neutron groups. BETA (1),..., BETA (ND) -- These quantities are the delayed neutron fractions, which may be obtained in Reference (24). SQ(l,l),..., SQ(ND,N) -- These are the fission spectra for each delayed neutron group.

-71V(l,l),..., V(N,REG) -- The group velocities. These are allowed to vary from region to region, but must be constant within each region. 2. Card Format The data required are in Fixed point, Floating point, and Integer forms. These are designated by the letters F, E, and I, respectively. E and I data must be punched to the extreme right of the alloted columns. F data may be punched anywhere within the alloted columns, but all decimal points must be included. For further information concerning the various number forms, see Reference (23). The contents of each input card are described below. First Card: Item Columns Form Comment M 1-2 I 1 < M < 20 N 3-4 I 1 < N < 4, but 1 < MN < 40 REG 5-6 I 1 <REG< 11 P 7-8 I KIN 9-10 I POINTS 11-14 I 1 < POINTS < 160 K1 15-24 F EPSK 25-34 F BUCK 35-44 F Next Card: Y(L),...,Y(REG) 1-10, F A maximum of seven values may be placed on 1 card. If REG>7, 11-20, etc continue on another card, placing Y(8) in Cols. 1-10, etc.

-72Item Columns Form Comment Next Card: XI(1),..., 1-10, F These are assumed independent XI(N) of position. 11-20, etc. Next Card: a: i(l),... 1-12, F If N = 1, do not include these a. *(REG) cards. The ordering of the J iG 13-24, etc. cards is as follows: i = 2, j = 1; i - 3, j = 1,2; i = 4, j = 1,2,3. Cards must be printed only for i = 2 to N. A maximum of 6 values may be punched on a card. If REG>6, continue on another card with the same format. See discussion of Q(1),...,Q(REG) for added details. Next Card: Dl(l1),..., 1-12, F A maximum of 6 values may D1 (REG) be punched on a card. If 13-24, etc. REG>6, continue on another card with the same format. See discussion of Q(l),..., Q(REG) for added details. Next Card: a (1),.., 1-12, 13-24, F See previous comment. a1 a(REG) etc. vaf1(l), 1-12, 13-24, F See previous comment....,vafl (REG) etc. The following cards contain values of Di, ai, vaf for i = 2 to N in the same order and format as is given for i = 1 above.

-73Item Columns Form Comment Next Card: ND 1-2 I This card and all succeeding ones are included only if KIN = 1, 1 < ND < 6 Next Card: BETA (1),..., 1-10, F BETA (ND) 11-20, etc. Next Card: SQ(1,1),..., 1-10, 11-20, F etc. SQ(1,N) The following cards contain values of SQ(i,l),..o,SQ(i,ND) for i = 2,...,N with the same format as is given for i = 1 above. Next Card: V(1,1),., 1-12, 13-24, E A maximum of 6 values may etc. be punched on a card. If V(1,REG) REG > 6, continue on another card with the same format. The V's are constant within each region. The following cards contain values of V(il),ll1.1V(i,REG) for i = 2 to N in the same order and format as is given for i = 1 above. 3. Output The information which is printed for each problem is as follows: (1) The input data. (2) keff after each iteration. (3) The harmonic coefficients of the group fluxes and adjoints (if desired).

-74(4) The group fluxes and adjoints (if desired) for points separated by Ax. (5) Ax. (6) After each iteration, the flux and adjoint checks, as discussed in Chapter I. The Fourier coefficients, the fluxes, and adjoints are normalized 1 1 so that a b1 = 1; that is, the first Fourier coefficients of the first flux and adjoint groups are set equal to unity. 4. Further details In addition to those given previously, some added comments regarding the code are in order. It has been pointed out earlier that if a coefficient of Equation (1) is not constant within a region of the reactor, a subroutine must be written to describe its space dependence. A brief description of the subroutines follows. The scattering cross sections a. i(x) are described by a subrouJ-~ - tine named SIGRE.; Di (x) (), and vaf (x) are described by the names Do, SIGA., and NUSIFG., respectively. These last three are different entries to the same subroutine. As a simple illustration, let us assume that a2 (x) in a certain a2 region which extends from x = 10 to x = 15 may be described by aa(x) = 1.15 +.05x -.04x2 Also assume that for another region a (x) =.04 +.005e53X for 20< x < 26 al and assume that all other coefficients are constants within regions.

-75The values assigned to aa and aa in the appropriate regions are assigned the value -1.0. The following subroutine is then included with the main program: EXTERNAL FUNCTION (I,x) INTEGER I ENTRY TO D. FUNCTION RETURN 0.0 ENTRY TO SIGA. WHENEVER x.G.10.0.AND. x.L.20.0.AND, I.E.2 FUNCTION RETURN 1.15 +.05*x-.04*x*x OTHERWISE FUNCTION RETURN.04 +.005*EXP.(-3.0*x) END OF CONDITIONAL ENTRY TO NUSIGF. FUNCTION RETURN 0.0 END OF FUNCTION Obviously some knowledge of MAD is required to write the necessary subroutines. The details of the language will not be discussed here because they are available, in great detail, in Reference (23). Note that the third, fourth, eleventh, and twelfth statements must be included although D and vaf are constant for the problem at hand. Entries must be provided for all subroutines even if they are not used. When all coefficients are constant within a region, one need not be concerned about the mechanics of MAD. If the subroutines are used} however, the binary IBM card marked EX1 (which describes SIGRE.) and/or the

-76cards marked EX2 (which describe D., SIGA., and NUSIGF.) must be removed and replaced by the appropriate subroutine(s) for the problem at hand. The main program consists of three sections, all of which are independent, except for some data and intermediate results which are common to all three. If it is desired to run more than one problem at a time, one simply stacks succeeding deck of data behind the first deck. 5. Timing As with most programs which have many options, it is difficult to assign an execution time estimate for a general problem. In this program the number of energy groups, the number of harmonics, the initial guess for keff which affects the number of iterations required, as well as several other factors affect the execution time. It might be useful to the reader, however, if one example is given. The fluxes, adjoints, and keff for illustrative reactor 1 were obtained in 46 seconds for 2 energy groups, 10 harmonics, and an initial guess for keff of 1.3. 6. Flow Diagram A flow diagram of the computer program follows:

-77FLOW DIAGRAM OF COMPUTER CODE -------------- / YYES READ M, N, REG, R I ( L,X"(_)....I PRINT....l. POINTS, P, KIN, INPUT NO ElJ = Y(1),...,Y(REG) -LY 4i = 1i,...,N-1 1)iDEPENDENT COEI B( = B(U) +,...,M YES U= PRINT E CIEN SP IGRE.(IJ,) | Y SIMPSON'S GRE.(I, V 0 NO R O/BR=1,..,EG U0.,2 -i C. -ANT S() BE.(IJx) IY R'')-"(~'SI+ 1tIl)](A U = 1. M; GVE 1,,M YES BAPEKDEMT =OEP- PRINT B(U) = B(U) = EXEUT - (). BA(I)EB = [ -)B(+v-l) REG ^*N y —— S NO ss^ i+l / [J =1,...-M - -1 (^( T^~=i IS1.,MYES. v = 1,...,M Cr = 1,...,M; v = 1,..,M

-78B(U) n 0 EXECUTE B(2M(i+l)+U) n B(U) U=O,..., 2M-1 MANIP. uO... 21. (NUSIGF.O), sj o LYES Sj1 NO L= 1 |N MNO - YES'-YE PRINT x (|L POINTS 5- 7 ~ 5 bi Sij as PRINT R, V i j v EXECUTE K 1K1' CEZEiaij j VAR. K3*Kl' -1 V V v ji j NO =20 1 — << K-Kil —-l ---- NO <EPSK YES,KK EXECUTE IS - _VAR. YESKIN YES COMPUTE COMPUTE i* NORM READ & PRINT -Leff by by Eq.(8. ~ bh SR Ey. (8~)....'1-I

-79SUBROUTINE MANIP.(F.,Z) ENTER I n /S Y PRINT I( I( |~AD~ Q~(R); I 1 11(R) EAD Q(R); -— ^Rriy_^ (R) ^ NO ^ Get subroutine R1,..,REG I T \ Idescribing spac -— \B *- PP )Idependent coefficient F.(I,x) By Simpson's rule get B(Z'M+U)- B(Z'M+U) + ( F (I,x) COS (U)dx R YES B(Z M+U) = B(Z.M+U) + B(Z-M) a B(Z-M L' (R)[SIN(U bY(?)-SIN(UY(Rl)) ]' Q(R)[Y(R)-Y(R-1U Ut L L Return to Main Program

-80o SUBROUTINE VAR. Tij EiJ + ^ SOLVE FOR A, BY GAUSS ELIMINATION,s p~ k ___ 14'Y k =SYSTEM OF EQUATIONS TH = W, WHERE E3p V p=1,,...-M TR IS T WITHPUT 1st ROW AND 1st l, r 3, l,....I'; COLUMN; A = (2,.,aN) AND W IS i - 1,..,N 1st COLUMN OF T WITHOUT Ti. J 1,..,,N OBTAIN TRANSPOSEOF T; SOLVE FOR B, BY GAUSS ELIMINATIION, T AiIJ Ti SYSTEM OF EQUATIONS TB = V, V WHERE TR IS T* WITHOUT 1st ROW *,. AND 1st COLUMN; B = (b,...,bM) V 1 M AND V IS 1st COLUMN OF ~ * l WITHOUT T*ll i l j =,N CALCUIATE AND / __________________. PRI -. Is T" zs~~ J 1'IX I CALCUIATE I I I —AND PRINT ~\ l2, l PRINT; i =1,...,N YES Ip l 1,..,M Mi AND GROUP FLUXES AT =P t -> INTERVAL Ax _I PR!NT_ bT; i=l,-,..,N l l,... M,M AND GROUP AD.JOINTS / AT INTERVAL Ax RETURN TO MAIN PROGRAM

APPENDIX COMPUTER PROGRAM The computer program used for the calculations described in the text is given on the following pages. The principal symbols have been defined previously. Those which have not been defined are temporary storage locations. -81

-82* COMPILE MAD, EXECUTE, DUMP, PUNCH OBJECTPRINT OBJECT CR1 PROGRAM COMMON MN,L,PIKlEPSK,XI POINTSEPVV,MN,VS 1,BYREGKIN,PIOL DIMENSION E(1599)*S(1599) XI(4) VV(40) V(4) Y(10) DX( 10) 1 BA(20),NMESH(10) B(159) Q(10) INTEGER VtMNVVtI,JMNMUNUPtR,US3,POINTStREGNMESHKIN Y(O)=0o READ FORMAT CARD1lMNREG,PKINPOINTSKl.EPSK PRINT FORMAT TITLEN,MREG READ FORMAT CARD3, Y(1)...Y(REG) PRINT FORMAT INT, Y(O)*..Y(REG) PRINT FORMAT WANT, K1, EPSK READ FORMAT CARD2, XI(1)...XI(N) PRINT FORMAT SPECT, XI(1).**XI(N) WHENEVER P.E.2, PRINT FORMAT CHAL WHENEVER KIN.E.1, PRINT FORMAT BOT PI3. 1415927 L=Y(REG) PIOL=pI/L MNUM*N THROUGH JOHNFOR R=1,l1,RG.REG NMESH(R) (Y(R)-Y(R-1))/L*200. NMESH( R ) 2* ( NMESH ( R ) /2 ) JOHN DX(R)=(Y(R)-Y(R-1) )/NMESH(R) THROUGH NAJ, FOR I=1,1, I*G.MN NAJ VV(I)= -1+(I-1)*MN THROUGH SFA, FOR I=1,1, I.G.N SFA V(I) ( I-1 )M WHENEVER N*E.*1TRANSFER TO REC THROUGH LAN.FOR I=1,1*I.GN-1 THROUGH LANFOR J=I+1,lJ.G.N THROUGH LAN, FOR MU1,L1,MU.G.M S3=VVtMU+V( I ) +V(J) THROUGH LAN, FOR NU=1t1,NU.G.M LAN E(S3+NU)=O.0 THROUGH CHUCK.FOR I=2,1,I.G.N THROUGH CHUCK,FOR J=1 1,JG.I-1 PRINT FORMAT SGROUPJI THROUGH SMITH, FOR U=O,1,U*G.2*M-1 SMITH B(U)=ot READ FORMAT CONSTy 0(1) ***Q(REG) THROUGH WMAN, FOR R=1,1,R.G*REG WHENEVER Q(R).E.O.O, TRANSFER TO WMAN WHENEVER Q(R).G.0.0, TRANSFER TO LIN Sl=O0 S2=0. THROUGH BAR, FOR XzY(R-1)+DX(R),2.*DX(R)~XeG*Y(R)-2.*DX(R) SIS1l+SIGRE.(I JX) BAR S2=S2+SIGRE.(IJX+DX(R)) B(O)=84O)+(SIGRE.( I JtY(R-1)+DX(R)/100.)+4,*Sl+2,*S2+ 1 4.*SIGRE.(IJY(R)-DX(R))+SIGRE.(IJy(R)-DX(R)/10.0))*DX(R) 2 /3. THROUGH MAN, FOR U=1lyUGo2*'M-1 Sl=O. S2=0. THROUGH ARF, FOR X=Y(R-1)+DX(R),2.*DX(R),X.G.Y(R)-2**DX(R) SI=Sl+SIGRE*( IJX)*COS.(U*PIOL*X) ARF S2=S2+SIGRE.(IJ,X+DX(R))*COS.(U*PIOL*(X+DX(R))) MAN B(U)aB(U)+(SIGRE.( I,JY(R-1)+DX(R)/100.)*COS.(U*PIOL*Y(R-1))

-831+4.*SI+2.*S2+4**SIGRE. IJY(R)-DX(R))*COS (U*PIOL*(Y(R)-DX(R 2 )))+SIGRE.(IJY(R)-DX(R)/100o)*COS.(U*PIOL*Y(R)))*DX(R)/3. TRANSFER TO WMAN LIN PRINT FORMAT SCATRQ(R) B(O)=B(0)+Q(R)*(Y(R)-Y(R-1)) THROUGH KAD, FOR U=1,1,U.G.2*M-1 KAD B(U)"B(U)+(SIN.(U*PIOL*Y(R))-SIN (U*PIOL*Y(R-1))*L*Q(R)/(U* 1 PI) WMAN CONTINUE THROUGH CHUCK, FOR MUz1,1, MU.G.M S3=W(MU+V( I )+V(J) THROUGH CHUCK, FOR NU=,1,1 NU.G.M CHUCK E(S3+NU)x B(MU+NU-1)+B(.ABS.(NU-MU)) REC THROUGHLUVE, FOR NU-11,NU.G*M LUVE BA(NU)=(NU-.5) PIOL THROUGH JONES. FOR I=11,I.G.N PRINT FORMAT PHAL9I THROUGH YIP, FOR U-O,19U.G.4*M-1 YIP B(U)=-0 EXECUTE MANIPo(D.,ODIFF) EXECUTE MANIP.(SIGA.,2ABS) THROUGH GAH, FOR MU=1*1, MU.G.M S3=VV(MU+V( I) )+V( I) THROUGH GAH, FOR NU=1~1, NU.G.M GAH E(S3+NU)= -B(2*M+MU+NU-1)-B(2*M+.ABS.(NU-MU))-BA(MU)BA(NU)* 1 (B(.ABS.(NU-MU))-B(MU+NU-1)) THROUGH GEE, FORU=O,1,U*G.2*M-1 GEE B(U)0o EXECUTE MANIP.(NUSIGF.,OFISS) THROUGH GLEN, FOR U0,1,U.Ge.2*M-1 GLEN B(2*M*( I- 1 )+U)B(U) THROUGH JONES, FOR J-ltlJ.G.N WHENEVER XI(J)*E.O., TRANSFER TO LUKA THROUGH ABC, FOR MU=1,1, MU.G.M S3=VV(MU+V(J) )+V( I) THROUGH ABC, FOR NU=1,1, NU.G.M ABC S(S3+NU) ( B(MU+NU-1 +B( ABS (NU-MU) ) )*XI (J) TRANSFER TO JONES LUKA THROUGH DEL, FOR MU=1,1,MUoG.M S3=VV(MU+V(J))+V( I) THROUGH DEL9 FOR NU=t,1,NU.G~M DEL S(S3+NU)= 0.0 JONES CONTINUE, S 3MN*MN-1 INTERNAL FUNCTION (F., Z, STAT) ENTRY TO MANIP. INTEGER Z. STAT READ FORMAT CONST9 Q(l) *..Q(REG) THROUGH RAM, FOR R=l1,lR.G.REG WHENEVER Q(R).E..O09 TRANSFER TO RAM WHENEVER Q(R).GO.O, TRANSFER TO ABE Sl=0O S2=0. THROUGH NAH, FOR X=Y(R-1)+DX(R),2.*DX(R)tXG-Y(R)-2.*DX(R) Sl=S1+F.( IX) NAH S2=S2+F (I X+DX(R)) B(Z*M)zB(Z*M)+(F*( I,Y(R-1 )+DXtR)/100 )+4*Sl+2. *S2+4,**F*( I 1 Y(R)-DX(R))+F.(I,Y(R)-DX(R)/100.))*DX(R)/3.

-84THROUGH AKA, FORU=11,U.G.2*M-1 Sl10. THROUGH JET, FOR X=Y(R-1)+DX(R),2.*DX(R),X.G.Y(R)-2.*DX(R) SlS1+F.( IX) *COS.(U*PIOL*X) JET S2=S2+F.( I X+DX(R))*COS(U*PIOL*(X+DX(R))) AKA B(Z*M+U) B(Z*M+U)+( F (ItY(R-1)+DX(R) /100.)*COS.(U*PIOL*YIR 1 -1))+4,*S1+2.*S2+4.* F o(I,Y(R)-DX(R))*COS.(U*PIOL*(Y(R)-DX 2'(R))+ F.(IY( DXR)XR)/100.)*COS.('U*PIOL*Y(R)))*DX(R)/3. TRANSFER TO RAM ABE PRINT FORMAT STAT, R, Q(R) B(Z*M)=B(Z*M)+O(R)*(Y(R)-Y(R-l)) THROUGH GRAS, FOR U=sl,1U.Go2*M-1 GRAS B(Z*M+U)=B(Z*M+U)+(SIN.(U*PIOL*Y(R))-SIN.(U*PIOL*Y(R-1) ) )*L* 1 (R)/(U*PI) RAM CONTINUE FUNCTION RETURN END OF FUNCTION VECTOR VALUES CARD1$512,vI4,2F12.6*$ VECTOR VALUES SGROUP=$26HOSCA CROSS-SECTIONS, GROUPI3,11H INT 10 GROUPI3*S VECTOR VALUES SCAT=$7HOREGIONI3,12H SIGMA SCAT=FlO*6*S VECTOR VALUES TITLE=$8H1GROUPS=I3,S10,6HMODES=I39S10S8HREGION1S I3*$ VECTOR VALUES CARD3S$(7F10.2)*$ VECTOR VALUES CONST=$(6F12 6)*$ VECTOR VALUES CARD2=S4F10.4*$ VECTOR VALUES INT=$.16HOBOUNDARY POINTS/(9F13.2)*$ VECTOR VALUES SPECT' $9HOSPECTRUM/4F20.7*S VECTOR VALUES WANT=$11HOINITIAL K=F10.7,S10t17HOESIRED ACCURA 1CYtFlo.7*$ VECTOR VALUES CHAL=$15HOADJOINT WANTED*$ VECTOR VALUES BOT: S16HOKINETICS WANTED*$ VECTOR VALUES PHAL=$6HOGROUPI3*S VECTOR VALUES DIFF=$26HODIFFUSION COEF FOR REGIONI3,3H ISFlO. 16*$ VECTOR VALUES ABS=S33HOREMOVAL CROSS-SECTION FOR REGIONI393H 1ISF10,6*$ VECTOR VALUES FISS=$19HONUSIGMAF IN REGIONI3,3H ISF10.6*SEXECUTE SEOPGM. END OF PROGRAM * COMPILE MAD, PUNCH OBJECT9 PRINT OBJECT EX1 EXTERNAL FUNCTION (I,J,X) ENTRY TO SIGRE.!'NTEGER I, J FUNCTION RETURN 0.0 END OF FUNCTION * COMPILE MAD, PUNCH OBJECT, PRINT OBJECT EX2 EXTERNAL FUNCTION (IX) INTEGER I ENTRY TO D. FUNCTION RETURN 0.0 ENTRY TO SIGA. FUNCTION RETURN 0.0 ENTRY TO NUSIGF. FUNCTION RETURN 0.0 END OF FUNTCTION * BREAK

-85* COMPILE MADO ExECUTE* DUMP, PUNCH OBJECTPRINT OBJECT CR2 PROGRAM COMMON M.NLPI,K1tEPSKXI.POINTS.EP.VVMN.V.S I,B"-T, REG', K tTPIOL, A C DIMENSION E(1599),S(1599),XI(4),VV(40),V(4).T(1599),A(40)t 1 C(40),BA(20),CC(40) DD(40) B(159) Y( 10) W(122).BB(40) EQUIVALENCE (W-BB),(W(41 )CC),(W(82),DD) INTEGER VMNVV, I,JMN,MU.NUP,RUS3,POINTSKIN LA:MN-1 THROUGH JIZq FOR NU1=lt, NUG.eM JI7 BA(NU)=(NU-.5)*PIOL DX=L/POINTS PRINT FORMAT SPACE, DX K=0, THROUGH PAM, FOR R=1llR.G.20 *OR..ABS.(K-K1).L*EPSK KEXECE VARK1 EXECUTE VAR. S2=O.0 THROUGH DUB. FOR Il=,l1 I.G*MN THROUGH DUB, FOR J=1,1 J.G.MN T=C(tI-1)*A(J-1) MU=VV( I )+J S1Sl+T*E (MU) DUB S2=S2+T*S(MU) K 1=-S2/S1 PAM PRINT FORMAT TEMP. RK1 K=K1 EXECUTE VAR. WHENEVER KIN.E.1 EXECUTE SEQPGM. OTHERWISE EXECUTE SELPGM*(1) END OF CONDITIONAL INTERNAL FUNCTION ENTRY TO VAR. THROUGH SOFT. FOR- MU=2,1. MU.G.MN I=WVVMU)+1 S3=VV(MU-1) T'(S3+MN) -E(I)-S(I)/K THROUGH SOFT, FOR NU=1*1, NU.E.MN J=I+NU SOFT T(S3+NU)=E(J)+S(J)/K Sl1SLINEQ (TLA. l.0ABB9CCDD) A(O)1.*0 THROUGH CLAR. FOR I=ll, I.E.MN CLAR A(I)= T(VV(I)+MN) THROUGH BILL, FOR I=11,lI.G.MN THROUGH BILL, FOR J=l,1J.G.*MN BILL T(VV J)+I )E (VV( I)+J)+S(VV( I )+J)/K THROUGH HARD, FOR MU=2.1. MU.G.MN I=VV(MU)+I S3=VV(MU-1 ) T(S3+MN)= -T(I) THROUGH HARD, FOR NU=1,1, NU.E.MN J=I+NU HARD T(S3+NU)=T(J) S1'SLINEQ. tTLAIl.OC.BBCCDOD) C(O)=1.O

-86THROUGH GLAT, FOR I=1,l1I.E.MN GLAT C(I)-T(VV(I)+MN) WHENEVER P.E.1, TRANSFER TO ONLY Sl1=Oo THROUGH MAF, FOR =0O,I1 I.eEMN J=VV(I+1)+l MAF SlsS1+C(I)*(E(J)+S(J)/K) PRINT FORMAT CHECK2, SI ONLY Sl=O.o THROUGH FAM, FOR I=0,1,I.E.MN FAM Sl=Sl+A(I)*(E(l)+S(I)/K) PRINT FORMAT CHEC, S1 WHENEVER P.E.1, TRANSFER TO SLIM THROUGH MYF, FOR I=0,1, I.E.MN MYF T(I)=C(I) EXECUTE FIN.(AFOUR, ADJ) SLIM THROUGH MYA, FOR I=0,1, I.E.MN MYA T(I)=A(I) EXECUTE FINo(FOURIE, FLUX) FUNCTION RETURN END OF FUNCTION INTERNAL FUNCTION (STT9STV) ENTRY TO FIN. THROUGH LASH, FOR I=191, I.G.eN LASH PRINT FORMAT STT, I, T(V(I)) * T(V(I)+M-1) IHROUGH BAR, FOR I=1,l1,IG.N Wo( )t0 THROUGH BUT, FOR MU'=I11MU.GeM BUT W(O)=W(0)+T(V(I)+MU-1) THROUGH ABZt FOR J=l,1,J.E.POINTS Sl=O. THROUGH GEM, FOR MU=1,1MU6G.M GEM Sl=Sl+T(V(I)+MU-1)*COS (BA(MU)*J*DX) ABZ W(J)=s1 BAR PRINT FORMAT STV,I-. W(O)...W(POINTS-1) FUNCTION RETURN END OF FUNCTION VECTOR VALUES TEMP=$1SHOITERATIONI3,4H, K=E15.7*$ VECTOR VALUES CHEC=$17HOCHECK PARAMETER=E15.7*S VECTOR VALUES FOURIE=$28HOFLUX COEFFICIENTS FOR GROUPI3/ 1 (7E17.7)*$ VECTOR VALUES SPACE=$16HOSPACE INTERVAL=E15.7*$ VECTOR VALUES FLUX= S6HOGROUPI3,6HFLUXES/(7E17.7)*S VECTOR VALUES AFOUR=$31HOADJOINT COEFFICIENTS FOR GROUPI3/ 1 (7E17.7)*$ VECTOR VALUES ADJ=$6HOGROUPI3,8HADJOINTS/(7E17,7)*S VECTOR VALUES CHECK2=$11HOADJ CHECK=E157*$ END OF PROGRAM * COMPILE MAD, PUNCH OBJECT, PRINT OBJECT EX3 EXTERNAL FUNCTION (ANN,MM. IIX,JXPVT) INTEGFR IJ,N,JXIXL,TTRKN MVERT110 BOOLEAN PVT MVERT115 ERASABLE JXMXLTDIVS,R,TMKNN ENTRY TO SLINEQ. MVERT140 KN = NN + *ABS.MM MVERT145 PVT = 1B MVERT185 N * NN MVERT225 THROUGH SET, FOR I=1,1,I *G. N MVERT245

-87Mti1 —'-Kl~i -1 i'-MVERT25O SET PVT(I) OB MVERT255 THROUGH ELMOL1 OtR l" t1~ L.GIN MVERTZ95 XMX=OO MVERT300 THROUGH MAXEL, FOR 1s,19.1.Go N MVERT330 WHENEVER PVT(I), TRANSFER TO MAXEL MVERT335 THROUGH NMAXELM. FOR J=,1,J.G. N — VERT340 WHENEVER PVT(J), TRANSFER TO MAXELM MVERT345 -T-I' r I) +J MvERKT3OWHENEVER.ABS.A(T) G.. *ABS.XMX MVERT355 XMX AT)} MVERT360 IX = I MVERT365 JX = J MVERT37O MAXELM END OF CONDITIONAL MVERT375 MAXrL CONTINUE MVERT38D WHENEVER IX *E. JX, TRANSFER TO DETRM MVERT405 THROUGH SWAP, FOR J11-'Z J.G6 KN -MERT410T=I ( I)+J MVERT415 R=I (JX)+J MVERT420 TM=A(T) MVERT425 A (T) A (R) — MlVERTT-s SWAP A(R)=TM MVERT435 DETRTM PVt (JX ) =B M-VERT460IX(L)=IX MVERT465 JX ( L )=JX MVERT40 — ( I (JX)+JX) 1*0 MVERT535 THROUGH DIV FOR Jtl,.J.G. KN -MVER T=I(JX)+J MVERT545 D V A(T)=A(T)/XMX MVERT550THROUGH ELMCOL, FOR I=1,1,I *G. N MVERT570 WHENEVER I.E* JX, TRANSFER TO ELMCOL -MVERT575 — T (I )+JX MVERT580 XMX A(T) MVERT5f — A(T)on0O MVERT590 THROUGH SUB, FOR J=ll,,J *G. KN MVERT"95 — TI (I +J MVERT600 StU- A(T)=AtT)-A( IJX)+J)*XMX MVERT605ELMCOL CONTINUE MVERT610 -.FUNCTION RETURN- 1.0 VERfT69END OF FUNCTION MVERT695 T —-BREAK COMPILE 4AD, EXECUTE, DUMP, PUNCH OBJECTPRINT OBJECT CR3 P-OGRAM COMMON MNtLPI,K1EPSKXI POINTSEtPVVtMN^tVr-S 1,BYREG,KIN,PIOLA,C -DMENSION E( 15"9,)SI1599),XI(4~,VV(40)t V(4)t Y<lOttI,'tl' 1 S2(6),BETA(6),SQ(24),Q(10) A(40),C(40 INTEGER V I,JMNMUNU,R,US3,REGNDS5,S6 READ FORMAT NDEL. ND PRINT FORMAT PDEL, ND READ FORMAT BET BETA(1 ) *BETA(ND) PRTNT FORMA'T'PBET, BETA(1)...BETA(ND) THROUGH-WAW, FOR U= 1,1,U.G.ND READ FORMAT DEL, SQ((U-1)*N+1)*.SQ(U*N) PRINT FORMAT DELPU SQO (U- 1)*N+1 ).SQ(U*N) WA"W 7S'Wsl=0.0 TO - SSUM.C 1F0 * N THROUGH SU FOR THROUGH SUM* FOR I=191~IoGoN

-88PRI NT - TOIRMA T - GRP, I THROUGH JIC, FOR U=sOl U.G.2*M-1 JT- E (U~~=0o READ FORMAT CONST, Q(1)**.Q(REG) THROUGH PHIL, FOR R* 11-. R.G.REG PRINT FORMAT PQ, R. Q(R)'E'(O0) -E-Y ro+ —lY( R)-Y(R-1) I/Q (R) THROUGH PHIL. FOR Ul=ol1 UG.2*M-1 PHi L TE(UtT E(U)+(SI f N. UPOL*Y ( R } )-S I rO*tOLYit - tl)'t:-tQ (R) U't 1 PI) THROUGH- S -UM- FOR MU 0 1, MU.E ~ M ADJaC V( I )+MU) THROUG'-SUMi — FOR NU= O 1, NU E.M S5S ABS (NU-MU)'S6aNU+MU+ IS 1SlI+AJ*A (V( )+NU)*( E(S5)+E(S6) ) -THROUGH SU^, — -FOR J=l 1 J G:.NS3=2*M *(J-1 ) T=ADOJ*A IVJ) i)*(B( 53+55):-B (S3+S6)NORM=NORM+T*XI (I) TfROGsu — p-OR UR 1,1, lU..Gl SUM 52(U) S2(U)+T*SQ((U-1)*N+I) -TROUGHR't-Fi-'-FOR - tU 1 tl i —,.t-ND RLF 52(U)=S2(U)*BETA(U)/NORM -PRtNT-'FORMAT BEFF — 52-( l tr-.-S2-t-ND PRINT FORMAT RECL, 51/NORM PRNT FORMAT MTEN, NORM VECTOR VALUES NDELZ $I2*S V-CTOR VAtUES- PDEL=S lHO i-i-,2'H- DEtAY-E -tEt-RON-GtS ftOUPS*VECTOR VALUES BETa$6F10O6*S VECT.OR VAt-UEt S PBET $ 26 HO Di-AY ED- UTRf-FA~Ti'ONS' F- c6 - VECTOR VALUES DEL=$4F10.4*$ -'ECTOR VAt:-DLP=$S12H ELAY GUPI39H SPECTRU. F15 $ -ftfeE-tMEt tw1; -! _ _ _ VECTOR VALUES GRP= $22HOVELOCITIES FOR GROUP I2*S -VECTOR VALtUES "COST= $ ( 6E12-i5-tS VECTOR VALUES PQ- $19HOVELOCITY IN REGION1392H aE155*S5 -VEC-TOR-VAt-LUE — B-EFF.= S1 $OHOBETA -tFF-. —FF6E-15.-6*$ VECTOR VALUES DEN=$12HONORMALIZER=E15.5*$ VECTOR VALt!S'-RECtS25.HOPROMPT s, iEUTR LlP icFET lME=E.ltL5.6* EXECUTE SEQPGM. -Et D- OF- PROGRAM* BREAK -D-ATA

BIBLIOGRAPHY 1. Goertzel, G. and Garabedian, H. L., "A Method of Solution of the Critical Mass Problem for a Thermal Pile with Slowing Down Properties Independent of Position," ORNL-30, (1948). 2. Edlund, M. C. and Noderer, L. C., "An Harmonics Method Applied to D20 Moderated Reactors," CF-54-3-120, (1954). 3. Garabedian, H. L. and Leffert, C. B., "A Time-dependent Analysis of Spatial Flux Distributions," Nuclear Science and Engineering, 6, 1 (1959). 4. Holway, L. H. Jr., "Perturbation Methods for Reactor Diffusion Equations, Nuclear Science and Engineering, 6, 191 (1959). 5. Foderaro. and Garabedian. L., "A New Method for the Solution of Group Diffusion Equations," Nuclear Science and Engineering, 8, 1 (1960). 6. Foderaro, A. and Garabedian, H., "Two Group Reactor Kinetics," Trans. Am. Nucl. Soc., 3-2 (1960). 7. Churchill, R. V., Fourier Series and Boundary Value Problems, McGraw-Hill Book Co., 1941. 8. Thrall, R. M. and Tornheim, L., Vector Spaces and Matrices, John Wiley and Sons, Inc., 1957. 9. Jeffreys, H. and Jeffreys, B., Methods of Mathematical Physics, Cambridge Press, 1950. 10. Meghreblian, R. V. and Holmes, D. K., Reactor Analysis, McGraw-Hill Book Co., 1960. 11. Morse, P. M. and Feshbach, Ho, Methods of Theoretical Physics, McGraw-Hill Book Co., Inc., 1953. 12. Hartree, D. R., Numerical Analysis, Oxford at the Clarendon Press, 1958o 13. Kohn, W., Journal of Chemical Physics, 17, 1949, p. 670. 14. Bodewig, E., Matrix Calculus, North-Holland Publishing Co,, 1956. 15. Churchill, R. V., Operational Mathematics, McGraw-Hill Book Co, 1958. 16. Gross, E. E. and Marable, J. H., "Static and Dynamic Multiplication Factors and Their Relation to the Inhour Equation," Nuclear Science and Engineering, 7, 281-291 (1960). -89

-9017. Sangren, W. C., Digital Computers and Nuclear Reactor Calculations, John Wiley and Sons, 1960. 18. Gyorey, G. L., "The Effect of Modal Interaction on the Xenon Instability Problem," Trans. Am. Nucl. Soc., 4-1 (1961). 19. Goertzel, G., "Minimum Critical Mass and Flat Flux," Journal of Nuclear Energy, 2, 193-201 (1956). 20. Wilkins, J. Ernest, Jr., "Minimum Total Mass," Nuclear Science and Engineering, 3, 229-232 (1959). 21. Ussachoff, L. N., Proc. Intern. Conf. Peaceful Uses of Atomic Energy, Geneva P/656 (1955). 22. Henry, A. F., "Computation of Parameters Appearing in the Reactor Kinetics Equations," WAPD-142 (1955) ~ 25. Arden, Galler, and Graham, "Michigan Algorithm Decoder," The University of Michigan Press (1960). 24. ANL 5800, Reactor Physics Constants. 25. Marlowe, O. R. et al., Wanda - "A One-dimensional Few Group Equation Code for the IBM-704," WAPD-TM-28 (1957). 26. Flatt, H. and Baller, D., "Notes on the AIM-6 Diffusion Equation Code," Unpublished. 27. LAMS-2255, "Neutron Cross Sections for Fast and Intermediate Nuclear Reactors." 28. Stevens, C. A. and Zweifel, P. F., "Multigroup Diffusion Equations with Arbitrary Space Dependent Coefficients," Trans. Am. Nucl. Soc., 4-1 (1961). 29. Varga, R. S. and Martino, M. A., "The Theory for the Numerical Solution of Time-Dependent and Time-Independent Multigroup Diffusion Equations," Intern. Conf. Peaceful Uses of Atomic Energy, Geneva P/1541 (1958).