A simplified analysis of the Chernobyl accident

Bertrand Mercier, Di Yang, Ziyue Zhuang, Jiajie Liang

To cite this version:

Bertrand Mercier, Di Yang, Ziyue Zhuang, Jiajie Liang. A simplified analysis of the Chernobyl
accident. EPJ N - Nuclear Sciences & Technologies, 2021, 7, pp.1. �10.1051/epjn/2020021�. �hal-
03117177�

HAL Id: hal-03117177

https://hal.science/hal-03117177v1
Submitted on 20 Jan 2021

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
EPJ Nuclear Sci. Technol. 7, 1 (2021) Nuclear
Sciences
© B. Mercier et al., published by EDP Sciences, 2021 & Technologies
https://doi.org/10.1051/epjn/2020021
Available online at:
https://www.epj-n.org

REGULAR ARTICLE

A simplified analysis of the Chernobyl accident

Bertrand Mercier1,*, Di Yang2, Ziyue Zhuang2, and Jiajie Liang2
1
CEA/INSTN, Gif sur Yvette 91191, France
2
Institut Franco-Chinois de l’Energie Nucléaire, Sun Yat-sen University, 519082 Zhuhai Campus, Guangdong Province,
P.R. China

Received: 13 July 2020 / Received in final form: 3 December 2020 / Accepted: 21 December 2020

Abstract. We show with simplified numerical models, that for the kind of RBMK operated in Chernobyl:
– The core was unstable due to its large size and to its weak power counter-reaction coefficient, so that the power
of the reactor was not easy to control even with an automatic system.
– Xenon oscillations could easily be activated.
– When there was xenon poisoning in the upper half of the core, the safety rods were designed in such a way that,
at least initially, they were increasing (and not decreasing) the core reactivity.
– This reactivity increase has been sufficient to lead to a very high pressure increase in a significant amount of
liquid water in the fuel channels thus inducing a strong propagating shock wave leading to a failure of half the
pressure tubes at their junction with the drum separators.
– The depressurization phase (flash evaporation) following this failure has produced, after one second, a
significant decrease of the water density in half the pressure tubes and then a strong reactivity accident due to
the positive void effect reactivity coefficient.
– We evaluate the fission energy released by the accident

1 Introduction The RBMK is a kind of thermal-neutron reactor

