E UNVERSIT OF MICHIGAN INDUSTRY PROGRAM OF TEE COLLEGE OF ENGINEERING ON THEE TEORY OF XENON INDUTCED INSTABILITIES IN NEUTRON, FLJX DISTR-IBUTION Geza Lo.- qorey A dissert~,ation submitted in partial fulfillmeint of the requirements for the degree of Doctor of Philosoplq in the University of Michigan 1960 June, 1960

TABLE OF CONTENTS Page LIST OF TABLES.............. o o o o o o o o o o o o..... o o... o o o o... iii LIST OF FIGURES........................................... iv NOMECLATURE o................. 000..................... v INTRODUCTION.................... 1 CHAPTER I PRELIMINARY INVESTIGATION OF AN ELEMENTARY REACTOR MODEL o......... o o o o o o o o o o o o. o. o o II THE FORMULATION OF THE PROBLEM..o...... 11 III THE STEADY STATE oo.o...oo..ooooo....... 16 IV THE EXPANSION OF THE EQUATIONS OF MOTION.....o.....o. 28 V THE DIMENSIONLESS FORM OF THE EQUATIONS OF MOTION oooo 32 VI THE MODAL COUPLING COEFFICIENTS....o...... oYoo.. O o... 36 VII THE DERIVATION OF THE STABILITY CRITERION WHEN THE EFFECTS OF THE NONLINEARITIES AND THOSE OF MODAL INTERACTION' ARE NEGLECTED................. 43 VIII THE EFFECT OF THE NONLINEARITIES....oo...oo o. 53 IX THE EFFECT OF MODAL INTERACTION o....o. o.......OOOOuOO 66 CONCLUSION..o.......................... 77 REFERENCES APPENDIX ~~0000................................. 79 BIBLIOGRAPHY...........................0........0..i......... 81 ii

Table Page 1 Maximum Steady State Flux Flattening Produced by Xenon-135 as a Function of Reactor Size..........00..O 27 2 Value s of 8 }T2 *... 0 a 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a 0 0 4~ 4 Values of 3-, 0 o a o a Q. 0 Values of H/M for which Sn 1 0.0000000000000-000000000 50 Values of S for H/M100 00 0000000. 50 iii

LIST OF FIGURES Figure Page 1 The Behavior of the Modal Coupling Coefficients Fij and Pij as Functions of_~............. 42 2 The Behavior of n (ePl _/) as a Function //l x +,n of o a * a a o a a e 0 a o a o ao o o oa 0 47 3 Solution Paths of the Second Harmonic for S2 = Jo9O0o 56 4 Solution Paths of the Second Harmonic for S2 = 180o 56 5 Periodic Motion of the Second Harmonic for S2 = o180..O......... ~ ~ a ~ ~ a a a 0 57 6 Limit Cycles of the Third Harmonic With and Without the Nonlinear Terms for S3 =.170..oo...oooooo.o 60 7 Solution Paths of the Third Harmonic for S3 = o170o 60 8 Periodic Motion of the Third Harmonic for S3 - 170.... a. o o D. o o. a a o a a a 00o ooe o o, o o o ~ o o oa a o o o o 62 9 Limt Cycles of the Third Harmonic or S3 1760 63 10 Limit Cycles of the Third Harmonic for S3 = o1760 o 63 10 Limit Cycles of the Third Harmonic for S = o.765. 63 3 11 Solution Paths of the Third Harmonic for S3 = o1770 ~..,.., 64 S 3e ~ 1770...o.................... o o 64 12 The Behavior of F///) V and of Fo _ as Functions of L a.... 72 13 The Effect of Modal Interaction of the Critical Value of S2 ~o-o................................U. o 73 14 Oscillations of N2, N4, and N6 for S2 = S4 = S6 = 409 o. o o o o o o o o o o 75 15 Oscillations of N2, N4, and N6 for S2 =.409, S4 = S6 =1 o.......... o....... a.... 76 iv

NOMECIATURE ajCog ) the time-dependent part of the i'th term in the expansion of!' (/t ) same as above in the expansion of /X(r/') same as above in the expansion of 5 (Pt ) h index index J' index the (infinite) multiplication factor of the reactor the mean neutron lifetime in an infinite medium J41 index index the resonance escape probability 9f the slowing down kernel r space variable time variable L7 the average speed of thermal neutrons X space variable A (iry) the absorption rate of neutrons by xenon-135 the material buckling of the reactor C the criticality factor of the reactor Cn the criticality factor of the n'th mode of the neutron flux distribution in the reactor -- Cn the effect of the xenon-135 poison on the criticality factor of the n'th mode ~scn ~ the excess criticality factor of the n'th mode (Ckl 1 ) v

D the thermal diffusion constant ~F(-, e) the fission rate F(IY) the Fourier transform of i (r) TV the modal coupling coefficient associated with the steady state neutron flux distribution F! the thickness of the slab reactor the concentration of iodine-135 io (v) and jC() the steady state iodine-135 distribution fIn (t) the reduced amplitude of the ntth mode in the expansion of 1iY(',t/) ~c(r, /) the deviation of the iodine-135 distribution from the steady state L the thermal diffusion length M/4 the migration length!%/h (~) the reduced amplitude of the n'th mode in the expansion of 6(rt,/) P:,j the modal coupling coefficient associated with the steady state xenon-135 distribution CQ et) the fast leakage escape probability of the n'th mode the reduced subcriticality of the n'th mode of the neutron flux distribution T the dimensionless time variable XCt,() the concentration of xenon-135 XoP_ ) and X(p) the steady state xenon-135 distribution Xn (Lid the reduced amplitude of the n'th mode in the expansion of'Xb / f) vi

~S(F/, X ) the deviation of the xenon-135 distribution from the steady state the fission yield of xenon-135 the fission yield of iodine-135 the fast fission factor Ax the decay constant of xenon-135 the decay constant of iodine-135 A - /\X / )1 the i'th characteristic value of the characteristic value problem underlying equation (22) the average number of neutrons produced upon fission the microscopic thermal neutron absorption cross section of xenon-135 Tr the neutron generation time in an infinite medium the effective neutron generation time of the n'th mode of the neutron flux distribution # (L~tv) the thermal neutron flux density -4 (k) and 4(P) the steady state thermal neutron distribution $) (k ir ) the deviation of the thermal neutron flux distribution from the steady state ~(fn ~_~ the n'th characteristic function of the wave equation ___ the variable of the Fourier transform vii

1n the maximum effect of the xenon-135 poison on the criticality factor of the w'th mode vJ ~n ~the n'th characteristic value of the wave equation AI the slowing down length the thermal absorption cross section the thermal fission cross section i_/ the thermal neutron flux distribution in the clean reactor h, the nonlinear modal interaction coefficient iJ, _5_- the reduced maximum steady state thermal neutron flux viii

INTRODUCTION Ever since the advent of nuclear reactors, the problem of reactor stability has been an important one. Although a great deal of time and effort has been expended on the study of nuclear reactor kinetics and stability, the knowledge of the subject is still far from satisfactory. The great majority of the published literature on the subject deals with the so-called space independent reactor kinetics in which it is assumed that the spatial and time variations of the neutron flux or of the reactor specific power can be treated separately. In the last few years, however, coinciding with the coming of age of large power reactors, there has been an increasing interest in space dependent nuclear reactor kinetics. One can, perhaps somewhat arbitrarily, classify types of reactor kinetic behavior according to the magnitudes of the time constants which dominate the behavior. According to such a classification one might speak about: 1. Short-term kinetics with time constants of the order of less than a few minutes such as behavior due to temperature effects, 2. Intermediate term kinetics with time constants of the order of hours such as the effects of the buildup and depletion of xenon-135 and samarium- 149, 3. Long-term kinetics dealing with problems such as fuel burnup and reactor lifetime with time constants greater than several days. -1

-2One might also classify reactor kinetic behavior in a different manner, and. speak about behavior induced by internal and external processes. The former refers to behavior induced by processes occurring within the core of the reactor such as poisoning and temperature effects, and the latter takes into account the effects of equipment associated with the reactor core such as the control and heat removal systems. At this point one must hasten to add that the types of behavior enumerated above are not distinct, but are interdependent to a certain degree. For this reason. the above classification of reactor kinetic behaviors may not be completely meaningful, but is certainly convenient. This thesis deals with the intermediate term, internal, spaceand time- dependent reactor stability problem of neutron flux shape variations due to the effect of xenon-135. The treatment is that of an intermediate-term reactor kinetic problem in that it is implicity assumed that the short-term behavior of the reactor is stable. The problem studied is internal in the sense that the study investigates the inherent stability of the spatial distribution of the neutron flux in a reactor, that is, it attempts to answer the question,'"is the shape of the neutron flux distribution self-correcting for perturbations from a steady state without interference from a control system?" Various treatments of the above mentioned problem have appeared in the literature (see references five through nine). The purpose of this thesis is to extend the treatment of the subject~ In the light of the amount of knowledge still necessary before satisfactory understanding of the problem is achieved, the forward step taken by this work is a

