EXTENSION OF XENON OSCILLATIONS SAFETY MARGINS USING WEAKLY NONLINEAR STABILITY ANALYSIS

Weakly nonlinear stability analysis is applied to study xenon oscillations in nuclear reactors using the approach of multiple time-scales method. This approach allows to characterize the dynamics of the system beyond the Hopf instability point. It provides important insight on the characteristics of the oscillations, namely if they diverge with time, or converge into a bounded limit cycle. Detailed derivation of the amplitude equation is presented. This equation is used to identify parameter ranges of bounded periodic oscillations, which may be allowed for safe operation. The influence of neutron generation time and the power feedback coefficients on the amplitude of limit cycles, as well as on convergence times, is discussed. The described method may be used to extend the safety margins required to prevent xenon unstable oscillations in reactor cores.


INTRODUCTION
The design of new and innovative nuclear reactors, as well as the improvement of present-day reactors, must include a detailed study of the dynamics and stability of such systems. In the presence of internal and/or external feedback processes, such as the formation of absorbing fission products, small deviation from nominal steady-state operation may increase in a diverging manner, driving the reactor core out of its safe operational zone. Such phenomena have been observed previously in several types of reactors [1,2] Various instability mechanisms have been identified and studied in the past decades [3][4][5]. Among these, xenon instabilities are of major importance. The formation and depletion of xenon-135, a strongly-absorbing fission product, have shown to yield space and time dependent flux oscillations. Such oscillations may either be bounded in a finite amplitude, or they may increase with time, yielding core instability [6].
In present-day reactors, control procedures are well established to prevent xenon unstable oscillations [7,8]. In addition, safety margins are typically applied, requiring that the operational conditions are considerably distant from the unstable region. Such margins are usually derived using complex reactor physics code platforms. While these tools are able to produce highly accurate results, they lack the ability to provide analytical insights into the governing feedback mechanisms in the core, which are of high importance when new reactor designs are considered.
In this work we use weakly nonlinear stability analysis in order to study the dynamics of xenoninduced oscillations in time. Stability conditions are derived analytically, and are used to distinguish between oscillatory instabilities of diverging amplitudes and instabilities that result in bounded periodic oscillations or limit cycles. The existence of such limit cycles can pave the way in future studies to extending operational limits and safety margins which are required for the suppression of xenon oscillations and to ensure reactor stability.

WEAKLY NONLINEAR STABILITY ANALYSIS
Based on previous studies [9,10], we assume a one-speed (thermal) neutron flux (φ) equation in an infinite homogeneous medium, coupled to number density equations of the isotopes iodine-135 (I) and xenon-135 (X ): where x and t are the space and time coordinates, v is the average velocity of thermal neutrons, ν is the average neutrons yield per fission event, σ X is the microscopic cross section of xenon-135 for neutron absorption, γ I , γ X , λ I , λ X are the direct yields and β-decay rates of iodine-135 and xenon-135, correspondingly, Λ is the neutron generation time of a reactor,ρ ∞ is the hot reactivity of an infinite reactor (excluding xenon-135 absorption) and α f is the power feedback coefficient of a reactor. Further description of the model, including boundary conditions and parameter values, is detailed elsewhere [9,10]. Notice that here neutron leakage is assumed to be negligible, hence the model corresponds to either an infinite or a fully reflective system.
Linear stability analysis, which describes the system response to small perturbations from equilibrium, was previously presented [9,10], and several uniform Hopf bifurcations were identified. However, beyond the instability thresholds, the perturbations grow exponentially fast and nonlinear terms cannot be neglected. Depending on these terms, the dynamics can either diverge or converge to a bounded limit cycle. To study this behaviour, a weak nonlinear analysis in the vicinity of the uniform Hopf bifurcation was carried out. This provides additional significant information, namely, the conditions under which the oscillations that develop beyond the instability thresholds are bounded and what parameters control their amplitudes' magnitude.
We consider small deviations from the neutral stability line by introducing a small parameter: whereρ ∞,0 is a specific value ofρ ∞ which corresponds to values of flux and power feedback on the neutral stability line [9,10]. Under such deviations, multiple time scales perturbation theory may be applied [11,12].
Solutions to Eqs. (1a)-(1c) can now be approximated as power series of The steady-state homogeneous solution is represented by the first term on the RHS of Eq. (3). The second term represents the leading order in the approximation. As we are in the vicinity of a uniform Hopf bifurcation, it may be written as a product of a slow varying amplitude and a fast oscillating Hopf frequency: Note that the amplitudeÃ is assumed to depend only on "slow" space and time coordinates: hence the derivatives ofÃ with respect to the "fast" space and time variables are of order √ λ and λ, respectively.
In the followings, the derivation of an amplitude equation, which describes the dynamics of the amplitudeÃ, is presented. First we substitute the expansion (3) in Eqs. (1a)-(1c), transform the time and space derivatives as in [11,12], and equate terms of equal powers of √ λ up to third order (the zeroth order yields the steady-state solution, as expected). At order √ λ we obtain: where the operator L is given by: Substituting Eq. (4) in Eq. (7) yields α and β in terms of the model's parameters, and also an expression for the Hopf frequency. At order λ we obtain, where C i and D i (i = 1, 2, 3) depend on the amplitudeÃ. The solution of Eq. (9) reads: where the quantities γ 1 , γ 2 , δ 1 , δ 2 , B 1 , B 2 depend also on the amplitudeÃ.
Finally, at order λ 3/2 we obtain: As with the previous order, the quantities E i and F i (i = 1, 2, 3) depend on the amplitude and on the parameters of the model. This equation contains secular terms of the form e ±iω H t , and according to Fredholm's theorem [11,13], a solvability condition should be imposed, namely, that the right hand side of Eq. (11) is orthogonal to the null space (kernel) of the adjoint operator of L, denoted as L † [11,12]. This kernel is defined by eigenvectors h that satisfy: We assume the form: and substitute it in Eq. (12) to obtain both and ζ. Now the solvability condition can be imposed: and using Eqs. (11) and (13), the amplitude equation is found: This equation has the form of the complex Ginzburg-Landau equation, the normal-form equation for a Hopf bifurcation [11,12], where A = √ λÃ. In Eq. (15), the quantities λη and λθ represent, respectively, the linear growth rate and the linear frequency correction to ω H . On the other hand, κ and ξ are the first order nonlinear corrections. These parameters determine the corrections to both the growth rate of the oscillatory mode and its frequency. In particular, a positive value of κ, leads to limit cycles, i.e. bounded periodic oscillations of finite amplitude.

