THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING THE EFFECT OF MODAL INTERACTION IN THE XENON INSTABILITY PROBLEM Geza Lo Gyorey December, 1961 IP -5 47

TABLE OF CONTENTS Page LIST OF TABLES AND FIGURES............................. oe.o.oo iii NOME NCLATUREo o l.o.o.o.o.o.o oo oc coo..oooo.oo.ooo.oooooooeoe iv INTRODUCTIONo................................................. 1 THE REACTOR MODEL.............................................. 3 THE MODAL EXPANSION.......................................... 5 THE MODAL COUPLING COEFFICIENTS.OOOOQOOOOOOOO..O e.o OOOOOOOe.Q 8 ANALYTICAL INVESTIGATION.....o o.....o.oo....BoooOoooeooe....ao 10 NUMERICAL RESULTS o o........................................ o 15 REFERENCES o........ o..o....... o................ o.......... 21 ii

LIST OF TABLES Table Page I Approximate Correspondence Between H/M and Sn....O.. 12 II The Values of the Constants Used in the Numerical CalculationsOO O...OOOOOo........ 16 LIST OF FIGURES Figure Page 1 The Behavior of the Modal Coupling Coefficients Fij and Pij as Functions of.................. 9 2 The Behavior of the Functions Fnn(Pnn 7x/7) X/Xx + Fnn and FnnPnm 11 and " oeooooooooooooooooooooooooooooooooooooo 11 Fnm 3 The Values of AS4 for the Condition AS2 = in the Region S2 > F2,2 P2,4/F2,4- -oo.............17 4 The Values of ZAS2 as a Function of S4 and..... 19 iii

NOMENCLATURE ai(t) the time-dependent part of the i~th term in the expansion of 0 (r,t) bi(t) same as above in the expansion of 6X(r,t) ci(t) same as above in the expansion of 6I(r,t) p the resonance escape probability q the slowing down kernel r space variable t time variable x space variable Cn the criticality factor of the n8th mode of the neutron flux distribution in the reactor D the thermal diffusion constant Fij the modal coupling coefficient associated with the steady state neutron flux distribution H the thickness of the slab reactor I(r,t) the concentration of iodine-135 Io(r) the steady state iodine-135 distribution In(t) the reduced amplitude of the nith mode in the expansion of 6I(rt) 6I(r,t) the deviation of the iodine-135 distribution from the steady state L the thermal diffusion length M the migration length Nn(t) the reduced amplitude of the n1th mode in the expansion of 8 (r, t) iv

P.j the modal coupling coefficient associated with the steady state poison distribution Q(Kn) the fast leakage escape probability associated with the n'th mode Sn the reduced subcriticality of the n'th mode of the neutron flux distribution T the dimensionless time variable X(r,t) the concentration of xenon-135 Xo(r) the steady state xenon-135 distribution Xn(t) the reduced amplitude of the n'th mode in the expansion of 6X(r,t) 6X(r,t) the deviation of the xenon-135 distribution from the steady state Yx the fission yield of xenon-135 YI the fission yield of iodine-135 7 Yx + YI bi~j Kronecker's delta the fast fission factor Xx the decay constant of xenon-135 kI the decay constant of iodine-135 X Xx + XI v the average number of neutrons produced upon fission Car the microscopic thermal neutron absorption cross section of xenon-135 T the neutron generation time in an infinite medium v

O(r,t) the thermal neutron flux _o(r) the steady state thermal neutron flux distribution 00(max) the maximum value of 0o(r) 50(r,t) the deviation of the thermal neutron flux distribution from the steady state tn(r) the n'th characteristic function of the wave equation fn(max) the maximum value of tnn(r) rn the maximum effect of the poison on the criticality factor of the n'th mode Kn the n'th characteristic value of the wave equation Za the thermal absorption cross section in the clean core Zf the thermal fission cross section the reduced maximum steady state thermal neutron flux vi

