Issue 
EPJ Nuclear Sci. Technol.
Volume 6, 2020



Article Number  48  
Number of page(s)  13  
DOI  https://doi.org/10.1051/epjn/2020009  
Published online  28 May 2020 
https://doi.org/10.1051/epjn/2020009
Regular Article
Modeling and control of xenon oscillations in thermal neutron reactors
^{1}
CEA/INSTN, 91191 Gif sur Yvette, France
^{2}
Institut FrancoChinois de l'Energie Nucléaire, Sun Yatsen University, 519082 Zhuhai Campus, Guangdong Province, P.R. China
^{*} email: bertrand.mercier@cea.fr
Received:
18
November
2019
Received in final form:
3
March
2020
Accepted:
12
March
2020
Published online: 28 May 2020
We study axial core oscillations due to xenon poisoning in thermal neutron nuclear reactors with simple 1D models: a linear onegroup model, a linear twogroup model, and a nonlinear 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 onegroup model, if we denote the migration area by M ^{2} 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 by M_{1}^{2}+M_{2}^{2} where M_{1}^{2} denotes the migration area for fast neutrons and M_{2}^{2} 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 wellknown 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.
© B. Mercier et al., published by EDP Sciences, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
A lot of publications can be found in the literature on the problem of xenon oscillations in nuclear reactors [1–10].
As reported in [1] xenon oscillations have occurred in Savannah River as early as in 1955.
The fact is that is the fission product with the highest capture cross section and that it can be produced either directly or via betadecay of another fission product which is .
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 μs), the delayed neutrons time scale (∼1 s), and the time scale for xenon oscillations (∼5 h).
A twogroup timedependent 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 timedependent 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 twogroup 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 twogroup model.
In the present paper, we show that it is possible to derive a onegroup model with the same behavior as our twogroup model provided that, in the onegroup 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 onegroup and twogroup 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 onegroup 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 test the wellknown Shimazu control law [13–15].
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.
2 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: (1)
where H is the height of the core. Φ 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:(2)
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 onegroup model, the solution Φ 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 Φ is defined up to a multiplicative factor.
2.1 Introduction of iodine and xenon
We have the following evolution equations:
where:

Y (resp. X) is the concentration of (resp. ) in the fuel evaluated in at/cm^{3} of fuel;

σ is the microscopic capture cross section for xenon in the thermal domain (evaluated in barns);

Σ_{f} is the macroscopic fission cross section of the fuel in the thermal domain;

γ the iodine fission yield (about equal to 0.064 for UO2 fuel);

γ′ the xenon fission yield (about equal to 0.001 for UO2 fuel);

λ (resp. µ) the time constant for the β − decay of (resp. ).
For the sake of simplicity, in the following, we shall neglect xenon generation from fission reaction directly, which is possible since γ ′ ≪ γ.
Note that in (3) and (4) Σ_{f} is evaluated in cm^{−1} and Φ in n/cm^{2}/s.
The core is assumed to be homogenized so that, in each cm^{3} of core, there are ω cm^{3} of fuel and (1 − ω) cm^{3} of moderator. In a PWR we have ω ∼ 1/3 so that the moderation ratio (1 − ω)/ω ∼ 2.
Since there are ωX at/cm^{3} of in the core, the xenon macroscopic cross section is equal to ωσX Once taking into account the neutron capture by equation (1) has to be replaced by(5)
where α = ωσ/ Σ and Σ is the absorption cross section of the core.
In view of (3) and (4) we shall have Y =Y(z, t), X = X(z, t) and then Φ = Φ(z, t) from (5).
2.2 Rescaling of the one group model
Following [16] we let
We also replace t by τ = λt, and we let Φ_{*} = λ/σ, a = μ/λ and . 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 onegroup 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.
2.3 Space approximation
We introduce a discretization step Δz = H / (N+1) where N is the number of discretization points. We now denote by ϕ the columnvector of values or the rescaled thermal neutron flux at the discretization points i Δz, 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(9) (10)
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.
2.4 Time discretization: semi implicit scheme
We denote the time step by δ. We let τ_{n} = nδ.
We use the following semi implicit scheme:
Which produces both a value for (the smallest eigenvalue) and a (nonunique) value for ϕ ^{ n }.
Then we replace ϕ ^{ n } by , where is a prescribed value for the average of the (normalized) flux.
2.5 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 , which is a rather high value since 3Φ_{*} = 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 τ.
We see that the oscillations are

damped out when M ^{2} = 100 cm^{2} and ;

stable when M ^{2} = 64 cm^{2} and ;

