Criticality and delayed neutrons

Criticality and delayed neutrons

Dr Chris Robbins, Grallator

Fissile nuclides have the property that when struck by a neutron they split into two smaller nuclides and, crucially, liberate more neutrons; typically two or three. These liberated neutrons can either be lost from the fissile system by leaking away at the periphery or by being absorbed by nuclides without leading to a fission even, or they can go on to produce another fission event. The net number of neutrons produced per fission is called the neutron multiplication factor, \(k\), and is given by the ratio of neutrons produced to neutrons lost with causing further fission,

\[ k = \frac{\mathrm{produced}}{\mathrm{absorbed \space + \space leakage}} \]

At each generation of production the number of neutrons present will change by a factor of \(k\). There are three regions of interest for \(k\):

• \(k<1\): The number of neutrons falls with each generation.
• \(k=1\): The number of neutrons stays the same with each generation.
• \(k>1\): The number of neutrons increases with each generation.

As fission events release a large amount of energy, the situation where \(k=1\) can sustainably produce a significant power output. The situation where \(k>1\) is very undesirable as not only can the power output be significant, it is also increasing exponentially - the essence of an explosion.

The evolution of the number neutrons (which is a proxy for power), can be described by the expression \[ \frac{dn}{dt}=\frac{(k-1)}{\ell}n \] where:
\(\quad\quad n \space \) is the number of neutrons
\(\quad\quad \ell \space \) is the mean neutron lifetime

The solution to this equation is \[ n=n_0 \space \mathrm{exp}\left(\frac{(k-1)}{\ell}t \right) \] where \( n_0 \space \) is the number of neutrons at time \(t = 0\). If \(k<1 \) the solution is an exponential decay. If \(k=1 \) the solution is constant at \(n_0\), and if \(k>1 \) the solution is an exponential increase. These three solution types are call sub-critical, critical and supercritcal respectively.

The above is all well and good if all the neutrons are a direct result of fission, or so-called prompt neutrons. In reality, the physics is a little more complicated as not all neutrons are prompt neutrons; a very small fraction of the neutrons are released as a result of decay of the fission products produced when the fissile atom is split. For example, one possible fission product is \(\mathrm{^{87}Br}\), which has a half life of about 55 seconds and which has as one of its decay modes to \(\mathrm{^{87}Kr}\) one which emits a neutron. There are other fission products that also decay by emitting a neutron that have shorter half life values. However, the important point is that these neutrons are produced a long time after the fission event relative to the prompt neutrons. The effects of these delayed neutrons is discussed below.

Modelling delayed neutrons

All the neutrons produced will be either prompt or delayed. Define the fraction of delayed neutrons of the \(i^{th}\) precursor with a decay constant \(\lambda_i\) as \[ \beta_i=\frac{\mathrm{delayed \space neutrons_i}}{\mathrm{delayed \space neutrons_i \space + \space prompt \space neutrons_i}} \] The total delayed neutron fraction is \(\beta=\sum_i \beta_i\) and the prompt neutron fraction is therefore \(1-\beta\). Using these the neutron evolution equation becomes \[ \begin{align} \frac{dn}{dt}&=\frac{(k(1-\beta)-1)}{\ell}n + \sum_i {\lambda_i p_i} \\ \\ \frac{dp_i}{dt}&=\frac{k\beta_i}{\ell}n - \lambda_i p_i \end{align} \] where \(p_i\) is the number of atoms of precursor \(i\). The first term in the upper equation gives neutrons produced from prompt neutron fission, while the second adds in the delayed neutrons produced by precursor decay. The first term in the lower equation gives the delayed neutron precursor production, while the second accounts for precursor decay.

