Preprint of a Paper Submitted to Nuclear Science and Engineering for Publication The Concept of Multiple-Input Zero-Power Describing Functions in Nuclear Reactors* by A. Z. Akcasu and C. M. Bost, Jr. Department of Nuclear Engineering The University of Michigan Ann Arbor, Michigan 48105 No. of Pages - 32 No. of Figures - 6 *Based on work performed under a grant from The Institute of Science and Technology, The University of Michigan

Send the galley proofs to: A. Z. Akcasu and C. M. Bost, Jr. The University of Michigan Department of Nuclear Engineering North Campus Ann Arbor, Michigan 48105

ABSTRACT The concept of multiple-input zero-power describing function is introduced and employed to obtain the describing function at high power using feedback circuit theory. Numerical application of the describing function to nonlinear stability analysis is presented and the various definitions of the describing function are discussed. iii

Introduction The describing function concept is introduced in stability analysis as an attempt to extend the methods of linear stability based on the transfer function concept to nonlinear systems. Essentially it is an amplitude dependent transfer function which accounts for the nonlinearity of the input-output functional relationship of a physical system. The first order correction to the reactor transfer function due to the nonlinearity in the neutron kinetic equations was obtained perhaps for the first time by Rumsey(1) in 1947. Nelkin(2) in 1955 generalized Rumsey's approach to include linear feedback processes, but his analysis did not take into account the change in the mean neutron density from the equilibrium value (3} in the presence of oscillations. Sandmeier(3) in 1957 presented a corrected (version of Nelkin's analysis. Akcasu(4) in 1958 showed that, in the firstorder approximation to the instantaneous reactor period, the logarithm of the neutron density plays the role that the neutron density plays in linearized kinetics. He also defined and calculated the "nonlinear gain" of a reactor as the ratio of the amplitude of the fundamental Fourier component (relative to the average of the periodic power oscillations) to the amplitude of the sinusoidal reactivity input. Smets(5) was perhaps the first to introduce the concept of "describing function" in reactor stability analysis. He calculated the magnitude of the zero-power describing function without accounting for the effect of the second harmonic and the variation of the average power level and consequently concluded that the magnitude was generally larger than thatof the low-power transfer function. This conclusion was at -1

-2variance with that obtained by Akcasu(4) for the nonlinear gain. The discrepancy was due to the oversimplification in Smets' analysis as clarified (6) by Wasserman) in 1962. Wasserman derived low and high power describing functions using a perturbation analysis which was equivalent to, but more systematic than, the procedures of Rumsey, Nelkin, and Sandmeier. He compared the theory to the results of an oscillator experiment on the SPERT Reactor and demonstrated the nonlinearities. He also pointed out the lack of uniform convergence in frequency when the low power describing function is expanded in powers of the amplitude of the reactivity insertion. In 1964 and 1965 Smets(7'8) extended his earlier analysis removing its shortcomings and presented a block diagram to study stability and limit cycles. The low power describing function he used in this diagram was obtained by performing a Fourier analysis in the logarithm of the reactor power in the presence of linear feedback and contained the effect of the second harmonic as well as the feedback parameters. In this sense his circuit diagram did not decouple the feedback effects from the nonlinearities of the neutron kinetics. Akcasu and Shotkin(9) in 1967 investigated limit cycles and their stability in various reactor models in terms of a describing function which was valid for larger amplitudes than those previously reported. In 1968(10) they proposed a zero-power describing function based on the 1967 analysis whose convergence did not depend on frequency. A note by Vaurio and Pulkkis(11) in 1970 extended the analysis of Akcasu and Shotkin by retaining the second harmonic in the logarithm of the reactor power. The zero-power describing function they obtained was expressed as a ratio of two infinite series of modified Bessel functions. They compared the various approximate expressions of the zero-power describing function numerically and determined the range of

-3frequency and reactivity insertions for the validity of these approximations. It can be observed from the above literature survey that the concept of describing functions in reactor dynamics has evolved continuously since 1947 with each paper extending and refining the previous one. However in 1967 Tan(l2) asserted that all the previous experimental and theoretical work was incorrect due to the definition of the describing function with respect to the average power level in the presence of power oscillation rather than the equilibrium power level prior to the periodic oscillations. A closer study of the development of the describing function concept reveals that it has not been possible in the past to obtain the high-power describing function by combining the zero-power describing function and the feedback transfer function in a closed loop and applying only feedback circuit theory. Such a derivation is desirable because it facilitates the calculations, clarifies the ideas by allowing the discussion of nonlinear effects and feedback effects separately, and demonstrates the utility of the zero-power describing function. The main objective of this paper is to show that such a derivation is possible if a new concept, the multipleinput zero-power describing function, is introduced. The application of the feedback circuit constructed by this approach will be demonstrated by plotting the amplitude-dependent Nyquist diagrams, and the existence and stability of limit cycles will be discussed on the basis of these plots. Another objective of the present paper is to clarify and hopefully (12) settle the controversy introduced by Tan's work(2) concerning the definition of the zero power describing function. In order to achieve these objectives, particularly the second one, we must review briefly the basic concepts and calculational techniques. For simplicity we shall neglect the delayed neutrons in these discussions and comment on the effect of

