WCAP-2o48 THE UN IV ER SIT Y OF MI CHI GAN COTLLEGE OF ENGINEERING Department of Nuclear Engineering Final Report THE DOPPLER EFFECT FOR A NON-UNIFORM TEMPERATURE DISTRIBUTION IN REACTOR FUEL ELEMENTS Jack E. Olhoeft ORA Project 04261 under contract with: WESTINGHOUSE ELECTRIC CORPORATION ATOMIC POWER DIVISION PITTSBURGH, PENNSYLVANIA For the U. S. Atomic Energy Commission Under USAEC-Westinghouse Electric Corporation Contract No. AT(30-1)-3064 administered through: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR July 1962

This report has also been submitted as a dissertation in partial fulfillment of the requirements for the degree of Doctor of Philosophy in The University of Michigan. 1962

ACKNOWLEDGIMIENTS It is a pleasure to acknowledge appreciation to Professor RoKo Osborn for his inspiring lectures on reactor theory and nuclear physics as a teacher and for his wise counsel 'and stimullating discussions as an advisor throughout the course of this investigation, in addition sincere gratitude is expressed to the following. To Professor W. Kerr for encouragement, administrative assistance, and enlightened discussions concerning this work, To Professor BA. Galler and the Computing Center Staff for helpful discussions and use of the IBM 704 and 709 digital computers, To Professors GL. Gyorey, F,G. Hammitt, T.B. Kammash and J. G Wendel, also members of my doctoral committee, for helpful discussions and comments, The latter two were on sabbatical leave during the final phases of this work, To Doctors GoHo Minton, W.H. Arnold, Jro, SM, Hendley and REA. Dannels, Westinghouse Electric Corporation, for their interest in jihis problem and discussions concerning the nLumerical phaseso To Doctor R.A. Boyd and Mr, SoJ. Plice, Project Representative s of the Office of Research Administration; for their dedicated admiristrative i ork~ TPo Mr. G.E, Prosser and staff of the Office of Research Administration for the final preparation and production of this manus 1ripT,o To my wife for her assistance and patience ~i preparing this and other manuscripts, The author is especially grateful to the Atomic Ernergy Co nission and to the Westinghouse Electrie Corpora ion for providing the necessary financial support to accomplish this work. ii

TABLE OF CONTENTS Page LIST OF ILLUSTRATIONS v ABSTRACT vi CHAPTER I. INTRODUCTION 1 A. Background 1 B. Scope of Research 6 ZI. INTEGRAL TRANSPORT EQUATION AND NEUMANN SERIES 11 A. Derivation of Integral Transport Equation from Boltzmann Equation 12 B. Formulation of Neumann Series and Convergence Criterion 16 C. Representation of Reactor Parameters by Neumann Series 23 D. Mathematical Modification of Kernel and Source 27 III. RANDOM WALK PROCESS 30 A. Formulation of Random Walk Chains for Analog Processes 31 B. Formulation of Random Walk Chains for NonAnalog Processes 42 IV. MONTE CARLO APPROACH 48 A. Specialization to Slowing Down Problems in a Reactor Lattice 49 B. Specialization to Doppler Coefficient and Resonance Escape Probability 58 V. MONTE CARLO PROCEDURE AND COMPUTER PROGRAM 72 A. Description of Mobnte Carlo Procedure 72 B. Brief Description of Program 75 VI. RESULTS 79 A. Program Constants and Input Parameters 79 B. Numerical Results and Discussion 81 iii

TABLE OF CONTENTS (Concluded) CHAPTER Page VII. CONCLUDING REMARKS 91 A. Conclusions 91 B. Future Work 94 APPENDIX A. DERIVATION OF GREEN'S FUNCTION 96 APPENDIX B. MATHEMATICAL DESCRIPTION OF BANACH SPACE 103 APPENDIX C. PROPERTIES OF DENSITY FUNCTIONS AND OTHER PROBABILITY FUNCTIONS 107 APPENDIX D. PROOF OF FINITENESS OF NON-ANALOG RANDOM WALK CHAINS DESCRIBED IN CHAPTER IV 111 APPENDIX E. SCATTERING PROBABILITY FOR STRUCK NUCLEUS IN MOTION 115 APPENDIX F. EVALUATION OF THE DOPPLER BROADENING FUNCTIONS AND THEIR DERIVATIVES 125 APPENDIX Go COMPUTATIONAL DETAILS USED IN THE MONTE CARLO PROGRAM, REPAD, FOR SELECTING RANDOM VARIABLES 138 APPENDIX H. A LIST OF SUPPLEMENTARY COMPUTER PROGRAMS WRITTEN FOR THIS INVESTIGATION 146 REFERENCES 148

LIST OF ILLUSTRATIONS Table Page I Values of Input Parameters for REPAD Program 80 II Results from 12,000 Monte Carlo Histories 83 Figure I The variation of resonance capture probability with fuel temperature; uniform temperature case, runs l, 2, 3. 84 2 The variation of the Doppler coefficient of reactivity and the derivative of resonance escape probability with fuel temperature; uniform temperature case, runs 1, 2, 3. 85 3 The variation of ln(-lnp) with fuel temperature; uniform temperature case, runs 1, 2, 3. 86 v

ABSTRACT The Monte Carlo method is used in this work to investigate the temperature effects on the Doppler coefficient of reactivity and resonance capture probability due to a non-uniform temperature distribution in the fuel rods arrayed in a reactCor lattice. Scattering is allowed in the fuel and moderator, and Doppler broadening of the resonances is taken into account. This investigation begins with the Boltzmann equation describing the transport of neutronsin steady state. By transforming this equation into an integral equation, the solution can be expressed as a Neumann series, for which a convergence criterion is established. It is shown that the Neumann series can be used to obtain series equivalents for general reactor parameters and their derivatives which depend on the neutron collision density, A Monte Carlo process that traces neutrons in spatial propagation and energy degradation through a reactor can be related to random walk chains governed by the source and kernel of the integral equation. Also, once the relationship between Monte Carlo histories and the Neumann series it established, one can determine the proper form for random variables which depend upon the paths of the neutron histories in order to evaluate reactor parameters and their derivatives by Monte Carlo. By application of this technique a new and direct method is devised for evaluat ing the Doppler coefficient. vi

As an example of this method a computer program was developed for a cylindrical cell and a fuel rod with arbitrary radial temperature dependence, Several cases for different fuel temperatures were investigated using water moderator and U02 fuel. Results illustrating temperature dependence of the Doppler coefficient of ractivity and the resonance capture probability are reported for uniform temperature and parabolic temperature distributions in the fuel. The results indicate that this direct method of calculating the Doppler coefficient by Monte Carlo is practical, and allows more efficient use of computer time in comparison to empirical and finite difference approximations previously used. It is also concluded that the Doppler coefficient is more sensitive than the capture probability to a temperature distribution appropriate to power reactors. This Monte Carlo program also provides the means of investigating these same effects for temperature distributions appropriate to reactors in excursion conditions. vii

CHAPTER I INTRODUCTION A. BACKGROUND Nuclear reactors may contain in the fuel certain nuclides, such as U-235, U-238, Pu-239, and Th-232, which have neutron resonances that lie within the neutron energy spectrum present in reactors. If the temperature of the resonance material were at absolute zero, then each resonance would have a shape, called 'natural line' shape, with very sharp peak and very narrow energy width. Because of the thermal motion of atoms at normal reactor temperatures, the relative energy distribution of the neutrons and target nuclei changes as the temperature of the material increases and causes the resonances to broaden and lower. This phenomenon is called 'Doppler broadening', and is named after Christian Dopplerl who noticed similar effects in his original work with light waves. The Breit and Wigner2 single level formula for neutron resonance cross sections is similar to that for the natural shape of an optical line. In 1937 Bethe and Placzek5 discussed Doppler broadening of a neutron absorption resonance caused by the thermal agitation of the resonance atoms and derived expressions describing the temperature dependence of the cross sections. The effect of the Doppler broadening of resonances on the neutron reaction rate is called the Doppler effect and the resulting influence on the reactivity of a nuclear reactor is called the Doppler temperature coefficient of reactivity. The dependence of the Doppler coefficient of 1

2 reactivity upon fuel temperature makes it important in reactor safety considerations because a decrease in reactivity with increasing fuel temperature helps to make a power excursion self-limiting. U-238 produces just such an effect whereas fissile resonance absorbers cause the opposite effect. Also the capture of neutrons in the resonances of U-238 and the ensuing conversion to Pu-239 is an important factor in long term temporal studies of reactor criticality. It is reported4 that as early as 1937, Meitner, Hahn and Strassmann first surmised that neutrons were absorbed in resonances of U-238, and that soon thereafter Fermi and Szilard recognized that resonance absorption could be decreased by mixing the fuel and moderator heterogeneously instead of homogeneously. At Szilard's suggestion, Wigner, Creutz, Jupnik, and Snyder5 performed experimental and theoretical studies in 1941 on resonance absorption of neutrons in heterogeneous fuel-moderator systems that included measurements of the temperature dependence in uranium dioxide (U02). Since that time many authors have reported extensive work, both of theoretical and experimental nature, that has contributed to the understanding and quantitative evaluation of the importance of the Doppler effect in the design and performance of nuclear reactors. The work for thermal reactors has been reviewed by Creutz, Jupnik, Snyder and Wigner4 in 1955, by Wigner6 in 1956, by Sampson and Chernick7 in 1958, by Nordheim8,9 in 1959 and 1961, and by Pearce10O in 1961. The importance of the Doppler effect in fast reactors has been studied by Feshback, Goertzel and Yamanchi,1 Frost, Kato and Butler,l2 and Nicholson.l3 The remainder of this discussion will deal exclusively with thermal reactors.

The temperature dependence of the Doppler effect has also been investigated with Monte Carlo techniques by Morton,14 and Arnold and Dannels.15 The latter authors used the REP code developed by Richtmyer, Van Norton and Wolfe,16 which treats the slowing down of neutrons from a "1/E" source in a rectangular cell with cylindrical fuel rod. The code allows use of the Breit-Wigner single level formula for the resonances. Pearce1O and Hellstrand, Blomberg, and Horner17 have indicated that the discrepancy among experimental and theoretical methods in evaluating temperature dependence of the resonance integral of uranium metal, uranium dioxide, thorium metal and thorium dioxide lies outside the limits of experimental error. However, recently Nordheim18 has discussed a new method of calculating the resonance integral which assumes isolation of the individual resonances and involves numerical integration to determine average collision densities in the fuel and moderator regions. He has obtained quantitative agreement of the resonance integrals for uranium metal and uranium dioxide and their temperature dependence with the experimental results of Hellstrand, Blomberg and Horner.17 In the references on thermal reactors cited above, the Doppler coefficient is defined as a= 1 (RI) (1-1) (RI) aT where T is the fuel temperature, and (RI) is the effective resonance integral, which is related to the resonance escape probability, p, by the

following expression: p = e(RI)/ s (1-2) where ~as is the average slowing down power per resonance atom. Since the temperature coefficient of reactivity is given by (l/k) ak/6T where k is the multiplication constant of the reactor, then the contribution to the temperature coefficient due to Doppler broadening of the resonances is given by 1 ap a n p (1-3) p aT and is called the Doppler coefficient of reactivity in order to make the verbal distinction between Eq. (1-1) and (1-3). In theoretical studies in the past it has been the practice to obtain the Doppler coefficient by approximating the differential in Eq. (1-1) or (1-3) as a finite difference or by fitting the temperature dependence of the resonance integral by an empirical curve functionally dependent on T1/2. Monte Carlo has also been used in conjunction with other calculational methods in reactor analysis in order to obtain more accurate design information and knowledge of temporal behavior. Ritsi and Minton19 have used Monte Carlo to treat the slowing down of neutrons through the resonance energy region and thence to provide energy group constants for the Muft and Sufocate codes, which are multigroup codes for the epi-thermal and thermal energy regions, respectively. In addition, Leshan, Burr, Temme,

5 Morrison, Thompson and Triplett20 have combined Monte Carlo with diffusion theory and kinetic equations for fuel life-time studies. All of the works mentioned previously have considered only a uniform temperature in the fuel, however in operating power reactors the temperature distributions are non-uniform. In fact Belle21 gives evidence that the ratio of the peak to surface fuel temperature (absolute temperature scale) for uranium dioxide could be greater than three depending on power density and size of fuel rods. Under certain cases he reports temperatures of about 2740~C and 400~C, respectively, in a simulation of an operating power reactor with uranium dioxide fuel. There has been very little published work discussing the effects of a non-uniform temperature distribution. The case of a 1 cm slab of uranium metal with surfaces-at O0K and mid-plane at 900~K was studied by Keane22 who used a simplified model of the U-238 resonance structure and assumed neutrons to be normally incident on the slab. In addition, Dresner23 performed a similar study for a single resonance on an isolated cylindrical fuel rod with maximum temperature ratio of two, using resonance parameters typical of U-238 resonances. Both authors conjectured that the resonance capture probability for a parabolic temperature distribution can be approximated by that of a uniform temperature at the mean temperature, however both authors- neglected scattering and neither discussed the validity of the approximation in obtaining the Doppler coefficient of reactivity. The neglect of scattering could be a serious error since the neutron width is comparable or greater than the capture width for at least half

6 of the resolved resonances of U-258 and the expected value of the neutron width is greater than the capture width for the unresolved resonances. Furthermore the absolute zero temperature used by Keane does not have physical significance for metal or oxide fuels as noted by Lamb.24 The Monte Carlo method is an attractive alternative to the analysis of Keane and Dresner, since the problem of a non-uniform fuel temperature can be solved by Monte Carlo for a realistic reactor simulation without making the assumptions and approximations introduced by them. For our purpose, Monte Carlo can be thought of as a mathematical experiment, involving the use of random sampling procedures, for simulating and studying a physical process or solving a deterministic mathematical problem. The physical process that we are studying is the slowing down of neutrons through the resonances of fuel in a heterogeneous reactor system. This process can be described by the Boltzmann transport equation, which is used to formulate the mathematical experiment. Alternatively, the Boltzmann equation could be solved by standard numerical techniques, but because of the resonance structure of the fuel and the geometry of the reactor lattice, there is a danger that the systematic errors or possible instabilities inherent in numerical analysis would mask the small but important effects that we are looking for in this problem. Monte Carlo avoids this difficulty and can be made more determinably accurate by increasing the number of histories processed. B. SCOPE OF RESEARCH The Monte Carlo method is used in this work to investigate the tem

7 perature effects on the Doppler coefficient of reactivity and resonance capture probability due to a non-uniform temperature distribution in the fuel rods arrayed in a reactor lattice. Scattering is allowed in the fuel and moderator, and Doppler broadening of the resonances is taken into account. A new and direct method of evaluating the Doppler coefficient by Monte Carlo is obtained and discussed in succeeding chapters. The investigation in Chapter II uses as a foundation the Boltzmann equation describing the transport of neutrons in a steady state condition. By use of the Green's function technique, the Boltzmann equation is transformed in this study to a transport integral equation, the solution of which is written in terms of the Neumann series. Although integral equivalents of the Boltzmann equation for specific conditions have been presented elsewhere25,26 the derivation presented in this investigation is kept general in order that the resulting integral transport equation be applicable for the same general conditions which apply to the Boltzmann equation. Furthermore this general form is required for our study of the problem of non-uniform fuel temperature. By using the mathematical properties of Banach abstract spaces, we can investigate the convergence of the Neumann series and establish a convergence criterion, which is required in the work of succeeding chapters. In view of the fact that the purpose of this study is to include investigation of reactor parameters, such as the resonance capture probability and the Doppler coefficient, a general parameter is represented by the integral of the product of a function and the collision density. The de

8 rivative of the parameter equals the derivative of the integral. It is then shown in Chapter II that the Neumann series can be used to obtain series equivalents for the integral and its derivative, and that these series converge for conditions that are usually satisfied in physical syst tems. A perusal of these Neumann series representations, which are given in Chapter II, indicates that closed form solutions are intractable because of the multiple integration involved in the terms of the series. However, these series can be evaluated by Monte Carlo. In 1953 Albert27 introduced a stochastic model for describing a random walk process, and in 1960 Spanier28 studied and generalized this model by using the rigorous mathematics of modern probability theory. Both authors discussed the application of Albert's model to the study by Monte Carlo of the solution of a Fredholm type integral equation by assuming uniform convergence of the Neumann series associated with the integral equation. Spanier recognized applicability to neutron transport problems under the assumption of convergence. Since the convergence of the Neumann series solution of the Boltzmann transport equation is established in Chapter II, the basic stochastic model of Albert for a general random walk process is used in Chapter III. Also in view of the convergence of the Neumann series representations of reactor parameters and their derivatives, the application of Albert's model to neutron transport is reformulated and extended to the evaluation of reactor parameters and their derivatives. The investigation presented in the second and third chapters results in the complete and straightforward formulation of a rigorous stochastic model

9 for evaluating reactor parameters and their derivatives by Mcnte Carlo. The work in this chapter is kept general; in fact, the results can be applied to the behavior of gamma rays for shielding problems. The stochastic model formulated in Chapter III is specialized in Chapter IV to the class of slowing down problems in a reactor. In this case the neutrons are emitted from a given high energy source and are processed by Monte Carlo according to random walk chains that are specialized in this chapter to specifically describe energy degradation and spatial propagation through the fuel and moderator regions. Since the type of collisions are restricted to those in which a neutron loses energy but cannot gain energy, this study is not applicable to the thermal energy region. However, this restriction is compatible with ts-he study of temperature coefficient and capture probability due to Doppler broadening of the resonances in uranium fuel. The Monte Carlo treatment for this problem is formulated for a general temperature distribution in the fuel rod. Previous to this work, the Doppler coefficient was obtained by an empirical fit to a curve of capture probability plotted versus temperature, or by a finite difference approximation to the differential, as described in Section I-A. This was done only in the case of a constant fuel temperature but could be applied to the case of a non-uniform temperature distribution. By applying the results of the investigation in Chapter III, a method of determining the Doppler coefficient directly by the Monte Carlo process is devised. In this way the Doppler coefficient is no longer dependent upon the inherent uncertainties of the two previous methods. Also

10 as discussed in Chapter VI, the direct method allows a very considerable saving of computer time in determining the Doppler coefficient to a given accuracy. As an example of the method given in Chapter IV, a computer program was written for a circular cylindrical cell geometry. This cell geometry was used instead of a rectangular cell geometry in order to save computer time. This savings in time is possible because neutron paths incident on the cell boundary are reflected along a chord similar to the incident chord. The results are presented for a parabolic temperature distribution, water moderator and U-238 oxide resonance material. This work results in a new method for determining the Doppler coefficient by Monte Carlo that is founded on a systematic deductive approach. Although the method could conceivably be formulated by intuition alone) evidently no one has done it. Furthermore the mathematical development avoids the introduction of assumptions that can be justified only by intuition. The method is shown to be applicable to physical systems of interest in reactor studies.

