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

1 Introduction

A better understanding of former nuclear accidents is obviously good for nuclear safety.

Many countries are developing nuclear reactor projects, and it is useful for them (and more generally for the world) to benefit from an adequate feedback from past events.

For teaching purpose, also, it is useful to show which conclusions can be obtained with simple models.

In the INSAG-7 report on the Chernobyl accident [1], the following design features of the RBMK have been identified as having a major responsibility:

• void coefficient of reactivity

• design and control of safety rods

• speed of insertion of the emergency control rods

• control of power

However, this does not explain why, when the operators pressed the AZ-5 button to command the reactor scram, nothing happened during 3 s and then the reactor exploded (see e.g. https://www.neimagazine.com/features/featurehow-it-was-an-operator-s-perspective/).

The RBMK is a kind of thermal-neutron reactor designed and built by the former Soviet-Union. It is an early Generation II reactor that dates back to 1954, some of which are still in operation in ex-USSR. In the RBMK reactor core, the fast neutrons generated from the fission of ²³⁵U are moderated by graphite so that they can easily actuate other fissions. Water is used as a coolant to transfer the heat from inside the core to outside. From a neutronics point of view, light water is both an absorber and a moderator. Since there is enough graphite to ensure moderation, the contribution of water to moderation is negligible, but this is not the case for absorption. This explains the positive void coefficient of reactivity which has been evaluated before by e.g. Kalashnikov [2].

In Section 2 we shall use TRIPOLI 4, which implements the Monte-Carlo method, to evaluate the safety rods efficiency. The Monte-Carlo method is not a simplified method, but we shall apply it on simplified geometries with a network of 3 × 3 fuel elements and lateral reflective boundary conditions as if the diameter of the core was infinitely large, or the radial reflectors fully efficient.

For the safety rod efficiency we find like in INSAG-7 that they may induce “a local insertion of positive reactivity in the lower part of the core.” And that “the magnitude of this ‘positive scram' effect depended on the spatial distribution of the power density and the operating regime of the reactor”.

In Sections 3–6, we shall study the stability and the reactivity control of the RBMK reactor under xenon poisoning by using 1D 1group axial simulations (see [3]) to prove that, with a weak power counter-reaction coefficient, such as the one we have on a RBMK of the 2_nd generation, the reactor is rather unstable.

At the same time, we shall propose a simple model to simulate the automatic control system of the core, which the operators have switched off before the disaster. If the core is unstable with the presence of the automatic system, then it would be highly possible that the manual steering of the staff has made the core more unstable and that possibly initiated the accident.

Then Section 7 is devoted to the hydrodynamic simulation of the coolant water in the pressure tubes which starts with a strong shock propagation and then, the depressurization process. We consider these two phases successively. The first phase is due to the fast release of a significant amount of energy in the liquid coolant which will be instantaneously transformed in (very) high pressure at the origin of a strong propagating shock in the pressure tubes. To simulate this phase we use a one-dimensional Lagrangian model to solve the system of conservation laws for mass, momentum and energy. During this phase we shall show that, if we start from a liquid state, the coolant will stay in the liquid phase, so that we can use a stiffened gas equation of state (see Corot-Mercier [4]). Our purpose is to estimate the shock wave velocity and to estimate the (short) duration of this phase. We conjecture that the shock will lead to the failure of the pressure tubes at their junction with the drum separators. Then this so-called steam explosion will be followed by a depressurization process in the 24 m long pressure tubes. Such a depressurization is also called flash evaporation. This is an isentropic process with liquid to steam phase change. We also use a 1D hydrodynamic model but we take into account the conservation of entropy to make the computational model simpler.

After the depressurization process is over, the water density in the pressure tubes will decrease a lot, so that the positive void coefficient of reactivity produces a strong reactivity accident. In Section 8 we evaluate the fission energy released by such a reactivity accident.

We use the Nordheim-Fuchs model [5,6], which is a 0D model and gives 165 GJ.

This is not the total energy released by the explosion, since we should also add the energy released through the depressurization process.

2 Neutron transport calculation: simplified model

Figure 1 extracted from Wikipedia represents a schematic view of the core of the Chernobyl-4 RBMK. The radial reflector (made of 4 rows of pure graphite blocks) is not represented.

211 control bars are represented: 32 are short and inserted from below whereas 179 are long and inserted from above.

The division of those 179 control rods between automatic and emergency protection groups was arbitrary; on Figure 1 the division is 167 + 12.

When control rods are inserted into the core, they should absorb lots of thermal neutrons effectively so that the chain reaction would be weaker, and the reactivity would be reduced accordingly. However, a graphite rod whose length is 4.5 m called “displacer” was attached to the end of the rod of B₄C. The graphite has a much smaller cross section of absorption than B₄C, and even than water. This absorber-displacer design has two advantages. One is to give more control range because when the displacer is in the core, it could give more reactivity, and the absorber has the opposite effect. The other is to prevent coolant water from entering the space vacated when the absorber is withdrawn, which enables a smoother thermal flux in the core to reduce local stress. Nonetheless, this design has a fatal flaw. As is well known, xenon oscillations induce periods where the neutron flux is higher in the lower part of the core and periods where it is higher in the top of the core. In the former case, after over-extraction of the control rod, when it is re-inserted into the core, the graphite tip would induce a lot of reactivity in the lower part of the core and thus increase the power significantly. No doubt that this is an essential contribution to the accident. (see GRS 121 [7], p. 44).

To check that this is the case, we shall use a Monte-Carlo code (Tripoli-4R) but on a simplified geometry. In view of Figure 1, it seems appropriate to consider a 3×3 network containing two kinds of components which are the fuel channel (FC) and channel of control protection system (CPS), both of which are axial models. Our 3×3 model will involve 7 FC and 2 CPS (see Fig. 2).

Our approach is similar to the one used in Kalashnikov [2] who considers also a 3×3 network but with 6 FC, 2 CPS and 1 AA (Additional Absorber). In fact, his target was to prove that, to reduce the void effect on a RBMK, it is a good idea to both increase the fuel enrichment and introduce additional absorber.

Our target is to identify the situations where the insertion of control rods can increase (and not reduce) reactivity.

The sectional view of a fuel channel is shown in Figure 3. FC consists of the fuel assembly and the pressure tube surrounded by the graphite. The assembly includes 18 fuel rods, which are arranged in the pressure tube in two concentric rings: the internal ring is formed by 6 rods (internal rods) and the outer ring is formed by 12 rods (outer rods).

The size of a graphite block is 25×25 cm. The pressure tube is situated between two co-axial cylinders (R = 4 cm and R = 4.4 cm). The radius of each fuel pin is 0.592 cm. We now describe how we modelize a CPS. The graphite block and the pressure tube have the same size as for FC. The B4C ring is situated between two co-axial cylinders (R = 2.52 cm and R = 3.28 cm). This is for the absorbing part, in the upper part of the control rod (see Fig. 4).

However, in the lower part, there is successively water (1.25 m), a graphite extension (4.5 m) and again water (1.25 m). When the rods are in a high position, their insertion begins by replacing the water in the lower part of the channel with graphite.

More precisely, Figure 5 represents a sectional view of the graphite displacer when it is located in the core.

An axial view of the control rod falling process is shown on Figure 6.

This process will bring about changes in reactor reactivity. We know (see INSAG-7) that there was xenon poisoning in the reactor at the time of the Chernobyl accident. In this section, we will study the effects of the insertion of the control rod on the reactivity of the reactor in two situations: normal operation and the presence of xenon poisoning.

In order to study the influence of the control rod insertion under different conditions, we shall use two models to simulate the states when the control rods are raised to the upper limit and when the control rods are initially inserted. For a reactor under normal operating conditions, the schematic diagrams of the two states are shown in Figures 7 and 8. Moreover, we also consider the xenon poisoning, which brings a heavy accumulation of ¹³⁵Xe in the upper half. In this case, the schematic diagrams of the two states are shown in Figures 9 and 10.

For the calculation with poisoning, we take 4.6 × 10¹⁵ atom of ¹³⁵Xe per cubic centimeter of fuel.

Our TRIPOLI-4 (R) simulation results give values of k _∞ which are shown in Table 1.

We find that when the upper half of the core is in a state of xenon poisoning, the insertion of the control rod will cause an increase of k _∞ then an increase in reactivity (+396 pcm with Xenon poisoning while it is −344 pcm without Xenon poisoning).

In the Chernobyl accident, when the operators saw the initial power increase, they pressed the emergency shutdown button AZ-5, causing the control rods to be inserted into the core. The simulation results above indicate that as the graphite followers displace water in the lower part of the core when the upper half of the core is in a state of xenon poisoning, they locally increase the reactivity, which led to a surge in energy and triggered the first explosion.

Fig. 1 Schematic view of the core (from Wikipedia).

Fig. 2 Unit representative of the core.

Fig. 3 Sectional view of FC.

Fig. 4 Sectional view of CPS with the absorber.

Fig. 5 Sectional view of CPS with the displacer.

Fig. 6 Control rod falling process.

Fig. 7 Control rods are raised to the upper limit without poisoning.

Fig. 8 Control rods are initially inserted without poisoning.

Fig. 9 Control rods are raised to the upper limit with poisoning.

Fig. 10 Control rods are initially inserted with poisoning.

2.1 Evaluation of the migration area

We know that (see e.g. [5]), when we use the one group 1D axial model, we have$$k_{eff}=k_{\infty }/\left(1+{\lambda }_0{M}^2\right)$$(1)

where λ ₀ = π ²/H ² is the smallest solution of the eigenvalue problem:$$-\frac{{\text{d}}^2\mathrm{\Phi }}{{\text{dz}}^2}={\lambda }_0\mathrm{\Phi },0<z<\text{H},\mathrm{\Phi }\left(0\right)=\mathrm{\Phi }\left(\text{H}\right)=0.$$

To evaluate the migration area M ² it is appropriate to perform two Monte-Carlo calculations : one with the leakage boundary conditions in z = 0 and z = H, which will give k _eff , and the other one with fully reflective boundary conditions which will give k _∞. Then we use (1). We find M ² = 190 cm².

The ratio M/H is an indicator to measure the stability of a nuclear reactor. According to reference [3], M/H in PWR is approximately 0.025, while M/H in the RBMK reactor obtained here is approximately 0.020, this means that in terms of xenon oscillation, the RBMK core is more unstable than the PWR.

According to [7], in the Chernobyl accident, the staff carried out inappropriate manipulation of the pilot bars, which led to a rapid decline in reactor power. To alleviate the power decline, they raised too many absorber rods. This move put the upper part of the core under xenon poisoning. It can be seen from the above conclusion that the instability of the RBMK core has been one reason for the accident.

3 Stability and reactivity control

It is important to know what happens to the reactivity of the core when the power of the reactor increases. In the book by Barré et al. [8], p. 55 it is explained that, for RBMK, there are at least 3 contributions to the power reactivity coefficient (a) the fuel temperature coefficient which is negative (b) the void effect in the coolant which is positive (c) the graphite temperature coefficient which is positive. Note that when there is a change of power there is a change of flow rate in the pressure tubes (see [9]). In normal operations, film boiling (DNB), and then significant variations of the steam quality x, should be avoided. At full power, we have x = 0.14 at the top of the core. Film boiling starts above x = 0.2. We have no detailed information on the flow rate tuning at intermediate powers, but we can assume that the void effect contribution to the power reactivity coefficient was limited in normal operations.

As specified in [8], in some (incidental) situations, the power coefficient could be positive so that a power increase could cause a power excursion.

As explained in WIKIPEDIA, the inlet temperature in the pressure tubes results from a mixture between feedwater, which is at a low temperature, and recirculating water which is at saturation temperature. At low power, the feedwater flow rate is low so that the inlet temperature is “dangerously close to the saturation temperature.”

As indicated in INSAG-7 [1], p. 105, in our study, we shall consider that the power coefficient is weakly negative, and this is sufficient, as we shall see, for the reactor to be difficult to control.

4 Safety test

The detailed sequence leading to the accident has been described in several papers like INSAG-7 [1] or GRS121 [7].

There are 2 × 4 = 8 main circulation pumps on the RBMK. At full power 2 × 3 = 6 are expected to deliver 8000 m ³/h each. As indicated by INSAG-7 [1], p. 8, fewer pumps are used during normal start-up or shutdown operations. The specific test which was planned required 2 × 4 = 8 pumps to run simultaneously. However for some reasons, the reactor power was much lower than expected, so that the feedwater flow rate (coming from the condenser) was low and the inlet temperature in the core was close to the saturation temperature of the recirculating liquid water. So, when 4 pumps out of 8, as expected by the safety test, suddenly stopped, there has been a sudden jump in reactivity due to the void effect. This may have also been the initiating event.

5 Basic equations describing the xenon poisoning

For simplicity, we shall consider 1D axial oscillations only.

Like in [3] we shall use the following one group model:$$\frac{dY}{\mathrm{d}t}=\gamma \sum \mathrm{\Phi }-\lambda Y$$(2) $$\frac{dX}{\mathrm{d}t}=\lambda Y-\left(\mu +\sigma \mathrm{\Phi }\right)X$$(3) $$-M^2\mathrm{\Phi }″+\left(1+\alpha X+ϵ\right)\mathrm{\Phi }=k_{\infty }\mathrm{\Phi },0<z<H$$(4) $$\mathrm{\Phi }\left(0\right)=\mathrm{\Phi }\left(H\right)=0$$(5)where:

• Y (resp. X) is the concentration of ${}_{53}^{135}I\left(resp.{}_{54}^{135}Xe\right)$ 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 .064 for UO₂ fuel)

