© B. Mercier et al., published by EDP Sciences, 2020

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 beta-decay 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 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: (1)

where H is the height of the core. Φ is the neutron flux of the core in n/cm²/s and M ² denotes the migration area in cm². 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 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:

 (3) (4)

where:

• Y (resp. X) is the concentration of (resp. ) in the fuel evaluated in at/cm³ 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⁻¹ and Φ in n/cm²/s.

The core is assumed to be homogenized so that, in each cm³ of core, there are ω cm³ of fuel and (1 −   ω) cm³ of moderator. In a PWR we have ω ∼ 1/3 so that the moderation ratio (1 −   ω)/ω ∼ 2.

Since there are ωX at/cm³ 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

(6) (7)

Finally, let , we get(8)

X_* has the following interpretation: it is the equilibrium value for the concentration of ¹³⁵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.

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

where we assume that ϕ ₀ = 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:

First step, solve

Second step, solve

Which produces both a value for (the smallest eigenvalue) and a (non-unique) value for ϕ ⁿ .

Then we replace ϕ ⁿ 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¹³ n/cm²/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 ² = 100 cm² and ;

• stable when M ² = 64 cm² and ;

• growing when M ² = 100 cm² and or M ² = 64 cm² 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². 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 cm2.

Fig. 2 Evolution of the axial offset AO p for various values of and M 2 = M2.

Fig. 3 Evolution of AOp for various values of height.

3 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¹³ n/cm²/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¹⁴ n/cm²/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¹³ fissions/cm³/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 Σ₂.

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³ per s in the core.

Note that the xenon capture cross section is negligible for fast neutrons.

Now, Σ₁Ψ is the number of fast neutrons which disappear per cm³ 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 ν ₁ ω Σ _f1Ψ where ν ₁ 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

Then, we rewrite the variables Φ = Φ_* ϕ, Ψ  =  Ψ_* ψ,  X = X _* x,  Y = Y _* y,  τ  =  λt.

And we obtain

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 Perron-Frobenius theorem to matrix (aI+B)⁻¹ 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 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%.

3.3 Time discretization: semi implicit scheme

We proceed as in Section 2 and introduce the following semi implicit scheme:

First step, solve

Second step, solve

which produces both a value for and a (non-unique) value for ϕ ⁿ .

Then we replace ϕⁿ 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.

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¹³ n/cm²/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 two-group 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 ₂ ≪H.)

In Figure 5 we compare one group and two-group 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 time-period of oscillations.

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

Fig. 5 One group vs two-group.

4 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:(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)

where

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 ϕⁿ and . 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 of the matrix A defined in (35).

Step 2: Knowing yⁿ , xⁿ we compute ϕⁿ by solving the N × N non-linear system

(For this, we use Newton's method).

Step 3: We compute yⁿ⁺¹ , xⁿ⁺¹ 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

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 ² =100 cm², 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 shut-down. 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 ² = 100 cm²).

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²; 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.

Fig. 10 Effect of our PID control strategy (38) with Kp = −75, Ti = 50, Td = 0.01.

Fig. 11 Effect of our PID control strategy (39) with Kd = 130.

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.

References

All Tables

All Figures

	Fig. 1 The oscillation of neutron flux in 1 group model with M 2 = 64 cm2.
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 AOp for various values of height.
In the text

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

	Fig. 5 One group vs two-group.
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 Kp = −75, Ti = 50, Td = 0.01.
In the text

	Fig. 11 Effect of our PID control strategy (39) with Kd = 130.
In the text