Often, the precursors are averaged in to groups with six groups being common. For this analysis, we'll take things further and average all the precursors in to a single item with decay constant \(\lambda\) (which we will use as \(\lambda = 0.077 \space s^{-1}\) ) so that the above becomes \[ \begin{align} \frac{dn}{dt}&=\frac{(k(1-\beta)-1)}{\ell}n + \lambda p \\ \\ \frac{dp}{dt}&=\frac{k\beta}{\ell}n - \lambda p \end{align} \] which can be written as the matrix differential equation \[ \begin{align} \frac{d}{dt} \begin{bmatrix}n\\p\end{bmatrix} &= \begin{bmatrix} \frac{(k(1-\beta)-1}{\ell} && \lambda \\ \frac{k\beta}{\ell} && -\lambda\end{bmatrix} \begin{bmatrix}n\\p\end{bmatrix} \\ \\ &=A\begin{bmatrix}n\\p\end{bmatrix} \end{align} \] where \(A\) represents the matrix. To solve this requires some initial conditions. We can arbitrarily take \(n_0 = 1\) and scale any final results if physical values are required. However, in this study we're only going to be looking at relative changes in power so relative changes in neutrons can be used as a proxy and absolute values are not required. The initial value of \(p\), \(p_0\), is also required. To calculate this we will assume the system is in a static equilibrium with a value of \(k=1\) so that \[ \frac{dn}{dt}=\frac{dp}{dt}=0 \Rightarrow \frac{k\beta}{\ell}n_0 - \lambda p_0 = 0 \Leftrightarrow p_0=\frac{k\beta}{\lambda \ell}n_0 = \frac{\beta}{\lambda \ell} \] and we will seek solutions for the evolution of the number of neutrons (remember this is being used as a proxy for power) with time when there is a sudden change in the value of \(k\). Such a change could arise from adding material, changing geometry, adding a reflector, adding a moderator etc.

Looking at the neutron part of the problem, the general solution will be of the form \(n=n_{1} \space \mathrm{exp}\left(\alpha_1 t \right) + n_{2} \space \mathrm{exp}\left(\alpha_2 t \right)\) where \(\alpha_1\) and \(\alpha_2\) are the eigenvalues of the matrix \(A\), and \(n_1\) and \(n_2\) are constants determined by the initial conditions (i.e. \(n_0\) and \(p_0\) above) and the eigenvectors of the matrix \(A\). Note, for consistency \(n_1 + n_2 = n_0\) etc.

The eigenvalues of \(A\) are found from the characteristic equation \[ \mathrm{det}\left( A - eI \right) = 0 \] where \(I\) is the identity matrix and \(e\) is a scaler - an eigenvalue. As \(A\) is a 2 x 2 matrix this determinant is easily evaluated as \[ \mathrm{det}\left( \begin{bmatrix} \frac{k(1-\beta)-1}{\ell}-e && \lambda \\ \frac{k\beta}{\ell} && -\lambda - e \end{bmatrix} \right) =\left( \frac{k(1-\beta)-1}{\ell}-e \right) \left( -\lambda - e \right) - \frac{k\beta\lambda}{\ell} \] which is quadratic in \(e\) indicating there will generally be two solutions, which will be the \(\alpha_1\) and \(\alpha_2\) constants in the general solution, and which give the timescale over which changes occur. Note, \(\alpha_1\) and \(\alpha_2\) will be real and may be positive or negative, indicating the solutions are exponential growths or decays. Once the eigenvalues are known, the eigenvectors can be constructed by substituting the eigenvalues in to \(\left( A - eI \right)E=0\) for eigenvector E. The important feature of an eigenvevtor is its direction and not its magnitude, and as we're looking at solutions for \(n\) it's convenient to write the eigenvectors in the form \[ E=\begin{bmatrix}1 \\ \eta \end{bmatrix} \]