• λ (resp. μ) the time constant for the $\beta -\text{decay of}{}_{53}^{135}I\left(\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{p}.{}_{54}^{135}Xe\right)$

• Φ = Φ (z, t) is the thermal neutron flux

• M ² is the migration area

• α = ωσ/∑, where ω is the fuel volume fraction in the core and ∑ the absorption macroscopic cross section of the homogeneous core without xenon poisoning

• ϵ = ϵ (z) is the additional absorption induced by the insertion of control bars.

We shall assume that $k_{\infty }=k_{\infty }^0.\left(1-\beta \mathrm{\Phi }\right)$ and β is the prompt power feedback coefficient. It takes the Doppler effect into account but also, in a simplified way, the moderator effect and the void effect. For stability it is necessary that β > 0.

When X = X (z)   and  ϵ   =   ϵ (z) are given, we notice that equations (4) and (5) is an eigenvalue problem which gives the infinite multiplication factor k _∞ at zero power.

To obtain a non zero power, we must select $k_{\infty }^0>k_{\infty }$ and the core power will depend on $k_{\infty }^0-k_{\infty }$.

The system (2)–(5) is a non linear system of differential equations with a constraint.

Both an implicit and a semi-implicit scheme are introduced in [3] for its solution.

In what follows, we shall use the semi-implicit scheme and the values given in Table 2 for the parameters of our model.

