# Precipitate Escalation of Reactivity in a Nuclear Reactor as its Criticality Proceeds from Delayed to Prompt [closed]

﻿The question actually stems in the first place from a question I put in the Physics forum about the Chernobyl nuclear accident - whether it was a true nuclear explosion. It is well known, and a concrete fact of nuclear engineeeing that the delayed neutron fraction is essential to reactor operation: as soon as the proportion excess of reactivity begins to exceed the delayed neutron fraction, the time-scale for exponential growth changes precipitately from the time delay of the delayed neutrons (10s of seconds) to the mean time a neutron is free in the reactor core (ms!). The longer timescale obtains within the narrow band of reactivity excess in which the system would be in exponential decay on the shorter timescale, and is only in exponential growth by reason of the small fraction of feedback that occurs on the longer timescale; and that within this band, howsoever narrow it might be, the total timescale is completely dominated by the longer one, rather than being some weighted mean of the longer & the shorter, as might be thought intuitively. There are various simple models that exhibit this property, showing that it is a very elementary property of that kind of exponentially growing system in which some fraction of the feedback is delayed. But it must be born in mind that these simple models are a long way from giving an accurate expression in detail for the behaviour of the reactor - the purpose of exhibiting them here is merely to show the properties under examination, that of the complete domination of the timescale by the longer one in the narrow band just described, and the very sudden flipping of the time scale with passing beyond this band, very robustly belonging to the exponential-growth nature of the system.

Three parameters introduced are $$\alpha$$ - the proportion of the total feedback that occurs on the longer timescale; $$\beta$$ - the deficit in the feedback occuring on the shorter timescale; and $$\zeta$$ - the ratio of the longer timescale to the shorter. The quantities $$\alpha$$ & $$\beta$$ are typically small fractions (maybe a fraction of a percent), whereas $$\zeta$$ is typically a large number, maybe hundreds or thousands.

A simple model is that of a sort of extreme Fibonnaci-like sequence

$$N_{k+\zeta+1}=(1-\beta)N_{k+\zeta}+\alpha N_k$$

If the solution $$N_k = N_0(1+\epsilon)^k$$ be inserted, $$\epsilon$$ becomes a root of the equation

$$(1+\epsilon)^\zeta = \frac{\alpha}{\epsilon+\beta}$$ ... ... ... *

I was wondering recently over what range approximating the high degree term as an exponential is valid, thereby enabling an expression for the root in terms of lambertw function. But doing it anyway for the present, we arrive at the solution

$$\epsilon = \frac{w(\zeta\alpha\exp(\zeta\beta))}{\zeta}-\beta$$

with $$w(x)$$ here being the solution in $$y$$ of

$$ye^y=x$$

& $$y$$ & $$x$$ both positive - that is, a branch of the lambertw function.

In the case of $$\alpha$$ & $$\beta$$ both positive, and $$\alpha>\beta$$ (the narrow band of operation with the long timescale), we can use the approximation

$$w(x)=\ln x - \ln\ln x(1-\frac{1}{\ln x})$$

as the argument to the lambertw function is obviously very large, and this gives

$$\epsilon ≈ \frac{\ln(\frac{\alpha}{\beta})(1 -\frac{1}{\zeta\beta})}{\zeta}$$

neglecting terms of higher order in $$\zeta$$.

The terms without $$\zeta$$ in the denominator vanish, leaving only those with at least the first power of $$\zeta$$ in the denominator.

At this stage we can also check the validity of using $$\exp(\zeta\epsilon)$$ as an approximation for $$(1+\epsilon)^\zeta$$. The proportion-error in this approximation is (unto quadratic terms) $$\zeta\epsilon^2/2$$, which in this case is of the order of atmost a few ÷ $$\zeta$$, so it is safe.

The approximation for the case of $$\beta\rightarrow 0$$ is obviously not suitable; but clearly

$$\epsilon=\frac{w(\zeta\alpha)}{\zeta}$$.

which, if we make the assumption that ζα is quite large even without the exponential in ζβ to multiply it,

$$\frac{\ln(\zeta\alpha)-\ln(\ln(\zeta\alpha))(1-\frac{1}{\ln(\zeta\alpha)})}{\zeta}$$

We can begin to see a hint of the precipitate rise in $$\epsilon$$ now, as whereas before, the numerator was the logarithm of atmost a few, now it is the lambertw of a smallish fraction of a large number; and although lambertw is actually sublogarithmic we can expect that it is now significantly bigger that it was before.

