Modeling and control of xenon oscillations in thermal neutron reactors

We study axial core oscillations due to xenon poisoning in thermal neutron nuclear reactors with simple 1D models: a linear one-group model, a linear two-group model, and a non-linear model taking the Doppler effect into account. Even though nuclear reactor operators have some 3D computer codes to simulate such phenomena, we think that simple models are useful to identify the sensitive parameters, and study the efficiency of basic control laws. Our results are that, for the one-group model, if we denote the migration area by M and by H the height of the core, the sensitive parameter is H/M. H being fixed, for the 2 groups model, there are still 2 sensitive parameters, the first one being replaced byM1 þM2 whereM1 denotes themigration area for fast neutrons and M2 the migration area for thermal neutrons. We show that the Doppler effect reduces the instability of xenon oscillations in a significant way. Finally, we show that some proportional/integral/ derivative (PID) feedback control law can damp out xenon oscillations in a similar way to the well-known Shimazu control law [Y. Shimazu, Continuous guidance procedure for xenon oscillation control, J. Nucl. Sci. Technol. 32, 1159 (1995)]. The numerical models described in our paper have been applied to PWR.

As reported in [1] xenon oscillations have occurred in Savannah River as early as in 1955.
The fact is that 54 135 Xe is the fission product with the highest capture cross section and that it can be produced either directly or via beta-decay of another fission product which is 54 135 I. Axial core oscillations have been first studied with two zones lumped models, and then by a one group diffusion model coupled to the evolution equations for xenon and iodine.
The first question is whether we should use a time dependent diffusion equation like in [11] or a stationary diffusion equation like in [1] and many other References With the former approach, we benefit from theoretical results on the Hopf bifurcation on non linear systems. With the latter it seems possible to prove theoretical results on the linearized system only.
More precisely, mathematical criteria have been found for the sign of the real part of the eigenvalues of the linearized system (see e.g. [4]).
Physically, there are at least three time scales in the core of a nuclear reactor: the prompt neutrons life time (∼20 ms), the delayed neutrons time scale (∼1 s), and the time scale for xenon oscillations (∼5 h).
A two-group time-dependent model taking the delayed neutrons into account is described in Ponçot's thesis [10].
Practically, in a PWR, criticality of the core is ensured by adjusting either the Boron concentration or by using the control rods. At the time scale of xenon oscillations, we can assume that the core is critical, so that there is no need to consider a time-dependent equation for the neutron flux.
Rather, at each time step, the xenon concentration being given, it is appropriate to solve an eigenvalue problem.
The second question is whether we should use a onegroup model or a two-group model (one group for fast neutrons, one group for thermal neutrons) like the one considered in [10].
Indeed, since the xenon capture cross section is significant for thermal neutrons only, one might think that it is more appropriate to introduce a two-group model.
In the present paper, we show that it is possible to derive a one-group model with the same behavior as our two-group model provided that, in the one-group model, the diffusion equation holds for the thermal neutron flux only and that the migration area is selected in an adequate way.
We show that for both one-group and two-group models, the height of the core and the magnitude of the neutron flux are sensitive parameters.
In [12], we have shown that, as far as xenon oscillations are concerned, introducing a reflector is equivalent to increasing the height of the core.
After showing that the one-group model is possible, we add a nonlinear term in the diffusion equation to take the Doppler effect into account: like in [1], we find that it helps damping out the oscillations.
Finally, we address the question, which is of primary importance, of how to better control the oscillations.
We check that it is efficient as expected, but we show that a PID control strategy can also be applied.
Our models are presented in a mathematical way, and we use adequate rescaling methods both to reduce the number of parameters and to facilitate their application (which is outside the scope of this work) to other thermal neutron reactors like CANDU or HTR. The results we show below are related to PWR.

