A look at recriticality during meltdown, part 1

The issue of recriticality in the damaged reactors at Fukushima pops up every now and then (a few examples link1, link2, link3, link4). Perhaps it is worth taking a look at what recriticality means, how likely it is and what it would mean if the cores goe critical. These posts will contain some maths and give some insight into basic reactor physics. Despite what most people think it is actually quite easy as long as one can follow the solution of some simple differential equations.

We will look at two different cases, in the first case the core has melted completely and is as a molten puddle or bed of "gravel" at the bottom of the vessel. In the second case the fuel rods are still mostly geometrically intact while the control rods have melted. If I have energy I might throw in a section about criticality in spent fuel pools as well at the end. We start with the completely molten core because it is easier and highlights all the relevant physics.

What exactly is criticality?

Fission is a reaction whereby a incoming neutron hits a nucleus, the nucleus then has a certain probability (depending on the energy of the neutron, what nucleus it is etc) of splitting into two roughly equally large pieces and in the process emit 2-3 new neutrons. Those neutrons can in turn hit new nuclei that causes more fissioning and voila, we have a chain reaction. If we assume we have a system where nothing is happening and we send in a burst of neutrons, those neutrons, that we will call the first generation, will cause an initial amount of fission reactions that produce a second generation of neutrons which goes on to create a third generation etc. Criticality is simply defined as the ratio between a subsequent generation with the one preceding it, it is usually designated by the letter K.

 \large K = \frac{Neutrons\ in\ generation\ x+1}{Neutrons\ in\ generation\ x}

If K is less than 1 we will get a chain of reactions that ultimately ends after a few generations. The total number of neutrons, N, that we get from a burst of initial neutron N_{i} in such a case is given by

\large N = N_{i}+N_{i}K+N_{i}K^2+N_{i}K^3.....

Which is equal to (when K less than 1)

\large N = \frac{N_{i}}{1-K}

This kind of system is called subcritical. While you get a multiplication of neutrons you do not get a self sustaining chain reaction and in order to get energy you constantly have to supply neutrons from an external source (this is the basis of accelerator driven systems or fusion/fission hybrid systems).

If K = 1 then the series above do not converge, in essence we can get an infinite amount of neutrons out of the initial neutron. We also see from the definition of K that it means that every generation of neutrons is as large as the preceding generation, in such a situation we have a stable chain reaction with a constant power being produced. This state is called a critical state.

if K is larger than 1 every generation will be larger than the preceding generation and the growth is exponential. If we have a K of 1.2 for instance then after 10 generations the neutron population will have increased by a factor of 6, after 20 generations its almost 40 times larger. The growth is very rapid. The time between neutron generations in a nuclear bomb for instance is on the order of 10^{-7}. In 10 microseconds (100 generations), with a K of 1.2,  the neutron population would go from 1 neutron to 83 million neutrons! This state is called supercritical.

The spontaneous though at this point would be "how on earth is a reactor controlled if the power can increase 80 million times in 10 microseconds". The secret is that some fission products decay and release neutrons with a time delay, those delayed neutrons makes the neutron generation lifetime much longer and instead of 10^{-7} the mean generation time turns into order of seconds. I am not going to delve into delayed neutrons however but it is good to know that they exist.

The job of a nuclear engineer, at least those of us who simulate how the neutrons behave in a reactor, revolves around the all powerful K. The K value is the end all, be all of reactor physics. We want to make sure that when all control rods are in that K is less than 1, we want to make sure K = 1 during operation of the plant and we want to make sure it never goes much above 1 regardless of what happens to the reactor. Keeping K = 1 is surprisingly easy because there are a lot of natural feedbacks that does the trick without us ever lifting our hands, but that is also a topic for another time. Lets now look at some different ways to calculate K, because that is what we are interested in here.

We will start from the easiest and most approximative method, the four factor formula and then move on to the most computationally intense monte carlo methods and between those extremes look at diffusion methods. First however we need to understand some simple concepts.

Cross sections and thermal neutrons

What does the word thermal signifies in reactor physics and what is a cross section? When we speak of thermal neutrons or a thermal spectrum we are talking about neutrons that have an energy that is about the same as energy the surrounding atoms have due to motion caused by heat. Heat makes everything vibrate and move and when the neutron has an energy about equal to this slight motion it is called thermal. A neutron is born in fission with an energy millions of times higher than thermal energies, but it loses energy by colliding with atoms until it has an energy equal to the thermal energy of the atoms it is colliding with. This slowing down process is called moderation and a material that moderates is called a moderator.