-3small one. The treatment is restricted to bare, homogeneous, solid fueled, thermal reactors- the word "homogeneous"' here refers to the cold clean reactor. The justification for restricting the treatment to solid fueled, thermal reactors is that they are the ones expected to suffer most from xenon trouibles. The rest of the above mentioned restrictions are justified by the fact that the asymptotic theory of homogeneous reactors is basic to all of reactor theory. There is one more major restriction placed on the problem.treated. It is assumed that the variation in the xenon-135 density is the only feedback (the word "feedback" here refers to any process by which a change in the neutron fl-ux affects the properties of the medium).effect influencing variations in the shape of the neutron flux distribution, the effects of other feedback processes, such as those of reactor temperature, are disregarded. The problem is treated using the method of harmonics, the spatial distribution of the neutron flux is expanded into an infinite set of orthogonal functions: the harmonics or modes of the distribution. By this expansion, the finite set of nonlinear partial differential equations describing the problem are replaced by an infinite set of ordinary nonlinear differ ential equat:ions. These equations are interdependent, so that the time behavior of each mode is affected by the time behaviors of all the other modes. In the treatments which have appeared in. the literature, the criterion for the stability of the shape of the neutron flux distribution is derived after making two major assumptions. These assumptions are:

-41i. That the interaction among the modes of the flux distribution is negligible, and therefore one can examine the stability of each mode by itself, and 2. That the stability criterion can be derived neglecting the nonlinearitieso The latter assumption requires some comments. In the case of a nonlinear system, one can, in addition to absolute (or nonlinear) stability, speak about linear stability, provided that the system has a linear approximation in a nonzero region about the steady state. Linear stability in this case means that the system is self-correcting to small disturbances about the steady state. For a linear system, of course, linear stability is synonymous with absolute stability. If stability means that disturbances die out, then in the case of a nonlinear system linear instability implies absolute instability, and absolute stability implies linear stability, but not conversely. In other words, linear stability is a necessary condition for absolute stability. Whether it is also a sufficient condition depends on the nature of the nonlinearitieso It is quite possible that the nonlinearities make a positive contribution to stability, so that linear stability is a sufficient condition. Nonlinearities can, however, also detract from stability and make the system unstable in spite of its linear stability. It must also be pointed out that an inherently stable system can be made unstable by tieing it to an improperly designed controller, just as an. inherently unstable one can be stabilized by a suitable control system. The purpose of this thesis can now be restated in more specific termso The purpose is~

- 51o To rederive the stability criterion when the effects of the nonlinearities and those of modal interaction are disregarded, from a formulation of the problem which includes the slowing down kernel and the neutron generation time, and to cast this criterion into a convenient dimensionless form; and 2. To investigate the effects of the nonlinearities and those of modal interaction on the stability. It is very important that the reactor designer have reliable information of the stability of the shape of the neutron flux distribution in the reactor, because if this is inherently unstable, he will be called upon to design a suitable instrumentation and control system to eliminate the instability. Therefore one must have some knowledge of the effects of the approximations involved in the analysis on the final predictions of stability, The original contribution of this thesis to the field of reactor kinetics is the investigation of the effects of the nonlinearities and those of modal interaction on the stability of the shape of the neutron flux distribution.

CHAPTER I PRELIMINARY INVESTIGATION OF AN EI1EMENTARY REACTOR MODEL Before embarking upon the formulation and treatment of the problem of xenon induced, space- and time-dependent neutron flux variations, it is instructive to examine an elementary nuclear reactor model which, in spite of its simplicity, gives a hint at the possibility of spatial instabilitiLes in the neutron flux. A hypothetical reactor will be examined in which the neutron flux can be represented by the one energy group equation r mtra hi~r,t) = D D (r,) + ( 92+-Ha) for +) (1) with the boundary condition that 5 vanish on the extrapolated boundary of the reactor, In this equation cr is some neutron speed, V is the average number of neutrons produced upon fission, D is the diffusion coefficien.t and,4 and - are the macroscopic fission and absorption cross section respectively. + can be expanded as follows: (2) where A is a solution of the characteristic value problem with the same boundary condition as the one that applies to equation (W-6

-7Furthermore, the's form. a complete, orthonormal set so that EA c rOe VOL 6/M1 One can now expand equation (1) into a set of ordinary differential equations by multiplying it by / and integrating over the volume of the reactor~ The resulting equations are: Defining the following quantities~ the set of equations can be written as 7-a, = (11 +Mao Z %4 @ _X2 A- a- (5) or as e, Iv/Z, 2 (8f2 _LL/rn2 lea) n, =,2, ~~(6) These equations describe the time behavior of all the modes, or harmonicas, of the neutron flux distribution i.n the hypothetical reactoro Now if on.e numbers the characteristic values such that do, < dep2 4 d e 4 -, i then the familiar resu.lts are apparent: the condition. to be satisfied if steady state is to be reached is l m = a2 > r1MaILe k

-8and in the steady state an'- 0 for n - 1o One may now ask- under what condition will a higher mode, say the second 4armonic, not go to zero? The answer is immediately clear from the above equations; the condition isoC22 Be g or 1 a' The answer makes the question academic because under this condition the fundamental mode is violently divergent, and the system will be destroyed~ One can ask a different questi.on however, which is not academic, namely~ what will happen if the reactor is poisoned in some manner, its size is increased in order to keep the fundamental mode just critical, and if in this process the condition where B2 refers to the clean reactor, becomes established? An answer to this question can be given by examining the criticality factor.* of the second harmonic after the addition of the poison. In the subsequent equations s ~c) and J,,(P) will refer to the eigenvalues of the clean and poisoned reactors respectively. For the clean, critical reactor one has ikr~n, -_ _ + h = and — ia < ld 4- tr(cr) i+ 1232 -- Ac ImoL a (c) * In this thesis the terminology and notation of Weinberg and Wigner (see Reference 1) are adopted. The term "critiicality factor"' and the symbol C are used in place of the more usual term "effective multiplication factor" and the symbol keff.

-9Now if the effect of the poison on the n'th mode is denoted by -AC, such that for the poisoned, critical reactor one has ~ h M2 se (P) then the criticality factor C2 of the second harmonic in the poisoned reactor is C, n,L 7 ( - - a/ dP) As long as (P)> B2 and is less than unity that is, the second harmonic is suband C2 is less than unity, that is; the second harmonic is subcritical, irrespective of / CJ If, however, 22 UP)4 13 -? c )z the subcriticality of the second harmonic is dependent on Z\ C2 t and variations in this quantity may cause divergent oscillations in the second harmonic. Similar arguments can be advanced for the higher modes It is apparent from the foregoing analysis that in order

-10to satisfy the last inequality mentioned above, the compensation for the addition of the poison must involve a relatively large change in reactor size. It is also apparent that the larger the clean, critical size of the reactor is in terms of its migration length M, the larger the change in its size will have to be to compensate for the poison. Therefore one would expect those reactors to be prone to spatial instabilities which are large with respect to their migration lengths.

CHAPTER II THE FORMULATION OF THE PROBLEM Having pointed out the possibility of spatial instabilities, the problem will now be formulated in terms of the equations of motion of the poison and of the neutron flux. The radioactive decay chain involving the production and decay of xenon-135 in a nuclear reactor core is shown by the diagram below.,3S FISSIO/N s/ Xe (s ) Te (Z.sF)-PI (6.7A) I (2e6,6.ro ) Xe (3. ) The numbers within the parentheses are the half lives of the nuclides, The following table summarizes the thermal fission yields of tellurium-135 and xenon-135 Fissioning nuclide U-233 U-235 Pu-2 39 Yield of Te-135 0 051 0.061 0.055 Yield of Xe-135 - 0.003 All the above mentioned data were taken from Reference 3. It is apparent from the diagram that the production and depletion of xenon-l35 in the decay chain are dominated by time constants of several hours. For this reason, the equations governing the behavior of the poison concentration can be written as if iodine-135 were produced directly from fission. Also. the branching that occurrs in the decay of iodine may be neglected since the xenon isomer decays to the ground state relatively fast. In addition to the radioactive decay, I-135 and Xe-135 are depleted - "burned out"' - by neutron -11

- 12absorption. The neutron absorption cross section of 1-135 is quite smaill however, so that its burnout rate is small compared to its radioactive decay rate, and therefore *the former will be neglected. The equations governing the poison concentration can now be written subject to the above mentioned approximations. They are~ T" X(r,~) =..I (r, —, ~'x,r- h(.t, ) -,,,(r,. ) _ A4 (', ) where X and I are the concentrations, AX and /\rthe radioactive decay constants, ~xand bthe fission yields of Xe-135 and 1-135 respectively F is the fission. rate, and A is -the absorption rate of neutrons by Xe-135. The last two quantities, of course, depend upon -the neutron flux. Now the formulation of the equation describing the rate of change of the neutron flux remai.ns~ Desregarding delayed neutrons and feedbacks (temperature, poison) for the time being, the thermal neutron flux in a bare homogeneous nuclear reactor can be described by the equation D VX/1A z o'B affl =A);7ly \ (! -,D i-2<,}(r>&)#SP2,cij+6' r t#tsE) ol t' (9) with the boundary conditi.cn-r that the neutron flux vanish at the extrapolated boundary (see References 1 and 2), The assumptions involved in this model are~

