ENGINEERING RESEARCH INSTITUTE THE UNIVERSITY OF MICHIGAN ANN ARBOR DIFFUSION OF NEUTRONS IN A HEAVY ELEMENTS MEDIUM AND APPLICATION TO THE MULTIGROUP DIFFUSION THEORY by,, E. Nasjleti Approved byi H. A. Ohlgren Project 2505 THE CHRYSLER CORPORATION DETROIT, MICHIGAN April, 1957

TABLE OF CONTENTS Page 1.0 INTRODUCTION 1 2.0 SCATTERING OF NEUTRONS WITH HEAVY NUCLEI 3 3.0 SLOWING DOWN OF NEUTRONS BY INELASTIC SCATTERING IN A HEAVY ELEMENT MEDIUM 8 4.0 MULTIGROUP DIFFUSION THEORY 16 APPENDIX I 21 APPENDIX II 22 LIST OF TABLES Page TABLE I 14 TABLE II 15

DIFFUSION OF NEUTRONS IN A HEAVY ELEMENTS MEDIUM AND APPLICATION TO THE MULTIGROUP DIFFUSION THEORY 1.0 INTRODUCTION The main nuclear interactions which occur during the diffusion of neutrons with energy between 0.1 mev and 10 mev in a heavy medium (A> 100) are elastic and inelastic scattering, capture and fission. Fast reactors are mediums of this type. The neutron spectrum "obtained from different reactors shows a broad maximum about 0.1 mev and extends up to the 3-5 mev region. The fission cross section of U235 and p239 are almost constant in this energy range with an average value about 1.5 barns. The U238 fission cross section has a threshold at about 1 mev'rising to about 0.55 barns so that it assumes considerable importance in the energy range under discussion. The capture cross sections (n,y) are rather small except for U238, and are well known experimentally. In this energy range, the inelastic scattering is the major mode of the neutron energy degradation'and is, as a consequence, one of the most important nuclear parameters. Elastic scattering is no longer an important source of neutrons degradation, hence its importance is lessened. The multigroup diffusion method has shown to give fairly accurate results for the flux distribution and critical mass of a fast assembly, so long as the radius is reasonably large. * * P/609 (International Conf.-Geneva) A survey of the theoretical and experimental aspects of the fast reactor physics. (ANL) -1-~

For small assemblies diffusion theory is less accurate since the neutrons mean free path is not small with respect to the geometrical dimensions of the assembly. In cases where more accurate results are required, the Sn method developed by B. Carlson * can be used. This method retains the angular dependence of the neutron flux. It has been programmed for the IBM-701 with S-4 in a multigroup form using spherical geometry using both three and ten energy groups.* For small assemblies (i.e., EoB.RoI where r = 10cm), the diffusion method gives a critical radius which is 9% higher than that given by the S4 calculations. For larger assemblies, the critical radii results of the two methods agree. In solving the multigroup diffusion problems, a numerical method can be used for the integration of the equations. Knolls Atomic Power Laboratory has provided a program for the IBM-650 automatic computing machine.** In order to use this program to solve a specific reactor design, group parameter involving the various cross sections must be computed. Using the inelastic slowing down treatment developed in this report, the inelastic scattering removal probabilities for U235 and p239 are computed. The numerical method used in the KAPL program, which is the same presented by Ehrlich and Hurwitz in Nucleonics, 54, is discussed in the final part of the report. i LA-1891 B. Carlson, Solution of the transport equation by Sn approximation. * KAPL-1415 One space dimensional multigroup for the IBM-650. Part I, Equations. KAPL-1531 One-space dimensional multigroup for the IBM-650. Part II, Machine program. -2