INTRODUCTION In the treatment of the space and time dependent reactor stability problem of neutron flux shape variations due to xenon-135, the method of harmonics is often used. Through this method, the finite set of partial differential equations describing the problem are replaced by an infinite set of ordinary differential equations which describe the time behavior of the various modes or harmonics of distribution. There is a finite subset of equations associated with each mode. It was pointed out by Kaplan (1) that particular sets of harmonics can be found such that each subset of these equations is independent of the other subsets, that is, such that the time behavior of each mode of distribution can be described independently of the behavior of the other modes. Often, however, it is desirable to choose a set of harmonics, such as the characteristic functions of the wave equation, which are simple, well known, and not dependent on the reactor power level. The use of such a set of harmonics will in general result in an interdependence of the subsets of ordinary differential equations mentioned above, such that the time behavior of each mode is influenced by the behavior of the other modes. It is convenient to attempt to predict the stability of the shape of the neutron flux distribution by examining the stability of each mode by itself, that is, when modal interaction is neglected. The results of such an attempt may not be very meaningful however, if it is possible that an unstable system is produced when several modes, each of which is by itself stable, interact~ This work deals with the effect of modal interaction on stability0 The characteristic functions of the wave equation are used 1

in the modal expansion. The reactor model used is described in a subsequent section. In this work, some of the terminology and notations of Weinberg and Wigner (2) are used, The term "criticality factor" and the symbol C are used in place of the more usual term "effective multiplication factor" and the symbol keffo In the modal expansion the first term will be called the first harmonic or fundamental mode, this in general should be distinguished from the steady state neutron flux distribution.

THE REACTOR MODEL The reactor model used in this work is a bare, thermal reactor with stationary fuel, it is homogeneous except for the xenon poison. It is assumed that variation in the xenon-135 density is the only process through which a change in neutron flux level affects the properties of the core medium, that is, xenon poisoning is the only feedback effect. It is assumed that the xenon affects only the thermal absorption cross section. Linear theory is used: the treatment is restricted to small deviations from the steady state. The attention is focused on the stablility of the shape of the neutron flux distribution, and therefore it is assumed that the fundamental mode of the flux distribution is held constant by a suitable control system which has negligible effect on the higher modes. All numerical computations and results are given for a slab reactor with effective boundaries at x = o and x = H, and in which no variations are allowed in the y and z directions in any of the variables. Subject to the enumerated assumptions, the reactor system can be described to a good approximation by the following equations: a rt) 2(r,t) - ZD(rt) - xX(rt)0(r,t) + vcpZf j 0(r',t)q(Jr-r'I)d3r' (1) all space a X(r,t) I(r,t) + xf(rt) - x(rt) - UxX(r,t)0(r t) (2) 6t a I(r,t) = 7ief(r,t) - XiI(rt) () at -3

-4with the boundary conditions that the variables are zero at the effective boundaries. The approximation that (1) is zero is reasonable for the higher modes of distribution, provided that the reactor dimensions are not more than a few hundred times the migration length.

THE MODAL EXPANSION It is convenient to separate the variables into steady state and time dependent parts: 0(r,t) = 0o(r) + 50(r,t) (4) X(r,t) = Xo(r) + 5X(r,t) (5) I(r,t) = Io(r) + 8I(r,t) (6) Substituting these into the equations of motion, subtracting the steady state equations from the resulting ones, and then neglecting the terms involving the product 506X, one obtains a set of equations describing the behavior of small deviations from the steady state. The time variables are now expanded in the following manner: co 80(rt) = Z ai(t)4i(r) (7) i =1 00 bX(r,t) = Z bi(t)>i(r) (8) i=l 8I(rt) = Z ci(t)4i(r) (9) where V72i (r) + Kii(r) = 0 (10) f *i(r)j(r)d3r =bi (11) reactor volume and the 8's satisfy the same boundary conditions as the variables. Multiplying the resulting equations by *n and integrating them over the reactor volume, one obtains the desired infinite set of ordinary -5