CHAPTER II INTEGRAL TRANSPORT EQUATION AND NEUMANN SERIES In this chapter an integral equation for the transport of neutrons is obtained from the steady-state Boltzmann equation. If one writes the integral equation in terms of neutron collision density rather than neutron flux, and represents a point in the 6-dimensional phase space (neutron energy, direction, and position) by the variable T. then the integral equation may be written in standard formo In Section II-B the Neumann series solution of the integral equation is determined by iteration. In investigating the convergence of the resulting Neumann series the more conventional techniques can not be relied upon since the kernel of the integral equation does not belong to Hilbert space. However the kernel does belong to the more general abstract space, namely Banach space. By using the appropriate definitions of norms and properties of transformations in Banach space the convergence for the Neumann series for neutron transport can be rigorously proven. The most important physical condition for the convergence is that the ratio of scattering to total cross section have a maximum value les's than one. This is a reasonable condition for realistic reactor or shielding systems. 11

12 A. DERIVATION OF INTEGRAL TRANSPORT EQUATION FROM BOLTZMANN EQUATION The Boltzmann equation describing the transport of neutrons serves as the foundation for this work. For the steady state case it is given by,26 V ~ Q _ ~(r,E,_) + ZL(r,E) ((r,E,Q) = S(r,E,_ ) (2-1) + EI(rET,,) 'Zs(E,'-E'Q,_';r) dQ' dE' E SQ. where the following definitions and identifications apply: n(r,E,Q) d3r dE dQ - The expected number of neutrons in d3r about r with energies in dE about E and going in the direction d2 about 2. t(_ (r,E,_) n(rE,)v(E) = Neutron "flux" v(E) = Speed of neutron of energy E S(r,E,Q) = Neutron source.t(r,E) = Total macroscopic cross section Zs(E,IjE',w';r) = Differential scattering cross section In order to derive the integral equation equivalent to (2-1), consider the Green's function G(r,E,Q2r',E',2') satisfying the differential equation