-4delayed neutrons at the end, without further explanation while presenting the results. Definitions of the Zero-Power Describing Function In the absence of delayed neutrons the point kinetic equations can be written as (i//(P/Jt) = k t) Prt) (,, where k(t) is the reactivity insertion measured in dollars, f is the prompt generation time, and B which is introduced here to convert the reactivity to dollars is the effective delayed neutron fraction. This equation is a linear, homogeneous equation with variable coefficient whose solution is determined within an arbitrary factor. To define the zero-power describing function, we consider a sinusoidal reactivity input as tkr)= _ t s$r xt (2) and determine the periodic power response. Since we are interested only in the stationary solution, the origin of time is irrelevant and can be chosen as desired. It is easily shown that the solution of (1) is P(t)-C c Mpp3-( k/0W) Cos a ] (3) where C is the arbitrary factor mentioned above. This constant can be calculated in terms of the value of P(t) at some t = to or in terms of the average power PAV defined by

-5PAv (l/T)fcJb P(t) 0 where T = (27/w), the period of the sinusoidal reactivity. In the latter case which will be used below, P(t) is given by Pt) = [PAV/T(s /eo )jJ XP/- (P k ciw) C0S W J (4) where Io(x) is the modified Bessel function of zeroth order defined by ~0~~~~~ o (X-(I/ ) 1 jdy exp ( Cosy) (5) If we just observe the periodic power oscillations induced by a sinusoidal reactivity insertion without any knowledge of the initial conditions, then we can use (4) to analyze the harmonic composition of the oscillations. Note that the equilibrium power level prior to the application of the sinusoidal reactivity insertion does not enter the analysis. We may also relate the arbitrary factor C to the initial conditions.an iniz/ial value if they are available. In this case we regard (1) a Vproblem in which the reactor power is assumed to have an equilibrium value Po prior to t = 0, at which time the sinusoidal reactivity input is turned on with a phase angle e, i.e., The solution is readily obtained as P(t)- PI eap[(p^k//os)cose] e.p[ (j /ew) oowt+&)] (6)

-6It should be noted that (4) and (6) represent the same periodic function, i.e., the same harmonic composition. The phase angle e in the time-dependent exponential factor in (6) does not affect the wave form of the periodic solution. Therefore by comparing (4) and (6) (or evaluating the average of both sides of (6) directly) we can relate the average power to the equilibrium power level as PV = Q(Swe) R (7a) where Q(S?,w9)=_Io( p 8/t) go ap[(BSk/z)c6O'S0 (7b) It should be noted for future reference that (PAV/Po) depends on the initial insertion as well as its magnitude and frequency. Consequently one may obtain the same periodic power oscillations (the same _.... power and harmonic composition) starting from different equilibrium power levels by adjusting the initial phase. The fact that PAV varies with the initial phase for a fixed Po implies that the reactor in the absence of feedback has a "perfect memory" and remembers the initial phase even after infinitely long times. We now focus our attention on the harmonic composition of the power oscillations using (4). By using a Fourier expansion we obtain: P(t) > V 2 ft SO + o CoS 2o t L T,(N, Sk Io ()N, k) ~'] -(8a)

-7where A/,- 3 /- )=1Z(Lw )| (8b) Z(iO) = (0/tia) being the zero-power transfer function in the absence of delayed neutrons. We can express this expansion in a more compact form by introducing the complex amplitude X of a sinusoidal function defined in general by A sin (c 4-+&)_ X p(pt~l )+ cc (9a) where -x (A/2 c) p ('( 0(9b) and cc denotes the complex conjugate of the first term. Introducing n [2 In (h S )/t Io (M so) (2 )] (lOa) with the convention 3Yn- Y-n I Soy 1(1Ob) we obtain from (8a) 400 [P (t)/Pa| X O(p(tnItI) (11) We note that yn denotes the complex amplitude of the nth Fourier component of the power oscillations relative to the average power level PAV and the actual complex amplitude of each component is PAV Yn. The complex amplitude of the input 6k sin wt is of course x1 = 6k/2i.

-8The conventional definition of the zero-power describing function is Dz ( Wx)i-(i /x,) 2I,(NSlk) kII(,Sk) (12) which is the ratio of the complex amplitude of the fundamental component of the power oscillations (relative to the average power) to the complex amplitude of the sinusoidal reactivity input. It should be noted that the describing function defined in (12) is a function of the magnitude and frequency of the reactivity input. When 6k - O, Dz(6k,w) approaches Z(iw) = B/tiw which is the zero-power transfer function in the absence of delayed neutrons. It can be measured in an oscillator experiment by observing the periodic power oscillations and comparing with the sinusoidal reactivity input. One does not need to know how the periodic oscillations are produced initially. Interpreting the experimental results using (12), one can for example determine the zero-power transfer function using N1 and (8b) without using the linearized kinetic equations. (This point will be discussed later.) We can also define a describing function measuring the complex amplitude of the fundamental relative to the initial power level which existed before the periodic power oscillations were induced. Using the relationship between PAV and Po in (7a), we obtain DT(SkC,B) Q(GKk6 ))0) DZ(kSW) (13) The describing function defined in this way not only depends upon the magnitude and frequency of the sinusoidal reactivity input but also on the Because of fhis dependenice on {he inil/cd va/ue, initial valueV'~T(6k,w,e) cannot be determined by observing only the periodic power and reactivity oscillations. The initial phase angle must be known as

-9well as the manner in which the final stationary periodic oscillations are reached. As an illustration, assume?(t)= ^t sn(utt+0),t it< 1k St' j+F e) ) t o I with yto = 6k. This represents a linear increase in the magnitude of the reactivity input during the time interval 0 < t < to. We can easily verify in the absence of delayed neutrons that the factor Q in (13) depends in this case upon the reactivity insertion rate y, the insertion period to, and the phase angle e. The value of Q given in (7b) corresponds to a sinusoidal reactivity insertion with constant amplitude for t > 0. It appears from the previous discussions that the describing function DT defined with respect to the initial power level does not characterize the periodic response of a reactor to a periodic input in terms of the kinetic parameters of the system, but rather it depends critically on the initial conditions. Since different initial conditions may give rise to the same periodic power oscillations with the same periodic reactivity input, it cannot be measured by observing only the periodic oscillations. Because (4) of these drawbacks, it was abandoned as early as 1958 in the literature(4) (12) In 1967, Tan 12 proposed this describing function as the "correct" definition and asserted that all the previous theoretical and experimental work which did not use this definition was "wrong". It seems that this assertion is based on the discrepancy between the zero-power transfer function and the zero-power describing function when the latter is calculated using the "linearized" point kinetic equations.

-10Connection Between the Describing Function and the Zero-Power Transfer Function The kinetic equations are linearized by introducing the incremental power p(t) from some reference power level PR as P(t) = PR + p(t) and then ignoring k(t)p(t) in the original kinetic equation. In the absence of delayed neutrons the "linearized" kinetic equation reads as p(t) = k(t)(I/e)PR. The various choices of PR will be discussed presently. Assuming that the reactor is at an equilibrium level Po prior to t = 0, the solution of this equation is obtained in the transformed domain as -P /R =- S) Z(s)+(/$s)[(o/P R-jX (14) where p(s) and k(s) are the Laplace transforms of p(t) and k(t) respectively and Z(s) = B/ts. In the presence of delayed neutrons (14) is still valid provided Z(s) is the conventional zero-power transfer function with delayed neutrons, i.e., Z(S)=[s( +-)1 (15) The response to a sinusoidal reactivity input of the form k(t) = 6k sin (wt+e) which is introduced at t = 0 follows from (14) after long times as fPE)/PR = 8,J?|l Z( )|]S<in(^ (9 B N + (J3iSk/e ) cos( +[(p,/p) -l] (16)

-11where (/ )=-[ (/)+, (at/y)] and Z= CS [Z(ci )] In (14) and (16) the reference power PR is still aribtrary. Conventionally PR is chosen as the equilibrium power level so that (14) reduces to p(s) = P Z(s) k(s) which yields the conventional definition of the zero-power transfer function as the ratio p(s)/Pk(s). In this case the last term in (16) vanishes so that the response to a sinusoidal input becomes after long times P()/ = I iZCc')Ii SZn( +ti +)n +(+ cs e/ew ) (17) The last term accounts for the change in the equilibrium power level in the presence of periodic oscillations. The average power level in the presence of oscillations is, = c+(sks )j (18) which depends on the initial phase angle e. We can compute the describing function from (17) using the definition (12) relative to PAV as Dz Z (Zf )/[ Sf cos e/e a)] (19) which differs from the zero power transfer function. One would expect the definition of the describing function to yield Z(iw) when it is evaluated from the linearized point kinetic equations. This seeming contradiction led Tan(12) to conclude that the conventional definition of the describing function must be "wrong". Tan also observed that DT = Z(iw) as it should if the describing function is defined relative to the initial equilibrium power level PO.

-12However, two points are overlooked in this conclusion. First, the linearized kinetic equations are valid only when |p(t)l<<Po which requires IZ(iw)16k<< 1 to hold for all w. This inequality becomes (6Ak/t*w)<<l for small frequencies. Therefore, the denominator in (19) must be replaced by unity to be consistent with the limitations of the linearized kinetic equations. If the term (A6k cose/t*w) in the denominator is not small compared to unity, the linearized treatment is not valid and one must use the original kinetic equations. The second and more important point which is perhaps overlooked in Tan's conclusion is that the choice of the reference power level PR is not restricted to P. If one is interested in the stationary periodic solution after long times rather than the behavior of P(t) for small times, it is more reasonable to choose the reference power level as the average power of the periodic oscillations. In this case, (16) reduces to ('t)/P,, = 81 I\Z(f)t \Sil (otLB6+ <) (19a) The average power level is related to PO by equating the last two terms in (16) to zero: I PAV POi- ol( FS1? e/e S)] (19b) It is obvious now that the definition of the describing function relative to the average power precisely yields DZ = Z(iw) whereas the definition relative to the initial equilibrium yields DT=Z(iw)[-I (pSk s e/wS)4 (20)

-13which differs from Z(iw). It is clear that both choices of PR are identical within the limitations of the linearized treatment which is the limit as 6k-*O. It should be clear from the above discussion that the discrepancy (12) observed by Tan 12 is due to approximations inherent in the linearized treatment, and emerges only if one attaches significance to second order terms in 6k which must be neglected in a consistent linearization procedure. It does not imply, as claimed by Tan, that the conventional definition of the describing function is "wrong". Actually it is meaningless to label any "definition" as "correct" or "incorrect", because definitions can only be "suitable", "convenient", "useless", etc. depending on the purpose. We shall show in subsequent sections that the conventional definition of the zero-power describing function arises naturally when constructing the describing function at power in the presence of feedback usinq the feedback circuit analysis. Describing Function in the Presence of Feedback In the presence of feedback which we assume to be linear, the kinetic equation (1) is replaced by (g/V) ( dP^~)/8t)-k) J<Tdu G(U) [Pt-u)-Po] (21) where G(u) is the feedback kernel, k(t) is the external reactivity which we again choose as k(t) = 6k sin At, and Po is the equilibrium power level determined by F^H Co~= k (22)

-14where H(s) is the Laplace transform of the feedback kernel, H(O) is the power coefficient, and ko is the externally introduced positive constant reactivity. This equation is nonlinear and inhomogeneous in contrast to (1) in the absence of feedback. Consequently the periodic response to a sinusoidal reactivity input does not contain any arbitrary factor. The average of the periodic power oscillations can easily be shown to be equal to the equilibrium power level for any periodic reactivity input with zero mean value (this result is true only in the absence of delayed neutrons as we have assumed, for simplicity, in these discussions). Since PAV = Po we can define the describing function unambiguously as the ratio of the relative complex amplitude of the fundamental of the output to that of the input. Assuming that the feedback is a low-pass filter (H(inw)? 0 for n >?) and ignoring the second harmonic in the logarithm one can obtain the describing function at power approximately as D(S P )= Dz(S )) -P H() DZ (23) where Dz(6k,w) is the zero-power describing function defined in (12). This can be interpreted as a closed feedback loop as in Figure 1 where yl = xlDz, sk xl =-2 + PAV Yl H(i), and ko = PAV H(o). We note that the describing function at power reduces to the zero-power describing function when Po is sufficiently small. Although simplified by various assumptions, the above analysis clearly shows what criteria should be used in defining the zero-power describing function. The actual purpose in introducing the zero-power describing function is to separate the nonlinearities of the neutron kine'tics from the feedback effects so that nonlinear stability and the response of the reactor

-15to a sinusoidal input in the presence of feedback can be investigated using the methods of linear feedback circuit analysis at given amplitude. If we judge the usefulness of a given definition of the zero-power describing function according to this criterion, then (23) indicates that the conventional definition of the describing function appears to be the suitable one to be used in feedback systems. The existence and stability of limit cycles (self-sustained periodic solutions) can be discussed by equating the denominator of (23) to zero i -P0 H(a ) Dz(s w) (24) which once more demonstrates the utility of the zero-power describing function concept. We note that the describing function at power cannot be obtained as easily as indicated in Figure 1 when the second harmonic is taken into account. This is possible only if we introduce the concept of multipleinput zero-power describing function which is the main contribution of this paper. Before we discuss this concept, the following remarks on delayed neutrons seem to be in order. Effect of Delayed Neutrons on the ZeroPower Describing Function When the delayed neutrons are taken into account, the kinetic equation cannot be solved exactly for a sinusoidal input. The presence of the delayed neutrons change drastically the character of the reactor response to periodic inputs. One finds that a negative bias in the periodic reactivity input must be provided to maintain a periodic output(4'13). Therefore, one has to consider a sinusoidal reactivity of the form

-16k6) Klko+S sln s t in constructing the periodic power oscillations, and to determine the appropriate negative bias ko as a function of 6k and w. In the lowest approximations ko is given by k, __(i, )2(I/2 ) Re ZRta))25a (25a) (14) Several approximate techniques, such as the prompt-jump approximation, (12) (e)7) WKB approximation, Fourier analysis of the logarithm of the power (, or Fourier analysis of the power response combined with a perturbation expansion in sk6 ) are available to analyze the power oscillations. As an example, we present the results of the last two methods for the zero-power describing function. Ignoring the second and higher harmonics in the logarithm of the power, Akcasu and Shotkin(1) obtained the zero-power describing function as D. -.......! (25b) 5k TOT(NI Sk) where N, e^ },)-Z, + [(st)[4Z2(Zl')- 2 ] \ ) When the delayed neutrons are neglected, (25b) reduces to (12). Vaurio and Pulkkis(ll) compared this result with a more accurate calculation including the second harmonic and found that the two agree with better

-17than 1% accuracy for frequencies above 102 cps and 6k up to 0.50. The perturbation analysis yields a power series expansion of D(6k,w) asC' Dz, (k,)u) _Z, I i^l [(w)/]4,a,-z(zl) i SX ) (26) As pointed out in the introduction this power series expansion is not uniformly convergent in w whereas the previous result is valid for the entire frequency range. However, the power series expansion gives the describing function correct to order (6k). The analysis of the power oscillations in the presence of feedback is more complicated when the delayed neutrons are taken into account. We shall construct the describing function in this paper using the multiple-input describing function and the closed loop feedback circuit theory and compare the results to that obtained by Smets(8) and Akcasu, et al.(4) Multiple Input Zero-Power Describing Function The power response of a reactor system with delayed neutrons but without feedback can be related to the reactivity input using the reactor kinetic equations (///?) Pt)= [krt)- l]P(t)+4JU" D ) c( -u (27) where D(2)- i a,, ep(-;i u) Upon substitution of infinite Fourier expansions for k(t) and P(t) as 400 t)-Z XnTp(n t)o )?= -OO

-18+00 +7-oo p(t)-PAVL ~ X Se8x(p'nt)W where xn = ( AV /2c ) exp( r,,7 ) J x _ x Xo= o into Equation (27) and equating terms of equal power in eiWt we obtain 0oo y _Zr Z Xn r-Ymn n=-O)t,+ +2' (28a) )= -00 where Zm = Z(iwm). Evaluating (28a) for m = 0, and noting that (1/Z ) = 0, we immediately obtain the negative bias required to maintain periodic power osci 11 ations: ko — 2Z Re xnynl 0k^_-2, Re.tn n] L(28b) This set of equations can be solved by a perturbation expansion(6 Yn = J /, n=+,l +2,, 0t tI to obtain

-19(,) Ym = Zr xr (29a) 4 w (z) 7 - ): =Zn _Z Xn Ym- n n-co (29b) +- 00+00 Y(3)-Zm =" -) -^ 0 X* /) (29c) [n — oo Xn-o o (n0 ojmn) (n o) This set of equations represents the basic equations of the Fourier analysis combined with a perturbation expansion. It enables one to construct the Fourier components of the power oscillations when the reactivity input is a periodic function of an arbitrary wave form. In the case of a pure sinusoidal reactivity input we have xl = Ak/2i and xn = 0 for n > 2. Equation (29a) indicates that only the fundamental component is excited in this case if we retain only the first order terms in the perturbation expansion, i.e., Y1) = 0 for all Iml > 2 and y = Z1 xl. This approximation corresponds to the "linearized" treatment mentioned earlier. Substituting this result into (28b), we obtain the negative reactivity bias in the lowest approximation as(4) X0 —2 IX, Re. (0) When the second order terms are retained in the computation of the complex amplitudes, we find from (29b) that only the first and second harmonics are

-20(2) f excited, i.e., y() = 0 for Iml > 3. In particular, the amplitude of the second harmonic follows as (2) 2 X31 y2 =zz x1 (2) 2 We may interpret y 2/ 2 Z X ZX as the internal harmonic generation in reactivity due to the nonlinearity of the neutron kinetics (see Figure 2). It follows that more higher harmonics are excited when the complex Fourier amplitudes are computed with increasing accuracy. Using (29) one can in fact show that the nth harmonic is given by n- =Z,2 Zn Xl (32) which is the n order in the reactivity amplitude. The describing function presented in (26) is obtained by calculating (yl/x1) using (29) correct to the second order. In a reactor with feedback, the higher harmonics of the power oscillations, generated internally by the sinusoidal reactivity input, are fed bad to the input through the feedback mechanism. Therefore, we must recalculate the zero-power describing function in the presence of these higher harmonics in the reactivity. In particular, we must consider the second harmonic if we desire to investigate the nonlinear response as well as the sustained oscillations of a reactor with feedback, correct to second order terms. This leads to the concept of "dual-input" zero-power describing function. It is defined as Dz(X1, x2, m) - (y1/x1) when the reactivity input contains first and second harmonics as (tf)= X~ [ x, eXp(itt)i CC I+ [X e^4p(2at)+ C 1 (33)

-21where x1 and x2 are independent and may be of the same order of magnitude. Using (29) and calculating the complex amplitude of the fundamental keeping the first three terms in its perturbation expansion, i.e., yl = yl) + l() + 3), one obtains Dz(xi) X2 =Z, 1+iX( Z,(Z2-2ReZ )+x (xP 7x,)(Z +ZZ 2 W' + X22 [^3 (ZZ2 )+Z1(z1+Z,2)-2Z, Re Z2] (34) Clearly, this reduces to the conventional single-input describing function Dz(x1, w) in (26) when x2 = 0. The dual-input zero-power describing function can be used in the interpretation of a pile-oscillator experiment with large amplitudes when the reactivity input is not a pure sine-wave. In such cases the amplitude of the second harmonic in the reactivity oscillations is usually an order of magnitude smaller than the amplitude of the fundamental. This is certainly true when the second harmonic is generated internally through the feedback mechanism as implied by (31). Referring to Figure 2 and performing a simple linear circuit analysis, one can in fact show that(4 x2 in this case is related to xl by X2= PAVZ H2 L2 x (35) where H2 = H(2iw), the feedback transfer function evaluated at 2w, and L2 is the power transfer function again evaluated at 2w, i.e., Z( -I (36) L(z Z )[1- V'^ Zr^ (36)

-22The average power in the presence of periodic oscillations increasesfrom the equilibrium value Po = ko/H(O)l in order to produce the appropriate reactivity bias calculated in (30), i.e., v= P,-i[ 2 2 IX Re /itI(o) 37) When x2 is known to be of the second order in xl as shown above, we can neglect the terms proportional to 1x2 2 in (34) and obtain the following simplified dual-input describing function: aD (x/, = 2, 1+ x, 2z (2-2 Re Z,) +x2 x x z z, The effect of the second harmonic input on the zero-power describingfunction has been investigated numerically using (34) for |x1/x21 = 5.0. The result is depicted in Figure 3 where we also present for comparison the zero-power transfer function and the single-input zero-power describing function calculated from (26). At a frequency w = 3.2 X 10-2 radians/sec. the phase difference in the single- and dual-input describing functions is about 20 degrees with a significant magnitude difference as well. Figure 4 shows the effect of the relative phases of the first and second harmonics on the dual-input describing function characteristics. Although in practical applications the ratio |x1/x21 would not be expected to be less than 5.0, the graphs are plotted for [xl/x21 between 1.0 and 5.0 to dramatize the effect of the second harmonic on the phase and amplitude of the zero-power describing function.

-23Construction of the High-Power Describing Function The response of a reactor with feedback to a sinusoidal reactivity insertion k(t) = 6k sin wt can be described by a high-power describing function. As in the case of no feedback, it can be defined relative to the average power PAV in the presence of power oscillations, although the equilibrium power level P could also be used in the definition without having the difficulties encountered in the zero-power describing function. Thus we adopt the following definition at present D(k ) PA^ -,,) - Y,,/( e /2 C j9) We note that no external negative reactivity bias is required to maintain periodic power oscillations when feedback is present because the bias is generated internally by an increase in the average power as discussed in (37). Using the block diagram shown in Figure 2 and feedback circuit analysis, one can easily obtain the high-power describing function as Dftb~,^Pa,)= D, I^ir^>w][\- P,v H^w) Di(x, ul] 4 D (S E ) ~, PAV') = DZ (X ) X2, L ) ) EiAV 14 (a 6()) (X )~2 X ]40) where x2 is to be substituted from (35) into DZ, i.e., Dz X',,X x )=Z|,lx- ZZ+L2(PAv!Z,)-2 ReZ] (4,) D2(x,) Z Z I (41) The average power PAV in (40) is determined from (37) which can also be written in terms of the externally applied positive reactivity as ko+ H(o)PAv — 21x, 1 ReZi (42)

-24We see from Figure 2 that xl in these equations denotes the complex amplitude of the fundamental of the net reactivity and is related to 6k by X, = (Sl/Z -PV HI DZ (43a) The right-hand side of (43a) contains Ix112. Although it can be solved exactly for xl, we can obtain an estimate of x1 by ignoring 1x112 in PAV and DZ, i.e., x, _ (Sb /2} ( L, /z,) ~(43b) where L1 is approximated by the high-power transfer function at the equilibrium power level Po rather than at PAV. A more accurate estimate of xl can be obtained from (43a) by iteration if desired. Once a value of xl is obtained from (43) for a given frequency and 6k, we can compute the high power describing function by first evaluating PAV and Dz from (42) and (41), respectively, and then substituting the latter in (40). The following remarks are important at this state: (a) The expression in the curly bracket in (41) is identical up to (6k)2 to M(w,jx!l) in Smets' analysis(8), which is referred to by Smets as the nonlinear frequency dependent gain. (b) In the absence of feedback, H(iw)= Q and the high-power describing function (40) reduces to the single-input zero-power describing function in (26). Furthermore, D(6k,w,PAV) reduces to the high power transfer function L(iO) when Ix1jilO.

-25(c) If we expand D(6k,w,PAV) in powers of 6k after we substitute x1 from (43b) into Dz(X1,x2,W), treat PAV constant (constant mean power operation, see (e) below), and retain terms up to (5k)2 we obtain 2, we~k. D( ) PAv)=L.s 2 i L(+ /+) LZL) -2 ReZi]+O[( Ji (45) This expression is identical to that obtained by a direct Fourier analysis and a perturbation expansion similar to equation (29) in the presence of feedback. Thus (45) represents the first two terms in the exact power series expansion of the high-power describing function. This result is the real proof of the validity of the feedback circuit analysis based on the dual-input describing function. (d) Since the relation between PAV and P represented by (37) does not involve initial conditions or the manner in which the periodic oscillations are started, we can calculate the high-power describing function relative to the equilibrium power level by multiplying D(6k,w,PAV) in (40) by (P AV/P)O The ambiguity associated with the various definitions of the zero-power describing function thus does not arise when feedback is present. (e) The average power level and the externally applied positive reactivity must satisfy equation (42). This relation may be satisfied in a pile-oscillator experiment in different ways. We may for example adjust ko externally in such a way to keep PAV constant

-26and equal to the equilibrium level Po at each frequency and magnitude of the sinusoidal reactivity input. This corresponds to the constantmean-power operation(8) which was mentioned in (c). One may also keep ko constant and allow PAV to vary according to (42). Since the magnitude. and frequency of the self-sustained oscillations are not subject to our control, we must adopt this mode of operation in the nonlinear stability analysis. Stability Analysis The nonlinear stability of a reactor with feedback can be investigated graphically by applying the Nyquist criterion to (40). The Nyquist diagram in this case is obtained by plotting /V( II)= - PAV (Wx) I, ) () D (lx 1i2) (46) in the complex plane as a function of w for various values of Ix 1. When |x1| = 0, the Nyquist plot yields the linear stability of the reactor depending on the encirclement of the -1 point. It is possible that the Nyquist plot may encircle the -1 point at large amplitudes even when the reactor is linearly stable. In this case, the reactor has an "unstable" limit cycle. On the other hand, the reactor may be linearly unstable with a Nyquist plot encircling the -1 point with |xlJ = 0 but not encircling the point when JxlJ exceeds a certain value. In this case, the reactor becomes stable for large amplitudes and a "stable" limit cycle arises. In either case, the amplitude and frequency of the limit cycle is determined by I- PA,(w)Ixi)VH4 w)Dz( IIXIJ) 2)=O (47) /- ~v C/~.~ Ix'12') kt(,'~) DzCu,j2 -o

-27which is the characteristic equation of the nonlinear feedback loop. If 2 this equation does not have a solution for a real w and positive 1|x12, then the system does not admi.t a limit cycle. There can be more than one limit cycle if this equation has several roots. In order to illustrate these ideas, we have examined a reactor model similar to the one employed by Smets(5). The feedback transfer function for the water-moderated system is L/is. A/ (s o.o)(is + O.Sj where the lags correspond to the heating time constant and the bubble formation time constant. In our model the constant A has been chosen so that, in the linear stability analysis, sustained oscillations occur for PAV o P = PC = 1.0. Therefore any Po > 1.0 will produce a linearly unstable system while any Po < 1.0 corresponds to a linearly stable system. For this model A =-10.37 and the critical frequency wc = 3.24 rad/sec. We have presented in Figure 5 the amplitude-dependent Nyquist plot in the vicinity of the -1 point for P = 1.1. The analysis shows that the linearly unstable system has a stable limit cycle as 6k is increased and eventually becomes stable as the Nyquist curve crosses the -1 point. The amplitude of the fundamental of the periodic power oscillations where the limit cycle occurs is obtained approximately using 2|xlJPAvDZIwith Ixll - 0.15 and w- 3.34 radians/sec. as 0.3 PAV. The average power level is found from (42) to be approximately equal to Po in this case. We have also investigated using (40) the magnitude of the describing function at high power as a function of frequency at different amplitudes.

-28The equilibrium power level was chosen to be Po = 0.1 so that the reactor was linearly stable. The internal closed loop amplitude xl was computed using (43b). The results are depicted in Figure 6a where we have also plotted the results obtained using (45) with PAV = Po. For values of 6k up to 0.5 the difference between PAV and Po was found to be less than 2% and insignificant. These curves describe the nonlinear response of a reactor with feedback. When 6k = 0, IDI is identical to the high power transfer function IL(iO)I which displays a resonance at max 1.05 rad/sec. max The height of the resonance peak IDMAXI = D(6k, WMAX' PAV)I as well as oMAX vary with sk as shown in Figure 6b. We observe that the results obtained with the circuit analysis using the dual input zero-power describing function (40) and those from direct Fourier analysis combined with a perturbation expansion differ appreciably for large values of 6k, i.e., 6k z 0.2 and higher. We believe, on the basis of our numerical analysis, that the results obtained with (40) are more accurate for larger values of 6k, particularly near the resonance frequency, because the zero-power describing function (41) is less sensitive to the amplitude variations represented by xl12 Z1 in (41) than the describing function at high power given by (45) in which the amplitude dependence is characterized by Ix112 L1. However the accuracy of the describing function approach in general should be checked by solving the point kinetic equations in the presence of feedback and computing the amplitude of the first harmonic as a function of the reactivity insertion and frequency. Such a calculation has not been attempted in this work because of the large computer requirement. In such numerical applications one must use caution to see that the correction terms due to nonlinearities, i.e. second term in (41) and (45), do not exceed unity because otherwise the expansion in powersof |xl1 may not be convergent.

-29Summary We have presented a new concept in describing function analysis which no longer places the low pass filter restriction on the feedback in a reactor system. The multiple-input describing function concept is a natural the extension from the perturbation analysis used to define4zero-power describing function. In addition, further extension can be made to analysis of any order harmonics arising from feedback in a similar manner to that described here for the second harmonic. For an N harmonic system, this could be represented by a block diagram similar to Figure 2 containing N inputs to the multiple-input describing function arising from N feedback mechanisms of the form H(iOn) where n = 1,... N. This extension has several significant features; first, the analysis of power oscillations in a nonlinear system in the presence of feedback can be made as accurate as needed by including higher harmonic feedback terms. Second, since the circuit analysis requires only a general form for the zero-power describing function, any acceptable form may be used which lends itself to a multiple input description The high power feedback circuit and stability analysis centered around the zero-power describing function illustrate the utility of the concept in the analysis of power oscillations. Decoupling the feedback effects from the nonlinearities of the kinetic equations allows us to treat both aspects of the nuclear reactor in a separate but systematic manner. Although the describing function concept has been used for many years in reactor dynamics, we believe the present work unifies the previous treatments of the subject and clarifies the definition of the describing function with regard to the average power level completely.

-30ACKNOWLEDGEMENT C.M.B. expresses his gratitude to the United States Atomic Energy Commission for support through the Special Fellowship Program. The authors would like to thank Dr. H. B. Smets for his helpful suggestions upon reading the manuscript.

-31References 1. Rumsey, V. H., "The Time Variation of Pile Power Due to Periodic Modulation of the Effective Multiplication Constant," CRP-333, PD-217 (1947). 2. Siegel, R. and Hurwitz, H., Jr., "The Effects of Temperature Coefficients on Reactor Stability and Reactor Transfer Function," KAPL-1138 (1955); Appendix E, by M. Nelkin. 3. Sandmeier, H. A., "The Influence of the Nonlinear Term KEX-n on Power Reactor Stability," Unpublished Memorandum, ANL, September 26, 1957. 4. Akcasu, A. Z., "General Solution of the Reactor Kinetic Equations without Feedback," Nucl. Sci. and Eng., 3, 456-467 (1958). 5. Smets, H. B., "The Describing Function of Nuclear Reactors," IRE Transactions on Nucl. Sci., NS-6, No. 4, 8-12 (1959). 6. Wasserman, A. A., "Contributions to Two Problems in Space-Independent Nuclear Reactor Dynamics," IDO-16755, Phillips Petroleum Co. (1962). 7. Smets, H. B., "Power Oscillations in Nuclear Reactors," Revue A, Belgian Society for Automatic Control, 6, 137 (1964). 8. Smets, H. B., "Derivation of the Describing Function for a Reactor at Power Level," Nucleonik, Band 7, Heft 7, 399 (1965). 9. Akcasu, A. Z. and Shotkin, L. M., "Power Oscillations and the Describing Function in Reactors with Linear Feedback," Nucl. Sci. Eng., 28, 72 (1967). 10. Akcasu, A. Z. and Shotkin, L. M., "Approximate Zero-Power Describing Function Valid at all Frequencies," Nucl. Sci. Eng., 32, 262 (1968). 11. Vaurio,i. K. and Pulkkis, G., "On the Validity of Some Approximate ZeroPower Describing Functions," Nucl. Sci. Eng., 39, 283 (1970). 12. Tan, S., "Sinusoidal Analysis of a Low-Power Nuclear Reactor by the WKB Approach," Nucl. Sci. Eng., 30, 436 (1967). 13. Bethe, H., "Reactor Safety and Oscillator Tests," APDA-117 (1956). 14. Akcasu, A. Z., Lellouche, G., and Shotkin, L. M., Mathematical Methods in Nuclear Reactor Dynamics, Academic Press (in press).7

-32CAPTIONS FOR FIGURES 1. Feedback Circuit Diagram with Single-Input Zero-Power Describing Function. 2. Feedback Circuit Diagram with Dual-Input Zero-Power Describing Function. 3. Amplitude and Phase Characteristics of Single- and Dual-Input ZeroPower Describing Functions. 4. Effect of Second Harmonic on Dual-Input Describing function Characteristics; = 104 sec., 6 = 0.0065, 6k = 0.50 Dollars and W = 3.2 X 10-2 rad/sec. The Second Harmonic was 6k(lx21/lxl1) sin 2wt ( —-), 6k(lx21/1x,1) cos 2wt ( —),and 6k(|x21/1|xl) sin (2wt + 7) (.*. ) for the First Harmonic 6k sin ot. 5. Amplitude-Dependent Nyquist Diagram for a Linearly Unstable Reactor Model 6. (a) Amplitude of the Describing Function at High Power; Eq. (40) -; Eq. (45) --- (b) Variation of the Resonance Peak and Frequency with Reactivity Input.

EXTERNAL REACTIVITY POWER RESPONSE..........- PAvDz (6k,w).. -....._. _.. iwt HARMONICS } H(if) _ HoPA + PAVH(iw)D (6k,)xle + cc Fig.1

2 2iwt: 2 iwt Z X21 + C PeV e + cc ~~~~~~~~11 / \ -- AV 2 SECOND HARMONIC GENERATI9N H2 2iwt EXTERNAL REACTIVITY X.-c\ X~e2" + cc /r, -POWER RESPONSE AVDZ (1'X2 ) 6k iwt k +-e + cc NET REACTIVITY 0 2 iwot Lot 6k0 + xe +cc PAV 1' + y1e +~ cc + y^ ^ + cc} iwt PAVHO + PAVH1lYle + Cc Fig.2

^ ~ ~~ -;\~~~ \ ^ ~=10- SEC 0?? X ^.0065 ~(3~~ "<, cfd~~~~k=0.50 DOLLARS 10 o — TRANSFER FUNCTION <4~~~ N^X "*O~~~''SINGLE INPUT DESCRIBING FUNCTION U --— DUAL INPUT DESCRIBING FUNCTION E-e H z 9 0O f3 En n ~~ ~ H~ " C4~~~ "~ ~< ~' P9 r 9 3 *}~~~~~~\. zi) I 0\ w-30O\ 2~ 0 * "*\. —:l*>^> —__, I 1 l \ i, ^ ^FREQUENCY? (RADIANS/S 4~~~~~~~~~~~~~~~^ ~

^ // /'2~ / / o8 - - 0.0//.05/.10 /.15 /20 /25 /.30 /.35 <^.!

10 9~~~~~~~~ / I ~* ~ /. \~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ III i~~~~~ 0 * IO!~~~~~ / * / / / /! / /r 0 2.^ j*ci~I~xZQ' ^t ~~~~~~~~~~~~~~ / ~~ ~ ~ ~ i~~~~~~~~~~~~~~~~~~~~~~C 7~~~~~~~~~~~ J''C ~~~~~~~~~~~~~~~~~~~~~~~~~ o~~^ / / "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" ~~~~~~~~~~~~/,

| /\ sk=0.3 1' \ 4.0/ - " /^I I 3I^ +, r-S-< 1; q7C1~~~~~~~~~~' -n 4 I \ \. (a) (R^AO/SeC) I Dnx\;, WMAX I 10.0 — i,..../ 5.0-.97 o.0 0-.2 0.4 S(o4o7aR) (b)