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 [Image omitted: See PDF.] 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 [Image omitted: See PDF.].
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 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 one-group 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 test the well-known 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:
[Image Omitted: See PDF.]
(1)
where H is the height of the core. Φ is the neutron flux of the core in n/cm2/s and M2 denotes the migration area in cm2. This differential equation should be completed by some boundary conditions (BC). In the following we shall choose the Dirichlet BC:
[Image Omitted: See PDF.]
(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 one-group 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:
[Image Omitted: See PDF.]
(3)
[Image Omitted: See PDF.]
(4)
where:
σ 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 xenon fission yield (about equal to 0.001 for UO2 fuel);
λ (resp. µ) the time constant for the β − decay of [Image omitted: See PDF.] (resp. [Image omitted: See PDF.]).
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/cm2/s.
The core is assumed to be homogenized so that, in each cm3 of core, there are ω cm3 of fuel and (1 − ω) cm3 of moderator. In a PWR we have ω ∼ 1/3 so that the moderation ratio (1 − ω)/ω ∼ 2.
Since there are ωX at/cm3 of [Image omitted: See PDF.] in the core, the xenon macroscopic cross section is equal to ωσX Once taking into account the neutron capture by [Image omitted: See PDF.] equation (1) has to be replaced by
[Image Omitted: See PDF.]
(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
[Image Omitted: See PDF.]
(6)
[Image Omitted: See PDF.]
(7)
Finally, let [Image omitted: See PDF.], we get
[Image Omitted: See PDF.]
(8)
X* has the following interpretation: it is the equilibrium value for the concentration of 135Xe 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.
[Image omitted: See PDF.]
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 column-vector 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
[Image Omitted: See PDF.]
(9)
[Image Omitted: See PDF.]
(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 {ya, xa} 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:
First step, solve
[Image Omitted: See PDF.]
(11)
[Image Omitted: See PDF.]
(12)
Like in Section 2, ωσX is the xenon macroscopic capture cross section : it should be added to Σ2.
In the same way, if Σf1 (resp. Σf2) 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 cm3 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 cm3 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
[Image Omitted: See PDF.]
(13)
[Image Omitted: See PDF.]
(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
[Image Omitted: See PDF.]
(15)
but it is also possible to consider the Fourier boundary conditions.
Remark: The number of neutrons produced by fast fission of [Image omitted: See PDF.] is actually 2.8, but, for the sake of simplicity, we choose ν = 2.4 both for [Image omitted: See PDF.] and [Image omitted: See PDF.]. 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 2-group model
Like in Section 2 the system (11)–(14) can be simplified.
At first, let
[Image omitted: See PDF.]
Then, we rewrite the variables Φ = Φ*ϕ, Ψ = Ψ*ψ, X = X*x, Y = Y*y, τ = λt.
And we obtain
[Image Omitted: See PDF.]
(16)
[Image Omitted: See PDF.]
(17)
[Image Omitted: See PDF.]
(18)
[Image Omitted: See PDF.]
(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
[Image Omitted: See PDF.]
(20)
[Image Omitted: See PDF.]
(21)
[Image Omitted: See PDF.]
(22)
The parameter k∞,1 being given, let B = B(x, k∞, 2) denote the 2N × 2N matrix defined by
[Image Omitted: See PDF.]
(23)
[Image Omitted: See PDF.]
(24)
[Image Omitted: See PDF.]
(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
[Image Omitted: See PDF.]
(26)
[Image Omitted: See PDF.]
(27)
[Image Omitted: See PDF.]
(28)
Algorithm : we notice that (26)–(28) is a particular case of the following differential system:
[Image Omitted: See PDF.]
(29)
[Image Omitted: See PDF.]
(30)
[Image Omitted: See PDF.]
(31)
where
[Image Omitted: See PDF.]
(32)
[Image Omitted: See PDF.]
(33)
[Image Omitted: See PDF.]
(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 [Image omitted: See PDF.], the reactivity is adjusted to zero not by the Doppler effect, but by the boron concentration, which means that the control system selects [Image omitted: See PDF.] as the smallest eigenvalue of matrix.
[Image Omitted: See PDF.]
(35)
where ϕ* = μnϕn and [Image omitted: See PDF.]. 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.
4.3 Time discretization: semi implicit scheme
Step 1: Like in the fully implicit scheme, we evaluate the smallest eigenvalue [Image omitted: See PDF.] of the matrix A defined in (35).
Step 2: Knowing yn, xn we compute ϕn by solving the N × N non-linear system
[Image Omitted: See PDF.]
(36)
After time discretization, we can obtain a new formula
[Image Omitted: See PDF.]
(37)
where:
U: variable that can be controlled
err (or e): error (that is AOp)
T: control period (in the following T = 0.02 in reduced units)
Kp: proportional coefficient.
In the following, we let [Image omitted: See PDF.] and [Image omitted: See PDF.].
Since U denotes the control rod's position and e the axial offset AOp, we can describe the position of control bars by this equation:
[Image Omitted: See PDF.]
(38)
Note that we have selected [Image omitted: See PDF.] as the initial position of our control bars.
Figure 10 shows the results obtained with PID when M2 = 100 cm2; Kp = −75, Ti = 50, Td = 0.01.
When we use PID regulator to control the system, we found that Td is a sensitive parameter. Let's just consider the differential item: Kd · [e(k) − e(k – 1)].
We let Ti = 0 and we neglect the proportional item, then we can describe the position of control bars by this equation:
[Image Omitted: See PDF.]
(39)
which means, practically, that when AOp increases, the control bars should be inserted.
We note that Kd is the more sensitive parameter. Figure 11 corresponds to Kd = 130.
We can conclude that an adequate PID regulator can also control the system quite well.
[Image omitted: See PDF.]
[Image omitted: See PDF.]
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 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.
Author contribution statement
Bertrand Mercier has provided the equations and the numerical methods. He has also implemented the one-group model (Sect. 2). Zeng Ziliang has implemented the two-group model and obtained the associated conclusions (Sect. 3). Chen Liyi has implemented the non-linear 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.
1. O. Gustaf, Spatial xenon instabilities in thermal reactors, Report 6910, 1969
2. J. Canosa, H. Brooks, Xenon Induced oscillations, Nucl. Sci. Eng. 26 , 237 (1966)
3. R.J. Onega, R.A. Kisner, Axial xenon oscillation model, Ann. Nucl. Energy 5 , 13 (1978)
4. 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)
5. 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
6. P. Reuss, Neutron physics (EDP Sciences, Les Ulis, 2008), p. 696
7. C. Lin, Y. Lin, Control of spatial xenon oscillations in pressurized water reactors via the kalman filter, Nucl. Sci. Eng. 118 , 260 (1994)
8. J.M. Pounders, J.D. Densmore, A generalization of λ-mode xenon stability analysis, PHYSOR 2014 − Kyoto, Japan, September 28–October 3, 2014
9. S. Massart, S. Buis, P. Erhard, G. Gacon, Use of 3D-VAR and Kalman Filter approches for neutronic state and parameter estimation in nuclear reactors, Nucl. Sci. Eng. 155 , 409 (2007)
10. A. Ponçot, Assimilation de données pour la dynamique du xénon dans les coeurs de centrale nucléaire, Thèse, Toulouse, 2008
11. S. Marguet, La Physique des Réacteurs Nucléaires (Lavoisier, 2013)
12. W. Xuejing, W. Jingyi, P. Yang, Un modèle simple pour comprendre les instabilités Xénon, thèse de Bachelor, IFCEN, 2017
13. Y. Shimazu, Continuous guidance procedure for xenon oscillation control, J. Nucl. Sci. Technol. 32 , 1159 (1995)
14. Y. Shimazu, Monitoring and control of radial xenon oscillation in PWRs by a three radial offset concept, J. Nucl. Sci. Technol. 44 , 155 (2007)
15. 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)
16. 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.archives-ouvertes.fr/hal-01625015, 2017
17. E. Thibaut, Les oscillations axiales du cœur par effet xenon. Mémoire, Cnam, 2018
18. G.S. Lellouche, Space dependent xenon oscillations, Nucl. Sci. Eng. 12 , 482 (1962)
Bertrand Mercier1, Zeng Ziliang2, Chen Liyi2 and Shao Nuoya2
1 CEA/INSTN, 91191 Gif sur Yvette, France
2 Institut Franco-Chinois de l'Energie Nucléaire, Sun Yat-sen University, 519082 Zhuhai Campus, Guangdong Province, P.R. China
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer
© 2020. This work is licensed under https://creativecommons.org/licenses/by/4.0 (the “License”). Notwithstanding the ProQuest Terms and conditions, you may use this content in accordance with the terms of the License.
Abstract
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 M2 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 M12+M22 where M12 denotes the migration area for fast neutrons and M22 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.
You have requested "on-the-fly" machine translation of selected content from our databases. This functionality is provided solely for your convenience and is in no way intended to replace human translation. Show full disclaimer
Neither ProQuest nor its licensors make any representations or warranties with respect to the translations. The translations are automatically generated "AS IS" and "AS AVAILABLE" and are not retained in our systems. PROQUEST AND ITS LICENSORS SPECIFICALLY DISCLAIM ANY AND ALL EXPRESS OR IMPLIED WARRANTIES, INCLUDING WITHOUT LIMITATION, ANY WARRANTIES FOR AVAILABILITY, ACCURACY, TIMELINESS, COMPLETENESS, NON-INFRINGMENT, MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. Your use of the translations is subject to all use restrictions contained in your Electronic Products License Agreement and by using the translation functionality you agree to forgo any and all claims against ProQuest or its licensors for your use of the translation functionality and any output derived there from. Hide full disclaimer