2.0 SCATTERING OF NEUTRONS WITH HEAVY NUCLEI 2.1 Elastic Scattering The elastic scattering is a reaction of the type n + ZNA -ZNA + n It is important to know the angular dependence. At low energies ( 0.1 mev) the angular distribution of the elastic scattered neutrons in U238, u235, p239 and others with A> 100, shows a single maximum in the forward direction. As the energy of the incoming neutron increases, the peak in the forward direction is more pronounced, but there also appears a secondary maximum.* As the energy of the incident neutrons gets higher than E = 2.5 mev, these secondary maxima are displaced toward smaller angles. With use of the continuum model, which requires that the spacing of the energy levels in the compound nucleus to be small compared to their widths, the angular distribution can be predicted and is in agreement with experimental data in the prediction of the strong for-ward peak shape. However, this method fails to explain the fact experiment- ally found that the total cross section of intermediate and heavy - elements measured as a function of the neutron energy shows a broad maxima and minima (between 0.5 > cos Q> 1) which shift slowly with atomic weight. To interpret this experiment, the strong interaction hypothesis was replaced by the assumption that the incident neutron interacts with * M. WaIt & H. Barschall, Phys. Rev., 93:1062 (1954) P/588 Angular distributions and non-elastic neutron scattering. (Los Alamos) (Internation Conf. - Geneva) -3

the average potential produced by the other nucleons.* The assumed potential is a complex square well of the form V(r) = -Vo (l+i~) r<R V(r) =0 r>R where R =.45 A1/3 x 103 cm; Vo = 19 mev and $ = 0.05. The calculated total cross sections are in good agreement with the experimental values for neutron energies up to 3 mev for heavy nuclei. If the differential elastic cross section for scattering through an angle 0 is known, the inelastic scattering and transport cross sections can be obtained from the total cross section. These cross sections are given by the following expressions gin = aT - jf(Q)dw atr = aT - Ja(Q) cos Qdw dw solid angle between Q and Q+dQ C(Q) diff. elastic cross section through Q aT total cross section. It can be said that a(@) is a parameter rather accesible by experimental techniques. In the last years, measurements have been performed for many nuclei and for several neutrons energies up to 14 mev.** * H. Feshbach & V. F. Weisskopf, Phys. Rev., 90:166 (1953) ** M. Walt & H. Barschall, Phys. Rev., 93:1062 (1954) P/588 Angular distributions and non-elastic neutron scattering. (Los Alamos) (International Conf. - Geneva) Interaction of n's (1.0;1.77;2.5;3.25 & 7 mev) with Nucl.-Phys. Rev., 104:1319 (1956),4

2.2 Inelastic Scattering The inelastic scattering is a nuclear reaction of the type A A n + NZ - NZ + n + 7 The inelastic scattering becomes important when the incident neutron energy is of the same order or higher than the first excited state of the target nucleus. For incident neutron energies below this threshold, the (n,y) is the more probable process. At neutron energies between 10 mev>E>0.2 kev for nuclei as U235, U238 and p239 the inelastic scattering is the more probable process. Nuclear reactions with proton or a emission will be very rare for these elements. The use of the statistical theory of nuclear reactions permit the computation of the inelastic scattering cross section and the prediction of the energy distribution function of the outgoing neutron. The cross section for a a(n,n) inelastic process may be written o(n,n) = an 1n where an is the cross section for the formation of the compound nucleus and In is the relative probability that a neutron will be emitted, and is given by Fn Bn = E c c n is the partial width of the compound state for the emission of a neutron, averaged over the compound states which are excited by the incoming neutron. The sum in the denominator should be extended over all the particles which are ejected. _-_