6 Simulation of the automatic control of the rods

Here is what INSAG-7 says “When the reactor was operated at low power (…) the operator had mainly to rely on the (…) ex-core detectors. However (…) they could not indicate the average axial distribution of flux, since they were all located in the mid-plane of the core.” To simulate the instability of the core at low power, we have introduced an appropriate feed-back law for the control rods.

In Section 2 we have specified that there were 32 short control rods inserted from below and an arbitrary number of (long) control rods inserted from above that were used for automatic control. In what follows, we shall consider two types of control:

• non symmetric control where only control rods inserted from above are used for automatic control.

• “Symmetric” control where both control rods inserted from above and control rods inserted from below are used for automatic control. They move from the same number of steps at each time, however, we give them different weight.

In the first case, we shall have:$$ϵ\left(z\right)=\{\begin{array}{c}\hfill 0,0<z\le z_0\hfill \\ \hfill {ϵ}_0,z_0<z<H\hfill \end{array}$$

whereas in the second one (we let z ₁ = H − z ₀) we shall have:$$ϵ\left(z\right)=\{\begin{array}{c}\hfill {ϵ}_1,z\le z_0\hfill \\ \hfill 0,z_0<z\le z_1\hfill \\ \hfill {ϵ}_2,z>z_1\hfill \end{array}\text{.}$$

The principle of our control is to adjust z ₀, that is, the insertion of controls rods, to stabilize the total power. However, for practical reasons, we shall assume that z ₀ can only take discretized values, namely an integer number m ∈ {0, 70} times 10 cm (note that for the core of a RBMK, H = 7 m).

The total power is proportional to the average thermal flux in our equations. Assume that at τ = τ _n we have obtained (Φ ⁿ ) (X ⁿ ) and (Y ⁿ ), and the position of control rods at the last time τ = τ _n−1 is m ⁿ⁻¹.

First, let p = ave (Φ ⁿ )/ave (Φ⁰) we adjust m according to p which is the ratio between the average flux at τ = τ _n and at τ = 0. Then we make m ⁿ  = m ⁿ⁻¹ + Δm by defining Δm as shown in Table 3 above.

The second step of adjustment is to verify if the control rod is in an extreme position. When the control rod is almost totally extracted, which means that m ⁽ⁿ⁻¹⁾ + Δm is small, the reactivity is insufficient for the core. Hence, we propose an increment Δk _∞ to the reactivity as if we extracted totally another control rod, then we reset m to a small value which means we insert the control rod deeply into the core. Moreover, when the control rod is almost totally inserted, we propose a similar method. The details of our model are shown in Table 4.

After these two steps of adjustment, we obtain m ⁽ⁿ⁾ for resolution in the next time step.

6.1 Simulation with non-symmetric control

We have simulated the temporal evolution of the core during 40 h starting from the initial values, with ϵ₀ = 0.01. Figure 11 ‑shows the evolution of the average thermal flux in the core, which also represents the power output of the reactor. We can see that at the beginning, the average flux has a sudden decrease and a following sudden increment. That is the consequence of the insertion of the control rods. The average flux starts to oscillate in a small amplitude but with a high frequency. This relatively stable state lasts for the first 17 h. Remarkably, between t = 7 h and t = 8.5 h the average flux evolves smoothly. Then at t = 17 h the oscillation becomes much more drastic than before. That is caused by our system of unsymmetric control rods.

To understand what happens, we should have a look to Figure 12 which gives the neutron flux profile with respect to time. We see that a xenon oscillation develops: before t = 17 h, the power is higher in the lower part. After, this is just the reverse.

But a one step insertion of the control rods from above (i-e 10 cm) has much more effect on the reactivity when the power is higher in the upper part.

During our whole simulation, the maximum and minimum of the average flux are 3.21 × 10¹³ n/(cm² ⋅ s) and 1.34 × 10¹³ n/(cm² ⋅ s) compared to the initial average flux which was equal to 2.30 × 10¹³ n/(cm² ⋅ s).

We can see that the maximum can reach 140% of the initial value and the minimum can drop to 58% of the initial value. With a little negative feedback, the power could still have such severe fluctuations although we have tried to maintain it as stable as possible. The inadequate operations of the staff at the night before the accident could very possibly make the situation worse. They tried to control the reactor manually but they did not pay enough attention to the instability caused by the accumulation of xenon, so they failed and possibly pushed the reactor to a state out of control.

Fig. 11 Evolution of the average thermal flux (unsymmetric control).

Fig. 12 Evolution of the thermal flux distribution. The colour code is from green (initially) to red (final time).

6.2 Simulation with “symmetric” control

We have selected ϵ₁ = 0.002 and ϵ₂ = 0.005, in other words the weight of the control rods inserted from below is 2.5 times less than the weight of the control rods inserted from above. However it is sufficient to give some stability as shown on Figures 13 and 14.

Using both control rods inserted from below and control rods inserted from above is then a key factor to stabilize the core power. Even so, we see that the frequency of control rods tuning is quite high and this can be achieved with an electronic control only.

Fig. 13 Evolution of the average thermal flux (“symmetric” control).

Fig. 14 “Symmetric” control. The colour code is from green (initially) to red (final time).

7 Hydrodynamic simulation of the coolant in the pressure tubes

In this section we try to explain why, when the operators pressed the AZ-5 button to command the reactor scram, nothing happened for a (short) while.

The reactor being unstable, its power was difficult to control. Moreover, when 4 of the 8 recirculating pumps were stopped, the (limited) void effect induced a reactivity increase so that the operators decided to command the reactor scram.

As we explained before, due to xenon poisoning in the upper part of the core, during the first second, the insertion of safety bars still increased (rather than decreased) the reactivity of the core.

We can assume that, at this point there was no prompt criticality yet, but that the fission power increase has been sufficient to transfer a significant amount of energy to the liquid coolant which has been instantaneously transformed in (very) high pressure state at the origin of a strong propagating shock in the pressure tubes. Assume that the pressure of the coolant was p = 8 MPa and its volumic mass was ρ   =  722  kg/m ³ (i.e. the volumic mass of the liquid at saturation) so that its specific internal energy is e ∼ 1310kJ/kg. Then the instantaneous energy deposition increases e but cannot change ρ until there is a significant expansion of the coolant, then it increases the pressure very much, say from 8 to 200 MPa. The coolant stays then at the liquid state, and the shock propagates in the liquid.

We conjecture that the shock will lead to the failure of the pressure tubes at their junction with the drum separators. Then the shock propagation will be followed by a depressurization process in the 24 m long pressure tubes. Such a depressurization is also called flash evaporation. We call “steam explosion” the process we just described.

7.1 Simulation of the shock propagation

To simulate the shock propagation in the liquid we shall use a one-dimensional Lagrangian model to solve the system of conservation laws for mass, momentum and energy.

During this period of time there is no mass transfer between the liquid and the steam phase, so that we can use a stiffened gas equation of state for both phases.$$p=-\gamma p_{\infty }+\left(\gamma -1\right)\left(e-q\right)\rho .$$(6)

For the liquid water, we select (see Corot [10] Tab. 6.1):$$\gamma =2.35,p_{\infty }=1.E9Pa,q=-1167E3J/kg\text{.}$$

Our purpose is to estimate the shock wave velocity caused by the energy deposition and to estimate its duration.

From B. Després [11], p. 20, in Lagrange coordinates (in 1D), we have to solve the following system (where a denotes the Lagrange coordinate):$$\frac{\partial \tau }{\partial t}-\frac{\partial u}{\partial a}=0$$(7) $$\frac{\partial u}{\partial t}+\frac{\partial p}{\partial a}=0$$(8) $$\frac{\partial E}{\partial t}+\frac{\partial }{\partial a}\left(pu\right)=0$$(9)where τ = 1/ρ and E = e + u ²/2 is the total energy per unit mass.

To solve this system numerically we shall use the so-called acoustic scheme (see [11]). We introduce a mesh of the Lagrange coordinate a such that$$\dots a_{i-2}<a_{i-1}<a_i<a_{i+1}<a_{i+2}\dots \text{.}$$

We also introduce a time discretization as follows$$0<\mathrm{\Delta }t<2\mathrm{\Delta }t\dots <n\mathrm{\Delta }t<\left(n+1\right)\mathrm{\Delta }t<\dots$$

Integrating equations (7)–(9) for a _i  ≤ a ≤ a _i+1 and nΔt ≤ t ≤ (n + 1) Δt and we get$$\frac{{\tau }_{i+1/2}^{n+1}-{\tau }_{i+1/2}^n}{\mathrm{\Delta }t}-\frac{u_{i+1}^{n+\frac{1}{2}}-u_i^{n+\frac{1}{2}}}{a_{i+1}-a_i}=0$$(10) $$\frac{u_{i+1/2}^{n+1}-u_{i+1/2}^n}{\mathrm{\Delta }t}+\frac{p_{i+1}^{n+\frac{1}{2}}-p_i^{n+\frac{1}{2}}}{a_{i+1}-a_i}=0$$(11) $$\frac{E_{i+1/2}^{n+1}-E_{i+1/2}^n}{\mathrm{\Delta }t}+\frac{{\left(pu\right)}_{i+1}^{n+\frac{1}{2}}-{\left(pu\right)}_i^{n+\frac{1}{2}}}{a_{i+1}-a_i}=0$$(12) where ${\tau }_{i+1/2}^n\left(\text{resp.}{u}_{i+1/2}^n\right)$ denotes the (assumed constant) value of τ (resp. u) at time nΔt on the cell a _i  ≤ a ≤ a _i+1.

To evaluate $u_i^{n+\frac{1}{2}}$ or $p_i^{n+\frac{1}{2}}$ we use the fact that for the (linearized) acoustic equations, p − αu where α ≡ ρ . c (and c is the sound speed) is a Riemann invariant propagating from right to left. Similarly p + αu propagates from left to right.

Let$$u_L=u_{i+1/2}^n,u_R=u_{i-1/2}^n\text{and}{u}^{*}=u_i^{n+\frac{1}{2}}$$ $$p_L=p_{i+1/2}^n,p_R=p_{i-1/2}^n\text{and}{p}^{*}=p_i^{n+\frac{1}{2}}.$$

We shall require u ^*, p ^* to be the solution of the 2 × 2 system$$p^{*}+{\alpha }_L{u}^{*}=p_L+{\alpha }_L{u}_L$$ $$p^{*}-{\alpha }_R{u}^{*}=p_R-{\alpha }_R{u}_R.$$

We proceed similarly for the boundary conditions.

From equations (10)–(12) we obtain ${\tau }_{i+1/2}^{n+1},u_{i+1/2}^{n+1}$ and $E_{i+1/2}^{n+1}$.

We then compute the new internal energy$$e_{i+1/2}^{n+1}=E_{i+1/2}^{n+1}-{\left(u_{i+1/2}^{n+1}\right)}^2/2$$

and we use the equation of state (6) to obtain the new pressure $p_{i+1/2}^{n+1}$ at time (n + 1) Δt

In the Chernobyl accident, we consider that there was 1 GJ of energy rapidly released into the coolant. Notice that such an energy corresponds to an increase for 1 s of 1 GW of the reactor power (that is the quarter of the nominal power). We know that the mass of coolant contained in the pressure tubes is about 32t.

We shall assume that the 1 GJ energy deposition increases e from 1310 to1510 kJ/kg for 5t of liquid water. From the equation of state this will increase its pressure from 8 to 200 MPa.

The section of each pressure tube is equal to 50 cm². However the coolant can flow in a section reduced to 22.5 cm² Taking into account the number of fuel channels (1661) we can assume that the 1 GJ energy is deposited on a limited height of 1.8 m that is about 25% of the core height.

We choose a simple model just for this basic study. The length of the fuel channel is 7 m in the core but we extend it to 12 m, by including out of core parts. We assume that the energy is deposited in its central part: by symmetry we shall consider a 6 m long tube with fixed pressure (8 MPa) in z = 0 and prescribed velocity (u = 0) in z = 6m (z is then the abscissa along the tube).

At t = 0 we shall start from ρ (z, 0) = 722kg/m ³ and:$$p\left(z,0\right)=\{\begin{array}{c}\hfill 8\text{MPa},0\le z\le 5.1\text{m}\hfill \\ \hfill 200\text{MPa},5.1\text{m}<z\le 6\text{m}.\hfill \end{array}$$

We take 600 cells each with Δz = 1 cm and we obtain the pressure profiles at t = 0.01ms and t = 0.02ms which are represented on Figure 15 and show that the pressure wave is propagating without attenuation at 2000  m/s.

We check that the pressure decreases behind the shock wave but this is not sufficient to obtain a phase change from liquid to steam. The numerical results show that the volumic mass remains everywhere above 700  kg/m ³.

On the other hand, we have a strong shock wave (100 MPa) which is sufficient to break the pressure tubes outside the core, namely at the junction with the drum separators.

Remark: assuming that the energy deposition takes place in the liquid water is conservative.

A simple calculation shows that if it takes place in a two-phase mixture liquid-with 14% steam, in the same 7.1  m ³ volume, then the pressure will climb to 35 MPa rather than 200 MPa. The fact is that steam is much more compressible than liquid water.

Fig. 15 Steam explosion: pressure profiles in the pressure tube.

7.2 Simulation of the depressurization phase

After the failure of the pressure tubes, the pressure at the open end of the pressure tube suddenly decreases from 8 to 0.1 MPa. The depressurization activates a phase change from liquid to steam which is an isentropic process. We also use a 1D hydrodynamic model but we take into account the conservation of entropy to make the computational model simpler.

After the depressurization process is over, the water density in the pressure tubes will decrease a lot, so that the positive void coefficient of reactivity produces a strong reactivity injection and then prompt-criticity.

To simulate the depressurization, we shall use the Wilkins' scheme (see Ref. in DESPRES [11]) and consider the conservation of mass and momentum only.$$\{\begin{array}{c}\hfill \frac{\partial \tau }{\partial t}-\frac{\partial u}{\partial a}=0\hfill \\ \hfill \frac{\partial u}{\partial t}+\frac{\partial p}{\partial a}=0\hfill \end{array}\text{.}$$

As before, we set the grid of the variable a like the following$$\cdots a_{i-2}<a_{i-1}<a_i<a_{i+1}<a_{i+2}\cdots$$

Then the middle of an element can be denoted by the subscript $i+\frac{1}{2}$ $$a_i<a_{i+\frac{1}{2}}<a_{i+1}.$$

Besides, a discretization of time in these instants is$$0<\mathrm{\Delta }t<2\mathrm{\Delta }t\cdots <n\mathrm{\Delta }t<\left(n+1\right)\mathrm{\Delta }t\cdots$$

A staggered grid is used in the Wilkins' scheme which writes$$\begin{array}{c}\hfill \frac{{\tau }_{i+\frac{1}{2}}^n-{\tau }_{i+\frac{1}{2}}^{n-1}}{\mathrm{\Delta }t}-\frac{u_{i+1}^{n-\frac{1}{2}}-u_i^{n-\frac{1}{2}}}{a_{i+1}-a_i}=0\hfill \end{array}$$(13) $$\begin{array}{c}\hfill \frac{u_i^{n+\frac{1}{2}}-u_i^{n-\frac{1}{2}}}{\mathrm{\Delta }t}+\frac{p_{i+\frac{1}{2}}^n-p_{i-\frac{1}{2}}^n}{a_{i+\frac{1}{2}}-a_{i-\frac{1}{2}}}=0.\hfill \end{array}$$(14)

We denote the position of the grid a _i at time nΔt as $z_i^n$. With the Lagrangian description, we know that the new position of this grid is given by$$\begin{array}{c}\hfill z_i^{n+1}=z_i^n+\mathrm{\Delta }t{u}_i^{n+\frac{1}{2}}\hfill \end{array}.$$(15)

Furthermore, ${\tau }_{i+\frac{1}{2}}^n$ being known we evaluate $p_{i+\frac{1}{2}}^n$ by using the isentropic process.

We assume that the pressure tubes will fail at the junction with the drum separator (Fig. 16).

We estimate that the length is 24 m, including a length of 7 m of the channel in the reactor core. According to the construction of the cooling system (Fig. 16), the whole channel we study is not straight. To simplify our simulation work, we still consider that the model is one-dimensional. As we have seen in Section 4 we can assume that half the pressure tubes will be filled of saturated liquid water. That's where the shock due to the fission energy increase will be the stronger. So these pressure tubes will eventually depressurize as follows. The initial distribution of the volumic mass is$$\rho \left(z,0\right)=705\text{kg}/{\text{m}}^3,0\le z\le 24\text{m}.$$

The pressure is also constant and equal to 7.5 MPa.

Using the same boundary conditions at both ends we get the results in Figure 17 where we represent the tube (24 m) on the right and 20 m of the outside world on the left.

After 1 s, the volumic mass has decreased to ρ =105 kg/m ³ which is sufficient for the reactivity accident to start and to break the other remaining pressure tubes.

Fig. 16 Cooling system of the RBMK reactor.

Fig. 17 Evolution of volumic mass along the 24 m pressure tubes.

8 Evaluation of the energy released caused by a reactivity accident

We shall use the Nordheim-Fuchs model [5]. Let P denote the fission power:$$\begin{array}{c}\hfill \frac{\text{d}P}{\text{d}t}=\frac{\rho -\omega }{\ell }P\hfill \end{array}$$(16)where ℓ is the mean prompt neutron lifetime, ρ the core reactivity and ω the fraction of delayed neutrons.

The increment of power heats the fuel, which brings negative feedback to the reactivity by the fuel temperature coefficient α $$\rho \left(t\right)-\omega =r-\alpha T\left(t\right)$$(17)where r = ρ (0) − ω. The physical interpretation of r is the initial excess of reactivity compared to the fraction of delayed neutrons. Concerning the temperature of the fuel, we can write that$$\begin{array}{c}\hfill Q\left(t\right)=m_{0}cT\left(t\right)={\int }_0^{t}P\left(s\right)\text{d}s\hfill \end{array}$$(18)where Q is the heat released, m ₀ is the mass of fuel and c is its heat capacity.

Following [5], we note $K=\frac{\alpha }{m_{0}c}$ and Q primitive of P, and $G\left(t\right)=Q\left(t\right)-\frac{r}{K}$.

Then equation (16) can be rewritten as$$\frac{d^{2}G}{\mathrm{d}{t}^2}=-\frac{K}{\ell }\cdot G\frac{dG}{\mathrm{d}t}.$$

By integration, there exists a constant ζ such that$$\begin{array}{c}\hfill G^{\prime }\left(t\right)=\frac{K}{2\ell }\left({\zeta }^2-G^2\left(t\right)\right)\hfill \end{array}.$$(19)

One can see that, if we choose b = ζ and η = Kζ/2ℓ, the function$$H\left(t\right)=b\tanh \left(\eta t\right)$$

is a solution to the differential equation (19).

By using the initial conditions of our system [5] we can derive that$${\zeta }^2=\frac{2\ell }{K}{P}_0+\frac{r^2}{K^2}.$$

Hence, we have obtained the analytical solution of our system. In particular the maximum value for $P_m=\frac{K{\zeta }^2}{2\ell }=P_0+\frac{r^2}{2K\ell }$ And the total energy released is$$Q_{\infty }={\displaystyle \underset{t\to +\infty }{\lim }}Q\left(t\right)={\displaystyle \underset{t\to +\infty }{\lim }}G\left(t\right)+\frac{r}{K}=\zeta +\frac{r}{K}.$$

Finally, we can apply the data to estimate the energy released.

In the RBMK reactor, the neutron mean lifetime is mainly determined by the diffusion time, which is about the order of 10⁻²s  [5] so we could choose ℓ = 0.01 s in our model. For the initial power, we choose P ₀ = 1000 MW for the reason that the operators started the test when the power was 200 MW and after a rapid increment of power they actuated the emergency shutdown system. The basic extra reactivity is caused by the void effect which is equivalent to 5  $    = 5ω of reactivity [7] so we choose r = 4 $ = 2400 pcm. The heat capacity of the fuel UO₂ is c = 216J ⋅ kg⁻¹ ⋅ K⁻¹ and the total mass of fuel in the core is m ₀ = 190 t. [2] As for the fuel temperature coefficient, we have α = 1.2 × 10⁻⁵ K⁻¹ according to the INSAG-7 report. Hence the total energy released is Q _∞ = 165 GJ and the maximum power is P _m  = 99 179 MW and it is reached at t _m  = 2.48s.

Some researchers have also estimated the energy released of the Chernobyl accident. For example, Pakhomov-Dubasov [12] estimated that the energy released is equivalent to 10t TNT, which means 42 GJ. And Malko obtained a result of 200t equivalent TNT, which means 837 GJ [13]. Our evaluation of the energy is closer to the first one, but there is still a significant difference. We can try to analyze the reason which may cause the disparity. First, the former research considers the impact of the reactor poison. This poison can effectively absorb the thermal neutrons, so it could prevent the increment of power, which results in less energy released. As for the latter research, it does not give more details about the estimation, but our model only estimates the energy released by the reactivity accident. In fact, there may be some energy released in another way, as the explosion of the pressure vessel in the core.

9 Conclusion

The subject of our work is to use simplified models to better understand the Chernobyl accident. We had the following five targets.

• To show that, from its very design, an RBMK is more unstable than a PWR: to keep constant power is not an easy matter and requires an automatic regulation. The initiating event, which leads to the accident, is when the operating staff decided to switch from automatic to manual control.

• To show that, in some situations, the design of emergency control rods can lead to reactivity (rather than antireactivity) injection when the emergency shutdown system is activated.

• To show that if the reactivity injected in the core is too high, this can produce a strong shock wave which will propagate in the fuel channels.

• Such a shock wave will eventually break at least half the fuel channels which induces depressurization of the coolant. Our target was to show that this phase takes a few seconds.

• To evaluate the fission energy released by the reactivity accident induced by the void effect.

Concerning our work, we have done the following researches.

• By using TRIPOLI-4 (R), we have evaluated the migration area M ². This is a contribution to target 1 since it is well known that the smaller M/H (where H is the height of the core) the more significant the instability is in the case of xenon poisoning. We have also shown that when there is xenon poisoning in the upper half of the core, the insertion of emergency absorbing bars leads to a significant reactivity increase.

• We have developed a one-dimensional axial model of the RBMK core showing its low stability and a strategy for controlling the overall reactivity of the core. Unfortunately, even with such a control the power is still volatile.

• As for the thermal-hydraulic properties of RBMK, we have taken two processes during the accident into account. We have developed a one-dimensional Lagrange model of the steam explosion, in which we have studied velocity of the shock wave. We have considered the stiffened gas as an approximation of the real situation. Concerning the following depressurization process, we have developed as well a one-dimensional model based on the Wilkins' scheme to simulate its duration.

• We have used the Nordheim-Fuchs model, which is a point reactor model, to evaluate the fission energy released by the reactivity accident.

Author contribution statement

Ziyue ZHUANG has implemented the TRIPOLI models to evaluate both the migration area and the reactivity induced in the core by the insertion of control rods with or without xenon poisoning.

Bertrand MERCIER has provided the general idea for this work, the equations and the numerical methods for the xenon poisoning, for the automatic control of the rods and for the hydrodynamic simulations.

Di YANG has implemented and tuned the automatic control of the rods with iodine and xenon evolution. He has also made the evaluation of the energy released by the reactivity accident.

Jiajie LIANG has implemented the hydrodynamic models, both for the steam explosion and for the depressurization phase.

References

All Tables

All Figures

	Fig. 1 Schematic view of the core (from Wikipedia).
In the text

	Fig. 2 Unit representative of the core.
In the text

	Fig. 3 Sectional view of FC.
In the text

	Fig. 4 Sectional view of CPS with the absorber.
In the text

	Fig. 5 Sectional view of CPS with the displacer.
In the text

	Fig. 6 Control rod falling process.
In the text

	Fig. 7 Control rods are raised to the upper limit without poisoning.
In the text

	Fig. 8 Control rods are initially inserted without poisoning.
In the text

	Fig. 9 Control rods are raised to the upper limit with poisoning.
In the text

	Fig. 10 Control rods are initially inserted with poisoning.
In the text

	Fig. 11 Evolution of the average thermal flux (unsymmetric control).
In the text

	Fig. 12 Evolution of the thermal flux distribution. The colour code is from green (initially) to red (final time).
In the text

	Fig. 13 Evolution of the average thermal flux (“symmetric” control).
In the text

	Fig. 14 “Symmetric” control. The colour code is from green (initially) to red (final time).
In the text

	Fig. 15 Steam explosion: pressure profiles in the pressure tube.
In the text

	Fig. 16 Cooling system of the RBMK reactor.
In the text

	Fig. 17 Evolution of volumic mass along the 24 m pressure tubes.
In the text