13 V QG + Zt(r,E) G(r,EdQlr't,E,Q' ) (2-2) =6(r - r') 6(E - E') 62( Q_') and the boundary conditions G = 0 at Ir - r = (2-3) The notation 6 is defined as the Dirac delta function and the notation 62 is defined as 1 2- ) 2-t '(Q -! ) (2-4) Let LL - Q2 * (2-5) then it is noted that (/ 2( Q ) dQ = 52( ~ _Q') dQ' Q Q2' (2-6) 1 2i = 1_ 6(L - 1) dpL dD = 1 2jT BL=- 1 e=O Equation (2-2) is solved in Appendix A by Fourier Transform techniques. The solution for the Green's function is found to be _, I I-'1 _ 7t(E,r' + R')dR G(r,E,q r',E', ') d3r' dE' d' e (2-7) d3r' dE' d2T - -- r - '1

14 r -r r- rt where R (2-8) K-ir R Let H(r',E',_') represent the right hand side of equation (2-1), then by superposition, ~(r,E,2) = ///G(r,E,Q2r',ETQ') r'E'Q' (2-9) H(r',E',Q') d3r' dE' dQ' is an integral equivalent of Eq. (2-1). The method of proving equivalence is given by Churchill29. Equation (2-1) can be obtained by multiplying Eq. (2-2) by H(r',E',2') and integrating over all r',E','. Substitute the right hand side of (2-1) for H(r',E',W') into (2-9). One can then integrate over E', and ', employing the properties of the Dirac delta function, and obtain lr-r'i S- - CZt dR' ir - r'1 2 2 R(2-10) + (rt,E',Q'_ ) 7s(E,_lEt',Qt;rt) dE" dQ" E" 2" Equation (2-10) can be put in the form of a standard integral equation by introducing the following definitions: = Collision density

15 Ir-r' 1 fl,,;. - Zt(E'"' -. RR' ) dR' T(rjr';E,_Q) d3r e -RR') dR Zt(E,r) -_.3_.._ 2(a. _R) d - r Transport kernel S-(rE,,Q) = S(r',E,Q) T(rr';E,Q') d3r' r' = Collision density source C(E,Q|E',Q';r')dEdQn - 1-, ----~ C(E.,QE',Q', r ) dEdQ Zt(rE') s Collision kernel ~K(rE,',E' -r C(E,Q|E',_;r') T(rlr';E,Q) Kernel of integral equation Finally the integral form of the Boltzmann transport equation can be written for the neutron collision density as,(r,E,a) = S(,E,a) (2-11) + | fr(r',E',') K(r,E,Qr',E',Q') d3r' dE' dQ' r t Q It is further noted that the following relations exist: /T(rIr';E,Q) d3r = 1 (2-12) r C(E,_E',';r') dEdQ = ) (2-13).t(Erl )

hence K(r,E,Q2r',E',2') d3r dEd2 =.(E — (2-14) ~~~r E 2~ ~,t(E,r') The integral transport equation for neutrons can also be obtained directly by physical arguments if one notes that the kernels be expressed as the following conditional probabilities: K(r,E,Q2r',E',Q') d3r dEdQ -The probability that a neutron with energy E' and going in direction 2' shall, upon suffering a collision at r', be scattered into dE about E and dQ about g after that collision and shall suffer its next collision in d3r about r. C(E,I;E',2';r') dE dQ - The probability that a neutron with energy E' an't going in direction Q' shall, upon suffering a collision at r', be scattered into dE about E and d~2 about 2 after that collision. T(rlr';E,Q) d3r - The probability that a neutron at r' with energy E and going in direction Q shall suffer its next collision in d3r about r. Zs(r,E)/Zt(r,E) _ The probability that a collision suffered by a neutron with energy E at r was a scattering collision. B. FORMULATION OF NEUMANN SERIES AND CONVERGENCE CRITERION The integral equation for neutron transport can be put in a more compact form for use in this section by letting T represent the 6-dimensional vector of phase space in Eq. (2-11), XV(T) = ((T) + K(TT) (T) dT 2-15)

17 where dT = d3r dEd2 (2-16) The Neumann series solution of Eq. (2-15) can be obtained by successive iteration.30 The first iteration is found by substituting the right hand side of (2-15) for $(T') appearing in the integral term of (2-15), obtaining = S() + K(TTO) 'S(TO) dTo T 0 (2-17) + K(TIT1) K(T 11T),t(To) dT0 dT1 1 0 The second iteration is formed by substituting the right hand side of (2-17) for V(T') into (2-15). Higher order iterations are formed by successively substituting the previous iterative value of $(T) into the integral term of (2-15). After m iterations there results, m V(T) = X $n(T) + Remainder (2-18) n=o where to(T) = S(T ) (2-19) and for n > 1 ( n-1) 1TK( 1IT) T S(T ) dTo **dTn (2-20) n- 1 o nT- T and also

18 Remainder. K(TITm)...K(TIjTO) *(To) dTo...dTm (2-21) Tm To The Neumann series is then formed by letting m approach infinity and neglecting the remainder term, 00 *(T) = n(-T) (2-22) n=o It is remarked that the nth term in (2-22)physically represents the contribution to the collision density at T due to source neutrons which have undergone exactly n scattering collisions. The question of convergence of the Neumann series has been thoroughly studied for the special case that Eq. (2-15) is a Fredholm integral equation of the 2nd kind and which has elements belonging to Hilbert L2 space.31'32 In that special case the Neumann series converges if the norm of the kernel is less than one, i.e., if JiKllIIJ I K( TIT ) dT Td < 1 (2-23) T T j Unfortunately, the kernel for the neutron transport integral equation is not square integrable, which is a necessary condition for the solution of Eq. (2-15) to belong to L2 space. However, the criterion for the establishment of convergence of the Neumann series can be obtained by generalizing to Banach space (LP space). Actually the functional space LP(1 < p < a) and the space of continuous

19 functions C are particular cases of the Banach abstract space. Briefly a Banach space is a set of abstract elements, and is a linear space, a complete space, and has a norm (positive) for every element of the abstract set. The elements of the set may be points, measurable functions,* transformed functions or functionals (operators). The mathematical description of a Banach space is given in Appendix B. For a Banach space the norm of a measurable function, f, defined in the region e is30 norm - fll - If(x)I ax (2-24) in LP space, 1 < p < o, or norm - lfll -max If(x)l (2-25) in C or L space. The norm of a linear transformation, T, belonging to Banach space is defined as norm ITII - Ma (2-26) where Ma is the smallest of the bounds M, and where |ITfll S MIIf l (2-27) *Functions which are continuous or have a finite number of finite discontinuities are measurable. Functions that concern us here, namely those which describe physical processes, are measurable. For a rigorous definition based on set theory the reader is referred to Halmos33 or Doob.34

20 The convergence criterion for the Neumann series for neutron transport can now be established within the framework of Banach spaces and linear transformations. Writing Eq. (2-15) as a functional equation: *(T) = S(T) + X8((T) (2-28) or (I- ) S = (2-29) where I is the identity transformation and kis the linear transformation, associated with the kernel K(TIT'), which transforms the measurable function ~ into the measurable function Hi, i.e., K( T) - K(T|T') (T') dT' (2-30) T the Neumann series in this notation becomes (T) X( -( -) = s+) ~+"-) + 3, +:... (2-31) n=o The elements of the integral Eq. (2-15) Or (2-29)) belong to Banach L1 space under the conditions that the norm of S exists and is finite and the norm of K is bounded. In fact the source is a non-negative quantity; therefore, for L1 space, = I |( T)I| dT = 3(T) dT (2-32) T T After referring to relation (2-12) and the definition of the collision

21 source, inverting the order of integration in r and r', and then integrating over all r, one can write the norm of the collision source as |S||= S(T) dT (2-33) T The norm equals one if the neutron source is normalized, and at least is finite for physical problems of interest, i.e., the class of slowing down problems in a reactor. Later it will also be required that the integral of the derivative of S be bounded. According to (2-26) and (2-27) the norm of the kernel is found from the norm of the function which was obtained by operating the kernel on the source term. I S = I WS' dT = J | K(TIT') S'(T ) dTdT (2-34) T T T < n JIK( TIT') S'(T')| dT' dT (2-35) T TI By employing Fubini's Theorem (page 84, Riesz and Nagy,30)the limits of integration in (2-35) can be interchanged. After integrating over all T and noting (2-14), there results s'I T ~|. 11 diT' < Zj~max * 111 (2-36) Therefore the norm of the kernel is no larger than the largest ratio of scattering cross section to total cross section.

22 1= Ma | S T < 1 (2-37) a- M Z<t(T) max Condition (2-37) would not be satisfied for a pure scattering system (e. g., pure Compton scattering for photons) but would be reasonable for true physical systems. In requiring conditions (2-37), the norm of the kernel is bounded. Also, as shown, the norm of the source exists, and therefore the solution of integral transport Eq. (2-15) belongs to Banach L1 space. This condition will now be shown to satisfy the criterion for convergence of the Neumann series. Later we will require the norm of the derivative of the kernel to be bounded. Consider the right hand side of Eq. (2-31) and denote the partial sum by Sn, then the norm of the partial sum is flS[1 = Is +S+ + (2 +8) I I | |=|S + /S + 32S + S 3S+S| (2-38) By using the properties of linear transformations(Eqs. (B-6) through (B-12)Appendix B), the norm of the partial sum becomes [s~n] _< lS~ll ~ (1 + t + Il +.2 + Iw 11n) (2-39) which converges as n approaches a. This can be seen by employing condition (2-37) and comparing the resulting series to the geometric series. Hence the Neumann series (2-31) is "strongly" convergent("strong convergence" in Banach space is a generalization of "convergence in the mean" in Hilbert space.30) Denote the limit of the partial sum Sn by j. By operating on the Neumann series with the transformation R and utilizing

23 the continuity property of linear transformations (Appendix B), one can show that Eq. (2-31) is indeed a solution of the integral Eq. (2-15) ' = S' ( + + ) -S + S + +: ' +2 + +.o (2-40) -S + S Q.E.D. C. REPRESENTATION OF REACTOR PARAMETERS BY NEUMANN SERIES In this work we are primarily concerned with evaluating reactor parameters such as resonance capture probability, resonance escape probability and Doppler coefficient. These parameters can be expressed in terms of the Neumann series solution to the integral transport equation. In this chapter let w(T) represent a general function which is bounded and has finite derivative, a.e., and which when "averaged" with respect to the neutron collision density, gives the desired parameters of interest. The quotation marks are used because, although the collision density can be thought of as a probability density function, it does not necessarily satisfy the requirement of normalization. Let I represent the integral I = w(T) (T) d (2-41) T then the true average of w(T) is < w(,) > = I (2-42) f A(T) dT T

24 Means of evaluating the integral I in terms of the Neumann series will now be considered. We will again rely on the mathematical theory of Banach space. Since it has been established that the elements of the integral Eq. (2-15)belong to L1 space, and consequently so does its solution *, it is necessary that the measurable function w(T) belong to Lw space in order that I represent the scalar product of w and V in Banach space.30 The norm of w must exist in order for w to belong to LX space. Then according to (2-25) the norm of w equals norm = Ilwj| = max I w(T)| (2-43) and must be finite. Now substitute the Neumann series for *(T) into the integral I (2-41), giving I = S ~w[S +y + + +....] dT (2-44) T According to Riesz and Nagy,30 strong convergence of the Neumann series implies "weak" convergence of the Neumann series, which means that (2-44) can be integrated term by term and the resulting series converges strongly to I. Therefore, I = w(T) S(T) dT T (2-45) 00 i J '"f W(Tn)K(T-nITnl) K( TT) (T T0) dTo... dTn n=l Tn T0 Let a be a variable upon which w, S, and K depend parametrically. (The

25 temperature of the fuel would be an example of a since neutron cross sections in general depend on the temperature of the medium.) Rewriting (2-45) as an explicit function of a, I() I = /w(T,,a) S(r,a) dT T (2-46) 00 + X J. W(,T a) KK Tn |n1)a) )..K(T 1To1a) S( OY do.dTn n=1 Tn To We will require that the conditions previously satisfied by the norms of w, S, and K hold for all values of interest for the variable ac Then by using the Weierstrass Test for uniform convergence3)'36 the Neumann series equivalent of the integral I (2-46) can be shown to be uniformly and absolutely convergent with respect to a over an interval in which the terms are continuous functions of a. II(,) wa. | S dT dT + dT. (2-47) T T < Iw(TNa)l max ( ' S dT + I dT + SIdT +..) (2-48) But according to the definition of the norm of the source and the kernel for L1 space, and by again using the properties of linear transformations reported in Appendix B, series (2-48) becomes I SI(C I Iw(T,)[max I SII (+ 1 + I [ |2 + *. (2-49) which converges since (2-39) converges as n approaches a. Therefore, Eq.

26 (2-46) converges uniformly since the function w(T,hc) is bounded and the Neumann series exhibits strong convergence. In later chapters the differential of I (2-41) with respect to the general variable will be needed. aI(e) _ = ~ (w(Tc-) (T,ea) ) dT (2-50) T assuming the limits of integration are independent of the variable a. In view of the uniform convergence of the Neumann series equivalent of I, Eq. (2-46) can be differentiated term by term if the resulting series converges uniformly. Performing term-wise differentiation, ic - wL - S-'w'("o) '(TO) dT. eI wo aa so a, JT TO + oan '+a + (2-51) n=l Tn To + K1 v + I o a, wn Kn —K1 So dTo ' 'dTn where the following notations are used Kn K(TnIT nl,0) (2-52) wn -- (Tn,a) (2-53) By a- ion t ) ( 2-4) By again using the Weierstrass Test and by imposing the requirements that

27 the derivative of w is bounded, that the integral over all T of the derivative of S is bounded, and that the integral over all rn (Tn = (rn,En,2n)) of the derivative of the kernel, Kn, is bounded for all n, then series (2-51) can be shown to be uniformly convergent with respect to C over an interval in which the terms are continuous functions of a. Under these conditions the absolute values of the terms of (2-51) can be shown to be less than the terms of the following geometric series: ~I(___ ) < M n S n=O which is convergent by the Ratio Test,35'36 for condition (2-37). M is a finite constant. D. MATHEMATICAL MODIFICATION OF KERNEL AND SOURCE The remainder of this chapter will consider the modification of the kernel in a mathematical sense, and the resultant effect on the Neumann series. It is noted that if one thinks of t as just a mathematical density function, then a new density term 4' can be introduced into the integral I (2-41) in the following mannerI = wE(T) M/ (T.) f (T) dT = JW(T) W(T),'(T) dT (2-55) T T where the weight factor W(T) is defined as W(T) - -(T)/t'(T) (2-56)

28 Hence the integral I represents in this second form the "average" of the product (w W) with respect to the "density" term r. Likewise the source and kernel in the Neumann series representation of I can be modified without effecting the sum if the correct weighting factor is introduced. Let S and K' be the new non-physical source and kernel respectively, then I = /W(T) Wo(T) S'(T) dT T (2-57) 00 X+ j j W(Tn)Wn(To...Tn)K'(TnITni) ' K'K(TlTO)Y'(To)dTo'dTn n=1 Tn TO where the weight factors are now defined as Wo(T) S(T)/'' (T) (2-58) and for n > 1, K(T nIT n)...K((TITo) '(To) Wn( To0' — n K= (2-59) K1(' nT n1)...K1'7T, Io) r (TO) Similarly the differential of I can be written in the form aI( ) /j 1wa + S a)w(T) Wo(T) S'(T) dT 00 ~+ — n + - +.. 9+ -- (2-60) wn -a Kn au K1 aa n=1 on TO.+ ~o ) WnWn(To,* Tn) K'n..K'1 'o dTo dTn where the notation Kn and ' are similarly defined as in (2-52) and (2-54)~ The new kernel and source may be completely arbitrary without

29 effecting the value of I or aI/6o as long as the weight factors are defined as above. However in practice it would be reasonable to require the norm of new source to exist and the norm of the new kernel to be less than one, and both norms to belong to L' space, so that the corresponding Neumann series would be convergent. This discussion on the modification of the source and kernel in this chapter serves as insight to the nonanalog random walk process that is examined in the next chapter.

CHAPTER III RANDOM WALK PROCESS In this chapter the mathematical description for a random walk process is formulated, and applied to the problem of neutron transport. The Neumann series and the random walk process are then related. In 1953 Albert27 presented a mathematical model based on stochastic process theory to describe a random walk process and discussed specialized applications to the problem of estimating the solutions of integral equations. Spanier28 reformulated and extended his work by using a clear and precise mathematical development based on modern probability theory. Both authors make the basic assumption that the Neumann series converges, which is needed to develop their mathematical theory rigorously. Since the criterion for convergence of the Neumann series solution to the integral equation for neutron transport was developed in Chapter II, the basic mathematical model of Albert can be used in this chapter and related to the Neumann series. Monte Carlo is a convenient method for processing random walk chains. This is accomplished by making decisions based on random variables selected from probabilities governing the random walk chains. These decisions determine the generation of the random walk process and in turn the evaluation of random variables that serve as estimators to desired parameters which depend upon neutron transport phenomena. As will be seen in this chapter, proofs of convergence that were given in Section II-C are needed 3o

31 to determine the proper random variables, called estimators, which provide an estimate of the quantities given by the integrals I(2-41) and aI/ca(2-50) in the preceding chapter. It will also be noted that Monte Carlo effectively provides a method of evaluating the individual terms of the Neumann series corresponding to the integral I(2-45) and aI/ao(2-51). Because this chapter draws upon probability theory as a background, using definitions of probability density functions (or frequency functions), distribution functions, conditional density functions, and expected values, all in the multivariate sense, these definitions are presented in Appendix C. A. FORMULATION OF RANDOM WALK CHAINS FOR ANALOG PROCESSES Random walk processes can be divided into two categories called analog and non-analog. Analog processes can be described as processes in which random walks are governed by probabilities simulating the physical laws of the system under consideration. Non-analog processes are those in which the probabilities are deliberately but carefully distorted. This distortion will then permit the adoption of special Monte Carlo techniques. A random walk chain can be described as follows~ Consider a particle propagating through a medium in a random manner and suffering collisions at successive points in phase space, T0, T1, T ~2,~. It may be absorbed at any one of the points, say Tn, with probability P thus terminating the random walk chain, or it may continue on to the next point Tn+1 with probability Qo The selective process is continued until termination occurs. The mathematical model can be developed by considering a random walk

32 chain consisting of a sequence of points [Tn) at which are assigned definite probabilities of continuing or stopping the chain. Let qn(To,' ~ Tn) = Probability of continuing the chain with sequence (Tn) beyond point Tno Pn(To, —.Tn) = Probability of terminating the chain with sequence (Tn) at point Tn. where Pn(To.-~Tn) + qn(To,-ooTn) = 1, for all n and (3-1) 0 <Pn(To..Tn) ~ 1, for all n Also let fO(To) = the probability density function of To, n=O and fn(Tn ITo,''Tnl) = the conditional probability density function of Tn, providing the particle was not absorbed at any preceding point of the chain, n > 0. It is noted that f (T nTc.. 1T ) is written as a density function condin n 0 n-1 tioned on the previous points of the chain and also that both Pn and qn are allowed to depend on all previous collision points. This dependence is actually more general than we need in later work. In accordance with these definitions a random walk chain can be constructed by the following process: 1. Select To from the density function fO(To). 2. Decide termination (with probability Po(To)) or continuation (with probability qO(To)).

33 3. If continuation is selected then select T from the conditional density function fl(T1IT ). 4. Decide termination (with probability p (T,T )) or continuation (with probability ql(To T1)). 5. If continuation is selected, continue the process by selecting decisions from appropriate probabilities of the next higher subscript until termination occurs. In general, if previous n decisions are of continuation type and the sequence (TO0,T,''T ) has been generated, then select Tn from the conditional density function fn(Tn To Tn_ 1). Let us assume that termination occurs at Tn with probability pn(To)...Tn), this random walk chain may then be represented by the notation Cn = (Tqo(T)), (T1, q (T1)) n(Tnqn l(TO, 'Tnl), (Tn-Pn(To,'Tn ) )T The length of the random walk chain is defined as the subscript corresponding to that phase point, among the sequence (T,T1,T2,.), at which a termination event was selected. Let L represent the chain length, which is equal to n for the random walk chain, Cn. It is noted that the chain must be of finite length in order to represent physical problems. The probability density function, fn(To,.'*Tn), for the chain of length n and sequence (T n can be determined by an iterative procedure from the n

conditional density function governing the chain. It is in fact defined by fn(To...Tn) = fn(TnTo'''n- fn-1(To...Tn- ) (3-2) By successive application of the definition, then fn(To'***Tn) = fn (TnITo. ) f n-(Tn-1 To, -'T ) (3-3) "f (T ITon) ~ f ( T1 o) ~ o(T) ( T fm(TmlTo'.. Tm-1 fo(T3) (3-4) m=l A random walk process consists of a sequence of random walk chains usually of variable length. Therefore, a random walk process is represented by a sequence of density functions (fn(To,' 'Tnn)) and a sequence of functions (qn(To,..Tn) l and (pn(To'...Tn)]. For a random walk process certain valuable probabilities can be obtained. Let P(L=nj -The probability that the chain length is exactly n. Descriptively this would be the probability that among all possible chains the first n-l collisions that occur are scattering collisions and that the nth collision is an absorption collision. Mathematically for the random walk process this probability is

35 P(L=n ] = f (T) q (T) dTod f (T To)q ((T T) dT TO T1 fn. (T n-1 To,l Tn.2'n-) qn (To'Tn-,) dT n- i (3-5) fn( nITo I 'T, *Tn) Pn(To ' 'n ) dTn Tn By using Eq. (3-3) this can also be written in shortened form as P([ L=n]} = / fn(To. Tn) Pn(To5-.. Tn) iT qi(To-.. Ti) dTo..dTn Tn TO i=o-6) (3-6) In accordance with the restriction that random walk chains must be of finite length, we must have P {L=oo) = 0 This is a natural condition to impose in order to process random walk chains by Monte Carlo. An analog random walk for the neutron transport problem can be generated by using the kernel and source term of the neutron integral transport Eq. (2-15). For convenience, assume normalization of the source and make the following identifications: Tq ZtT (..)=/K Z|T)ts(Tn) n(To. Tn) K( TTn) dent of chain points To n and In this case, q~(T0,''T ) iS independent of chain points T,'' ~T _, and

36 physically represents the probability of the neutron suffering a scattering collision at T given that the collision occured there, and therefore n the probability of continuing the chain. Furthermore from Eq. (3-1), pn(To... Tn) represents the probability of an absorption collision. pn(To, Tn) = 1 - qn(TO.- 'T n) (3-8) The conditional density function is given by ( ) K(TnITnln) t(Tn-1) f (T IT )0**TKn(To ''T 3)9 n n o'n- 1 qn(T. T) Z (T n n-1 0n-1 n s n-1 which is normalized as can be seen by integrating over all Tn. The probability density function for n=0 is related to the neutron source as follows, assuming that the source is normalized, fo(TO) dTo = S( To) dTo (3-10) Upon substituting (3-7) through (3-10) into (3-3), the probability density function for n > 0 is given by fn(ToT, -.Tn)dTo..dTn: S (Ti) K( T o)dTo'' d Tn (3-11) i=o 0 (Ti) Eqs. (3-10) and (3-11) are also normalized as can be shown by integrating over all To through Tn and under assumption that the source is normalized. The procedure previously mentioned for constructing the random walk chains can now be used for this analog problem. However in order to use

37 Monte Carlo for processing these random walks, it is necessary to show that this analog problem produces chains of finite length with probability one. Equivalently we can show that the probability of a chain of infinite length is zero by using the convergence criterion of the Neumann series. From Eq. (3-6) the probability that the chain is of infinite length is just m P(L=oo) lim... fm(To,'' Tm) 77 qi(To, 'Ti) dTo **dTm (3-12) Tm TO After substituting Eqs. (3-11) and (3-7) into Eq. (3-12) one obtains, after simplification P(L=oo) im.. K(T1 (T ) dT0 idTm Tm To (3-13) m lim + K(Ti+ 1 i S(T0) dT m+l i=o Since the source and kernel are non-negative quantities the inequality of Eqs. (2-3 ) and (2-37), can be used in Eq. (3-13), obtaining P(L=a)] < lim i S( T) dT (3-14) = 0 It is thus noted that the analog random walk chain is of finite length in the applications where the Neumann series converges. As mentioned in the previous chapter, we are concerned with evaluating

38 integrals of the form I = W(T) t(T) dT (2-41) T where 4 is a solution of the integral Eq. (2-15). It is desirable to obtain an estimator which when used in conjunction with the random walk process will give an unbiased estimate of I. Suppose I represents the random variable (or estimator) which is evaluated as the random walk chains are processed by Monte Carlo. After N histories of random walk chains are processed, the average of the estimator is given by N I= N1 X I(j) (3-15) j=l where N = total number of histories = total number of random chains The requirement that I be an unbiased estimate of I is just that the expected value of I equals I, E[I) = I (3-16) Let it be assumed that the estimator is evaluated only at the termination point of each chain. In order to identify the chain length with each history a special function,called a characteristic function (or indicator), is introduced. It is defined as Xn(j) = 1, if jth history terminated at Tn;(L=n) 0, otherwise

39 The estimator I can now be defined as k A EI (j = lim Xn(J) In(To''*Tn) (3-17) n=o and the Eq. (3-15) can be written as N k -k I X() In(T on ) (3-18) j=1 n=o A where In can be called the partial estimator corresponding to histories terminating at Tn. For unbiasedness, we must now determine In such that E(I) = E(I) = I (3-19) (The first equality in (3-19) is known from probability theory.37) It is convenient to employ conditional expectations to show unbiasedness. Using the property of the characteristic function, Eq. (3-19) becomes z kk E(I] = Ei l im L Xn(J) A (T ''' k-*00 n n n=o (3-20) k l= im E(InL=n) ~ P(L=n) n=o Because of (3-14) it is noted P(L < } = lim P(L=n = 1 (3-21) n=o The conditional expectations are of the form:

40 E(iL=n) = n h (T,'..TnIL=n) dT...dT (3-22) ' In o n0 f Tn To where hn is the conditional probability density of chain Cn, conditioned on a chain of length exactly n, h (TL=O) =f(To) PO(To) (-2) (3-23) h~(ToL=0) P[L=O) and h (T -.* T nL=n) =fn(To'' ''Tn) | i=O qCi(To*'*'Ti) n(To T n)(3-24) h (T.... n (-4) P(L=n) Using Eqs. (3-22), (3-23), and (3-24) in Eq. (3-20), ~~~~~k ~n-1 A E =IA) lim In fn(To, ' 'Tn) qi( TO,' Ti ) k+ 1-o i=o n=o Tn To (3-25) * Pn(To' **Tn) d'To. —dTn Upon substituting Eqs. (3-7) through (3-11) into (3-25) there results, E(I ) PO( T) d TO k + lim I (T,' 'T ) p (T T ) K(T IT ) (3-26) n=l Tn T0 '''K(T1ITo) S(To) dTo O'dTn Equation (3-26) is now in the form to determine the proper estimator In(To,' ' Tn) SO that I is an unbiased estimate of I. Comparison of Eq.

41 (3-26) to the Neumann series solution of I (2-45), leads to the following choice of the partial estimator In, T n ) W(Tn) t(Tn) (27) In(To,' n) = ToW(Tn) ) (3-27) The division of w by pn can be understood physically by recalling that the estimator is averaged with respect to the termination of the random walk chains which is analogous to capture density, however in order to obtain I from the integral equation, w was averaged with respect to neutron collsion density. The proper derivative estimator for the integral aI/aa (2-50) can be determined in a manner similar to the foregoing analysis for the estimator of the integral I (2-41). Then after referring to the Neumann series solution of 6I/ac (2-51), the proper choice for the partial differential estimator is A F 1 w 1 _ 1O W(TO) DIo(To) wt) + n= W( To) Sa 'i70) ha Po(To) and nA i wn L1, (3-28) DIn(To'''n ) = n + ''' Kn ~ 1 ~ 1 ~~] W(Tn) K1 a S '(To) J p (To a*Tn) where the notations are the same as in Chapter II (2-52) through (2-54). Because of the derivatives involved in estimator (3-28), the derivative

42 estimator is explicitly dependent on previous collision points of the random walk chain whereas the estimator (3-27) is dependent only on the terminal point of the chain for neutron transport phenomena. The estimate of the differentiated integral is given by the average of the partial derivative estimators. = E D~zI(j) (3-29) j=l N k j= k n=o where Xn(j) is the same characteristic function defined previously. By using the Monte Carlo method to generate the individual random walk chains of the analog process and thus averaging the partial estimators (3-27) and (3-29), estimates of the individual terms of the Neumann series solution of the integrals I(2-45) and aI/6a(2-51) can be obtained. B. FORMULATION OF RANDOM WALK CHAINS FOR NON-ANALOG PROCESSES The techniques and mathematical descriptions of the random walk process developed in the previous section serve as basis for this section in which we consider the use of non-analog chains. These random walk chains may be completely arbitrary as long as the condition of finite length with probability one holds. In Chapter II it was discussed that the kernel of the integral equation and of the Neumann series could be modified without mathematically

effecting the solution, if the correct weighting factor was introduced to compensate for the modification. Similarily, non-physical kernel and source terms will be used in this section. If this modified kernel is then used to govern the random walk chain, then a non-analog process would be produced. Let Cn' = o (Tn,'(To Tn)) ( n ( Tn represent a random walk chain of length n that is determined by the arbitrary kernel K'(TIT' ) and arbitrary source S'(T). The source and kernel may be completely arbitrary except that they must have the properties required in Chapter II and in addition determine finite random walk chains. The continuation and termination probabilities and the conditional density function governing the chain Cn' are defined as follows: qn (To, * Tn) - K'(TITn) dT (3-31) Pn (T ~O' Tn) = 1- qn'(T,''T ) (3-32) K' (nlTn n l) fn' (lTnsI Tos ''' aTn (3-33)) q n-i 0 n-i Tahe probability density functions governing the chain C 'are obtained as in Section III-A. fo'(To0) dTo = S'(To) dTo; n=0 (3534)

44 fn '(To '. T )dTo* *dTn = fn (nTon' *** Tn )I ' f 'i( TI T )fo'(To)dTo... dTn n-1 Kt ( Ti+11 Ti) (3-35) 1=0 qK (To0"" i.) S (To)dTo *dT n; n > i-Q 1 0 i Equations (3-33), (3-34), and (3-35) are normalized as can be shown by integrating over all Tn, Toy and To through Tn, respectively, and under the assumption that the source S'(T) is normalized. The random walk process for this non-analog problem is now represented by a sequence of density functions (fn 1(To..Tn)) and a sequence of functions (qn' (To, ''Tn) } and Pn T(To0 ' Tn) ) As in the previous section, the conditional probability density, hnt, conditioned on a chain Cn' of length exactly n will be needed for this non-analog process. It is given by h '(To Ift(TL=) p (T n = 0 (3-36) 0 0 P[L=n) and n-1 fn,(To...*Tn ) i=o qi,(To,...Ti) Pn (T oT. -n) hn'(To...TnlL=n) =- -; n> 1 P(L=n] (3-37) In order to estimate the integral 1(2-41) by the use of non-analog chains processed by Monte Carlo, the proper estimator must first be determined. As in Section III-A, an estimator is acceptable only if it gives an unbiased estimate of the integral I, i.e., EC[)' = I (3-38)

45 whe re k I'(j) = lim L X (j) In (T.n) (3-39) k-*o n n n=o where the characteristic Xn(j) assumes the same definition as in Section III-A and where In' is the partial estimator corresponding to histories terminating at Tn for the non-analog random walk chain Cn' After N histories have been generated by Monte Carlo the estimate of I is given by N I' = N \ '(j) (3-40) j=1 The proper estimator is again determined by the use of conditional probabilities and the density functions governing the non-analog random walk chain, k E([I' = lim j EIn' IL'=n) ~ P(L'=n] (3-41) n=o where the conditional expectations are equal to Elin'IL'=n = n hn (To-'' Tn IL'=n) d0 dTn (3-42) Tn To After using Eqs. (3-31) through (3-37) in Eq. (3-42), and then substituteA ing the result into Eq. (2-41), the expected value of the estimator I' can be written as

46 EAI'] = ' t (T0) po'(To)dT0 T 0 n-1 + lim In (T.. *Tn)Pn (To,''' Tn) 7 K'(Ti+llTi) ken- i=o n=lTn To (3-43) * '(To)dTo '' d Tn 0 0 n The selection of partial estimator I'tn such that the expected value of I' equals I, is obtained by comparing Eq. (3-43) with (2-57), this gives I'A( ) w(Tn)..( T )44) I'n(To,'' Tn) p Wn(To, Tn) (3-44) where the weighting function is given by "(To) Wo(T) = n = 0 (3-45) K(T )ln...K(, -.ITo) 'S(To) Wn(Ton Tn) K 0) n >1 (3-46) KT(ITn K'('TITo) '(To n or expressed in terms of the probability density functions governing the analog and non-analog random walk chains, by n-1 fn(To...T) o qi(To...Tn) Wn(To'..Tn ) = fWn( o' n) F q i(T0 * Tn)(3-47) i=o By using the same techniques as above, the proper choice for the partial derivative estimator in the non-analog process can be obtained. Referring to Eq. (2-60), it is given by

47 I1 WoW 1 W+o j( a ]oJ () W(T)o); n=0 (3-48) I TW(T ) Saa + aa K a s( ( F a___~wn Kn _K1 __ o W(Tn) P (TO' Tn) Wn(To ''Tn); n >1 where the notations are the same as in Chapter II, (2-52) through (2-54), and where the weighting factor is given by (3-45) and (3-46). The estimate of the differentiated integral (2-50) is given by the average of the partial derivative estimators, averaged over N histories of the non-analog process. N DIT = DI'(j) (3-50) j=1 where k DIA(j) = limr Xn(j) DI n(To -.Tn) (3-51) n-o Referring to Eqs. (3-44), (3-48) and (3-49), it is noted that both the partial estimator and partial derivative estimator for the non-analog process are explicitly dependent on all the phase points of the given random walk chain. Furthermore the probability of termination P'n(To*.. Tn) occurs in the denominators since the integrals I and aI/y represent averages with respect to neutron collision density and not with respect to capture density.

CHAPTER IV MONTE CARLO APPROACH In this chapter the Monte Carlo chains are specialized to the neutron slowing down problems in a reactor. Implicit in this case is the fact that neutrons undergo only energy degrading collisions. Usually a neutron is introduced at some high initial energy and followed by random walk process until it is absorbed or its energy has decreased below some cutoff energy, which is above thermal energies. Certain properties can be obtained by averaging reactor parameters over the points of the random walk chain. In the first section a non-analog process is formulated, which forces all collisions to be scattering collisions. It is then necessary to prove that the chain is finite with probability one and to determine the proper estimator that gives an unbiased estimate of the desired property. In order to obtain the proper estimator, the use of expected values is evolved from this non-analog process, i.e., the probability of terminating the chain by absorption is analytically calculated and also the quantity to be averaged is evaluated at each chain point not at just the terminal point. In practice, when the contribution of succeeding chain points to the averaged quantity becomes negligible, the chain is terminated even though the neutron energy may be greater than the cutoff energy. In Section IV-B the non-analog slowing down problem is further specialized to obtaining information about resonance absorption probability and Doppler temperature coefficient. A heterogeneous reactor with infinite 48

49 lattice array is assumed; therefore, a single cell can be considered into which a neutron source is introduced and in which the neutrons undergo perfect reflection at the cell boundaries, i.e., angle of reflection equals angle of incidence. In order to utilize experimentally obtained parameters describing the resonances, the fuel cross sections are expressed as a function of the BreitWigner formula. By neglecting any temperature expansion of fuel materials, all the temperature dependence for this problem is contained in the Doppler broadened cross sections, which are expressable in terms of the Doppler broadening functions. The estimators as well as the transport kernel used in generating the random walk chains are functions of the temperature dependent Doppler broadening functions. In addition the derivative estimators are related to the derivatives with respect to temperature of the Doppler broadening functions. A. SPECIALIZATION TO SLOWING DOWN PROBLEMS IN A REACTOR LATTICE For this type of problem it is convenient to use a high energy source and a non-analog process which forces all collisions to be scattering collisions, but being careful that an unbiased estimate of the desired reactor property is obtained. The chains are terminated only when the neutron energy is less than some arbitrary cutoff energy. Since all collisions are forced to be scattering collisions, the collision kernel in Chapter II is accordingly modified by replacing the ratio Zs/t by one. In this case H is normalized and is given by

50 H(T T1n-1) ( T ITn 1) (4-1) Zs( Tn_ )/Zt(Tn_ l) T(r Irn E 'EnQ ) s( EnnIEn -n r (4-2) Zs(En- 1-nl) The random walk chains are governed by the following probability density function, termination, and continuation probabilities, respectively: dn(To..*'Tn), p*(Tn) q*(Tn) where d(nIT - *..Tn) = H(Tn ln-_ ) (4-3) 0, if En < Ecutoff q*(Tn) = ) (4-4) 1, otherwise p*(Tn) = 1 - qC*(Tn) (4-5) and therefore the probability density function for Tn is given by dn(To Tn) = d(TnIT ''' Tn) dn- (T ITo n-2)..ddo( (To) (4-6) H(TnITn_) '''.H(T1ITo) S'(TO) For convenience the source is assumed normalized in this section. It is noted that the values of p* and q* do not depend on the sample sequence but rather are determined by the given phase point in the chain. As mentioned in Chapter III, the non-analog chains may be completely

51 arbitrary, provided the chain has probability one of being finite. Therefore, as in Section III-A, it is necessary that the probability of generating a chain of infinite length is zero. This probability is given by P(L=o)] = lim dm(TO...Tm ) T q*(Ti) dT'~ *dTi mmo Tm T i=o (4-7) m-1 im q*(T) "(T) f H(Ti+llTi) q*(Ti) dT0o mdT mo*0m i=o T To Because of the nature of q* in this non-analog problem for a heterogeneous lattice, it is difficult to simplify this expression and show that the limit is zero directly. However it is shown in Appendix D by physical arguments that the above process is finite with probability one. Since elastic collisions in the epi-thermal region are energy degrading, it seems reasonable that the neutron energy will decrease below the cutoff energy in a finite number of collisions. By combining the use of expected values with this non-analog chain, it is possible to obtain more information from a given chain and hence a more efficient determination of the desired reactor property. Briefly, the use of expected values involves analytically calculating the expected value of a random variable rather than using a sample value. For the problem of interest the use of expected values is applied by calculating at each chain point the true probability of absorption, Za/Zt, and multiplying this by the quantity to be averaged, w, at each chain point. The estimator for this problem will now be determinedo

52 The random walk process governed by Eqs. (4-1) through (4-6) has the following conditional density function, conditioned on the chain length exactly n. n-1 bn(To,...TnIL=n) _ dn(To,... Tn) P*(Tn) 1=o q*(Ti) (4-8) P(L=n] n-1 p*(Tn) gS(To) i=0o H(Ti+1Ti) q*(Ti) (4-9) P[L=n) But the termination of the chain is determined only by the energy of the last point, i.e., the chain terminates only when a phase point has energy less than the cutoff energy. In the previous chapter the estimator contributed to the knowledge of the integral 1(2-41) only at the termination point of the chain, although it was allowed to depend on all previous points. However since each chain point can contribute to the estimate of I in this chapter, the estimator is written N k I = lim Xn(j) In(T'' Tn) (4-10) N k+o n j=l n=o where * n In(To~ Tn. ) = I*m( To ~ *Tm) (4-11) m=o *Spanier28 has shown that by using the integral equation adjoint to (2-15) and (2-41), a random variable determined by summing over all the phase points in the chain can be an unbiased random variable for an analog random walk process. In the investigation of slowing down problems, we obtain the proper estimator for I(2-41) for a non-analog process without relying on the adjoint equation.

55 As in Chapter III k I(j) = lim Xn(j) In(To' Tn) (4-12) k+o n=o where Xn(j) is the same as defined in Chapter III. According to the reA quirement of unbiasedness, In must be determined such that E(I = I (4-13) Employing conditional expectations with the non-analog process, n E(I) = E Xxn(i) I*m(TO,' 'Tm) (4-14) IE * I( m(To. *.Tm)IL=n| PtL=n) (4-15) n=o =o But by Eq. (4-9), E I*m(To ~*Tm)IL=n (m=o n = ' (ZI'm( To,' Tm)) bn( To, ~~TnL=n)dTo*' dTn (4-16) Tn To m=o n-1 p*(Tn) S(TO) i=o (Ti+Ti) (Ti) I*M(TO, ~Tm) -.. H(..I. dTo -'*dT m=o T T n o The next step is to substitute E4s. (4-16) into (4-15) and employ the Cauchy product rule, which is given by

54 00 n o0 00 00 00 1 F(n,m) = ' 7F(n+m,m) = F(n,m) n=o m=o n=o m=o m=o n=m Equation (4-15) becomes A E (I (4-17) 00 n n-1 X o / J moI* (T m. nP() () iT H(T lTi)q*(Ti)dTo * *dT n=o m=oTn To 1=0 00oo CO '. same integrand (4-18) m=o n=m Lj I'*m( T,' m) -[ H( Ti+l Ti) q*(Ti) S(To) dT dT m=o Tm To (4-19) where n-1 = p*(T + p*(Tn) H(Ti+lTi)q*(Ti) dTm+1 ndT n=m+1 Tn Tm+ i= (4-20) but since p*(T ) = - q*(T) and H(T i+lTi) dT i+ T. 1+1 then

n-1 T. p n*(T) [ H(T. jTi.) q*(Ti) dT * dT 1+1 1 1 m+l n Tn Tm+l n-1 = *( Tm ) f q*( i) H(Ti I Ti_ ) dTm+l..dTn 1 (4-21) T T+ =m+ 1 n-1 m+i n q*( T) ~ q*( T i) H(Ti I Ti+) dTm+ *T m m+ 1 n Tn Tm+l 1 =m+i Upon using Eq. (4-21), (4-20) becomes + q *(Tm) 1 q*( Tn) H(TnTn) dTn (4-22) Tm++1 q*(T m- q*(T n ) H(T nI-) (4-22) n=m+2 Tn- Tm+l Tn n-1 q*(T.i ) H(TiT.i) dTM+ n1 11-1 m+i n-i After cancellation of like terms, { } = 1 - lim q*(Tm) q*(Ti) H(Ti Ti 1) dm+l...dT n n m+ i=m+l (4-23) = 1 (4-24) Upon comparing the integral term in (4-23) to Eq. (4-7), the limit is zero

56 under the conditions of non-zero source and finiteness of the random walk chains. Equation (4-15) now becomes 00 = j.. I*'m( To)*** Tm) H(Ti+Tij)qg*(Ti)S((T0)dTo'.'dTm i=o m=o Tm T (-25) (4-25) I*O(To) S(T) dT T 0 (4-26) 00 + I*m( To ' ' ) 'Tm H(Ti+llTi) S(To) dTo *dT m=1 T T 10 m o The q*(Ti) are all equal to one since the integral was conditioned on a chain of length given by the index. Comparison of this series to the Neumann series for I in Chapter II (2-57) leads to the following estimator which satisfies the condition of unbiasedness: I (T...) = W(Tm) Wm_(1(To''T ) (4-27) where Zs(TM) ZS(To0) Wm(To''' *Tm) (T t(TO) (4-28) Therefore N 0 n Xn( )I*m(ToT) (4-29) j=i n=o m=o

57 is the estimate value of I(2-41) as determined by a Monte Carlo process generating non-analog random chains which terminate in a reactor only when the corresponding energy of the chain becomes less than some cutoff energy and in which each collision point of the chain contributes knowledge to the average value. In practice the chain may actually be terminated when the contributing knowledge becomes negligible (when the product of the ratios of Es/At becomes << 1). Because of the uniform convergence of the Neumann series representation of I(2-45) and aI/6a(2-50), and the finiteness of the non-analog random walk chains in this chapter, the partial derivative estimator corresponding to aI/aC is determinable. By combining the methods of Chapter III and the preceding analysis of this chapter and referring to Eq. (2-60), the corresponding partial derivative estimator is found to be DI*m(To,... Tm) = W(Tm) Wm-_( T70.'Tm-_) DWm(To).. Tm) (4-30) where 1 6Wo 1 6~o DWo(To) 1 w 1S( (4-31) w(T) T c T() ac and for m > 1 DW(1 Tm 1 wm 1 Km 1 S ~DW m(To ''Tm -w aCl K 601 K &C S C(4-2) m m 01 and where again the following notations are used: w m (m) (4-33) Km - K(TmlTm-_) (4-4)

58 SO~ o(~~ O~ ~(4-35) The derivatives are evaluated at the chain points. Thence N cO n n N = ti ZZ X (j) DI*m(T~, -T' (4-36) j=l n=o m=o is the Monte Carlo estimate of aI/8a (2-50). B. SPECIALIZATION TO DOPPLER COEFFICIENT AND RESONANCE ESCAPE PROBABILITY In this section the background provided by the previous sections are applied to the specific problem of resonance absorption probability and temperature coefficient in a heterogeneous reactor having a non-uniform temperature in the fuel rod. In order to determine the absorption probability, the random variable w(T) is set equal to the ratio of neutron absorption cross section to the total cross section. Recalling that * represents the collision density, then I = w(T) *V(T) dT = ZTa(.(T) dT (4-37) R R Zt(T) represents the number of neutron absorptions per unit time in region R. The resonance absorption probability can be thought of as the ratio of the number of absorptions per unit time in the phase space with energy greater than the thermal cutoff to the number of neutron births per unit time in the same space. For numerical results this definition will be specialized

59 to the case of no absorptions in non-fuel regions. By using the nomenclature p = Resonance escape probability 1-p = Resonance absorption probability then j ', *(T) dT p = (4-38) S(T) dT R where R = region of the whole space T with E > E cutoff = number of source neutrons per unit time (or birth rate) = S(T) dT /S( dT R. R The derivative with respect to temperature of the resonance absorptioh probability is given by 6T S 66T- S R That part of the temperature coefficient of reactivity due to absorptions of neutrons in the resonances of the fuel will be called the Doppler coefficient of reactivity, and is then given by 1 6p -1 aI!..~p~ _ = -l~ ~~ ~(4-4o) p aT s - I aT R -~T(Z~~~) dT (4-41) S dT - dT R R

60 Since this work is concerned with the problem of temperature effects due to the Doppler broadening of resonances in the fuel, any effects due to changes in geometry or material density with temperature are purposely neglected. Specifically, U-238 is considered in oxide form (U02) because of its use in power reactors, although the methods proposed are valid for the metal form. In order to exhibit the temperature dependence of the fuel cross section explicitly, two assumptions are made, (1) the motion of the fuel atoms which are bound in a crystal are described by the perfect gas model, i.e., their velocity distribution is approximated by the Maxwell-Boltzmann distribution for free particles in aidealgas, and (2) the resonances of the fuel are spaced far enough apart that there is no interaction among neighboring levels and therefore the single-level BreitWigner formula2,38 can be used to describe the microscopic cross section in the center of mass coordinate system. Lamb24 studied the resonance capture of slow neutrons by atoms bound in a crystal and developed a criterion for using the ideal gas model in calculating the resonance cross section. The criterion depends on the Debye temperature of the crystal, the natural width of the resonance, the mass of the absorber nucleus, the neutron energy, and the environmental temperature. According to Lamb's criterion the gas model is a good approximation for the lowest resonance of U-238 (6.68ev~) when the absorbing atoms are bound in the metallic crystal at room temperature, and the approximation improves in accuracy with increase in temperature or with the higher resonance levels. The uncertainty of the Debye temperature for uranium

61 dioxide does not allow reliable application of Lamb's criterion. The large discrepancy among measurements and calculations of the Debye temperature of uranium dioxide has been noted by Belle21 who reported estimates varying from 160~K to 800~K. He does indicate that there is more evidence for correctness of lower values, which would fulfill Lamb's criterion. In 1960 Jackson, Bollinger, and Cot359 performed an experiment for determining the shape of the effective neutron cross section over the 6.68 ev resonance peak of U-238. They observed no significant crystalline effect in uranium metal or U308 at room temperature. However at liquid nitrogen temperature a large discrepancy was noted between the observed shape of the cross section in U308 and that predicted by the gas model corresponding to the effective temperature given by Lamb. The authors used a Debye temperature of about 4500K for U308 which corresponds to about the mean in the Debye temperature range reported by Belle for U02. Since reactors used for power purposes have fuel temperatures considerably above room temperature, it appears that the gas model approximation is acceptable for determination of the U-238 resonance cross section in both uranium metal and U02o The second assumption is substantiated by the experimental work of Rosen, Dejardins, Rainwater, and Havens40 in which they were able to resolve the individual resonances of U-238 in metal foils up to about 1000 ev. All but a few of the resonances were observed to be clearly separable. Furthermore the resonances were sufficiently isolated that resonance parameters could be determined from the experimental data by a theoretical

62 analysis using the single level Breit-Wigner formula and the Maxwell-Boltzmann distribution. In fact the authors use the same expressions for the Doppler broadened cross sections as those given by Eqs. (4-42) through (4-5o). The derivation of the Doppler broadened cross sections, using the single-level Breit-Wigner formula to describe the center of mass interaction probability and using the Maxwell-Boltzmann distribution to describe the motion of the struck nuclei, has been extensively investigated. Bethe and Placzek3'41 and Dresner42 make a basic approximation in their method, namely that the term involving the square of the ratio of atom speed divided by the neutron speed can be neglected in determining the kinetic energy of the neutron relative to the struck atom. This approximation is known as the Bethe approximation. Solbrig,43 and Osborn (footnote in reference 44) both derive a more exact expression for the cross section, in which Solbrig later introduces approximations which he shows to be valid for incident neutron energies above.01ev and for heavy absorbers at temperatures not less than room temperature. Nordheim9 has also investigated this latter method. For the conditions present in this work, i.e., incident neutron and resonance energies above thermal energies, and fuel temperatures above room temperature, both methods produce the same expressions for the Doppler broadened neutron resonance cross sections. Since the experimentally determined resonance parameters are obtained from measurements in the laboratory coordinate system, the parameters in the usual expressions of Dop

63 pier broadened cross sections should be appropriately modified by the ratio of neutron mass to reduced mass to take account of this fact. However for fuel nuclei this correction amounts to less than j? and will accordingly be neglected. Therefore all parameters in the following expressions for cross sections represent values in the laboratory coordinate system. The Doppler broadened microscopic cross section for the capture of neutrons in the neighborhood of an isolated resonance is ac(E,T) = a() + aco /(xe) (4-42) The Doppler broadened microscopic cross section for the elastic scattering of neutrons in the neighborhood of an isolated resonance is Cs(ET) p + a sO Y (x,@) + J7 sX (x,o) (4-43) 1 1 where a(v) and p represent the background "-" absorption cross section and potential scattering cross section respectively and where acO = aO r (4-44) co o T aso (4-45) so 0 The unbroadened cross section for the formation of the compound nucleus, with neutron energy corresponding to the resonance peak, is 0 _ 4~keg 0_r (4-46) ~o F

64 where E = energy of neutron in laboratory system. Eo = energy of neutron corresponding to resonance peak. t = h2, kois the neutron wave length at resonance peak. 2J+l g 2(2I+1) I,J are the spins of target and compound nucleus respectively. r = rno+ry = total resonance width, full width at half the maximum peak. rno = neutron width with E = Eo. rF = capture width. The Doppler broadening functions, also called spectroscopic integrals because of their original use in optics, are defined as 1 - (x-y)2 /')/ (x,a) J e 49 (4?47) 1 + y2 + Y2 (X-y)2 00 4E t (x,@) _ _1 1 e 2 (4-48) so_ i001 + y where E - Eo x = (4-49) 2 = 4mkTE (m + M) F2 m = mass of neutron M = mass of target nucleus T = Temperature of resonance medium

65 It is noted that the temperature dependence of the microscopic cross sections appear solely in the Doppler broadening functions. The potential scattering cross section is temperature independent for neutron energies above thermal. In accordance with the restriction that there is no interaction between neighboring resonance levels and with the approximations discussed above, the total macroscopic cross section for a neutron of energy E is given by a simple sum of contributions from all the resonances, i.e., Tt(ET) N A p + r(l/v) + cr +- ) 41 (xiOi) + a G5 ( xiqi)i; (4-51) where i represents the index of the resonance levels and NA is the atom density of the resonance absorber. For the sake of convenience in the mathematical description of the estimators and kernels in the random walk process, the cross sections will be written explicitly for only one resonance, the nearest dominating resonance, although the effect of the summation of resonances will be implied. In practice, for a given neutron energy the one resonance of U-238 or Th-232 with Eo nearest the neutron energy is dominant and all the rest of the resonances make a negligible contribution to the total cross section. In the case where the fuel is a mixture of resonance atoms and nonresonance atoms, e.g., UO2, let oa(l/v) and cp be respectively the back

66 ground absorption and potential scattering cross sections per atom of resonance absorber. Also let the total cross section in a reactor lattice composed of separate fuel and moderator regions be represented by the following notation: It(T) = NA(OA + cA (l/v)) + ~a | V(x,@) (4-52) + y (x, ) + i (X, ) (T) + M(oM + (/V)) V(T) rs ri (C +Ca v where Zo A (AA) ra = N co (4-53) O e o aA (4-54) rs so Fri NAgop a A (4s55) where the superscripts M and A denote the moderator and the resonance material respectively, and where B and v are heterogeneity factors defined as ( 1, if the neutron position r lies in the fuel. ( T) 0, otherwise. ( 1, if the neutron position r lies in the moderator. V(T) 0, otherwise. Recall that T is the phase space vector of the neutron and thus denotes (rEn ). Also it is remarked that the temperature T is a function of r, and hence of T, since a general temperature distribution is allowed in the fuel. Thence the derivatives with respect to temperature of the cross sections are

a 4(T) -Za (T) + - S(T) (4-56) T a T (4-57) Z5(T)- - (=-ps ) aS S(T) L [-AS ai aTI( (4-58) Methods of evaluating the derivative of the Doppler broadening functions are discussed in Appendix F. One further approximation is necessary before writing explicitly the kernels governing the random walk chain. Once a scattering collision occurs, the scattering angle and new neutron energy must be selected from the collision kernel. When the neutron collides with a non-resonance nucleus there is no difficulty since the well known scattering frequency (4-60) for elastic collisions with struck particle at rest can be used for the neutron energies of interest, i.e., above thermal and below inelastic threshold energies. For collisions with resonance nuclei, this same frequency (4-60) is used for the scattering distribution. The assumption of struck particle at rest is harder to justify for collisions involving resonance nuclei and is investigated in Appendix E, where the resonance scattering frequency is formulated for the case of the target nucleus in motion, (E-8). Although the resonance cross section is effected by the motion of the struck particles, once a collision occurs and the decision is made that it is of a scattering type, then the dynamics of the scattering are similar to that in which the target is at rest since

68 the incident velocity of the neutron is much greater than that of the nucleus. The collision kernels defined in Chapter II can now be written, EZs(r' El) Id& C(E,_IE',_';r') dEd = F(E'+E;r') dE 5(Q ~ ' - g(E,E')) (4-59) Also F(E'+E,r') becomes F(E'+E), once a decision is made as to the kind of nucleus involved in the scattering collision. The scattering frequency for an elastic collision with the struck nucleus of mass number A at rest is F(E'+E) dE when ( + d E' < E < E' 4A E' A + 1 - - (4-60) =; otherwise Also g(E,E') A + 1 EA - 1 (4-61) 2 El' 2 E By utilizing the definition of the transport and collision kernels in Chapter II and Eqs. (4-59) and (4-60), the kernel H (4-5) for the nonanalog process described in Section IV-B becomes - >Zt(E, r-QRR')dR' t(E,r) H(T ' T') dT —e O 62(s. R) - d3r ' - (4 126 (4-62) F(E'+E,r') dE b(_'Q' - g(E,E')) 4 Thence the kernel K defined in Chapter II becomes

69 -J Z t(E,r-2RR') dR' Zt(E,r) K( IT')dT = e- r2 *2(2-'2-R)!r - r' -2 (4-63) s(' E')Et) Z_ ( F(E'+E,r') dE 5(2s.' - g(E,E')) -- ~t(j' E') -- 4| The derivative with respect to temperature of the kernel K is then -Kr(TIT' it) t aZt (T) dR 6_T - t(T) (T (4-64) + 1 — ES (T') - m i a (T') J K( TIT') In practice the initial phase point of the random walk chains are selected from the conventional 'birth' source S. rather than the collision source S, and the position of the first collision is selected from the transport kernel. Since the energy of the source is established above the resonance energy region, the cross sections are temperature independent-at source energy. Therefore the derivative of the source S with respect to temperature appearing in (4-31) and (4-32) is set equal to zero. With the random variable w defined as Z'(T) w(T) - t_( (4-65) the problem of estimating the resonance capture probability and Doppler coefficient in a heterogeneous lattice can now be solved by using the Monte Carlo method to process the random walk chains described in Sec

70 tion IV-B. Referring to Eqs. (4-59) through (4-64), the partial estimator (4-27) and partial derivative estimator (4-30) can be explicitly described as follows: I*rm(To.'Trn) - W (T **Tr ) (4-66) Z(Tm) and Za(Tm) DI* (TO,, *TM) Z (Tm) Wm-_ 1 (To.T) DWm(To...Tm) (4-67) where Zs(Tm) s( TO) W (T...T (4-68) Wi( ~ "'' Tm) Zt (Tm) Ft (TO) and DWo(To) = 0 (4-69) and for m > 1 DWm(To~'' TM) Za(Tm) 1(Tm) dR Za(Tm) aT o T (4-70) + Km- +... '' - 1 where Za Z (Tm) — J a- Zt(T) dR (4-71) ms (Tin) 6T s and Rm r I~m ~ xmn-iI (4-72)

It is remarked that derivatives with respect to temperature of the total cross sections do not appear outside the integrals in (4-70) even though w is defined as (4-65). This is because all factors of Z(T) cancel from the product of w(T ) K...K S m m 1 0 Estimates by the Monte Carlo method of capture probability and Doppler coefficient of reactivity are given by the following: aP- n - I (4-73) and i ~ -DI (4-74) where N N =) %,A Xn(J) I*m(mo' rm) (4-75) j=l n=o m=o and N oo n DI = XXn(j) DI*m(To,' Tm) (4-76) j=l n=o m=o

CHAPTER V MONTE CARLO PROCEDURE AND COMPUTER PROGRAM A. DESCRIPTION OF MONTE CARLO PROCEDURE The procedure for processing Monte Carlo histories and for evaluating the Doppler coefficient and the resonance capture probability on a digital computer is described in this chapter. The procedure utilizes the mathematical theory developed in the preceding chapters. A semi-infinite lattice cell with a heterogeneous mixture of fuel and moderator is considered for this problem. The fuel is contained in a cylinder which has its axis colinear with that of the cell. Stated in brief, Monte Carlo is used to generate a history by tracing neutron propagations through the lattice cell as the neutron is degraded in energy by scattering collisions from source energy to cutoff energy. The value of the cutoff energy is preset above thermal energies and below the energy of the fuel resonances. Decisions regarding the behavior of neutrons in a reactor lattice are determined by random variables selected from the probability density functions governing the random walk chains of Chapter IV. Data is collected as the neutron histories are generated and is used to estimate values of the desired reactor parameters, i.e., resonance capture probability and Doppler coefficient of reactivity (Eqs. 4-73 and 4-74). The application of sequences of pseudo-random numbers to select random variables from prob72

ability density functions which depict physical processes has been discussed by many authors.45,..48 The assumption that the pseudo-random numbers are uniformly and independently distributed on the interval between zero and one is a very good one for numbers generated by the Method of Congruences,49'''.53 which is also used here. The mathematical details regarding the selection of random variables from the appropriate distribution functions are given in Appendix G. The verbal description of the Ybnte Carlo procedure for this problem follows. 1. Each history is started by selecting a pseudo-random number from the first random number generator (G-l), which is used to initiate the sequence of pseudo-random numbers generated by the second random number generator (G-2). This latter sequence of pseudo-random numbers is used to determine the behavior of the neutron for the given history. Since the number of random numbers used in each history depends upon its length and the paths taken in the cell, the use of two random number generators allows one to know beforehand the initial random number of each history. 2. The initial energy, position, and direction of the neutron is selected from the source distribution functions, Appendix G-3. The neutron source is taken to be isotropic in direction, uniformly distributed in the cell, and a line source in energy. The initial weight of the neutron is set equal to one (4-68). 3. Upon knowing the neutron energy, direction, and position at the source point (or last collision point), the distance, measured in terms of optical thickness, to the next collision point is selected from the

74 exponential distribution law of attenuation, Appendix G-4. The actual position of the next collision is determined by knowledge of the total macroscopic cross section, taking into account the crossing of boundaries separating dissimilar media, the variation of the resonance cross section due to the temperature distribution in the fuel, and the reflection of the neutron path on the outer cell boundary. At the outer cell boundary, the angle of reflection is set equal to the angle of incidence. 40 Upon knowing the position in the cell at which the collision takes place, the species of the interacting nucleus is selected by random variables according to the contribution to the total cross section at the phase point of the collision. A scattering type collision is then forced to occur by virtue of the modified kernel H (4-62) in Chapter IV and the weight of the neutron is multiplied by the ratio of scattering to total cross section (4-68). The derivatives with respect to temperature of the Doppler broadening functions are integrated over the neutron path in the fuel, and the partial estimators for p and ap/)T are calculated by Eqs. (4-66) and (4-67). If the neutron weight is less than a preset constant, say.001,* the history is terminated and the procedure is continued as in paragraph 6, otherwise the history is continued as in paragraph 5. 5. The scattering angle in the C.M.C.S. is selected from an isotropic *The use of Russian roulette for deciding termination of a history was provided for in the program and probably should be used if the input constant is made considerably larger than.001. Russian roulette was tried but did not produce an advantage in time or statistics over the straightforward use of the preset constant.001.

75 distribution, Appendix G-5. The new direction in the L.C.S. is determined by a transformation of coordinate systems, and the neutron energy after the scattering collision is calculated, according to the species of interacting nucleus, Appendix G-5. If the neutron energy is now greater than the preset thermal cutoff energy the history is continued by again following the procedure in paragraph 3, otherwise the history is terminated and the procedure in paragraph 6, is followed. 6. Once a decision has been made to terminate the history, either by conditions in paragraphs 4 or 5, the contributions of the history to p and ap/6T are added to the sums in Eqs. (4-75) and (4-76). If the total number of histories processed is less than a preset maximum, another history is processed by again starting with the procedure in paragraph 1. Otherwise the Monte Carlo estimate is obtained for the resonance capture probability and Doppler coefficient of reactivity by evaluating Eqs. (4-73) through (4-76). In practice the histories are processed and the above estimates obtained for subgroups of histories so that trends can be observed and statistical variations noted. Also estimates of the variance of capture probability and its derivative are obtained by first calculating the average values of the square of the partial estimators (4-66) and (4-67). B. BRIEF DESCRIPTION OF PROGRAM Since part of the purpose of this work is to provide a computer program which utilizes the theory presented in the preceding chapters and

76 which is capable of handling the case of a non-uniform temperature distribution in the fuel, a brief description of the program is appropriate here. Most of this program was written in the MAD language (Michigan Algorithm Decoder54) for the IBM 709 digital computer. The Random subroutine discussed in paragraph 6 was written in symbolic machine language. The program requires a machine with 32,768 words in fast access memory. Because of the large size of the program, it was necessary to break it into component parts that could be separately translated by the MAD compiler. The division of the program according to separate functions has facilitated debugging and will make future modifications discussed in the succeeding chapters much easier. The program, called REPAD (Resonance Escape Probability and its Derivative with respect to temperature) includes the following routines: 1. The Main routine generates neutron histories and follows them in energy degradation and spatial propagation through the fuel and moderator regions of the reactor lattice. When the option of a non-uniform temperature distribution is selected, the following routine is used for the neutron propagation through the fuel. 2. The Non-uniform Temperature subroutine computes the optical thickness of the neutron path (and its derivative with respect to temperature) through the temperature zones of the fuel, using cross sections obtained from the Cross Section subroutine, and determines the next collision point. If the neutron energy is such that resonances need not be taken into account, this subroutine is bypassed.

77 3. The preparation subroutine prepares input data, cross sections of resolved resonances evaluated at the peak energy, and macroscopic crosssections for the program. The radial temperature distribution can be completely arbitrary, however if the parabolic temperature distribution is selected, this subroutine calculates the value of the temperature in each of the desired temperature zones. 4. The Cross Section subroutine determines the appropriate macroscopic scattering, absorption, and total cross sections in the fuel for the current value of neutron energy. In the neutron energy region of unresolved resonances, the neutron width is selected from the Porter-Thomas distribution,40 and the unbroadened resonance scattering, absorption, and interference scattering cross sections are calculated. In the region of resolved resonances the decision is made whether the neutron energy is in a limited range between resonances where the resonance cross sections are negligible with respect to the potential scattering cross section or not. If not, the unbroadened cross sections are determined for the appropriate resonance or resonances. The Doppler broadening of resonances is accounted for in the next subroutine. The "l/v" absorption cross section is calculated if it is important with respect to potential scattering cross section for the given neutron energy. 5. The Psi-Chi subroutine evaluates the Doppler broadening functions for the appropriate resonance parameters and fuel temperature to a prescribed accuracy, and calculates their derivatives with respect to temperature, Mathematical details of calculating these functions are given in Appendix F.

6. The Random subroutine generates pseudo-random numbers and variables used in the problem. Selection of random variables from uniform, exponential, polar cosine, azimuthal sine and cosine, and parabolic distributions are required. 7. The Random Input/Output subroutine accepts the initial random numbers for input and prepares initial, final, and intermediate random numbers for output. This routine is necessary for running a series of many histories on the digital computer with the same physical conditions but different random numbers. 8. The Edit subroutine computes averages, and estimated variances and standard errors of the reactor lattice parameters for subgroups of histories and for the completed problem. The parameters include resonance capture probability and its derivative with respect to temperature, resonance escape probability and Doppler coefficient of reactivity. The first two parameters are also given for three separate energy groups. This subroutine also punches output cards to be used as input for succeeding problems. The program REPAD is written for a cylindrical fuel rod and lattice cell, both of arbitrary dimensions. Both cell and fuel rod are circular and have colinear axes.

CHAPTER VI RESULTS A. PROGRAM CONSTANTS AND INPUT PARAMETERS For the total time available on The University of Michigan IBM 709 computer for this problem, it was decided that for purposes of comparison it would be better to run only a few cases with 12,000 Monte Carlo histories rather than more cases with fewer histories and larger statistical error. The production runs were for a.25 in. diameter rod of U23802 fuel and a 1/1 fuel to moderator volume ratio. These and other input parameters, such as fuel surface temperatures and non-resonance cross sections, were selected to be compatible with those appearing in practice and consistent with those in WCAP-1572 in order to make a comparison with the results of Arnold and Dannels15 for the case of uniform fuel temperature. Values of input parameters that are common for all the runs are given in Table I. The recent resonance parameters determined at Columbia University by Rosen, Dejardins, Rainwater and Havens40 are used for U-238. Since the resonance parameters used in REP16 were older measurements which included only 28 resolved resonances, this number of resonances was also used in the production runs with REPAD. Unresolved resonances were not included because a considerable savings of computer time could be realized by omitting the unresolved resonances, and because it seemed incongruous to select unresolved resonances parameters from a probability distribution when actually another 27 resonances are 79

TABLE I VALUES OF INPUT PARAMETERS FOR REPAD PROGRAM Parameters Value Description vf/vm 1.000 Fuel-Moderator volume ratio Rf e3175 cm Radius of fuel rod (.25 in.) Rm.4490 cm Radius of moderator cell.416 cm-l Potential scattering cross section of fuel-U02 1.468 cm-' Potential scattering cross section of moderator-H20 o 1.468 cm-l Total cross section of moderator-H20.0522 cm-1 "l/v" absorption cross section of fuel at 1 ev ERN(1) 100,000 ev Reference value - source energy ERN(2) 5530 ev Reference value - not used ERN(3) 545.75 ev Reference value - max energy of resolved resonance region ERN(4) 250.00 ev Reference value - cutoff energy for "l/v" cross section ERN(5).625 ev Reference value - thermal cutoff energy

81 now resolved. According to the results presented in WCAP-1572 for a uniform fuel temperature, the 28 resonances of U-238 below 545 ev, which also include the strong resonances, contribute more than 80% to both the resonance capture probability and the Doppler coefficient; it is, therefore especially important to investigate the effects of a non-uniform temperature distribution due to these resonances. Furthermore, for fuel temperatures above 800'K, unresolved resonances corresponding to the average reduced neutron width (< Fn > <= <no/ Eo>) are broadened sufficiently so that at peak energies the average neutron path through the fuel rod is of the order of 1 mean free path, and consequently neutrons incident on the fuel rod have a good chance of passing through most of the temperature regions in the fuel. It is expected therefore that the concept of an average fuel temperature would have more meaning for the unresolved resonances than the resolved resonances. B. NUMERICAL RESULTS AND DISCUSSION The running time of REPAD is about.01 min/history for uniform fuel temperature and.02 min/history for non-uniform fuel temperature on the IBM 709. The latter time is for 24 temperature regions in the fuel. Runs of 12,000 histories each were made for the following conditions: Uniform fuel temperature Run Ts,~K 1 84o 2 1700 3 2800

82 Parabolic fuel temperature Run Ts,~K Tp,~OK Ts/Tp (radially) 4 840 2800 3.34 5 1820 3780 2.08 where Ts = surface temperature Tp = peak temperature at centerline Also let Ts + Tp Tav 2 Trm = ( XFTs +3 2 As mentioned previously one of the features of IREPAD is that the derivative with respect to temperature of the capture probability is directly evaluated by Monte Carlo. This then allows calculation of that part of the temperature coefficient of reactivity due to the resonance capture probability and also of the Doppler coefficient. Given these results one may determine D in the empirical T1/2 law for resonance integrals,10,17 RI(T) = RI(To) ~ (1 + P (AfT = fWTo)) (6-1) From Eqo (1-2), and Eqo (6-1), f can be expressed as - ap/6T 2fT (6-2) p. in (p) The results of the Monte Carlo runs are given in Table II and are summarized in Figs. 1 through 3. (The error indications shown in the figures are for probable error.) Results from WCAP-1572 are also included

TABLE II RESULTS FROM 12,000 MONTE CARLO HISTORIES Flat Temperature Distribution, sK Parabolic Temperature Distribution, K _Ls Ts Ts Ts Tp Ts T 840 1700 2800 840 2800 1820 3780 (l-p).21329.22777.24180.22533.24094 +.6745 a +.00207 ~.00213 +.00218 +.00211 +.00224 -ap/l T (10-4).2134.1666.1355.1796.1348 +~6745 a (10-4) ~. 185 ~.0106 ~.0073 +.0135 +.oo85 P.78671.77223.75820.77467.75906 -En p.23990.25848.27681.25532.27587 -2n( -n p) 1.428 1.352 1.283 1.365 1.289 -1 -p ( 10-4).271.216.179.232.178 p 6T 1 a" (10-) 1.131.834.646.909.644 penp aT (10-2) ~655.687.683.744.671 T' =.675 x 1O.2.030 X 10-2

84.30 1.29.28 -.27 DATA FROM WCAP 1/572 5/25 H/STORIES — ~26~~~~ 4~~ REPAD ~.26 1 /iO/2,000 H/STORIES c.25 / RUN 5.24 52.1 (2730 OK).23 r 7 ff /RUN 42.65.22 ( 820 t soOK) v./ 39.0 (/520 "K) Treff.2 i (300~OK) (840K) (1/700~K) (2800~K) (6000 0A /1732 2898 41.23 52.92 773.20 I I I i I L I I 10 20 30 40 50 60 70 80 1/2 1/2 T =(TEMPERATURE, K) Fig. 1. The variation o' resonance capture probability with fuel temperature; uniform temperature case, runs 1, 2, 3. Values for parabolic temperature distribution indicated on lower curve, runs 4, 5..25 in. dia. U02 fuel, 28 resolved resonances.

.28 REPAD.26 1{ /2,000 H/STORIES.26.24 -RUN RIUN 4~.22 -.20- (/400S) ';ff 8 -/ p (-p/dT) RUN 4 - T RUN 5 o.16 374 (14000K) 7e.14 Teff -dp/aT (/820 0) RUN.12 -To0 (Run 4) Teff. 10.8 -.6.4.2' (840 OK) (/700~K) (2800~K) 28.98 41.23 52.92 0OI II 1i I I l 10 20 30 40 50 60 70 80 T 1/2= (TEMPERATURE, OK )i/2 Fig. 2. The variation of the Doppler coefficient of reactivity and the derivative of resonance escape probability with fuel temperature; uniform temperature case, runs 1, 2, 3. Values for parabolic temperature distribution indicated on curves, runs 4, 5..25 in. dia. U02 fuel, 28 resolved resonances.

86 1.6 REPAD 1/2,000 H/STORIES 1.5- t DATA FROM WCAP /572 5/25 HIS TORIES SLOPE CORRESPONDING 1.4\ TO, =.675 x /O 1.30 \wSL OPE B =.6/x/0 1.2 2 SLOPE = =.66 x /0 - I.! 1.0.9-.8 (300 OK) (8 40 K) (1/700 OK) (2200 OK) (6000 oK /732 28.98 41.23 52.92 773.? II I I I I I I 10 20 30 40 50 60 70 80 T /2 (TEMPERATURE, OK) )/2 Fig. 3. The variation of ln(-lnp) with fuel temperature; uniform temperature case, runs 1, 2, 3..25 in. dia. U02 fuel, 28 resolved resonances.

87 for comparison. Values of capture probability were obtained by using data in WCAP-1572 from the resolved resonance energy region only. The differences between the results of WCAP-1572 and REPAD for flat fuel temperature are attributed to the difference in outside cell boundaries, i.e., rectangular and circular, respectively. The average value of P obtained by Eq. (6-2) from the three flat fuel temperature runs of REPAD was found to be.675 x 10-2, which is close to the value obtained from WCAP-1572 data (see Fig. 3). By representing the slope corresponding to P = 0675 x 10-2 by the dashed line superimposed on the upper curve in Fig. 3, one notices that the dashed line lies within the error limits, and thus helps to substantiate relations (6-1) and (6-2). Figure 2 illustrates the dependence upon temperature of ap/6T, the derivative of capture probability, and(l/p). ap/aT, the temperature coefficient due to Doppler effect. The curves are for a flat temperature in the fuel. It is noted that both values decrease about 35% as the temperature increases from 840 to 2800~Ko The curves of (l-p), ap/6T, and (l/p)oap/6T, and the average value of f for the flat temperature are useful for illustrating the effect of a nonuniform temperature distribution and the determination of an effective temperature. Let the effective temperature, Teff be defined as that temperature at which the case of a flat temperature gives the same value of the parameter of interest as the case of a parabolic temperature distribution in the fuel. A characteristic effective temperature applies to only one parameter under a given set of conditions. Since the functional dependence on temperature of the capture probability and its derivative are dif

88 ferent, it is not expected that their effective temperature would be the same. By superimposing the values of capture probability for the parabolic temperature distribution cases on the lower curve in Fig. 1, the effective temperatures are obtained. Similarly the effective temperatures for ap/6T and (i/p) (6p/3T) are obtained in Fig. 2. Values for Teff) Tav and Trm are given below: Run Ts,~K Tp, OK Tav, K Trm, K TeffOK ( 1520 for (1-p) 4 840 2800 1820 1670 140o for )p/6T 1400 for (l/p) ap/)T ( 2730 for (l-p) 5 1820 3780 2800 2720 2800 for ap/6T 2800 for (l/p) ap/6T For the higher temperature run, the parabolic temperature distribution can be characterized by the average temperature since Tav, Trm, and Teff appear to give the same values of (l-p), 6p/6T, and (l/p) ap/)T within limits of the statistical error. However for the lower temperature run, which corresponds more nearly to practical operating fuel temperatures in a power reactor, the average temperature is a poorer approximation to the effective temperature. Use of the average temperature to characterize the parabolic distribution in this case would cause the following errors:

89 % Error caused by replacing Teff by Ts,~K Tp, K Parameter Tav Trm 1 - p 1.7 1. 840 2800 ap/aT -10. -7. (/p) ap/6T - 9. -6.5 For the case of (l-p) the discrepancy in using Tav is outside the limits of probable error, however for the other two parameters there is a slight overlap in the error limits. The concept of a constant 5 (Eq. 6-1) for a large temperature range also provides determination of the effective temperature. For the first run of parabolic temperature, the effective temperature is found from Eq. (6-2) to be 1380~K, which is within the statistical error of that given for the Doppler coefficient of reactivity. The reasons that the average temperature appears to be a much better approximation to the effective temperature for the second run at the high temperatures than the first may be attributed to two factors. 1) The peak to surface temperature ratio is smaller in the second run. 2) At higher temperatures the resonances become lowered and broadened by the Doppler effect and therefore the self-shielding factor is decreased. At the lower temperature, neutrons incident on the fuel surface near the resonance peak do not travel far before being absorbed, and hence experience cross sections corresponding more nearly to the surface temperature. However in the limit of extremely high temperature the resonances are com

90 pletely broadened and the self-shielding factor vanishes, hence incident neutrons could penetrate further into the fuel and thus experience crosssections over a larger range of the temperature distribution. At this limit the volume average of the cross section would equal the effective cross section. For a parabolic temperature distribution this averaged cross section corresponds to the cross section at the average temperature. The error in using the average temperature concept is due to an attempt of replacing the average of a function by the function of an average, however as the rod temperature becomes higher and the temperature distribution becomes less extreme, the results indicate that this error decreases.

CHAPTER VII CONCLUDING REMARKS A. CONCLUSIONS Pearce10 has indicated that there is a definite need for investigating the effect of a non-uniform temperature distribution in the fuel on capture probability and Doppler coefficient. He states that "there are two distributions of temperature within fuel elements that are of particular interest: (1) the distribution present in an operating reactor, where the temperature is highest away from the coolant channels, and (2) the distribution present in a reactor runaway, where the temperature is highest in the region of highest thermal flux. Both these distributions are difficult to realize in either activation or reactivity experiments and there have been no experiments reported on the Doppler effect in these distributions. The theory of Doppler effect in non-uniform distributions is also in an unsatisfactory state." The current theory reviewed in his article "shows that one cannot simply take a surface Doppler coefficient with a surface temperature and a volume coefficient with a volume temperature." It is hoped that the program written in this work partially fulfills this need, since REPAD provides the first means of evaluating the effects of an arbitrary, radially dependent temperature distribution without introducing approximations that are difficult to justify. The numerical 91

92 results are for a parabolic temperature distribution since it closely approximates the distribution existing in a U02 fueled reactor.21 The results of this work indicate that the substitution of a parabolic temperature distribution in the fuel by a uniform distribution at the average temperature as suggested by Keane22 and Dresner23 can be done for the conditions considered in Section VTI-B only if one is satisfied with an error of about 2% in resonance capture probability. However the results further indicate a trend toward a larger error as the peak to surface temperature ratio becomes larger than 353 and/or as the average fuel temperature becomes lower than 1820'K. Both of these tendencies are quite probable in operating power reactors using U02 fuel. The error caused by using a uniform average temperature as a substitute for the parabolic temperature distribution is much greater for the Doppler coefficient of reactivity, being of the order of 10% for a parabolic temperature distribution in a.25 in. U02 fuel rod with surface and peak temperatures of 840'K and 2800'K respectively. The results indicate that the Doppler coefficient has the same trend toward a larger error as indicated above for capture probability. This 10% error is significant when one considers that the roppler coefficient is of prime importance in reactor safety because it reacts promptly with changes in fuel temperature. Furthermore, as pointed out by PearcelO the temperature coefficient of reactivity is more sensitive to the Doppler coefficient than to changes in thermal utilization and resonance capture due to thermal expansion of the fuel in power reactors.

95 The root mean square temperature (Section VI-B) represents a better approximation than the average value as a substitute to a parabolic temperature distribution for purposes of evaluating the capture probability and Doppler coefficient, but it is still not satisfactory for the usual operating conditions of a power reactor. The Monte Carlo program REPAD includes a new and direct method of evaluating the derivative with respect to temperature of the capture probability that was rigorously derived from the Neumann series solution of the Boltzmann equation by utilizing the mathematical model of Albert27 for random walk chains. This formulation for determining the Doppler coefficient circumvents the approximations inherent in the other two methods discussed in Chapter I, and the results indicate that this new method compares favorably with these other methods. For the special case of a flat fuel temperature, a direct comparison was made with the method of describing the temperature dependence of the resonance integral by an empirical T1/2 law. Both methods agree within experimental error. The use of a finite difference approximation to the derivative would give a large error for our results since the temperatures were widely spaced, and therefore this approximation was not used. However Mortonl4 used this technique in his Monte Carlo program for a uniform temperature distribution and obtained a probable error of about 25% for 4096 histories. For 12,000 histories one would then expect the finite difference method to give a probable error of about 14%. This compares with an average of about 7% for the runs reported in the results using the direct method, and about 12% using the empirical

94 method. Furthermore both the empirical fitting and finite difference techniques require at least two runs at different temperatures in order to obtain an estimate for the Doppler coefficient. Therefore for a given probable error, the new method allows a savings in computer time by a factor between 6 and 8 in comparison to the two other methods. In summary, it is felt that this work fulfills two purposes. 1) It illustrates the application of the method of constructing random walk chains from the Boltzmann equation to develop a mathematical formulation for determining the Doppler coefficient directly in a Monte Carlo process. The method developed produces a considerable savings in computer time. 2) On the basis of this work the quantative influence of a fuel temperature distribution on the resonance capture probability and Doppler coefficient of reactivity is evaluated by means of the present digital computer program. B. FUTURE WORK There are many interesting problems that can be investigated with the Monte Carlo program REPAD as presently written; these include the following: 1. One could evaluate the temperature dependence of the resonance capture probability and Doppler coefficient for the 54 resolved resonances and the unresolved resonances of U-238 and Th-232,8 in order to compare with the experimental results of Hellstrand, Blumberg, and Hornerl7 for a uniform temperature.

95 2. One could investigate further the apparent divergence of the effective temperature for a parabolic distribution from the average temperature as the average fuel temperature decreases and the peak to surface temperature ratio increases. Also one could investigate the effect of fuel rod diameter and fuel-moderator volume ratio on this trend. 3. One could use temperature distributions more appropriate to a reactor in a runaway condition where the temperature is highest in the fuel region of highest flux. With some modification of REPAD, the following features could be included: 1. Doppler broadened U-235 resonances in fuel of low enrichment using the Breit-Wigner single-level formula and Doppler broadening functions; 2. variation of fuel composition in the separate radial fuel zones; (This is important for fuel burnup studies.) 3. provision for the determination of spatial dependence of Doppler effect;55 4. provision for p wave and higher angular momentum neutron interaction for the unresolved resonances; 5. provision for a more generalized geometry to include a rectangular lattice cell; 6. use of the Breit-Wigner multi-level resonance formula for U-235 at high enrichments; 7. provision for a spatial and energy dependent source to investigate source effects on Doppler coefficient and capture probability.

APPENDIX A DERIVATION OF GREEN'S FUTNCTION The solution of Eq. (2-2) can be determined by Fourier Transform techniques, however in order to gain insight to the method, the special case of spatially uniform cross section will be considered first. Case I: -<- L (E); function of E only. Equation (2-2) then becomes v*-C G C((E)G~rE| r'ltE' 2) - k SF db J S E ~) J2 ~ ~s )~ (A-1) The definition of the symbols and the boundary condition remains the same as in Chapter II, Let the Fourier Transform of the Green's function be CG T r, r E (A-2) then by definition of the inverse transform 2 -)eX (7)3(A-3) also let 96r t k tv (A-4) 96

97 In this case the transform of Eq. (A-i) with respect to the variable r becomes _Q _k G + S(tE)G = e -- Y(E-'E ')g ( )~ (A-5) hence e(k r -' E' (.r n'J (A-6) define: a~'- g ( E-E'3,(-~l') In integral form, the solution is then obtained by the property of the inverse transform and is f_; k. tr-r') G (Jr} EaR I rJE'JQ) 72';)3 'ty,,) -. (A-7) k I- Z;/E) — iR~rG (A-7) define _ - _ r -r' - (A-8) -_~ = ~ (A-9) Let 2 be along x axis for convenience, then

98 8z -,' ex +J + y +h RY where i, j, h are unit vectors along x, y, z axis respectively, then (A-7) becomes +4 Jie -t+ Lkk (A-10) k.x k,k but one definition of the Dirac delta function is /Heoc &(>) / e -= (A-ll) Hence 2G = z J (A-12) kX but according to Volume 1 of the Bateman Manuscript Project,56 -17 ikXG, x 1;,if g,, o; (_.(- r-") ) ie ~~(A-1 let (A-14) now

99 but note ( R 2 12 j- (e, 7 6e~(A-Z1) therefore Since(E)R~$ (A-16) Since d3 =| R|Z d/R4, then ~G~,E,/rr E' '_Q Jdr= G- J3R (A-17) 5.(E) -a 3 'e Z (,_e)S'(nn. ' Or(-E') - is a solution to (A-l) and is the Green's function for the case of spatially constant cross section. Substituting (A-17) in (2-9) and integrating over all E', p2', the integral equivalent of Eq. (2-1) can be written in this case as d3r - S'() [L-(r-"- -)I (A-18) t ", I E, 2 E'r J _el, j

100 Since J 3r ' _ 13? = ~R 1Q IL-r I |oP-r | (A-l9) and r* r R = _s then an alternate form is aa, 6,(E)R iRi ~-lr e, E~ D (A-20) is the integral equivalent of the Boltzmann transport equation with cross section spatially constant. Case IIT 9 4 (E r); function of both E and r. With insight afforded by Case I, Eqo (2-2) can be solved with the cross section a function of E and r, Physically the solution represents the transport of particles with energy E in direction Q from r' to r only if Q = QR (A-9). We can impose this condition immediately on Eq. (2-2) by utilizing the delta function; 52(QOR). Let s be in direction Q then ' V s (A-21)

101 Equation (2-2) then becomes d- G + (E,C)G (_,,,nl',F, ) ds (A-22) Define: a - 4.JSIC2( g) s ) S(6 -E') Multiply Eq. (A-22) by and note e,_ -a ~( r) Since I d for non-trivial solution then (A-22) becomes (A-23) Let E -~ e (_A-24)

102 then (A-23) becomes Ids (A-25) Taking the Fourier Transform with respect to R of (A-25) v+ J -t dR ~ -r_ =.o(Re) de a (A-26) hence ~E -,k and by definition of the inverse transform GE =2r- -e k = y/(R) z (A-27) -00 where again (CR) i~.P\ t~his effect included 7 - (R) = t fU > D 0 in 82(_._R) ) Therefore (Cc E,| E, L( E_'. n,) ) a (E) r (Q.Q) 2, BQ Q) /( - AE2) Hence (r, J_,Ez n') 3r' - Cd3R te (A-29) is the Green's function for Eq. (2-2).

APPENDIX B MATHEMATICAL DESCRIPTION OF BANACH SPACE This section presents the mathematical description and some useful properties of the Banach abstract space, and is summarized from the book by Riesz and Nagy.30 BACKGROUND The fundamental spaces LP (where 1 < p ao) and the space of continuous functions C are particular cases of the Banach abstract space (generalization of Hilbert space). A Banach space is a set B of abstract elements f, g,...which satisfies the following conditions: (a) B is a linear space, i.e., the operations of addition and multiplication by real or complex numbers are defined for its elements, and these operations follow the rules of vector algebra. (b) To every element f of B there is a norm with properties: norm I i/ II| - if and only if f = O, a.e. Il F +9 C + C1 11 1 1 - { a = constant (c) B is a complete space in the following sense, if [fn) is a sequence of elements in B that satisfies 103

104 I/1K7I(t o + > ~ for k-n n oo then there exists an element f of B 3 fI( - +- I( -* C for fl - c(Riesz-Fischer Theorem) Remark: the elements of the set may be points, measurable functions, transformed functions, or functionals (or operators). Measurable functions can be simply defined as a function which is the limit a.e. (almost everywhere) of a sequence of step functions (see Reference 30, p. 44). The norm of the measurable function defined in a measurable region e, for which f IP is summable is defined in Banach space as norm j ( If~ i)JIt; in LP space (B-l) a n C or L space (B-2) A linear operation or linear functional belonging to Banach space is an operation which assigns to every element f of B a numerical value Af and which is 1: additive: Af,+ ) = Af, 4 A -F 2: homogeneous: A (cf) = cA f, where c is arbitrary constant 3: bounded: there exists a constant M such that for all f....

105 IA |4 ~ 1Iff/ (B-3) Denote the smallest of the possible bounds M by MA or IlAll called norm of the linear functional Ao Hence IA I i ^ jAi(f /I (B-4) hence IJA 11 rVA (B-5) A linear transformation belonging to Banach space is a tranformation of the space B into itself which transforms the element h into the element Th and which is 1: additive: 7(, g ) - T 1, + T 1 2: homogeneouso T(c 4 cT I 3o bounded: there exists a constant M such that for all h /ITh11 M! " Ii if MA is the smallest of (B-6) these bounds then norm = -J T4I - MA (B-7) Multiplication of linear transformations is defined as (T, T,) h T' (T, h) (B-8)

106 and = TT = T(T) (B-9) (B-9) for iterated transforms, Their norms have the properties: I/~J ' I-~~[ II ~, + T. Ii ~1 s T T, IIT,1 1 1c#7i) = ( cj| || 7 || T 11 iJ+Tz || — < ||~~ T ||(B-10) I| T, Tj I <~ J! Ti | II Tz (B-ll) hence for iterated transformations: iiT2|i II2 | T; |i | ' T|| (B-12) Continuity of linear transforms: Every linear transformation is continuous in the sense that if the sequence fhn] converges strongly to h in B, then the sequence jThnj also converges strongly to Th. Proof: (lh-ITbwTHlII= -IIT'll( -))j IrI iVn-ki-l0 (B-13)

APPENDIX C PROPERTIES OF DENSITY FUNCTIONS AND OTHER PROBABILITY FUNCTIONS Consider a probability density function, or in shortened nomenclature, density function or frequency function, f(x),which satisfies the condition: _ U(x)dx = 1 (C-1) then for the random variable, PX~ '5: x j = x)dx (C-2) The cumulative distribution function (or distribution function), F(x) for the random variable, is x.:,, — x 00 -00 hence) Pit a -POI F (hio 10 if 4(x) =O far X S 0) then p i g G?= F(o) = o The expected value of a function g(x) w/r to the density function 107

108 f(x) is -oo —.00 This one-dimensional concept can be generalized to many dimensions. Let l, ~ 2... k be k random variables on sample space, a2with density function Pf_,~X; _ -x. sk; -=: ~,i,,'2,.. xa ) lx, d,:. - xt (C-5) The multivariate distribution function is then f',; /'; "Xk ) ),.-Y), ) /x,. JXk (C-6) - ---- o-i The normalization condition is then If ko If k = 2, with density function, f2(x,y), then - ds( Jqf<)yJ Y tJ l' =J (,X, C/X ~(c-_8) are called margina density funtions for x and y respectively This are called marginal density functions for x and y respectively. This concept can be generalized by considering the above frequency functions and normalization conditions: 00 00 $-C/)=~~ x,,- Jx,) (x ~ ex* (c-9) - oO., x,,, ).. (r,..,,j.,, Jr...(-00

109 where n k The corresponding distribution function is then Xt X, ~C oo ( s (A - j Xh Idx, 4, txr, (C-10)....- -oe Conditional density functions are defined as fk (xt;.xXk) k _r (x,...x ) probability of random variables xn+l,..xk lying in intervals dxn+l,...dxk, respectively, upon given xl, *..Xn. Conditional probability is defined as P fEZ EIl}=e, (C-12) P{EIZ

110 where E1 and E2 are events in the sample space 2, and P E1, E2] is the probability for both events to occur. For mutual independence we have for probabilities P[E2jE1J = PjE2J since PfE1l E2} = 0 and P{E1lE2} = PE1)} P(E2} and for random variables, i PfX,,t, s t%+bsc; L=1,~f 3 -n 1tIX;C~rgt~lx,4 (C-13) =) i (x;)dx4 (C-14) L=/

APPENDIX D PROOF OF FINITENESS OF NON-ANALOG RANDOM WALK CHAINS DESCRIBED IN CHAPTER IV In a heterogeneous reactor containing an infinite array of lattice cells, the neutrons slow down by striking moderator or fuel atoms. Under the assumption of elastic collisions and struck nuclei at rest, the longest chain from birth to death results if the neutron were to slow down to the cutoff energy entirely within the fuel regions. The expected number of collisions with one kind of nucleus sustained by a neutron in slowing down from source energy Eo to cutoff energy Ec is approximated by -B g h j (e~E (D-1) where S is the average logarithmic energy decrement for the given kind of nucleuso Therefore the expected length of the chain for the non-analog ramdom walk process described in Section IV-A, in which all collisions are scattering collisions, is n In ( C_ (D-2) In measure theory language using Lebesque integration 111 - ph{; to b)(D-3) 111

112 where p(c) is the probability measure of the set o Q is the complete sample space 4(Q) = L1 Consider the set A and its complement defined by (D-4) A' - Ix "(X) = b 1 }= oo Thus cr<>/O nc/L)> 4r ) { t hl/W) fh/ y) A1tX) (D-5) /4 A' But by the definition of the set At ~h h -frdw)&dw(X) + bl4 (D-6) A Since <n> is finite and b is unbounded then 4(A') must be zero in view of Eq. (D-2). But 4(AT) is just the probability measure of the set containing random walk chains of infinite length. therefore the probability of a chain of infinite length is zero and therefore Eq. (4-7) becomes piL=m? = (D-7) Osborn has proposed using the hydrogen model to prove (D-7) for the case of non-analog random walk chains in which all collisions are forced to be sc attering collisions

113 In lethargy units the scattering frequency for elastic neutron collisions with the struck hydrogen nucleus at rest is given by / II F(IA'a'a) jgJ I _ — u - /, (D-8) 8= 0) otherwise Let Pn be the probability of crossing lethargy u in exactly n collisons, then Pi L n =P With the source at lethargy zero, then 00oo 00 p I e J~ = e (D-9) ~LA.t~~ P3 - J F(~, ) F(u, )" d>t': - {,~/ Le (D-l1) (/_ ) -t8 (D-12) ag /A /& U o<~ The nth term is determined by induction to be P n- { -(D-13) /.-~

114 The limit: /JM PP o 0 ~(D-14) can be shown by taking the logarithm of (D-13) and using Sterling's formula for the factorial with large n. And therefore p r3=P 0 Q(D-15) It is also noted that 2 Rn = 1 (D-16) by substituting (D-14) into the summation and recognizing the resulting sum to be an expansion formula for eU.

APPENDIX E SCATTERING PROBABILITY FOR STRUCK NUCLEUS IN MOTION In the previous sections, the Monte-Carlo process was developed under the assumption that the collision kernel ~C(E QI E/, s; 1) could be defined for elastic collisions of neutrons with the struck nuclei at rest. The motion of the nuclei, which produces the Doppler broadening of the resonances, was considered in determining the total cross section in the transport kernel and the ratio of scattering to total cross section. However once it was determined that a scattering collision occurred, the new energy and direction of the neutron were selected as if the struck nucleus were at rest, i.e., it was assumed that at the energies of interest (above 10 kT) the motion of the nucleus had negligible effect on the collision dynamics. For resonance scattering, this last assumption is somewhat questionable, and it is the purpose of this section to investigate this assumption in some detail. The method employed for this discussion is to derive the correct scattering frequency for the case of elastic collisions and the struck nuclei in motion. The correct scattering frequency and the assumed scattering frequency can then be compared and perhaps some evaluation made of the assumption. The work is an extension of the method of Osborn57 in which the scattering frequency is 115

given in integral form for an arbitrary differential scattering cross section and in closed form for a constant scattering cross section in the CMCS. Since details are given in Reference 57, the steps are only outlined here. The desired probability is P(vv'I; A,T) = the probability that a neutron of initial velocity v', shallupon suffering a scattering collision with a nucleus of mass number A and whose velocity distribution is Maxwellian and characterized by a temperature T, have a final velocity in d3v about v. This scattering probability can be derived by formulating the scattering rate density of neutrons from all possible initial velocities into d3v about, v. After a coordinate transformation this can be written as V/ where n(rv') is the initial neutron density in phase space. P(lI( T A).____'J __, (E-2) Primes indicate precoilision variables vI = initial neutron velocity V' = initial atom velocity m(V') = Maxwellian distribution of atom velocities.

117 = vV -, _ c Upon evaluating the Jacobian, assuming that the differential scattering cross section is isotropic in the CMCS, and performing the azimuthal integration, the scattering probability can be written as? (Vff)Pi) 2ZA\I k ()Js (E-3) (__)4_(v_ _ v ) 44J~' where N = number scattering atoms per cm3 JrO = zero order Bessel function of first kind D 2 kVV COS B - - D _- vt' It has been assumed that Gc is the center-of-mass scattering angle, then the angle 7, and the relative speed R', are related by ~ - - e ____c =1

118 The usual form for the scattering cross section is related to the form appearing in Eq. (E-3) under the assumption of isotropy in the center-of-mass coordinate system. ls- ( (2c,) 4 For the case of potential scattering, for which the microscopic crosssection is constant in the center-of-mass system, the scattering probability can be written p~y(I'j TA) =,),_ where __4- V- / $VIZd (' t ) Z (E-6) V v' - _ ')/ In the limit as hv' becomes large (v', and upon integrating over angles, the scattering probability for potential scattering reduces to the common probability for struck nuclei at; rest, JFv-;V/ n. -- pv | A,T)-v -- (A-i+-)2 Zvclv v AV1 ycY 4A 4 v,; A- - FW(A-~E E )' E ' -Iz

119 For the case of resonance scattering, the Breit-Wigner single-level formula is introduced into Eq. (E-3) before integration. Straight forward integration does not seem possible in this case but by expanding the cross section in Tchebicheff polynomials of the second kind, integrating the individual terms of the series as Hankel transforms and thereby producing Laguerre polynomials, and then performing a sum of some of the terms, the resonance scattering probability can be written as F, (f v-U7A) _RM / J JhX / ^)*Z i e 'A Notes I + US 4- (7r 7r LVI 5) e (E-8) where C 2 -k C, /. z. Zh ( / {a gE<'-:, Oo;'So-C,._h=o ___ c, _14 X (EC,_ - [I =L f G,< - e8) -) C

120 Un is the Tchebicheff polynomial of the second kind Ln is the Laguerre polynomial. In the resonance energy range of U-238, Eq. (E-8) can be simplified by replacing (1+BS) by 1. This causes a negligible effect for all scattering angles except backward scattering, in which case the maximum error is about.8% at room temperature and 8% at 3000~C. The term in the exponential equals zero when the scattering angle and the final speed are determined for the case of struck nucleus at rest and becomes large rapidly for other scattering angles and final speeds. Hence the motion of the struck nuclei cause the probability function for the final speeds to be distributed in Gaussianl form about the speed corresponding to struck nuclei at rest. and therefore, the exponential term causes Eqo (E-8) to approximate the Dirac delta function,and qualitatively the effect of the motion of the struck nuclei is small, However it is difficult to evaluate quantitative effects due to the complexity of the resonance scattering probability. Fuirther insight can be obtained by considering special cases of scattering. Since for heavy nuclei, scattering is almost isotropic in the laboratory system for the energy range of interest, the average scattering angle is j/2(;os 0 = 0). Change the scattering probability from velocity to energy and direction coordinates. PB (E, / i T, A)JE J (E-9)

121 _ COS 9 9 = scattering angle in LCS. Let E* be the final energy of the neutron for a given scattering angle for the case of struck nuclei at rest, i.e., for a given scattering angle in CMCS, E 'L' C K = 'r + 4 +r C0r- _ (E-10) where l -/ A-/ (IA ) or for a given scattering angle in LCS, -E 1;' z z (E-ll) --- = F + a OP + Z 79 Let the final energy of the neutron due to motion of the nuclei be written as a=E (I+) (E-12) where A~ - E E -~E X E - E_ The approximation that A << 1 can be satisfactorily made since the mass number of the struck nucleus is large. The scattering probability for large scattering mass and it = 0(0 = i/2) becomes Pa (I- a(zkT: A), t,lo< ~ z e /V a /\' < a ' ( (E-l3); (2a7;T/ '- v /.~', z

122 where The coefficient of A2 in the exponential is about 7000 at room temperature for E' = 6ev, hence A must be very small, and the final energy very close to E*, But since the ratio of the neutron half-width to the "peak" energy of the resonance can also be quite small, A can be comparable to r/E. Therefore near a resonance let A be proportional to some multiple, y, of the half-width, AE r E_? Also let r be the reduced total width, i.e., 0 _ iTherefore 8AT AT v Equation (E-13) in approximated form becomes P 0 (%.o E7,- A) ~ _ C(A p... A (E-14)

123 Intuitively, one would expect that if the motion of nuclei in the collision dynamics is an important consideration for the resonance escape probability and Doppler coefficient, then the motion must manifest importance for the strong resonances which contribute a large portion to the resonance absorption probability. However it would seem that a definite trend in modification of the collision dynamics would have to be established in order to realize an important effect, e.g., if the motion of the nuclei could with reasonable probability cause the neutron to lose enough energy in a single collision to completely jump over the Doppler broadened resonance, then surely the effect would be important. For the strong resonances the practical width (resonance width at which the resonance cross section equals the non-resonance cross section) is about 50 half-widths (x = 50) and the reduced total width is about 10 mv/ 4Tv. For the strong resonances the exponential term is then approximately, eA Z /yZ b T 02 = 3 e(v. roo temp) e - / AT -,zs3 e v. (/3$oooK) For a unique case in which the expected energy loss for the struck nucleus at rest is 25 half-widths and if E' corresponds to x = 50, then E* - Eo. Furthermore if y = -50 then the neutron the neutron could lose enough energy to jump over this hypothetical resonance, However the exponential term would be

124 A?02 ( kT =.o273s ev e-3 k 7 -f 3 iV, and therefore the scattering probability, Eq. (E-14), would be so small that the event would be extremely rare. In reality it takes on the order of 10 collisions with U-238 nuclei for the neutron to skip over the strong Doppler broadened resonances of U-238. Therefore the motion of the nuclei probably cannot cause an effect as important as jumping over the resonances. By studying the combination of exponential and denominator terms, it seems that on the average the nuclei motion can cause the neutron to lose more energy when E' > Eo and less energy when E' < Eo, than it would lose under the assumption of struck nuclei at rest. Intuitively these small changes would not necessarily change the number of collisions required to skip over the resonance. And since the ratios and remain constant throughout a given resonance, under the assumption of negligible potential and interference scattering, then the probability that a neutron suffers a sufficient number of scattering collisions to skip across the resonance does not change.

APPENDIX F EVALUATION OF THE DOPPLER BROADENING FUNCTIONS AND THEIR DERIVATIVES Various methods of evaluating the Doppler broadening functions, which were introduced in Chapter IV, will be discussed in this appendix. Among these methods, a new one is presented which utilizes the Moebius Inversion Technique for evaluating the inverse of a Fourier transform. In Chapter IV, the Doppler broadening functions were defined as _ (F-) X (x, a) e (F-2) which satisfy the following partial differential equations and boundary conditions: 9(4,o) /..~ 2 xz _ i/ a'n (0d 4 = X e e r - e and ~aX ~ ~ X~~ (x,o) = X w_(oe) - o 125

126 where -X erfc0j) = -& ( I= /. er#y) It can also be shown that thie Doppler broadening functions also satisfty the following differential equations: = e.........:i/ - 7 - I -- x4) 40 464 Equations (F-3) and(F-4) offer a convenient way of evaluating the derivatives with respect to temperature of the Doppler broadening functions. Recall that values of these derivatives are needed to determine the Doppler broadening function in the Morz:te Carlo program. Equations (F-3) and (F-4) are useful fow-r finding higher order derivatives-of Psi and Chi. By recognizing the exponential and non-exponential parts of the integrands of the Doppler broadening functions,(G-l) and (G-2),as inverse Fourier transforms and recognizing a resultant form of the Dirac delta function the Psi and Chi can be written as ) (F-5)0 t~x9) f cosw Y (aZ c~~~~~~~et

127 Xl~a = i -- v e in dv (F-6) which are the inverse transforms of, 00 Q(2 e) -- =. Y '(xie)J cosxv d (F-7) a 0O 12 (ve) -z Xi =e(x,3) x/ xV rx (F-8) Therefore Psi and Chi can be thought of as the inverse Fourier cosine and sine transform of Omega. 1. MOEBIUS INVERSION METHOD Goldberg and Varga58 and, independently, Duffin59 have developed a method of numerically evaluating the inverse Fourier transform by employing the Moebius inversion formula from number theory and the Poisson formula from Fourier analysis. In the notation of Goldberg and Varga, li't) = F/A) C tL A/ (F-9) is the inversion of F )-4 (F-10) o

128 Upon being given an arbitrary sequence of numbers ak let the sequence (b/ be defined by 6 = /J (FTll) k / / O _ 3S 4 - - (all divisors of m) Then the Moebius inversion formulae are given by Gl(t)= - k (F-12) m(t) = jj_ 4, G (nt) (F-15) n=1 Suppose the ak are given by k, (F-14) [ = k ~=k- - / 2,3 2,,, then bn are defined by ( 27-~ / a~ - 6, -- '- o i B ) ce~(F 15) I -/% ~=/1zZ3 -,, Where the I-t are Moebius nu-mbers, 1 m = 1 = (-l) if m is the product of s distant primes 01 O, if m is divisible by a square~

129 If we let G. tr~+ (-i F( 7 (F-16) k and use the sequences (F-14) and (F-15) and Eq. (F-10), then Oo k -(F-l7) and hence oo. -7 = 4/, is a solution of Eq. (F-9). The fact that the result (F-18) is a double infinite series is somewhat disconcerting, however in applying the method of Moebius inversion to the functions Psi and Chi certain techniques can be employed to obtain more rapid convergence. This will be illustrated for the Psi function. By defining.- 1V- e f'ev) =J - ~ (F-19) then Eq. (F-5) can be written in the form 00 (xW _,e) ( F (e, v) co. x v ('-2o) 0T

130 hence 's (E@) ='Yi~et ) cos xv dv (F-21) Let m = 2n-1 L _ r ' DX = | } S )S -' e?_, x-/~ -- @ x )(F-22) Z.n X Defining = _,', -(x,6) becomes v. C iz ' I,~ J IL,l-Ok F~ m x (F-23) zX G z, ~ / 2 +' Since it is desired to evaluate Psi for a set of values of x and ~, it is convenient to evaluate the exponentials on the computer only once and store the values —this of course requires that the two series be truncated, The non-linear transformation method of Shanks60 was employed to speed up the convergence of the k-series. In order to obtain rapid convergence of the n-series the terms were compared to the terms of a series that can be written in closed form. Insight into this method was obtained using the Euler-Maclaurin

131 sum formula on the k-series, expanding the resulting functions in Taylor series, and recognizing the n-sum of the first few terms of the k-sum. Employing the Euler-Maclaurin sum formula, Eq. (F-22) can be written after some manipulation as + I-E) +rn w (u( ) - r H,(- e 2-)) (F-24);Zt X B()- 3 \2) 3! 3 where Bi = Bernouli numbers Hi = Hermite polynomials. After expansion for mo and mx greater than 1, ($(Zk) _ / {-(l ); _; ~+!4 ('1(F-25) + 7 ( 7 ( Now consider the sum Avf e expansinfr 4~ an x m )m 4ra ohaX n,-/ 00 6/ 2,7-; G; (>; - xJ ~~~(F-26) Nrtah 62x~~~~~~~) "~~o

132 where we define 6Y tX (d g ) ' 02 ) (F-27) The sum is performed by using the Moebius inversion formula and the Fourier cosine series for and letting Another useful sum can also be obtained, namely 04.7 - 12mX-/ r n= where i -/ Z- + - -- 2+ +, ~ (F-29) In this case the k-sum involving x alone can be written in the closed form: G2. Therefore, for large x and I and by comparing Eqso (F-27) and (F-29) with (F-25), 6 /m/} C78+z,

133 hence This same result can be obtained in the limit of large x and 0 by considering (F-1) as a Weierstrass transform and formally evaluating by the method of Hirschman and Widder.61 For small x and, but large m, we still have C- ~x) L G, S~G Hence the sum (F-23) can be written y -x'1) -p = w + tt-, I n 1 An G (X) - (, (b/,)) Zml)] L + 6 4 / fX 2 ~ (ZK)2 (F-30) '"- " L 42r '~......................- )- '+ oo — '-/ 1K 4 I Zmx 4J - Z -' -n7 x - 2m -6

134 In the program written for the computer and using this method, it was desired to evaluate Psi for x and i quite small, therefore the subtracting of G1 and G2 was delayed until the 6th term of the n-series so that mx and mi were less than one, The \K 1 and Y2 terms were correspondingly modified. For this case the series for Psi becomes _ /. o 00o /s7 S7?o 88 V(,/} -Z 0) - r /41 I) + + - ' (F-31) where the indicator In is defined as B s t,v e ) By this method values of Psi were calculated to at least 5-figure accuracy and were compared to those obtained by methods described in the next section. 2. OTHER METHODS Other methods of evaluating the Doppler broadening functions have been tried, The Runge-Kutta technique of numerically evaluating dif

ferential equations on a digital computer was applied to (F-3) and (F-4). The numerical evaluation of the partial differential equations satisfied by Psi and Chi is used to build up a table in the REP Code.l6 Furthermore the derivation of asymptotic and convergent series for Psi and Chi has been done by others.20'42 None of these methods (including Moebius Inversion) of evaluating Psi and Chi over their full range of arguments are entirely satisfactory for use in the program REPAD, either because of inaccuracies for certain values of x and theta, or because they require too much time to evaluate on the computer for this Monte Carlo program, or both. In consideration of these various possibilities of evaluating the Doppler broadening functions and the necessity of evaluating their derivatives, it was decided to use a combination of convergent and asymptotic series. Define: Convergent series: Y(x o = e 5~! e) (F-32) 2 o

136 where:rO = 7Y(oweJ - } e er1-c( (40).+, rn) - =1 + Z 4- - (4-)n. Asymptotic series: '(xroe) = Y.~: (x,e>= rn o () ( X (X, g) _ e W (Z~.i)! (-,) n (F-35) g~ (z.c () 2 l (z.,7) + )n/,t.- x~z Z (F-57) Itt — Otl=

137 where The convergent series were used in REPAD to evaluate Psi and Chi for values of x and G such that 0 < 3. The derivatives of Psi and Chi with respect to theta was determined by (F-3) and (F-4), when the convergent series was used. In order to obtain a 1% accuracy in the derivatives it was found that the convergent series should be truncated only when the contribution of the last term in Chi to the sum was less than.01%. The asymptotic series were used for values of x and G such that 0 > 4. The asymptotic series for the derivatives, (F-36) and (F-37), were obtained by differentiating (F-34) and (F-35). The resulting series were checked for accuracy before using in REPAD. For x and G such that 3 <, ~ 4, it was found that by adding the first term of each convergent series to the asymptotic series for Psi and Chi, satisfactory accuracy resulted. Likewise the derivatives with respect to G of these first terms were added to (F-36) and (F-37) for this same range of B. The number of terms to be used in the asymptotic series was determined by a careful comparison of evaluating techniques; the decision for specifying the proper number of terms is built into REPAD.

APPENDIX G COMPUTATIONAL DETAILS USED IN THE MONTE CARLO PROGRAM, REPAD, FOR SELECTING RANDOM VARIABLES 1o RANDOM NUMBER GENERATORS The two pseudo-random number sequences used in the Monte Carlo program, REPAD, are obtained from generators using the Method of Congruences, On binary computers with 35 numeric bits and one sign bit, this method involves multiplying two numbers together and taking the 35 least significant bits of the product as the random number. The pseudo-random numbers for the first generator are given by 2/ 3 an - 3 a rood 2 3s-,3-% (G-l) for n 1,22,o. The initial integer, ao, can be read in as an arbitrary number, however provision is made in REPAD to use ao = 235-1 if desired. The pseudo-random numbers for the second generator are given by 14. b. _ o 5 L. azCb_ o (G-2) for i l, 12,.. (G-l) is used to generate the initial random number of each history and (G-2) is used for successive elements in the chain; therefore, bo = an for the nth history. In practice, the random numbers, bi, are divided by 235 to obtain fractions on the interval (0,1). 138

139 2. GENERAL METHOD Given the probability density function, f(x), Pjy <4y/d = rxy) (G-3) its cumulative distribution function, F(x), defined as Pf So! =F0) - fWdJ (G-4) The random variable x is selected from its probability function by first selecting the random number, RN, RNC (0,1) then determine x such that x A) = F/X) = f f{i} 1y (G-5) If the inverse of F(x) is known, then x = F-1 RN. (G-6) 3. NEUTRON SOURCE POSITION AND DIRECTION For a source uniformly distributed in the circular, semi-infinite, lattice cell with circular fuel rod, only the distance, r, from the axis is required because of symmetry. Let R = radius of cell

140 then the density function is Pr r Ztda g 7, for 0 r R = R2 o, otherwise and r PIjr} = F(r) = | _ _ Now select RN C (0,1), find r such that RN = F(r) = r2/R2 i.e., r R= R vo (G-7a) Alternatively,48 select RN1 and RN2 C (0,1); let RN = max (RN1,RN2,); find r such that RN = r/R r = R~RN = Romax(RN1,RN2) ~ (G-7b) The neutron direction is selected from an isotropic distribution. The density function for the polar angle, y, measured with respect to the radial line extending from the source point r to the axis, is Y A y'd: = B 8-, for 0wi y< e o1, otherwise

141 The distribution function is then Now select RN2 (0,1), find cos y such that i.e., (G-8) Then, Y, - V / C-cos (G-9) The density function for the azimuthal angle, J, measured with respect to the same reference line as above, is:yi v=' y, -for 0'.2) r =, ~ otherwise Then elc E -v Y. v) _ sc Now select RNC(O,1), find $ such that z r

142 i.e., 2- Zt RA (G 10) Then evaluate sin9 and cos). Alternativelyl6 45 select RN1 and RN2 C (O1). Test: is (2RN2-1) + RN12 < 1? If no, select two new random numbers and test again. If yes, then = / - (G-ll) (A RI - + RN, rsn ~, = -ziz cN2 l -/)R!N, (G-12) (2 RN, -,)Z R/ Z 4. OPTICAL THICKNESS The density function for the optical thickness, K, traveled by a neutron, between successive collision points, is given by the exponential law of attenuation; Pl --- yoy P <,d}e, forO. - 0, otherwise Therefore, Now select RNC (O,1), find M such that R - I - e

143 i.e., = -ln(l-RN). (G-13) Alternatively, the evaluation of the logarithm may be avoided by using a rejection technique devised by von Neumann.45 Upon selecting the optical thickness the actual distance traveled by the neutron, t, is determined from t r= ( <( r('4, R') djR (c-14) where St is given by Eq. (4-52). When the cross section is temperature independent (G-14) can be integrated, otherwise t is found by trial and error. 5. NEUTRON SCATTERING The conditions of scattering are the same as established in Chapter IV, namely, elastic collision, struck nucleus at rest, scattering isotropic in the CMCS. Let Ge = polar angle, CMCS G = polar angle, LCS = azimuthal angle, CMCS = azimuthal angle, LCS E' = neutron energy, precollision E = neutron energy, postcollision. The azimuthal angle is uniformly distributed on the interval (0,2X) for

144 both CMCS and LCS. Select the random variables, cos 0c, sin Gc, cos ~, and sin I as in Section G-3 for the polar and azimuthal angles. For the conditions given above, the scattering angle, G, and neutron energy E can be determined from Gc by the following relations: E = /b (I t * -(1) Cal 8t ) D (G-15) A-I E' costar A+I 4 4 A d Z e Z- D ) A(G-16) S/,. = J 1- oG7 e (G-17) where /A- /) Hydrogen scattering is treated as a special case since a = 0, eoC e = (14 cesec ) D (G-18) Iose VJDG (G-19) But by Eq. (G-8), for the polar angle of scattering, Go, we have RA- =i(/IAoCO)) (G-20)

145 Therefore, (G-21) and Csca e 9 RN (G-22) The evaluation of the square root of RN can be avoided by using the same trick as in Section G-3 (Eq. G-7b), giving cos e= vnax( )N RN2) (G-23) And then E =E'cos (G24) slrg = v/-cos (G-25) The direction of the neutron after the collision, measured with respect to the radial line extending from the collision point r to the axis, can now be determined since the precollision angles, y' and 'I (Section G-3), are known and the scattering angles, 9 and 0 are known, cos = 5/, cos0.S/k Cos cos ' (G-26) sin X = /-dcosz~ (G-27) ~os9 = = / F (SmnerCsocs+ ')CeSV/-(s/ e 99sinr)s/7v' +-cs- e 9ca " (G-28) =- Fs~ (;oqSc. s,- ri '.i (snes,,n)cosv.' (Coes, 5r'),; Y' J (G-29)

APPENDIX H A LIST OF SUPPLEMENTARY COMPUTER PROGRAMS WRITTEN FOR THIS INVESTIGATION In addition to the main digital computer program discussed in Chapter V, several programs were developed for preliminary information and for checking subroutines of the main program. The following programs were written to investigate methods of evaluating the Doppler broadening functions, Psi and Chi, and derivatives, and to study the behavior of these functions with respect to their arguments: 1. evaluation of Psi by the method of Moebius Inversion, 2. evaluation of Chi by the method of Moebius Inversion, e. evaluation of Psi and Chi by the method of Runge-Kutta, 4. evaluation of Psi and Chi by the use of convergent series, 5. evaluation of Psi and Chi by the use of asymptotic series, 6. evaluation of Psi for argument x=O, using the approximation of Hastings62 for the error function, 7. approximation to Psi by extension of Hastings' approximation to the error function with complex argument, 8. check on other approximations to Psi and Chi, 9. evaluation of first and higher order derivatives of Psi and Chi by application of differential equations, Eqs. (F-3) and (F-4), 10. evaluation of first derivatives of Psi and Chi by asymptotic 146

147 series plus empirical functions, and comparison -to results of item 9. The following programs were written in order to study the behavior of the Doppler broadened resonance cross sections of U-238 with varying neutron energy and fuel temperature: 1. calculation of peak cross sections from resonance parameters, 2. e-valuation of practical widths of resonances of U 238 metal and oxide for various fuel temperatures for two cases: 1) neglect of interference scattering, 2) inclusion of interference scattering, 3. comparison of the Doppler broadened cross sections of U-238 with the U02 potential scattering cross section for various energy increments and fuel temperatures. Provision for overlapping of resonances included. The following programs were written to check the subroutines of the main Monte Carlo program: 1. checking of' random number generators and of random variables obtained from the various distributions discussed in Appendix G, 2. generation and listing of random numbers used to initiate each hi story, 3. other auxiliary programs for checking the various subroutines discussed in Chapter V.

REFERENCES 1. M. Born and E. Wolf, Principles of Optics, Pergamon Press, p. XXVI, 1959. 2. G. Breit and E. Wigner, "Capture of Slow Neutrons," Phys. Rev., 49, 519 (1936). 3. H. A. Bethe and G. Placzek, "Resonance Effects in Nuclear Processes Phys. Rev., 51, 450 (1937). 4. E. Creutz, H. Jupnik, T. Snyder, and E. P. Wigner, "Review of Measurements of the Resonance Absorption of Neutrons by Uranium in Bulk," Jo Appl, Phys., 26, 257 (1955) 5. E. P. Wigner, E. Creutz, H. Jupnik, and T. Snyder, "Resonance Absorption of Neutrons by Spheres,," J. Appl. Phys., 26, 260 (1955); Eo Creutz, H. Jupnik, To Snyder, and E. PO Wigner, "Effect of Geometry on Resonance Absorption of Neutrons by Uranium," J. Appl. Phys,, 26, 271 (1955); E. Creutz, H. Jupnik, and E. P. Wigner, "Effect of Temperature on Total Resonance Absorption of Neutrons by Spheres of Uranium Oxide," J. Applo Phys., 26, 276 (1955) 6. E. P. Wigner, '"Results and Theory of Resonance Absorption," ORNL 2309, p. 59 (1956), Confo on Neut. Phys, by Time of Flight, Gatlinburg. 7. J. B. Sampson and JO Chernick, "Resonance Escape Probo in Thermal Reactors," Prog. in Nuclo Energyr Phys. and Math. II, 223 (1958). 8. L. W. Nordheim, t"Theory of Resonance Absorption," GA-638 (1959). 9. L. W, Nordheim, "ResonanceAbsorption of Neutrons," Neutron Physics Conference, Michigan Memorial Phoenix Project, The University of Michigan, June 1961. 10. R. M. Pearce, "The Doppler Effect in Thermal Reactors," Reactor Science, J. of Nuclear Energyo Part A, 13, No. 3/4, January 1961. llo H, Feshback, G. Goertzel and H, Yamanchi, "Estimation of Doppler Effect in Fast Reactors," Nuclear Sci, and Engo,, 1, 4 (1956). 12. Ro To Frost, W. Yo Kato and D. Ko Butler, "Measurement of Doppler Temperature Coefficient in Intermediate and Fast Assemblies," P/ 1777, Proc. Second Intern, Confo Peaceful Uses Atom, Energy (1958). 148

149 REFERENCE (Continued) 13. R. B. Nicholson, "The Doppler Effect in Fast Neutron Reactors," APDA-136 (1960). 14. K. W. Morton, "The Calculation of Resonance Escape Probability by Monte Carlo Methods," P/19, Proc. Second Intern. Conf. Peaceful Uses Atom. Energy (1958). 15. W. H. Arnold, Jr. and R. A. Dannels, "A Monte Carlo Study of the Doppler Effect in U02 Fuel," WCAP-1572 (1960). 16. R. D. Richtmyer, R. Van Norton, A. Wolfe, "Monte Carlo Calculation of Resonance Capture in Reactor Lattices," P/2489, Proc. Second Intern. Conf. Peaceful Uses Atom. Energy (1958). 17. E. Hellstrand, T. Blomberg and S. Horner, "The Temperature Coefficient of the Resonance Integral for Uranium Metal and Oxide," Nuclear Sci. and Eng., 8, 497 (1960). 18. L. W. Nordheim, "A New Calculation of Resonance Integrals," Nuclear Sci. and Eng., 12, 457 (1962). 19. H. A. Risti, G. H. Minton, P. W. Davidson and J. D. Cleary, "Microscopic Lattice Parameters in Single- and Multi-Region Cores: A Comparison of Theory and Experiment," SCAP-1434 (1961). 20. E. J. Leshan, J. R. Burr, M. Temme, R. Morrison, G. T. Thompson and J. R. Triplett, "RBU-Calculation of Reactor History Including Details of Isotropic Concentration. Part I, Method," ASAE-34 (1958); ATL-A-101 (1959). 21. J. Belle, Ed., "Uranium Dioxide: Properties and Nuclear Applications," U.S. Ao.E.C, Gov't, Printing Office, July 1961. 22. A. Keane, "Resonance Absorption in a Slab with a Parabolic Temperature Distribution," AERE-R/M-198 (1958). 23. L. Dresner, "Some Remarks on the Effect of a Non Uniform Temperature Distribution on the Temperature Dependence of Resonarce Absorption," Nuclear Sci. and Eng., 11, 39 (1961). 24. W. E. Lamb, Jr., "Capture of Neutrons by Atoms in a Crystals" Phys. Rev., 25, 190 (1939).

150 REFERENCES (Continued) 25. B. Davison, Neutron Transport Theory, Oxford at the Clarendon Press, London, 1957. 26. R. V. Meghreblian and D. K. Holmes, Reactor Analysis, McGraw-Hill, 1960. 27. G. E. Albert, "A General Theory of Stochastic Estimates of the Neumann Series for the Solution of Certain Fredholm Integral Equations and Related Series," ONRL-1508 (1953); Symposium on Monte Carlo Methods, ed., H. A. Meyer, pp 37-46, Jo Wiley, 1956. 28. J. Spanier, "Monte Carlo Methods and Their Application to Neutron Transport Problems," WAPD-195 (1959). 29. R. V. Churchill, Operational Mathematics, McGraw-Hill, 1958. 30. F. Riesz and B, Nagy, Functional Analysis, Frederick Ungar Publ, Co., 1955. 31. WO Vo Lovitt, Linear Integral Equations, Dover, 1950. 32. F, Go Tricomi, Integral Equations, Interscience Publ., Inc., 1957. 33. P. Halmos, Measure Theory, Van Nostrand, 1950. 34. J. L. Doob, Stochastic Processes, J. Wiley, 1953. 355 K. Knopp, Theory and Application of Inf;inite Series, Hafner Publ. Co., Fourth Edition, 1950, 36. Ivan S. Sokolnokoff, Advanced Calculus, 9 McGraw-Hill, 1939. 37. D.AoSo Fraser, Statistics: An Introduction, J. Wiley, 1958. 38. J. M. Blatt and V. Fo Weisskopf, Theoretical Nuclear Physics, J. Wiley, 1952, 39. H. E. Jackson, L. Mo Bollinger and R. E. Cote, "Capture of Slow Neutrons by Nuclei Bound in Crystals," Phys. Rev. Letters, 6, No. 4, 187 (1961). 40. J. Lo Rosen, J. S0 Desjardins, J. Rainwater and W. WO Havens, Jr., "Slow Neutron Spectroscopy I. U 238," Physo Rev_., 118, 687 (1960).

151 REFERENCES (Continued) 41. H. A. Bethe, "Nuclear Physics: B. Nuclear Dynamics, Theoretical," Rev. Mod. Phys., 9, 69 (1937). 42. L. Dresner, "Resonance Absorption of Neutrons in Nuclear Reactors," ONRL 2659 (1959); Pergamon Press, 1960. 43. A. W. Solbrig, Jr., "Doppler Effect in Neutron Absorption Resonances," Am. J. Phys., 29, 257 (1961); "Doppler Broadening of Low-Energy Resonances," Nuclear Sci. and En., 10, 167 (1961). 44. E. C. Smith, G. S. Pawlicki, P.E.F. Thurlow', G. W, Parker, W. J. Martin, P. M. Lantz and S. Bernstein, "Total Neutron Cross Section of Xe-135 as a Function of Energy," Phys. Rev., 115, 1693 (1959). 45. J. von Neumann, "Various Techniques Used in Connection with Random Digits," Monte Carlo Method, NE.B.S. Applied Math. Series, 12, 36 (1951). 46. H. Kahn, "Applications of Monte Carlo," AECU-3259, Revised April 27, 1956. 47. G. Goertzel and M. Karlos, "Monte Carlo Methods in Transport Problems," Progress in Nuclear Energy, Phys. and Math, II, 315 (1958). 48. E. Cashwell and C. Everett, A Practical Manual on the Monte Carlo Method for Random Walk Problems, Pergamon Press, London, 1959. 49. J. Moshman, "The Generation of Pseudo-Random Numbers on a Decimal Calculator,," J3 Assoc. Comp. Mach., 1 88 (1954). 50. D. L Johnson, "Generating and Testing Pseudo-Random Numbers on the IBM Type 701," Math, Tables Aids Comp., 10, 8 (1956). 51. J. Todd and 0O Taussky Todd, "Generation of Pseudo-Random Numbers," Symposium on Monte Carlo Methods, ed., H, A, Meyer, pp. 15-28, J. Wiley, 1956. 52. Eve Bofinger and V. J. Bofinger, "On a Periodic Property of PseudoRandom Sequences," JO Assoc. Compo Mach,, 5, 261 (1958). 53. R~ R. Coveyou, '"Serial Correlation in the Generation of PseudoRandom Numbers," J. Assoc, Comp. Mach,, 7, 72 (1960).

152 REFERENCES (Concluded) 54. B. Galler, B. Arden and R. Graham, Michigan Algorithm Decoder, Computing Center, The University of Michigan, Feb. 1962. 55. R. M. Pearce, "Radial Dependence of Doppler Effect in Bars of Uranium and Thorium," Reactor Science, J. of Nuclear Energy: Part A, 11, No. 2/4, Feb. 1960. 56. A. Erdelyi, Tables of Integral Transforms, Vol. I, II, Bateman Manuscript Project, McGraw-Hill, 1954. 57. R. K. Osborn, "Some Characteristics of the Thermal Neutron Scattering Probability," Nuclear Sci. and Eng., 3, 29 (1957). 58. R. R. Goldberg and R. S. Varga, "Moebius Inversion of Fourier Transforms," Duke Math. J., 24, 553 (1959). 59. R. J. Duffin, "Representation of Fourier Integrals as Sums, Parts I, II, III," Bull. Amer. Math. Soc., 51, 383 (1941); Proc. Amer. Math. Soc., 1, 250 (1950); Proc. Amer. Math. Soc., 8, 272 (1957). 60. D. Shanks, "Non-Linear Transformation of Divergent and Slowly Convergent Sequences," J. Math. and Phys., 34, 1 (1955). 61. I. I. Hirschmann and D. V. Widder, The Convolution Transform, Princeton Univ. Press, 1955. 62. C. B. Hastings, Approximations for Digital Computers, Princeton Univ. Press, 1955.

UNIVERSITY OF MICHIGAN 3 9015 03483 21731111111 3 9015 03483 2173