The application of the statistical model for the compound nucleus permit drawing the important conclusion for reactor physics that the inelastic scattering is spherically symmetric.* The distribution function obtained for the (n,n) reaction is** f(E'-E)= a E, exp (-r E) where a is a numerical parameter that is adjusted from experiment, Eis the incoming neutron enegy, and is the outgoing neutron energy. 2.3 The (n,2n) Reaction When the neutron energy is above many mev (E> 2.5 mev) the (n,2n) reaction becomes relevant. In this case, the residual nucleus is left after the emission of the first neutron in an excited state with an excitation greater than the binding energy of a neutron. Therefore, in absence of competition from other modes of decay, the emission of a second neutron is then possible. The cross section for this process is o(n,2n) = a f f(E/E) dE where cn = E-Et Et = threshold energy for the (n,2n) process. The energy distribution function of the outgoing neutrons (including the first and second neutron) is given by the following expression* NYO-636 Final report of the fast neutron data project (NDA) **Phys. Rev., 52,295 (1937) -6

F(E'- E) = T (x e -+2 exp - i )) + x ((-l)e y dy 1i) x where x = E; = x-xb = (E ) Eb = energy required to remove a neutron from the nucleus. Experimental information about inelastic scattering cross sections can be found in recent publications.* Theoretical values can also be obtained using the statistical model.** The (n,2n) reaction's cross sections for reactor material are quite small for the energies between 0.1<E <3 mev (fast reactors). Hence, it can be neglected in fast reactor calculations without introducing an appreciable error. * M. Walt & H. Barschall, Phys. Rev., 93:1062 (1954) P/588 Angular distributions and non-elastic neutron scattering. (Los Alamos) (International Conf. - Geneva) Interaction of n's (1.0;1.77;2.5;3.25 and 7 mev) with Nucl. Phys. Rev., 104:1319 (1956) ** NYO-636 Final report of the fast neutron data project (NDA) -7

3.0 SLOWING DOWN OF NEUTRONS BY INELASTIC SCATTERING IN A HEAVY ELEMENT MEDIUM 3.1 Neutron's Energy Loss After an inelastic scattering the energy of the neutron is reduced by the amount of the excitation energy of the residual nucleus and by the recoil energy loss. This last is exceedingly small for heavy nuclei. Let 0(u') be the neutron flux at lethargy u' and ZAs (u') the corresponding macroscopic inelastic scattering cross section. The number of inelastic scattering collisions per cc per sec experienced by the neutrons in the lethargy element du' are E. (ut )(ut )du The number of neutrons that will emerge in the lethargy element du will be Eis(u')0(u')fn(u' ->u)du' du where fn(u' - u) is the normalized distribution function. (See Appendix 1) If the lethargy range is divided into intervals Aui (i = 1,2,3,...,n) the number of neutrons that will emerge in du originated by inelastic scattering in Au' will be f Eis(U')(u') fn(u' -4u) du' du (3.1.1) Aut The following average values will be defined Au' Au' f (u')du' =<.>,. Au Au' Au' -8