where \(\eta\) is calculated from \(\left( A - eI \right)E=0\) when the first component of \(E\) is set to 1. The general solution for \(n\) and \(p\) is now given as \[ \begin{bmatrix}n\\p\end{bmatrix}=c_1 \begin{bmatrix}1 \\ \eta_1 \end{bmatrix} \mathrm{exp}\left(\alpha_1 t \right) + c_2 \begin{bmatrix}1 \\ \eta_2 \end{bmatrix} \mathrm{exp}\left(\alpha_2 t \right) \] for coefficients \(c_1\) and \(c_2\), which are determined from the intial conditions at \(t=0\) (for which \(\mathrm{exp}\left(\alpha_1 t \right)=1\) etc ), i.e. \(n=n_0=1\), \(p=p_0= \frac{\beta}{\lambda \ell}\),i.e. \[ \begin{align} n_0=1&=c_1+c_2 \\ p_0=\frac{\beta}{\lambda \ell}&=c_1\eta_1+c_2\eta_2 \end{align} \]

Example effects

The above theory shows that when starting from critical, \(k=1\) and there is a sudden change of multiplication factor there are, for this simplified case, two timescales that determine how the number of neutrons changes with time. Recall, we are using number of neutrons as a proxy for fission power, so these timescales are relevant to the rate of change of power after a change in \(k\), which may be caused by adding material, changing geometry, adding a reflector, adding a moderator etc. In all cases below, where there are delayed neutrons it will be assumed that the production decay constant is \(\lambda = 0.077 \space s^{-1}\).

No delayed neutrons and short neutron lifetime

In this case we will use \(\beta = 0\) and \(\ell = 10^{-8}\), and an instantaneous change of \(k\) to \(k=1.001\). The results are plotted below, and show that in this case the system rapidly increases in power, reaching ten times the initial power in 23 microseconds and one hundred times the initial power in 46 microseconds. Such an excursion would clearly be uncontrollable as the power increase occurs very quickly.

Delayed neutrons and long neutron lifetime

In this case we'll assume that there are a small but significant fraction of delayed neutrons and that prompt neutron lifetimes are longer by using \(\beta = 0.06\) and \(\ell = 10^{-4}\), and keep an instantaneous change of \(k\) to \(k=1.001\). These values are representative of a commercial thermal nuclear reactor. The results here show how delayed neutrons vastly increase the timescale over which change takes place. In this the ten times power level is reached in 140 seconds and the one hundred times power level is reached after 290 seconds. These values do not change appreciably even if the prompt neutron lifetime is reduced by several orders of magnitude as the speed of change is determined by the delayed neutrons. Such change times allow for control intervention to prevent the power build up.

A larger excursion

Using the above values of \(\beta = 0.06\) and \(\ell = 10^{-4}\), there will be a reduction in the time it takes the power to grow as the instantaneous change in \(k\) is increased. For example, when \(k=1.005\), the ten times and hundred times power levels are reached after 1.7 and 8.1 seconds respectively. Note, for a power reactor this would be an extremely large reactivity event.

The figure above also shows the effect of two timescales for change. Zooming in and looking at the ten times power region shows an initial rapid increase in power from the increase in \(k\), however, the increase is soon slowed down and becomes an increase with a speed determined by the delayed neutrons. It should be noted that the increase in \(k\) in the last two examples is lower than the delayed neutron fraction. Increasing \(k\) above this level means that the system will be prompt-critical and the delayed neutrons will not slow down the rate of power increase, which will change on a timescale determined by the prompt neutron lifetime.

Reducing the reactivity

The previous cases have all used \(k>1\) to see how power increases. Also of interest is what happens when control is applied to reduce the value of \(k\) in terms of how fast the power output from fission falls. The left figure below shows the response for \(\beta = 0\) and \(\ell = 10^{-8}\), and an instantaneous change of \(k\) to \(k=0.995\). As with the increase in power case, the response is very rapid, taking less than 5 microseconds to fall to 10% of the initial value. The right figure below shows the response for \(\beta = 0.006\) and \(\ell = 10^{-4}\), and an instantaneous change of \(k\) to \(k=0.995\). The sudden drop associated with reducing the number of prompt neutrons is apparent, however, the delayed neutrons extend out the time to reduce to 10% of initial power, in this case to 48 seconds.

Contacts

Dr Chris Robbins can be contacted at this website: Grallator.