© B. Mercier and V. Borysenko, Published by EDP Sciences, 2023

1. Introduction

On the web, we find https://www.nuclear-power.com/nuclear-power/reactor-physics/reactor-dynamics/reactor-stability/ which states:

“The response of a reactor to a change in temperature (i.e., the overall reactor stability) depends especially on the algebraic sign of α_T. A reactor with negative α_T is inherently stable to changes in its temperature and thermal power, while a reactor with positive α_T is inherently unstable.”

In the present paper, we consider the case of a reactor for which α_T is positive at low power and negative at large power. This case is of interest since this was how the RBMK behaved at the time of the Chernobyl accident.

Let k be the neutron effective multiplication factor.

Let ρ = k − 1/k denote the reactivity.

Obviously when ρ > 0 (resp. ρ < 0 the reactor is supercritical (resp. subcritical). In other words, the linear theory predicts that the reactor is unstable. It also tells us that when the reactor is critical (i.e. ρ = 0) to increase its power, one has to make a positive reactivity injection ρ₀ > 0, e.g. by moving up the control rods, and then, when the power has reached the desired power level, a negative reactivity injection −ρ₀ e.g. by moving down the control rods where they were before.

This is not, as we shall see, how it works in a real nuclear reactor where there is a reactivity feedback.

Let n denote the normalized reactor power such that n = 1 at nominal power, the linear theory predicts that

(1)

where τ is the effective neutron lifetime, which means that n(t) increases or decreases exponentially.

In a PWR, when the power increases, the fuel temperature increases inducing the Doppler effect. The moderator temperature also increases, inducing a moderator effect.

In other words, there is a negative reactivity feedback α(n) < 0 such that

(2)

Therefore ρ is not a constant and we should replace (1) by

(3)

Note that a point reactor model is used for n; τ takes into account the delayed neutrons; α takes into account both moderator and Doppler effects and is supposed to depend only on n.

In view of

we obtain

(4)

Finally, we have to study the 2 × 2 nonlinear differential system:

(5) (6)

Note that we can exclude the particular case n₀ = 0 which is not of interest.

2. Reactor stability

If, for a while, we consider ρ as a function of n, we have the differential equation (2) which should be complemented by the “initial” condition

(7)

So that we have

(8)

which holds both for n > n₀ and for n < n₀.

The reactivity injection ρ₀ can be positive (in which case the reactor power increases) or negative (in which case it decreases).

The reactor stability requires, as we shall see, that there exists α_* such that:

(9)

Indeed, with such an assumption, for all starting point {n₀, ρ₀}, the curve:

will cut the horizontal axis once and only once in the plane {n, ρ}. In other words, the equation

(10)

has a unique solution. Note, however, that we may have n_∞ < 0.

The particular case where α(n) is constant has been considered in Reuss-Mercier [1].

In such a case, we have

so that the curve n → ρ = ρ(n₀, ρ₀, n) is a straight line.

In this section, we consider the general case described by (9).

Proposition 1: Under the above assumption (9), then n(t) → max(0, n_∞) when t → ∞.

Proof. We divide (4) by ρ and (3) by n so that we get

(11) (12)

1°/Case where ρ₀ > 0.

Equation (12) shows that n(t) increases, at least in the neighborhood of t = 0, but can reach a maximum only if ρ vanishes which means n(t) = n_∞.

As long as n(t) < n_∞, we have ρ(t) = ρ(n₀, ρ₀,n(t)) > 0.

Since we have

we deduce that Logρ(t) → −∞ if t → ∞ so that ρ(t) → 0.

The curve n → ρ = ρ(n₀, ρ₀, n) being continuous, we deduce that n(t) → n_∞.

2°/Case where ρ₀ < 0.

In an analogous way, (12) shows that n(t) must be decreasing, at least in the neighborhood of t = 0 but can reach a minimum only if ρ vanishes, then if n = n_∞.

Let us assume first that n_∞ > 0.

We have n_∞ ≤ n(t) ≤ n_o and then α(n(t)).n(t) < α_*.n(t) < α_*.n_∞ < 0.

Therefore Log|ρ(t)| → −∞ and then |ρ(t)| → 0 so that n(t) → n_∞.

Let us now assume that n_∞ < 0.

The curve n → ρ(n₀, ρ₀, n) cuts the vertical axis n = 0 in ρ = ρ_* < 0.

We deduce that ρ(t) increases and that ρ(t) → ρ_*.

Since ρ₀ ≤ ρ ≤ ρ_* < 0 we deduce from (12) that Log n(t) → −∞ and then n(t) → 0.

With this result, we have shown that for a stable reactor, due to the reactivity feedback, to increase the power of the reactor, it is sufficient to make an appropriate reactivity injection ρ₀, but it is not necessary to make a negative reactivity injection equal to −ρ₀ after this positive reactivity injection.

3. Evaluation of α(n) for the RBMK

It is well known that, before the Chernobyl accident, there was a significant void effect on the RBMK. This void effect originated from the fact that graphite is used as a moderator and water as a coolant. Since light water absorbs thermal neutrons, if it disappears, the core reactivity increases.

More precisely, in [2, p.55] we find that the Doppler effect brings a negative reactivity of −1000 pcm and that the void coefficient may be as high as +2500 pcm. (This was the case just before the accident).

We assume that the Doppler effect is proportional to n(e.g. −100 pcm for n = 0.1.

To evaluate the coolant effect, we assume that

• the coolant pressure is 7 MPa, then the saturation temperature is 285.8 °C.

• As for the enthalpy, we have h_l,sat = 1267 kJ/kg for saturated liquid h_v,sat = 2772 kJ/kg for saturated steam.

• With obvious notations, we have ρ_l,sat = 739 kg/m³ and ρ_v,sat = 36 kg/m³.

• The associated specific volumes are τ_l,sat = 1.35 L/kg and τ_v,sat = 27.8 L/kg.

• The mass flow rate is q = 10417 kg/s.

The inlet enthalpy is adjusted in such a way that the saturation level is at z_sat = 1.855 m (note that H = 7 m is the height of the core).

To obtain the values given in Table 1, we proceed in the following way.

Given the reactor power P, we calculate the enthalpy increment per meter Δh = P/(H.q).

We calculate h_in = h_{l,sat − zsat} * Δh and h_out = h_l,sat + (H − z_sat) * Δh.

Then we calculate the steam mass fraction

x_out = (h_out − h_l,sat) /(h_v,sat − h_l,sat), τ_out = (1 − x_out) τ_l + x_outτ_v and ρ_out = 1000/τ_out.

We assume that ρ(z) is constant for 0 ≤ z ≤ z_sat and then varies linearly from ρ_l,sat for z = z_sat to ρ_out for z = H.

We calculate its average ρ_ave in the pressure tubes and then the coolant effect (in pcm).

Coolant effect = 2500 * (ρ_l,sat − ρ_ave) /ρ_l,sat

(In the limit case ρ_ave = 0 we would get that the coolant effect would be equal to the void effect).

Starting from criticality at n = 0, the core reactivity for n ≠ 0 is computed as the algebraic sum of the (positive) coolant effect and the (negative) Doppler effect. It is computed in column 6 of Table 1. The derivative in pcm/MW is computed in column 7 and plotted in Figure 1.

Fig. 1. n → α(n) for the RBMK before the accident.

Finally, our α(n) would be obtained by multiplication of column 7 by 3200.

What it means is that rather than being always negative, α(n) is positive for n < n_* and α(n) < 0 when n > n_* with n_* = 0.238 that is P = 762 MW.

Clearly, condition (9) is violated.

In Section 4, we shall prove that, in such a case, the reactor is unstable and we shall study how it is unstable.

4. Behavior of the unstable reactor

In this section, we shall assume that there exists n_* such that

0 < n_* < 1 and α(n_*) = 0

α(n) > 0 for n < n_*

α(n) < 0 when n > n_*.

The curve n → ρ = ρ (n₀, ρ₀, n) that we have defined in (8) is then increasing when n < n_* and strictly decreasing when n > n_* as indicated in Figure 2. Equation (10) may then have 2 solutions: one smaller than n_* and the other one greater. We shall prove that only the latter (which we denote by n_∞) is stable. We shall prove that

• if ρ₀ > 0 then n(t) → n_∞ when t → ∞,

• if ρ₀ < 0 and n₀ > n_*, then n(t) → n_∞ when t → ∞,

• if ρ₀ < 0 and n₀ < n_*, then n(t) → 0 when t → ∞.

Fig. 2. Example of curve n → ρ = ρ (n0, ρ0, n).

Let us first consider the case where ρ₀ > 0.

Let us choose n₁ > n_* such that ρ(n₀, ρ₀, n₁) ≥ ρ₀ (see Fig. 2) and then

From (12) we deduce that

and then, as long as n(t) ≤ n₁

Then the set of real numbers t such that n(t) ≤ n₁ is bounded. Then, there exists a time t₁ such that n(t₁) = n₁.

Let ρ₁ = ρ(n₀, ρ₀, n₁); we have ρ(t₁) = ρ₁.

Let α_* = min {ρ(n):n > n₁} with the above assumptions, we have α_* < 0.

Since α(n) < α_* < 0Ɐn > n₁ we can apply proposition 1 which tells us that n(t) → n_∞.

The case where ρ₀ < 0 and n₀ < n_* is illustrated in Figure 3.

Fig. 3. Curve n → ρ = ρ (n0, ρ0, n) obtained for n0 = 0.05 and ρ0 = −26.6 pcm.

Equation (12) shows that n decreases and hence ρ also decreases. We see that Log n → −∞ hence n → 0.

To be complete we should also consider the case where n₀ > n_* and ρ₀ < 0.

We can just apply proposition 1, which shows that ρ(t) → 0 and n(t) → n_∞.

We have thus proved the desired result.

Note that we have proved that the starting points n₀ < n_* with ρ₀ = 0 are unstable:

• if rather than being equal to zero, ρ₀ is slightly positive we shall have n(t) → n_∞ > n_*.

• If ρ₀ is slightly negative we shall have n(t) → 0.

Note also that with 3 starting points n₀₁ = 0.05, n₀₂ = 0.12 and n₀₃ = 0.5 and the same reactivity injection ρ₀ = 25 pcm, as can be seen in Figure 4 the final power n_∞ may be larger when the initial power n₀ is smaller.

Fig. 4. 3 starting points n01, n02, n03 with the same ρ0 = 25 pcm.

5. Conclusion

With our simple model, we have shown that reactivity feedback is an important phenomenon.

When the reactor is stable, to increase its power, it is sufficient to make an appropriate reactivity injection ρ₀, but, as the operators know it, it is not necessary to make a negative reactivity injection equal to −ρ₀ after this positive reactivity injection.

We have also studied one case of a reactor that is unstable at low power, like the RBMK.

We have shown the practical effects of such an instability:

• 1)

  with the same reactivity injection, the lower the initial power then the final power may be higher,

• 2)

  the line ρ₀ = 0 0 < n₀ < n_* is a bifurcation line.