- 131o That Fick's rule provides an adequate representation of the neutron current, 2, That the extrapolated boundary has the same position for neutrons of all energy, and 3. That the d2/aZ term of the rigorous time-dependent diffusion equation may be neglected. In equation (9) r is the average speed of thermal neutrons, D is the diffusion constant Z5and Zf are the thermal absorption and fission cross section respectively, V is the average number of neutrons produced upon fission, is the fast fission factor, p is the resonance escape probability, and -''I-) is the slowing down kernel, that is, the probability per unit volume and time -that a neutron. produced at a point r' and time t' will become thermalized at r and t in the absence of fast fission and resonance captureo The slowing down kernel used here is normalized to unity, that is, its inter gral over all space and time in an infinite medium is unity. In order to take into account the feedback effect of Xe-135, an additional term is now added to the above equation. This term involves the space- and. time-dependent poison concentration X(rt) and the thermal neutron absorption cross section of the poison ac. In adding this term, the approximations are used that Xe-135 absorbs on.ly thermal neutrons and therefore does not alt er p, and also that vLu, Dc, and E are not influenced by the poison. Subject to these

-14approximations, equation. (9) now becomes C I (0) where is as in equation (9) the thermal absorption cross section of the where ~is, as in equation (9). the thermal absorption cross section of the clean core. Now if the changes in b are slow with respect to the slowing down time of the neutrons, this equation can be written as follows- ( / ) = D V ( - ) (11)'T —, pif No~) )'(it _// k One should at this point observe that the time constants governing the variations in Xe-135 and I-135 densities are very large compared to the mean neutron life time*,f and even compared to the longest delay time of the delayed neutrons. This means that the response of the neutron flux distribution to variations in the poison distribution will be very fast compared to the speed of the variations in the poison distribution, and if the short-term behavior of the neutron flux distribution is stable, the flux distribution will at all times be very nearly in equilibr.ium with the poison distribution. Therefore, the flux distribution can * In this thesis the term "mean neutron life timer refers to the average time elapsed between the appearance and disappearance (by means of absorption or leakage) of the neutrons, and the term Itmean generation time" will refer to the average time elapsed between the occurrance of a fission and the disappearance of the neutrons (promt and delayed) produced by that fission.

-15adequately be described by the following equation. (12) where ~ is a time constant of the order of the mean generation time of the neutrons, about 01ol second. In order to demonstrate in what manner the short-term kinetics enters the problem, the ~/} term of equation (12) will be carried along up to a certain point in the development. To summarize~ subject to the enumerated approximations, the following equations describe the interaction of the neutron flux and xenon poison distributions in the absence of other feedbacks. o T r;r j) y - (13) ~ 15

CHAPTER III THE STEADY STATE The model. described by the final equations of motion of the previous Chapter will now be examined in the steady state. This examination is not only instructive, but is also a necessary forerunner of the dynamic analysis. In the steady state the final equations of the previous section become D7+(t) -Z (-u) -dX(r ) #)) + - ae o (16) (17 A (D) 6 v E +r) - A4 Xr )- X()4(t =,1) From these equations the expressions for the steady state distributions of iodine-135 and xenon-135 can be derived quite simply as functions of the neutron flux distribution. They are~:Z7(Ptr~~~) =a ~ $ $~(Sr~)~, (19) ~X(t) = }rJ, 2(20) where = — ~z +,. It is apparent from these equations that the iodine concentration is everywhere proportional to the neutron flux. The xenon concentration is approximately proportional to the flux if A/x - 4 (- ) For Ax 4 r ~-x (r) the xenon concentration saturates and becomes - 16

-17practically independent of flux level. The assumption of the former condition is called the low flux approximation, while that of the latter is the high flux approximation. The break point between these two approximations occurs at 9= ) (r). For the values Ax = 2.1x10O- sec and - 3x108 cm x x the flux level at the break point is about 7x1012 neutrons per cm2 per sec. It is interesting to note that thermal power reactors usually operate with an average thermal flux of about this order of magnitude, so that neither the high nor the low flux approximation is valid for an appreciable portion of their cores. One must realize that even if the high flux approximation is valid near the center of the reactor, it will not be valid near the boundaries where the neutron flux goes to zero. Substituting expression (20) into equation (16)2 one obtains the equation for the neutron flux distribution in the steady state poisoned reactor. The equation is: Dv x2 + x q (ek ) (21) The effect of the poison on the criticality factors of the different modes of flux distribution, and also its effect on the steady state flux distribution will now be calculated by comparing equation (21) with the corresponding equation of an imaginary critical reactor which is identical in size and compositionr except that it contains no poison and that the A of its fuel has a different value:o )/ 2 The equation of' the

-18-:maginary reactor is~ +7 ( t) t ~ (k) 4)~ (22) with the same boundary conditions as those of equation (9)~ It is one of the fundamental theorems of reactor theory (see Reference 1) that in a reactor model described by equation (22) the spatial distributcion of the steady state neutron flux is the fundamental solution of the wave equation, V72S,j) + X4, ( 1) 0) with boundary condit;ions the same as above, so that ()=,o (-k where Ao is an arbitrary constant. It is also a fundamental theorem that if d(r) satisfies the wave equation then i)) ) (23) whereO(J,)i.s the Fourier transform of'(iw-e';) p" Q) 4c(Q)sn:,ir2 12 (24) The proof of these theorems i.s shown in the appendix.

-19 - Using the results of these theorems, equation (22) can be written. as follows: D V (r) 2 - pf Q (4)#r)= 0Q (25) or 5- D4 -2 r p2 QC()]r ir) (= o so that the familiar criticality condition is apparent: C' 1 = ) where and D Since is everywhere proportional to 1 the latter function also satisfies equation (22), and in general one can write D7 7 ttr - R, v 4 ) 94'") J a'- 1k t 1 (26) or since VZ (J) = -- L /:' () (27) and k e ta) t o/ (i = - e (h) (28) where the/'s are constants, )A4, being zero. One can now write down the