For slow varying functions in an interval, it can be assumed that the average of the product is the product of the average, and making use of the mean value theroem for integrals, expression (3.1.1) can be written < Ais ~u',' fn (u' -4u) du' du Hence, the total number of neutrons that will emerge in the interval Au will be <iaU,' <>u'<, au,,t fn (u'-u) du' du = S(u'-u) (3.1.2) The double integral can be analytically computed. If the index "i" is associated with the group Au', and "n" with the group Au, expression (3.1.2) can be written S(in n) is. <>i P(i- n) where P(i- n) = fi f n f(i-n) du du' (3.1.3) P(i- n) can also be expressed as P(i- n) = p(i- n) Ai from where p(io-n) =P(i- n) fi fn fn(i- n) du du' p(i-4n) (3.1.4) p(i- n) can be interpreted as the average probability that a neutron which made an inelastic scattering in the group "i" will emerge in the group "n." Expression (3.1.2) as is written, can be introduced in the diffusion equation as the inelastic scattering source term. -9

3.2 The p(i- n) Function The normalized inelastic scattering distribution function is* E exp E dE fn(u'- u) du = -a, 1-+ aE') exp - E} where the normalization factor Ci =1 ( + aET) exp {.IT) (see Appendix 1) will be slowly varying and nearby unit for most of the relevant energies because of the rather large value of the parameter "a" for U235, p239 and i238, and can be treated as a constant. Hence, from expression (3.1.3) is obtained aE P(i- n) = (-1)Ci 1i fn E exp - E dE du' (3.2.1) The index i and n corresponds to the intervals Aui and Aun of the lethargy range. U. and ui are the extremes of the interval Aui, hence, Aui =Ui - ui u > u. Corresponding in the energy range will be + X+ ui -E1 u. E u7 - E7 1 1 therefore E1 > E+ it is also true that + _ u = Ui+l and correspondingly + - Ei = Ei+l E. U. I t path of integration Ei ui * NYO-636 Final report of the fast neutron data project. (NDA) AERE-T/R-1500 Ms. M. E. Mandl, Multigroup theory with an application to the inelastic scattering in uranium. -10

let f + r (3.2.2) and x = E (3.2.3) from the last expression dx = ~ dE, expression (3.2.1) becomes P(i ->N) - (-1)Ci.i j x exp 1-x. d u' (-l)Ci Ji -x exp (-x} - exp {-x] dut n therefore P(i -n) Ci fi (X+ exp -X3+ exp -x ) - (3.2.4) -(Xn exp &xn; + exp xnJ )J du Because of the symmetry of the expression, the integration over the interval "i" will be carried out in detail only for the x4 part of the above expression. Integrating over interval "i," (u ), two integrals are immediately obtained, namely I= f i x+ exp -x du' I1 = ~i n and =+ =f exp -n du' 2 i { u' setting du' = d(lg Eo dE and from (3.2.2) and (3.2.3) 3/2 d = - (-1/2) ra E' dE' (E' corresponds to Ei) Integral Il becomes 1'2E^;' exp E I = (2 - exp { -En 11= (2) exp { En7- (2) exp {- En E -11

In a similar manner and using relations (3.2.2) and (3.2.3) as before, the following can be obtained for I1 I i E exp i -E+n dEi = (2) exp E I+ (2) E + (2) E E E i i where Ei(-x) is the exponential function integral defined for negative values of the argument by the following expression o0 t -Ei(-x) = dt x -t For values of the argument larger than x>10, the exponential integral function defined above can be approximated by the following asymptotic expansion ( )_ x)l =. K1 1 2... 3 This function is tabulated for positive and negative values of the argument. Hence, the expression regarding x+ in (3.2.4) becomes [ + I - (2) exp n - (2) exp En + (2) Ei {En - (2) Ei{ E n i n (2) exp {E (2) E -E - (2) exp {- En 3 (2) Ei En If it is defined a function B(Ei -E )= (2) exp "En - (2) Ei - En E -12

-13therefore I + I - E - E) - B(Ei - n El Using the above results, expression (3.2.4) can now be written P(i -*n) = C + 1 - I - 1 (3.2.5) Ci B, E+ >.. B(E+ -+E) = C LB(EI - En) - B(E 4 En) i n i - B(E ->E-) + B(En )E] (See Table I) Remembering expression (3.1.4) p(i - n) will be given p(i - n) =(i - n) (See Table II) (3.2.6) u. - u'.' Ui 1i The normalization constant Ci can be computed as an average constant over the group i. In the simplest approach may be used E- Ei+ 2 The function B(Ei - En) has a useful property, namely B(ti+2Eo -_ tn+lEo) = B(tiEo-tnEo) (See Appendix II) This property simplifies calculations in multigroup diffusion problems if the lethargy range is divided into equal intervals. 3.3 Inelastic Scattering Slowing Down The number of neutrons which are slowed down from i to group n, using expressions (3.1.2), (3.1.3) and (3.1.4), can be expressed as S(i ->n) = iS) i p(i -n) A. i-4n if it is set 4s>)' p(i -rn) = qs S(i -n) = -Zisn<>i. (3.3.1) This expression can be used in the multigroup diffusion theory as the source term of the diffusion equation for the group "n."

TABLE I P(i- n) for U235 and P239 (a =13.4 mev-l) i intervals P(i- n) =Ci B(E En)-B(E En)+B-(EtE )-B(Ei- E)E * ui u. n=i+l n=i+2 n=i+3 n=i+4 n=i+5 n=i+6 1 0 1.5.3176.6755.1732.088.069.0012 2 1.5 2.249.1032.0556.0264.0176 3 2 3.288.2392.1283.0145 4 3 3.5.123.2217.0439 5 3-5 4.284.0685 6 4 5.371 7 5 6.3087 8 6 7.2955 9 7 8.56 10 8 oo * The value of "a" was obtained from NY0-636 from a graph for odd-nuclei and correspond to U235. The graph was made up using experimental values of "a" for many nuclei. Between the heavy elements only the value for Th was experimental. The extrapolated "a" for U235 and UJ239 differ so slightly and because of the uncertainty in their values, the P(i-*n) obtained above can be extended to Pu239 ~ B(Et-En) = (2) exp - Enf a - (2) EilEn J Ci constant averaged over group "i" -14

TABLE II p(i- n) for U35 and Pu239 lp (i-n) P(i- n) p(i->n) U -+U1 i ui Ui n=i+l n=i+2 n=i+3 n=i+4 n=i+5 n=i+6 1 o21173.45033.11546.05866.o46.oo0008 2.498.2064.11120.0528.0352 3.288.2392.1283.0145 4.246.4434.0878 5.568.1370 6.371 7.3087 8.2955 9.56 10 -15

4.0 MULTIGROUP DIFFUSION THEORY 4.1 The Diffusion Equation for a Heavy Element Medium (Fast Reactor) The processes which have to be considered are fission, radiative capture, inelastic scattering and elastic scattering. The neutron slowing down is attributed only to the inelastic process. If i(u,r) is the neutron flux at position r, and lethargy u the diffusion equation can be written as 3Zt V2 0(r,u) + (a + Eis) s (u,r) Qis (r.u' u) + (4.1.1) + Qf(ru) where Za(u) - Macroscopic absorption cross section included fission. Zis(u) - Macroscopic inelastic scattering cross section. D~r(u) - Macroscopic transport cross section. Qis(r,u'- u) - Inelastic scattering source at r. Qf(r,u) - Fission source at r. Can be set Qis =u is(u) i (r,u') fn(u - u) du' and Qf = vX(u) (u ) Et (u) dut +,6(k) Z(k] where k - last lethargy group and f, X(u) du = 1 A new function for symmetrical cases can be introduced, namely O(u,r) = r.o(u,r) and setting Sf(u,r) = rQf -16

the diffusion equation becomes in the particular case of spherical symmetry o2(D (u _ i -D(u) @2 u r + [Za(u) + Zis(u)] O(u,r) - (4.1.2) - f is(ut) * (u',r) fn(ut'4n) du' _ Sf(r,u) 4.2 Multigroup Approximation In this part, the procedure given by Ehrlich and Hurwitz, will be closely followed for converting equation (4.1.2) into a set of equations, one for each group.* The variable u will be divided into groups of width Au. Au=u n+ -u (uu) n n -n U (u > n) un = Un+l + un corresponds to the incoming neutrons and un to the outgoing neutrons of the group n. After defining f un (r,u)du =un. AV(r,u) and integrating equation (4.1.2) over Aui,. the following equation is obtained ( ) nAyv un is + - ) nAv = dr nAv i - = Z ^T f f iUn s (u')t (r,ut). fn(u'-u) du du' + n=l i n + u' Aui Zis(u') $ (r,u') fn(U'-u) du du' + fAn S(r,u) du. * R. Ehrlich and. H. Hurwitz, Jr., Multigroup methods for neutron diffusion problems. Nucleonics, page 23, Feb. 1954. -17

Making use of expression (3.3.1) results -AUn (D r nAv n + 4 nAv i-1 = E S(i- n-1) + S(n-n) + AunS(r,u) du i=ln 4.3 Numerical Method Expressions similar to (4.2.1) can be written for each group n in which the lethargy range is divided. The equation will be treated by a simple finite difference method. The spatial region will be divided into sub-intervals Arm and the index m will denote the radius rm at the mth space point. With the use of a simple difference approximation, the second derivative can be written (d2FM) F= +l + Fm_ -, 2Fm dr2 r ( 2 setting this equation into relation (4.2.1) will be found the recurrent relation for the Fn s. Fn n n n (4.3.1) m+l = kn F - Fm-1 - I (31) where kn 2 + an(Za +Zis)n - (is)n P(n-4n) (^ Dn - n-1 (Lrm)2 n-l Z (Zis)t P(t-n) Ft + Rn S(rm,') du m aDn a -n n (Arm) It is convenient to choose Ar so that the boundaries will occur at space points. The boundaries conditions for the current into the medium have to be modified so that the conservation of neutrons is not altered by the substitution of the differential equations by difference equations. Hence, the expression used for the neutrons current is then uniquely determined by the spatial integration. -18

The formula that will be used is r2 / ^I~f_'' dr =r = m rm Ar -( r + r2 + ~ 0r2 - r 1r) rl m-rl m=r1 This formula is rigorous if f is a linear function of r between rl and r2 since r2=(-~+ r r2 r r r (re O)drr ) r1 The neutron current can be evaluated by the l.h.s. of the above equation, and is obtained 2 rm n n 2 n + F1 ( 4.3.2) rm =-Dn (Fm+l - FMo') 3' (F+l + Fn ) (4 3 2) Equation (4.3.1) will be solved by R. H. Stark's method that is convenient for solving the equation on fixed point automatic computing machines. Two new variables will be introduced Am+l m Ij ^1-^ ancL ^ A,^ (4.3.3):~m+1- Am and.m = Z Aj (433) C6+1 Am d~mj=o J Am am and. m satisfies the following recursion relations 1 am-1 am+l = kn - and mm = + I Proof will now be given that for any linear boundary conditions which involve Fm-.i Fm and Fm+l can be reduced by the use of equation (4.3.3) to a form involving only two values of F.* Assume that Fm m+l = Fm+l + ~m then m FmAm+1 = Fm+ Am + Z Aj I j=o m (kn Fm - Fm1- Im) Am + E A Ij j=o * AERE-T/R-1500 Mrs. M. E. Mandl, Multigroup theory with an application to the inelastic scattering in uranium. -19

from where the following expression can be obtained m-l Fm_1 Am = Fm Amil + Z Aj Ij j=o or Fm1 am = Fm + 3m-1 Hence, by induction, Fm can be calculated by the following recurrent relation Fm = m+l m (4.34) 4.4 Boundary Conditions At the mesh point rn, the physical conditions can be approximated using any one of the following ways a. Fi = O m b. - rm 2r -l m m rm+j m-1 6) 1I 3Ztr Case a. - making use of (4.3.1)and (4.3.4), and conditiona' is obtained. i im+2 =ki m+l = Im+l Case b. - using (4.3.1) and condition'b'is obtained. ~i+l = rr Dmi [im Irm 2 i -r m m 3+ -20

APPENDIX I Normalization of f(E-> E) f(E'- E) = aE aE therefore fn(E' E) = f(E'-E) (I.1) Ef (E' I- E)dE and fo fn(Et- E)dE = 1 introducing a new variable E - = dx = E dE The integral in the denominator of expression (I.1) will be called I and after integration by parts there is obtained I = f x exp {-x} dx =1 - (1 + aE') exp - {aE from where E exp I-Jr 5E} dE f (E'tE)dE = a E exp - E dE fn E'E 1-(1 + IaE') exp \- a'} -21

APPENDIX II B(tn+2 Eo - tJ+l Eo) = = (2) exp {t t Eo itn2 E- i {o H } t~+ atn Fo E hence B(+n+2 Eo -ttj+l Eo) = B(tn Eo -Atj Eo) -22