The 1 group diffusion model
According to the neutron transport theory and the diffusion theory, the one group diffusion equation for a homogeneous reactor can be written, in the uniaxial case, as: where H is the height of the core. F is the neutron flux of the core in n/cm 2 /s and M 2 denotes the migration area in cm 2 . This differential equation should be completed by some boundary conditions (BC). In the following we shall choose the Dirichlet BC: We refer the reader to [3] for the Fourier BC.
In a one group model, the migration area is arbitrarily evaluated by taking into account both fast and thermal neutrons.
In the present section which is dedicated to the one-group model, the solution F to equation (2) will be assumed to be the thermal neutron flux in the core.
We remind the reader that (1) and (2) is an eigenvalue problem and then that F is defined up to a multiplicative factor.

Introduction of iodine and xenon
We have the following evolution equations: where: -Y (resp. X) is the concentration of 53 135 I (resp. 54 135 Xe) in the fuel evaluated in at/cm 3 of fuel; s is the microscopic capture cross section for xenon in the thermal domain (evaluated in barns); -S f is the macroscopic fission cross section of the fuel in the thermal domain; g the iodine fission yield (about equal to 0.064 for UO2 fuel); g 0 the xenon fission yield (about equal to 0.001 for UO2 fuel); l (resp. m) the time constant for the b À decay of 53 135 I (resp. 54 135 Xe).
For the sake of simplicity, in the following, we shall neglect xenon generation from fission reaction directly, which is possible since g 0 ≪ g.
Note that in (3) and (4) S f is evaluated in cm À1 and F in n/cm 2 /s. The core is assumed to be homogenized so that, in each cm 3 of core, there are v cm 3 of fuel and (1 À v) cm 3 of moderator. In a PWR we have v ∼ 1/3 so that the moderation ratio (1 À v)/v ∼ 2.
Since there are vX at/cm 3 of 54 135 Xe in the core, the xenon macroscopic cross section is equal to vsX Once taking into account the neutron capture by 54 135 Xe equation (1) has to be replaced by where a = vs/ S and S is the absorption cross section of the core.

Rescaling of the one group model
Following [16]  We also replace t by t = lt, and we let F Ã = l/s, a = m/l We get the simpler system X Ã has the following interpretation: it is the equilibrium value for the concentration of 135 Xe atoms when the neutron flux is infinitely large. Equations (6)-(8) completed with Dirichlet BC is the one-group system which we shall solve.
We could have also changed the space variable z → Z/M to prove that the main parameter is H/M.
Typical values of the parameters to be selected for PWR reactors are shown in Table 1.