-20relationship between /p and oLl. The relationship is: /e = - DoJ0Z - Z, + V ep ~ Q Co~ )v (29) or / _ k _ _ (_ _ _ ) where k' (X)) is the criticality factor Cl of the i'th mode ~[(e)of the neutron flux distribution in the hypothetical reactor. The effect of the poison on the criticality factor of the i'th mode can now be derived by determining the change in V necessary to compensate for the addition of the poison. This change can be determined by comparing the equation of the iith mode in the poisoned reactor with that in the hypothetical reactor. The equations areo Ax + - qc(r) (30) +- 2G f zk i()t-' - / 4' Lk ) and op'~q,': (t) -'. x <,,')rl-)d'' (.). (31) Making use of relationships (27) and (28), multiplying the equations by and integrating them over the volume of the reactor, and then subtracting one from the other, one obtains _ _ _ _ _ _ _ _ _ _ -r. _ _ _ _ X (r)&a)d 3 (32) z_ vJ66e2)s ( c ZlP:Ldr, ~ ~rc cpV

-21The left side of this equ.atio;n is clearly the change in the criticality factor that had to be introduced by raising V to 11 in order to counteract the effect of the poisons Therefore the effect of the poison on the criticality factor is: ___1 ____ _ z ~ r) ( 33) Now if in the first approximation then ACT - 2, (qALf20 Z4) i 2 dt Z~ (34) and in the first approximation the change in the criticality factor of the reactor is *LCAC, t(6-,) XOne must be very careful to distinguish between the criticality factor of the reactor and that of the fundamental mode ( of the wave equation 7 t+k +4, =-0, In the case of the hypothetical reactor, where these two quantities are of course the same. In the case of the poisoned reactor, however, the former terrm refers to the criticality factor of the

-22flux distribution + that is, to that of the persisting flux distribution, and the two quantities mentioned will be equal only in the first approxirlation. The accuracy of the above mentioned approximations can be investigated by computing the approximate difference between the flux distributions (r) and i(r). The computation of this difference is based on equations (21) and (22). Introducing the expansion into equation (21), making use of the realtionships (27) and (28), multiplying the equations by IL, integrating them over the volume of the reactor, and then subtracting one from the other, one obtains +-' >') z?Ae (36) According to equation (32) and. therefore equation (36) may be written in the following form~ Ai JVI 3x Ct1e r=D(37j Ai /O

-23Now if in. the first approximation then the series in equation (37) will be dominated by the n= _ term, and therefore in the first approximation A A~~~~~x 13 3 __ ____ -~ Kr (38) V and since in the steady state the flux and poison distributions are even functions about the center of the reactor, A =O for i even*. According to equation (29) because/~ =O, and therefore /wM Ole -,I L Ji e C.2e 31 Q~deL) (dP 7 Ce~ (3.9) so that finally 9 c X 3,4d A1 ZG~JI#~L2QIJ) 0&,LeL) Ic112~La0L _ ~141AJ-C: 7 (40)' (*&,) J - Q (U~) In this thesis the fundamental mode of the wave equation is called the first harmonic, so that the even numbered harmonics are odd functions about the center of the reactor.

-24One should note that according to equation (33) the quantity is the maximum effect that the poison can have on the criticality factor of the itth mode. It should also be noted that this analysis leaves An undefined. One must at this point remind himself that the results calculated on the basis of the approximations used can be trusted without further analysis only if the predicted changes in the criticality factors and the neutron flux distribution are relatively small. A simplification is possible if the above results are applied to large reactors. From equation (24) it is apparent that Q I() can be written as the following infinite series: E7 (_Wkh Q/ R ~ 1 - < Xe -' + gag (41) ZL (n2(41) where e2 z ~-'i I and 4- 4i[ 2 /)2aV Pso that P. is the mean square distance traveled by a neutron during moderation. Now if o. is small, Q (, ) can be adequately expressed by the first two terms of the series: _ 4 ( P )i a (42)

_25where -A= 16 M Using this approximation, the criticality condition and the expression for the perturbed neutron flux distribution become C____ t 2 _ be 2s 0, —_. (_43). -"~1.4 7ao0L zv(/Ax + where l/ -/= 2 and AI'1 i = Z a(d-#-LW7,) fi 2(. M 14') J _ _f ( z From the last equation it is apparent that the magnitude of the distortion in the shape of the neutron flux distribution caused by the poison is governed by the maximum effect of the poison of the criticality factors, the size of *the reactor, and the values of the integrals The values of these intergrals will be examined in detail in a subsequent Chapter. At this point it is sufficient to note that the value of 7-K

-26and therefore the values of the integrals, will be small when K //)x ~. and that when >a i/ >x 2> o~x I 4 4L and therefore the values of the integrals will again be small since / Ar a Oe L. The absolute values of the integrals reach maxima for (maximum) of the order of 10 neutrons per cm per sec. A measure of the distortion in the neutron flux distribution can now be obtained by computing A3/A 1 For a slab reactor with boundaries at x=o and x=H, k/p=l2. I(maximum)=1013 neutrons per cm2 per sec. the maximum effect of the poison on the criticality factor of any mode may not exceed furthermore 5 -4 X2, 2 itj dZff+ / -*i~ a 6 4-12 -.06) Ax 4

-27and the calculations yield the results exhibited in Table 1o TABLE 1. HD 25 50 75 100 A3 016.059.13.23 Al Maximum steady state flux flattening produced by xenon-135 as a function of reactor size. From this it is apparent that the reactor has to be very large compared to its migration length in order that the flux shape distortion caused by Xe 135 be more than a few percent. Since A3 is positive, the flux distortion is such that q(r) is a flatter distribution than (r), and the ratio of the maximum to average flux in the poisoned reactor is less than that in the clean reactor.

CHAPTER IV TIE EXPANSION OF THE EQUATIONS OF MOTION In order to investigate the response of the flux shape described by the equations of motion (13), (14), and (15) to disturbances about the steady state, the quantities +, X, and I will each be separated into two components- one standing for the steady state value, and another denoting the deviation from the steady state. Accordingly: T r,) c Substituting these into the equations of motion and subtracting the steady state equations from the resulting ones, one obtains the following seto ~~ y t cA- 4- Lxo6XCaXf]~i?} e3_ (45) wt SX=Ii1 slcti~eZf +-/)k X- F~j~o i*+@o itxSt I (46) One now has a set of nonlinear, partial differential equations describing the motion of the deviations from the steady state~ In order to reduce these to ordinary nonlinear di.fferential, equations, the - 28

-29deviations will be expanded in the following mannero oat where, as in the previous sections, and Substituting these into the equations, making use of relationships (27) and (28), multiplying the resulting equations by Jt' and integrating them over the volume of the reactor, one obtains the following result~ cZ- -' o;6 d I~, jt ~d /- +* Q 2'2 50) -l Cx" + VC~O ev+,,r h r _7 O, - ~ir~o) LX I"1~i/b~i~il~drea~C~;3II +a50