SOLUTIONS OF THE AMPLITUDE EQUATION
It is possible to obtain analytical solutions of the amplitude equation (15), by expressing it in terms of the polar variables, ρ = |A| and φ = arg (A) (i.e. A = ρe iϕ ), obtaining: The equation for ρ does not depend on ϕ and can be readily integrated. Substituting the solution for ρ in the equation for ϕ and integrating, we find the solution: where ρ 0 and ϕ 0 are the initial conditions for ρ and ϕ, respectively.
The asymptotic behaviour of ρ can be analyzed using the limit t → ∞ in Eq. (17a)), yielding: The parameters η and κ must have the same sign in order for the solution to exist. In addition, η is always positive in linearly unstable regions of the steady-state solution. Hence, κ must also be positive for the solution (18) to exist.
In Fig. 1, the temporal dynamics of the neutron flux is presented in respect to the steady-state flux value. The flux was calculated using numerical integration of Eqs. It can be seen that the numerical solution indeed converges with time to the analytical asymptotic solution. As depicted in the figure, for shorter generation times, the amplitude of the limit cycle increases. In addition, for weaker power feedback, the amplitude also increases. However, these parameters yield only slight modifications to the time required for the convergence to the limit cycle behaviour. In Fig. 2, the phase plane spanned by φ and I, calculated using a numerical integration of the model's equations, is presented. As depicted in the figure, a small deviation from the steady-state solution results in convergence to a constant-amplitude oscillatory solution represented as a closed trajectory around the original steady-state solution, i.e. a limit cycle.
Figure 2: Convergence to a limit cycle in the phase plane spanned by the neutron flux and the iodine number density. The limit cycle, illustrated as a black curve, represents stable oscillatory reactor dynamics. The black dot corresponds to the steady-state solution (from [10]).
The existence of limit cycles suggests that flux dynamics in regions of linear instability may in fact be bounded, due to nonlinear feedback mechanisms between neutron flux and xenon. In terms of reactor operation, such regions may not necessarily be excluded for safe reactor operation, provided that the amplitude of oscillations (ρ) is small with regard to the steady-state solution of the neutron flux.

CONCLUSIONS
Nonlinear stability analysis was applied in order to study xenon oscillations in a nuclear reactor. Using multiple-time scales analysis, an equation for the amplitude of oscillations in the vicinity of the newly discovered Hopf bifurcations was derived, and was used to study the nature of this bifurcation as a function of operational parameters, i.e. whether it is a super-critical bifurcation leading to small-amplitude oscillations, or a sub-critical bifurcation resulting in unbounded oscillations or large-amplitude oscillations.
Numerical integration of the model's equations in the vicinity of the neutral stability line shows the influence of different generation times and power feedback values on the amplitude of the limit cycles. The nonlinear analysis demonstrates the utility of the amplitude equation in identifying flux and power feedback values for which safe, small-amplitude oscillations may be achieved.
Xenon oscillations are typically avoided in current reactors by either using advanced control systems or by limiting the reactor operation to stable flux and power feedback ranges. The analysis presented here presents deviations from the neutral stability line that may be acceptable as they induce oscillations which are weak and stable. Once derived from a more general description, using similar methods to those presented here, this approach may therefore be used to extend the operational limits and safety margins which are required for the suppression of xenon oscillations and to ensure reactor stability.