Then we have cross sections. I steal this reasoning straight from Lamarsh book "Introduction to nuclear reactor theory"". If we shot a burst of neutrons into a slab of material nuclear reactions will take place. Lets assume we have a thin slice of a material with a thickness x and area A and that we shoot a ray of neutrons on it with intensity N. How many reactions will happen in the material? If we double the area we will obviously get double the reactions, so it must be proportional to A, if we double the thickness we also double the reactions, assuming the thickness is small enough not to consume to much of the beam. If we double the amount of atoms in the material, N, then we again obviously get double the amount of reactions. One can express this as the formula

 \large Reaction\ rate\ =\ \sigma INAX

NAX is obviously the number of atoms in the target (thickness, times area is volume and N is density of atoms, density times volume is total number of atoms). I is the beam intensity, lets assume it is 1. What is left is the quantity \sigma, this is some kind of proportionality constant and it can be interpreted as the likelihood that a neutron in the beam will react with an atom in the material. Simply put, number of atoms in the target, times this constant, times number of neutrons in the beam gives the number of reactions. \sigma is called the cross section of the atom and could be tought of as the size of the atom that the neutron experiences when it goes through the material. Its the same as if you randomly throw darts on a wall with holes in it, then the probability of hitting a hole is dependent on the total area of the holes. The curious property of cross sections however is that they are not constant, it depends on how fast the neutron is. It would be like a wall with holes that shrink if you throw the dart hard and enlarges if you throw the dart softly. Below is an example of how the seize of the cross section of one uranium isotope changes with the energy of the incoming neutrons.

As we see the cross section is much larger at low energies and that is why thermal neutrons are important, if we slow down neutrons to thermal energies (0.025 eV) the probability of nuclear reactions are much higher!

\sigma almost always appears together with the atom density of the material, N. To simplify notation one often multiplies them together into another constant \Sigma.

\large \Sigma = \sigma N

The four factor formula