-30One now has, in place of the three nonlinear, partial differential equations, an infinite set of ordinary nonlinear differential equations. There is a set Of-three equations for each mode of the flux distribution, but the sets are interdependent so that the time behavior of each mode is influenced by the behaviors of all the other modes. The influence of the i'th mode on the m'th one is represented' by the expression within the second pair of brackets in equation (48). The same expression appears in equation (49). The integrals these expressions act as coupling coefficients. Since the /{/5 are orthonormal functions, one can expect that the absolute values of the integrals within the second pair of brackets are less than those of the corresponding ones within the first pair of brackets in equation (48). Furthermore, one can expect that the greater the difference between i and m, the smaller the absolute values of the integrals become. Now then if the deviation of the neutron flux distribution from the steady state flux shape is dominated by one harmonic, say the n'th, so that the total effect of the other harmonics on it is negligible, then the time.behavior of -this mode can be described adequately by the following set of three equations. (51) fSE [X 6f -<JXO a3 ig CL~ - [o"3 3 ~ a(52), /lt Dr YI 6 W g ad -z I

-31The validity of these equations depends on the degree of the modal interaction, that is, on the size of the modal coupling coefficientsand on the degree of predominance of the n'th mode. The last terms of equations (51) and (52) are nonlinear. If the deviations from the steady state are smalli these terms are negligible with respect to the other terms, and the system behaves approximately as a linear one. It must be pointed out that not all the nonlinearities involved in this problem are represented by the last terms of equations (51) and (52). The requirement thatq( / ) X (_p /), and I (etj) be non-negative, places a limit on the amplitudes of the higher harmonics. The derived equations of motion of the system will therefore hold only if the ats, b's, and c's are within certain limits.

CHAPTER V THE DIMENSIONLESS FORM OF THE EQUATIONS OF MOTION Before beginning the examination of the equations of motion, it is convenient to cast them into a dimensionless formo From the definitions of amn, bm, and Cm, the amplitudes of the m'th modes of the deviations of the neutron flux, xenon,and iodine distributions from the steady state are dw,,~ (max), -, (max), and C 4 (max) respectively, where, (max) is the maximum value, that is, the amplitude of ~ (~). Since.(),4)X(Kj)j and I(_,)must be non-negative, |ct*t76K (wax) X0 (V~ 0J-X x 4- (max) k&4ot! (/d/h~ ($4x) 4 X0 (ThX =x) / + mx ) if the deviations from the steady state are purely in the m'th mode. If one now defines 2 4 - (L d X) I (f) - r I n/4I ( Xj ( x- 1k1 (s x) -32

-33then under the same condition Ixi4 4.IY) are dimensionless and could be called reduced amplitudes. In terms of the new variables, the equations of motion can now be written in the following form o ~C>1 + 4L2)L/ R 2Ii?,, ) + — / 4/,..,a ix, -V _- /Y5 __ - [ J4-, fX (a:) 2~ - bL IV. t j r X'j4c -d4 4 4/,f L 0 ) (5 6) oihiH) /: —:. IV~ -Z,. j,... 0- d 3,ON /I _Ex c. (:d)<) 3

-34The terms of these equations are dimensionless. The coefficient of, in equation (54) should be recognized as _g~, =,,- - i the excess criticality factor of the m'th mode, and the expression multiplying the rest of the terms on the right hand side of the same equation should be recognized as the maximum effect of the poison on the criticality factor of the mt.h mode. In order to reduce the equations of motion to a tractable form, a number of new quantities will now be defined0 the effective generation time of the m'th mode, the dimensionless time variable, r 6_ the maximum effect of the poison on the criticality factor of the mtth mode, *the reduced subcriticality of the mith mode, Q-fL- - X(n)

-35 - the reduced maximum steady state flux; j AR (MYaX) fV k Hi i and so that 2' P., is the effect of the poison on the criticality factor of the mgth mode. In terms of the newly defined quantities, the equations of motion can be written in the following form. A Sm Alk]~, *- mXim L+ ew1\/XN- (57) =IT /t1, 1,2,..... 4... 5 ol T VA/

CHAPTER VI THE MODAL COUPLING COEFFICIENTS The behavior of the modal coupline coefficients FijYij04 and Pij will now be examined. According to the definitions in the preceeding Chapter, these syrmbols stand for the following expressions. (60) ~jh~ k (M~x) tv di,(61) (nc (62) The examination of these expressions leads to the following observations~ 1. Fij and ijh are linear functions of the reduced maxi.mum flux L o 2o P, is a nonlinear function of the same quantity0. 1ja 3o For T40/( -ax)9 = //!/ (max), that is ~o=i Fij= j ij 4. As JL- * o R FQj 5. As &L -~,?J -. for =, and Pj 0 for z' J since k and /, are orthonormal functions. -36

- 376. As qo approaches a flat distribution FF> f Fj )j Fe j W =j * 0 F o# J F:; _ for _j > p6j w a Xor ci j. that is3 the flattening of the steady state flux distribution tends to decrease the interaction among the modes. In order to illustrate the behavior of the modal coupling coefficients and to provide data for subsequent numerical calculations3 the val.ues of some of the coefficients were calculated for a slab reactor with boundaries at x=O and x=H, using the approximation (max): ( r) $ll (max), and the requirement that no variations be allowed in. tke y and z directions in any of the variableso As it was shown in the Chapter on stead state, So is a flatter function than, and therefore the calculation based on the above mentioned approximation will yield the maximum modal interaction. For this reactor -, = SLZ ir X so that; (tW.X)= C S and _Stn Introducing the following change in variable~ 77-"'X

-38the expressions for the modal coupling coefficients become j~ =l t t S~h & 5 a {( S) 5tt (ji) 2&2(635 gir = iL 2 JD Sin tK6B) six/ (ji) stz (AS^JOS (64) ~' 2LJ = L2 ns SL'& 6)sii ((J) L (65) Since, about the center of the reactor, i is an even function for i odd, and an odd function for i even,;.= _pj O if i j is odd, and =j h 0 if:J4 is eveno Therefore, in the linear range of oscillations, where the nonlinear terms are negligible, there is no interaction between. an even and an odd harmonic, Also, if modal interaction is neglected, the equations of motion of an even numbered harmonic contain no nonlinear terms, since /- = O for m even, and therefore the only nonlinearities influencing the behavior of such a mode will be the limits placed on the magnitudes of Nml Xm. and Imo It should also be noted that the following relationships hold, and P P - -P, where m=e, +J)' y=-(-,') ~ro _'.t od,= t fI_ (j- o,

- 39Since j is directly proportional to 4I g / -is a costanto The general. expression for this quantity is shown below. ji =- T ( j)_2[ (_J)_' -) (66)'. j + h not even. Some of the values of this expression are shown in. Tables 2, 33 and 4, As it was pointed out previously, j =:,,j, and therefore the tabulation of the latter function includes the values of the former. The behavior of Pi j a nonlinear function of -1-" is shown in Figure lo

-40TABLE 2 Values of 3T- 7T 37 e&LrL, 8zL J 1 2 3 4 5 6 7 j>> [ 1 _ __ 1 0 1 3 5 35T5 35 1 27 5 7 9 35 0 21 5 3 4 8 0 16 0 8 0 12 15 2 - 1 33 3 - 3 d o - -- - - - 44.. 1 8 4 612~ TABLE 3 Values of 1 2 3 4 5 6 7 j >> L' 4 4 4 4 12 2 3 0 0 1 0 4 8 208 18 3 0 7 0 15 0j o8 0 8 0 40 0 56 24 35 ~ 15 77 195 j

- 41TABLE 4 Values of 3_ 3 8JL J3,,J] 1 2 3 4 5 6 7 j>>d+3 1 -~~27 1 7 9 _5 35 O ~ o 0 4 8 1o8 18 7 15 o35 Do 27 1 27 0 27 27 4n 15 0 8)1 55 16|216 -3 9 - o 5 o 55 5 o L 31 4 c 84 2- 3Z,42

(xow)p o~x 0~:y 73 (ww)O X.0 001 08 09 0 Oz 01 O 9 tp I a' 9' t'' I' 10' OD 4~U ~ I 0 0 4-,,+U'Ud-UI!~~~~~~~~~~~~~~~~~~~~t d ra 0)'-I _.r~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~-I /,.,, //~go' Pg0 C: /~~~~~~~~~~~~~~~~~~~P 0 O" 001- ~~~~~~~~90' 60',d o 0I' // / Z+U'Ud-W!'-X,,i': " -4 0~' 0 H 09',/ // 0; / o'.,////~~~~~~~~~~~~~~~~~~~~~~~~e'./ 0~: O// o m. 00, o~ uw,-D /, 0 / r, d / 09,,,4, - ~'~O " -~ %, u 08''~~1~~~~~~~~~ II = oz._,.,; — 77~~ —-7-~~~~~~~ OD- —;, - -- V —7-/Tt*+Ulu -U4 tur E.0 ODN& ~ +u'u d- Wl,', 00%0

CHAPTER VII THE DERIVATION OF THE STABILITY CRITERION WHEN THE EFFECTS OF THE NONLINEARITIES AND THOSE OF MODAL INTERACTION ARE NEGLECTED The investigation of the stability of the system described by the equations of motion (57), (58). and (59), will now be begun. As it has been pointed out in Chapter IV, if the effects of modal interaction and of the nonlinearities are disregarded, the behavior of each mode can be described by a set of three linear differential equations. In this section the stability criterion for such a simple system will be derived. The effects of modal interaction and those of the nonlinearities will be investigated in the subsequent Chapters. The system of equations describing the behavior of the n'th modes of the deviations of the neutron flux, xenon, and iodine distributions from their steady state values, in the absence of modal interaction and nonlinear effects, are shown below. id / -AI7// - Fs v X ] (67) d14 I = 4- I t g -&Ix /V F 7 (68) -z _- (69) -43

-44Eliminating Xnand Inn one obtains the following third order differential equation. ~_r A L 12 Lt fli~hI z l 4kda where /`- / I /Ix, An examination of the coefficients of this equation reveals that all of themi except that of the first derivative term, are at all times positive for S' 0.~ The coefficient of the first derivative term, under the same condition, is also positive if A/! > Po o Furthermore, all of the coefficients, except that of the third derivative term, involve /I S~, the suberiticality of the n'tth mode~ The order of magnitude of the coefficient of the third derivative 10-6 term is 10-6 For = -: n S O. 09 but for ~t> it will be of the order of 10-2, unless the reactor is of extremely great size. this will be shown later in the hapter in the chapter Therefore as expected, for the higher modes the coefficient of the third derivative term is very much smaller those of the others, Since th those of the others Since this hess deals with the stability of the shape of the flux distribution, the interest is focused on the higher modes, The assumption that the short-term behavior of the reactor is stable implies that the fundamnental mode is

-45stabilized by some process which has negligible effect on the higher modes. The stability of this system, for kZ, will now be examined by the application of the Hurwitz-Routh stability criterion. According to this criterion, the differential equation e3 ~e t ex oL3+ d~x + e dt + 2ox =o represents a stable system if and only if e3>o, e >o, o eo?o, and eae - 3 e3 ~o. To satisfy the last condition, one must have r3r [- /) nSH 2F ) — x Fn, / - z) — A / / K) XI1 OKg^F" (71) If one now disregards th. ~-last two terms of this inequality, which involve U*Ax 9 this condition becomes (/hriz)Sn F,, 4-v ( P,-,t) % ~> (72)

-46which is precisely the condition for the stability of the equation +0,,~'- a -, + L, ( _,,, - o' + (F,) F (P)]n - (73) which in turn results from equation (70) when. the terms involving 7n 2\x are neglected. It should be observed that, unless Fnn is very large, the presence of the short-time constant An detracts from stability. The physical reason for this is that taking into account the fact that a change in the neutron flux level does not instantaneously follow a change in xenon concentration introduces additional time delay into the problem. This means that there is an increase in the time required for a decrease in xenon concentration to produce an increase in flux level and in turn an increase in iodine concentration and in xenon production rate, which will reverse the original decrease in xenon concentration. The approximate stability criterion for the n'th mode can now be cast in the following form. S1 > (74) The behavior of the right hand side of this inequality is shown in Figure 2. Since S >, this criterion is always satisfied when P' n < YVf/6. The physical meaning of this becomes clear

1.0.vrl l l l l l I rI~~~~~~~~~~~~~~~~~~~~~~~~i I I 1 111 l l l l l~~~~~~~~iI' i iI i ii i,f'/)F:.95 ______ _......8~~~~~~~~~~~~9 n:3.6 - - -.4 ----.2 5 1 2 5 10 20 50 100 Q Fnn (Pnn YX/- ) Figure 2. The Behavior of/x+ F as a Function of Q.

if one recalls that P, is associated with the burnout rate, and ~x/d with the direct production rate of the poison by the n'th mode of the flux distribution. If RPg x/, then an increase in the n'th mode of the neutron flux distribution will result in a prompt net increase in the production rate of the poison, which in turn will tend to reverse the increase in the flux. If, however,., > /', then. an increase in the flux will result in a prompt decrease in the production rate of the poison, which in turn will tend to increase the flux level further, and possibly result in instability. The steady state neutron flux level above which the burnout rate of xenon exceeds its direct production rate is about 3x101 neutrons per cm2 per seco The right hand side of inequality (74) is an increasing function of &L, and approaches dI/ as L — oc. Therefore, according to the approximate stability criterion, if Sn,~~> I 4~~/,~~ _2,~ ~~(75) that is, if the subcriticality of the n'th mode h go is greater than the maximum effect of the indirectly produced poison X [/~ on the criticality factor of the same mode, then the n'th mode is stable. Or, in other words, if the ntth mode is subcritical irrespective of the effect of the pqison, then it is stable. One should now recall the very similar result reached in the preliminary investigation of an elementary reactor model in Chapter I.

-49The stability criterion is often exhibited in a form involving the difference between the squares of the characteristic values o n and OY, instead of Sno This form will now be derived in order to show the approximations involved in it. Recalling that, if the reactor is large with respect to its migration length, and if n is not very large, then 25 1 - r(76) and that for a large critical reactor from which ( 2 + - / ) then F S ( 0 gg- l A m(2 _ (77) and if one now ignores the first term on the right hand siade, 5 ^ _ A 2( 2 _ ) (y8) One can now perform some rough calculations and get the approximate reactor size at which flux shape instabilities become a possibility. For a slab reactor with r/ 7- fi t. ~, ILI1)

-50 - For the first few higher harmonics, Table 5 gives the values of H/M for which S lo TABLE 5 no 2 3 4 5 6 H/M: 32 52 71 90 109 Values of H/M for which 5g,., According to Table 5, in a reactor whose dimensions are larger than 32 migration lengths, unstable second mode oscillations will result if the steady state flux level is higher than some critical value. It is also interesting to calculate the values of Sn for a reactor with dimensions as large as 100 migration lengths. Table 6 summarizes the results of such calculations based. on equation (79)~ TABLE 6., n 2 3 4 5 6 S, eoo0 7 ~ 51 1 1 2 Values of Sn for..M x 1,00 Now since,7 is of the order of percents, the subcriticality of the second mode even in this large reactor is greater than 10-3, and therefore one is justified in neglecting the terms involving the generation time if H/M < 100o

-51 Recalling the definition of rP 2 Y(14 )2 2) one should note that using a one energy group analysis, this quantity will be underestimated by writing it as Tin2o (1~~/d12L6n) Consequently, Sn and the stability of the reactor will be overestimated by one energy group approach. The calculations in this Chapter show that when modal interaction and the terms involving the neutron generation time are rinegleced, the time behavior of each mode can be described by a second order differential equation. The examination of this equation in the linear range yields a stability criterion which states that in order that the nTth mode be stable, its subcriticality must be larger than an increasing function of the maximum neutron flux level. Since the value of this function is always less than the maximum effect of the indirectly produced xenon-135 on the criticality factor of the ntth mode, if the subcriticality factor of the ntth mode exceeds this value, the n'th mode will be stable irrespective of the neutron flux level in the reactor. It was shown that the subcriticality of the higher modes decreases as the reactor size increases, and therefore as the reactor is made bigger and bigger, first the second harmonic becomes unstable, then

-52the third, and so on, It is of course still possible to have a reactor in which the second harmonic is stable due to the action of a suitable control system, and the third harmonic is unstable. It was shown that the terms involving the neutron generation time can be neglected when the reactor size is not extremely large; their effect, however, detracts from stability and will become important for extremely large reactors. It was also pointed out that a one energy group analysis will overestimate the stability of the shape of the neutron flux distribution

CHAPTER VIIIL THE EFFECT OF THE NONLI'NEARITIES In the previous Chapter the stability of the n'th mode was studied with the assumptions that the oscillations are small about the point where Nn = Xn = In 0, and therefore the effect of the nonlinearities can be neglected, and that the oscillations are predominatly in the n'th mode, and therefore modal interaction is negligible. In this section the effect of the nonlinearities will be studied while retaining the second assumption mentioned above, The terms involving the relatively small time constant An will hereafter be neglected. The system of equations to be examined, therefore, is ~= 5n //4 - FnhX n / -n n A//A X< (80) of XI2~~~zT =F~ 7x,-Av~ r~ T(82) with appropriate limits placed on the variables, The method of analysis will be a topological one, based on the study of the solution paths of the system in the phase plane, mapped by point singularities and singular lines, or limit cycles, dividing the phase plane into topological domains in which the behavior of the solution paths of the system can be investigated by relatively simple methods. As a general reference on this subject, see Reference 4. -53

A few definitions will now be made before beginning the analysis, The phase plane will be defined by axes Nn and I, so that the solution paths, or phase trajectories, will consist of plots of Nn against In as time progresses, A limit cycle is a closed solution path corresponding to periodic motion. Point singularities correspont to points of equilibrium. of the system. On the basis of topological considerations one can state that any solution pathin the phase plane approaches either a singularity or infinity for both f --- + - and + -— o. A stable singular point is one that is not approached by solution paths as - - —,- all other singular points will be called unstable. A stable limit cycle is one which is approached by solution paths from both sides as + —o +; for an unstable limit cycle the same holds true as i -o — A half-stable limit cycle is one which is approached by solution paths from. one side as,-4-,og and from the other side as i — oo The phase plane will how be explored for the presence of singular pointso This will be followed by the search for the location and stability of limit cycles. Since singular points correspond to points of equilibrium of the system, their locations will be given by the solutions of equations (80), (81), and (82) when the derivatives are set equal to zero, if for the time being the limits placed on the variables are ignored. One should note that at any singular point Nn = I n If n is even, the nonlinear terms vanish, and the only solution is at Nn = Xn = In.= 0.

The stability of precisely this singular point was studied in the previous Chapter. Since, if modal interaction is negligible, the equations of motion are linear for n even, linear stability is a necessary and sufficient condition for stability for the even numbered modes. The effect of the limits placed on the variables produces a limit cycle about the singular point at the origin of the phase plane so that if the system is unstable, the solution paths approach the limit cycle, that is, the oscillations build up to some periodic motion. This is illustrated in Figures 3 and 4. Using the particular values of the coupling coefficients at -L 2 calculated previously, the equations of motion, together with the limiting conditions, were solved with an analog computer. Marginal stability occurs for S2 =.186, Firgure 3 shows the solution paths of the system for S2 = 190 (stable), and Figure 4 for S2 =.180 (unstable). The former shows that the oscillations damp out, and the latter that the paths approach a limit cycle from both sides, resulting in periodic motion. Figure 5 shows the time behavior of the variables for this periodic motionr In this case the size of the oscillations is determined by the limits placed on N2, the other variables do not reach their limits. Limits of +.70 were used for N2 in this case, As S2 is decreased, that is, as the system is made more and more unstable, the limit cycle moves outward until the limiting on-the other variables comes: intuo pl.ay. For n odd the nonlinear terms do not vanisho It is interesting.-to compare:; the ma-gnitmd...e of the nonlinear term with that of the

-56r2 N2 0.8 START 0.8 START /S ATART 0.4 2.4 Figure 3. Solution Paths of the Figure 4. Solution Paths of the Second Harmonic for Second Harmonic for S2 =.190. S2 =.180.

-o09T = oS ao$ $EfUOUU0H puoa' a9Z 9 o uo= o J~poW Psa anT suTnoH ib, q z sS ys / Z H "357IVS 31Wll 2~~~~~~~~~~~~~~~~~~~ ~~~~~~~~~~~~~~~ - -01 - - -- - -r - --- ------ X i Xt tX X ~~~I IF - -- - - - - +llllltllflIXlTIF Ti T ~~~~~~~~~~~~~~~~t Iliiillll~lli iii- ---------- E E lIl t E i!11 Ili.llT-I I I 11 [II Illlll lLllmll ~jj~]~_: ~ _ll.llll ff fg1 fgfg fiNfr I Ff:i I —'_T_4-fI_ 1_ 11_l111II I II1IIII — -0 r lI TI DNr 1LII III _sI II 1 II~ II ILl L lli IIX s|X ]Z^ III I —-IT -- II -- -- --. It, llTT- 11 llllli-llllIIIlG

-58linear term preceeding it. The ratio of the two terms is S nn Xn_ I;A n xn Using the particular expression computed in the section on the modal coupling coefficients. — ann (4vl2 <for n odd. Fnn 3 n3 For n~l >31 - and ln~nr -O as n- o~. Furthermore,!Nnl 1, and therefore the effect of the nonlinear term is not expected to be very large even for n = 3 and will be negligible for large n.,. For n odd. the equations of motion predict, in addition to the singular point at the origin, a sin.gular point at x2=X~ — - ( (83) ItH S +j (84) One should observe that both of the terms within the last pair of brackets are positive, and therefore, since a/, >l /Vn DBut becaulse of the limits placed on the variables, Nn must. remain less than unity at all times, so that the system can never reach this point of

59equilibrium predicted by the equations of motion. One should also note that at this point all the variables are negative. In the case of small oscillations about the singular point at the origin, Nn and Xn are of opposite sign, as one can see from equation (80) when it is solved for Nn (85) /F'X Xt= In order to reach the new point of equilibrium, Nn would have to pass the discontinuity at Sn = - Jnm X It has been pointed out previously that nonlinear terms may contribute to or detract from the stability of a system. The effect of the nonlinear terms in this problem is demonstrated in Figure 6. This figure shows the limit cycles reached by the unstable third harmonic with and without the nonlinear terms. Here again the values of the coupling coefficients atL_=2, were used; the limits placed on N3 were +1 and -~5. The inner limit cycle was reached when the nonlinear terms were left out, that is when;3,33 was taken to be 0. When the nonlinear terms were included, the solution paths converged to the outer limit cycle, as is shown in Figure 7, indicating that the nonlinear terms promote instability. On the basis of this:fact one must conclude that if linear analysis predicts marginal stability, for n odd, the system is actually unstable. If Sn is slightly larger than the critical value predicted by the approximate stability criterion, then there must be a region about the stable singular point at the origin of the phase plane in which the solution paths approach the singular point, If however,

-60N3 N3 START 1.0 1.I0, - LIMIT CYCLE, if 0~~'.7333,//C LIMIT \ CYCLE.0 /'~~~'"; i LIMIT CYCLE / 7 O~~~~~~~~~~~~~~~~ 0.5 / 05 / $ ~~~~'3a~~~~~~~~;' O 5 li /,/ Figure 6. Limit Cycles of the Third Figure 7 Solution Paths of Harmonic With and Without the Third Harmonic the Nonlinear Terms for f r S 170. S3 =.170. 3~ ~ ~ ~~I i

the system is displaced. sufficiently so that its solution path leaves this region, then the nonlinaar terms will be'sufficiently great to make the system unstable and the solution paths will approach the stable liTmit cycle established by the limits placed on the variables. One therefore has in this case a stable singular point surrounded by an unstable limit cycle which in turn is surrounded by a stable limit cycle. One can expect that as Sn,increases, the region of stability about the origin becomes larger, that is, the unstable limit cycle moves outward, and that if S is sufficiently large, the unstable n limit cycle reaches the stable one, and the system becomes absolutely stable. The above described expectations were verified by observing the behavior of the soluti.on paths of the third harmon.ic for a number of values of JL and S3. The results at _ — =, which are representative. are shown in Figures 7, 8, 9, 10 and 11. The critical value of S3 predicted by the approximate stability criterion is.175. The stable limit cycle reached, for S3 =.170, and the t-ime behavior of the variables corresponding to the periodic motion. on t;his limit cycle are shown in Figures 7 and 8 respectively. One should note that in this case N3 is the only variable that reaches its limits. Figures 9 and 1.0 show the locations of the stable and unstable limit cycles for S3 =1760 and.1765 respectively. It can be seen that as S3 increases the unstable limit cycle moves outward and the stable one inward. Figure 11 shows the solution paths of the system for S3 =.1770o Here the two limit cycles are superimposed, resulting in an inwardly

-62!~~~~~~~~~ 0? l l W r S, I II ~11 W - - - - - - - - - 3 --- - --- LLL i ~ ~ ~ ~ ---- --- I1 Z ~ ~ ~ ~ ~ ----- I 1~1 ~~~~~~~Tl TIESAE2AFe.2. OR FiueZ eidcMto fteTidHroi

-63N3 N3 I.OT I.O.61.6 I ~~~~~~~~.57. 5 Xr x /STABLE LIMIT CYCLE UNSTABLE LIMIT CYCLE 0.5 0.5 Figure 9. Limit Cycles of the Figure 10. Limit Cycles of Third Harmonic for the Third Harmonic S3 =.1760. for S3 = 1765.

-64N3 1.0 START 0.5 Figure 11. Solution Paths of the Third Harmonic for S3 =.1770

unstable half-stable limit cycle. The solution path approaches a periodic motion, and then spirals away from it toward the stable singular point at the origin.. For S3 >.1770 the system is stable. One should note that only a small change in S3, from.175 to.1779 is necessary to change linear stability to absolute stability. The calculations in this Chapter show that if modal interaction is neglected, linear stability is a necessary and sulfficient condition for the stability of the even nlmnbered modes, but is only a necessary condition for that of the odd numbered modes. This means that if the approximate stability criterion derived in the previous Chapter predicts that an odd numbered mode is marginally stable, it is in fact un.stable. The numerical calcuiations show that only a small increase in the subcriticality of a particular mode is required in order to go from linear stability to absolute stability.

C HAPTER IX THE EFFECT OF MODAL INTERACTION The approximate stability criterion derived in a previous Chapter predicts the stability of each mode by itself, that is, for the case when. oscillations occur in one mode only. The effect of the interaction among the modes on stability will now be investigated. It is reasonable to expect that if there are several modes which are predicted to be unstable by the approximate stability criterion, then the system consisting of these modes coupled together will also be unstable. This will be shown in the course of the analysis, The basic question to be answered, therefore, is as follows: if there is one mode, the n th, which is predicted to be marginally stable when modal interaction is disregarded, and all other modes are predicted to be stable under the same condition, then, if modal interaction is taken into account, will the system be stable or unstable? When the n'th mode begins to oscillate, it forces all the other modes to oscillate. The effect of the n'th mode on the mtth decreases as In-mi increases since the magnitudes of the coupling coefficients fall off rather rapidly. Now then, if the infinite series in equation (57), which takes into account modal interaction, converges, then one can estimate the -total effect of modal interaction oh the ntth mode by investigating the effect of the adjacent modes on it. This infinite series resulted from the expansion of a finite term, and therefore it must be convergent. -66

- 67Since it was found in the previous section that the effect of the nonlinear terms is small, modal interaction will be studied in the linear range only. The general equation describing the behavior of two coupled modes, the ntth and the m'th, was derived by combining the equations of motion of the two modes with the assumption that Fij = Pij = 0 except when i = n or i = m and j = n or j = mo This equation is shown on the following page. If the n'th mode is predicted to be marginally stable by the approximate stability criterion, the underlined terms in equation (86) vanish. The subsequent analysis will be carried out under this condition. The stability of the system represented by equation (86) can now be studied by applying the Hurwitz-Routh stability criterion. According to this criterion2 the differential equation d4x ~n-~e d~ dX e4 029+ e 1 ee X represents a stable system if and only if 2 e3 6 e* t v 6 O and. A3 - - e a- eo e3 >Ge The examination of the stability of the individual modes yielded a very simple stability criterion for5Llarge, namely Sn > XI/. The stability

-68 - - / {M Sn S - /nm } d': )'(I?4 I ~7 7~ ~ J L- )z" kLP 4- F/ IV tH) 5z Ttr X -< r() Ao Pn h X 1 g * SS, j + g?q 0 & Ps HA_ (e *? + 15? & _ & )4 L7 ( Ar?~ t? 54+i P )' d~~' (L~op 1? 1-P d

-69of equation (86) will now be examined for _large. As fL P, —. 2- ) t,,(l_-Pn — ) 0-P and if Sh - / that is, the n'th mode is predicted to be marginally stable by the approximate stability criterion, equation (86) reduces to the following form. WI ) C/ )/n 4 (87) Sn 5 T 4 L o_ T3 + AF4 )S^(S- Ty)0 ~ E(F __ F /5 _ 1/T3 ~O) so that and

-70Recalling that / />I/4 and /,I>/| /|, the condition for the stability of this equation is Sm > /, which is precisely the criterion for the stability of the m'th mode by itself. Therefore, for _L_ large, if the m'th mode is by itself stable it contributes to, and if it is by itself unstable it detracts from the stability of the combined system of the ntth and mtth modeso It should be observed that as S becomes large so that 5,m- ~r/~,, it can be factored out of equation (87). This means that as C-A, the effect of the m th mode on the stability approaches a nonzero limit. An examination of the coefficients of equation (86) reveals that this is true not only for ~L large, but for the whole range of ~ o The implication of this fact is very significanto It states that even if Sm is very large, that is, even if no variations are allowed in the m1th mode of the neutron flux distribution, there will be some interaction between the ntth and the mlth modes. The reason for this is that the large Sm does not keep the n7th mode from inducing oscillations in the m'th modes of the xenon and iodine distributions. Therefore, even if the magnitude of, say, the fundamental mode of the flux distribution is held constant by a suitable control system, in examining the stability of, say, the third mode, interaction between it and the fundamental must be taken into account:. The effect of modal interaction for the whole range of Q-L is not as simple as that for ~L large. A rather tedious examination of the coefficients of equation (86) in the light of the HurwitzRou.th stability criterion reveals that if n d Fn PnZ //v- I

- 71interaction with the m'th. mode always detracts from the stability of the n'th modeo For 5, >/h/~, the interaction promotes stability if the mrth mode is stable by itselfO Since an. increase in S promotes stability, it can'be expected that 5,= - will give the least effect if modal interaction detracts from stability, and that it will result in the largest possible effect if the interaction promotes stability. The critical values of Sng as predicted by the approximate stability criterion, together with, the function' FP / Fml are shown in Fl.gu~re 12 as functions of _2o In order to obtain an indication of the magnitude o.f the effect of modal interaction, and in order to extend the analysis'to the interaction among three modes, the equations of motion. of the second, fo'urth., and sixth modes were solved simultaneously with an analog computer using the particular values of the modal coupling coef:ficients calculated in a previous section~ Solutions were obtained for a number of valules of _L, S4, and S6 and 6, s'the difference between the actual critical value of S2 and the one predicted by the approximat;e stability criterion was noted in each case. The results are summarized. in Figure'3~ Curve 2 on Figure 1.3 shows the effect of the fourth. mode only. The difference between curves 1 and 2, therefore, shosit the effect of the sixth. mode when the subcriticalities of all the three modes were taken to be equal. For 56 >/ 9 tihe effect of t.he s ixth: mode is so small.X that it is unobservable even on the seale used in Figure 13. It should be noted that for 54>d, the effect of the fouirth mode is also qu.ite smalL~ According to Table 5, the dimensions

1.0 Fi(X x - F2Q8 F22 F26 6;_00 F44 P44 Y(X/ X/ Xx-+4 00nf 0.4;~~~~~~~~~~4 F,,~~~~I~ 24Pz 1 ~ ~~ ~ ~ ~ ~ ~ ~~~5 1 0est 20 Ii Figure 12. ~~~~~~Fn (nn ) Fnn Fiur 1.The Behavior -of the Functions and -- Pnm. xx~ nn

-73+.10 CURVE I: S6: S4: S2 CURVE 2 Se S2 +.08 1= t CURVE3: Si 54- 1 CURVE4: S6= S4= I +.06 A S2 +.O 4 +.02 -.02 2 5 10 20 Figure 13. The Effect of Modal Interaction on the Critical Value of S2.

of the reactor hbave to exceed about 71 migration lengths in order that S4 become less than unity. In order that S6 become less than unity, the reactor dimensions have to exceed about 109 migration lengths. Therefore, unless the size of+ the reactor is extremely large, modal in:teraction has only a small. effect on stability~ Figures 14 and 15 show the relative magnitudes of N 2 N4, and N6 when the system consisting of the'three coupled modes is marginally stable. The calculations in this Chapter show that if the effects of the nonlinearities and the terms involving the neutron generation time are neglected, for FO 1 0n0PAVI' i.nteraction with the m'th mode detracts from tne stability of the nt.th mode irrespective of the stability of the m'th. For 5n >,/I interaction with the m'th mode contributes to the stability of the ngth provided that the mlth mode is predicted to be stable by the approximate stability criterion0 For very high values of the maxi-, mrnum neutron flux the effect of' modal interaction, can be described by the simple statement,~ interaction with the -malth mode contributes -to -the stab.i.lity of +.the nTth. if the myth mode is by itself stable, ot:herwise it detracts from. it~. The numerical, calculations show that unless the reactor size is extremely ilarge, the effect of modal interaction on the stability of a mode is small~

* 608* = 9S -s = Zs oi 9H pu'e "tlx JO suoT.4GTTTSO * T aalUTA sunOH b 9Z.0 *aS -x/ z3 TVOS 3W11 00i~jfmm~ 9 F1XG XfF~~I IIIII th il W iLk0H1 0t WF - FESSE IIX I I-ZfttS~t W4Xt-f~glWX3 Att Xttl Ek S X - W i t t T - h L ff- Wt0Wt4H M XXIX A t-# til r; I VT -i fill \-F i M Mitt t~~~~~~~~~~t!thtetjr HPtt'. oi t~~e 2 t~i ii ii i j ItfV -T F- t H: tt x i,, T j! rti t i 1 It t tli i ii I f i- t I - t 1r i t t i ifl-t 1 1;.I i 1 1!. 4.4 t titi i ejit t'' t t: -I J ri i~ lc -ni I 44 41TT — Tt. i ti;rf ii I~~~~~~~~~i Hi I, I II, ti W IM i. t 1 T- 1 Hli t1 1 1 It t J' [ ititif r i if J,., fit I fNtIift fI iVi 1 ~ AIV A 1 11 tt t tli i:t i i t~~~~~~~~~~~~~~~~~~~~~~i 41 f f f I f iit: t N~:t t;it ffil liUKH- 1-1 11 iif ~ ~ i tlt t f i ti i I~~t ii f i ~t t~~f-f 1 1 It tit, - I I 1-1 A If t tt t i 1 i i,~~~~~~~~~~~~~~~~~~~~~~~Ii. i I 1

-76gir i I i i i ~~~ ~ ~~~~I I I F-f O -f.02O —-.02 -- Ne~~~~~~~~~~~~~~~A 4+ f i Ollatol rN, 4 adN for S~~~~~~~~ =.409, S4 = S6 = 1~~~~~~~~~~~~~~~~~~-9 1

-77CONCLUSION The cal.clatlions in the last two chapters show that the effects of the nonlinearlties and those of modal interaction on stability are sm.all, unless the reactor size is in excess of about. 100 migration lengths. Since the dimensions of present day reactors do not exceed even 50 times their migration lengths, the approximate stability criterion derived in Chapter VII can be expected to give good predicti.ons of their stability. This means that if the approximate stability' criLterion predicts th.at periodic oscillations will occuar followring a disturbance in the shape of'the neutron flux distribution, then immediately followi.ng t+he distulrbance the oscillations will inr fact be nearly periodic, and will converge or diverge only slightl.y The fact that the effect s of the nonlinearities and those of modal interaction are small permits one to assess the effects of the former in'the absence of -the latter, and vice versa. To sun.marize the results linear stability is a sufficient condition for the stability of the even numbered modes, but i.s only a necessary condition for that of the odd numbered ones. If the sibcriticality of a mode is less than a certain value, modal interaction detracts from its stabilityo If its suberiticality is greater than this value, modal interacti.on will promote its stabi.lity provided that the other modes are stable. An increase in the reactor size decreases the subcrit;icaly onf the higher mnodes and th:ereby resul.ts in a decrease Jn -the stability of the shape of the neutron flux distribution.

Although the calculations performed in this thesis indicate that the approximate stability criterion can be expected to give good predictions of stability, the reactor designer must bear in mind that the effects of the terms involving the neutron generation time, the effects of the nonlinearities, and those of modal interaction may all detract from the stability of the shape of the neutron flux distributiono Reactor stability studies, therefore, should include an assessment of the total effect of these factors on the stability of the reactor,

APPENDIX THE PROOF OF THE FUNDAMENTAL THEOREMS OF REACTOR THEORY USED IN CHAPTER III Taking the Fourier transform of equation (22) one obtains where Q and /are the Fourier transforms of c, and i respectively, being the variable of the transform. For nonzero i, and therefore for nonzero F -zD-2 _ r + ltp - o ( ) = 2> or By the physical nature of the problem Q (wZ)must be associated with the probability that a neutron will not leak out of the reactor during slowing d.own, and therefore it must be a decreasing function of w2 The right hand side of the last equation is an increasing function of 5a and therefore the equation has a solution for only one real value of SJ~ which will be called SJo From this it follows that in the steady state ~(r) is the solution of the wave equation VP /t) 4 J2 (L) = because by taking the Fourier transform of this equation one obtains which has a solution for nonzero FCprecisely when, t2_ Now -79

-80othen if We —=, a constant, then the Fourier transform of which is Q (wg) Fai) equals which in turn is the Fourier transform of and therefore C EbJA O r/)L

BIBLIOGRAPHY 1. Weinberg, A. M. and E. P. Wignero The Physical Theory of Neutron Chain Reactors. The University of Chicago Press, 1958. 2. Krieger, T. JO and P. F. Zweifelo "tTheory of Pulsed Neutron Experiments in Multiplying Media't, Nuclear Science and Engineering, 5, No. 1, January 1959. 3o'"Reactor Physics Constants', Argonne National Laboratory, United States Atomic Energy Commission Document ANL-5800. 4. Minorsky, N, Introduction to Nonlinear Mechanics, J. W. Edwards, 1947, Ann Arbor. 5. Ward, A. G., "rThe Problem of Flux Instability in Large Power Reactors"t, Atomic Energy of Canada Ltd., Document CRRP-657, July, 1956. 6. Henry, A. F. and J. D. Germann, "Oscillations in Power Distribution Within a Reactor"t, Nuclear Science and Engineering, 2, No, 4, July 1957. 7. Randall, D. and D. S. St.John,'FXenon Spatial Oscillations"t, Nucleonics, 16, No. 3, March 1958. 8. Wick, R. S., "Space- and Time-dependent Oscillations (and instability) in Thermal Reactors due to Non-uniform Formation and Depletion of Xenont', United States Atomic Energy Commission Document WAPD-TM-138, August 1958. 9. Stewart, J. C.,'"A Generalized Criterion for Xenon Instability'., United States Atomic Energy Commission Document KAPL-M-JS-4, -81