It gets really interesting when $$\beta$$ goes negative; beause then the argument of the lambertw goes rapidly very small by reason of that whereby it went large before - that it's argument has an exponential in $$\zeta$$ × something. In this régime, we have

$$\epsilon ≈ \alpha\exp(\zeta\beta)-\beta$$

which is now dominated by $$\beta$$ itself. Quite frankly, the approximation of exponential for finding this root has certainly by now become abrogate ... but it doesn't matter anymore, because the error is on a term that only makes a miniscule contribution now.

I find it fascinating that such a simple analysis with an approximation for finding the root yields a solution that reptoduces the behaviour so very well. But in this case you don't even need to get an analytical solution. You can actually demonstrate qualitatively the behaviour purely graphically, by plotting the functions in equation *. When β>0, the extremely steep (1+ε)^ζ curve obviously intersects the 1/(ε+β) part at a point of the order of (atmost a few)/ζ from the point ε=1; but when β<0, the 1/(ε+β) part now lies a distance -β to the right of the line ε=1, and the solution for ε only departs from -β by the amount by which the intersection is now to the right of this line, which is miniscule by reason of the extreme stepth of the (1+ε)^ζ curve, combined with the way the 1/(ε+β) curve cleaves to the vertical axis with increasing heighth up it.

It still seems unsatisfactory though that the negligible part is so much in error. You could improve it by noting that in the region of ε=-β the exponential is now overestimating by a factor of approximately $$\exp(-\zeta(\beta +\ln(1-\beta))$$, and because it is a simple reciprocal with which the intersection is being made, the departure-from-beta must be multiplied by this same factor: whence we get, by what is effectively the first step in a polishing of the solution by iteration

$$\epsilon = \frac{w(\zeta\alpha\exp(\zeta\beta))\exp(-\zeta(\beta +\ln(1-\beta))}{\zeta}-\beta$$

$$≈\alpha\exp(-\zeta\ln(1-\beta))-\beta$$

$$≈\frac{\alpha}{(1-\beta)^\zeta}-\beta$$

(remembering that β is negative in this part!) but it scarcely matters really, and is just really an expedient for removing the irkity arising from knowledge that one of the terms is grossly in error - negligible a term though-be-it. I do find it appealing, though how pleasant results keep dropping out.

A second method for setting up a model is that of the delay defferential equation.

﻿Defining α & β as before, the model would be a differential equation, with $$t_0$$ being the mean time between a neutron's being emitted & its causing a new fission.

$$\frac{d}{dt}N(t)=\frac{\ln(1-\beta)}{t_0}N(t)+\frac{\ln(1+\alpha)}{t_0}N(t-t_1)$$

This is a very flawed & very obviously flawed model of the behaviour of the neutrons in a reactor core, as it premises that the delayed neutrons are suddenly just dumped from a previous generation, $$t_1$$ earlier; but the behaviour of a delay equation is interesting in its own right, and does show forth the quality being queried here.

Insert parametrised solution $$N(t)=N_0\exp(\lambda t)$$.

If $$\alpha$$ & $$\beta$$ be assumed sufficiently small, this results in the equation for $$\lambda$$

$$\lambda t_0=-\beta + \alpha\exp(-\lambda t_1)$$.

Introducing parameter $$\zeta = \frac{t_1}{t_0}$$:

$$\lambda t_0=-\beta + \alpha\exp(-\zeta\lambda t_0)$$

$$\lambda t_0-\alpha\exp(-\zeta\lambda t_0)=-\beta$$.

Let $$x =\alpha\exp(-\zeta\lambda t_0)$$:

$$\frac{\ln(\frac{x}{\alpha})}{\zeta} + x = \beta$$

$$\ln(x)+\zeta x=\zeta\beta+\ln(\alpha)$$

$$\ln(\zeta x)+\zeta x=\zeta\beta+\ln(\zeta\alpha)$$

$$\zeta x=w(\zeta\alpha\exp(\zeta\beta))$$ .

This results in

$$\lambda t_0=\frac{\ln(\frac{\zeta\alpha}{w(\zeta\alpha\exp(\zeta\beta))})}{\zeta}$$.

This shows how wonderful the lambertw function is, in that it can solve $$x+e^x = y$$ as well as the function $$xe^x=y$$, which it is by definition the inverse of! You just have to exponentiate the argument & then take the logarithm of the answer!

This function might look different from - & yet similar to - the one obtained from the Fibonacci-like sequence method ... but remarkably it exhibits very similar behaviour with respect to changes in β: most particularly blowing up & becoming about equal to -β when β becomes ever so slightly negative. Examing the behaviour again with α & β both >0 & α > β, we can use the approxumation for the lambertw that was used before (but remembering that there's no ζβ to be subtracted this time!) & get

$$\epsilon ≈ \frac{\ln(\zeta\alpha)-\ln(\zeta\beta+\ln(\frac{\alpha}{\beta})(1 -\frac{1}{\zeta\beta}))}{\zeta}$$

$$\epsilon ≈ \frac{\ln(\frac{\alpha}{\beta})(1 -\frac{1}{\zeta\beta})}{\zeta}$$

which is exactly the same result as before, even unto the term in 1/ζβ!

The approximation in the case of β=0 is now

$$\epsilon=\frac{\ln(\zeta\alpha)-\ln(w(\zeta\alpha))}{\zeta}$$

which, if we make the assumption again about ζα being quite large even without the exponential in ζβ to multiply it,

$$\frac{\ln(\zeta\alpha)-\ln(\ln(\zeta\alpha))(1-\frac{1}{\ln(\zeta\alpha)})}{\zeta}$$

again the same as before!

Remarkably again, if we let β go substantially negative, and we take the _second term in the expansion of the lambertw function for small arguments, we recover

$$\epsilon ≈ \alpha\exp(\zeta\beta)-\beta$$

There is clearly an extremely close correspondence between these two methods for modelling the escalation of reactivity when the reactivity excess exceeds the delayed neutron fraction. It seems to stop at this point though, because there does not appear to be an analogue of that little trick for polishing the solution; that was shown before; but it scarecely matters really, as at this point the solution is utterly dominated by the -β 'undiluted' by any ζ in the denominator.

However, this assumption of the delayed neutrons being 'dumped' from the generation $$t_1$$ earlier is not a very satisfactory one. If instead we assume that they are from the decay of a suite of nuclides having mean life $$t_1, t_2, ... , t_n$$ and fractions of the total delayed neutron suite of nuclides $$\upsilon_1, \upsilon_2, ... , \upsilon_n$$ (with

$$\sum_{k=1}^n\upsilon_k=1$$),

having mean lives relative to $$t_0$$ $$\zeta_1, \zeta_2, ... , \zeta_n$$; and weight the size of the population from which the current generation of delayed neutrons coming 'online' proceeds by the decay curves of this suite of nuclides, then the function of $$\lambda t_0$$ multiplying $$\alpha$$ in the resulting equation becomes, instead of $$\exp(-\lambda t_d)$$,

$$\int_0^\infty\exp(-\lambda t_0\theta)\sum_{k=1}^n\frac{\upsilon_k}{\zeta_k}\exp(-\frac{\theta}{\zeta_k})d\theta$$

$$=\sum_{k=1}^n\frac{\upsilon_k}{\zeta_k\lambda t_0+1}$$

whence we would get the equation for $$\lambda t_0$$

$$\lambda t_0-\alpha\sum_{k=1}^n\frac{\upsilon_k}{\zeta_k\lambda t_0+1}=-\beta$$.

which is a polynomial of degree n+1!

If however we just assume there is a single nuclide (unindexed), we get

$$\lambda t_0-\frac{\alpha}{1+\zeta\lambda t_0}=-\beta$$.

$$\zeta(\lambda t_0)^2 +(1+\zeta\beta)\lambda t_0 = \alpha - \beta$$

which has solution

$$\lambda t_0 = \frac{1+\zeta\beta}{2\zeta}(\sqrt{1+\frac{4\zeta(\alpha - \beta)}{(1+\zeta\beta)^2}}-1)$$

This function has similar behaviour again ... but a little different in detail, as might be expected, as for $$\exp(-\zeta\lambda t_0)$$ we have substituted $$1/(\zeta\lambda t_0+1)$$ - a rational function for an exponential. But in the range α & β > 0, & α>β, the solution becomes

$$\lambda t_0 ≈ \frac{\alpha-\beta}{\zeta\beta+1} ≈ \frac{\frac{\alpha}{\beta}-1}{\zeta}$$

This resembles the solutions derived above containing $$\ln{\alpha/\beta}$$, inthat that function is roughly approximated by $$\alpha/\beta-1$$.

When β=0 we have

$$\lambda t_0≈\sqrt{\frac{\alpha-\beta}{\zeta}}$$

which does not really much resemble the previous solutions when β=0: but when β goes significantly negative, we now have

$$\lambda t_0 ≈ -\frac{\zeta\beta+1}{\zeta}(1+\frac{\zeta(\alpha-\beta)}{(\zeta\beta+1)^2}) ≈ -\beta$$

... again!

There is no trouble passing through the point β=-1/ζ ... but it must be kept in mind that when β passes through that point and the standard quadratic-equation-solution is factorized in that way, the sign attaching to the square-root symbol must flip from positive to negative in order to keep track of the same solution!!

And I think it's pretty obvious here that the higher-order corrections (ie in 1/ζ) aren't going much to resemble the ones in the previous two scenarios (though they showed a remarkable resemblence to each other) - which is not surprising by reason of the function according to which the delayed neutrons are phased in according to the size of earlier populations being of a fundamentally different nature.

What I have not explored yet is whether in that very hardly tractable case (actually the most realistic one) inwhich the contribution to the delayed neutrons was averaged over n species of nuclide resulting in a degree n+1 polynomial to solve ... does that actually morph into the exponential case to which the first two scenarios considered here were two complementary approaches? ... with a single somehow-averaged delayed neutron fraction & a single somehow-averaged delay ... as is indeed generally represented as being so in physics treatises.

When I say 'morphs into', what I have in mind is something similar to, say, what happens when you have many events, each of which alone might have an ever-so strange probability distribution for some one of it's indices; but when it is the probability distribution of some aggregate of all those indices that is under consideration (which might well be a simple sum of them), and the number of events increases without bound, the individual probability distribution gets washed out, and the distribution of the aggregate converges to a Gaussian ... always a Gaussian, unless there is some very special pathological reason for it not to do so. I wondered whether something similar might happen in this setup - corresponding to reality, as the delayed neutrons in a nuclear reactor originate in a wide range of fission products with a wide & random range of decay constants.

Update

I've just realised something about the similarity in behaviour of the first two solutions. I see now that they absolutely the same solution as setting them equal is tantamount to saying

$$\ln x=w(x) +\ln(w(x))$$

which is the very definition of w(x). This is only to be expected, as the two formulations are to each other essentially what

$$e^x$$ & $$\lim_{\zeta\rightarrow\infty}(1+\frac{x}{\zeta})^\zeta$$

are to each other.

Further Update

Actually this problem has become transparent now, through working it out. The question kept changing as I was writing this and the matter thereby became incrementally clearer. The final question was going to be "how is a single delay abstracted as a parameter from a medley of different nuclides, each being present in its own proportion, deaying exponentially, & having its own mean life?".

Basically you get an equation of the form

$$\lambda t_0 + \beta = \alpha f(\lambda t_0)$$

where f(x) is a function that begins at (0,1), descends with a certain steepness (ζ) towards the x-axis, then levels-off to tend asymptotically to the x-axis. These factors $$\zeta_k$$ are large numbers, being the ratio of the mean life of a neutron-delaying nuclide to the mean time a fission-inducing neutron is free in the reactor core, and typically are of the order of 100s or 1000s. Precisely because these are so large, the function f is extremely steep at its beginning; and precisely because it is so steep at the beginning, the precise form of it is nigh-on immaterial: to the extent that you could make it two straight line segments - one from (0,1) to (1/ζ,0), & another from (1/ζ,0) to (+∞,0); and it would make negligible practical difference to the outcome: and by the same token ζ might as well simply be

$$\sum_{k=0}^n\upsilon_k\zeta_k$$

I'm sure a proper mathematician would have seen that immediately!

## closed as unclear what you're asking by Servaes, user10354138, max_zorn, José Carlos Santos, ArsenBerkNov 14 '18 at 9:44

Please clarify your specific problem or add additional details to highlight exactly what you need. As it's currently written, it’s hard to tell exactly what you're asking. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

• 1. Is there a question buried in there somewhere? 2. However did you manage to express yourself before italics were invented? – Gerry Myerson Nov 13 '18 at 1:53
• Yes - the question is right at the end of the part before the update. It is whether or not when there is a mixture of many different proportions & lifetimes of neutron-delaying species, the overall scenario morphs into a delay-differential-equation or Fibonacci-like-sequence (as I have called the latter - it isn't standard usage!) scenario; which, as you will see if you look at the update, I now realise are essentally the same - which I ought really to have realised all along, but was distracted into missing by the seeming difference in the expeessions obtained for the rate. – AmbretteOrrisey Nov 13 '18 at 2:04
• $\it Careful$, there's an $\it impending\ worldwide\ shortage\$ of $\it italics$. – Gerry Myerson Nov 13 '18 at 6:11
• We can just use nuclear power to generate more! – AmbretteOrrisey Nov 13 '18 at 12:44