• If we start from ρ₀ slightly positive we have n(t) → n_∞ > n_*;

• if we start from ρ₀ slightly negative we have n(t) → 0.

This was the case for the RBMK before the Chernobyl accident: the operators were obliged to frequently use the control rods, and, practically, remaining at low power was not comfortable at all for them. This is one explanation for the accident, but of course, it is not the only one [3, 4].

Note that our 2 × 2 model can only give a qualitative behavior of the reactor; to get accurate kinematics, we would have to adjust the effective neutron lifetime τ or use the 8 × 8 model with delayed neutrons, as described in [1].

As has been shown in [1], the asymptotic behavior is the same for both models.

Our results explain better why Anatoly Dyatlov has written in https://www.neimagazine.com/features/featurehow-it-was-an-operator-s-perspective/:

“The nuclear safety department of the power plant (…) measured the fast power coefficient at power levels close to nominal full power, ie in the region where it was negative. The results were used by the operators in their everyday work. The latest data prior to the accident gave a value of minus 1.7 × 10⁻⁴ ß/MW”.

Conflict of interests

The authors declare that they have no competing interests to report.

Funding

This research did not receive any specific funding.

Data availability statement

This article has no associated data generated and/or analyzed.

Author contribution statement

1. Introduction B. Mercier – 2. Reactor stability B. Mercier – 3. Evaluation of α(n) for the RBMK V. Borysenko – 4. Behavior of the unstable reactor. B. Mercier.

References

All Tables

All Figures

	Fig. 1. n → α(n) for the RBMK before the accident.
In the text

	Fig. 2. Example of curve n → ρ = ρ (n0, ρ0, n).
In the text

	Fig. 3. Curve n → ρ = ρ (n0, ρ0, n) obtained for n0 = 0.05 and ρ0 = −26.6 pcm.
In the text

	Fig. 4. 3 starting points n01, n02, n03 with the same ρ0 = 25 pcm.
In the text