Now lets look at the simplest method to estimate K. Despite its simplicity its a very powerful way to estimate the criticality in a system. It is called the four factor formula. In the four factor formula we simplify the system and only care about the materials present, we do this by assuming the system is infinite in all directions. The effect of this is that neutrons can not leak out of the system, they are all created and absorbed in the material and nothing else can happen to them. The four factor formula, as the name implies, contains 4 quantities that that I will go through. (The inf in K_{inf} signifies that the system is infinite and no leakage occurs).

 \large K_{inf} = \eta_{T}fp\epsilon

  • First we have \eta_{T} which is the number of neutrons emitted by fission per thermal neutron absorbed. It can be expressed as \large \eta_{T} = \nu \frac{\Sigma_f}{\Sigma_a} Where \nu is the number of neutrons produced per fission, \Sigma_f is the fission cross section and \Sigma_a is the total absorption cross section (absorption includes fission as well as other reactions that doesn't cause fission).
  • Second is f which is called the thermal utilization factor. f gives the probability that a neutron gets absorbed in the fuel instead of in non fuel materials in the material mixture. It is defined as \large f = \frac{\Sigma_{aF}}{\Sigma_{aF}+\Sigma_{aM}} Where \Sigma_{aF} is the absorption cross section of fuel and \Sigma_{aM} is absorption cross section of all other materials.
  • p is a quantity called the resonance escape probability. If we look again at the cross section pictured earlier we see that in the energy region 10 - 1000 eV there is a lot of spikes in the cross section, these are called resonances. They exist because if a neutron with a specific energy hits a nucleus and the energy matches precisely what is needed for the nucleus to enter a excited state then the probability of absorption of the neutron increases. Nuclear excited states are analogues to when an atom absorbs light (see wikipedia for an explanation). Because the resonances are so large it is quite likely that an atom will be caught there before becoming thermal. Because of the complexity of the resonances this value is not easy to calculate and one relies on experimental data.
  • \epsilon is the fast fission factor and it gives the ratio between the number of neutrons produced by all fission (including those caused by fast neutrons) and neutrons produced by fission in thermal energies.

That was quite an handful, thankfully most of the values of the 4 quantities exists in tables and one can estimate quite a lot from the equation without needing to care about the details. A lot of the behavior of a reactor can be understood conceptually from how these 4 factors behave during different circumstances. One drawback of the four factor formula as described above is that it assumes the material is mixed homogeneously. In other words there is no physical separation between fuel, moderator, structural materials and so on, this is obviously not the case for a real reactor where one has fuel pins that are clearly separated from the moderating water. But after a meltdown all the fuel, much of the steel structure, the cladding of the fuel etc has melted and formed a approximately homogeneous melt at the bottom of the reactor vessel (or at the bottom of the containment if unlucky). The four factor formula can then be used to roughly asses the K_{inf} of the molten core and see if it can go critical.

To simplify everything further we can make some assumptions. Lets assume there is no resonance capture by putting the resonance escape factor, p, to 1 (in an equal mixture of uranium and hydrogen it would be about 0.8, it would get increasingly smaller with more uranium and less hydrogen). The justification behind that conservatism is that we will smear out the neutron spectrum and resonance capture and normal capture will become the same thing. This is because we have so little moderation. The fast fission factor, \eta will also be put to 1 for the same reason.

Now we only need to figure out the values of \nu\Sigma_f, \Sigma_a,  \Sigma_{aF} and \Sigma_{aM}. In the fuel we have both Uranium-238 and Uranium-235, the factor

\large\eta_{T}=\nu\frac{\Sigma_f}{\Sigma_a}

then turns into (I have separated the macroscopic cross sections into microscopic onces in order to introduce enrichment in the next step)

\large\eta_{T}= \nu\frac{N^{235}\sigma^{235}_f}{N^{235}\sigma^{235}_a +N^{238}\Sigma^{238}_a}

Here we need to think of enrichment. N^{235} and N^{238} are related by the enrichment E as, where N^{U} is uranium atom density.

\large E = \frac{N^{235}}{N^{235}+N^{238}} = \frac{N^{235}}{N^{U}} => N^{235}=EN^{U}

\eta_{T} can then be rewritten as.

\large\eta_{T}=\nu\frac{E\sigma^{235}_f}{E\sigma^{235}_a +(1-E)\Sigma^{238}_a}

With little moderation the cross sections will be about the same as one finds in a fast spectrum reactor, \sigma^{235}_f=2.12\sigma^{235}_a=2.84\sigma^{238}_a=0.41. All cross sections are given in the unit barn which is equal to 1\ b=10^{-24}\ cm^{2}. Someone might now object and state that on top of the melt there might be a layer of water, this means the spectrum must be more thermal than what one can see in a fast reactor. To account for that and be conservative let us double the U-235 fission cross section to 4.24 barn.

Now what remains is the thermal utilization factor f. In the core melt mixture we have a lot of junk, there are a whole bunch of fission products, there is molten control rods, steel etc. All of those things absorb neutrons and lower the K. Once again, to be conservative, lets ignore all that junk just to see what K we get in a nice melt that is only uranium dioxide and nothing else. In other words all neutrons gets absorbed in the fuel and f=1.

To summaries we now have:

  • \large\eta_{T}=\nu\frac{E\sigma^{235}_f}{E\sigma^{235}_a +(1-E)\Sigma^{238}_a}
  • f =1
  • p = 1
  • \eta = 1

The four factor formula reduces to one factor.

K_{inf}=\large\eta_{T}=\nu\frac{E\sigma^{235}_f}{E\sigma^{235}_a +(1-E)\Sigma^{238}_a}

When a core is loaded fresh fuel has an enrichment between 3-5% and when a fuel bundle is unloaded it has enrichment of about 1% and another % or so of fissile Pu (we can assume the fissile Pu is the same as one % more U235 in this case). Average enrichment in the core is from 2'ish % to 5'ish %. What is the K value if we assume all fuel assemblies are fresh and nice with an enrichment of 5%? Turns out to be K = 0.76, way below criticality. Lets look at how it varies with enrichment.

We need to have an enrichment (or fissil content) of nearly 8.3% to go critical given the assumptions above. Factoring in the lowering effect of fission products, control rods etc we would end up with a much higher enrichment needed. If we then also take into account that the real life molten core is not infinite we also end up with a lot of neutron leakage that further lowers K. The conclusion we can draw from this simple example is that if the core has melted down into a nice homogeneous blob at the bottom of the vessel then there is no way it can go critical.

There is one possible oversight in the above reasoning however. If there is a thick layer of water on top of the melt, then by the upper edge of the melt a lot of neutrons will go out into the water and some of them will return back nicely moderated. In the worst case the upper layer of the melt might become critical due to this reflection of thermal neutrons. The four factor model is not refined enough for that and in order to calculate such a situation we need to turn to a model that can take into account geometry and regions with different materials in it. Next refinement will be a diffusion model and that will be described in part 2.