THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING THERMAL UTILIZATION AND RESONANCE INTEGRALS IN HETEROGENEOUS REACTORS Muhammad Ao Mannan A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Science in the University of Michigan Department of Nuclear Engineering 1963 January, 1964 IP-655

Doctoral Committee: Professor Paul Fo Zweifel, Chairman Assistant Professor Geza Lo Gyorey Professor William Kerr Professor Richard K. Osborn Associate Professor Franklin H. Westervelt

ACKNOWLEDGMENTS Many individuals and organizations have contributed to the completion of this work, It is with great pleasure that the author takes this opportunity to express his sincere gratitude: To Professor Paul F. Zweifel, the chairman of the doctoral committee, for his invaluable guidance, supervision and encouragement throughout the course of this investigation; To Professor William Kerr without whose help it would have been impossible to obtain financial support in order to complete this work, and for his helpful discussion and comments; To Professors G. L. Gyorey, R. Ke Osborn and F. H. Westevelt, also members of the doctoral committe, for their helpful suggestions and comments; To Dr. Sidney Yip for his enlightened discussions; To the members of the staff of the Computing Center for their assistance with the programming; To the staff of the Industry Program of the College of Engineering for the preparation and production of this manuscript; To the International Cooperation A ministration (now, Agency for International Development) for granting a fellowship to cover the expenses in the United States of America; To the Pakistan Atomic Energy Commission for granting a leave of absence during these studies and for paying the expenses for international travel; Last, but by no measn the least, to many friends whose encouragements have been a constant source of inspiration. 111

TABLE OF CONTENTS Page ACKNOWLEDGEMENTS................................ *,*...0 iii LIST OF TABLES............................................ vi LIST OF FIGURES...................v........... vii CHAPTER I INTRODUCTION......................................... 1 A. Purpose and Statement of the Problem.............. 1 B. Review of Previous Works........S.............. 3 II SOLUTION OF AGE-DIFFUSION EQUATION.................. 12 A. Rectangular Parallelepiped Cell with a Spherical Fuel Element at the Center of the Cell............ 12 B. Rectangular Cell with a Fuel Rod at the Center of the Cell..................................... 19 C. One Dimensional Slab Cell..................... 24 III THERMAL UTILIZATION................................... 29 A. Rectangular Parallelepiped Cell................... 29 A'. Equivalent Spherical Cell.......5............ 32 B. Rectangular Cell............................ 34 B'. Equivalent Cylindrical Cell..................... 35 C. Comparison of Rectangular Cell with the Equivalent Cylindrical Cell....................... 37 IV COLLISION PROBABILITIES IN A SLAB LATTICE............ 41 A. Solution of Age-Diffusion Equation with Delta Function Source and Sink of Neutrons.............. 43 B. Calculation of First-Collision Probabilities...... 47 V INTERFERENCE BETWEEN RESONANCES IN A SLAB LATTICE..... 54 A. Number of Neutrons Absorbed in the Fuel Element per Unit Volume per Second by the First Resonance.............................. 55 I. Interference between Two Closely Spaced Resonances... 59 C. Computational Details and Results................ 65 VI CONCLUDING REMARKS............................o..... 76 iv