-6differential equations. These equations are cast into a convenient dimensionless form by defining the following quantities.. n(max) Nn(t) _0n(max an() (12) Xx0(max)ant Xn(t) Xx7En(max) bn(t), (13) Z~fio (max) In(t).x!n(max) cn(t), (14) YI fefo (max) the reduced amplitudes of the ntth modes of the neutron flux, xenon, and iodine concentration deviations from steady state; -x - 0o(max), (15) ~X the reduced maximum steady state flux, Fij ax f 0o(r)ti(r)*j(r)d~r (16) ijXx reactor volume Pij _ x X0 (r) (r) * (r) d3 r (17) reactor volume the modal coupling coefficients associated with the steady state flux and poison distribution respectively; 7ECf,YZf, (18) Za (l+L2K2) E Kn the maximum effect of the steady state poison on the criticality factor of the n'th mode (so that rnPnn is the actual effect on it); vcZp;a(K+-n) Cn m~~~~~~~nn t ~~~~~(19)

-7the criticality factor of the n'th mode; S 1-Cn (20) rn the reduced subcriticality of the n'th mode, and finally T -- Xxt, (21) the dimensionless time variable, In terms of the newly defined quantities, the equations of motion can be written in the following form, 0 ~SnNn + FnnXn + n[PinNi + FinXi], (22) dXn _ I In + (x _ Pnn)Nn - (l+Fnn)Xn dT 7? - Z [PinNi + FinXi] (23) iAn dIn = XI (NnIn), n = 1, 2, 53 e cX (24) dT Xx The influence of the i'th mode on the n'th one is represented by the expression within the brackets in Equations (22) and (23).

THE MODAL COUPLING COEFFICIENTS The values of the modal coupling coefficients Fij and Pij, as indicated in (16) and (17), are influenced by the shapes of the steady state flux and poison distributions, the latter distribution being uniquely determined by the former. The shape of the flux distribution depends on the maximum flux level as well as on the reactor size. In the poisoned reactor the flux distribution is flatter than the fundamental mode of the expansion used here, but for the combinations of reactor sizes and flux levels considered in this work, it is not much flatter. For numerical computation of the coupling coefficients, therefore, the shape of the steady state flux distribution is approximated by the fundamental mode I1o A flattening in ~o, and therefore also in Xo, increases the values of Fij and Pij for i=j (which might be called the self-coupling coefficients), and decreases the values of the coefficients for iAj, because the I's are orthogonal functions. The approximation mentioned above, therefore, results in a slight apparent strengthening of modal interaction. The values of some of the coupling coefficients are shown in Figure 1 as functions of the flux level0 When the nonlinearities are neglected, as is done here9 there is no interaction between an even and an odd harmonic0 The absolute magnitudes of the coefficients fall off quite rapidly as li-jl increases, and. therefore one may expect the strongest interaction to appear between the i'th and j Ith modes when j-=i+-2. -8

2.00 ~.50..200 Lim' F, 5F \ n-. - F,2,l'oo...0 f2,2.,- / /m-, n n+1.00 5 * FD P,., Lim Pn,..50 00~ n a.20 le ~~~~~~~~~Lim -Pn~~.10.05 F.02 an1Pa3 nP2 to of "1J/~' r-~-P "J.*Jln;:pnn4.01.10.20.50 1.00 2.00 CD 5.00 I0.00 20.00 50.00 I00.00 I ii I 12 -2 I $ 13 14 7x 10 7x I01 Qo(max)CM2 SEC 7x 10 7x I0 Figure 1. The Behavior of the Modal Coupling Coefficients Fij and Pij as Functions of (.

ANALYTICAL INVESTIGATION An approximate stability criterion can be derived for the higher modes when modal interaction is neglected. Under this condition, the system represented by Equations (22), (23) and (24) is stable if Sn> Fnn(Pnn -, n = 2, 3,.=, o (25) X/Xx + Fnn This criterion predicts marginal stability for the n'th mode when Sn equals the right hand side of (25); this value will be referred to as the critical value of Sn as predicted by the approximate stability criterion, This critical value is shown as a function of flux level in Figure 2. It is desirable to have a rough indication of what this criterion means in terms of reactor size rather than in terms of subcriticalityo For nAl, but not very large, and for the range of reactor sizes of interest here, one may write (3) as a rough approximation S M2 ( - ) e (26) rn For a slab reactor with 2n -.03: S n ~2 (nZ1)(M)2 (27).03 Table I shows this approximate correspondence between reactor size and Sno It should be noted here that the exact relationship between Sn and reactor size involves terms which are dependent on the power level. The analytical investigation carried out was directed at answering the following question. If the n'th mode is predicted to be marginally stable by the approximate stability criterion given above, then will -10

1.0 0.8 X/X + F44 n.oo X/XE + F 0.6.P Lim 4 - Pn n+2 0.4 I 0.2 0.5 1 2 CD; 5 10 20 50 100, 12 I2,~_ c' I, 1 13 2.1 x 10 7xI0 o (max)CM SEC 7xIO 7 x Figure 2. The Behavior of the Functions Fnn(Pnn - 7x/Y) and FnnPnm /~x, + Fnn Fnm

-12TABLE I APPROXIMATE CORRESPONDENCE BETWEEN H/M AND Sn H/M: 31 35 41 50 70 99 140 S2: 1.8.6.4 e2 1.o05 S4: 5 4 5 3 2 1.5.25

-13-the combined system of the n'th and m'th modes be stable or unstable? In order to attempt to answer this question, the equations of motion (22), (23) and (24) of the n'th and m'th modes were combined under the assumption that Fij = Pij = 0 except when i=n or i=m and j=n or j=m. The resulting fourth order differential equation, which is too long to exhibit here, was subjected to a rather tedious examination in the light of the Hurwitz-Routh stability criterion (4) under the condition that the n'th mode is predicted to be marginally stable by the approximate criterion (25). For the reactor model used here, this examination yielded the following results. For Sn < FnnPnm/Fnm the combined system is unstable, irrespectively of the stability of the m'th mode. The behavior of the right hand side of this inequality as a function of flux level is shown in Figure 2. For Sn > FnnPnm/Fnm the combined system is stable if the m'th mode is by itself sufficiently stable, otherwise it is unstable. Except for a very small region where Sn, FnnPnm/Fnm, in order that the m'th mode be by itself "sufficiently stable" in this sense, the requirement Fmm(Pmm - jx/) (28) X/Xx + Fmm is necessary, but not sufficient. In other words, the value of Sm must be larger by some amount than that necessary for marginal stability as predicted by the approximate criterion (28). A necessary and sufficient condition simple enough to be useful was not found for this "sufficient stability" of the m'th mode. The numerical computations discussed in the next section, however, indicate that the value of Sm does not have to be significantly greater than that required by condition (28)~

-14For very large values of Sm, the equation of motion of the combined system becomes independent of Sm and the effect of modal interaction approaches a nonzero limit. This fact is quite significant. According to Equation (22), an infinitely large Sm means that Nm, the m th mode of the flux distribution is not allowed to vary. Nevertheless, it was found that interaction with the m'th mode has a nonzero effect on the stability of the n'th mode. The reason for this is that although the m'th mode of the flux distribution is not allowed to vary, the m'th mode of the xenon distribution will be forced into oscillation by variations in the n'th mode of the flux distribution. In other words, even if the oscillations in the flux distribution are purely in a single mode, the oscillations in the xenon distribution will involve all the other modes. This result agrees with those reported by Pearce (5), This means that even if the magnitude of the fundamental mode of the flux distribution is held constant by some control system, in examining the stability of the third mode, interaction between it and the fundamental must be taken into account,

NUMERICAL RESULTS Since the least stable higher mode is the second harmonic, the numerical calculations were centered on it. An exploration, by means of analog computer simulation, of the behavior of the system consisting of the second, fourth, and sixth harmonics coupled together indicated that the effect of the sixth harmonic on the second is negligible unless the reactor size is in excess of a few hundred migration lengths. Such large reactor sizes are not considered here because they are beyond the range of validity of the approximations used. Therefore, in further numerical computations, which were done on a digital computer, the system consisting of only the second and fourth modes was considered. Table II shows the values of the constants which were used in the calculations, For a given reactor composition the power level and the reactor size uniquely determine the values of S2 and S4, or alternately, the reactor power level and the value of S2 determine the value of S4o In the computations, however, the value of S4 was varied in the full range S2 4 S4 co (effectively). Although this led to some computations which may not be physically realizable, the results are quite interesting and do contribute to the understanding of the problem. In the previous section, it was stated that in the region where S2 > F2,2 P2,4/F2,4, interaction with the fourth mode will enhance stability provided the value of S4 is sufficiently large. Figure 3 shows the results of the search for values of S)L in this region such that the critical value of S2 as predicted by the approximate stability -15

-16TABLE II THE VALUES OF THE CONSTANTS USED IN THE NUMERICAL CALCULATIONS ax = 3 x 10-18 cm2 Xx = 2.1 x 10-5 sec-1 XI = 2.9 x 10-5 sec-1'x =.003 7I =.061

.10....08.06 04 H.02 24 I 2 5 10 20 Z 50 100 200 500 1000 I12 I13 -2 -1 14 7 5 7xIO 7x1O 4o(max)CM SEC 7x 10 7x 10 Figure 3. The Values of AS4 for the Condition AS2 = 0 in the Region S2 > F2,2 P2,4/F2,4.

-18criterion (25) is also the critical value for the combined system of the second and fourth modes, The plot shows ZS4: the difference between the values of S4 found by this search and the critical values of S4 as predicted by the approximate stability criterion. This difference is relatively small everywhere in this region. Table I shows that in this region the value of S4 is considerably larger than the critical values found by the search described above. Therefore, in this region, if the second mode is predicted to be marginally stable by the approximate stability criterion, the combined system of the two modes considered here will be stable. The remainder of the numerical calculations consisted of the search for values of S2 for which the combined system is marginally stable. The flux level and the value of S4 were varied over wide ranges. Figure 4 shows the computed values of AS2: the difference between the values of S2 found by this search and the critical values of S2 as predicted by the approximate stability criterion. The results indicate that if the reactor dimensions are not very large, less than about 60 times the migration length for the numerical values used here, second mode instability will set in at a high flux level: above 1013 neutrons/cm2/seco Interaction of the second mode with the fourth in this case enhances stability. The actual critical value of S2 is percentagewise only slightly different from the critical value obtained when modal interaction is neglected. If the reactor dimensions are large, 70 times the migration length or more, second mode instability sets in at flux levels below 1013 neutrons/cm2/sec. Interaction with the fourth mode in this region

.05.04...03.02 /-s4 S s-,, /.OI IS4... 2 -01 1 I / \v I hII 1I I t 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 0 __ _ S4 =D.01 2 F24, F2~ 4 -4 I.02.1.2.5 I 2 CD 5 10 20 50 100 200 500 I00 I Oji I 12 -2 -I x I013 14 I15 7x10 7x10 o(max) CM SEC 7x10 7x10 Figure 4. The Values of ZS2 as a Function of S4 and 0.

-20detracts from stability. As the reactor dimensions are made larger, the difference between the actual critical value of S2 and that predicted when modal interaction is neglected becomes percentagewise quite large. Thus the stability of a large reactor is seriously overestimated if modal interaction is neglected~

REFERENCES 1. S. Kaplan, Nuclear Sci. and Eng., 9, 357-361 (1961). 2, A. M. Weinberg and E. P. Wigner, The Physical Theory of Neutron Chain Reactors. The University of Chicago Press, Chicago, 1958. 3. Go L. Gyorey, On the Theory of Xenon Induced Instabilities in Neutron Flux Distribution, PhoDo Thesis, The University of Michigan, Ann Arbor, 1960. 4. E. A. Guillemin, The Mathematics of Circuit Analysis. John Wiley and Sons, New York, 1949. 5. R. M. Pearce, Nuclear Sci, and Engo, 11, 328-337 (1961). -21