growing when M ^{2} = 100 cm^{2} and or M ^{2} = 64 cm^{2} and .
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]).
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.
Fig. 1 The oscillation of neutron flux in 1 group model with M ^{2} = 64 cm^{2}. 
Fig. 2 Evolution of the axial offset AO _{p} for various values of and M ^{2} = M2. 
Fig. 3 Evolution of AO_{p} for various values of height. 
3 The 2group 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 Ψ (resp. Φ) denote the fast (resp. thermal) neutron flux, we shall use the following system of coupled diffusion equations(11) (12)
Like in Section 2, ωσX is the xenon macroscopic capture cross section : it should be added to Σ_{2}.
In the same way, if Σ_{f} _{1} (resp. Σ_{f} _{2}) denotes the macroscopic fission cross section of the fuel in the fast (resp. thermal) neutron domain, ω Σ _{f1}Ψ (resp. ω Σ _{f2}Φ) 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, Σ_{1}Ψ 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(13) (14)
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(15)
but it is also possible to consider the Fourier boundary conditions.
Remark: The number of neutrons produced by fast fission of is actually 2.8, but, for the sake of simplicity, we choose ν = 2.4 both for and . However, it would be easy to change the term ν ω Σ _{f1}Ψ in (11) by ν _{1} ω Σ _{f1}Ψ where ν _{1} is some appropriate average of 2.4 and 2.8.
3.1 Rescaling of the 2group model
Like in Section 2 the system (11)–(14) can be simplified.
At first, let
Then, we rewrite the variables Φ = Φ_{*} ϕ, Ψ = Ψ_{*} ψ, X = X _{*} x, Y = Y _{*} y, τ = λt.
To simplify the above equations, we choose X _{*} = Y _{*} and
If we choose
.
We can simplify the system one more time and we have(16) (17) (18) (19)
3.2 Space approximation
We introduce the tridiagonal matrix A as in Section 2.
Let {y _{∞,} x _{∞,} ϕ _{∞,} ψ _{∞}} be a steady solution of (16)–(19) satisfying(20) (21) (22)
The parameter k _{∞,1} being given, let B = B(x, k _{∞, 2}) denote the 2N × 2N matrix defined by
The solution {y _{∞,} x _{∞,} ϕ _{∞,} ψ_{∞}} of (20), (21) and (22) cannot be unique: does not vanish only if Ker B(x _{∞,} k _{∞,2}) is not reduced to , which means that the core is critical, in which case
(NB. Applying the PerronFrobenius theorem to matrix (aI+B)^{−1} we see that dim(Ker B) = 1).
After elimination of ψ_{∞} 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 onegroup 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%.
3.3 Time discretization: semi implicit scheme
We proceed as in Section 2 and introduce the following semi implicit scheme:
which produces both a value for and a (nonunique) value for ϕ ^{ n }.
Then we replace ϕ^{n} by , 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.
Typical numerical data for our 2group model.
3.4 Numerical results
The power of the reactor will be maintained constant which means the average neutron flux is fixed. so that Φ_{*,ave} = 3Φ_{*} = 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 (which is the case here) the results given by the one group and the twogroup models look pretty much the same.
(We have observed in our results that ϕ ≅ ψ, but this is no surprise when we consider (17): this results from α _{*} x ≪ 1 and from M _{2} ≪H.)
In Figure 5 we compare one group and twogroup results. We plot AO _{p} as a function of the reduced time τ. 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 timeperiod of oscillations.
Fig. 4 The oscillation of neutron flux in the 2group model with . 
Fig. 5 One group vs twogroup. 
4 The nonlinear 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:(23) (24) (25)
β 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.
4.1 Space approximation
Proceeding like in Section 3, we obtain the following system(26) (27) (28)
Algorithm : we notice that (26)–(28) is a particular case of the following differential system:(29) (30) (31)
which we should complement with initial conditions:
When we do not take the temperature effect into account (i − e when β = 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 (β > 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.
To identify (26)–(28) and (29)–(31) we need to choose
4.2 Time discretization: fully implicit scheme
We have used successfully the following fully implicit scheme.(32) (33) (34)
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 as the smallest eigenvalue of matrix.(35)
where ϕ ^{*} = μ_{n} ϕ^{n} and . We leave to the reader the details of implementation of Newton's method to solve the nonlinear system (32)–(34) of 3N equations for 3N unknowns.
4.3 Time discretization: semi implicit scheme
Step 1: Like in the fully implicit scheme, we evaluate the smallest eigenvalue of the matrix A defined in (35).
Step 2: Knowing y^{n} , x^{n} we compute ϕ^{n} by solving the N × N nonlinear system
(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 (β = 0) is basically the same, except that we avoid step 2;

the results we obtained with both algorithms are almost the same.
4.4 How to choose β?
In PWR 1300 MWe, the rated reactor power corresponds to .
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 β Φ_{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.
Fig. 6 The oscillation of neutron flux in 1 group model with the Doppler effect β = 0.01/3. 
Fig. 7 Evolution of AOp with/without Doppler effect (β = 0.01/3). 
4.5 Shimazu's representation
Like Shimazu [13] we introduce another axial offset for iodine
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}, and α _{*} = 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.
Fig. 8 Shimazu's representation with/without Doppler effect. 
Fig. 9 Effect of Shimazu's control method. 
5 Control laws
5.1 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 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 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:
ϕ _{i} et ϕ _{p}: 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.
γ: stability index of the xenon oscillation ω: 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}).
5.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:(36)
After time discretization, we can obtain a new formula(37)
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 and .
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:(38)
Note that we have selected 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:(39)
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.
6 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 twogroup 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.
Author contribution statement
Bertrand Mercier has provided the equations and the numerical methods. He has also implemented the onegroup model (Sect. 2). Zeng Ziliang has implemented the twogroup model and obtained the associated conclusions (Sect. 3). Chen Liyi has implemented the nonlinear model (Sect. 4). Shao Nuoya has implemented Shimazu's method and found an adequate PID control law (Sect. 5).
Acknowledgments
The authors would like to thank the referee for clever remarks and careful reading.
References
 O. Gustaf, Spatial xenon instabilities in thermal reactors, Report 6910, 1969 [Google Scholar]
 J. Canosa, H. Brooks, Xenon Induced oscillations, Nucl. Sci. Eng. 26 , 237 (1966) [CrossRef] [Google Scholar]
 R.J. Onega, R.A. Kisner, Axial xenon oscillation model, Ann. Nucl. Energy 5 , 13 (1978) [CrossRef] [Google Scholar]
 M.H. Yoon, H.C. No, Direct numerical technique of mathematical programming for optimal control of xenon oscillation in load following operation, Nucl. Sci. Eng. 90 , 203 (1985) [CrossRef] [Google Scholar]
 G. Mathonnière, Etude de problèmes neutroniques liés à la présence du xénon dans les réacteurs à eau pressurisée, Thèse, 1988 [Google Scholar]
 P. Reuss, Neutron physics (EDP Sciences, Les Ulis, 2008), p. 696 [Google Scholar]
 C. Lin, Y. Lin, Control of spatial xenon oscillations in pressurized water reactors via the kalman filter, Nucl. Sci. Eng. 118 , 260 (1994) [CrossRef] [Google Scholar]
 J.M. Pounders, J.D. Densmore, A generalization of λmode xenon stability analysis, PHYSOR 2014 − Kyoto, Japan, September 28–October 3, 2014 [Google Scholar]
 S. Massart, S. Buis, P. Erhard, G. Gacon, Use of 3DVAR and Kalman Filter approches for neutronic state and parameter estimation in nuclear reactors, Nucl. Sci. Eng. 155 , 409 (2007) [CrossRef] [Google Scholar]
 A. Ponçot, Assimilation de données pour la dynamique du xénon dans les coeurs de centrale nucléaire, Thèse, Toulouse, 2008 [Google Scholar]
 S. Marguet, La Physique des Réacteurs Nucléaires (Lavoisier, 2013) [Google Scholar]
 W. Xuejing, W. Jingyi, P. Yang, Un modèle simple pour comprendre les instabilités Xénon, thèse de Bachelor, IFCEN, 2017 [Google Scholar]
 Y. Shimazu, Continuous guidance procedure for xenon oscillation control, J. Nucl. Sci. Technol. 32 , 1159 (1995) [CrossRef] [Google Scholar]
 Y. Shimazu, Monitoring and control of radial xenon oscillation in PWRs by a three radial offset concept, J. Nucl. Sci. Technol. 44 , 155 (2007) [CrossRef] [Google Scholar]
 Y. Shimazu, Xenon oscillation control in large PWRs using a characteristic ellipse trajectory drawn by three axial offsets, J. Nucl. Sci. Technol. 45 , 257 (2008) [CrossRef] [Google Scholar]
 X. Wang, B. Mercier, P. Reuss, A simple 1D 1 group model to study uniaxial xenon instabilities in a large PWR core and their control. HAL CCSD : https://hal.archivesouvertes.fr/hal01625015, 2017 [Google Scholar]
 E. Thibaut, Les oscillations axiales du cœur par effet xenon. Mémoire, Cnam, 2018 [Google Scholar]
 G.S. Lellouche, Space dependent xenon oscillations, Nucl. Sci. Eng. 12 , 482 (1962) [CrossRef] [Google Scholar]
Cite this article as: Bertrand Mercier, Zeng Ziliang, Chen Liyi, Shao Nuoya, Modeling and control of xenon oscillations in thermal neutron reactors, EPJ Nuclear Sci. Technol. 6, 48 (2020)
All Tables
All Figures
Fig. 1 The oscillation of neutron flux in 1 group model with M ^{2} = 64 cm^{2}. 

In the text 
Fig. 2 Evolution of the axial offset AO _{p} for various values of and M ^{2} = M2. 

In the text 
Fig. 3 Evolution of AO_{p} for various values of height. 

In the text 
Fig. 4 The oscillation of neutron flux in the 2group model with . 

In the text 
Fig. 5 One group vs twogroup. 

In the text 
Fig. 6 The oscillation of neutron flux in 1 group model with the Doppler effect β = 0.01/3. 

In the text 
Fig. 7 Evolution of AOp with/without Doppler effect (β = 0.01/3). 

In the text 
Fig. 8 Shimazu's representation with/without Doppler effect. 

In the text 
Fig. 9 Effect of Shimazu's control method. 

In the text 
Fig. 10 Effect of our PID control strategy (38) with K_{p} = −75, T_{i} = 50, T_{d} = 0.01. 

In the text 
Fig. 11 Effect of our PID control strategy (39) with K_{d} = 130. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.