TABLE OF CONTENTS (CONT'D) Page APPENDIX A The Integral I(i,j,k,l,m,n)......................... 78 B The Integral Hz(ijk,l,m,n).................. 83 C The Integral I(k,l,m,n).............................. 85 D Derivation of the First-Collision Probability in a Fuel Lump...........................88 E Derivation of the First-Collision Probability in the Moderator Lump.......................96 F Computer Programs........................... 103 v~~~~~~~~~0

LIST OF TABLES Table Page III-1 Thermal Utilization of a Natural Uranium-Graphite Cell for Various Flux Approximations,................... 38 III-2 Thermal Utilization of a Natural Uranium-Heavy Water Cell for Various Flux Approximations.............. 38 III-3 Comparison of a Rectangular Natural Uranium-Graphite Cell with an Equivalent Cylindrical Cell.........5..... 39 III-4 Comparison of a Rectangular Natural Uranium-Heavy Water Cell with an Equivalent Cylindrical Cell......... 40 III-5 Comparison of a Square Natural Uranium-Heavy Water Cell with an Equivalent Cylindrical Cell.............. 40 IV-1 Roots of Determinants of Several Natural UraniumGraphite Slab Lattices................................ 45 V-1 Potential Scattering Cross Sections..................... 70 V-2 Resonance Parameters for U238.......................... 70 V-3 Resonance Parameters for Th232..,.,,.,,,,,,,,,,,,,, 71 v-4 Resonance Interference between the 116.5 ev and 102.5 ev Resonances for U238-Graphite Cells............. 71 V-5 Resonance Interference between Several Pairs of Resonances in a U238-Graphite Cell.................. 72 V-6 Resonance Interference between Several Pairs of Resonances in a U02-Heavy Water Cell................. 72 V-7 Resonance Interference between Several Pairs of Resonances in a Th232-Graphite Cell................ 73 V-8 Interference between Resonances in Th232-Heavy Water Cell................................. 74 V-9 Resonance Interference in a (U238Th232)-Graphite Cell... 75 V-10 Resonance Interference in a (U238Th232)-Heavy Water Cell......................................... 75 vi

LIST OF FIGURES Figure Page 1 Sketch of a Rectangular Parallelepiped Cell3.... 13 la Semiquadrant of a Rectangular Parallelepiped Cell..... 13 2 Cross Section of a Rectangular Cell................... 21 3 Cross Section of a Slab Lattice2................ 25 vii

CHAPTER I INTRODUCTION A. Purpose and Statement of the Problem In this thesis, a mathematical technique has been developed for solving the neutron (age) diffusion equation in "mixed" geometries. Thus, for example, we are able to obtain rigorous solutions to the diffusion equation for rectangular cells containing cylindrical fuel rods, cubical cells containing spherical fuel elements, etc. The reason for developing this technique is that in this way a rigorous* expression for thermal utilization and resonance escape in heterogeneous reactors be obtained. Most calculations(l-3) of these quantities are based on approximations which attempt to avoid the complexities introduced by mixed geometry. Thus, for example, equivalent cylindrical cells replace the actual cells, and we attempt to evaluate the error introduced in this way. More seriously, in resonance escape calculations, it is assumed that the flux recovers completely to a constant in space and lethargy between resonances. Our method permits an accurate calculation of resonance escape without this approximation, and permits us to study the interference between closely spaced resonances. The various approximations which we are attempting to evaluate are described in more detail below. * Rigorous in the context of the diffusion approximation. -1

-21. It is generally assumed, in thermal utilization calculations, that rectangular cells with fuel rods at the center can be replaced by cylindrical cells for the purpose of flux calculations. The assumption that the reactor can be divided into cells, with a zerocurrent boundary condition at the surface, is clearly correct only for a very large system. However, we accept this assumption, but test the additional approximation that the actual cell shape may be replaced by a cylinder. We have done this by calculating the thermal utilization exactly for a system composed of rectangular cells, and comparing it with the same calculation carried out in the "equivalent cylindrical cell" approximation. The methods of our calculations are described in Chapters II and III. The numerical results, obtained via a code described in Appendix F, are presented in Chapter III. Our rigorous calculations have also been carried out for one- and three-dimensional systems, although no comparison with approximate methods are made for those cases. To anticipate the numerical results, we find that the "equivalent cell" approximation is rather good. 2. In resonance escape calculations, it is usually assumed that the flux feeding neutrons into the resonance is constant in space and lethargy. In case of closely spaced resonances the flux may not recover completely before reaching the next lower energy resonance and thus the flux may depend on both space and lethargy. When a resonance is in the transient region of the flux due to the presence of a higher energy resonance, the resonance integral of the lower energy resonance may differ from that calculated with the flat flux assumption. We investigate this effect,

-3the resonance interference, by quasi-rigorous calculations described in Chapters II, IV and V. Although our numerical results are restricted to the case of one-dimensional (i.e., slab) lattice, the two- or threedimensional cases can be carried out using the expressions derived in Chapter II with proper first-collision probabilities. We find that resonance integral changes due to the presence of a nearby resonance at higher energy. Therefore, the flat flux assumption may not be valid for close resonances. In order to understand better the purpose and methods of calculations described here, it is important to review the previous work which has been carried out. B. Review of Previous Works The thermal utilization, defined as the fractional number of thermal neutrons absorbed in the fuel, for a homogeneous system depends upon the composition of material, but in a heterogeneous system it depends also upon the size and shape of the unit cell, The effect of cell shape on thermal utilization has been considered by several authors.(4-11l) The equivalent unit cell approximation of Wigner and Seitz,(l2) which was originally conceived for the calculation of wave functions of crystal lattices, are usually used in calculating the flux in the unit cell. In this approximation it is assumed that the unit cell may be replaced by a cylinder (or a sphere, in case of a three-dimensional cell), although it is clear that a collection of cells cannot occupy the exact volume of any system unless the cells have straight sides. Nevertheless this assumption has become popular for the simplicity of the calculations. The heterogeneous method of Feinberg(13) also gained general acceptance.

-4Weinberg(4) investigated the thermal utilization of a square lattice and compared with that of a hexagonal lattice. He observed that the difference is negligibly small when the radius of the equivalent cell is small compared with the moderator diffusion length. Clark and Newmarch(5) considered a finite-sized fuel rod in a square cell with a capturing moderator. They concluded that the thermal utilization calculated with the cylindrical cell approximation differs from their calculations only by a few tenths of a percent even for water lattices. Cohen(6) considered a square cell with a non-absorbing moderator and replaced the fuel rod by a line sink. He observed that the flux pattern is cylindrically symmetric to a radius of about one-fourth the lattice spacing and showed that a sufficiently accurate value for the thermal utilization is given by the equivalent cylindrical cell approximation. Galanin(7) assumed a line fuel rod and small moderator absorption to find thermal flux patterns and thermal utilization of several lattices. He summed over lattice arrays of modified zero order Bessel functions of the second kind (i.e., the Ko, according to Watson's(l4) notation). Galanin's approach is conceptually the same as the small source theory of Horning,(l5) who showed that the small source theory holds good if the ratio of the fuel rod radius to the diffusion length is small. Neumann(8) combined the lattice sum technique of Galanin(7) with the harmonic development of Clark and Newmarch(5) to obtain expressions for the thermal flux and the thermal utilization for several cells. He observed that the propagated effect of the fuel rod in a lattice will

-5be cylindrically symmetric when the lattice is highly symmetric, when the ratio of lattice spacing to moderator diffusion length is small, or when the ratio of fuel diameter to the lattice spacing is small. Bailly du Bois(9) replaced the fuel rod by a cylindrically symmetric line sink, and solved the diffusion equation for a nonabsorbing moderator. He expressed the thermal utilization in terms of Jacobian theta functions. Pazy and Goshen(lO) extended his method to the absorbing moderator case. All the works mentioned above used diffusion theory and assumed a flat source distribution in the moderator. Horning and Galanin showed that these approximations should be satisfactory for all shapes. Amouyal et al. (11) have developed an improved method for computing the thermal flux disadvantage factor and the thermal utilization in lattice cells in which the diffusion theory assumption was relaxed. They assumed that diffusion theory is applicable a few mean free paths away from the rod, and integral transport is used in the fuel. The method may break down if the moderator region is only one or two mean paths thick. Amouyal et al, concluded that the method gives results as good as those obtained from the P3 or P5 approximations. Our approach has been somewhat different from all of those described above. We assumed diffusion theory to be valid in both fuel and moderator, and expand the flux in the cell in a Fourier series. In this sense, our work is similar to that of Clark and Newmarch,(5) although it is basically simpler. For example, we are able to satisfy the zerocurrent boundary condition exactly because of the form of our expansion. Thus, we do not need to consider, as Neumann(8) does, the contribution

-6of other cells. In this way, we are led to an infinite determinantal equation to solve for the expansion coefficients of the flux, although we find that a reasonably small number of terms represent a good approximation. The same technique is used both in the thermal utilization and resonance escape problems. The basic procedure is as follows: In Chapter II we solve the age-diffusion equation for rectangular parallelepiped, rectangular, and slab cells. Fluxes have been expanded in terms of the solutions of the Helmholtz equation with lethargy dependent coefficients. In Chapter III expressions have been derived for the disadvantage factors and the thermal utilization for a rectangular parallelepiped cell with a spherical fuel element at the center, and a rectangular cell with a cylindrical fuel rod at the center. The disadvantage factors and the thermal utilizations have been computed for several rectangular cells and compared with those computed for the equivalent cylindrical cells. A flat source distribution of thermal neutrons in the moderator and no source of thermal neutrons in the fuel have been assumed. Turning to the question of resonance escape, previous works have been along the following lines. Basically, we represent p, the resonance escape probability, by I' where Ii is the resonance integral of the i-th resonance. where ~l is the potential microscopic scattering cross section of the absorber, and 0S2 is the microscopic scattering cross section of the moderator per absorber nucleus. Therefore, we are interested in evaluating resonance integrals.

-7In order to simplify the resonance integral calculations Wigner(l6) suggested two approximations: (i) the narrow resonance approximation (NR), in which it is assumed that the resonance line is extremely sharp and that a single collision with the absorber is enough for the neutron to pass the resonance region; (ii) the infinite mass (IM) approximation in which it is assumed that neutrons are not slowed down by the absorber, The distinction between these two extreme approximations is made with the help of the concept of practical width, defined as the interval over which the resonance cross section is larger than the potential cross section. When the practical width is small compared to the maximum energy change per collision, the NR approximation is used. The IM approximation is used in those cases where the practical width is larger than the maximum energy loss per collision, But there are cases where none of these approximations gives a satisfactory result. First order corrections to the NR and IM approximations have been made by Spinney,(l7) Chernick, and Vernon, (18) and Rothenstein.(l920) Nordheim(21-23) has developed a numerical method of evaluating resonance integrals without the NR or IM approximation. Recently Goldstein and Cohen(24) developed an analytical method (X-method) of evaluating resonance integrals for homogeneous media. The constant X, which depends upon the resonance parameters only, has a value between 0 and 1 and the extreme values correspond to the usual NR and IM approximations. Hill and Schaefer(25) have extended the X-method for the heterogeneous case. Sumner(26) has written a code based on the X-method. In this work we have considered the NR and IM approximations only. However, the — method may

-8be used in our calculations for the resonances in which the maximum energy loss per collision of the absorber and the practical width have nearly the same value. The usual procedure for calculating resonance integrals in heterogeneous systems has been to assume a flat neutron source (both spatially and lethargy wise) feeding into resonance. Then the resonance integral can be expressed simply in terms of constant source first collision probabilities such as those calculated by Case, et al.(29) While the flat source assumption is valid for resonances spaced sufficiently widely in lethargy, we are interested in finding corrections for closely spaced resonances. We are also interested in finding out whether the spacing in lethargy which may be considered wide depend upon the size of the cell. Schermer and Corngold(27) investigated the interference between resonances in an infinite homogeneous media using a variational technique and concluded that the interference is negligibly small. Brown, et al. (48) performed an experiment to measure the interference between the resonances of gold, indium and rhenium. They placed a spherical lump of indium mixed with gold or rhenium in water. They observed a decrease in resonance integral of indium resonance due to the presence of interfering resonances. Foell, et al.(49) calculated the resonance absorption in gold and indium lumps using a Monte Carlo code(50) and compared with those calculated using the ZUT-MOD-3 code,(51) which considers two overlapping resonances. In our calculations we consider two non-overlapping resonances only. In a recent paper Iijima investigated the resonance disadvantage factor due to (i) the incomplete recovery of flux at off resonance

-9energies and also due to (ii) failure of narrow resonance approximation in the moderator. The error, due to the first cause, of about two percent or less for rectangular uranium-graphite cell was observed. The error in the surface absorption due to the second cause may be sizeable. In this calculation Iijima used the Wigners rational approximation for the escape probability. No correction was made for non-flat source distribution. In our investigation we study the interference between close resonances in a heterogeneous lattice. Our procedure has been to first find the neutron flux, as a function of energy and position, in the region of a given resonance due to the "negative source" contribution of a higher resonance. This is done by utilizing the solutions of the age-diffusion equation described in Chapter II. The results are given in Chapter V. In order to evaluate the resonance integral it is necessary to find the first-collision probabilities for spatial distributions of the neutrons above the resonance. Case, et alo(29) have calculated the flat source collision probabilities for various geometries. Wigner's(l6) rational approximation gives consistently higher value for the collision probability. Sauer(30) has developed an expression for collision probability based on Wigner's rational approximation and has shown that it gives a better result. Di Pasquantonio(1) used a different approach in developing analytic expressions for the collision probabilities for lattice cells composed of different media. But in all these calculations the flat source approximation has been used. Corrections for dense lattices, which were first investigated by Dancoff and Ginsburg,(32) and were later developed by

-10several other investigators,(33,34) are usually made in the calculation of resonance integrals. A simple equivalent expression for escape probabilities in dense lattices has been given by Bell(35) Levine(36) computed resonance integrals by the Monte Carlo method and compared with those calculated from analytical expressions and concluded that the flat source collision probabilities are inapplicable in thick lumps. In Chapter IV we derived expressions for first collision probabilities in a slab lattice due to flux distribution which is a sum of a spatially constant term and a cosine term. The cosine term appears due to the presence of first (higher energy) resonance absorption which has been replaced by a negative source. (We have shown that the next higher terms in the expansion of the flux may be neglected.) If the separation between resonances is large, the coefficient of cosine term vanishes. Since we have to compute the collision probabilities numerically, we cannot prove the reciprocity theoremo(29) In Chapter V we have calculated the resonance integrals of two resonances separately, and also in the presence of one another, All resonances are assumed to be resolved. The unresolved resonances may be treated by the methods developed by Nordheim,(21) and Adler, et al.(37) The interference between resonances has been computed for several selected pair of resonances of U-238, Th-232, their oxides and their mixtures. We find that the resonance interference increases with the increase of thickness of fuel lump. We also observe that the resonance integral of the lower energy resonance increases for the graphite moderated lattices but it decreases for the heavy water moderated lattices,

-11although the total absorption by the lower energy resonance decreases due to the presence of higher energy resonance in both the cases. Hayes, Luming, and Zweifel(38) investigated the cases where the flux should take its asymptotic value below a resonance. This work is a further extension of their method. A part of our results was reported earlier.(39)

CHAPTER II SOLUTION OF AGE-DIFFUSION EQUATION In this chapter we solve the age-diffusion equation in (A) a rectangular parallelepiped cell with a spherical fuel element at the center of the cell, (B) a rectangular cell with a cylindrical fuel rod at the center of the cell and (C) a one dimensional slab cell. These cells are assumed to be one of many similar cells of infinite extent so that zero current boundary condition can be assumed at the surface of the cell. (In a reactor this assumption is valid for the cells not close to the boundary of the reactor.) It is also assumed that: i. Age-diffusion theory holds both in the fuel and the moderator region. ii. The cross sections are independent of energy and are region-wise constant. iii. A uniform monoenergetic neutron source S1 is present in region 1 (fuel) and S2 in region 2 (moderator). A. Rectangular ParallelepipedCell with a Spherical Fuel Element at the Center of the Cell Let us consider a rectangular parallelepiped cell of length 2a along the x-axis, 2b along the y-axis, and 2c along the z-axis with a spherical fuel element of radius p at the center of the cell (Figure 1). Let us choose the center of the cell as the origin of the coordinate system. For the steady state condition the age-diffusion equation can be written as -12

-135z 2b Figure 1. Sketch of a Rectangular Parallelepiped Cell. C Figure la. Semiquadrant of a Rectangular Parallelepiped Cell.

_ I Aft Af a;+D`at) BD aQ d a + Za -Zs;- + S(,x') 6(U) (2-1) (x,y,z,u) is the usual neutron flux, E c(xy,z) and E s(x,y,z) are macroscopic absorption and scattering cross sections respectively, D(x,y,z) is the diffusion coefficient, ~(x,y,z) is the average logarithmic energy decrement per collision and 6(u) is Dirac's delta function. Thus S(x,y,z) 6 (u) is a source of monoenergetic neutrons of lethargy u = 0, and is assumed constant in each region (as are the various material parameters). Since we propose to enforce a zero current boundary condition at the surface of the cell, we expand the flux in the following way: Q oC / 00I c (xYZ~u) = -~ot )O+ If >(l) =-o j-o k=o CoS TX s J7TY COS (2-2) cos COC The summation is over all values of i,j and k except the case i = j = k = 0. We insert Equation (2-2) in Equation (2-1), multiply by cosT i, divide by EIS and integrate over the cell. After taking a Laplace transform with respect to u, we obtain (A = transform variable)

lct23.. LOOO() [bc (k + SO )s) ).orn sort, + 6I / S 2.. IS ) ~~~~where ~ / OC / / i-O j=o k=o [L't ( t> ) + (JbLT. (O0O, O, n m, n) otherwise + ( Ij ( i) rL) +xabl\;;o i, t3Skj = S(ulo' ) (2-3) where I(O,O,O,l,m,) = (2-4) I'(O,O,0,1,m,n), otherwise,ij etc. are Kronecker deltas, l,m,n = 0,1,2,3,... etc. Aoo = 8 A o- =Ao0o ==Aoo, = 4, for i,j,k O0 Aijo =A [0 =Aojk = 2, for i,j,k A 0 nAiJj= 1, for i,j,k A 0 (2-5) and I(i,jkl,m,n) = COS CoS Cn'O C So CO ce c OS b J (2-6)

-16~Hx(i,~jkl~.")7,mn) D sin cos J rcos kirr cUiL 5Urx mrry c TTo -dydE co x Cos CoS (2-7) Hy(i,J,k,l,m,n) = Co s-D COS i t LX sIn co5 CoS COS m Y CIOS rL7r F Y0 k (2-8) Hz(i,J,k,l,m,n) = 1 ) C-og a TX c ry *ir COS — rt os - — g - LTrx CC''Vo iTTE d,Xj (2-9) S(1,m,n) = _ (Xy, ) X M t -Y ceZll Cos C.. 7 X xy A~ c, (2-10) The integrations defining I(i,J,k,l,m,n), Equation (2-6), have been carried out by dividing the cell into several parts such that in each

-17part is constant and can be taken out of the integral, These integrals are evaJluated in Appendix A. We find I(i,~,klrn) ( j L- r; F(ipj llm n) + I(ijk,#1n,mn) + (it,j,k,l,-m,n) + F(i,j,k,l,wm,n) + (ijk,-l,m,)-n) + F(i,j,k,-Zl,-m,n) + F(i,,k,l,-m,-n) + F(i, k, -l, -m,-n)] + b2.L abcAijk i jjm lkn (2-11) Here F(i,j,k,l,m,n) is given by F(i,J,k,lm,n) = 3/ / + J3/2 /.s ae J3/2 is a Bessel function. To carry out the integrations defining Hx, Hy and Hz let us note that is a product of a step function

-18and a delta function oiS DB x 3 (D-D2 ) (x- - D (2-13) Thus, these integrations which contain,and can be carried out easily. The integrations are carried out in Appendix B. We find Hx(ipjJk,1,mpn) -. + DT2) - D1)[ ( t)(i, j Ik mn) + F(i,j,k,l,m,-n) + F(i,j,k,l,-m,n) + F(i,Jk,-l,m,-n) + F(i,J,k,-l,-m,n) + F i Jk 1 m-n ](2 -14) Similarly we obtain Hy(i,,klmn) - +. (D2 - D) [ ( ) jF(ij,k,l,m,n) + F(i,J,k,l,m,-n) + F(i,J,k,-l,m,n) + F(i,J,k,-l,m,-n) + (j-L ) 7r F(iJ,kl,-mn) + F(i,j,k,1,a-m,-n) + F(i,J,k,-1,-m,n) (2-15) + F(i, j,k, -!,-m,-n) ](2-15) and Hz(i,j,k,l,m,n) = c tI + (D2 - D1) [( C Fijklmn + F(i,,k,l,-m,n) + F(i, J,k,-l,m,n) + F(ij,k,1,-mn) +F(i,,k,,m,-n) + F(i,j,k,-l,m,-n) + F(i,j,k,l,-m,-n) + F(ij,k,-1, -m,-n)J ('-) (2-16)

-19The integrations containing the source term can be performed in the same way as the integrations containing Thus we get S(1,mn) 8 I S F(0 O,O,,l0,,m,n) +'s'.- 8abcsb0.S4 o ( Si S. 4 7T am sxs 8 /~ SOL so son (2-17) where 0, when 1 =m = n =0 F(0,0,0,l,m,n) = (2-18) F(0,0,0,,lm,n), otherwise I'(i,j,k,l,m,n) is exactly same as I(i,J,kl,m,n), (Equation (2-6)), except that D1 and D2 should be replaced by Yal and Xa2 respectively. We can solve the matrix Equation (2-3) and obtain the expansion coefficients aiJk(u). In order to calculate the thermal utilization and disadvantage factor, we may assume that S = O, ES = 1 and X = 0. In other words it is assumed that monoenergetic neutrons are produced uniformly in the moderators only and there is no slowing down. Absorption and diffusion are the processes by which neutrons are removed from the system. The energy of the neutrons is the average thermal energy and the measured valued of material parameters clan be used in the calculations of disadvantage factor and thermal utilization. B. Rectangular Cell with a Fuel Rod at the Center of the Cell Let us consider a rectangular cell of length 2a along the x-axis and width 2b along the y-axis with a fuel rod of radius /0 at the center

-20of the cell (Figure 2). The cell is assumed to be infinite along the z-axis. For the steady state condition the age-diffusion equation is written as -D(x,y) + - - ) ]' + ~ )~( ) ( O X );3Y a y = g (X. vy U( S(XOY) S(Nl (2-9) All the quantities are defined as in the previous section. In order to satisfy the boundary conditions,'we expand the flux C (xy,u) in the following way: q. (xYyu) ZL + ak,(u)cos - coSs Ab (2-20) a —o La* The sunation o is over all values of k and 1 except the pair 1 = k = 0. We insert Equation (2~-20) in Equation (2-19), multiply by Cos WICrXCos- and divide by Then upon integrating over the cell, we find after taking Laplace transform ( ~ transform variable) the following equations for the expansion coefficients akl *:o ( _),( n nr) + (.)H?(kLtfnn) S ~ 2 t- (2- 21,) + Q1. It. $.~ = $(~. ~) (L- z)

-21y _'I' Fgr 2 Cosetoo2a g Figure 2. Cross Section of a Rectangular Cell.

-22where _, (O.0, when m = n = 0 I(O,O,m,n) = (2-22) I'(O,O,m,n), otherwise m,n = 0,1,23,ooo.. etc. AO0 =4 A ko =Aot 2, for k,1l 0 At =1, for k,l O 0 (2-23) and = 4 _,tx ~ {.} + c,. - -i +a ~~$~'~ ](2-24) O -D ITTX LrTrY Hx(k,l,m,n) = I s COS C/OS V'L~TX rL7TtY aX 7Y (2-26) -Hy(kJ l mn) = e cos tL cX s; Y[ SZsl r:~ -- COS- Ob t jLZS2-dxdy a. S b y(2-26)

S-rvr. /z 5 Ss.,mY) Cos Cos b' 7T S dx Ay =s cos t k (2-27) After performing the integrals in Equation (2-24) (see Appendix C), we find I(k,~l,m,n) =2 ( - i) [F(k,l,m,n) + F(k,l,-m,n) + F(k,l,m,-n) + F(k,l,-m,-n) ] + abA kl Smk Snl (2-28) where F(k,l,m,n) = J + ( 2[ L (2-29) Jl(.') is a first order Bessel function of the first kind. As before let us note that i4 - D is a product of a step function and a delta function. S A1 (52 S 5 - (2-30) Here is the unit step function and D _~ (Dl - D2) (x 2 ) (2-31) Evaluating the integrals in Equation (2-25) we find Hx(klmn)= (S + (D2 - Dl) [ a (F(klmn) + F(k,l,m,-n)} + (~ - ) I ~F(k,l,-m,n) + F(k,l,-m,-n) }] (2-32)

Similarly, we get / 1 _'[ 6 Hyki)=1 + (2 - O1) (k,l,m,n) + F(k,l,-m,n)} + (L-) F(k,l,m,_n) + F(kl,-m,-n) 3 (2-3) I'(k~l,m,n) has the same expression as I(kl,m,n), Equation (2-28), except that D1 and D2 should be replaced by Eal and Ea2 respectively. S(mn) = 4 (F(0,0,m,n) + - X l +4, SA }- om S.. (2-34) where f0, when m = n = 0 F(O,O,m,n) = (2-35) LF(O,O,m,n), otherwise Solving the matrix Equation (2-21) and taking inverse Laplace transform, we get the expansion coefficients akl(u). C. One Dimensional Slab Cell Let us consider a slab cell of width 2b with a central fuel lump of width 2a (Figure 3). For the steady state condition we write the age-diffusion equation as - + SX)C (xu)) X+ (x6, --- -~ (x) a;*x. + scx) S(u) (2-~6)

1' IX I I I' 1 1 I I I FI 1~~~~~~~~~~~~~~ I II I ~~~MODERATOR I i I 1 I I III~~~~~~~~~~~~~~~~~~~~~~~~~~~ I I I I I I III~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1 I I I II I I I I III I I I I I I I I I I I I I I I I I I 1 I I I I III~~~~~~~~~~~~~~~~~~~~~~~~~~~ I I I I ~~~~~~~~~~~~~~I'I I I I I I I'I I I i~~~~~ I _ _ _ _ I I + ~~ ~~ ~~~~~~~~~~I -I I I I I I I I Figure~~~~~ 5.CosScIono lbLtie

-26The one-dimensional case can be treated either by the method used in the two and three dimensional cases, i.e., by integrating over the singularity in or alternately, by measuring the distance in units of diffusion coefficient, D(x). In the present section we adopt the latter approach. Introducing a new variable X / X 1 a x (2-37) we write Equation (2-36) in the following form _______ Il.'LA)(y,u) + +(,L) "a L S. () $(1) (2-38) Here we use the relation Di = (2-39) where / is the average cosine of the scattering angle per collision in the laboratory coordinate.system" The boundaries x = a and x = b correspond to y = ( and y =. They are given by C( = XL Ti (2-40) - __

-27As before we expand the flux OC (Y ) Glo(U- a (=, u 0 cos " Tr2-Y 1Y... ) (2-41) n) -I We insert Equation (2-41) in Equation (2-38), multiply by cos n rfa and integrate over y from -- to + o Taking the Laplace transform, we get ac.OS[( gs_,) $; (nI h7Tr ((3 + +A\ i +, a)~ =1 + a_ e,- l2 + (VJI -WZ) )~i~Ts }~3 (D"' 6,n ~14g +wz( ~r27i~Tj Fs + \- g C) (2-42) where m = 0,1,2,3, etc. gi 5; ZsC (2-43) 3CW- t).t

-28and F(m) = s ( - )+ dL sOm (2-44) For the particular case when there is no absorption during slowing down, Fa = 0 and Equation (2-42) becomes ao(%N) Aom + jann( ) [ Am + mn] (2-45) where Amn T T( 5) [(wl( W2)i +) + < rn In + uW2r T rn (2-46) We can take account of resonance absorption during slowing down by replacing the resonances by delta function sinks. Solving Equation (2-45) for an( A) and taking the inverse Laplace transform we obtain the expansion coefficients for the case in which there is no non-resonance absorption during slowing down. Since at high energy resonance absorption is the most dominant absorption term during slowing down, the neglect of non-resonance absorption is a good approximation.

CHAPTER III THERMAL UTILIZATION In this chapter we obtain expressions for the thermal utilization and disadvantage factors for (A) a rectangular parallelepiped cell, (A') an equivalent spherical cell having the same volume as the rectangular parallelepiped cell, (B) a rectangular cell, and (B') an equivalent cylindrical cell having same cross sectional area as the rectangular cell. For the cases (A) and (B) the fluxes obtained in Chapter II are used in the calculations of the average fluxes. But for the cases (A') and cB') the diffusion equation is solved separately in the fuel and the moderator regions with the boundary conditions that the current and the flux of thermal neutrons are continuous at the fuel-moderator interface. In all cases a flat source of thermal neutrons in the moderator is assumed. In Section C numerical values of the thermal utilization and disadvantage factors are computed for a rectangular cell and compared with those computed for the equivalent cylindrical cell. The disadvantage factor is defined as the ratio of the average flux in the moderator region to the average flux of thermal neutrons in the fuel region of the cell. A. Rectangular Parallelepiped Cell The flux for this cell is given by Equation (2-2) oC oc, o0 (XYi, _A)= %, +II oL c iorX Ls-o _o lao.co~s:~no, cos ia TC S ~(3-1) -29

-30 - In order to calculate the thermal utilization aijk are obtained by solving the matrix Equation (2-3) with A\ = 0, 1 Z5 = 1 and S = 0 When the expansion coefficients, aijk, are known, the average fluxes in the fuel and the moderator region can be found by integrating over and dividing by the volume of fuel and the moderator respectively. Thus, the average flux,4l, in the spherical fuel element of radius P is 4M- 8' + rr. """ =+ 3ZZ R jJ COS Tr O l=o o o xCos 7 Y Cos kt,.X 4sy Aj b C 2ooo' = 7ro + F3 ktf(;ip, 06) (3-2) L=o j=3 k=o Here F(i,j,k,O,O,O) is defined as in Equation (2-8). The average flux, LP, in the moderator is I 6( i c 4 ) E ) C; 7 3 C o j=o k-co ( C0o0 o05 j 0xd0 od: ~~mos f o o~bc 6 g- j-o =o

-31The average flux on the surface of the fuel element is Cre?. o~'tltt cos;'x L-o $-o k-o =so.wr LsCer cos. bY cos k~__ c C /~ X on _ aoo zi = am >>'t'o j o l> o (2 tX [ Ict F (b)~( C ) + (3-4) The above integral can be evaluated easily. The thermal utilization, f, is the ratio of thermal neutrons absorbed in the uranium to the total number absorbed in the cello Let us consider that the fuel element contains natural uranium only. The thermal utilization is usually expressed by the following relation YI __=VI t 4 1 V1 and V2 are the volumes occupied by the fuel and the moderator in the cell. F is the ratio of the average flux on the surface of the fuel sphere to the average flux in the fuel element, i.eo, F = E is characteristic of flux distribution in the moderator, and is given by E - = VzT 2 vt..V... (3-6) V/,1..,

-32The absorption in the fuel, Vit o I~, is equal to the net current into the spherical fuel element. (E-l) is indeed that part of the absorption in the moderator which is due to the excess of the neutron flux in the moderator over the value of the flux at the fuel-moderator interface, divided by the absorption in the fuel, This is sometimes called excess absorption. The disadvantage factor, d, is defined as a = 4/. (3-7) The thermal utilization, f, can be expressed in terms of the disadvantage factor, d. - =1 +.: (3-8) Y_0~! V3 The thermal utilization depends upon the atomic ratio and the volume ratio of the fuel and moderator. It also depends upon the size and shape of the fuel and moderator separately and is introduced by the presence of disadvantage factor, d, in Equation (3-8) A'. Equivalent Spherical Cell The rectangular parallelepiped cell is approximated by an equivalent spherical cell having equal volume. The normal derivative of flux is assumed to vanish on the surface of the cell. This assumption, which underlies the method for computing the cohesive energy in a crystal lattice,(l2) is introduced to simplify the mathematics. The radius, o, of the equivalent cell is - = ( 6c>i C)23 (5-9)

-335In order to find the flux in the cell, we must solve the following diffusion equations: Div 1,(h) -- c-1 1,(t) = (in fuel) (3-lOa) D2 V'2T('t) - Dz a2z( 4 ) + O (in moderator) (3-lOb) The boundary conditions are 4 (e (3 -11) DI 1= PDC t-1 (3-12) Solving Equations (3-lOa) and (3-10b) with above boundary conditions, the average fluxes in the fuel and the moderator can be found easily. Thus, the expressions for E and F can be obtained(2) in the following form:,( X( tvL~ Xt P r {~;~totxp Xi(3-12) E x..-&..F. - X..cottX2.(r-.) 1 (3-13) t3~. e a _ -xaC-X:.(_) cotkx.( _)) where Xl= - 9 ___ (3-14)

-34B. Rectangular Cell For the rectangular cell the flux has been shown in Chapter II to be ((X, Y) = a + C,. cos -X Cos -- 3-15) 4 + (3515) akl are obtained by solving the matrix Equation (2-17) with,. = 0, Es = 1 and S1 = 0. Thus, the average flux in the fuel rod is given by N o k o t o Ps, - + -0', os Cos 4 7T pz I, O. " ~ k=O t= (3-16) F(k,l,O,O) is defined as in Equation (2-22). The average flux in the moderator is x CO +Z5ab_ 0q) X AJyCOS J kxCQS Lrrvdr 4 VOj ii =- 4. - c*o

-35The average flux, 4s, on the surface of the fuel rod is cP' = k Z c L cos krrx c Los-2 Y,=o tLo ct rcle G= too +Y X a A (" Ad (Ct L bV ) _- a1- o (3 -18) The above integral can be evaluated easily following the method as shown in Appendix C. BE Equivalent Cylindrical Cell In this case the rectangular cell is approximated by an equivalent cylindrical cell having the same cross-sectional area. The diffusion Equations (3-10a) and (3-lOb) are solved with flux and current continuity condition at the fuel-moderator interface. It is also assumed that the normal derivative of flux vanishes at the surface of the cell. The radius, 9, of the equivalent cylindrical cell is g ( f (3-19) Solving Equations (3-10a) and (3-lOb) in cylindrical coordinate system it is found(37) that the flux ) (h), in the fuel rod is a(k), =A Io ( t (-20) and the flux, C2(-), in the moderator is

-36h) = C[ii(X2 ) K(x) +K.(X2. )I()X +)] + (3-21) Io0 Il, Ko and K1 are zero and first order, modified Bessel functions of the first and second kind respectively.(l4) The constants A and C may be evaluated by means of the boundary conditions. A1( Xt) =X c[ Ixt X(2 )+K o(X'2C),(x,)] 4- (3-22) and A,r, (xi ) =[ 1,(xO K(X t - (xI4.(X2,) ] (3-23) Tie thermal utilization, f, in the equivalent cylindrical cell is given by. _I = V2 F Ya F +.E-i (3-5) where F 10(Xe) (3-24) and X____ _) Io( XP) K(X ) +K o(XZ j)I1(X2~) 1 E a =X ) Ki(x(Xf ) - K(X ) II(XX,) J (3-25) Xi and X2 are defined as in Equation (3-14).

-37C. Comparison between a Rectangular and an Equivalent Cylindrical Cell A computer program (Appendix F) has been written for IBM-7090 using MAD(42) language to solve Equation (2-21) for akl and to compute the average fluxes in the fuel and moderator, the disadvantage factor and thermal utilization. Let us note that Equation (2-21) is an infinity by infinity matrix. Tables III-1 and III-2 show the average fluxes, disadvantage factors, and thermal utilization of two different cells calculated keeping different terms in the expansion of the flux, From the tables we find that the expansion coefficients converge rapidly. We get an excellent approximation for disadvantage factor and thermal utilization if we keep the expansion coefficients for k, 1 = 0,1,2,3,4o In Tables III-3, 4 and 3 we have calculated the disadvantage factors and the thermal utilizations for three different cells varying the radius of the fuel rod and compared with those calculated with the equivalent cylindrical cell approximation. The computer program written for the equivalent cylindrical cell is given in Appendix F. From these tables we observe that the equivalent cylindrical cell approximation is accurate except for very tight lattices. The thermal cross sections and the diffusion coefficients used in these calculations has been taken from ANL-5800.

-38TABLE III-1 THERMAL UTILIZATION OF A NATURAL URANIUM-GRAPHITE CELL FOR VARIOUS FLUX APPROXIMATIONS a = 2 cm, b = 2.5 cm, p 1 cm Average flux* Average flux* in the fuel in the moderator Thermal Utilization k 2 f 1 62/1 f 1 1 14.039043 16.524909 1,177068 0.995510 2 2 14.060333 16.703461 1.187985 0.995469 3 3 14.075512 16.756275 1.189035 0.995465 4 4 14,090440 16.762192 1.189615 0.995463 5 5 14.099829 16,774762 1.189714 0.995462 6 6 14.105524 16.782078 1.189752 0.995462 7 7 14.111894 16,789801 1,189762 *0995462 * Arbitrary unit TABLE III-2 R UTILIZATION OF A NATURAL URANIUM-HAV WATER CELL FOR VARIOUS FLUX APPROXIMATIONS a-2 }ri b -2 cm p 1 cm Average flux* Average flux* ih tthe fu1 irn thejn@rator A Thermal Utilization k. 2 Q1 2/Ql f 1 1 9.95 61-20 l25057761 1. 2 528 0 999550 2 2 9, 89984 12, 55554, 2187 O, 999548 3 S 9856714 12.059424 1.221444 Oa999547 5 9,840412 12,0025 2.81,209 O, 999547 6 6 9*837569 12.025747 la22241 0, 999547 7 7 9*855 622 12022149 1222556 0999547 4Arbi1trary niLt

-39TABLE III-3 COMPARISON OF A RECTANGULAR NATURAL URANIUM-GRAPHITE CELL WITH AN EQUIVALENT CYLINDRICAL CELL a = 2 cm, b = 2.5 cm, k = Q = 5 p ( cm) ((2/Fl)RC (2/F1)ECC fRC fECC 1.8 1.435236 1.300985 0.999012 0,998699 1.6 1.353071 1.273241 0.998566 0.998568 1.4 1.292296 1.242156 0.997930 0.998404 1.2 1.239316 1.208218 0o996982 0,998189 l o 0 1189714 1.171716 0.995462 0.997893 0.8 1.o 142149 lo 133084 0.992756 0.997454 0.6 11.096617 093305 0o. 987104 0.996727 0.4 1.054568 1.054131 0.971620 0.995273 RC = Rectangular cell ECC = Equivalent cylindrical cell

TABLE III-4 COMPARISON OF A RECTANGULAR NATURAL URANIUM-HEAVY WATER CELL WITH AN EQUIVALENT CYLINDRICAL CELL a = 2 cm, b = 2.5 cm, k = a = 5 p (cm) (72/1)RC (~2/~)ECC fRC fECC 1.8 1.595902 1.325893 0.999861 0.999832 1.6 1.484790 1.303616 0.999800 0o999814 1.4 1.401965 1.272824 0.999714 0.999792 1.2 1.328663 1.238422 0.999588 0.999764 1.0 1.259559 1.199773 0.999388 0,.999726. 8 1.193202 1.157289 0.999033 0.999669 0.6 1.130091 1.111750 0.998294 0.999576 0.4 1.072665 1.065726 O.996243 0.999391 TABLE III-5 COMPARISON OF A SQUARE NATURAL URANIUM-HEAVY WATER CELL WITH AN EQUIVALENT CYLINDRICAL CELL (cm) (~2/D1)RC (~2/~1)ECC fRC fECC 0.8 1.085149 1.070877 0o.999903 0.999863 0.6 1.061513 1.053869 0.999756 0.999820 0.4 1.039283 1.035871 0.999345 0.999735 0.2 1.016171 1.015347 0.997169 0.999481

CHAPTER IV COLLISION PROBABILITIES IN A SLAB LATTICE The resonance integral, I, is defined as =I - 4O l(u)aa(u)du (4-1) Res where r10 is the flux in the fuel element which would exist at the resonance energy in the absence of the resonance. 41(u) is the average flux in the fuel element in the resonance region. The lethargy integral in Equation (4-1) is over the resonance interval. In order to get an expression for l1(u) in the resonance we write the integral equations for +(u) in two media, following Chernick,(43) in the form given below: XPl(u) VLZsijl(ut)eU1u i i U~ti =( 251\1Ju s2e du' i i u.-Ai u'-u (4-2a) E V2Esj2(u')e + (l1-P2(u)) V1i du' u-Aj -41 -

The summations. and. are over the atoms in the fuel and moderator respectively. Zj and Esi are the total and scattering macroscopic cross sections respectively of the i-th atom. l(u) and +2(u) are the fluxes in the fuel and moderator respectively. The maximum fractional energy loss of a neutron in an elastic collision with an atom of mass, Ai, is (1 - ai), where (i (Ail) 24-3) Ai = ~n(1/ai) (4-4) The collision probabilities Pl(u) and P2(u), which enable us to relate the collision densities 1()(L(u) and Z (u)j2(u) to the previous collisions at higher energy, i.e., at lower lethargy, are defined for the lattice in the following way: Pk(u) is the probability that the neutrons of lethargy u originating in region k (k = 1 for fuel, k = 2 for moderator) will make their firstcollision in the same region k (probably after traversing other regions). The integrals on the right hand sides of Equations (4-la) and (4-2b) can be evaluated easily in the NR and IM approximations. The details of the calculations of the resonance integrals will be shown in Chapter V. In this chapter we derive expressions for the first-collision probabilities. Where two resonances are close together the flat-source collision probability is inapplicable in the calculation of second (lower energy) resonance. In our calculations we replace the first (higher energy) resonance by a delta function sink and solve the age-diffusion equation

-43to obtain the energy and spatial dependent flux above the second resonance. In Section A we obtain the spatial and energy dependent flux. In Section B we derive expressions for the first-collision probabilities using the flux obtained in Section A as the source distribution of neutrons. A. Solution of Age-Diffusion Equation with Delta:Function Source and,Sink of Neutrons The age-diffusion equation is exactly the same as Equation (2-29) with Z = 0 and a modified source term. A delta function negative a sQurce, S, at lethargy u w uo is added to the source term. -'- (D(x) O (xu)). -t(x)s(x) a + S(x)5(u)-so(x)5(u-uo) (4-5) As in Section C, Chapter II, the flux is expanded in the following way 4(y,u) - () + an(u)cos C (4-6) n=1 All the terms have been defined in Chapter II. The matrix equation for the expansion coefficients becomes o00 Bo(X)k bom +Z n(X)[Amn + X6mn] = (m)-o(m)e-o (4-7) n=l Equation (4-7) is exactly the same as Equation (2-38) except for the modified source term. In order to find the expansion coefficients, an(u), we have to solve Equation (4-7) and take the inverse Laplace transform. Let us

note that the roots of the following determinant, Equation (4-8), multiplied by u or (u-uo) will appear in the expansion coefficients as negative exponent. X A01 A02... An.. 0 All+X A12 Aln o A21 A22+X... A2n...ee O (4.8) o Aml Am2 Am+X.. Amn m O Anl An2... Amn+X Equation (4-8) is an infinity by infinity determinant. Taking different values of m=n we have solved Equation (4-8) for X. Table IV-1 shows that the roots converge rapidly and the smallest root may be approximated by All, i.e., by the solution of a 2 x 2 determinant. The actual value of the smallest root, X1, is about three percent less than All. Let us also note that the next higher root, X2, is about five times as large as X1l Since Xi appear in the expansion coefficients as negative exponent, we can neglect all the terms involving Xi except for i = 1. Solving Equation (4-7) for an(X) and then taking the inverse Laplace transform we obtain the expansion coefficients an(u) for different approximation. Solving a (n+l) x (n+l) matrix equation we find the expansion coefficients ai(u), where i < n, have the following form

TABLE IV-1 ROOTS OF DETERMINANTS OF SEVERAL NATURAL URANIUM-GRAPHITE SLAB LATTICES Half-width of fuel = a cm Half-width of Moderator = (b-a) cm aU = 10375 barns aC = 4.5 barns P P a b 1 2 3 X4 1126 798.o5 1o 1114.54 5905.7 1093o 19 5905. 1 10999.9 1092o 75 5865.9 10954.9 21876 0 281.7 1.0 2 278.6 1476.4 273~3 1476,2 2749 9 273 2 1466 4 2738.7 5469 0 125.2 1.5 3 123.o 8 656.2 121.5 656.1 1222.2 121.4 651.7 1217.2 2430.7 70.4 2.0 4 69.6 369.1 68.3 369.0 687.5 68. 29 366.6 684.7 1367. 2 33.1 2.5 6 31.7 143,.5 31.6 140o7 301.5 31.2 140o6 298.1 571.8

ai(u) = Cnl e-l(u-u~)+ Cn2 e-2(u-u)+...+Cnn e-n(u-uo) _ Cn1 e-kl(U-Uo) where Cnl...,Cnn are constants. We found in the one-speed case, i.e., the thermal energy group, that the first two terms in the expansion of the flux are correct to about two percent. The expansion coefficients for the higher terms are very small, We, therefore, keep only the first two terms in the expansion of the flux ~(y,u), Equation (4-6), and write (y,u) = a(u) + al(u) cos (4-9) 2 where |x/D1, for -a < x < a = (4-10) a/Dl + x-a/D2, for a < x < b The expansion coefficients are ao(u) = S( ) -01 ( -elU) Al1 - (u-uo) so() lS(All (l-e-l(uuo) (4-11) and al(u) = S(l)e-klu - (u-uo)So(l)) e-Al(u-uo) (4-12)

-47where 1 for u > uo (U-lo) (4-13) 0 for u< u (4-1 ) is a unit step function. The terms containing e can be neglected. The terms containing e 1(uuo) will be important only when Xl(u-uo) g 1. For widely separated resonances %(u-uo) at the second resonance would be large and consequently e -l(u-uO) 0, i.e., al + 0 and as expected flux would be constant in lethargy and in space. Therefore, we find that the flux above the second resonance becomes spatially and lethargy dependent due to the presence of a close first (higher energy) resonance only. B. Calculation of First Collision Probabilities The first collision probabilities are usually calculated assuming a flat source distribution as given by Case et alo(29) We derive expressions for the first collision probabilities for source distributions represented by the spatially dependent flux, C(xu), shown in Equation (4-5). The single collision kernel at x' due to a plane source of one neutron per cm2 per sec located at x' is given,.by 2 E1(Zlx-x'l) (4-14) where En(x) is defined as(44) En(x) ext tndt 11 n-2 (4-15) 0

-48for positive values of x ~ 7 is the total macroscopic cross section of the medium. For different media we have to take the optical distance between the plane source and the point under consideration. Therefore, the probability, Pl(u), that the neutrons of lethargy u originating in the lump (-a < x' < a) of region 1 will make their first collision anywhere in the same region 1 of the lattice is given by (see Figure 3) Y'l a a a Pl(U) = 1 dx (xv) El(ilx-x' I)dx (X')dxv 2 -a -a a 0a 2nb+a + P _ adxE(x') EL[n2(b-a) (Z2-l) n=l -a nb-a + Z1(x-x')]dx O(xj)dx' -a ~p~ 0 a -2nb+a + dx-(x-) dx El n2 b-a)(Z-4) n=l -a -2nb -a a + lxx)/ (x')dx (4-16) -a The above integrations have been evaluated in Appendix D. Finally for Pl(u) we get ao(u) = -E a(()( nK 7 -2( Pl(U) 1..5 - (E3(k(n)) - 2E3(2(n)) + E5(3(n))j n=O

-49al(u) [sin CCa { 1j t Z +? Z(cos 1 + c2) - tan + log 00 + a1(u) [zilcos'afjl(n) - 2 12(n) + W(n) 3 n=O + alsin'ala i!(n) - i 3(n) (4-17) where dl = 2a d2 = 2(b-a) (4-18) e-i (n)/ttdt'-p (n) =0 |t+d2t2'i = 1,2,3 (4-19)' (n ) = e i /tt 1^ dt (4-20) +2 2 tD2 at D2a+Dl ( )4-21) 1 = 2(x)dx' -a

-5oand el(n) = (n+l)d22 + ndlC1,2~(n) = (n+l)(d22 + dl1) 33(n) = nd2Z2 + (n+l)dl,1 (4-22) Similarly, we find the probability, P2(u), that the neutrons of lethargy u originating in the moderator lump (a < xg < 2b-a) of region 2 will make their first collision in the same region 2 anywhere in the lattice. 12b-a 2b-a 2b-a P2(u) = / dx'(x') ELz 2 x-x'|]dx d(x')dx' a k,. —2b-a 2(n+l)b+a.2b-a + / dx'(x'?) El[2na(Zl-Z)+(x-x'),2 dx (x')dx' 2 J n=l a 2nb+a a 00,2b-a — 2(n-1)b-a -2b-a + 2 dx'4(x') E1 2na(Zl-Z )+(xV-x)Z2 dx |(xv)dx' n=1 a -2nb+a a (4-23) The above integrals have been evaluated in Appendix E o We find

-5100 P2(u) 1..- -.....2 3 _d2~2 2 ~ 5- X E3(~l(n))-2E3( 2(n))+E3( 3(n)) n=O 2d2$2 2Pt22 + tan (sin'(2b-a)-sin l'a a1+(u)sin loge 1 + j(sin P'(2b-a)+sin p'a) 2dA2 2I 2' + (1 - tan >(cos'C(2b-a)+cos'a) +/ - a(u)c os (' (2b-a)+cos B'afFl(n) -2I2((n) 2d242 n=O + II3(n)}+ c'{sin P'(2b-a)-sin 5'a) [IIl(n)-Il3(n)] 00 - al(u) z2sin in (2b-a)+sin p'a)] fl(n)-2I12(n)II3(n)) n=O -'(cos P'(2b-a)-cos p'a) IIli(n)-l3(n))] (4-24) where ita (1 1 ),-= (4-26) 1DD2

- 2b-a 2= 2(b (x )dx (4-27) a f'l (n) e tdt i = 1,23 (4-28) 1 -.Pi(n)/ 2 e t dt i (n) = (4-29) 0 l(n) = (n+l)dl71 + nd2Z2.(n) (n+1) (d-i1l + d2 ) 3 (n) = ndlZi + (n+l )d2Z Let us note that the first two terms of P-1 (u) (also P2(u)) are exactly the same as derived by Rotherstein>l9)9 (20) with a flat source assumption. The remaining terms are due to the presence of the cosine term in the flux:, $(xu)> Hence, for t-he flat source assuiilption we get 1 ) PO (u)I 1 ) 2E3 () )+E3 (123 ( n ) n=O (4-31) and P2(u) > P.(u>) = 1 - L- 3(~(n))-2Es( 2(n))+E3(53(n)

-53 - From the reciprocity theorem it can be proved29 that ZlVl(l - P(u)) = 2V2(1 - P2(u)) (4-33)

CHAPTER V INTERFERENCE BETWEEN RESONANCES IN A SLAB LATTICE In this chapter the interference between two close resonances is investigated. It is assumed that the flux takes its asymptotic value before reaching the first (higher energy) resonance. The flat source collision probability is used in the calculation of the first resonance integral. The second resonance (lower energy) integral is evaluated in the absence and also in the presence of first resonance, The lethargy and spatially dependent flux, due to the presence of the first resonance, and the corresponding first-collision probabilities are used in the calculations of the second resonance integral, In Section A, the first resonance integral and the strength of the sink due to resonance absorption are obtained, In Section B, the second resonance integral is calculated in the absence of first resonance and in the presence of a delta function sink due to the higher energy resonance. The NR and IM approximations are made in these calculations. It has been assumed that the Breit-Wigner one-level resonance cross-section formula can be used for all resonances. In Section C, the computational details are discussed and the numerical values of several pair of resonance integrals are shown, Throughout these calculations it has been assumed that the change of temperature is inconsequential in the evaluation of resonance integrals,

A. Number of Neutrons Absorbed in the Fuel Element Per Unit Volume Per Second by the First Resonance In this case we assume that the flux takes its asymptotic value before reaching the resonance, i.e., it is constant in space and in lethargy within at least one collision interval above the resonance, It is also assumed that there is no absorption during slowing down. The asymptotic flux in a reactor lattice is given(45) by Oasy = q(Vl+V2) (5-1) Vl 1Z1+V2 2p2 where q is the slowing down density, 9l and 2 are the average logarithmic energy decrement per collision in the fuel and the moderator region respectively. 4 1 and 2 are the macroscopic potential scattering cross sections in the fuel and the moderator respectively. The average number of neutrons absorbed in the fuel lump per unit volume per second, So, is(46) So = No, Oasy (5-2) Where No is the number of absorber nuclei per unit volume of the fuel (region 1), I is the resonance integral defined in Equation (4-2). In this case 0 10 =.asy In order to evaluate the resonance integral, I, we make the NR or IM approximation and solve the integral Equation (4-2a) to obtain 4l(u) at the resonance lethargy. For the higher resonance the escape probabilities given by Equations (4-31), (4-32) and (4-33) can be used.

1. Narrow Resonance Approximation When the maximum energy loss in a collision with an absorber atom is greater than the practical width of the resonance, the narrow resonance (NR) approximation is applied. Usually the absorber is the heaviest atom in a reactor system. If the NR approximation can be used for the absorber atom, it is expected that the NR approximation is also applicable for the atoms mixed with the fuel element and for the moderator atoms. In other words we can use the asymptotic values of $1(u') and 02(u') inside the integrals in Equation (4-2a) and evaluate the flux, 0l(u), at the resonance. The asymptotic value of the flux is given by Equation (5-1), i.e., 0L(u') = 02(u') asy (5-3) From Equations (4-2a) and (5-3), we find i 1 J (5-4) Let us note that / z' l(U) and /i = ZS are the macroscopic total and scattering cross sections respectively in the fuel element; and E.. - s 2 is the macroscopic scattering cross section in the moderator lump. Let us assume that the fuel element may contain U238, Th232 and O16-atoms. From Equations (4-1) and (5-4) the resonance integral for the i-th absorber, INR, is obtained

asy 1l(u)fu) i( du Oasy Res Z h Cr'(fu)P (uP u du(u)PO(u)du + +1 C u)PO u) a 1 Res Res + |o1(u) du (5-5) Res (Here i = Th or U), the relation, Equation (4-33) between PO(u) and P0 (u) has been used to obtain Equation (5-5). 2. Infinite Mass (IM) Approximation When the practical width is larger than the maximum energy loss in a collision with the absorber atom, the IM approximation is employed. It is assumed that the scattering by the absorber atoms changes only the direction of motion of the neutron but not its energy. Therefore, for this particular case, we assume OU''TOh + 1, which make the integrals due to the absorbers delta functions. But the NR approximation is used for the moderator and the moderating material present in the fuel element, for example, the oxygen in the uranium oxide. As before, we assume that the fuel element contains U,Th and 0. Thus, in the infinite mass approximation Equation (4-2a) becomes

+ = VlZPo(u)l(u) + L VlPO(u) l(u) + VjlasysPl(U) + Zs2V20asy(l-PO(u)) (5-6) Assuming that there is no absorption in the moderator and using the reciprocity relation, Equation (4-33), we obtain the resonance integral for the i-th absorber in the following form: ITM - (z U Th 0 Th U Res ( + + s)) - PT(u)( + ZS) ( + Zh)i a(u)Plu)du J+ ( 4+h + -P _ (u)(5 + Zsh) (5-7) here i - U or Th, Let us note that, and when i = Th, U U E,= l, the potential scattering cross section of uranium, because we assumed that there is no non-resonance absorption during slowing down. Th Now, from Equations Similary, when i = uranium, T Now from Equations ~~ ~ Z

(5-1), (5-2), (5-5) and (5-7) we can compute the strength of the sink So, due to the higher energy resonance at lethargy uo B. Interference Between Two Closely Spaced Resonances In this sectionwe calculate the interference between two close resonances, The absorption rate by the first (higher energy) resonance has been calculated in Section A. In order to calculate the second (lower energy) resonance integral in the presence of the first resonance, we consider the first resonance as a delta function sink, The solution of the age-diffusion equation with this approximation has been obtained in Chapter IV. The average fluxes 0l(u) and ~2(u) in the fuel and moderator respectively are used in the calculations of the second resonance integral. The average first-collision probabilities Pl(u) and P2(u), derived in Chapter IV, due to spatially dependent source distribution are used in these calculations. In the absence of the first resonance integral, the uninterfered second resonance integral will have the same form as shown in Section A, The NR and IM approximations have been used in these calculations, We write the neutron balance Equation (4-2a) assuming that the fuel element is a homogeneous mixture of U, Th and 0 atoms and that the moderator may contain one or more light atoms in the following way:

u V1(u)(4 ++ + ) = V ZP(u) j 1(u?) - du' l-ao' ~ dU0 u- A0 + AVZP1(%U) 1(u') e du' u-AU u + V ZsPj(u) | (u) U - u- ATh UfU du-u + V2(1- P2(u)),j.Z 2(U,) e du u-Aj The average fluxes l1(u') and 42(u') above the second (lower energy) resonance is (u)= < + f e A(U-Uo for u' < u (5-9a) Q 1(u) = All('-u~) 2(U) = 0 + 2 e-All( o), u' < uo (5-9b) ~o' 1' ~ ~2 are defined in Equations (D-2), (D-3), and (E-7) respectively. In order to solve Equation (5-8) we make the usual NR and IM approximations.

1, Narrow Resonance Approximation In the NR approximations we replace 0l(u') and 42(u') in Equation (5-8) by their values in Equations (5-9) and (5-10) when ii = n (1/Qci) is less than or equal to Au = Ur-Uo. Here Ur is the lethargy at the second resonance. And when in(l/oi) > (Ur-Uo), = 02 X= asy in the interval 2n(l/i) - (UrUo) and 41 and 42 are given by Equations (5-9) and (5-10) in the interval (Ur-Uo). Let us evaluate a typical integral. Igu u'-u Ri(u) = / tl(u,) e du' (5-10) u —n( l/oi) For the case Qn(1 /i) < (Ur-Uo) i((i) = ( + ~ [ e(All-l) n(l/ai) 11 e All(u-uo) (-11 RiAu) 00 + f e (5-11) For the case Yn(l/li) > (Ur-Uo) u Ri(u) = 1() eu'-u 01(uf) e du' l-czi u-sn( 1/+ i' ~ ~seu-n(l dui) 11y1-1i "u-Y~n(l /~ )1

-62asy [e (u-uo) - eUn( l/ani)] + [ -(u) 1- 0 I4-e 1 + Le(u-uo) - e-All(U-Uo) (1-ai) (All-l) (5-12) Similarly Rj (u) = + If [e(A1l-1 n(1/oBj) 1 e -(u-u0) (5-13) for the case In(l/c ) < (Ur-Uo) and Rj(u) = asy [e-(Uuo) - e~n(1/j) + o [-(U-Uo)] + -(u-uo) _ e-All(u-uo) +' e2 - ( t-hcj ) (All- 1) (5 -14) for the case ln(l/Uj) > (Ur-Uo)

-63The changes of Ri(u) and Rj(u) are very small in the resonance region compared to the change of cross-sections, So, we take Ri(u) and Rj(u) out of the integral and put Ri(u) = Ri(Ur) and Rj(u) = Rj(Ur) when the integration is over the resonance, Finally, we obtain the interfered resonance integral in the NR approximation for the i-th absorber as (INRint ) (u)a(u)du Res T iF vLsRu (Ur) + ES RTh (Ur) + RO(Ur) F Pl(u)ai(u)du 1 U Th 0 (Ur V Z+Zr V1 Res Rs (l-P2(u) )a(u)du (ZT +Z +Z) Res

2. Infinite Mass Approximation Now we set UU, Th + 1 in Equation (5-8) and use the narrow resonance approximation for the other integrals, i.e., we replace $1(u') and ~2(u') by their values in Equations (5-9a) and (5-9b). Finally we get (IIM)int ()(u)u (Ur) Res 0 F_' P (u)i(u)du (RlO(Ur) ( hI(+Z +ZT )- Pl(u)(Zs +Zs) Res + V2 s Rj(Ur) (1-P2(u))a(u)du(5-16) V1 / ~ (U ) J (MU Th 0 Th U ), j1(Ur) (Zr+Zph+q Z) - Pl(U)(Zs +4s)

-65C, Computational Details and Results The cross sections at the resonance have been calculated using the Breit-Wigner one-level formula, The absorption and the scattering cross sections are aa(E) = ao 4(EE)2/22+1 (5-17) 4(E-Eo)2/F2+ n 1 + (Cn1/2 2(E-EO)/r +s(E) g a, 4(E s 7L +0)2/F'2+1 + ap (5-18) The first term in Equation (5-18) represents the resonance scattering, the second term is due to interference between the resonance and potential scattering, and the third term is the potential scattering cross section of the absorber, rr, rn and r are the radiation, neutron and the total width at the resonance respectively. o = g4-li2 n- (5-19) where Xo is the wavelength associated with a neutron of energy Eo Here Eo is the energy at the peak of the resonance. A simple relation between Xo and Eo is ~o _ 4.55, x 10 ElO / (5-20)

where Eo is in electron volts and Xo is in cm. The statistical weight factor, g, is given by 2J+1 (5-21) 2(2I+1) where I and J are the nuclear and channel spins respectively, For resolved resonances in U238 and Th232 we may take I = 0 and J = 2 But in the computer program (Appendix F) provisions have been made for selecting any value of I and J The integrations over the resonances have been performed numerically using Simpson's rule.(47) The range of integration has been taken six times the practical width of the resonance, No correction has been made for the absorption in the wings of the resonance, The error due to the finite interval of integration is negligible. It has been shown(21) that even when the range of integration is five-times the practical width, the wing correction is less than a few percent of the main part of the absorption. In the calculations of collision probabilities the values of the function E3(x) have been interpolated from the tables given in Reference 44. The integrations in Equations (4-19), (4-20), (4-28) and (4-29) have been evaluated numerically using the Simpson's ruleo The potential scattering cross sections and the resonance parameters used in the calculations are shown in Tables V-l, V-2, and V-3o The practical width, rprac, has been calculated using the relation

-67rprac = r07p (5-22) where all the terms in the above equation have already been defined. In the Tables V-4 through V-10 we have shown the resonance integrals of several pairs of resonances for different lattices. The second (lower energy) resonances have two values for the resonance integral, one assuming no absorption of neutrons during slowing down before reaching the resonance and the other assuming a delta function sink above it. In order to find the rate of absorption of neutrons by the lower energy resonance in the absence of the first (higher energy) resonance, we have to multiply the resonance integral by the density of the absorber atoms and also by the lethargy dependent flux, which is a constant and has the same value as the flux above the first (higher energy) resonance when it is present. Therefore, the absorption rate by the second resonance in the absence of the first resonance is p = ~10 Na(I),unint (5-23) where (10 is the asymptotic flux above the first resonance. The absorption by the second resonance in the presence of the first resonance is P2 = ~2 Na(I)int (5-24)

-68where 02 is the average flux which would be present at the second resonance if the second resonance were not there; (I)int is the interfered resonance integral given by Equation (5-15) or Equation (5-16), Let us define p = pPi (5-25) Pt p/ is the ratio of capture due to interference to non-interference capture. p' has been computed for all pairs of resonances under consideration. In Table V-4 we have shown the interference between the 116.5 ev and 102.5 ev resonances in several U238 - graphite lattices. It is found that the interference is larger in thicker lumps. In Table V-5 we have shown the interference between several selected pairs of resonances of a U238 _ graphite lattice having 2 cm thick fuel lump and 18 cm thick moderator lump. It is found that the interference depends upon the resonance integral of the higher resonance, which is proportional to the strength of the sink, and also upon the lethargy separation, Au, between the first and second resonances. In Table V-6 we have shown the interference between several pairs of resonances in a U02 - heavy water lattice; the thickness of the fuel being 1 cm and that of the moderator 3 cm. In Tables V-7 and V-8 we have shown the interference between pairs of resonances of a Th232 - graphite and a Th232 - heavy water lattice respectively. In Tables V-9 and V-10 we have shown the interference between resonances when the absorber lumps contain a homogeneous

-69mixture of U238 and Th232 having graphite and heavy water as moderator respecticely. From all the tables mentioned above we notice that the interfered resonance integral is larger than the resonance integral with a flat flux approximation in the graphite moderated cells, i.e., defining AI = (I)int - (I)unint (5-26) AI is positive for the graphite moderated lattices. But for the heavy water moderated lattices A I is a negative quantity. But this should not mislead us. The presence of a high energy resonance, of course, decreases the absorption by the second (lower energy) resonance in both the graphite moderated and heavy water moderated lattices. The negative value of p' shows this effect clearly.

-70TABLE V-1 POTENTIAL SCATTERING CROSS SECTIONS ap Elements (barns) D 3.4 c 4.5 0 3.8 Th232 12.0 U238 10.75 TABLE V-2 RESONANCE PARAMETERS FOR U238(In Electron Volts) Eo rn rr rp Approximation (Calculated) 6.7 1i52 x 10-3 2.46 x 10-2 1.080339 IM 21.0 8.9 x 10-3 2.46 x 10-2 1.853616 IM 36.9 3.25 x 10-2 2.46 x 10-2 3.488660 IM 81.3 2.1 x 10-3 2.46 x 10-2 0.408538 NR 90.0 9.0 x 10-5 2.46 x o10-2 0.077299 NR 102.5 6.5 x 10-2 2.46 x 10-2 3.708180 IM 116.5 1.5 x 10-2 2.46 x 10-2 1.110817 NR 189.6 1.35 x 10-1 2.46 x 10-2 5.244164 IM 208.5 5.50 x 10-2 2.46 x 10-2 2.254223 NR 264.5 2.3 x 10-4 2.46 x 10-2 0.072286 NR 274.0 2.7 x 10-2 2.46 x 10-2 1.109286 NR 398.5 1.00 x 10-2 2.46 x 10-2 0.458391 NR 411.0 1.7 x 10-2 2.46 x 10-2 0.645301 NR * Taken from Reference 21.

-71TABLE V-3 RESONANCE PARAMETERS FOR Th232* (In Electron Volts) Eo rn rr rprac Approximation (Calculated) 21.84 2.40 x 10-3 3.0 x 10-2 0.878572 IM 23.48 4.0 x 10-3 3.0 x 10-2 1.120588 IM 113.15 1.10 x 10-2 4.2 x 10-2 1.056990 NR 121.00 1.85 x 10'2 4.1 x 10-2 1.404356 NR 128.50 1.00 x 10-4 4.1 x 10-2 0.083271 NR 129.40 3.10 x 10-3 4.1 x 10-2 0.478585 NR 196.8 1.3 x 10-4 4.10x 10-2 0.076747 NR 199.8 9.0 x 10-3 4.10x 10-2 0.698770 NR 203.00 1.88 x 10-2 4,10x 10-2 1.095739 NR 366.70 3.5 x 10-2 4.10x 10-2 1.255238 NR 371.00 2.4 x 10-2 4.10x 10-2 0.954775 NR 402.8 9.0 x 10-3 4.10x 10'2 0.492138 NR 413.4 2.68 x 10-2 4.10x 10'-2 0.976164 NR 456.4 2.82 x 10'2 4.1 x 10-2 0.962788 NR 465.0 4.5 x 10-2 4.1 x 10-2 1.343243 NR * Taken from Reference 21. TABLE V-4 RESONANCE INTERFERENCE BETWEEN THE 116.5ev AND 102.5ev RESONANCES FOR U238 - GRAPHITE CELLS Thickness of the Absorber Lump = 2a Thickness of the Moderator Lump = 2(b-a) 2a 2(b-a) E0 (I)unint (I)int (I)int-(I)unint x 100 p (cm) (cm) (ev) (barns) (barns) (I)int 116.5 0.136624 1.0 9.0 2.900018 -0.043485 102.5 0.175556 0.180799 116.5 0.117069 2.0 18.0 6.197581 -0.072165 102.5 0.134082 0.142941 116.5 0.117076 2.0 28.0 6.592241 -0.108977 102.5 0.134096 0.143560 116.5 0.106046 4.0 26 9.220431 -0.109864 102.5 0.126165 0.138980 116.5 0.106046 4.0 36.0 11.637581 -0.105127 102.5 0.126164 0.142780 116.5 0.103780 5.0 45.0 18.543231 -0.044865 102.5 0.150339 0.184563

-72TABLE V-5 RESONANCE INTERFERENCE BETWEEN SEVERAL PAIRS OF RESONANCES IN A U238 - GRAPHITE CELL Thickness of the Absorber = 2cm Thickness of the Moderator = 18cm E0 A u (I)unint (I)int (I) nt(I)unint P (ev) (barns) (barns) (int -O 90.0 0.015188 0.101664 0.493352 -0.014496 81.3 0.088780 0.089220 208.5 0.068054 0.095022 4.862141 -0.043569 189.6 0.059805 0.062861 274.0 0.038607 0.035287 1.393043 -0.059245 264.5 4.757314x10-3 4.835669x10-3 411.0 0.01833 0.030886 0.801507 -0.028683 398.5 0.015957 0.016086 TABLE V-6 RESONANCE INTERFERENCE BETWEEN SEVERAL PAIRS OF RESONANCES IN A U02 - HEAVY WATER CELL Thickness of the Absorber = 1.Ocm Thickness of the Moderator = 3.0cm Eo Au (M)unint (I)int (I)int-(I)UnintxlO0 P' (ev) (barns) (barns) (I)int 90.0 0.022161 0.101664 -3.796641 -0.037345 81.3 0.138441 0.133377 208.5 0.102702 0.095022 -4.759987 -o.o48970 189.6 0.112919 0.107788 291.6 0.048700 0.062255 -2.251708 -0.023812 274.0 0.058837 0.057541 411.0 0.028091 0.030886 -1.085112 -o0.012130 398.5 0.024579 0.024315

-73 - TABLE V-7 RESONANCE INTERFERENCE BETWEEN SEVERAL PAIRS OF RESONANCES IN A Th232 - GRAPHITE CELL Thickness of the Absorber = 20cm Thickness of the Moderator = 18.0cm Eo Au ( I) unint (I)int (I)int- (I)unintx100 P (ev) (barns) (barns) (I)int 23.48 0.650028 0.072406 43.312819 -0.383695 21.84 0.568972 1.003704 121.0 0.189870 0.067253 4.763815 -0.155140 113.15 0.174218 0.182933 129.4 0.079632 6.979465x10-3 7.213855 -0.047351 128.5 0.012501 0.013473 203.0 0.087732 0.015889 3.601960 -0.088312 199.8 0.067155 o0.0o69664 413.4 0.033996 0.025975 1o 23086 -0.032289 402.8 0.023351 0.023642 465.0o 0.033140 0.o018668 1.277404 -0.032649 456.4 0.029775 0.030160

-74TABLE V-8 INTERFERENCE BETWEEN RESONANCES IN A Th232 _ HEAVY WATER CELL Thickness of the Absorber Lump = 2.0cm Thickness of the Moderator Lump = 6.0cm Eo Au (I)un(I)int (I)int (I)unintx P (ev) (barns) (barns) () int 23.48 0.630786 0.072406 -2.837988 -0.057473 21.84 0.551959 0.536727 121.0 0.187481 0.067253 -1.723887 -0.025939 113.15 0.171990 0.169076 129.4 0.078597 6.979465x10-3 0.483585 -6.309174 x 10-3 128.5 0.012371 0.012432 203.0 o0.0o86629 0.015889 -0.410150 -0.011506 119.8 0.0o66294 0.066023 371.0 0.038097 0.013569 -0.270201 -6.34224x10-3 366.0 0.043781 0.043663 413.4 0.033575 0.025975 -0.728074 -9.293795x10-3 402.8 0.023053 0.022887 465.0 0.032739 0.018668 -0.522162 -7.695124x10- 3 456.4 0.029407 0.029254

TABLE V-9 RESONANCE INTERFERENCE IN A (U238Th232) - GRAPHITE CELL Thickness of the Absorber = 2cm Thickness of the Moderator = 18cm Number of U238 atoms per unit vol: 0.023915 x 1024 Number of Th232 atoms per unit vol: 0.014650 x 1024 Absorber Eo Au (I)unint (I) int ( ()unintx1.. (ev) (barns) (I)int u238 116.5 0.157977 3.233714 0.029177 Th232 113.15 0.265604 0.274480 -0.133932 U238 208.5 0.090544 0.026733 2.059386 Th232 203.0 0.133127 0.135926 -0.075256 TABLE V-10 RESONANCE INTERFERENCE IN A (U238Th232) - HEAVY WATER CELL Thickness of the Absorber = 2cm Thickness of the Moderator = 18cm Number of U238 atoms per unit vol: 0.023915 x 1024 Number of Th232 atoms per unti vol: 0.014650 x 1024 Absorber Eo Au (I)unint (I)int (I)int-(I)unintx100 p (ev) (barns) (barns4 (I)int u238 116.5 0.157975 9.929177 -0.444601 Th232 113.15 0.265600 0.264425 -0.013552 u238 208.5 0.090542 0.026733 Th232 203.0 0.133125 0.132400 -0.010764 -0.547137

CHAPTER VI CONCLUDING REMARKS Solving the age-diffusion equation in reactor lattices by a Fourier series expansion of the flux, and observing that the series converges rapidly, we have been able to treat several problems of reactor physics satisfactorily. From the investigations made in this thesis we can make the following observations. 1. The equivalent cell approximation used in the calculation of the disadvantage factor and the thermal utilization in diffusion theory is good for most practical lattice cells. Similar observations have been made earlier by several investigators. They used different methods of calculation. This confirms the usefulness of the method of calculation employed here. 2. In the case of a slab lattice, it has been observed that a resonance is influenced by the presence of another resonance above it, i.e., at higher energy. Since there is depletion of flux due to absorption by the higher energy resonance, the absorption of neutronsby the lower energy resonance decreases. The resonance integral does not depend upon the magnitude of the flux if the flux is constant spatially and lethargywiseo But in the case of closely spaced resonances, the resonance integral is different from that calculated by the flat flux approximation. The difference between the resonance integrals is a measure of the resonance interference which arises due to the incomplete flux recovery after the absorption by the higher energy resonance. In the transient region the flux is spatially and lethargy dependent and -76

the flux in the moderator is higher than the flux in the fuel lump. Depending upon the excess flux in the moderator above the lower energy resonance the resonance integral may be greater or smaller than that calculated by the flat flux approximation. The lethargy interval of the transient region depends upon the strength of the higher energy resonance and also upon the size of the fuel lump. We have observed that the ratio of the interfered resonance integral to the flat flux resonance integral is greater than unity in the case of graphite moderated cells but it is less than unity in the case of heavy water moderated cells. The resonance interference in the two- and threedimensional lattices may be treated in the same way. In the calculation of resonance interference, the higher energy resonance has been replaced by a delta function sink. This approximation is expected to be quite good for narrow resonances. The delta function sink approximation may not be satisfactory for a wide resonance. In these calculations, we have assumed two non-overlapping resonances. But there may be cases where two resonances may overlap each other, for example, in the case of a mixture of two absorbers. The 21.84 ev and 23.48 ev resonances of Th232 are wide resonances and also they overlap each other. The results obtained in Table V-7 and V-8 for those cases are not reliable. A Monte Carlo calculation is desirable in order to investigate those overlapping cases0 Also, ZUT-MOD-3 code(51) can be used to treat the overlapping resonances0 The technique of solving the age diffusion equation may also be used to treat lattices with hydrogeneous moderator. In that case we can solve the Selengut-Goertzel equations0

APPENDIX A THE INTEGRAL I(i,j,k,l,m,n) Here we evaluate the integration shown in Equation (2-6) I(i, j,k,l,m,n) = Jo C iT/ 6 o T. 5L- C - O C osj k los f cos cos 7 cos COS os cos --- dxdydz (2-6) b C. Since it is assumed that D is constant in each region of the cell, namely the spherical fuel element (region 1) and the outer moderator lump (region 2), we divide the entire volume of the cell into several parts such that in each part D is constant and integrate over the parts separately. Let us note that the integrand has eight-fold symmetry so that we can integrate over one-eighth of the cell only, (see Figures 1 and la). The integration over the cell can be carried out in the following way: C C OS - COS Cos COS ---------- COS c b C CL Cell p -Vp x,O-x-ye cos l 7k C dxdydz = 8 dy dz I__ Gt s lC + d- x dy dz b27S 0 DEi P S -ZL l t-x-y (A-l) - dx dy dz (A-l) -78

Note that 1__ __x(i-)EX (~m) 7y = - ( COS + Cos )OS 4+ cos (-m)y C)os (ktn) (Z + cos (k (A-2) Let us consider the first term of the above expression after performing the multiplications. The remaining seven terms will give similar results except that the signs of l,m,n will be different in each case. Let us define F(i j,k.l,m,n) as F(iyjk lm n) dx dy dzcos (osj 0 o0 /OP F/- -X cos (ktn)1 TZ 1 dx dy cos(c~x) cos(3y) sin(' 2Xz ) (A-3) where t _; - (j+m)- and -y=(k+n)7i (A-4) Let y =zx sin e dy = / - cos edO x = p sin dx = p cosdC{P (A-5)

80 - We now expand the sines and cosines of Equation (A-3) in terms of its arguments and substitute the values of x and y from Equation (A-5) and get C oc p+q+r 2P R2q 2r 2(P q + r)+3 F(i,j,k,l,m,n) = ZI 2 - ( ) p=o q=o rzo 2 2 2P 2(q+r)+3, 2q 2(r+I) sin o COS p(q ) sin2q Cos Od (A-6) The definite integrals in Equation (A-6) are known. Thus, we get C C I plq+r0r 2P a q 2r 2(p+q+r)t3 F(iy j kl mn) = /,- (2_p)1 (.2q) (2r-l)! P P=o q:o r-o ~2p+l Z-t'~2 2q32 2r(p+q+r+) 2r(+)2) Z-,c Olc 2( p+qr 2P (A-7) P=-o q=o rzo o (,p+qtr+ ) 22(p*t1r)+3 p! q! r we know ( —2 _F( f )IF- ( 2). - (A-8) (2p)!(2q)! (2r+I)! -22(Ptq+r)+ P!q! r!

-81Let us make the following substitutions in Equation (A-7) s =p + q + r t qq+ r and obtain F(i,j,k,l,m,n) =Zi-()f P _ ( (Y) S! 1 (s=o + Z) 22'3 (S-t) (t -)= r! s 3T 2 25+3 72 (/O Z Z + ) (A-9) -o (s 5) 22+3 S! We know(40) that $) Zn+2S Jn(z) Z > n+25 S 2, where n is any number. (A-10) 5-3 2 S! [(n+S+I) Using above relation we can express Equation (A-9) as F(i,j,k,l,m,n) = J3/ (V +) ( )/. (Z2+( )z ); 3/2 2 iJs/ 2 a Bessel function / 2 ~ 2] 2 (A-ll) J3/2 is a Bessel function.

-82The second integral of Equation (A-l) is relatively simple. gdx dy dzcos oT cos J z cos cos cos os a b C b c Co A ijkS il6jmSkn (A-12) A ijk has been defined in Equation (2-5). The third integral of Equation (A-l) is exactly the same as the first one. Taking all the terms of Equation (A-2), we finally get I(i,j,k,l,mn) = (I Di D )F(ij,k,l,m,n) + F(i,j,k,l,m,-n) + F(i,jk,l,-m,n) + F(i,j,k,-l,m,n) + F(i,j,k,-l,-m,n) + F(i,j,k,-l,m,-n) + F(i,j,k,l,-m,-n) + F(i,j,k,-l,-m,-n) Dz abcAijsi lSjmkn (A-13) 2+ Sijk k S~S

APPENDIX B THE INTEGRAL Hz(ij,k,l,m,n) Let us perform an integration containing terms like I D We know that is constant in each region of the cell except at the $3iD fuel-moderator interface where it is a step function. Whereas is a delta function at the fuel-moderator interface but zero everywhere else. Therefore, we can take the average value of at the interface. e XD (D2 - D1)(x 2 y2- z2) (B-l) D = (D2 - D1)S (y - V2 - x2 _ z2) (B-2) ay and D- = (D2 - D () (z o2 -2 y2) (B-3) Let us consider the following integral Hz(i,j,k,l,m,n) =.D cos a cos sin cos ce_/') Zs' at cos o cos n dxdydz (B=4) b C The term inside the bracket can be written as.i...8 =Cos + COS COS + os ) ( sin (k + sin C ) (5) Let us take the first term of Equation (B-5), i.eo, (5+K) g, (j+m)M) - (k -n) Z cos cos sin ct b C and carry out the integration shown in Equation (B-4)o The remaining -83 -

-84seven terms of Equation (B-5) will give similar results except for the change of sign of 1, m, and n. - Cos9 9 cos cos sin4Z dxdydz cell/y -)X-y -IE = 8- ( + (D2 - D) d dy dz cosox cosDy sin (z _ 2 2 _ y2 /0 - / x -xx 2 8.1 + (D2 - DJ) dx dy cosX cosVy siny fP-x-' 0 0 (B-6),f,P ( are defined as in Equation (A-4). Observing that the above integral has exactly the same form as in Equation (A-3), we can easily write i i cos o x cosjy sin= y dxdydz Taking all the terms of Equation (B-5), we finally get kHz(ijgkl,m,n) = 4 (-, )(D2 - D1)[ (k+n) IT F(i,9j,k,l,m,n) + F(i, j,k,-l,m,n) + F(i,j,k,l,-m,n) + F(ij,k,-l,-m,n)} (k-n) ( T F(i,j,k,l,m,-n) + F(i,j,k,l,-m,-n) + F(i, j,k,-l,m,-n) + F(i,j,k,-l,-m,-n) (B-8) Similarly we can derive expressions for Hx(i,j,k,l,m,n) and HyTi/i*k*l1min)

APPENDIX C THE INTEGRAL I(k1,lm,n) The Integration in Equatioh (2-24) is very similar to the integration carried out in Appendix A except that here it is a surface integral. In this case also the entire area of the cell has been divided into several parts in each of which s is constant (see Figure 2) and integrations have been carried out separately. I(k,l,m,n) = L cos kT cos b cos cos dxdy = +4 Di dx + dx {y ax dy { — jO 0 — + dx dy31 (2-20) We know kqx LI)y' nmx ny (mc rk) x cos cos cos Cos b = (cOS = Ot ~ b c) b o (k- m) x ((+ ) y (k-n) Y y + OS G( ) (Cos b + cos (C-1) Let us consider the first term on the right hand side of Equation (C-l) and perform the integrations shown in Equation (2-24) 0 X j dx dy cos o/x cosQy = F(k,l,m,n) where =- (k+rnm)1r and (+n) (c-2) F(kl,m,n) = dx(cos qX) sins p (c-3)

As in Appendix A let us expand the sine and cosine in Equation (C-3) and then make the following substitutions: x = p sin (c-4) dx = p cosGdG cos O(x = (-)P (c-5) p=-o ~~102+ sin/3 - (C-6) Li (2q +)!( Inserting Equation (C-4), (C-5) and (C-6) in Equation (C-3) we get C — cC 2p 2q 2( p+q+) (q F(k,,1,m,n) = - Sirn COS Ode (2 P) (2q+ 2+ 2 0 2 / 3 2? P=-dqo (2P)! (2q)! 2F(P+q+Z) _ (-Z)P4~ O(q ot2 p _ ~(C-7) p=Oqao,~- 2 (P+q+2) _2(Ptq)+l p! q! let s = p + q, and get o< 2S+2 S- ) F(k,l,m,n) =Z 2s+2 r/' 0 (C q) j)) s! s:= 2 [(s+2) 2=o (S-q)! q! s! 2S+1 vZ- u 11/0 o ] P+ / <so 27S s! F(S+2) 2(c+/3)~ Yz

-87J J. 22 (C-8) We used the relation in Equation (A-10) in the above equation. Substituting the values of K( and 3, we get F(k,l,m,n) J+ (km (b~J + k j( ) ( cn-9) In the same way we integrate the second and third term of the right hand side of Equation (2-24) and finally get Di D2 I(k,l,m,n) =, ( s2 - F(kl,m,n) + F(k,l,m,-n) + F(k,l,-m,n) + F(k,l,-m,-n) + C bAk SmT kL (C-10) A/kt has been defined in Equation (2-23)o Following the same method as has been shown in Appendix B we can obtain the expressions for Hx(kl,mn) and Hy(k,l,m,n) easily.

APPENDIX D DERIVATION OF THE FIRST-COLLISION PROBABILITY IN A FUEL LUMP Pl(u) is the probability that the neutrons of lethargy u originating in the fuel lump (-a < x' < a) of region 1 will make their first collision in the same region 1 (probably after crossing other regions such as the moderator in this case) anywhere in the lattice. From Equation (4-16) we write Pl(u) = dx'4(x') E(Zl(x' -x))dx 4(x,u)dx' -a -a -a_ x-a =a a a 00 -a -2nb+a + cLk dx'(x',u) E1( 2n( b-a )Z2-1 +Z1 x-x) du)dx n=l — aa a -2nb+a + LJ dx'4(x',u) E1[2n(b-a) (Z2-Z1) 2 j -a -2nb-a a + El(x'-x) ]dx 4(x',u)dx' (D-l) -a -88

-89The flux, 4(x',u), in the fuel lump is given by Equations (4-9) and (4-10) ao(u) $(x',u) = 2 + al(u)cos cX'x, for -a < x' < a a +c(x',u)dx' = 2al1(u) a where a $L(u) 2a 4(xV,U)dX' -a ao(u) 1 2 + al(u)'a sin(Oaa) 1 (0) AOlS(1) A01So(i) 2 K A11 (u-uo)js (o) - A,1 Jj) + 2 2(u-uo) I- A~013(1) -Ai (u-uO) 2 All - (u-uo)So(l) a1 sin(cOa)eAj1( o) If <+ e -All ( u-u0) e(D-2) Here = 2 - A) _1, uuo{So(O) - soA1)Jjl (D-3) 1 - A( S(1) A(u-uo) s (0) AS~(t)] ~o IS ]] (D-3) Ko2SoA1_ So j (U-UO) I + sin(Z'a) (DI4) 2All C'a

-9o - By definition En(x) = e x n d4 for x > 0 and 1 1 n, a non-negative integer = e/e t tn2 dt (D-5) 0 We also know(44) that En(Ax+B)dx = - A En+l(Ax+b) (D-6) where A and B are constants. Integrating Equation (D-l) with respect to x, we obtain a 14a J dx' 2 + a1(u)cos'x E2 () -E2 (x'+a -a + 1 dx -U ) + ai(u)cos'xJ LE2(O)-E2(Ela-ZlxI a c0 a + 1i dx, 2 + ai(u)cos X'xf E2 {nd2Zn n=l -a + (2n-1)aZ1-E1x' - E2{nd2ZI + (2n+1)aZi-Ex'hj 00 a + 4 a dx' 2 + al(u)cos'x] LE2fnd22 n=1 -a + (2n-1)sZI+71x' - E2{nd2Z + (2n+1)a- l+Zlx'] (D-7

-91where d2 = 2(b-a), the thickness of the moderator. We use the integral form, Equation (D-5), for E2 ~ o f when it is multiplied by cos U'x' and perform the integrations in Equation (D-7) and find Pl(u) = 1 + 4 - E3(2aI) - E3(0 a(u) 00 + 1 ~( {E3(nd2Z2+(n-l)dEZ1) - E3(nd2e+ndl71) Lz 4a4l 2,1 n=l - E3(nd2Z2+ndll1) + E3(nd2e2+(n+l)dlZj)j 00 + 4a1 ao(u) IE3(nd2Z2+(n-1)dlzl - 2E3(nd2Z2+ndlZ1) n=1 + E3(nd2~+(n+1)dllZ)} a - al(u) dx'cos a'x' eat 4a~i -a 0 + |dx'cos CQ'x' e dl(a-x)/t dt -a O 00 a 77 dx'cos x au - end2ae+(2n-1)asl-lx'] /t - a~u- j dx'cos GZx' e dt, 4a4 L n=l -a O dx' cos a'x' e [nd2e+(2n+1)aZL-lx ] /t (D-8) -a 0

-92a + dx'cos'x' e dt -a - dx'cos 0'x' e-nd2 +(2n+)a+x]/t (D-8) -J~~~~~~~a Jn ~~~~~~(Cont'd) here dl = 2a, the thickness of the fuel lump. In Equation (D-8) we have integrations of the following types a 1 [-K l - [K+F x' /t e -K/t cos C'x'e dx' 0 -a [a-K/t ='sin ~l e / ldt ex' /tt2dt 2 2 e ) - CI'sin Ct'a O OI - [K+la] /tt 2t e- [ K+la] /tt + Ocos I'a K Z+'2t2.'+ct' tJ(D-9) 0 0

-93 - and cosZ al'x1/t -K+Elai/tt 2 C~os cx' e LK-j( dt ='sin U'a e t2dt ~1 ~2+ 2t2 -a O 7-LK-Zijai/t - IK-'laI/tt2dt I e tdt E2 122 2 a~j Z'+c'2t2 0 0 r110 - [K+laO /t e tdt (D-10) 21+a'2t2 0 Let us also note that a a.1 dx'cosX el(a+x )/tdt = j dx'cos Ceax' e dt -a 0 -a O C cOsin a. e -2a/ t2dt + t2dt ] Z +o'2t I J+w'2t2 tdt 1 a+1/t2 dt + Zi os ce f a Z1+a)'2t2 O 1(-11~ tt) 0 1 2 t 2

We know that 1 t2dt 1; 1 tanl - (D-12) and 1 tdt 1 loge + _ _ Z+2t2 2,2 loge (D-13) Evaluating all the integrals of Equation (D-8) and rearranging the terms we finally get Pl(u) = 1 - do5 E (n))-2E3(2(n))+E3((n)) n=O ai(u) Ksin Cla' -ZicosaIa ie1 +' d,41 O' 1 ~ge 00 + X [LL1cos alTjE(n) 2- 2 *(n) + T3(n) n=O + Ca'sin C'a 1(n) - 3(n) (D-14) where \Yl(~)1 _i(n)= Te /ttdt n)e ~~)ti = 1,2,3 (D-15)

1 -( i(n)/ttadt ~| aj(n)/(It ji = 1,3 (D-16) \ nd 2+'=t 1 and 11(n) = ndlj1 + (n+1)d2Z2 a2(n) = (n+l)(dlZl+d22) 23(n) = (n+l)dlC 1 + nd2aZ (D-17)

APPENDIX E DERIVATION OF THE FIRST-COLLISION PROBABILITY IN THE MODERATOR LUMP P2(u) is the probability that the neutrons of lethargy u originating in the moderator lump (a < x' < 2b-a) of region 2 will make their first collision in the same region 2 (probably after traversing other regions) anywhere in the lattice. From Equation (4-28) we know 2b-a x' 2b-a P2(u) = dx'4(x') ELZ2 (x'-x)]dx 4(x')dx' S a a b-a 2b-a + 2 dx'E(x') ElZ d(x-x']dx (b-a)f2 a x' 00 g2b-a 2(n+1)b-a 2 dx'(x') ElL2na(Zl- Z)+(x-x' )Z]dx 2(b-a)T2 n=l a 2nb+a 00 2b -a IO rl.2b-a -2(n-l)b-a + ) -- dx'~(x') El 72na(Z1-Z2)+(x'-x)Z]2dx 2(b-a)42 n=l1 a -2nb+a 2b-a.j = b-(x')dx' (E-2) -96

-97From Equations (4-9) and (4-10) we know that 4(x') (u) + a (u)cos + / for x < 2b-a Therefore, we find b$2(U) 2(b-a ) + a1(u)cos 0 cos i'x' a (E-3) sin G sin'x'l dx' where i ~a (1 1 O c = D1, D2 = <~ It ~~~~(E-4) 5D2 Performing the integrations in Equation (E-3) we find P2a(u) = 2 + al(u) sin C (2D-a)-sin I (E-5) a and D have been defined in Equation (2-40). 42(u) can also be written as c2 (U) = 4o~+ 2 e Aj((u-) (E-6) where $o has same value in Equation (D-3) and 2 = - (U-Uo) FAoSo(1) - so( ) sin ( ] (E-7)

-98Integrating Equation (E-1) with respect to x we find 2b 2b-a P2(u) 2E2() 4(x')dx' - dx' ao(u) a a ~a~~~~a + al(u)(cos Q cos B'x'-sin 0 sin P'x')t E2 (2b-a)-exf 00 2b-a + dx' ) +bao(u ))(cos cos'x'-sin 0 sinx)E2(a) a (nd an ao(u)- E (n X+ 2b dxi' 2 + al(u)(cos 0 cos 5'x'-sin 0 sin 5'x'l) +al(u) n=l a ~E2(ndaZi+(n-l)d2-az-2+zx') - E2(ndjZ1+nd2Ze-aeZ+x'), (E-8)

-99As in Appendix D we express E2(. o ) in integral form, Equation (D-5), when it is multiplied by cos i'x' or sin i'x'and get ao(u) P2(U) = 1 - 2d2T2 o.5-E3(d2Z2) - 3(ndlZl+(n-1)d2) n =1 2E3(ndll+nd2Z2) + E3,(ndlZl+(n+l)d2)tl al(u)cos G 2b -a 1 cos 5'x' e dt a1(u)sin 8 2b-a 1 al(u) si n 5'x' e[ i (2b a) ~2x ]/tdt 2d2T2 a 0 i 2b-a.1 al(u)cos Gsi x' eZx-a]/tdt I cos 5'x' a 0 2b-a 1 al(u)sin 2 I b (1 X-Za]/td 0 0 al(u) J [ c e dndFli+nd2Ze+aZe-ex A/t 2d2T2 YI —i a 2b-a /I - -ndll+nd2Z2+aZ -xx'] /t dt + sin @ sin'x' e - (E-9) a O 2b-a 1 + sin 9 sin'x'/ et ndl21+(n+1)d24+a4-Zx']/tdt (E-9) i

-100 - + Cos @ cos'x' P e-1[nddll+1(n -l) d2Ze-aZe+2x't/tdt J a 0 2b-a 1 - sin G sin B'x' endt a 0.-2b-a 1 - cos OJ cos'x' e- [ndlZl+nd2Ze-a4+Zx~ /t dt a 0 -2b-a -1 endiZ,+nd2e_-a~e+ex] /t + sin | sin I'x'J e dt (E-9) J0n -JB /t dt (Cont'd) a O In Equation (E-9) we have integrations of the following types: 2b-a 1 ~~~COS Xf- [K+-X2 /t cos 5'x' e dt a 0 1g -_ [K+Ze (2b -a ) /tt2t e- K+Zeaj/tt2dt ='sin f' (2b-a) e +'t -sin ~ a _ 2+a/2t2 0 0 / a K+a/t K(2b-a)] /t + Zecos l'a - Zcos P'(2b-a) K+(2 /ttdt Z2~_,2t2 "S — 2t2 0 O (E-10) and 2b-a 1 1 -2b-a L [K+-2Xij /t -K/t sin x' e dt = dt e- sin'x' e dx a 0 0 a -K/t" e- Z(2b-a)/tt2 {j. sin 5'(2b-a)-5 cos (2b-a L~+8,2ta _ e t ( t-) sin 2'a-1'cos 1' dt (E-ll) T2~,2+2in~y~~co

-1011 _ K+-/ttdt _ ~esin 5(2b-a) K2-a)Z = Zsin a e - [K+2 A ]tdt - e2sin 0I (2b-a) e tdt 1 e_ [K+ a]/tt2dt + 5'cos h'a / I f ~+5'2t2 - 5'cos 5'(2b-a) i efj+ 2ajt t (E-1l) j Z+ t(Cont'd) Evaluating all the integrals of Equation (E-9) we finally get P2(u) = 1 - ~ ~0. 5 -: E3(' (n))-2E3( 2(n))+E3(i3(n)) n=O a(u)cs loe (loge + (cos 5'(2b-a)+cos h'a) + ( - tan' (sin P'(2b-a)-sin'a) ai(u)sin a Z2 / 2~' a(u)sn 2' loge 1 + 2) (sin 5'(2b-a)+sin h'a) (E-12)

-102+ (1 - tan-' )(cos'(2b-a)+cos 5'a) + Z2{cos'(2b-a)+cos' a} Il(n)-2iI2(n)+n3(n) 2d242 n=O +'[{sin' (2b-a)-sin p'a} iII(n)-zI(n)3 00 a,( ) Z; [ iusin I' (2b-a)+sin 5'a} )}(n)-212(n)+n3(n n=O - %'cos P'(2b-a)-cos I'a II(n)- (n)- (n) (E-12) (Cont'd) where n)e (nl)/tdt 1 (n)1e i = 1,2,3 (E-13) 0' / e- 4(n)/tt2dt (E-14) a(n) = (n+l)dlEl+nd2a l (n) = (n+l)(dlZl+Zed2) Qs(n) = ndla1 + (n+l)d2e (E-15)

APPENDIX F COMPUTER PROGRAMS The computer programs used in obtaining the tables of results given in the text are listed in the following pages. The programs, written in MAD,(42) are for the IBM-7090 computer. They are written mainly for checking the calculations used in this thesis. A brief description of the input and output variables are given below: Program No. 1 This program calculates the average fluxes in the fuel and the moderator lumps, the disadvantage factor and the thermal utilization of a rectangular cell with a cylindrical fuel rod at the center of the cell. Input HKHLHMHN = the highest value of k,l,mn respectively SA = half-width of the cell (in cm) along the x-axis (a) SB = half-width of the cell (in cm) along the y-axis (b) SOR2 = source density in the moderator (S2) SA1 = macroscopic absorption cross section of the fuel element (Cal) for thermal neutrons SA2 = macroscopic absorption cross section of the moderator material (Za2) for thermal neutrons Dl - diffusion coefficient of thermal neutrons in the fuel element (D1) -10o3 -

-104D2 = diffusion coefficient of thermal neutrons in the moderator (D2) DF = an arbitrary constant by which each term of the matrix are divided so that the determinant does not become too large for the machine Output AFLUX1 = average flux in the fuel rod (1) AFLUX2 = average flux in the moderator (F2) RATIO = the disadvantage factor (d) THERU = the thermal utilization (f)

-105PROGRAM NO. 1 $ COMPILE MADEXECUTEPRINT OBJECTPUNCH OBJECT,DUMP DIMENSION A(6650,ADIM),DEL(6650,DIM)*PHI(O100PDIM)I 1 PHIF(100,PDIM),PHIM(100,PDIM) VECTOR VALUES PDIM = 1,1 VECTOR VALUES ADIM s 2,1,0 VECTOR VALUES DIM =4,90,0,00 START READ DATA PRINT RESULTS HM,HNtHK,HLSA,SBSA1*SA2D1,tD2*RODFSOR2 LC = HM+HN+HM*HN+2 LR = (HK+1)*(HL+l) ADIM(2) = LC DIM(l) 1 (HK+1)*(HN+1 )*(HL+1)+(HN+)*(HL+1)+HL+2 DIM(2) = HK+1 DIM(3) =HN+1 DIM(4)= HL+1 PI = 3*142592 RP = PI*R D= D1-D2 S= SA1-SA2 R R CALCULATION OF FIRST AND LAST COLUMN OF THE MATRIX R PRINT COMMENT $0 FIRST AND LAST COLUMN OF THE MATRIX$ A(1.1) = (S*'RP*R/4.0 + SA2*SA*SB)/DF A(1,LC) = (4.0*SA*SB-RP*R)*SOR2/DF PRINT RESULTS A(1,1),A(1LC) I 1 THROUGH LOOP9,FOR L1,1tL.G.HL I = 1+1 Z9 = L/SB X9 = RP*Z9 Y = BSL1.(X9,11,,B9,0) WHENEVER Y.E.2.0,TRANSFER TO COMM A(I,1) = (S*R*B9/(Z9*2.0))/DF A(I,LC) " (-2,0*R*B9/Z9)*SOR2/DF LOOP9 PRINT RESULTS A(I,1),A(ILC) I = HL+1 THROUGH LOOP1,FOR KI1i1. K.G.HK THROUGH LOOP1,FOR L=0-1. L.G*HL I = I+1 Z = SQRT*((K/SA).P*2+ (L/SB)*Po2) X = RP*Z Y = BSL1.(X,I1,,B,O) WHENEVER Y.E.2.0QTRANSFER TO COMM A(I,1) = (S*R*B/(Z*2))/DF A(I.LC) = (-2*R*B/Z)*SOR2/DF LOOP1 PRINT RESULTS A( I 1),A(ILC) I = 0 THROUGH LOOP2, FOR K=O,1, K*G.HK THROUGH LOOP2, FOR L=0.1, L.G.HL I I+1 J = 1 THROUGH LOOP3, FOR N=1,1, N.G.HN

-106J = J+1 WHENEVER N.E.L *AND. K.E.O DEL(OKNL) a SA*SB*2 OTHERWISE DEL(OKtNL) =0 END OF CONDITIONAL Z1 a SQRT(({K/SA)*P*2 + ((N+L)/SB).P*2) X1 =RP*Z1 Y = BSL1.(X191,l1Bl0) WHENEVER Y.E.2.O0TRANSFER TO COMM Z2 = SQRT*((K/SA).P.2 +((N-L)/SB).P.2) X2 = RP*Z2 Y = BSL1*(X2,1,1,9B20) WHENEVER Y.Ee2.0,TRANSFER TO COMM A(ItJ) = (((N*PI/SB)*P.2*(D*R*(B1/Zl+B2/Z2)+D2*DEL(OKNL)) 1 +N*PI/SB*(-D)*R*PI/SB*((N+L)*B1/Zl+(N-L)*B2/Z2)+S*R*(Bl/Z1+ 2 B2/Z2)+SA2*DEL(O0KNL))/2)/DF LOOP3 CONTINUE LOOP2 PRINT RESULTS A(I,2)...A(ItHN+1) I a 0 THROUGH LOOP4, FOR K=O1, K.G.HK THROUGH LOOP4, FOR L=O,1t L.G.HL I I+1 J = HN+1 THROUGH LOOP5sFOR M=u1.1.,iG.HM J a J+1 WHENEVER L.E.O.AND. M.E.K DEL(MK,0OL) =SA*SB*2 OTHERWISE DEL(MKO,L) NO END OF CONDITIONAL Z3 a SQRT*(((M+K)/SA),P.2 + (L/SBI*P.2) X3 a RP*Z3 YV BSL1l(X3,191.B3O0) WHENEVER Y.E.2.O.TRANSFER TO COMM Z4 3 SQRT*(((M-K)/SA)}P.2 + (L/SB )P.2) X4 a RP*Z4 Y a BSL1X(X4ItlslB4*0) WHENEVER Y.E*2.0,TRANSFER TO COMM A(I.J) a I((M*PI/SA)eP.2.0*(D*R'*IB3/Z3+B4/Z4)+D2*DEL(MKOL) 1 )+ (M*PI/SA)*(-D)*R*PI/SA*({M+K)*B3/Z3+(M-K)*B4/Z4)+S*R*(B3/ 2 Z3+84/Z4)+ SA2*DEL{M,K,OL))/2)/DF LOOP5 CONTINUE LOOP4 PRINT RESULTS A(IHN+2)*,.A(I,HN+HM+1) I = 0 THROUGH LOOP69FOR K=01,9 K.G.HK THROUGH LOOP69 FOR L=O,t1 L*G*HL I = 1+1 J 4 HM+HN+1 THROUGH LOOP7,FOR M1,.1, M*G*HM THROUGH LOOP79FOR N1ltl, NeG.HN J 3 J+1 WHENEVER LeEeN *AND. K.E.M DEL(MKtNL) a SA*SB OTHERWISE DEL(M*K*N#L) = 0 END OF CONDITIONAL Z5 = SQRTe(((M+K)/SA).P,2 +((N+L)/SB).P.2) X5 = RP*Z5

-107Y = BSL1(X5,1,1,B5 0) WHENEVER Y.E*2.O0TRANSFER TO COMM Z6 = SQRT.(((M+K)/SA).P.2 +((N-L)/SB)*Po2) X6 = RP*Z6 Y = BSL1*IX6,,lslB60) WHENEVER Y.E.2.o0TRANSFER TO COMM Z7 = SQRT*(((M-K)/SA)IP.2 +((N+L)/SB)*P*2) X7 = RP*Z7 Y = BSL1I(X7,1,iB7.0) WHENEVER Y.E.2.O0TRANSFER TO COMM Z8 = SQRT*(C(M-K)/SA)*P*2 +((N-L)/SB)*P*2} X8 = RP*Z8 Y a BSL1o(X81,1iB8tO) WHENEVER YVE.2e0ATRANSFER TO COMM A(IJ) a (((M*PI/SA).P.2,O +(N*PI/SB).P.2.O)*(D*R/22.*(B5/Z5 1 +B6/Z6+87/Z7+B8/Z8)+D2*DEL(MKNL))+ M*PI/SA*(-D)*PI/SA*R/ 2 2.0*((M+K)*(B5/Z5+86/Z6)+(M-K)*(B7/Z7+B8/Z8)) + N*PI/SB* 3 (-D)*PI/SB*R/2eO*((N+L)*(B5/Z5+B7/Z7)+(N-L)*(B6/Z6+B8/Z8))+ 4 S*R/20*(B5/Z5+B6/Z6+B7/Z7+B8/Z8) + SA2*DEL(M,KINL))/DF LOOP7 CONTINUE LOOP6 PRINT RESULTS A(IIHM+HN+2)...A(I.LC-1) PRINT COMMENT $0 THE MATRIX$ PRINT RESULTS A(1,1).*.A(LRtLC) R R CALL FOR GJRDT SUBROUTINE R PRINT COMMENT $0 CALCULATION OF COEFFIENTSS G a GJRDTI(LR*LC,A(1,1)gDE) WHENEVER G *E. O.09TRANSFER TO OUT PRINT RESULTS DE THROUGH LOOP8. FOR I=olI *G. LR LOOP8 PRINT RESULTS A(I.LC) PRINT COMMENT $0SO AVERAGE FLUX IN FUELS I = 1 PHIF(I) = A(ILC)/4 THROUGH LP1. FOR N=1,1, N.G.HN I I+1 S1 = N/SB T1 = RP*S1 Y * BSL1 (T11II,W190) WHENEVER Y*EE2,0#TRANSFER TO COMM LP1 PHIF(I) = PHIF(I-1)+ A(ILC)/RP*W1/S1 I a HN+l THROUGH LP2s FOR Ml,l, M.GHM I = I+1 S2 = M/SA T2 * RP*S2 Y = BSL1*(T2,11,W2*0) WHENEVER Y.E.2.Q0TRANSFER TO COMM LP2 PHIF(I) = PHIF(I-1) + A(I,LC)/RP*W2/S2 I - HM+HN+1 THROUGH LP3, FOR M=1.1, M*G*HM THROUGH LP3, FOR N=1,1, N*G.HN I = I+1 S3 = SQRT*((M/SA)*.P2 + (N/SB).P*2) T3 = RP*53 Y = BSL1i(T3,1,1,W3,0) WHENEVER Y.E.2.OTRANSFER TO COMM LP3 PHIF(I) = PHIF(I-1) + 2*A(ILCI)/RP*W3/S3

AFLUX1l PHIF(I) PRINT RESULTS AFLUX1 PRINT COMMENT S0 AVERAGE FLUX IN MODERATORS F - 4.0*SA*SB - PI*R*R I 1 PHIM(I) a A(ILC)/4 THROUGH LP49 FOR Nlo,1. NeG.HN I a 1+1 S4 a N/SB T4 a RP*S4 Y a BSL1i(T4,1l,1W490) WHENEVER Y.E.2.0,TRANSFER TO COMM LP4 PHIM(I) = PHIM(I-1) -R*A(I,LC)/F*W4/S4 I = HN+1 THROUGH LP5, FOR M=1,1, M.G*HM I a I+1l 55 = M/SA T5 = RP*S5 Y = BSL1*(T5,itW5S,0) WHENEVER Y.E.2.0*TRANSFER TO COMM LP5 PHIM(I) = PHIM(I-1) -R*A(I,LC)/F*W5/S5 I = HM+HN+i THROUGH LP69 FOR M=1,l1 M.G.HM THROUGH LP6, FOR Ni.l1, N.G.HN I I1+1 S6 = SQRT*((M/SA)*,P2 + (N/SB).P.2) T6 = RP*S6 Y = BSL1,(T6,1.1.W6,0) WHENEVER Y*E*2O.0TRANSFER TO COMM LP6 PHIM(I) = PHIM(I-1) -2,O*R*A(ILC)/F*W6/S6 AFLUX2= PHIM(I) PRINT RESULTS AFLUX2 RATIO a AFLUX2/AFLUX1 INVDIS a= 1+(4*SA*SB-PI*R*R )*SA2/(PI*R*R*SA1 )*RATIO THERU = 1/INVDIS PRINT RESULTS RATIO.THERU INTEGER HM,HNHKtHLLCLR I K L J,.N,M TRANSFER TO START COMM PRINT COMMENT $0 ARGUMENT OF BSL1 TOO LARGES TRANSFER TO START OUT PRINT COMMENT $0 ERROR IN GJRDT OPERATIONS TRANSFER TO START END OF PROGRAM $DATA

-og09Program No. 2 This program calculates the average fluxes in the fuel and the moderator lump, the disadvantage factor, the thermal utilization and flux at different points in the equivalent cylindrical cell. Input A = half-width (in cm) of the rectangular cell whose equivalent cell calculation is need, along the x-axis (a) B = half-width (in cm) of the rectangular cell along the y-axis (b) R1 = radium (in cm) of the fuel rod (p) SAl = macroscopic absorption cross section (Cal) in the fuel for thermal neutrons SA2 = macroscopic absorption cross section (Za2) in the moderator for thermal neutrons D1 = thermal diffusion coefficient (D1) of the fuel material D2 = thermal diffusion coefficient (D2) of the moderator material SEP = intervals of the radius at which the fluxes are calculated Output APHIl = average flux in the fuel rod (1) APHI2 = average flux in the moderator lurp p(2) RATIO = the disadvantage factor (d) THERU = thermal utilization (f)

-110$COMPILE MADEXECUTE,DUMP,PRINT OBJECTI/O DUMP,PUNCH OBJECT START READ DATA PRINT RESULTS ABR1*SA1*SA 2,DDl9D2SEP K1 = SQRT~(SA1/D1) K2 = SQRT*(SA2/D2) R2 = 2*5SQRT*(A*B/3*141592) Zl = K1*Rl Z3 = K.*R1 Z4 = K2*R2 PRINT RESULTS ZltZ39Z4 R R ZERO ORDER MODIFIED BESSEL FUNCTION OF FIRST KIND R L = BSL1*i(Zi,20tBl#O) WHENEVER L *Ee 2*OTRANSFER TO OUT L z BSL1*(Z3.2tOB2.O) WHENEVER L *E. 2.0OTRANSFER TO OUT R R FIRST ORDER MODIFIED BESSEL FUNCTION OF FIRST KIND R L = BSL1.(Z192,1B390) WHENEVER L *E~ 2*0,TRANSFER TO OUT L = BSL1.o(Z3=92.B4,0) WHENEVER L.E. 2.OTRANSFER TO OUT L = BSL1 (Z492,1,B590) WHENEVER L *E* 2*O,TRANSFER TO OUT R R ZERO ORDER MODIFIED BESSEL FUNCTION OF SECOND KIND R L = BSL1*(.Z3.3o0B6,0) WHENEVER L *EE 2*.,TRANSFER TO OUT R R FIRST ORDER MODIFIED BESSEL FUNCTION OF SECOND KIND R L = BSL1e(Z3,31,87.TO) WHENEVER L *E. 2.0.TRANSFER TO OUT L = BSL1.(Z4,3t1,B8,0) WHENEVER L *~E 2.0,TRANSFER TO OUT R R CALCULATION OF Q/(SA2*A) R CONS = Bl-Dl*Kl*B3*.(B2*B8-B6*B5)/( D2*K2*(B4*B8-B7*B5) ) PRINT RESULTS CONS R R CALCULATION OF C/A R CBYA = (81-CONS)/(B2*B8+B6*B5) PRINT RESULTS CBYA R R CALCULATION OF FLUX IN FUEL R PRINT COMMENT $0 FLUX IN THE FUEL$ THROUGH LOOP1,FOR R=0,SEP. R.G. Ri

F = R*K 1 L = BSL1.(F,2,O,F190) WHENEVER L *.E 2.OTRANSFER TO OUT PH1BYA = F1 LOOP1 PRINT RESULTS R,PH1BYA R R CALCULATION OF FLUX IN THE MODERATOR R PRINT COMMENT $0 FLUX IN THE MODERATOR$ THROUGH LOOP2, FOR R =R1,SEP, R *G. R2 M = R*K2 L = BSL1*(M,2*,0M1.0) WHENEVER L *.E 2.0,TRANSFER TO OUT L = BSL1(M,3 O,M2.O ) WHENEVER L *E. 2.09TRANSFER TO OUT PH2BYA =CBYA*(M1*B8+M2*B5)+CONS LOOP2 PRINT RESULTS R,PH2BYA PRINT COMMENT $0 AVERAGE FLUXES$ APHI1 = 2*B3/(K1*R1) APHI2 = 2*CBYA/(K2*(R2*R2-Rl*Rl))*(B8*(R2*B5-Rl*B4) +B5*(R1* 1 B7-R2*B8 ))+CONS PRINT COMMENT $0 DISADVANTAGE FACTOR$ RATIO = APHI2/APHI1 INDIS = 1+ R2*SA2/(Rl*SA1)*RATIO THERU = 1/INDIS PRINT RESULTS APHI,APHI2,RATIO,THERU TRANSFER TO START OUT PRINT COMMENT $0 ERROR IN BSL1 SUBROUTINES TRANSFER TO START END OF PROGRAM $DATA

-112Program No. 3 This program calculates the roots of the determinant, Equation (4-8). Input SA = half-width of the fuel lump (a) SB = half-width of the one-dimensional cell (b) SIGS1 = macroscopic scattering cross section of the fuel material (ZsE) SIGS2 = macroscopic scattering cross section of the moderator (ZS2) MEU1 = the average cosine of the scattering angle per collision in the fuel in the Laboratory Coordinate System (1-1) MEU2 = the average cosine of the scattering angle per collision with the moderator atom in the Laboratory Coordinate System (g2) CIl = the average logarithmic energy decrement of neutrons per collision in the fuel element (1) CI2 = the average logarithmic energy decrement of neutrons per collision in the moderator lump (~2) S1 = source density in the fuel region (S1) HN = the maximum value of n in the summation of the flux. (HN = 4 in this program) Output R(1),R(2), = are the roots of the determinantal equation. R(1) + iR(2) etc. is the first root, but R(2) = R(4) =..o = O, i.e., the imaginary part of the root becomes zero, the roots appear with a negative signo

-113PROGRAM NO. 3 SCOMPILE MADEXECUTE,PRINT OBJECT,DUMPPUNCH OBJECT DIMENSION A(100,ADIM) S(10),8(20) R(20) VECTOR VALUES ADIM = 29,00 READ READ DATA PRINT'OMMENT $0 INPUT DATA FOR NO ABSORPTION CASES PRINT RESULTS SASBSIGS1,SIGS2,MEUl1MEU2tCI1,CI2,SltHN ADIM(1) = HN+2 ADIM(2) = HN+1 Dl = 1*0/(3*SIGS1*(1-MEU1)) D2 = 1.0/(3*SIGS2*(1-MEU2)) W1 = 3e0*(1-MEU1)/CI1 W2 = 3.0*(1-MEU2)/CI2 W = W1-W2 ALP = SA/D1 BETA - ALP+(SB-SA)/D2 PI! 3.1415926 T = PI*ALP/BETA A(O0,) = 0 THROUGH LOOP1,FOR MOl1. M.GiHN A(M90) = 0 THROUGH LOOP1.FOR Nal1l, N.G.HN WHENEVER M.E.N A(M.N) =PI*(N/BETA)e.P2*(W*(SINe((N+M)*T)/(M+N)+T)+W2*PI) OTHERWISE A(M.N) = PI*(N/BETA),P.2*W*SINC((M+N)*T)/(M+N) END OF CONDITIONAL WHENEVER M.E.0 S(M) = 2*S1*T/(PI*CI1*SIGS1) OTHERWISE S(M) = 2*S1*SIN*(M*T)/(PI*CI1*SIGS1*M) END OF CONDITIONAL LOOP1 CONTINUE PRINT RESULTS A(OO).,*A(H'N,HN) PRINT RESULTS S(O)***S(HN) WHENEVER HN*G.2, TRANSFER TO HN3 B(0) = 1 B(1) = 0 B(2) = A(1,1) + A(2.2) 8(3) = 0 8(4) = A(1,1)*A(2,2) - A(2,1)*A(192) B ('5) = 0 Q = ZER2.(2,B(O),R(1)) WHENEVER Q*EE2,Q, TRANSFER TO OUT2 WHENEVER Q*E*30, TRANSFER TO OUT3 PRINT RESULTS R(1).*.R(4) TRANSFER TO READ HN3 WHENEVER HN.G.3tTRANSFER TO HN4 B(O)=1 B(1) = 0 B(2) = A(191)+A(292)+A(3.3) 8(3) = 0 8(4) = A(1,l)*A(292)+A(1,l)*A.(3,3)+A(2,2)*A(3,3)-A(3,2)* 1 A(2,3)-A(1,2)*A(2,1 i-A( 1,3)*A 3,1)

-114B(5) = 0 B(6) = A(1,1)*(A(2,2)*A(3,3)-A (3,2)*A(2,3))-A(1,2)*(A(291)* 1 A(3,3)-A(3l1)*A({23))+A(193)*(A(291)*A(392)-A(292)*A(391)) 8(7) = 0 Q = ZER2*(3,B(O)RR(1)) WHENEVER Q.E*2.09TRANSFER TO OUT2 WHENEVER Q.Eo3.OTRANSFER TO OUT3 PRINT RESULTS R(1).e.R(6) TRANSFER TO READ HN4 WHENEVER HNeGe49 TRANSFER TO HN5 B(0) = 1 B(1) = O0 B(2) = A(l1,) + A(2.2) + A(3,3) + A(4,4) 8(3) = O0 B(4) a A(ll,)*(A(2,2) +A(3.3) +A(494)) +A(2,2)*(A(3,3)+ 1 A(494)) +A(3,3)*A(4,4) -A(293)*A(3,2)- A(2,4)*A(4,2)-A(3,4)* 2 A(493)-A(1,2)*A(21l)-A(193)*A(3,1)-A(194)*A(491) B(5) = 0 B(6) = A(1,1)*(A(2,2)*A(3,3) +A(2,2)*A(4,4)+A(3,3)*A(4,4)1 A(2,3)*A(3,2) -A(2,4)*A(4,2)-A(394)*A(4,#3) +A(2,2)*(A(3,3)* 2 A(4,4)-A(394)*A(4,)3)) -A(23 )*(A(3,2)*A(4,4)-A(3t4)*A(4,2) )+ 3 A(294)*(A(3,2)*A(4,3) -A(4,2)*A(3,3)) -A(192)*(A(291)*A(3,3) 4 +A(2,t)*A(494)-A(2t3)*A(3,1) -A(294)*A(491))+A(1,3)*(A(2,1)* 5 A(392)+A(394)*A(4,1)-A(2,2)*A(391)-A(3,1)*A(494)) -A(194)* 6 (A(2,2)*A(491)-A(3,1)*A(4,3)+A(4.1)*A(393)-A(291)*A(492)) B(7) O0 B(8) =A(l1l)*(A(2,2)*(A(393)*A(4,4)-A(394)*A(4,3)) -A(2,3)*( 1 A(3,2)*A(4,4)-A(3,4)*A(4,2)) +A(2,4)*(A(3,2)*A(4,3)-A(492)* 2 A(3,3))) -A(1,2)*(A(2,1)*(A(3,3)*A(4t4)-A(3,4)*A(4,3)) - 3 A(2,3)*(A(31 )*A(494)-A(3,4)*A(4,l).) +A(2,4)*(A(3,01)*A(4,3)4 A(4,1)*A(3,3))) +A(1,3)*tA(2,1)*(A(32)*A(4,4)-A(3,'4)*A(4,2) 5 ) -A(2,2)*(A(3,1)*A(4,4)-A(3,4)*A(4,1)) +A(2,4)*(A(391)* 6 A(4t2)-A(3,2)*A(4,1))) -A(1,4)*(A(2,1)*(A(392)*A(4,3)-A(4,2) 7 *A(3,3)) -A(292)*(A(391)*A(4,3)-A(491)*A(3,3)) +A(293)*( 8 A(3,1)*A(4,1)-A(3,2)*A(4,1))) B(9) = 0 Q = ZER2.(4B,(O)9R(1)) WHENEVgR Q*.E.2.0 TRANSFER TO OUT2 WHENEVER Q.E.3.0, TRANSFER TO OUT3 PRINT RESULTS R(1).*.R(8) HN5 TRANSFER TO READ OUT2 PRINT COMMENT SO ARGUMENTS OF ZER2 ARE OUT OF RANGES TRANSFER TO READ OUT3 PRINT COMMENT $0O IMPOSSIBLE TO LOCATE THE ROOTS$ TRANSFER TO READ INTEGER MN,HN END OF PROGRAM $DATA

Program No. 4 This program computes the resonance interference between a pair of resonances in a slab lattice. The fuel lump contains U238, Th232 and 016 but the moderator may contain any monoatomic moderator. Input SA = half-width of the fuel slab (in cm) (a) SB = half-width of the cell (in cm) (b) NU = number of U238 atoms in fuel/barn-cm NTH = number of Th232 atoms in fuel/barn-cm NOX = number of oxygen atoms/barn-cm SIG2 = total macroscopic cross section of the moderator ( 2) MEU2 = the average cosine of the scattering angle per collision with the moderator atom in the Laboratory Coordinate System (~iM) CI2 = the average energy decreasement per collision (~2) SORS = source density at higher energy SPINJ = channel spin (J) SPINI = nuclear spin (I) ALPM = (JA )2 AM+l EFR = energy of the first (higher energy) resonance (Eo) GNF = the neutron width (in ev) at the first resonance (rn) GGF = the radiation width (in ev) at the first resonance ( Fr) 0, NR approximation for the first resonance IMF =, IM approximation for the first resonance

-116{0, when the first resonance is due to Th232 SURAN =, when the first resonance is due to U238 HMIHN = the highest value of the m,n in the expansion of the flux, here HM = HN = 2 will do NFR = one half the number of divisions used in integrating over the first resonance (the range of integration is six times the practical width) ESR = energy at the peak of the second (lower energy) resonance (0, when the second resonance is due to Th232 SURAN = 1, when the second resonance is due to U238 0, when NR approximation is used for the second resonance IMS =, when IM approximation is used for the second resonance GNS = the neutron width (in ev) at the second resonance ('n) GGS = the radiation width (in ev) at the second resonance (Vr) NSR = one half the number of divisions used in integrating over the second resonance (range of integration is six times the practical width) Output GAMP = the practical width (in ev) of the first resonance (pra prac FRI = resonance integral of the first resonance (I) GAMPS = the practical width (in ev) of the second resonance (Frac) UNIRI = the uninterfered second resonance integral, (I) unint DELU = the lethargy difference between the first and second resonance, Au = ur-uo

-117RI = the interfered second resonance integral (I)int PCERR = (I)int - (I)unint x 100 (I)int P: - Pl PCINA = P1

-118 - PROGRAM NO. 4 $COMPILE MAD, EXECUTE, PRINT OBJECT, DUMP, PUNCH OBJECT DIMENSION A(81,ADIM),S(10) AS(10),XT(282).YT(282),F(250), 1 FN(250),IN(100),JN(100) EXECUTE FTRAP. VECTOR VALUES ADIM = 2,0,0 READ DATA SIGPU = 10.75 SIGPTH = 12.0 SIGPO = 3.8 MEUU = 0.0028 MEUTH = 0.0029 MEUO=0.0317 CIU = 0.0084 CITH = 0.0086 CIOX = 0.12 START READ DATA PRINT RESULTS SA, SB, NU, NTH, NOX, SIG2,MEU2.CI2.SORS, 2 EFR, NFR, SPINJ, SPINI, GNF, GGF, IMF, FURAN, HM,HN, 3 ESR, NSR, GNS, GGS,ALPM, SURAN,IMS ADIM(1) =HM+2 ADIM(2) = HM+1 SIGSU = NU*SIGPU SIGSTH = NTH*SIGPTH SIGSO = NOX*SIGPO SIGS1 = SIGSU+SIGSTH+SIGSO MEU1 =(MEUU*SIGSU + MEUTH*SIGSTH + MEUO*SIGSO)/SIGS1 D1 = 1*0/(3*SIGS1*(1-MEU1)) D2 = 1.0/(3*SIG2*(1-MEU2)) CI1 = (CIU*SIGSU+CITH*SIGSTH+CIOX*SIGSO)/SIGS1 W1 = 3*0*(1-MEU1)/CIl W2 = 3.0*(I-MEU2)/CI2 W = W1-W2 ALPH = SA/D1 BETA = ALPH + (SB-SA)/D2 PI = 3.1415926 PIT = PI*CI2*5IG2 PIC = PI*CI1*SIGS1 T = PI*ALPH/BETA TM = 2*(SB-SA) TF = 2*SA THROUGH LOOP, FOR M=0.1, M*G.HM A(M.0) =0 WHENEVER ME.0C S(M) = 2*SORS*T*(1/PIC-1/PIT) + 2*SORS/(CI2*SIG2) OTHERWISE S(M) = 2*SORS*SIN.(M*T)/M*(1/PIC-1/PIT) END OF CONDITIONAL THROUGH LOOP, FOR N=ll, N.G.HM WHENEVER M.E.N A(M,N) = PI*(N/BETA).P.2*(W*(SIN.((M+N)*T)/(M+N)+T)+W2*PI) OTHERWISE A(M.N) = PI*(N/BETA).P.2*W*SIN.((M+N)*T)/(M+N) END OF CONDITIONAL

LOOP CONTINUE PRINT RESULTS A(OO0),..A(HM,HM) PRINT RESULTS S(0)...S(HM) R CALCULATION OF THE STRENGTH OF THE SINK R EFR IS THE ENERGY AT THE PEAK OF THE FIRST RESONANCE PHIASH = SORS*SB/(SA*CIl*SIGSi + (SB-SA)*CI2*5IG2) PRINT RESULTS PHIASH G = (2*SPINJ+I)/(2*(2*SPINI+i)) GTF = GNF+GGF LAMH = (4.55E-10)/SQRT.(EFR)*(.1OE12) C1 = 4*PI*L*LA*LAMH*GNF/GTF*G C2 = Cl*GGF/GTF*SQRT.(EFR) C3 = Cl*GNF/GTF PRINT RESULTS C1,C29C3 WHENEVER FURAN *E.l GAMP = GTF*SQRT.(C1/SIGPU) PRINT COMMENT $OFIRST RESONANCE ATEFR IS DUE TO URANIUM$ OTHERWISE PRINT COMMENT $OFIRST RESONANCE AT EFR IS DUE TO THORIUM$ GAMP = GTF*SQRT*(C1/SIGPTH) END OF CONDITIONAL FEL = 3.0*GAMP PRINT RESULTS GAMP,FEL HE = EFR+FEL LE = EFR-FEL H =FEL/NFR R R2FEL IS ENERGY INTERVAL OF INTEGRATION OF FIRST RESONANCE R 2NFR IS NUMBER OF DIVISION USED IN INTEGRATING FIRST RESON WHENEVER IMF.Eo1 PRINT COMMENT $0 INFINITE MASS APPROX FOR FIRST RESONANCE$ OTHERWISE PRINT COMMENT $0 NARROW RESONANCE APPROX FOR FIRST RESONANCES END OF CONDITIONAL THROUGH LOO 0 THROUGH LOOP19 FOR E=LEH, EoGoHE X = 2*(E-EFR)/GTF RT = 1+X*X SIGA - C2/(SQRT,(E)*RT) SIGSR = C3/RT WHENEVER FURAN *E1i SIGSI = SQRT.(C3*G*SIGPU)*X/RT SIGT = NU*(SIGA+ SIGSR+ SIGSI+ SIGPU) SIGT1 = SIGT + SIGSTH + SIGSO TRANSFER TO KHAL OTHERWISE SIGSI = SQRT.(C3*G*SIGPTH)*X/RT SIGT = NTH*(SIGA+ SIGSR+ SIGSI+ SIGPTH) SIGTl = SIGT+ SIGSU+ SIGSO END OF CONDITIONAL R CALCULATION OF FIRST COLLISION PROBABILITY IN FUEL P1 KHAL Y =0 THROUGH EDA. FOR N=091, N.G.HN L1 = (N+i)*TM*SIG2 + N*TF*SIGT1 L2 = (N+1)*(TM*SIG2 + TF*SIGT1) L3 = N*TM*SIG2 + (N+i)*TF*SIGT1 WHENEVER LlGo100 Y1 = 0 OTHERWISE

-120Y1 = TAB.(LlXTYT9,1,1,3, 281,SW) END OF CONDITIONAL WHENEVER SW.E.209,TRANSFER TO OUT1 WHENEVER L2.G.100.O Y2 = 0 OTHERWISE Y2 = TAB.(L2,XTgYT9191,3,281,SW) END OF CONDITIONAL WHENEVER SW.E2.09,TRANSFER TO OUT1 WHENEVER L3.G.10.0 Y3 = 0 OTHERWISE Y3 = TAB.(L3,XT,YT,19193t281,SW) END OF CONDITIONAL WHENEVER SW.E.2.OTRANSFER TO OUT1 EDA Y = Y+Y1+Y3-2*Y2 P1 = 1-(1-2*Y)/(2*SIGT1*TF) WHENEVER IMF*E.1 F(J) = (SIGT1-(SIGT1-SIGSO)*P1)*SIGA/(E*(SIGT1-Pl*(SIGSU+ 1 SIGSTH))) J = J+1 OTHERWISE F(J) = SIGA/E - (SIGT1-SIGS1)*P1/(SIGT1*E)*SIGA J = J+l END OF CONDITIONAL LOOP 1 CONTINUE F1 = 0 THROUGH JAHAN, FOR N=1,1, N.G.(J-1)/2 JAHAN F1 = Fl + F(2*N-1) F2 = 0 THROUGH ARA9 FOR N=1ll, N.G.(J-3)/2 ARA F2 = F2 + F(2*N) R CALCULATION OF THE STRENGTH OF THE SINK WHENEVER FURAN.E.1 NF = NU OTHERWISE NF = NTH END OF CONDITIONAL F.RI = H/3*(F(O) + 4*F1 + 2*F2 + F(J-1)) PRINT RESULTS FRI PRINT COMMENT $OSTRENGTH OF SINK AT HIGHER ENERGY RESONANCE$ SINK m PHIASH*NF*FRI PRINT RESULTS SINK THROUGH RAJA, FOR M=01, M.GoHM WHENEVER M.E.O AS(M)= 2*SINK*T*(1/PIC) OTHERWISE AS(M)= 2*SINK*SIN.(M*T)/M*(1/PIC) END OF CONDITIONAL RAJA CONTINUE PRINT RESULTS AS(O).*.AS(HM) R R CALCULATION OF UNINTERFERRED SECOND RESONANCE R GTS = GNS + GGS LAML = (4.55E-10)*(*1OE+12)/SQRT. (ESR) CS1 = 4*PI*LAML*LAML*G*GN.S/GTS CS2 = CS1*SQRT.(ESR)*GGS/GTS R

-121R ESR IS THE ENERGY AT THE PEAK OF SECOND RESONANCE R CS3 = CS1*GNS/GTS PRINT RESULTS CS1,CS2,CS3 WHENEVER SURAN.Ee 1 PRINT COMMENT $0 SECOND RESONANCE AT ENERGY ESR IS URANIUM$ GAMPS = GTS*SQRT.(CS1/SIGPU) OTHERWISE GAMPS = GTS*SQRT.(CS1/SIGPTH) PRINT COMMENT $0 SECOND RESONANCE IS DUE TO THORIUM$ END OF CONDITIONAL SEL = 3.0*GAMPS PRINT RESULTS GAMPSSEL HE = ESR + SEL LE = ESR-SEL H = SEL/NSR R R 2SEL IS THE ENERGY INTERVAL FOR THE INTEGRATION OF SECOND R RESONANCE R 2NSR IS THE NUMBER OF DIVISIONS USED IN INTEGRATING THE R SECOND RESONANCE R WHENEVER IMS.E*1 PRINT COMMENT $OINFINITE MASS APPROX* FOR SECOND RESONANCE$ OTHERWISE PRINT COMMENT $ONARROW RESONANCE APPROX FOR SECOND RESONANCE$ END OF CONDITIONAL J = 0 THROUGH LOOP2. FOR E=LEH, E.G.HE X = 2*(E-ESR)/GTS RT = 1+X*X SIGA = CS2/(SQRT.(E)*RT) SIGSR = CS3/RT WHENEVER SURAN E 1 SIGSI = SQRT.(CS3*G*SIGPU)*X/RT SIGT = NU*(SIGA + SIGSR + SIGSI + SIGPU) SIGT1 = SIGT + SIGSTH + SIGSO TRANSFER TO PRi OTHERWISE SIGSI = SQRT.(CS3*G*SIGPTH)*X/RT SIGT = NTH*(SIGA + SIGSR + SIGSI + SIGPTH) SIGT1 = SIGT + SIGSU + SIGSO END OF CONDITIONAL R R CALCULATION OF COLLISION PROBABILITY PR1 Y = 0 THROUGH YURA, FOR N=0,1, N.G.HN L1 = (N+1)*TM*SIG2 + N*TF*SIGT1 L2 = (N+1)*(TM*SIG2 + TF*SIGT1) L3 = N*TM*SIG2 + (N+1)*TF*SIGT1 WHENEVER L1G.10O0 Y1 = 0 OTHERWISE Y1 = TAB.(L1, XTtYTl,31*39281,SW) END OF CONDITIONAL WHENEVER SW.E.2.0, TRANSFER TO OUT1 WHENEVER L2.G.10.0 Y2 = 0 OTHERWISE

-122Y2 = TAB.(L2, XTgYT191,,3,281,SW) END OF CONDITIONAL WHENEVER SWE.2o0, TRANSFER TO OUT1 WHENEVER L3.G,10.0 Y3 = 0 OTHERWISE Y3 = TAB.(L3, XT.YT,11i,3281,SW) END OF CONDITIONAL WHENEVER SW*E.2'.0, TRANSFER TO OUT1 YURA Y = Y+Y1+Y3-2*Y2 P1 = 1-(1-2*Y)/(2*SIGT1*TF) PRINT RESULTS EY',P1 WHENEVER IMS.E.1 F(J) = (SIGT1-(SIGT1-SIGSO)*Pl)*SIGA/(E*(SIGT1-Pl*(SIGSU+ 1 SIGSTH))) J = J+1 OTHERWISE F(J) = SIGA/E-(SIGT1-SIGS1)*Pl/(SIGT1*E)*SIGA J = J+1 END OF CONDITIONAL LOOP2 CONTINUE F1 = 0 THROUGH PIEM, FOR N=19ll N.G.(J-1)/2 PIEM F1 = Fl+F(2*N-1) F2 = 0 THROUGH SUWAN9 FOR N='11, N*G.(J-3)/2 SUWAN F2 = F2+F(2*N) PRINT COMMENT $0 UNINTERFERRED SECOND RESONANCE INTEGRAL$ UNIRI = H/3*(F(0)+4*F1 + 2*F2 + F(J-1)) PRINT RESULTS UNIRI R R INTERFERRED SECOND RESONANCE INTEGRAL R ES IS THE ENERGY OF SOURCE NEUTRONS R ALPOX = 15.0*15,0/(17.0*17*0) ALPU = 237.0*237.0/(239.0*239.0) ALPTH = 231.0*231,0/(233,0*23300) ES = 2.00E6 US = ELOG.(ES/ESR) UF = ELOG.(ES/EFR) DELU = US-UF PRINT RESULTS USUFDELU WHENEVER US.G.UF STEP: 1 OTHERWISE STEP = 0 PRINT COMMENT $0OFIRST RESONANCE DO NOT HAVE HIGHER ENERGY$ END OF CONDITIONAL CB = STEP*(AS(0)-A(O,1)*AS(1)/A(1,1))/2 CD = STEP*(S(O)-A(0,1)*S(1)/A(1,1))/2 PHIZ = CD-CB PHIPI = -STEP*AS(1)*(A(091)/(2.0*A(1,1))+ SIN.(T)/T) PHIP2 =-STEP*AS(1)*(A(0,1)/(2.0*A(1,1)) - BETA/(PI*(BETA1 ALPH))*SIN.(T)) AZ = 2.0*PHIZ-STEP*A(Ol)*AS(1)/A(1,1)*EX EX = EXP.(-A(1,1)*DELU) Al = -STEP*AS(1)*EX PHIBi = PHIZ+PHIP1*EX PHIB2 = PHIZ+PHIP2*EX

-123CON1 = AZ/(2*TM*PHIB2) CON2 = A1/(2*TM*PHIB2)*COS*(T-T*D1/D2) CON3 = A1/(2*TM*PHIB2)*SIN*(T-T*D1/D2) SBP = SIN.(BETAP*(2*SB-SA)) CBP = COS.(BETAP*(2*SB-SA)) CBA = COS.(BETAP*SA) SBA = SIN. ( BETAP*SA) ALPHP = PI/(BETA*D1) BETAP = PI/(BETA*D2) J = 0 THROUGH LOOP3, FOR E=LE,H, E*G.HE X = 2*(E-ESR)/GTS RT = 1+X*X SIGA = CS2/(SQRT.(E)*RT) SIGSR = CS3/RT WHENEVER SURAN*.E1 SIGSI =SQR T. (CS3*G*S I GPU)*X/RT SIGT = NU*(SIGA+SIGSR+SIGSI+SIGPU) SIGT1 = SIGT+SIGSTH+SIGSO TRANSFER TO PR2 OTHERWISE SIGSI = SQRT.(CS3*G*SIGPTH)*X/RT SIGT = NTH*(SIGA+SIGSR+SIGSI+SIGPTH) SIGT1 = SIGT+SIGSU+SIGSO END OF CONDITIONAL R R CALCULATION OF COLLISION PROBABILITIES R PR2 PSI = 0 PSIP = 0 CPI = 0 CPIP = 0 YM = 0 YF = 0 THROUGH KENJI, FOR N=0,1, N.G.HN L1 = (N+1)*TM*SIG2 + N*TF*SIGT1 L2 = (N+L)*(TM*SIG2+TF*SIGT1) L3 = N*TM*SIG2+(N+1)*TF*SIGT1 LP1 = N*TM*SIG2 + (N+1)*TF*SIGT1 LP2 = (N+1)*(TM*SIG2+TF*SIGT1) LP3 = (N+1)*TM*SIG2 + N*TF*SIGT1 WHENEVER L1 *.G 10.0 YF1 = 0 OTHERWISE YF1 = TAB.(L1,XT.YT,11,3q281,SW) END OF CONDITIONAL WHENEVER SWE.2.0, TRANSFER TO OUT1 WHENEVER L2 *.G 10.0 YF2 = 0 OTHERWISE YF2 = TAB.(L2,XT,YT,1l1l3v281,SW) END OF CONDITIONAL WHENEVER SW.E.2.0. TRANSFER TO OUT1 WHENEVER L3 *G. 10.0 YF3 = 0 OTHERWISE YF3 = TAB.(L3,XT,YT,11,3tq281,SW) END OF CONDITIONAL WHENEVER SW.E.2.O0, TRANSFER TO OUT1

-124YF = YF+YF1+YF3-2*YF2 WHENEVER LP1.G. 10.0 YM1 = 0 OTHERWISE YM1 =TAB.(LP1,XTYT,1,1,3, 281,SW) END OF CONDITIONAL WHENEVER SW.E.2.o0 TRANSFER TO OUT1 WHENEVER LP2.G. 10.0 YM2 = 0 OTHERWISE YM2 =TAB. (LP2,XT,YT,i1,13,281,SW) END OF CONDITIONAL WHENEVER SW.E.2.09 TRANSFER TO OUT1 WHENEVER LP3.G. 10.0 YM3 = 0 OTHERWISE YM3 =TAB. ('LP3,XT,YT 1,1,3,281,SW) END OF CONDITIONAL WHENEVER SW.E.2.0, TRANSFER TO OUT1 YM = YM + YM1 + YM3 - 2*YM2 R R CALCULATIONS OF THE INTEGRALS PSI AND PSIPRIME T1 = L1/15.0 WHENEVER T1iLe1.0.AND. T1.G.O.5 H1 = (1*0-T1)/20.0 OR WHENEVER T1.LE.O.5 *AND. T1,G*.00000i1 H1 = (1i0-T1)/40.0 OTHERWISE PSI1 = 0 PSIPI= 0 TRANSFER TO SI2 END OF CONDITIONAL K = 0 THROUGH RUSTAM, FOR ST' T1,H1, ST.G*1.0 IN(K) = ST*EXP.(-Ll/ST)/(SIGTi*SIGTi+(ALPHP*ST)oP2O) JN(K) = IN(K)*ST RUSTAM K = K+I JJ = 0 II = 0 THROUGH ALI, FOR M=li,. M.G.(K-1)/2 II = II+IN(2*M-1) ALI JJ = JJ+JN(2*M-1) IK = 0 JK =0 THROUGH GOLAM, FOR M=-'I,l M.G.(K-3)/2 IK = IK+ IN(2*M.) GOLAM JK = JK + JN(2*M) PSI1 = H1/3.0*(IN(O)+4*II+2*IK+IN(K-1)) PSIP1 = H1/3.0*(JN(O)+4*JJ+2*JK+JN(K-1)) SI2 T2 = L2/15.0 WHENEVER T2.L.i.0 *AND. T2.G.O.5 H2 = (1.0-T2)/20*0 OR WHENEVER T2.LE0..5 *AND. T2G.OOO000001 H2 = (1.0-T2)/40.0 OTHERWISE PSI2 = 0 TRANSFER TO SI3 END OF CONDITIONAL K =O

-125THROUGH FAROOQ, FOR ST =T2,H2, ST.G.1.0 IN(K) = ST*EXP.(-L2/ST)/(SIGT1*SIGT1+(ALPHP*ST).P.2.O FAROOQ K = K+1 II = 0 JJ = 0 THROUGH ANO, FOR M=1l,1, MoG.(K-3)/2 ANO II = II+IN(2*M-1) IK = 0 JK =0 THROUGH WARA, FOR M=1,1, M.G.(K-3)/2 WARA IK = IK+ IN(2*M) PSI2 = H2/3.0*(IN(O)+4*II+2*IK+IN(K-1)) SI3 T3 = L3/15.0 WHENEVER T3*L.1*0 *AND. T3.G.0.5 H3 = (1.0-T3)/20.0 OR WHENEVER T3.LE.O05 *AND. T3*G.0.000001 H3 = (1.0-T3)/40'0O OTHERWISE PSI3 = 0 PSIP3= 0 TRANSFER TO NAHID END OF CONDITIONAL K = 0 THROUGH BEGUM, FOR ST:T3,H3,ST.G-1.0 IN(K) = ST*EXPo(-L3/ST)/(SIGT1*SIGT1+(ALPHP*ST)Po2.O) JN(K) = IN(K)*ST BEGUM K = K+1 II = 0 JJ = 0 THROUGH RABEYA, FOR M=1,1, M.G.(K-1)/2 II = II + IN(2*M-1) RABEYA JJ = JJ+JN(2*M-1) IK = 0 JK =0 THROUGH MALEKA, FOR M=1,1 M.G.(K-3)/2 IK = IK+ IN(2*M) MALEKA JK = JK + JN(2*M) PSI3 = H3/3.0*(IN(O)+4*II+2*IK+IN(K-1)) PSIP3 = H3/3.0*(JN(0)+4*JJ+2*JK+JN(K-1)) NAHID PSI = PSI + PSI1 + PSI3 - 2*PSI2 PSIP = PSIP+PSIP1-PSIP3 R R CALCULATIONS OF INTEGRALS CAP-PI AND CAP-PI PRIME R TP1 = LP1/15.0 WHENEVER TP1.L1lO0 *AND. TP1.G.O.5 HP1 = (10.-TP1)/20.0 OR WHENEVER TP1.LE.O.5.AND. TP1.G.O.OO001 HP1 = (1*0-TP1)/40.0 OTHERWISE CPI1 = 0 CPIP1 = 0 TRANSFER TO PI2 END OF CONDITIONAL K = 0 THROUGH MATIOR. FOR ST=TP1,HP1, ST.G.1.O IN(K) = ST*EXP.(-LP1/ST)/(SIG2*SIG2+(BETAP*ST)*P.2) JN(K) = IN(K)*ST MATIOR K = K+1

-126THROUGH RAHMAN. FOR M=1,1 M.G.(K-1)/2 II = II + IN(2*M-1) RAHMAN JJ = JJ + JN(2*M-1) IK = 0 JK =0 THROUGH KHALEQ, FOR M=l,1i M.G.(K-3)/2 IK - IK+ IN(2*M) KHALEQ JK = JK + JN(2*M) CPIi = HP1/3.0*(IN(O)+4*II+2*IK+IN(K-1)) CPIP1 = HP1/3.0*(JN(O)+4*JJ+2*JK+JN(K-1)) PI2 TP2 = LP2/15.0 WHENEVER TP2.L.1.0 *AND* TP2.G.O5 HP2 = (1*0-TP2)/20,0 OR WHENEVER TP2.LE.05 *AND. TP2*GO00001 HP2 = (1O-TP2)/40*0 OTHERWISE CPI2 = 0 TRANSFER TO PI3 END OF CONDITIONAL K = 0 THROUGH KERR, FOR ST=TP2.,HP2, SToG.1,O IN(K) = ST*EXP. (-LP2/ST ) / (SIG2*SI G2+( BETAP*ST ) *P2) KERR K = K+1 II - 0 THROUGH BOHR, FOR M:1,1, M.G.(K-1)/2 BOHR II = II + IN(2*M-1) IK = 0 THROUGH PAUL, FOR M=11, M.G.(K-3)/2 PAUL IK = IK+ IN(2*M) CPI2 = HP2/3.0*(IN(O)+4*II+2*IK+IN(K-1)) PI3 TP3 = LP3/15.0 WHENEVER TP3.L.1*0 *AND. TP3G.O0*5 HP3 = (1*0-TP3)/20,0 OR WHENEVER TP3.LE*O0*5 *AND. TP3.oG.0*0001 HP3 = (1*0-TP3)/40,0 OTHERWISE CPI3 - 0 CPIP3 0= TRANSFER TO YIP END OF CONDITIONAL K = 0 THROUGH ZWE, FOR ST=TP3,HP3, ST*G*1*0 IN(K) = ST*EXP. (-LP3/ST)/(SIG2*SIG2+(BETAP*ST ) P-2) JN(K) = ST*IN(K) ZWE K = K+1 II = 0 JJ = 0 THROUGH IFEL, FOR M=1,1, MG ( K-1)/2 II = II + IN(2*M-1) IFEL JJ = JJ+JN(2*M-1) IK = 0 JK =0 THROUGH OSBORN,FOR M=l,l, M-G.(K-3)/2 IK = IK+ IN(2*M) OSBORN JK = JK + JN(2*M) CPI3 = HP3/3.0*(IN(O)+4*II+2*IK+IN(K- 1)) CPIP3 = HP3/3.0*('JN(O)+4*JJ+2*JK+JN(K-1)) YIP CPI = CPI+CPI1+CPI3- 2*CPI2 KENJI CPIP = CPIP+CPIP1-CPIP3

-127P1 = 1-(AZ/SIGT1*(0.5-YF)/2+Al*(SIN.(T)/ALPHP*(l-SIGT1/ALPHP* 1 ATAN.(ALPHP/SIGT1))+SIGT1*COS.(T)*ELOG.(l+(ALPHP/SIGT1).P.2) 2 /(2*(ALPHP).P*2)) -Al*(SIGT1*COS.(T)*PSI+ALPHP*SIN.(T)*PSIP) 3)/(TF*PHIB1) P2 = 1-CON1/SIG2 *(0O5-YM )-CON2*(SIG2/(2*BEETAP*BETAP)* 1 ELOG.(l+(BETAP/SIG2)oPo2)*(CBP+CBA) + (1-SIG2/BETAP*ATAN. 2 (BETAP/SIG2))*(SBP-SBA)/BETAP-SIG2*(CBP+CBA)*CPI-BETAP* 3 (SBP-SBA)*CPIP) + CON3*(5IG2/(2*BETAP*BETAP)*ELOGo(1+(BETAP 4 /SIG2).P.2)*(SBP+SBA) + (CBP+CBA)*(1-SIG2/BETAP*ATAN.(BETAP 5 /SIG2))/BETAP - SIG2*(SBP+SBA)*CPI+BETAP*(CBP-CBA)*CPIP) PRINT RESULTS E, P1, P2 WHENEVER IMS.E*1 F(J) = SIGA/((SIGT1-Pl*(SIGSTH+SIGSU))*E)*Pl FN.(J) = (1-P2)*SIGA/(SIGT1-Pl*(SIGSTH+SIGSU.))/E J = J+1 OTHERWISE F(J) = SIGA*Pl/(E*SIGT1) PN(J) SIGA*(1-P2)/(E*SIGT1) J =.J+l END OF CONDITIONAL LOOP3 CONTINUE Fl: U FN1 = 0 THROUGH KILL. FOR M=1919 M.G.(J-1)/2 Fl = Fl+F(2*M-1) KILL FN1 = FNi+FN(2*M-1) F2 = 0 FN2 = 0 THROUGH EEN. FOR M=1,1, M.G.(J-3)/2 F2 = F2+F(2*IM EEN FN2 = FN2+FN(2*M) RI1 = H/3*0*(F(O)+4*F1+2*F2+F(J-1)) Ri.: H/3.0*(FN(O)+4*FN1+2*FN2+F(J-1)) PRIN. RESULTS RI1,RI29PHIZ9PHIBlPHIP19PHIP2 R INTERFERED RESONANCE INTEGRAL AMU = (1-ALPU)*(A(1,1)-1) AMTH = (1-ALPTH)*(A(1,1)-1) AMOX = (1-ALPOX)*(A(11l)-l) AMM = (1-ALPM)*(A(191)-1) LNU = ELOG.(1/ALPU) LNTH = ELOG.(1/ALPTH) LNOX = ELOG.(1/ALPOX) LNM = ELOG ( 1/ALPM) EXDU = EXP.(-DELU) WHENEVER LNU *G*DELU IU = CD/(1-ALPU)*(EXDU-EXP (-LNU)) + PHIZ/(i1-ALPU)*(1-EXDU) l +PHIP1* (EXDU-EX)/AMU OTHERWISE IU = PHIZ + PHIP1*(EXP.((A(l,l)-l)*LNU-A(191)*DELU)-EX)/AMU END OF CONDITIONAL WHENEVER LNTH.G-DELU ITH = CD/(1-ALPTH)*(EXDU-EXP,(-LNTH)) + PHIZ*(1-EXDU)/ 1 (1-ALPTH) + PHIP1*(EXDU - EX)/AMTH OTHERWISE ITH = PHIZ + PHIP1* (EXP*((A(1,l)-l)*LNTH-ACll)*DELU)-EX)/ 1 AMTH END OF CONDITIONAL WHENEVER LNOXI.GDELU IOX = CD/(1-ALPOX)*(EXDU-EXP.(-LNOX)) + PHIZ*(1l-EXDU)/(1

-1281 ALPOX) + PHIP1*(EXDU-EX)/AMOX OTHERWISE IOX = PHIZ + PHIP1*(EXP.((A(ll)-l)*LNOX-A(191)*DELU)-EX)/ 1 AMOX END OF CONDITIONAL WHENEVER LNM *G.DELU IMOD = CD/(1-ALPM)*(EXDU-EXP.(-LNM)) + PHIZ*(1-EXDU)/(1-ALPM) 1+PHIP2*(EXDU-EX)/AMM OTHERWISE IMOD = PHIZ +PHIP2*(EXP.((A(1,1)-l)*LNM-A(1,1)*DELU)-EX)/AMM END OF CONDITIONAL WHENEVER IMS.Eo1 RI = (SIGSO*IOX*RI1 + (SB-SA)/SA*SIG2*IMOD*RI2)/PHlB1 OTHERWISE RI = ((SIGSO*IOX + SIGSU*IU +SIGSTH*ITH)*RI1 +(SB-SA)/SA* 1 SIG2*IMOD*RI2)/PHIB1 END OF CONDITIONAL DELRI = RI-UNIRI PCERR = DELRI*100.0/RI PRINT COMMENT $OINTERFERED RESONANCE INTEGRAL$ PRINT RESULTS RI,UNIRI,DELRI tPCERR PCINA =(RI*PHIB1 - UNIRI*CD)/(UNIRI*CD) PRINT RESULTS PHIB2,PCINA TRANSFER TO START OUT1 PRINT COMMENT $OERROR IN TAB SUBROUTINE$ TRANSFER TO START INTEGER NFRNSRHM,HNM,NJ,K END OF PROGRAM $DATA XT=.00,.01, *02, *03,,04,.05,.06, *07,.08,.09, *10,.11,.12,.13, *14, *15,.16, *17,.18, 019,.20,.21,.22,.23,.24, *25, *26,.27,.28. *29,,30,.31,.32,.33,.349.35, *36,.37, *38,.39,.40, *41, *42,.43,.44,.45,.46,.47,.48,.49,.50,.51,.52, *53,.54,.55,.56,.57,.58,.59,.60,.61,.62,.63,.64,.65, 966,.67, -68, *69..70, *71,.72,.73,.74,.75, *76,.77, *78..79,.80,.81, *82,.83,.84, *85,.86, *87,.88, *89,.90,.91,.92,.93,.94,.95,.96,.97,.98,.99, 1.00, 1.01, 1.02, 1.03. 1.04, 1.05, 1.06, 1.07, 1.08. 1.09, 1.10, 1.11. 1.12, 1.13. 1.14, 1.15, 1.16, 1.17, 1.18, 1.19, 1.20, 1.21, 1*22, 1.23 24, 1. 25, 1.26, 1.27, 1*28, 1.29, 1.30, 1.31, 1.32, 1.33, 1.34, 1.35, 1.36, 1.37, 1.38, 1*39, 1*40, 1.41, 1.42, 43 44, 1.45, 1946, 1.47. 1.48, 1449, 1.50, 1.51, 1.52, 1.53, 1.54, 1.55, 1.56, 1.57, 1.58, 1.59, 1.60, 1.61, 1.62, 1.63, 1.64, 1.65, 1.66, 1167, 1.68, 1.69, 1.70, 1.71, 1.72, 1.73. 1.749 1.75, 1.76, 1.77, 1.78, 1.79, 180 1, 182, 1834, 1.85, 1.886 1.87, 1.88, 1.89, 1.90, 1.91, 1*92, 1.93, 1.94, 1.95 1196, 1-97, 1*98, 1.99, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2*7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3. 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.3, 4*4, 4.5, 4.6, 4.7, 4.8, 4.9, 5.0, 5.1, 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.3, 6.4, 6.5, 6.6, 6.7, 6.8, 6-9, 7.C, 7.1, 7.2, 7.3. 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8.0, 8.1. 8.2, 8.3, 8.4, 8-5, 86, 8.7, 8.8, 8.9, 9.0, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9,10.0, YT=.5000,.4902766,.4809683,.4719977,.4633239,.4549188,.4467609,.4388327,.4311197, *4236096,

-129*4162915,.4091557, ~4021937, ~3953977, ~3887607,.3822761,.3759380, *3697408, ~3636795, *3577491,.3519453,.3462638, *3407005, ~3352518, ~3299142, ~32468419, 3195585, ~3145343, *3096086, ~3047787,.30004189.2953956, ~2908374, *2863652, ~2819765,.2776693,.2734416, *2692913, *2652165, ~2612155,.2572864,.2534276, *2496373, *2459141, *2422563, ~2386625,.2351313,.2316612,.2282508, ~2248990,.2216044,.2183657,.2151818, ~2120516,.2089739, o2059475,.2029715,.2000448, ~1971664, ~1943353,.1915506,.1888114, ~1861166, ~1834656, o1808573,.1782910,.1757658, ~1732810, ~1708358, ~1684294,.1660612, ~1637303, ~1614360, ~1591778, ~1569549,.15476679, 1526125,.1504917, ~1484037, *1463479,.1443238,.1423307, ~1403681, ~1384355, ~1365324,.1346581,.1328122, ~1309943, ~1292037, ~1274401,.1257030, ~1239919, o 223063, *1206459, ~1190102, *1173988,.1158113, ~1142472, ~1127063, ~1111880,.1096920, ~1082179, ~1067654, ~1053342, *1039238,.1025339,.1011643, o0998145, ~0984842, *0971731,.0958809,.0946074,.0933521, *0921149, *0908953,.0896932, ~0885083,.0873402, o0861888, o0850537,.0839347,.0828315, o0817439, ~0806717, ~0796146, *0785723, ~0775447, ~0765316, ~0755326, ~0745476,.0735763, o0726186,.0716742, *0707429, o0698246,.0689191,.0680260, 90671453, o0662768, o0654203,.0645755,.0637424, *0629207, o0621104, o0613111,.0605227,.0597452, ~0589782, ~0582217, ~0574755,.0567395, ~0560135, o0552973, *0545908, o0538939, ~0532064, o0525283, o0518592, 90511992, o0505481, ~0499057,.0492720, ~0486467, ~0480299, ~0474213, 90468209,.0462284., o0456439, o0450672, o0444982, o0439367,.0433827, o0428361, o0422967, *0417645,.0412393,.0407211, ~0402097, *0397051, o0392071,.0387157,.0382308,.0377522, ~0372800,.0368139, o0363540,.0359001, ~0354521,.0350100,.0345737, *0341430,.0337180, o0332986, ~0328346, o0324759,.0320727,.031'6746,.0312817, *0308939, o0305112, 3.01334E-2, 2.66136E-2, 2o35207E-2, 2o08002E-2, 184054E-2, 1.62954E-2, 1.44349E-2, 1.27932E-2, 1.13437E-2, 1.00629E-2, 0.89306E-2, 0.79290E-2, 0.70425E-2, 6~25744E-3, 5-56190E-3, 4.94538E-3, 4o39865E-3, 3o91360E-3, 3o48310E-3, 3.10087E-3, 2.76136E-3, 2.45969E-3, 2*19156E-3, 1~95315E-3, 1.74110E-3, 1.55244E 3, 1.38454E-3, 1.23507E-3, 1.10197E-3, 0-98342E-3, 0.8778CF-3, 0.78368E-3, 6o9978E-4, 6o2498E-4, 5-5827E-4, 4.9877E-,+ 4o4569E-4, 3.9832E-4, 3.5604E-4, 3.1830E-4, 2~8460E-4, 2o5451E-4, 2*2763E-4, 2o0362E-4, 18217E-4, 1.6300E-4, 1.4586E-4, 1@3055E-4, 1@1685E-4, 10461E-4, 0*9366E-4, 0o8386E-4, 7o5100E-5, 6~7261E-5, 6~0247E-5, 5~3970E-5, 4o8352E-5, 4@3323E-5, 3o8821E-5, 3.4790E-5, 3~1181E-5, 2o7949E-5, 2o5054E-5, 292461E-5, 20138E-5, 1.8057E-5, 1.6192E-5, 14521E-5, 1.3024E-59 11682E-5, 1.0479E-5, 0o9400E-5, 8.4335E-6, 7.5668E-6, 6o7896E-6, 6.0927E-6, 5. 4677E-6, 4.9071E-6, 4,4044E-6, 3@9533E-6,3.5488E-6*

REFERENCES 1. A. M. Weinberg and E. P. Wigner, The Physical Theory of Neutron Chain Reactors, The University of Chicago Press, 1958. 2. R. V. Meghreblian and D. K. Holmes, Reactor Analysis, McGrawHill Book Co., New York, 1960. 3. A. D. Galanin, Thermal Reactor Theory, Pergamon Press, New York, 1960. 4. A. M. Weinberg, "Thermal Utilization in a Square Lattice," CP-330, 1942. 5. A. C. Clark and D. A. Newmarch, "The Diffusion Theory of Fine Structure in Thermal Neutron Reactor Assemblies Consisting of a Cylindrical Fuel Element Set in Square Lattice Array," AERERP/R 1657 (1955). 6. E. Richard Cohen, "The Thermal Neutron Flux in a Square Lattice Cell," Nucl. Sci. and Eng., 1, (1956) 268-279. 7. A. D. Galanin, "The Thermal Coefficient in a Heterogeneous Reactor," The Proceedings of the International Conference on the Peaceful Uses of Atomic Energy at Geneva(PICG), P/666, 5_, (1956) 477-483. 8. H. Neumann, "Thermal Utilization, Thermal Flux, and Lattice Cell Shape," HW-52048, 1957. 9. B. Bailly du Bois, "Influence de la forme des cellules sur le laplacien et structure fine du flux thermique dans une pile heterogene," Rapport C.E.A. No. 740, 1957. 10. A. Pazy and S. Goshen, "The Thermal Utilization in a Rectangular Cell," Nucl. Sci. and Eng., 11, (1961) 357-361. 11. A. Amouyal; P. Benoist, and J. Horowitz, "Nouvelle Methode de determination du Facteur d'Utilisation Thermique D'une Cellule," Journal of Nuclear Energy, 6, (1957) 79-98. 12. F. Seitz and E. P. Wigner, "On the Constitution of Metallic Sodium," Physical Review, 43, (1933) 804-810. 13. S. M. Feinberg, "Heterogeneous Methods for Calculating Reactors: Survey of Results and Comparison with Experiments," PICG, P/669, 5, (1956) 484-500. -130

-13114. G. N. Watson, Theory of Bessel Functions, Cambridge University Press, 1958. 15. W. A. Horning, "A Summary of Small Source Theory Applied to Thermal Reactors," HW-34021, 1954. 16. E. P. Wigner, E. Creutz, H. Jupnik, and T. Snyder, "Resonance Absorption of Neutrons by Spheres, " Journal of Applied Physics, 26, (1955) 260-270. 17. K. T. Spinney, "Resonance Absorption in Homogeneous Mixtures," BNL-433(C-24), (1956) 103-109. 18. J. Chernick and R. Vernon, "Some Refinements in the Calculation of Resonance Integrals," Nucl. Sci. Eng., 4, (1958) 649-672. 19. W. Rothenstein, "Collision Probabilities and Resonance Integrals for Lattices," BNL-563(T-151), 1959. 20. W. Rothenstein, "Collision Probabilities and Resonance Integrals for Lattices," Nucl. Sci, and Eng., 7, (1960) 162-171. 21. L. W. Nordheim, "A Program of Research and Calculationsof Resonance Abiorption, " GA-2527, 1961. 22. L. W. Nordheim, "A New Calculation of Resonance Integrals," Nucl. Sci. and Eng., 12, (1962) 457-463. 23. G. F. Kuncir, "A Program for the Calculation of Resonance Integrals" GA-2525, 1961. 24. R. Goldstein and E. R. Cohen, "Theory of Resonance Absorption of Neutrons," Nucl. Sci. and Eng.o, 13, (1962) 132-140. 25. J. G. Hill and G. W. Schaefer, referred to in Reference 26, (to be published). 26. H. M. Sumner, "ERIC-1, A Fortran Programme for Calculating Resonance Integrals, and some examples of its use," AEEWM304, 1963o 27. R. Schermer and N. Corngold, "Resonance Capture of Neutrons in Infinite Homogeneous Media," Proc. Phys. Society, Volo LXXIII, (1959) 561-571. 28. S. Iijima, "Resonance Absorption and the Resonance Disadvantage Factors," Nucl. Sci. and Eng., 17, (1963) 42-46.

-132290 K. M. Case, F. de Hoffman, and G. Placzek, Introduction to the Theory of Neutron Diffusion, Vol. I, Superintendent of Documents, U. S. Government, 1953. 30. A, Sauer, "Approximate Escape Probabilities," Nucl. Sci, and Eng., 16, (1963) 329-335. 31. F. Di Pasquantonio, "Calculation of First-Collision Probabilities for Heterogeneous Lattice Cells Composed of n-Media," Reactor Science and Technology (J. of Nuclear Energy, Parts A/B), 17, (1963) 67-51o 32. S. M. Dancoff and M. Ginsburg, "Surface Resonance Absorption in a Close-packed Lattice," CP-2157, 1944, 33. I. Carlvik and B. Pershagen, "The Dancoff Corection in Various Genometies, " AE-16, 1959. 34. J, A. Thie, "A Simple Analytical Formulation of the Dancoff Correction," Nucl. Sci. and Eng., 5, (1959) 75-77. 35. G. I. Bell, "A Simple Treatment for Effective Resonance Absorption Cross Sections in Dense Lattices," Nucl. Sci. and Eng,, 5, (1959) 138-139. 36. M. M. Levine, "Resonance Integral Calculations for U-238 Lattices," Nuclo Sci, and Eng., 16, (1963) 271-279. 37. F, T, Adler, G. W. Hinman, and L. W. Nordheim, "The Quantitative Evaluation of Resonance Integrals," PICG, P/1988, 16, United Nations, New York, (1958) 155-171. 38, W. T. Hayes, M. Luming, and P. F. Zweifel, "Resonance Escape, Flux Disadvantage, and Collision Probabilities," Trans. Amer. Nucl, Society, 3, (1960) 366, 39. Mo A. Mannan and P. F. Zweifel, "Solution of Diffusion Equations in Rectangular Cells," Trans. Amer. Nucl. Society, 6, (1963) 4. 40o E. T, Whittaker and G. No Watson, A Course of Modern Analysis, Cambridge University Press, (1962) 359. 41, S. Glasstone and M, C. Edlund, The Elements of Nuclear Reactor Theory, D, Van Nostrand Co., Inc., Princeton, New Jersey, 1952, 42. B, Arden, B, Galler, and R, Graham, Michigan Algorithm Decoder, University of Michigan, 1963.

-133 - 43. J. Chernick, "The Theory of Uranium-Water Lattices," PICG, P/603, 5, United Nations, New York, (1956) 215-228. 44. G. Placzek, "The Functions En(x) = exuu-ndu, I MIT-l(NRC No, 1547). 45. J. H. Ferziger, The Theory of Neutron Slowing Down in Nuclear Reactors, Thesis, University of Michigan, 1962. 46. L, Dresner, Resonance Absorption in Nuclear Reactors, Pergamon Press, London, 1960. 47. F, B. Hildebrand, Introduction to Numerical Analysis, McGrawHill Book Co., New York, 1956, 48, H. L. Brown, Jr., T. J, Connolly, and W. K. Foell, "Effective Resonance Integrals and Interference between Resonances of Indium-115, Gold, and Rhenium-185," Trans. Amer. Nucl. Soc., 5, (1962) 375-376. 49. W. K. Foell,'R. A. Grimesey, and S. Tong, "A Monte Carlo Study of Resonance Absorption in Gold and Indium Lumps," Trans. Amer. Nucl, Soc., 6, (1963) 272-273. 50. J. Ro Triplett, et al., "The RBU Reactor Burnup Code," HW70049, 1961. 51. H. L. Brown, and T. J. Connolly, (personal communication to Foell et al., Reference 49).