Space approximation
We introduce a discretization step Dz = H / (N+1) where N is the number of discretization points. We now denote by ' the column-vector of values or the rescaled thermal neutron flux at the discretization points i Dz, 1 i N. (Our numerical results will be obtained with N = 39, but this is not a sensitive parameter. ) We introduce the tridiagonal matrix A such that where we assume that ' 0 = 0 and ' N+1 = 0 which corresponds to the Dirichlet boundary conditions. We shall make the following convention : (x') i = x i ' i , which means that, depending on the context x may denote either a vector in ℝ N or the N Â N diagonal matrix containing the components of x on its main diagonal.
Let {y ∞ , x ∞, ' ∞ } be a steady solution of (6)-(8) satisfying We note that the solution {y ∞ , x ∞, ' ∞ } of (9) and (10) cannot be unique: in fact x ∞ being given, (10) is an eigenvalue problem (which means that the core is critical) so that ' ∞ is defined up to a multiplicative constant and we fix this constant so that the core power is given.
In the case where our steady solution {y ∞ , x ∞, ' ∞ } is not stable, if we start from {y a , x a } close to {y ∞ , x ∞ } we shall activate an oscillation. This means that the reactor core is unstable.

Time discretization: semi implicit scheme
We denote the time step by d. We let t n = nd.
We use the following semi implicit scheme: First step, solve Which produces both a value for k n ∞ (the smallest eigenvalue) and a (non-unique) value for ' n .
Then we replace ' n by ' n times ' = averageð' n Þ, where ' is a prescribed value for the average of the (normalized) flux.

Results with the one group model
When the average flux is sufficiently high the flux oscillations are first growing and then stabilize.
A typical result is given in Figure 1 for which ' ¼ 3, which is a rather high value since 3F * = 4.38 Â 10 13 n/cm 2 /s.
To illustrate the evolution of the oscillations we shall use the Axial Offset AO p defined by AO p = (P t À P b )/ (P t + P b ) where P t (resp. P b ) is the average power (or neutron thermal flux) in the upper (resp. lower) half part of the core. Figure 2 shows the evolution of AO p with the reduced time t.
When the oscillations are growing they reach a limit cycle, and this is because we ensure criticality at all times. Now from the operator's point of view these oscillations are instabilities.
On a standard PWR we know that the migration area is more like 64 cm 2 . Note that the period of oscillations depends both on the migration area and on ': it can range between 28 and 36 h, which is consistent with the measurements made in actual PWR [17] (∼30 h for a 900 MWe [5]). Table 1. Typical numerical data for our one-group model (see [11] p. 382).
In Figure 3 we show the influence of the height of the core on the oscillations.
When the reactor core is lower than 250 cm, no oscillation is activated. Once the height is more than 300 cm, a growing oscillation can be produced. Moreover, the higher the height is, the bigger AO p is.

The 2-group model
There are many neutrons with different energy in a reactor core. It is usually considered that the fast neutrons possess an energy more than 0.625 eV, and on the contrary, that the thermal neutrons have an energy less than 0.625 eV.  Applying the data in [11], the thermal neutron flux is equal to 3.7 Â 10 13 n/cm 2 /s in average in the fuel pins of a PWR. In addition, the fast neutron flux is almost 7 times as high as the thermal neutron flux. It's then approximately equal to 2.7 Â 10 14 n/cm 2 /s in average in the fuel. Even though the highest flux in a PWR core corresponds to the fast domain, more than 80% of the fissions happen in the thermal domain. For a 1300 MWe PWR, there are about 1.1 Â 10 13 fissions/cm 3 /s in the fuel at nominal power (3800 MWth).
Let C (resp. F) denote the fast (resp. thermal) neutron flux, we shall use the following system of coupled diffusion equations Like in Section 2, vsX is the xenon macroscopic capture cross section : it should be added to S 2 .
In the same way, if S f1 (resp. S f2 ) denotes the macroscopic fission cross section of the fuel in the fast (resp. thermal) neutron domain, v S f1 C (resp. v S f2 F) denotes the number of fast (resp. thermal) fission per cm 3 per s in the core.
Note that the xenon capture cross section is negligible for fast neutrons. Now, S 1 C is the number of fast neutrons which disappear per cm 3 per s. Some of them are physically absorbed but a fraction p of them reappear, after slowing down, in the thermal neutron group.
The factor p appears to be the resonance escape probability from Enrico Fermi's four factors formula.
For the sake of simplicity, we shall assume that the diffusion coefficients do not depend on the space variable z.
Even though (11) and (12) are stationary diffusion equations, we may have X = X(z, t).
The xenon evolution is coupled to the iodine evolution through Such an evolution is rather slow at the scale of one hour, so there is no need to consider the neutron evolution at the timescale of delayed neutrons. The 2 groups model will be complemented by some boundary conditions for the neutron flux.
We shall mainly consider the Dirichlet boundary conditions but it is also possible to consider the Fourier boundary conditions. Remark: The number of neutrons produced by fast fission of 92 235 U is actually 2.8, but, for the sake of simplicity, we choose n = 2.4 both for 92 238 U and 92 235 U. However, it would be easy to change the term n v S f1 C in (11) by n 1 v S f1 C where n 1 is some appropriate average of 2.4 and 2.8.

Rescaling of the 2-group model
Like in Section 2 the system (11)- (14) can be simplified.
At first, let Then, we rewrite the variables And we obtain To simplify the above equations, we choose X * = Y * and We can simplify the system one more time and we have

Space approximation
We introduce the tridiagonal matrix A as in Section 2.
Let {y ∞, x ∞, ' ∞, c ∞ } be a steady solution of (16)-(19) satisfying The parameter k ∞,1 being given, let B = B(x, k ∞, 2 ) denote the 2N Â 2N matrix defined by The solution {y ∞, x ∞, ' ∞, c ∞ } of (20), (21) and (22) cannot be unique: is not reduced to 0 0 , which means that the core is critical, in which case (NB. Applying the Perron-Frobenius theorem to matrix (aI+B) À1 we see that dim(Ker B) = 1). After elimination of c ∞ thanks to (22) we find that k ∞,2 is the smallest eigenvalue of And this is the way we compute k ∞,2 . As in the one-group model, if our steady solution is not stable, even if we start from {y a, x a } close to {y ∞, x ∞ } we may activate an oscillation. Like above this means that the reactor is unstable.
To simulate such a perturbation, the xenon concentration will be increased by 10% in the half lower part of the core. In the half upper part, it will be decreased by 10%.

Time discretization: semi implicit scheme
We proceed as in Section 2 and introduce the following semi implicit scheme: First step, solve y n À d zc nÀ1 þ ' nÀ1 À y n À Á À y nÀ1 ¼ 0; x n À d y n À ða þ ' nÀ1 Þ À Á x n À x nÀ1 ¼ 0: which produces both a value for k n ∞ and a (non-unique) value for ' n .
Then we replace ' n by ' n times '=average ð' n Þ, where ' is a prescribed value for the average of the (normalized) flux.
Typical values of the parameters to be selected for PWR reactors are shown in Table 2.

Numerical results
The power of the reactor will be maintained constant which means the average neutron flux is fixed. ' ave ¼ ' ¼ 3 so that F Ã,ave = 3F * = 4.38 Â 10 13 n/cm 2 /s. (This is somewhat too high for a PWR reactor).
In Figure 4 we give the thermal flux profile as a function of time. This should be compared to Figure 1. We notice that when M 2 ¼ M 2 1 þ M 2 2 (which is the case here) the results given by the one group and the two-group models look pretty much the same.
(We have observed in our results that ' ≅ c, but this is no surprise when we consider (17): this results from a Ã x ≪ 1 and from M 2 ≪H.) In Figure 5 we compare one group and two-group results. We plot AO p as a function of the reduced time t. We see that the behavior of one and two group model is very much the same. This is because we have selected the same (Dirichlet) boundary for both fluxes.
We check again that the migration area has a significant influence on the time-period of oscillations.

The non-linear model (power feedback effect)
Most water moderated reactors have a negative moderator temperature coefficient, which means that when the moderator temperature increases, the reactivity of the core decreases (see e.g. [6]).
For lightly enriched Uranium fuel there is also a negative fuel temperature coefficient which is due to the  Doppler effect which leads to a broadening of isotope capture resonances of U238.
Controlling the chain reaction is obviously of primary importance for nuclear safety.
Such negative temperature effects are then very useful.
Like in [3] we shall use the following nonlinear model: b is the prompt power feedback coefficient. It takes the Doppler effect into account and is determined by the temperature difference between the actual reactor shutdown and full power operation, independent of other parameters (time, boron, flux etc.). This model accurately describes what happens with the fuel (Doppler) temperature effect.
To take into account the moderator temperature effect, we should compute the moderator temperature locally, which would mean coupling neutron diffusion equations and thermal hydraulics.
To keep simple models, we shall not do that and we refer the reader to [18] for a simplified approach.

Space approximation
Proceeding like in Section 3, we obtain the following system Algorithm : we notice that (26)-(28) is a particular case of the following differential system: where which we should complement with initial conditions: When we do not take the temperature effect into account (i − e when b = 0) the RHS of (28) is just equal to k ∞ ' which means that (28) is an eigenvalue problem. Therefore, it has a solution ' ≥ 0 only if k ∞ is equal to the smallest eigenvalue of matrix In the nonlinear case (b > 0) we can hope (and that's what we see in our simulations) that if k ∞ is larger than this smallest eigenvalue, (28) will have a solution ', and that this solution is larger if k ∞ increases.

Time discretization: fully implicit scheme
We have used successfully the following fully implicit scheme.
x n À x nÀ1 À Á =d ¼ F y n ; x n ; ' n ð Þ ; ð33Þ This seems quite obvious, however we note that in a nuclear reactor, when the power is given, it means the average value of the flux is equal to ', the reactivity is adjusted to zero not by the Doppler effect, but by the boron concentration, which means that the control system selects k n ∞ as the smallest eigenvalue of matrix.
where ' Ã = m n ' n and m n ¼ '=averageð' n Þ. We leave to the reader the details of implementation of Newton's method to solve the non-linear system (32)-(34) of 3N equations for 3N unknowns.

Time discretization: semi implicit scheme
Step 1: Like in the fully implicit scheme, we evaluate the smallest eigenvalue k n ∞ of the matrix A defined in (35).
Step 2: Knowing y n , x n we compute ' n by solving the N Â N non-linear system G y n ; x n ; ' n ð Þ¼0: (For this, we use Newton's method).
Step 3: We compute y n+1 , x n+1 with: Remarks: the algorithm we use in the linear case (b = 0) is basically the same, except that we avoid step 2; the results we obtained with both algorithms are almost the same.

How to choose b?
In PWR 1300 MWe, the rated reactor power corresponds to ' ∼ 3.
The "power defect" in the EDF terminology is of the order of 2000 pcm. But we're only going to take the Doppler effect which is of the order of 1000 pcm. (In fact, differing from the Doppler effect, the moderator effect is not a local effect because of the water circulation in the core, and therefore an increase in its temperature certainly reduces the reactivity of the core, but over much of its height).
This means that we consider b F moy ≅ 0.01. Even with the Doppler effect, some oscillations can be observed in Figure 6. From Figure 7, we see that, with the Doppler effect, the amplitude of the oscillation decreases.

Shimazu's representation
Like Shimazu [13] we introduce another axial offset for iodine And a different one for xenon NB. Note that at constant power the normalized xenon level should be < 1 so that we have both X t < 1 and X b > 1, which proves that Shimazu's ellipses are a special representation of the results where we plot AO Y À AO X with respect to AO p À AO X . Here is what we get with M 2 =100 cm 2 , ' ¼ 3 and a Ã = 0.0386. The initial point starts from the horizontal axis. It turns in the trigonometric sense. Note that the curve looks more like a divergent spiral than an ellipse in case there is no Doppler effect in Figure 8. However, as can be seen in Figure 9 it converges to a limit ellipse if we follow it on a sufficient period of time. When we take the Doppler effect into account the spiral is convergent to origin which means that the oscillations are killed.
Comparing with the diverging curve of the model without the Doppler effect, we can draw a conclusion: the Doppler effect effectively weakens the oscillation. Actually, the Doppler effect plays an important role in the inherent safety of reactor.
However, even if the oscillations are more limited than without Doppler effect, they still need to be controlled, and that is what we study in the next section.

Shimazu's control method
Shimazu has been the first to use the axial offsets AO p , AO x and AO y to control the xenon instability.  As has been seen previously, in some circumstances, the power may oscillate from the lower part to the upper part of the core. Now, boron concentration is not the only way to control reactivity of the core, we also have the control rods. A natural idea is to insert the control rods when the neutron flux is high in the upper part of the core, and at the same time, decrease the boron concentration to avoid shutdown. Now we should not wait for the time when the power is maximum in the upper part of the core to insert the control rods. Rather we should anticipate. Shimazu's representation gives a way to do so.
When we drop the control rods, the representative point will move to the left because AO p will decrease, however, if we raise the control bars, it will move right.
Note that we use lightly absorbing control rods, which means that when they are fully inserted, they reduce the reactivity of the core by 1100 pcm.
As in [16], we select 10 cm for a step.
Here are the criterions: -If AO y À AO x > 0 and AO y ÀAO x AO p ÀAO x < tanð' R Þ we raise the control rods for a step.
To improve the result, if AO y À AO x > 0.2 we raise the control rods for two steps.
-If AO y À AO x < 0 and AO y ÀAO x AO p ÀAO x < tanð' R Þ we drop the control rods for a step.
To improve the result, if AO y À AO x < À0.2 we drop the control rods for two steps.
To calculate ' R , we should use the following formula: phase at the initial time of AO y and AO p respectively. AM p , AM y , AM x : amplitude of AO p , AO y and AO x respectively.
g: stability index of the xenon oscillation v: angular velocity of oscillation.
At the initial moment, it has an oscillation because there is the disturbance and we used an application named Origin (https://www.originlab.com/) to get the sinusoidal functions of AOP and AOI, then we know ' i et ' p and finally ' R .
We can see in Figure 9 that Shimazu's method eliminates the xenon oscillation quite well. (The computation is carried out with M 2 = 100 cm 2 ).

PID method
We wanted to try a more standard PID regulator.
PID is a proportional-integral-derivative controller which is widely used in industrial control systems.
To regulate a process it has been noticed that the control should not depend only on the gap between the observation and the target. Here the control is the rod's insertion U(t), the observation is the axial offset AO p , and the target is zero. In the following, we shall say that the error is AO p − 0 = AO p .
In the real system we can only compute the control U(t) at some discrete times.
Following the PID principle, to define U(t) we shall combine a proportional term, a derivative term and an integral term: After time discretization, we can obtain a new formula where: U: variable that can be controlled err (or e): error (that is AO p ) T: control period (in the following T = 0.02 in reduced units) K p : proportional coefficient. In the following, we let K i ¼ K p T T i and K d ¼ Since U denotes the control rod's position and e the axial offset AO p , we can describe the position of control bars by this equation: x 1 ðkÞ ¼ 10 30 þ K p :AO p þ K i S k n¼0 eðnÞ þ K d : eðkÞ À eðk À 1Þ ½ : Note that we have selected x 1 ¼ 300 as the initial position of our control bars. Figure 10 shows the results obtained with PID when M2 = 100 cm 2 ; Kp = À75, T i = 50, T d = 0.01.
When we use PID regulator to control the system, we found that T d is a sensitive parameter. Let's just consider the differential item: K d · [e(k) À e(k -1)].
We let T i = 0 and we neglect the proportional item, then we can describe the position of control bars by this equation: which means, practically, that when AO p increases, the control bars should be inserted. We note that K d is the more sensitive parameter. Figure 11 corresponds to K d = 130.
We can conclude that an adequate PID regulator can also control the system quite well.

Conclusion
Even though Shimazu's control law has been designed by using a simple 2 nodes lumped model of the core, it works well also when applied to a more detailed one group finite element model of the core. This is well known by nuclear reactor operators like EDF who developed the COCCINELLE code [10,17].  In our paper, we study in Section 3 a two-group model and show that it gives basically the same trend as the one group model, but the choice of physical parameters is more obvious with the two group model.
Both models show that xenon oscillations are significant on large PWR cores.
The oscillations are reduced when some nonlinear effects like the Doppler effect are taken into account (see Sect. 4).
However, the xenon oscillations still need to be controlled.
In Section 5, we show that Shimazu's method is efficient, but more standard methods like the widely used PID regulator are also efficient.
We have applied our numerical models to PWR. They could be easily applied to BWR since we adjust criticality at each time by means of a homogeneous type of control of reactivity (Boron concentration adjustment for a PWR, recirculation pumps for a BWR). For other reactors like HTR where criticality is adjusted by means of control bars only, our python programs should be adapted, but this is feasible.