designed and built by the former Soviet-Union. It is an
A better understanding of former nuclear accidents is early Generation II reactor that dates back to 1954, some of
obviously good for nuclear safety. which are still in operation in ex-USSR. In the RBMK
Many countries are developing nuclear reactor projects, reactor core, the fast neutrons generated from the fission of
235
and it is useful for them (and more generally for the world) U are moderated by graphite so that they can easily
to benefit from an adequate feedback from past events. actuate other fissions. Water is used as a coolant to transfer
For teaching purpose, also, it is useful to show which the heat from inside the core to outside. From a neutronics
conclusions can be obtained with simple models. point of view, light water is both an absorber and a
In the INSAG-7 report on the Chernobyl accident [1], moderator. Since there is enough graphite to ensure
the following design features of the RBMK have been moderation, the contribution of water to moderation is
identified as having a major responsibility: negligible, but this is not the case for absorption. This
– void coefficient of reactivity explains the positive void coefficient of reactivity which has
– design and control of safety rods been evaluated before by e.g. Kalashnikov [2].
– speed of insertion of the emergency control rods In Section 2 we shall use TRIPOLI 4, which imple-
– control of power ments the Monte-Carlo method, to evaluate the safety
rods efficiency. The Monte-Carlo method is not a
However, this does not explain why, when the operators simplified method, but we shall apply it on simplified
pressed the AZ-5 button to command the reactor scram, geometries with a network of 3
3 fuel elements and
nothing happened during 3 s and then the reactor exploded lateral reflective boundary conditions as if the diameter of
(see e.g. https://www.neimagazine.com/features/feature the core was infinitely large, or the radial reflectors fully
how-it-was-an-operator-s-perspective/). efficient.
For the safety rod efficiency we find like in INSAG-7
that they may induce “a local insertion of positive reactivity
* e-mail: Bertrand.mercier@cea.fr in the lower part of the core.” And that “the magnitude of

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0),
which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
2 B. Mercier et al.: EPJ Nuclear Sci. Technol. 7, 1 (2021)

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 2nd 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 Fig. 1. Schematic view of the core (from Wikipedia).
(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 211 control bars are represented: 32 are short and
the failure of the pressure tubes at their junction with the inserted from below whereas 179 are long and inserted from
drum separators. Then this so-called steam explosion will above.
be followed by a depressurization process in the 24 m long The division of those 179 control rods between
pressure tubes. Such a depressurization is also called automatic and emergency protection groups was arbitrary;
flash evaporation. This is an isentropic process with on Figure 1 the division is 167 + 12.
liquid to steam phase change. We also use a 1D When control rods are inserted into the core, they
hydrodynamic model but we take into account the should absorb lots of thermal neutrons effectively so that
conservation of entropy to make the computational the chain reaction would be weaker, and the reactivity
model simpler. would be reduced accordingly. However, a graphite rod
After the depressurization process is over, the water whose length is 4.5 m called “displacer” was attached to the
density in the pressure tubes will decrease a lot, so that the end of the rod of B4C. The graphite has a much smaller
positive void coefficient of reactivity produces a strong cross section of absorption than B4C, and even than water.
reactivity accident. In Section 8 we evaluate the fission This absorber-displacer design has two advantages. One is
energy released by such a reactivity accident. to give more control range because when the displacer is in
We use the Nordheim-Fuchs model [5,6], which is a 0D the core, it could give more reactivity, and the absorber has
model and gives 165 GJ. the opposite effect. The other is to prevent coolant water
This is not the total energy released by the explosion, from entering the space vacated when the absorber is
since we should also add the energy released through the withdrawn, which enables a smoother thermal flux in the
depressurization process. 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
2 Neutron transport calculation: simplified the core and periods where it is higher in the top of the core.
model In the former case, after over-extraction of the control rod,
when it is re-inserted into the core, the graphite tip would
Figure 1 extracted from Wikipedia represents a schematic induce a lot of reactivity in the lower part of the core and
view of the core of the Chernobyl-4 RBMK. The radial thus increase the power significantly. No doubt that this is
reflector (made of 4 rows of pure graphite blocks) is not an essential contribution to the accident. (see GRS 121 [7],
represented. p. 44).
 B. Mercier et al.: EPJ Nuclear Sci. Technol. 7, 1 (2021) 3

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

Fig. 2. Unit representative of the core.

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

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
Fig. 3. Sectional view of FC. 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
To check that this is the case, we shall use a Monte- part, in the upper part of the control rod (see Fig. 4).
Carlo code (Tripoli-4R) but on a simplified geometry. In However, in the lower part, there is successively water
view of Figure 1, it seems appropriate to consider a 3
3 (1.25 m), a graphite extension (4.5 m) and again water
network containing two kinds of components which are the (1.25 m). When the rods are in a high position, their
fuel channel (FC) and channel of control protection system insertion begins by replacing the water in the lower part of
(CPS), both of which are axial models. Our 3
3 model will the channel with graphite.
involve 7 FC and 2 CPS (see Fig. 2). More precisely, Figure 5 represents a sectional view of
Our approach is similar to the one used in Kalashnikov the graphite displacer when it is located in the core.
[2] who considers also a 3
3 network but with 6 FC, 2 CPS An axial view of the control rod falling process is shown
and 1 AA (Additional Absorber). In fact, his target was to on Figure 6.
prove that, to reduce the void effect on a RBMK, it is a This process will bring about changes in reactor
good idea to both increase the fuel enrichment and reactivity. We know (see INSAG-7) that there was xenon
introduce additional absorber. poisoning in the reactor at the time of the Chernobyl
Our target is to identify the situations where the accident. In this section, we will study the effects of the
insertion of control rods can increase (and not reduce) insertion of the control rod on the reactivity of the reactor
reactivity. in two situations: normal operation and the presence of
The sectional view of a fuel channel is shown in Figure 3. xenon poisoning.
FC consists of the fuel assembly and the pressure tube In order to study the influence of the control rod
surrounded by the graphite. The assembly includes 18 fuel insertion under different conditions, we shall use two
rods, which are arranged in the pressure tube in two models to simulate the states when the control rods are
concentric rings: the internal ring is formed by 6 rods raised to the upper limit and when the control rods are
(internal rods) and the outer ring is formed by 12 rods initially inserted. For a reactor under normal operating
(outer rods). conditions, the schematic diagrams of the two states are
4 B. Mercier et al.: EPJ Nuclear Sci. Technol. 7, 1 (2021)

Fig. 6. Control rod falling process.

Fig. 8. Control rods are initially inserted without poisoning.

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

shown in Figures 7 and 8. Moreover, we also consider the

xenon poisoning, which brings a heavy accumulation of
135
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
1015
atom of 135Xe 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 Fig. 10. Control rods are initially inserted with poisoning.
 B. Mercier et al.: EPJ Nuclear Sci. Technol. 7, 1 (2021) 5

Table 1. Results of the insertion effect of control rods a change of flow rate in the pressure tubes (see [9]). In
(k∞ values). normal operations, film boiling (DNB), and then significant
variations of the steam quality x, should be avoided. At full
CPS state Reactor state 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
No xenon poisoning Xenon poisoning the flow rate tuning at intermediate powers, but we can
Raised to the 1:2509 ± 2:2
104 1:2370 ± 2:1
104 assume that the void effect contribution to the power
upper limit reactivity coefficient was limited in normal operations.
As specified in [8], in some (incidental) situations, the
Inserted 1:2466 ± 2:1
104 1:2419 ± 2:2
104
power coefficient could be positive so that a power increase
could cause a power excursion.
As explained in WIKIPEDIA, the inlet temperature in
in a state of xenon poisoning, they locally increase the the pressure tubes results from a mixture between
reactivity, which led to a surge in energy and triggered the feedwater, which is at a low temperature, and recirculating
first explosion. water which is at saturation temperature. At low power,
the feedwater flow rate is low so that the inlet temperature
2.1 Evaluation of the migration area is “dangerously close to the saturation temperature.”
As indicated in INSAG-7 [1], p. 105, in our study, we
We know that (see e.g. [5]), when we use the one group 1D shall consider that the power coefficient is weakly negative,
axial model, we have and this is sufficient, as we shall see, for the reactor to be

 difficult to control.
keff ¼ k∞ = 1 þ l0 M 2 ð1Þ
4 Safety test
where l0 = p2/H2 is the smallest solution of the eigenvalue
problem:
The detailed sequence leading to the accident has been
described in several papers like INSAG-7 [1] or GRS121 [7].
d F
2
 ¼ l0 F; 0 < z < H; Fð0Þ ¼ FðHÞ ¼ 0: There are 2
4 = 8 main circulation pumps on the
dz2 RBMK. At full power 2
3 = 6 are expected to deliver
8000 m3/h each. As indicated by INSAG-7 [1], p. 8, fewer
To evaluate the migration area M2 it is appropriate to
pumps are used during normal start-up or shutdown
perform two Monte-Carlo calculations : one with the
operations. The specific test which was planned required
leakage boundary conditions in z = 0 and z = H, which will
2
4 = 8 pumps to run simultaneously. However for some
give keff, and the other one with fully reflective boundary
reasons, the reactor power was much lower than expected,
conditions which will give k∞. Then we use (1). We find
so that the feedwater flow rate (coming from the condenser)
M2 = 190 cm2.
was low and the inlet temperature in the core was close to
The ratio M/H is an indicator to measure the stability
the saturation temperature of the recirculating liquid
of a nuclear reactor. According to reference [3], M/H in
water. So, when 4 pumps out of 8, as expected by the safety
PWR is approximately 0.025, while M/H in the RBMK
test, suddenly stopped, there has been a sudden jump in
reactor obtained here is approximately 0.020, this means
reactivity due to the void effect. This may have also been
that in terms of xenon oscillation, the RBMK core is more
the initiating event.
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
5 Basic equations describing the xenon
the power decline, they raised too many absorber rods. poisoning
This move put the upper part of the core under xenon
poisoning. It can be seen from the above conclusion that For simplicity, we shall consider 1D axial oscillations
the instability of the RBMK core has been one reason for only.
the accident. Like in [3] we shall use the following one group model:

dY P
3 Stability and reactivity control ¼g F  lY ð2Þ
dt f

It is important to know what happens to the reactivity of

the core when the power of the reactor increases. In the dX
¼ lY  ðm þ s FÞX ð3Þ
book by Barré et al. [8], p. 55 it is explained that, for dt
RBMK, there are at least 3 contributions to the power
reactivity coefficient (a) the fuel temperature coefficient M 2 F00 þ ð1 þ aX þ eÞF ¼ k∞ F; 0 < z < H ð4Þ
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 Fð0Þ ¼ FðH Þ ¼ 0 ð5Þ
6 B. Mercier et al.: EPJ Nuclear Sci. Technol. 7, 1 (2021)

Table 2. Typical values for the parameters introduced in Table 3. Adjustment rules of the position of control rod.
our model.
p Dm p Dm
l 2:92
105 s1 v
P
0.0317
p < 0.95 +5 1.01  p < 1.02 1
m 2:12
105 s1 0.24 cm−1
Pf 0.95  p < 0.97 +3 1.02 < p  1.03 2
g 0.064 0.0129 cm−1
0.97  p < 0.98 +2 1.03 < p  1.05 3
s 2.106 barns
0.98  p < 0.99 +1 p > 1.05 5

where:

– Y (resp. X) is the concentration of 135 135
53 I resp: 54 Xe in of (long) control rods inserted from above that were used
the fuel evaluated in at/cm3 of fuel for automatic control. In what follows, we shall consider
– s is the microscopic capture cross section for xenon in the two types of control:
thermal
P domain (evaluated in barns) – non symmetric control where only control rods inserted
– f is the macroscopic fission cross section of the fuel in from above are used for automatic control.
the thermal domain – “Symmetric” control where both control rods inserted
– g the iodine fission yield (about equal to .064 for UO2 from above and control rods inserted from below are used
fuel) for automatic control. They move from the same number
– l (resp.
m) 135
the time
 constant for the b  decay of of steps at each time, however, we give them different
135 weight.
53 I resp: 54 Xe
– F = F (z, t) is the thermal neutron flux In the first case, we shall have:
– M is the
2
Pmigration area
– a = vs/ (
P , where v is the fuel volume fraction in the core 0; 0 < z  z0
and the absorption macroscopic cross section of the eðzÞ ¼
homogeneous core without xenon poisoning e0 ; z 0 < z < H
– e = e (z) is the additional absorption induced by the
insertion of control bars. whereas in the second one (we let z1 = H  z0) we shall
have:
We shall assume that k∞ ¼ k0∞ : ð1
 b FÞ and b is the
prompt power feedback coefficient. It takes the Doppler 8
effect into account but also, in a simplified way, the < e1 ; z  z 0
>
moderator effect and the void effect. For stability it is eðzÞ ¼ 0; z0 < z  z1 :
>
:
necessary that b > 0. e2 ; z > z1
When X = X (z) and e = e (z) are given, we notice that
equations (4) and (5) is an eigenvalue problem which gives The principle of our control is to adjust z0, that is, the
the infinite multiplication factor k∞ at zero power. insertion of controls rods, to stabilize the total power.
To obtain a non zero power, we must select k0∞ > k∞ However, for practical reasons, we shall assume that z0 can
and the core power will depend on k0∞  k∞ . only take discretized values, namely an integer number
The system (2)–(5) is a non linear system of differential m ∈ {0, 70} times 10 cm (note that for the core of a RBMK,
equations with a constraint. H = 7 m).
Both an implicit and a semi-implicit scheme are The total power is proportional to the average thermal
introduced in [3] for its solution. flux in our equations. Assume that at t = t n we have
In what follows, we shall use the semi-implicit scheme obtained (Fn) (Xn) and (Yn), and the position of control
and the values given in Table 2 for the parameters of our rods at the last time t = tn1 is mn1.
model. First, let p = ave (Fn)/ave (F0) we adjust m according
to p which is the ratio between the average flux at t = tn
and at t = 0. Then we make mn = mn1 + Dm by defining
6 Simulation of the automatic control Dm as shown in Table 3 above.
of the rods The second step of adjustment is to verify if the control
rod is in an extreme position. When the control rod is
Here is what INSAG-7 says “When the reactor was almost totally extracted, which means that m(n1) + Dm is
operated at low power (…) the operator had mainly to small, the reactivity is insufficient for the core. Hence, we
rely on the (…) ex-core detectors. However (…) they could propose an increment Dk∞ to the reactivity as if we
not indicate the average axial distribution of flux, since extracted totally another control rod, then we reset m to a
they were all located in the mid-plane of the core.” To small value which means we insert the control rod deeply
simulate the instability of the core at low power, we have into the core. Moreover, when the control rod is almost
introduced an appropriate feed-back law for the control totally inserted, we propose a similar method. The details
rods. of our model are shown in Table 4.
In Section 2 we have specified that there were 32 short After these two steps of adjustment, we obtain m(n) for
control rods inserted from below and an arbitrary number resolution in the next time step.
 B. Mercier et al.: EPJ Nuclear Sci. Technol. 7, 1 (2021) 7

Table 4. Adjustment rules of the control rod at extreme

positions.

m(n1) Dk∞ m(n)

m(n1) + Dm < 4 +0.01 64
4  m(n1) + Dm  65 0 m(n1) + Dm
m(n1) + Dm > 65 0.01 5

Fig. 12. Evolution of the thermal flux distribution. The colour

code is from green (initially) to red (final time).

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
Fig. 11. Evolution of the average thermal flux (unsymmetric possibly make the situation worse. They tried to control the
control). 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
6.1 Simulation with non-symmetric control control.

We have simulated the temporal evolution of the core 6.2 Simulation with “symmetric” control
during 40 h starting from the initial values, with e0 = 0.01.
Figure 11 −shows the evolution of the average thermal flux We have selected e1 = 0.002 and e2 = 0.005, in other words
in the core, which also represents the power output of the the weight of the control rods inserted from below is 2.5
reactor. We can see that at the beginning, the average flux times less than the weight of the control rods inserted from
has a sudden decrease and a following sudden increment. above. However it is sufficient to give some stability as
That is the consequence of the insertion of the control rods. shown on Figures 13 and 14.
The average flux starts to oscillate in a small amplitude but Using both control rods inserted from below and control
with a high frequency. This relatively stable state lasts for rods inserted from above is then a key factor to stabilize the
the first 17 h. Remarkably, between t = 7 h and t = 8.5 h the core power. Even so, we see that the frequency of control
average flux evolves smoothly. Then at t = 17 h the rods tuning is quite high and this can be achieved with an
oscillation becomes much more drastic than before. That electronic control only.
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
7 Hydrodynamic simulation of the coolant
to time. We see that a xenon oscillation develops: before in the pressure tubes
t = 17 h, the power is higher in the lower part. After, this is
just the reverse. In this section we try to explain why, when the operators
But a one step insertion of the control rods from above pressed the AZ-5 button to command the reactor scram,
(i-e 10 cm) has much more effect on the reactivity when the nothing happened for a (short) while.
power is higher in the upper part. The reactor being unstable, its power was difficult to
During our whole simulation, the maximum and control. Moreover, when 4 of the 8 recirculating pumps
minimum of the average flux are 3.21
1013n/(cm2 ⋅ s) were stopped, the (limited) void effect induced a reactivity
and 1.34
1013n/(cm2 ⋅ s) compared to the initial average increase so that the operators decided to command the
flux which was equal to 2.30
1013n/(cm2 ⋅ s). reactor scram.
8 B. Mercier et al.: EPJ Nuclear Sci. Technol. 7, 1 (2021)

e ∼ 1310kJ/kg. Then the instantaneous energy deposition

increases e but cannot change r 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
Fig. 13. Evolution of the average thermal flux (“symmetric” between the liquid and the steam phase, so that we can use
control). a stiffened gas equation of state for both phases.

p ¼ gp∞ þ ðg  1Þðe  qÞr: ð6Þ

For the liquid water, we select (see Corot [10] Tab. 6.1):

g ¼ 2:35; p∞ ¼ 1:E9 P a; q ¼ 1167E3 J=kg:

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

∂t ∂u
 ¼0 ð7Þ
∂t ∂a

∂u ∂p
þ ¼0 ð8Þ
∂t ∂a

∂E ∂
þ ðpuÞ ¼ 0 ð9Þ
∂t ∂a
Fig. 14. “Symmetric” control. The colour code is from green
(initially) to red (final time). where t = 1/r and E = e + u2/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
As we explained before, due to xenon poisoning in the the Lagrange coordinate a such that
upper part of the core, during the first second, the insertion
of safety bars still increased (rather than decreased) the . . . ai2 < ai1 < ai < aiþ1 < aiþ2 . . . :
reactivity of the core.
We can assume that, at this point there was no prompt We also introduce a time discretization as follows
criticality yet, but that the fission power increase has been
sufficient to transfer a significant amount of energy to the 0 < Dt < 2Dt . . . < nDt < ðn þ 1ÞDt < . . .
liquid coolant which has been instantaneously transformed Integrating equations (7)–(9) for ai  a  ai+1 and
in (very) high pressure state at the origin of a strong nDt  t  (n + 1) Dt and we get
propagating shock in the pressure tubes. Assume that the
pressure of the coolant was p = 8 MPa and its volumic nþ1 nþ1
tnþ1
iþ1=2  t iþ1=2
n
u 2  ui 2
mass was r = 722 kg/m3 (i.e. the volumic mass of the  iþ1 ¼0 ð10Þ
liquid at saturation) so that its specific internal energy is Dt aiþ1  ai
 B. Mercier et al.: EPJ Nuclear Sci. Technol. 7, 1 (2021) 9

nþ1 nþ1
iþ1=2  uiþ1=2
unþ1 n
p 2  pi 2
þ iþ1 ¼0 ð11Þ
Dt aiþ1  ai

nþ1 nþ12
iþ 1=2  E iþ1=2
Enþ1 n
ðpuÞiþ12  ðpuÞi
þ ¼0 ð12Þ
Dt aiþ1  ai

where t niþ1=2 ðresp: uniþ1=2 Þ denotes the (assumed constant)

value of t (resp. u) at time 1nDt on the cell ai  a  ai+1.
nþ1 nþ
To evaluate ui 2 or pi 2 we use the fact that for the
(linearized) acoustic equations, p  au where a ≡ r . c (and
c is the sound speed) is a Riemann invariant propagating
from right to left. Similarly p + au propagates from left to Fig. 15. Steam explosion: pressure profiles in the pressure tube.
right.
Let
shall consider a 6 m long tube with fixed pressure (8 MPa)
 nþ1 in z = 0 and prescribed velocity (u = 0) in z = 6m (z is then
uL ¼ uniþ1=2 ; uR ¼ uni1=2 and u ¼ ui 2
the abscissa along the tube).
At t = 0 we shall start from r (z, 0) = 722kg/m3 and:
pL ¼ pniþ1=2 ; pR ¼ pni1=2 and p ¼ pi
nþ12
: (
8 MPa; 0  z  5:1m
pðz; 0Þ ¼
We shall require u*, p* to be the solution of the 2
2 200 MPa; 5:1m < z  6m:
system
We take 600 cells each with Dz = 1 cm and we obtain the
p þ aL u ¼ pL þ aL uL pressure profiles at t = 0.01ms and t = 0.02ms which are
represented on Figure 15 and show that the pressure wave
p  aR u ¼ pR  aR uR : is propagating without attenuation at 2000 m/s.
We check that the pressure decreases behind the shock
We proceed similarly for the boundary conditions. wave but this is not sufficient to obtain a phase change from
iþ1=2 ; uiþ1=2 and
From equations (10)–(12) we obtain tnþ1 nþ1
liquid to steam. The numerical results show that the
nþ1
Eiþ1=2 . volumic mass remains everywhere above 700 kg/m3.
We then compute the new internal energy On the other hand, we have a strong shock wave
(100 MPa) which is sufficient to break the pressure tubes
 2 outside the core, namely at the junction with the drum
iþ1=2 ¼ E iþ1=2  uiþ1=2 =2
enþ1 nþ1 nþ1
separators.
Remark: assuming that the energy deposition takes
and we use the equation of state (6) to obtain the new place in the liquid water is conservative.
iþ1=2 at time (n + 1) Dt
pressure pnþ1 A simple calculation shows that if it takes place in a
In the Chernobyl accident, we consider that there was 1 two-phase mixture liquid-with 14% steam, in the same
GJ of energy rapidly released into the coolant. Notice that 7.1 m3 volume, then the pressure will climb to 35 MPa
such an energy corresponds to an increase for 1 s of 1 GW of rather than 200 MPa. The fact is that steam is much more
the reactor power (that is the quarter of the nominal compressible than liquid water.
power). We know that the mass of coolant contained in the
pressure tubes is about 32t. 7.2 Simulation of the depressurization phase
We shall assume that the 1 GJ energy deposition
increases e from 1310 to1510 kJ/kg for 5t of liquid water. After the failure of the pressure tubes, the pressure at the
From the equation of state this will increase its pressure open end of the pressure tube suddenly decreases from 8 to
from 8 to 200 MPa. 0.1 MPa. The depressurization activates a phase change
The section of each pressure tube is equal to 50 cm2. from liquid to steam which is an isentropic process. We also
However the coolant can flow in a section reduced to use a 1D hydrodynamic model but we take into account the
22.5 cm2 Taking into account the number of fuel channels conservation of entropy to make the computational model
(1661) we can assume that the 1 GJ energy is deposited on simpler.
a limited height of 1.8 m that is about 25% of the core After the depressurization process is over, the water
height. density in the pressure tubes will decrease a lot, so that the
We choose a simple model just for this basic study. The positive void coefficient of reactivity produces a strong
length of the fuel channel is 7 m in the core but we extend it reactivity injection and then prompt-criticity.
to 12 m, by including out of core parts. We assume that the To simulate the depressurization, we shall use the
energy is deposited in its central part: by symmetry we Wilkins’ scheme (see Ref. in DESPRES [11]) and consider
10 B. Mercier et al.: EPJ Nuclear Sci. Technol. 7, 1 (2021)

the conservation of mass and momentum only.

8
>
> ∂t ∂u
<  ¼0
∂t ∂a
:
>
> ∂u ∂p
: þ ¼0
∂t ∂a
As before, we set the grid of the variable a like the
following

⋯ ai2 < ai1 < ai < aiþ1 < aiþ2 ⋯

Then the middle of an element can be denoted by the

subscript i þ 12

ai < aiþ1 < aiþ1 :

2

Besides, a discretization of time in these instants is Fig. 16. Cooling system of the RBMK reactor.

0 < Dt < 2Dt ⋯ < nDt < ðn þ 1ÞDt ⋯

After 1 s, the volumic mass has decreased to r =
A staggered grid is used in the Wilkins’ scheme which 105 kg/m3 which is sufficient for the reactivity accident to
writes start and to break the other remaining pressure tubes.
tn  tn11 n
1 1
n 2
1
iþ 2 iþ 2 u 2  ui ð13Þ 8 Evaluation of the energy released caused
 iþ1 ¼0
Dt aiþ1  ai by a reactivity accident
1
nþ 2
1
n 2 pn  pn We shall use the Nordheim-Fuchs model [5]. Let P denote
ui  ui 1
iþ 2
1
i 2
the fission power:
þ ¼ 0: ð14Þ
Dt aiþ 1  ai 1
2 2 dP r  v
¼ P ð16Þ
We denote the position of the grid ai at time nDt as zni . dt ℓ
With the Lagrangian description, we know that the new
where ℓ is the mean prompt neutron lifetime, r the core
position of this grid is given by
reactivity and v the fraction of delayed neutrons.
1 The increment of power heats the fuel, which brings
znþ1 ¼ zni þ Dtui
nþ 2 : ð15Þ negative feedback to the reactivity by the fuel temperature
i
coefficient a
Furthermore, tniþ1 being known we evaluate pniþ1 by rðtÞ  v ¼ r  aT ðtÞ ð17Þ
2 2
using the isentropic process.
We assume that the pressure tubes will fail at the where r = r (0)  v. The physical interpretation of r is the
junction with the drum separator (Fig. 16). initial excess of reactivity compared to the fraction of
We estimate that the length is 24 m, including a length delayed neutrons. Concerning the temperature of the fuel,
of 7 m of the channel in the reactor core. According to the we can write that
construction of the cooling system (Fig. 16), the whole Z t
channel we study is not straight. To simplify our simulation QðtÞ ¼ m0 cT ðtÞ ¼ P ðsÞds ð18Þ
work, we still consider that the model is one-dimensional. 0
As we have seen in Section 4 we can assume that half the
pressure tubes will be filled of saturated liquid water. where Q is the heat released, m0 is the mass of fuel and c is
That’s where the shock due to the fission energy increase its heat capacity.
will be the stronger. So these pressure tubes will eventually Following [5], we note K ¼ ma0 c and Q primitive of P,
depressurize as follows. The initial distribution of the and GðtÞ ¼ QðtÞ  Kr .
volumic mass is Then equation (16) can be rewritten as
d2 G K dG
rðz; 0Þ ¼ 705 kg=m3 ; 0  z  24m: ¼  ⋅G :
dt 2 ℓ dt
The pressure is also constant and equal to 7.5 MPa. By integration, there exists a constant z such that
Using the same boundary conditions at both ends we
get the results in Figure 17 where we represent the tube K
2 
G0 ðtÞ ¼ z  G2 ðtÞ : ð19Þ
(24 m) on the right and 20 m of the outside world on the left. 2ℓ
 B. Mercier et al.: EPJ Nuclear Sci. Technol. 7, 1 (2021) 11

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

One can see that, if we choose b = z and h = Kz/2ℓ, the Some researchers have also estimated the energy
function released of the Chernobyl accident. For example,
H ðtÞ ¼ b tanhðhtÞ Pakhomov-Dubasov [12] estimated that the energy
released is equivalent to 10t TNT, which means 42
is a solution to the differential equation (19). GJ. And Malko obtained a result of 200t equivalent
By using the initial conditions of our system [5] we can TNT, which means 837 GJ [13]. Our evaluation of the
derive that energy is closer to the first one, but there is still a
significant difference. We can try to analyze the reason
2ℓ r2
z2 ¼ P0 þ 2 : which may cause the disparity. First, the former research
K K considers the impact of the reactor poison. This poison
Hence, we have obtained the analytical solution of can effectively absorb the thermal neutrons, so it could
our system. In particular the maximum value for prevent the increment of power, which results in less
2
r2
energy released. As for the latter research, it does not
P m ¼ Kz
2ℓ ¼ P 0 þ 2Kℓ And the total energy released is give more details about the estimation, but our model
r r only estimates the energy released by the reactivity
Q∞ ¼ lim QðtÞ ¼ lim GðtÞ þ ¼zþ : accident. In fact, there may be some energy released in
t!þ∞ t!þ∞ K K another way, as the explosion of the pressure vessel in
Finally, we can apply the data to estimate the energy the core.
released.
In the RBMK reactor, the neutron mean lifetime is
mainly determined by the diffusion time, which is about the 9 Conclusion
order of 102s [5] so we could choose ℓ = 0.01 s in our
model. For the initial power, we choose P0 = 1000 MW for The subject of our work is to use simplified models to better
the reason that the operators started the test when the understand the Chernobyl accident. We had the following
power was 200 MW and after a rapid increment of power five targets.
they actuated the emergency shutdown system. The basic – To show that, from its very design, an RBMK is more
extra reactivity is caused by the void effect which is unstable than a PWR: to keep constant power is not an
equivalent to 5 $ = 5v of reactivity [7] so we choose easy matter and requires an automatic regulation. The
r = 4 $ = 2400 pcm. The heat capacity of the fuel UO2 is initiating event, which leads to the accident, is when the
c = 216J ⋅ kg1 ⋅ K1 and the total mass of fuel in the core is operating staff decided to switch from automatic to
m0 = 190 t. [2] As for the fuel temperature coefficient, we manual control.
have a = 1.2
105 K1 according to the INSAG-7 report. – To show that, in some situations, the design of emergency
Hence the total energy released is Q∞ = 165 GJ and the control rods can lead to reactivity (rather than
maximum power is Pm = 99 179 MW and it is reached at antireactivity) injection when the emergency shutdown
tm = 2.48s. system is activated.
12 B. Mercier et al.: EPJ Nuclear Sci. Technol. 7, 1 (2021)

– To show that if the reactivity injected in the core is too Di YANG has implemented and tuned the automatic
high, this can produce a strong shock wave which will control of the rods with iodine and xenon evolution. He has
propagate in the fuel channels. also made the evaluation of the energy released by the
– Such a shock wave will eventually break at least half the reactivity accident.
fuel channels which induces depressurization of the Jiajie LIANG has implemented the hydrodynamic
coolant. Our target was to show that this phase takes a models, both for the steam explosion and for the
few seconds. depressurization phase.
– To evaluate the fission energy released by the reactivity
accident induced by the void effect.
Concerning our work, we have done the following References
researches.
– By using TRIPOLI-4 (R), we have evaluated the 1. International Nuclear Safety Advisory Group, The Chernobyl
migration area M2. This is a contribution to target 1 accident: updating of INSAG-1 (International Atomic
since it is well known that the smaller M/H (where H is Energy Agency, Vienna, 1992)
the height of the core) the more significant the instability 2. D. Kalashnikov, APOLLO2 calculations of RBMK lattices.
is in the case of xenon poisoning. We have also shown that France, 1998:34
3. B. Mercier, Z. Ziliang, C. Liyi, S. Nuoya, Modeling and
when there is xenon poisoning in the upper half of the
control of xenon oscillations in thermal neutron reactors, EPJ
core, the insertion of emergency absorbing bars leads to a
Nucl. Sci. Technol. 6, 48 (2020)
significant reactivity increase. 4. T. Corot, B. Mercier, A new nodal solver for the two
– We have developed a one-dimensional axial model of the dimensional Lagrangian hydrodynamics, J. Comput. Phys.
RBMK core showing its low stability and a strategy for 353, 1–25 (2018)
controlling the overall reactivity of the core. Unfortu- 5. P. Reuss, Fission nucléaire, réaction en chaine et criticité
nately, even with such a control the power is still (EDP Sciences, 2016, 91940 - Les Ulis France)
volatile. 6. T.V. Holschuh, W.R. Marcum, Modified Fuchs-Nordheim
– As for the thermal-hydraulic properties of RBMK, we model for use in reactor pulse measurements, Ann. Nucl.
have taken two processes during the accident into Energy 116, 314–318 (2018)
account. We have developed a one-dimensional Lagrange 7. GRS-121, The accident and the Safety of RBMK-Reactors, 1996
model of the steam explosion, in which we have studied 8. B. Barré, P. Anzieu, R. Lenain, J.B. Thomas, Nuclear reactor
velocity of the shock wave. We have considered the systems: a technical historical and dynamic approach (EDP-
stiffened gas as an approximation of the real situation. Sciences, 2016, 91940 - Les Ulis France)
Concerning the following depressurization process, we 9. A. Kaliatka, E. Ušpuras, R. Urbonas, Verification and
have developed as well a one-dimensional model based on validation of RELAP5 code against Ignalina NPP and E-108
the Wilkins’ scheme to simulate its duration. test facility data // Trans. of International Meeting on
– We have used the Nordheim-Fuchs model, which is a “Best-Estimate” Methods in Nuclear Installation Safety
point reactor model, to evaluate the fission energy Analysis BE-2000, 16–18 November 2000. Washington, USA.
released by the reactivity accident. 27 p
10. T. Corot, Numerical simulation of shock waves in a bi-fluid
flow: application to steam explosion. Conservatoire national
Author contribution statement des arts et metiers  CNAM, 2017
11. B. Després, Numerical methods for Eulerian and Lagrangian
Ziyue ZHUANG has implemented the TRIPOLI models to conservation laws (Springer International Publishing, Cham,
evaluate both the migration area and the reactivity 2017)
induced in the core by the insertion of control rods with 12. S.A. Pakhomov, Y.V. Dubasov, Estimation of explosion
or without xenon poisoning. energy yield at Chernobyl NPP accident, Pure Appl.
Bertrand MERCIER has provided the general idea for Geophys. 167, 575–580 (2010)
this work, the equations and the numerical methods for the 13. M.V. Malko, The Chernobyl reactor Design features and
xenon poisoning, for the automatic control of the rods and reasons for accident 11–27 (2016). https://inis.iaea.org/
for the hydrodynamic simulations. search/search.aspx?orig_q=RN:48080457

Cite this article as: Bertrand Mercier, Di Yang, Ziyue Zhuang, Jiajie Liang, A simplified analysis of the Chernobyl accident, EPJ
Nuclear Sci. Technol. 7, 1 (2021)