© V. Urbanevych et al., published by EDP Sciences, 2020

1 Introduction

The idea of autowave (or traveling-wave) regime in the medium with diffusion and multiplication of neutrons was suggested in the paper [1] more than 30 yr ago. In [2] authors give the equations which describe a physical system with a possible nuclear burning. These equations take into account among others the multi-stageness of burn-up reactions, the neutrons energy distribution and the possible dependence of the diffusion coefficient on the energy, time, and coordinate. These equations are too complicated for the exact analytical solutions being looked for, so one has to seek either for the simplified systems of equations or the numerical simulations. In a number of papers [1,3–8] authors regard significantly simplified expressions for this system. Using it, they try to analyze physical effects and look for the traveling-wave solutions to this problem. In present work, we consider the system of equations, similar to one from [1] and [8], but we suggest a different way to analyze it. Substantial difference between models of [1,4,7,8] and [3,5,6,9] consists lies in describing the combustion processes. In [3] author regards only one differential equation for the neutrons diffusion, what implies that the nonlinearity and nonlocality of the system are regarded as an effect of the fuel’s combustion, while the equations which describe these processes are not considered at all. In [5,6,9] all the kinetic equations for the fuel’s combustion do not have terms which are not proportional to the neutron flux and could provide the system with a contribution of the β-decay. Neglecting of these terms leads to the appearance of the integrals, which are not instilled in the original system. Considering these terms makes impossible to apply the methods, similar to that in [5,9]. There is an important difference in [1,8,10] and [5,11] and it is that in the differential equations we have some terms which are not proportional to the neutron flux and are responsible for the β-decay. Considering these terms makes impossible to apply the methods, similar to that in [5,11].

In [8] the numerical simulation was performed only with time-constant external neutron flux. Since the system is not left to itself in such a mode, one can question the interpretation of these results as autowave. In addition, the numerical calculation does not give an answer to two important, from our point of view, questions. Firstly, what are the conditions, which would lead the system of equations to the autowave solution? And secondly, where the speculations about conditions for the wave velocity come from? In [1,4,7] authors used an analogy with an eigenvalue problem in the famous Schrödinger equation in order to answer these two questions. In our opinion, the inconsistency in [1,4,7] arises from the fact that in the Schrödinger equation, potential does not depend on the eigenfunction, which is being looked for, while the analog of this potential in the system from [1] does depend by means of burn-up equations Besides that, the eigenfunction of the Schrödinger equation can take a negative value, in spite of the neutrons density. We would like to demonstrate that the analysis of the problem like in [1,4,7] can be performed using the analogy with the dynamic equation from the classical mechanics instead of the quantum equation which was used there. In such a new analogy the coordinate is substituted by the integral over the neutrons density, and time – by the autowave variable. In such case, the “force” (here and below the quotation marks mean, that it is not this physical quantity, but just its analog within our model) in the equation is consisted of the potential and dissipative parts. It comes to light that the traveling-wave mode is possible only when the “potential energy” function has the minimum point, which in classical mechanics led to the regime with oscillations. If the oscillations appear than in some part the velocity should become negative. As in our case, “the velocity” is proportional to the neutrons density, the solutions with the negative velocity leads to the unphysical negative neutrons density. The presence of dissipation can prevent the system from the negative velocity values if all initial stored “energy” will be dissipated before the minimum point is reached. Similarly, the ball in the elastic pendulum, which is brought out of the equilibrium in a very viscous medium will be able just to return to the equilibrium point, but not to oscillate. It will be shown below, that the “dissipative” contribution is directly connected to the β-decay processes, which are not taken into account in [5,9], and it is appeared to be proportional to the autowave velocity. Out of it, we come to the condition for this velocity - it has to be such, that “dissipation” would not allow oscillations. Hence the restrictions for the eigenvalues spectrum in the regarded model appear not out of the boundary conditions, but out of the non-negative solution requirement. We perceive the main difference of our work from another in this area, as e.g. in [1,4,7] the spectrum comes out from the boundary conditions.

Different kinds of nonlinear equations are being analyzed in [12]. In the paper, authors find the traveling-wave velocity from the condition that there is some phase trajectory in such a system, that connects points, consistent with the boundary conditions. A similar way of finding the wave velocity is used in [9], but using only the part of the trajectory, which corresponds to the positive (that is physical) values of neutron flux. But this path is being built with the help of integral, which can lead to an unphysical branch with negative values because of the square of neutrons’ flux in it. These two branches are connected in the origin of thecoordinates and the transition between the branches can happen at that point. Finally, the nonlinear eigenvalue problem in [11] is not concerning the velocity, but the effective multiplication coefficient. Again, the authors find the eigenvalues requiring that there is a nontrivial solution in such a system, which satisfies the boundary conditions.

We find that the traveling-wave velocities spectrum, which is determined by the condition of the non-negativity of the solution, is continuous and limited from below. Looking at the results of [12], such a situation is typical for the nonlinear autowaves. Next, we will discuss the above considerations in more detail. In Section 6 we will analyze the possible connection between the obtained results and more realistic models.

2 System of kinetic equations

Let us start with the system of kinetic (balance) equations for neutrons and nuclides, which describe the process of wave neutron-nuclear burning of the uranium-plutonium fissile medium. We use the work by L.P. Feoktistov [1] as a basis. Note that the system considered further on has a simplified form. As in [1], the equations describing the fission fragments are omitted, and the fuel consists only of ²³⁸U and there are no fissile nuclides like ²³⁵U, as in [13].

Let us first consider the equation for neutrons. With the mentioned simplifications in the diffusion approximation it looks as follows: $$\begin{array}{lll}\frac{\partial n\left(t,x\right)}{\partial t}\hfill & =\hfill & D\frac{{\partial }^{2}n\left(t,x\right)}{\partial x^2}+(\left(N_8^0\left|v\right|{\sigma }_{f,Pu}\right)\left(\nu -1\right)\hfill \\ \hfill & \hfill & -\left(N_8^0\left|v\right|{\sigma }_{c,Pu}\right))n\left(t,x\right)n_{Pu}\left(t,x\right)\hfill \\ \hfill & \hfill & -n\left(t,x\right)(n_8\left(t,x\right)\left(N_8^0\left|v\right|{\sigma }_{a,8}\right)\hfill \\ \hfill & \hfill & +n_9\left(t,x\right)\left(N_8^0\left|v\right|{\sigma }_{a,9}\right)),\hfill \end{array}$$(1)

where n, n₈, n₉ and n_Pu are the dimensionless concentrations of neutrons, ²³⁸U, ²³⁹U and ²³⁹ Pu respectively (all the concentrations are dimensionless and normalized by the initial concentration of ²³⁸U–N${}_0^8$); σ_c,i is the neutron capture cross-sectionfor the ith nuclide, σ_f,i is the fission cross-section, σ_a,i is the neutron absorption cross-section (σ_a,i = σ_f,i + σ_c,i); D is the neutron diffusion coefficient; ν is the average number of neutrons produced per one ²³⁹Pu fission reaction.

Quantities $N_8^0\left|v\right|{\sigma }_{f,\text{Pu}}$, $N_8^0\left|v\right|{\sigma }_{c,\text{Pu}}$, $N_8^0\left|v\right|{\sigma }_{a,8}$, and $N_8^0\left|v\right|{\sigma }_{a,9}$ have the dimension of inverse time. So one can treat them as the inverse mean free times for the neutrons with respect to the corresponding nuclear reaction with relevant nucleus, given the nuclei concentrations $N_0^8$. Therefore, they may be denoted as follows: $$\begin{array}{llll}{N}_8^0\left|v\right|{\sigma }_{f,\text{Pu}}\hfill & =\frac{1}{{\tau }_{f,\text{Pu}}},\hfill & N_8^0\left|v\right|{\sigma }_{c,\text{Pu}}\hfill & =\frac{1}{{\tau }_{c,\text{Pu}}},\hfill \\ N_8^0\left|v\right|{\sigma }_{a,8}\hfill & =\frac{1}{{\tau }_{a,8}},\hfill & N_8^0\left|v\right|{\sigma }_{a,9}\hfill & =\frac{1}{{\tau }_{a,9}}.\hfill \end{array}$$(2)

However, the introduced mean free times are approximately equal. At least, the difference between them is much less than their difference from the overall time scale of the problem – the characteristic time of β-decay τ_β. So in our work we apply a simplification, that $${\tau }_{f,\text{Pu}}\approx {\tau }_{c,\text{Pu}}\approx {\tau }_{a,8}\approx {\tau }_{a,9}=\tau .$$(3)

In addition, one can switch to the dimensionless time coordinate, which is normalized by the τ and spatial coordinate x, defined as follows: $$x=\sqrt{D\tau }y.$$(4)

All these simplifications bring the equation (1) to the next form: $$\begin{array}{lll}\frac{\partial n\left(t,y\right)}{\partial t}\hfill & =\hfill & \frac{{\partial }^{2}n\left(t,y\right)}{\partial y^2}\hfill \\ \hfill & \hfill & +\left(\left(\nu -2\right)n_{\text{Pu}}\left(t,y\right)-n_8\left(t,y\right)-n_9\left(t,y\right)\right)n\left(t,y\right)\text{}.\hfill \end{array}$$(5)

Next, we would like to formulate the kinetic equations which we mentioned in the introduction. These equations will be responsible for the fuel combustion with a next chain of reactions: $${}_{92}^{238}\text{U}+n{\to }_{92}^{239}\text{U}{\to }_{94}^{239}\text{Pu}.$$(6)

We suppose that ²³⁸U can only burn out and cannot accumulate in any way. Then its kinetic equation is $$\frac{\partial N_8\left(t,x\right)}{\partial t}=-{\sigma }_{a,8}\left|v\right|N_8\left(t,x\right)\tilde{n}\left(t,x\right).$$(7)

Using the same dimensionless densities and coordinates one obtains: $$\frac{\partial n_8\left(t,y\right)}{\partial t}=-n_8\left(t,y\right)n\left(t,y\right).$$(8)

²³⁹U is produced during the neutron capture by ²³⁸U. Its amount decreases due to absorption of neutrons (neutron capture or nuclear fission) and β-decay with characteristic time τ_β. So the equation for ²³⁹U is $$\begin{array}{lll}\frac{\partial N_9\left(t,x\right)}{\partial t}\hfill & =\hfill & \left({\sigma }_{a,8}{N}_8\left(t,x\right)-{\sigma }_{a,9}{N}_9\left(t,x\right)\right)\left|v\right|\tilde{n}\left(t,x\right)\hfill \\ \hfill & \hfill & -\frac{1}{{\tau }_{\beta }}{N}_9\left(t,x\right).\hfill \end{array}$$(9)

After scaling, the equation (9) takes the form: $$\frac{\partial n_9\left(t,y\right)}{\partial t}=\left(n_8\left(t,y\right)-n_9\left(t,y\right)\right)n\left(t,y\right)-\frac{\tau }{{\tau }_{\beta }}{n}_9\left(t,y\right).$$(10)

Since ²³⁹Pu is produced by the β-decay of ²³⁹ U, and burns out absorbing neutrons, the kinetic equation for ²³⁹Pu is $$\begin{array}{lll}\frac{\partial N_{\text{Pu}}\left(t,x\right)}{\partial t}\hfill & =\hfill & \frac{1}{{\tau }_{\beta }}{N}_9\left(t,x\right)\hfill \\ \hfill & \hfill & -\left({\sigma }_{f,\text{Pu}}+{\sigma }_{c,\text{Pu}}\right)N_{\text{Pu}}\left(t,x\right)\left|v\right|\tilde{n}\left(t,x\right)\text{}. \hfill \end{array}$$(11)

After simplification and combining together with equations (5) (8) and (10) one can obtain the next system of equations: $$\{\begin{array}{l}\frac{\partial n}{\partial t}\left(t,y\right)=\frac{{\partial }^{2}n}{\partial y^2}\left(t,y\right)\hfill \\ +\left(\left(\nu -2\right)n_{\text{Pu}}\left(t,y\right)-n_8\left(t,y\right)-n_9\left(t,y\right)\right)n\left(t,y\right),\hfill \\ \frac{\partial n_8\left(t,y\right)}{\partial t}=-n_8\left(t,y\right)n\left(t,y\right),\hfill \\ \frac{\partial n_9\left(t,y\right)}{\partial t}=\left(n_8\left(t,y\right)-n_9\left(t,y\right)\right)n\left(t,y\right)-\frac{\tau }{{\tau }_{\beta }}{n}_9\left(t,y\right),\hfill \\ \frac{\partial n_{\text{Pu}}\left(t,y\right)}{\partial t}=\frac{\tau }{{\tau }_{\beta }}{n}_9\left(t,y\right)-\frac{1}{{\tilde{n}}_{\text{Pu}}}{n}_{\text{Pu}}\left(t,y\right)n\left(t,y\right)\hfill \end{array}$$(12)

where $${\tilde{n}}_{\text{Pu}}=\frac{{\sigma }_{c,8}}{{\sigma }_{f,\text{Pu}}+{\sigma }_{c,\text{Pu}}}.$$(13)

System (12) is almost the same as in [8], but in the diffusion equation we use explicitly the contribution of the combustion of fuel (namely ²³⁸U and ²³⁹U).

As one can see, a very small parameter $\frac{\tau }{{\tau }_{\beta }}$ emerged in the problem – the ratio of the neutron mean free time τ (about 10⁻⁷ s) to the characteristic time of β-decay τ_β ≈ 3.3 days. We shall denote this parameter as: $$\epsilon \equiv \frac{\tau }{{\tau }_{\beta }}\approx {10}^{-14}[\text{s}].$$(14)

Let us convert the resulting system of kinetic equations into the autowave form. For this purpose, we make the variable substitution according to the relation z = y + ut. With such substitution, the burning wave goes in the negative direction of the coordinate axis of the autowave variable z. The derivatives change as follows: $$\begin{array}{l}\frac{\partial n}{dt}=u\frac{\partial n}{dz},\\ \frac{{\partial }^{2}n}{d{y}^2}=\frac{{\partial }^{2}n}{d{z}^2}.\end{array}$$

System (12) in autowave form reads $$\{\begin{array}{l}u\frac{dn}{dz}\left(z\right)=\frac{d^{2}n}{d{z}^2}\left(z\right)\hfill \\ +\left(\left(\nu -2\right)n_{\text{Pu}}\left(z\right)-n_8\left(z\right)-n_9\left(z\right)\right)n\left(z\right),\hfill \\ u\frac{d{n}_8\left(z\right)}{dz}=-n_8\left(z\right)n\left(z\right),\hfill \\ u\frac{d{n}_9\left(z\right)}{dz}=\left(n_8\left(z\right)-n_9\left(z\right)\right)n\left(z\right)-\epsilon n_9\left(z\right),\hfill \\ u\frac{d{n}_{\text{Pu}}\left(z\right)}{dz}=\epsilon n_9\left(z\right)-\frac{1}{{\tilde{n}}_{\text{Pu}}}{n}_{\text{Pu}}\left(z\right)n\left(z\right).\hfill \end{array}$$(17)

The second equation in system (17) may be integrated. Than we obtain the following expression for the ²³⁸U concentration: $$n_8\left(z\right)=n_{8,-\infty }\exp \left(-\frac{1}{u}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)d{z}_1}\right),$$(18)

where n_8,−∞ is the ²³⁸U concentration when $\left(z\to -\infty \right)$. Since we use the values scaled to theinitial ²³⁸U concentration, n_8,−∞ = 1.

Applying the variable substitution and some basic mathematical transformations, one can integrate the rest of the equations from (17). The system of equations then takes the form: $$\{\begin{array}{l}\frac{1}{u}\frac{dn}{dz}\left(z\right)=n\left(z\right)-n_9\left(z\right)+\left(\nu -2\right){\tilde{n}}_{\text{Pu}}{n}_{\text{Pu}}\left(z\right)\hfill \\ -\frac{\epsilon }{u}\left(1+\left(\nu -2\right){\tilde{n}}_{\text{Pu}}\right){\displaystyle {\int }_{-\infty }^z{n}_9\left(z_1\right)d{z}_1}\hfill \\ +2\left(1-\exp \left(-\frac{1}{u}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)d{z}_1}\right)\right),\hfill \\ n_8\left(z\right)=\exp \left(-\frac{1}{u}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)d{z}_1}\right),\hfill \\ n_9\left(z\right)=\frac{1}{u}\exp \left(-\frac{1}{u}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)}d{z}_1\right)\hfill \\ \times \left(\exp \left(-\frac{\epsilon }{u}z\right){\displaystyle {\int }_{-\infty }^{z}n\left(z_2\right)\exp \left(\frac{\epsilon }{u}{z}_2\right)d{z}_2}\right),\hfill \\ n_{\text{Pu}}\left(z\right)=\frac{\epsilon }{u}\exp \left(-\frac{1}{u{\tilde{n}}_{\text{Pu}}}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)d{z}_1}\right)\hfill \\ \times {\displaystyle {\int }_{-\infty }^z{n}_9\left(z_2\right)\exp \left(\frac{1}{u{\tilde{n}}_{\text{Pu}}}{\displaystyle {\int }_{-\infty }^{z_2}n\left(z_1\right)d{z}_1}\right)d{z}_2}.\hfill \end{array}$$(19)

Let us consider the last equation of the system (19) in more detail. We then direct argument to plus infinity in this expression, i.e. $\left(z\to +\infty \right)$. Then, we have the following: $$\begin{array}{lll}{n}_{\text{Pu}}\left(z\to +\infty \right)\hfill & =\hfill & \frac{\epsilon }{u}\exp \left(-\frac{1}{u{\tilde{n}}_{\text{Pu}}}{\displaystyle {\int }_{-\infty }^{+\infty }n\left(z_1\right)d{z}_1}\right)\hfill \\ \hfill & \hfill & \times {\displaystyle {\int }_{-\infty }^{+\infty }{n}_9\left(z_2\right)\exp \left(\frac{1}{u{\tilde{n}}_{\text{Pu}}}{\displaystyle {\int }_{-\infty }^{z_2}n\left(z_1\right)d{z}_1}\right)d{z}_2}.\hfill \end{array}$$(20)

This expression may be conventionally split into three multipliers: the constant $\frac{\epsilon }{u}$, the exponent, and the integral. The integrand is the exponent multiplied by ²³⁹U concentration, and the interval of integration is from $\left(-\infty \right)$ to $\left(+\infty \right)$. The exponent is strictly more than zero throughout the entire integration interval (the properties of this function). The ²³⁹U concentration is a non-negative value, and it has to be greater than zero at some points because otherwise there would be no reaction andno wave of nuclear burning at all. So the integral is of the product of two positive functions, and consequently, is also a positive quantity which is certainly not equal to zero. Since the constant (the quotient of two positive values) and the exponent (property of the function) are also positive values, then the entire product is certainly more than zero. Thus it may be concluded that the plutonium cannot burn out completely. Instead, it tends to some constant level: $n_{\text{Pu}}\left(z\to +\infty \right)\equiv n_{\text{Pu},+\infty }\ne 0$. This is also confirmed by the results of numerical simulation of the burning kinetics, e.g. [14,15]. The fact, that plutonium can not completely burn out indicates, that in the active zone of traveling-wave nuclear burning the reproduction coefficient is greater than unity, which the boundary condition $n_{\text{Pu}}\left(z\to +\infty \right)\equiv n_{\text{Pu},+\infty }\ne 0$ is corresponds to.

Let us assume the following: $${\displaystyle {\int }_{-\infty }^{z}d{z}_2}n\left(z_2\right)\exp \left(\frac{\epsilon }{u}{z}_2\right)\approx n\left(z\right){\displaystyle {\int }_{-\infty }^{z}d{z}_2}\exp \left(\frac{\epsilon }{u}{z}_2\right).$$(21)

Applying the approximation (21), we proceeded with the following considerations. The integral in (21) is ${\displaystyle {\int }_{-\infty }^{z}n}\left(z_2\right)\exp \left(\frac{\epsilon }{v}{z}_2\right)d{z}_2$. If we substitute the function $n\left(z_2\right)$ under the integral by its value at the upper limit $n\left(z_2\right)\approx n\left(z\right)$, then for the points z₂ near the upper limit of integration, such substitution it is possible because of the approximate equality $n\left(z_2\right)\approx n\left(z\right)$ for the small deviation of z₂ from z. On the other hand, for the points far from the upper limit of integration such substitution becomes possible due to the fact that $\exp \left(\frac{\epsilon }{v}{z}_2\right)\ll \exp \left(\frac{\epsilon }{v}z\right)$ for z₂ ≪z. So the concentration $n\left(z_2\right)$ in this case will enter the integral with a coefficient that is exponentially small in comparison with that of n(z). Therefore, one can expect that the value of the integral will be much less sensitive to the accuracy of the approximation in the region z₂ ≪z as compared to the region, where z₂ is close to z, where the approximation $n\left(z_2\right)\approx n\left(z\right)$ should be a good one.

With this approximation, the system of equations (19) will take the form: $$\{\begin{array}{l}\frac{1}{u}\frac{dn}{dz}\left(z\right)=n\left(z\right)-n_9\left(z\right)+\left(\nu -2\right){\tilde{n}}_{\text{Pu}}{n}_{\text{Pu}}\left(z\right)\hfill \\ -\frac{\epsilon }{u}\left(1+\left(\nu -2\right){\tilde{n}}_{\text{Pu}}\right){\displaystyle {\int }_{-\infty }^z{n}_9\left(z_1\right)d{z}_1}\hfill \\ +2\left(1-\exp \left(-\frac{1}{u}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)d{z}_1}\right)\right),\hfill \\ n_8\left(z\right)=\exp \left(-\frac{1}{u}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)d{z}_1}\right),\hfill \\ n_9\left(z\right)=\frac{1}{\epsilon }\exp \left(-\frac{1}{u}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)}d{z}_1\right)n\left(z\right),\hfill \\ n_{\text{Pu}}\left(z\right)=\frac{{\tilde{n}}_{\text{Pu}}}{1-{\tilde{n}}_{\text{Pu}}}[\exp \left(-\frac{1}{u}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)}d{z}_1\right)\hfill \\ -\exp \left(-\frac{1}{u{\tilde{n}}_{\text{Pu}}}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)d{z}_1}\right)].\hfill \end{array}$$(22)

Since the burning wave is propagating in a negative direction, the boundary conditions are set at minus infinity and have the following form: $$\begin{array}{llll}\hfill & n\left(z\to -\infty \right)=0,\hfill & \hfill & n_8\left(z\to -\infty \right)=n_{8,-\infty }=1,\hfill \\ \hfill & n_9\left(z\to -\infty \right)=0,\hfill & \hfill & n_{\text{Pu}}\left(z\to -\infty \right)=0,\hfill \end{array}$$(23)

where n_8,−∞ is the initial ²³⁸U concentration.

3 The analogy to Newton’s second law

Let us consider the first equation of the system (22). As we have already derived the expressions for ²³⁹U and ²³⁹ Pu from other equations, we can substitute them here. After substitution we can collect terms and integrate some parts of this equation. Finally we obtain the following expression: $$\begin{array}{lll}\hfill & \hfill & \frac{1}{u}\frac{dn}{dz}\left(z\right)=\left[1-\frac{1}{\epsilon }\exp \left(-\frac{1}{u}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)}d{z}_1\right)\right]n\left(z\right)\hfill \\ \hfill & \hfill & +\left(1-\left(\nu -2\right){\tilde{n}}_{\text{Pu}}\right)\hfill \\ \hfill & \hfill & +\frac{(\nu -1){\tilde{n}}_{\text{Pu}}-1}{1-{\tilde{n}}_{\text{Pu}}}\exp \left(-\frac{1}{u}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)}d{z}_1\right)\hfill \\ \hfill & \hfill & -\frac{\left(\nu -2\right){\left({\tilde{n}}_{\text{Pu}}\right)}^2}{1-{\tilde{n}}_{\text{Pu}}}\exp \left(-\frac{1}{u{\tilde{n}}_{\text{Pu}}}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)d{z}_1}\right)\text{}\text{}. \hfill \end{array}$$(24)

So we got integro-differential equation. Let us introduce a new variable N(z) into it, as follows: $$N\left(z\right)\equiv \frac{1}{u}{\displaystyle {\int }_{-\infty }^{z}n\left(z_1\right)}d{z}_1.$$(25)

The physical meaning of this value is neutron fluence divided by the wave velocity or called neutron density accumulation over time. Let us treat N(z) as an analog of coordinate. Note that the speed, i.e. a derivative of this coordinate, cannot be negative. In our case the derivative is a neutron concentration, thus a negative sign at the speed means the same sign at n(z), which takes us out of the physical region.

Taking into account (25), (24) becomes $$\begin{array}{lll}\hfill & \hfill & \underset{\text{resultant}}{\underbrace{\frac{d^{2}N\left(z\right)}{d{z}^2}}}=\underset{\text{viscous force}}{\underbrace{u\left(1-\frac{1}{\epsilon }\exp \left(-N\left(z\right)\right)\right)\frac{dN\left(z\right)}{dz}}}\hfill \\ \hfill & \hfill & +\underset{}{\underbrace{\left(1-\left(\nu -2\right){\tilde{n}}_{\text{Pu}}\right)+\frac{(\nu -1){\tilde{n}}_{\text{Pu}}-1}{1-{\tilde{n}}_{\text{Pu}}}\exp \left(-N\left(z\right)\right)}}\hfill \\ \hfill & \hfill & -\underset{\text{potential force}}{\underbrace{\frac{\left(\nu -2\right){\left({\tilde{n}}_{\text{Pu}}\right)}^2}{1-{\tilde{n}}_{\text{Pu}}}\exp \left(-\frac{N\left(z\right)}{{\tilde{n}}_{\text{Pu}}}\right)}}.\hfill \end{array}$$(26)

Let us note that equation (26) looks like this because we did not neglect the derivative over time in the system of kinetic equations, when switching to the autowave variable z = y + ut in Section 2, as it was done e.g. in [1,3,5,11,16,17].

One could draw a parallel between equation (26) and Newton’s second law. If we consider the N(z) an analog to coordinate, its second derivative is the acceleration. So one has a resultant force scaled to mass on the left in equation (26). On the right, there is a sum of a viscous force (a term including velocity) and a force caused by some potential of interaction (which is a gradient of potential energy with a negative sign). It is easy to build an expression, the negative derivative of which would coincide with this term and which can be regarded as an analog to the potential energy. One can easily obtain the equation (26) in the form: $$\begin{array}{l}\frac{d}{dz}\left(\frac{1}{2}{\left(\frac{dN\left(z\right)}{dz}\right)}^2+U\left(N\right)\right)\\ =-u\left(\frac{1}{\epsilon }\exp \left(-N\left(z\right)\right)-1\right){\left(\frac{dN\left(z\right)}{dz}\right)}^2,\end{array}$$(27)

where $$U\left(N\right)\equiv k_{1}N+k_2\exp \left(-N\right)+k_3\exp \left(-\frac{N}{{\tilde{n}}_{\text{Pu}}}\right),$$(28) $$\begin{array}{lll}{k}_1\hfill & \equiv -\left(1-\left(\nu -2\right){\tilde{n}}_{\text{Pu}}\right),\hfill & k_2\equiv \frac{(\nu -1){\tilde{n}}_{\text{Pu}}-1}{1-{\tilde{n}}_{\text{Pu}}},\hfill \\ k_3\hfill & \equiv -\frac{\left(\nu -2\right){\left({\tilde{n}}_{\text{Pu}}\right)}^3}{1-{\tilde{n}}_{\text{Pu}}}.\hfill & \hfill \end{array}$$(29)

A derivative of the sum of kinetic and potential energy in the expression (27) equals to the viscous force. Thus this equation may be considered as the energy conservation law, and the total energy may change only through the work of the viscous force.

We are interested in finding cases of equilibrium, which corresponds to the autowave fission mode. This happens when the total force consisting of the “potential” force and “dissipative” viscous force, is zero. For this, the potential energy must have a point of minimum, because the derivative (the potential force) is zero at this point. The speed (and hence the kinetic energy) must also be zero at this minimum point, meaning the dissipative viscous force is zero. Analyzing different values of coefficients k₁, k₂ and k₃ depending on the ñ_Pu parameter, we found that in order to obtain the potential energy (28) with a minimum point, it is enough that the ñ_Pu > 1. If this condition is satisfied, than the potential energy has a form, as shown in Figure 1. There is an apparent point of minimum in this graph – a stable stationary point. The derivative of potential energy (the potential force) is obviously zero at this point. So this is a possible stationary state, if the kinetic energy is zero, and corresponds to the autowave fission mode.

Remind that the derivative of N cannot be negative, because it would lead us to negative neutron concentration, which is non-physical. If we treat N as an analog of some coordinate, then its derivative is the speed. So this speed cannot be negative – as non-physical. If the object starts from the point N = 0 and we want it to stop at the point of equilibrium (minimum), it should not overshoot this point. Otherwise the derivative of the potential energy will be non-zero, the force will act on the body, and the body will keep moving. It is known that the potential energy cannot be greater than total energy, so at some moment, when the kinetic energy completely transforms into potential, the object will start moving in the opposite direction. This means that the speed of the body will become negative. And we cannot allow this as non-physical. Therefore it is necessary for the entire kinetic energy to dissipate due to viscosity by the time when the potential energy reaches a point of minimum. This way point the object stops at a stationary point. Graphically this path is shown by the arrow in Figure 1.

As ñ_Pu and ν (via Eqs. (13) and (43), respectively) depend on the neutrons energy, in the Section 5 we will find all possible energy values, at which the potential function has a minimum point.

Fig. 1 A single possible form of potential energy that allows the existence of an auto-wave mode is shown in black. The arrow showsthe path of total energy along which this mode can be set.

4 Finding a criterion for the wave velocity

We determined that for the autowave mode to be possible, the coefficients k₁, k₂ and k₃ (taking into account the value ñ_Pu) must be chosen in such a way, that the potential energy (or its analog) has a point of minimum. The potential force (the derivative of U(N) with opposite sign) in this case, will look like it is shown in Figure 2 (solid line).

Since in the autowave form the wave propagates over an infinite interval, and there is no explicit dependence on variable z in equation (27), it is always possible to perform a variable substitution that shifts z relative to the function N. Let us thus switch to a new variable so that z = 0 at the point $\frac{N_0}{2}$. We replace the function of force by a triangular one (dashed line in Fig. 2).

Next, we try to replace the function of force with two straight lines – one on the negative part of the coordinate axis z, and the other on the positive part. If we solve the equation (27) for N with the functions of force in the form of two lines, we obtain two unknown constants for each one of solutions (because the equations are of second-order). We adjust these constants to satisfy the boundary conditions. We also hope that the constants adjustment will end up with zeroing out one of them for each solution since the functions at which they appear will diverge to infinity. Then, if we join the functions and their first derivatives at the point $\frac{N_0}{2}$, we would obtain one spare equation, which can be used to find the wave speed.

So, let us first consider the interval $z\in \left(-\infty ,0\right)$. Here we make the following substitution: $$\begin{array}{l}F\left(N\right)=qN,\\ U\left(N\right)=U_0-q\frac{N^2}{2},\end{array}$$

where q is the slope, and is equal to: $$q=-{\frac{d^{2}U\left(N\right)}{d{N}^2}|}_{N=0}=-\left(k_2+\frac{k_3}{{\tilde{n}}_{\text{Pu}}{}^2}\right)=1.$$(32)

By substituting (31) into (27) we obtain a homogeneous differential equation of second order, whose solution looks like: $$\begin{array}{l}N\left(z\right)=C_1\exp \left({\chi }_{1}z\right)+C_2\exp \left({\chi }_{2}z\right)\\ {\chi }_{1,2}=\frac{1}{2}\left(-\frac{u{n}_{8,-\infty }}{\epsilon }\pm \sqrt{{\left(\frac{u{n}_{8,-\infty }}{\epsilon }\right)}^2+4q}\right).\end{array}$$

χ₁ is apparently positive, and χ₂ is negative. So in order for the N(z) to converge to (−∞), C₂ must be zero. So $$N\left(z\right)=C_1\exp \left({\chi }_{1}z\right).$$(35)

Now let us consider the interval $z\in \left(0,+\infty \right)$. In this case we make the following substitution: $$\begin{array}{l}F\left(N\right)=-\left(N-N_0\right),\\ U\left(N\right)=U\left(N_0\right)-\frac{1}{2}{\left(N-N_0\right)}^2.\end{array}$$

If we substitute the potential energy again into (27), we obtain a similar solution. $$\begin{array}{lll}N\left(z\right)\hfill & =\hfill & N_0+D_1\left(u\right)\exp \left({\lambda }_1\left(u\right)z\right)\hfill \\ \hfill & \hfill & +D_2\left(u\right)\exp \left({\lambda }_2\left(u\right)z\right),\hfill \end{array}$$(38) $$\begin{array}{lll}{\lambda }_{1,2}\hfill & =\hfill & \frac{1}{2}[-u\frac{n_{8,-\infty }}{\epsilon }\exp \left(-N_0\right)\hfill \\ \hfill & \hfill & \pm \sqrt{{\left(\frac{u{n}_{8,-\infty }}{\epsilon }\exp \left(-N_0\right)\right)}^2-4}].\hfill \end{array}$$(39)

As we can see, both λ₁ and λ₂ are negative, so at (+∞) both constants are at the “good” exponents (the converging ones), and we cannot zero out any of them. So this way we do not obtain an additional equation, which could be used to find the wave speed.

There is an expression that may be negative under the root in equation (39). Thus, the solution will be the sum of sine and cosine with some coefficients. Since they have a finite period, and the wave propagates on an infinite interval, the derivative will become negative at some point. As we noted above, this would lead us to the negative neutron concentration, which is non-physical. So we require that $${\left(\frac{u}{\epsilon }\exp \left(-N_0\right)\right)}^2-4\ge 0.$$(40)

Considering this requirement relative to u, we obtain the following restriction: $$u\ge 2\epsilon \exp \left(N_0\right).$$(41)

It means, that the wave velocity may take a continuous spectrum of values greater than a certain minimum value. In the framework of the analogy to the energy conservation law, one might say that the viscosity cannot be less than a certain value. It makes sense, because with higher viscosity the kinetic energy will still dissipate completely to the minimum point, but if the viscosity is not high enough, the body will possess some nonzero speed at the stationary point, which eventually leads to the reverse motion. And as we noted above, the reverse motion is not allowed as non-physical.

Let us check if the speed can take different values. To do this, we choose the parameters so that the potential energy has the point of minimum, and solve the equation (26) numerically. These parameters are: ñ_Pu = 0.74, ν = 2.8, N₀ = 1.8, where N₀ is the point of potential energy minimum. With such values, the condition (41) has the form: u ≥ 15ε. Figure 3 shows the dependences of N on z for different velocities. Comparing it to the potential energy (Fig. 1a), one cannot notice the following: the coordinate is initially zero, and then it goes to the point of the potential energy minimum. The change in speed, i.e. dissipative force, affects the form of this transition and shifts it in time.

Fig. 2 A solid line shows a schematic representation of the potential force F(N) (a derivative of U(N) taken withthe opposite sign), when the potential energy has a minimum point (the type of potential from Fig. 1c). A dashed line shows an approximate force which we use in our work for while looking for the wave velocity criterion.

Fig. 3 Solutions of equation (26) as a function N(z) when a minimum exists in potential energy. The curves are given for different velocities that satisfy the requirement (41).

5 The spectrum of neutron energies, suitable for an auto-wave mode

In order to make possible the existence of the traveling-wave mode, the potential energy (or its analog) must have a point of minimum, and also the wave speed must be such that the kinetic energy is completely dissipated when the potential energy reaches this minimum point. So now as we realize that, it is interesting to study how the fulfillment of the first condition (existence of the minimum point) depends on the neutrons energy. According to (29) and (28), U(N) depends on ν and ñ_Pu. These parameters, in turn, depend on the neutron energy. We find the parameter ñ_Pu using equation (13) and examine it depending on the neutron energies.

According to [18–20], the mean number of instantaneous neutrons ν produced bysingle fission has the following dependence on energy: $$\begin{array}{lll}\nu \left(Z,A,E_n\right)\text{}\text{}\hfill & \hfill & =2.33+0.06\left(2-{\left(-1\right)}^{A+1-Z}-{\left(-1\right)}^Z\right)\hfill \\ \hfill & \hfill & +0.15\left(Z-92\right)+0.02\left(A-235\right)\hfill \\ \hfill & \hfill & +\left(0.13+0.006\left(A-235\right)\right)\left(E_n-E_{\text{threshold}}\right)\text{}.\hfill \\ \hfill & \hfill & \hfill \end{array}$$(42)

Since we are interested in the average number of neutrons ν for Pu²³⁹, we choose the following parameters for the equation (42): $$A=239,Z=94,E_{\text{threshold}}=-0.89\left(\text{MeV}\right).$$(43)

So now we know the dependence of the coefficients in the potential energy on the energy of the neutrons. We have to find the entire spectrum of neutron energies at which the potential energy has a point of minimum. It may be done as follows. Break some interval of N (from zero up to a certain maximum value) into sufficiently small segments, and compare the values of U(N) at three adjacent points. If the value at the middle point is less than those at both ends, then the minimum of the function exists. If we do not find such three points, then there is no minimum at this energy. Having analyzed the entire spectrum of neutron energies for which it was possible to find a cross-section values in [21], that is from 0 to 15 MeV this way, we find all the energy intervals with minima.

We show all the neutron energies, for which we found the minimum of potential energy, as a set of points in the graph, where the neutron energy is along the X-axis and the parameter ñ_Pu is along the Y -axis. Only the values, for which there is a stationary point in the potential energy, are shown.

Figure 4 shows the values of the parameter ñ_Pu for neutron energies in the range from 0 to 100 eV. We used the dependence of the cross-sections of nuclides on the neutron energy from the database [21]. As seen from Figure 4 (according to the considerations in Sect. 3), there are seven regions at about 6–7 eV, 19–21 eV, 34–40 eV, 66–67 eV, 81–82 eV and 90–100 eV, having a minimum of potential energy, which makes the autowave of nuclear fission possible. The spectrum of ñ_Pu values, where the minimum of potential energy was found, looking among the all the neutron energies which were available in [21] (from 0 to 15 [MeV]) is shown in Figure 5.

It is worth mentioning, that results described in Figure 5 for fast neutrons region have a preliminary nature. This is due to the fact the calculating with one-group approximation admits a significant error in the fast neutron spectrum. Therefore it is necessary to use the multigroup approximation for quantitative calculations.

For all the previous considerations we assumed that in equation (27) the right-hand side is always negative. We did it as the dominant negative term in the brackets is proportional to $\frac{1}{\epsilon }~{10}^{14}$. Nevertheless if N(z) is large enough, than the exponent can make the term very small and the total viscous force will become positive, what is difficult to interpret within our Newton’s law analogy. As we are interested in a potential of Figure 1 - kind, than the maximum possible value of N(z) is a N₀. If in that point the viscous force in equation (27) is still zero, than all are speculations above are valid. So we checked the value of the coefficient $\left(\frac{1}{\epsilon }\exp \left(-N_0\right)-1\right)$ for the different energies of neutrons, but requiring that the potential equation (28) has a minimum point and N₀ is this point ofminimum. The values of this factor are shown at Figure 6 and one can see that all of them are much greater than zero, what means that viscous force in equation (27) takes always only the negative values.

Fig. 4 The dependence of the parameter ñPu on the energy of neutrons, for the energies, when the potential energy has the point of minimum, and the autowave may exist. The image shows the energy range from 0 to 105 eV.

Fig. 5 The dependence of the parameter ñPu on the energy of neutrons, for the energies, when the potential energy has the point of minimum, and the existence of an autowave is possible. This picture embraces all possible points, as the entire region of neutron energies, where it was possible to obtain an experimental data for the cross-sections [21]. Namely, the energy range is (0, 15) MeV.

Fig. 6 Coefficient at the viscous force equation (27) taken at the point of Potential energy minimum as a function of the energy of neutrons, where there is a minimum point.

6 Numerical solution of the non-stationary system of equations

To confirm the conclusions of the theoretical model, by the numerical simulation we selected the epithermal region of neutrons energy. Also, we will try to confirm the relation between automodel regime existence and the form of potential energy analog. As follows fromFigure 5, the form of potential energy analog in this region (about 5.9 eV) is appropriate for the automodel regime’s existence. Note, that research of autowave nuclear burning regime existence is the foreground for authors of this article. e.g. [13].

The system of equations studied above is simplified and does not fully reflect the physical processes that take place in the reactor core.Let us consider the system of equations given in [13], which is more accurate (though not ideal), and consists of 18 equations. The kinetic equation for the neutron concentration $$\frac{\partial n\left(y,t\right)}{\partial t}=D\frac{{\partial }^{2}n\left(y,t\right)}{\partial y^2}+q\left(y,t\right),$$(44)

where q(y, t) is the internal source of neutrons, which has the form: $$\begin{array}{lll}q\left(y,t\right)\hfill & =\hfill & \left[{\nu }^{\text{Pu}}\left(1-p^{\text{Pu}}\right)-1\right]n\left(y,t\right)V_n{\sigma }_f^{\left(\text{Pu}\right)}{N}_{\text{Pu}}(y,t)\hfill \\ \hfill & \hfill & +\left[{\nu }^5\left(1-p^5\right)-1\right]n(y,t)V_n{\sigma }_f^{\left(5\right)}{N}_5(y,t)\hfill \\ \hfill & \hfill & +\ln 2{\displaystyle \sum _{i=1}^6\left[\frac{{\tilde{N}}_i^{\left(\text{Pu}\right)}(y,t)}{T_i^{\left(\text{Pu}\right)}}+\frac{{\tilde{N}}_i^{\left(5\right)}(y,t)}{T_i^{\left(5\right)}}\right]}\hfill \\ \hfill & \hfill & -n(y,t)V_n{\displaystyle \sum _{5,8,9,Pu}{\sigma }_c^{\left(i\right)}{N}_i(y,t)}\hfill \\ \hfill & \hfill & -n(y,t)V_n{\displaystyle \sum _{i=1}^6\left[{\sigma }_c^{i\left(\text{Pu}\right)}{\tilde{N}}_i^{\left(\text{Pu}\right)}(y,t)+{\sigma }_c^{i\left(5\right)}{\tilde{N}}_i^{\left(5\right)}(y,t)\right]}\hfill \\ \hfill & \hfill & -n(y,t)V_n{\sigma }_c^{\text{eff}\left(\text{Pu}\right)}{\overline{N}}^{\left(\text{Pu}\right)}(y,t)\hfill \\ \hfill & \hfill & -n(y,t)V_n{\sigma }_c^{\text{eff}\left(5\right)}{\overline{N}}^{\left(5\right)}(y,t).\hfill \end{array}$$(45)

In equations (44) and (45) n(y, t) is the neutron concentration; V_n is the neutron velocity; ν^Pu and ν⁵ are the mean number of instantaneous neutrons per ²³⁹Pu and ²³⁵U fission respectively.N₅, N₈, N₉ and N_Pu are the ²³⁵U, ²³⁸U, ²³⁹U and ²³⁵ Pu concentrations respectively;${\tilde{N}}_i^{\text{Pu}}$ and ${\tilde{N}}_i^5$ are the concentrations of fragments with neutron excess, produced due to ²³⁹ Pu and ²³⁵U fission respectively;${\overline{N}}^{(\text{Pu})}$ and ${\overline{N}}^{(5)}$ are the concentrations of all the other fission fragments of ²³⁹Pu and ²³⁵U respectively;σ_c and σ_f are the microscopic cross-sections of the neutron capture and nuclear fission; the parameters p_i $\left(p={\displaystyle \sum _{i=1}^6{p}_i}\right)$ and T_1∕2 characterizethe groups of delayed neutrons. They are well known and presented in [20]. ${\sigma }_c^{\text{eff}}$ is some effective microscopic cross-section of neutron capture by fragments.

One can find the whole set of kinetic equations, together with the constants used there in [13]. We use here dimensionless concentrations and coordinates, normalized in the same way as we did for the simplified system of equations. Now we assume, that initial medium consists not only from the ²³⁸U, but also from ²³⁸U and initial conditions are: $$\begin{array}{llll}{N_8\left(y,t\right)|}_{t=0}\hfill & ={\eta }_8,\hfill & {N_5\left(y,t\right)|}_{t=0}\hfill & ={\eta }_5,\hfill \end{array}$$(46)

where η₈ and η₅ are constants adjusting the fuel enrichment (η₈ + η₅ = 1). The remaining elements are initially absent.

In contrast to [8,13] we want to find the traveling-wave mode, so when there are no external neutrons come to the system and it burns in itself. In order to obtain it, we want firstly to switch the external neutron flux on and then, after some time, switch it off. We want to check if the system can set up the equilibrium with a soliton-like burning wave. The initial and boundary conditions for neutrons we set as $$\begin{array}{llll}n\left(y,0\right)\hfill & =0,\hfill & n\left(0,t\right)\hfill & =n_0\left(t\right),\hfill \end{array}$$(47)

where $n_0\left(t\right)$ is some function specifying the number of external neutrons and which satisfies our statements above.

For the numerical simulation, we use Wolfram Mathematica. We perform the calculations for the energy of neutrons E_n = 5.923 (eV). During the manipulation with a full system of equations we obtain a similar parameter as in equation (14). Unfortunately, at this stagewe cannot solve the system with a real value of this parameter, so let us take some approximation. We take ε = 10⁻³. This value is not realistic (ε is about 10⁻¹⁴), but it is difficult to carry out the numerical calculations with more precise values for our model.

We want to set the function n₀(t) in such a way that the external neutrons flux will be switched on for some time and then it drops to zero. In particular, we set the dependence of the boundary neutrons density on the time as it is in the equation (48). With such a function the neutrons come only for a time period around t =10 (in our dimensionless units) and then it goes to zero. As we want to see the setting of the autowave, we will show the results which are much beyond the time of switching off a source and out plots starts from t =10⁵. $$n_0(t)=3000\times t\times \exp \left(-0.5\times t\right).$$(48)

Figures 7–10 present the results of numerical simulation (with the neutron energy E_n = 5.923 (eV)) with the external source of neutrons being switched off already. In Figures 7–10 are presented the concentrations of neutrons, ²³⁸U, ²³⁹U and ²³⁹ Pu, respectively.The setting of the autowave mode of nuclear burning is confirmed by the fact that starting from a certain moment in time (when the external source of neutrons is already switched off), each subsequent curve differs from the previous only by a shift along the y-axis.

We also studied the regions of neutron energies for which there was no minimum in potential energy. In such cases, the autowave mode did not establish, and the wave faded away after the external source of neutrons was switched off.

The results of numerical modeling confirm the obtained theoretical conclusions. According to these results, one of the possible areas of autowave burning in the uranium-plutonium medium is the epithermal region of neutron energies (as in [13,22]). Meanwhile in the region of fast neutrons, the wave burning of the uranium-plutonium medium requires the constant supply of neutrons from an external source and the wave fades out when the source is switched off.

The discovered burning mode with external support can be used for implementing a traveling-wave reactor, e.g. to stop the burning at any time by switching off the external source of neutrons.

Fig. 7 The concentrations of neutrons as a function of the spatial coordinate y in the different moments of time (from t = 200 000 to t = 1 000 000) obtained from the numerical solution. The different lines correspond to the particular moments of time.

Fig. 8 The same as in Figure 7 but for the concentration of the 238U.

Fig. 9 The same as in Figure 7 but for the concentration of the 239U.

Fig. 10 The same as in Figure 7 but for the concentration of the 239Pu.

7 Conclusions

A kinetic system of equations, describing the wave mode of neutron-nuclear burning in the uranium-plutonium medium is formulated. Its autowave form is also obtained. In contrast to many other papers like [1,3,5,11,16,17], in which it was also considered in autowave form, we do not neglect the derivative of neutron concentration with respect to time in the neutron diffusion equation. This way we study the non-stationary traveling-wave burning mode.

For the first time, a kinetic equation for neutrons is obtained in the form of the energy conservation law for a mechanical system withdissipation, and the mechanical analogy for the fission traveling wave is developed. This allowed us to formulate the conditions for the existence of the autowave burning mode, and determine the possible values of the wave speed. The values of the burning wave rates have a continuous distribution and are bounded at the bottom, that is, they have a continuous spectrum. We also succeeded in determining the regions of neutron energy, for which the autowave burning is possible. We conclude that in other energy ranges, in which the autowave mode is impossible, the wave of nuclear burning may still be established using the support of an external neutron source. It should be noted, that these results are preliminary and further refinement is planned with the aid of more accurate calculations.

To confirm the theoretical conclusions, we performed a numerical 1D simulation of the neutron-nuclear burning in the uranium-plutonium medium in a single-group diffusion approximation (E_n = 5.923 (eV)). The results of numerical modeling confirm the obtained theoretical conclusions. According to these results, one of the possible areas of autowave burning in the uranium-plutonium medium is the epithermal region of neutron energies (as in [13,22]).

We regarded a significantly simplified model in the paper. Nevertheless, obtained results can be used as a first approximation in more realistic models. In particular, we did not count the energy distribution of neutrons and cross-sections (as it was done in [2]), but the limiting stage of the investigated processes was a β-decay, which does not depend on the energy of the neutrons. If one will use a multigroup approximation (as in [2]), then different groups of neutrons will have their typical times for the processes which happen to them. So during the problem solution, these times will be included as the ratios between them and their ratio to the typical β-decay time. Despite its small value, the ratio to the β-decay cannot be neglected as setting it to zero leads to the loosing of the neutrons multiplication channel and consequently, the traveling-wave solution is lost. Therefore, one cannot get the autowave iteratively, starting from an approximation, when all of them are equal to zero. But, as we showed in the work, one can obtain a traveling-wave solution putting them all equal to each other. Similar speculations can be applied to the ratios of typical neutrons times to each other. Putting them all equal to one, we get an initial solution for the iterative procedure. Having the wave at first approximation, it is possible to move on and count the multigroup effects by sequential iterations for discussed parameters. The fact that the β-decay plays a crucial role in the wave formation and it does not depend on the energy of the neutrons, allows us to expect that such an iterative approach will change a form and speed of the wave, but will not ruin the autowave solution.

One more approximation is that we work one the one-dimensional system within plane geometry. The possible improvement can be performed in a way as in [9,11], that is solving a system for a wave spreading in a cylindrical tube. In a diffusional approximation, the boundary conditions represent the continuity of the neutrons flux that is a continuity of the radial component of the gradient. It is possible to regard an ideally reflective or ideally absorbing tube walls, and neutrons will not be able to get inside the walls. In such a case the radial flux will be equal to zero, so this is a requirement for the boundary conditions. Such a requirement is also satisfied by the particular solution when the radial component is zero along all the tube surface. In that case, the problem is reduced to the one, regarded in the present work. So it is possible to take our results as a first approximation and move towards more realistic problems.

Author contribution statement

V. Urbanevych: performed analytical derivation of the model and numerical calculations in order to prove it. Prepared the figures and wrote the manuscript draft. I.V. Sharph: conceived of the model, supervised the work. V.A. Tarasov: co-supervised the work, worked on the text of the article. V.D. Rusov: discussions about the model, useful pieces of advice to improve the model. All authors put their efforts in order to obtain the final version of the article.

References

All Figures

	Fig. 1 A single possible form of potential energy that allows the existence of an auto-wave mode is shown in black. The arrow showsthe path of total energy along which this mode can be set.
In the text

	Fig. 2 A solid line shows a schematic representation of the potential force F(N) (a derivative of U(N) taken withthe opposite sign), when the potential energy has a minimum point (the type of potential from Fig. 1c). A dashed line shows an approximate force which we use in our work for while looking for the wave velocity criterion.
In the text

	Fig. 3 Solutions of equation (26) as a function N(z) when a minimum exists in potential energy. The curves are given for different velocities that satisfy the requirement (41).
In the text

	Fig. 4 The dependence of the parameter ñPu on the energy of neutrons, for the energies, when the potential energy has the point of minimum, and the autowave may exist. The image shows the energy range from 0 to 105 eV.
In the text

	Fig. 5 The dependence of the parameter ñPu on the energy of neutrons, for the energies, when the potential energy has the point of minimum, and the existence of an autowave is possible. This picture embraces all possible points, as the entire region of neutron energies, where it was possible to obtain an experimental data for the cross-sections [21]. Namely, the energy range is (0, 15) MeV.
In the text

	Fig. 6 Coefficient at the viscous force equation (27) taken at the point of Potential energy minimum as a function of the energy of neutrons, where there is a minimum point.
In the text

	Fig. 7 The concentrations of neutrons as a function of the spatial coordinate y in the different moments of time (from t = 200 000 to t = 1 000 000) obtained from the numerical solution. The different lines correspond to the particular moments of time.
In the text

	Fig. 8 The same as in Figure 7 but for the concentration of the 238U.
In the text

	Fig. 9 The same as in Figure 7 but for the concentration of the 239U.
In the text

	Fig. 10 The same as in Figure 7 but for the concentration of the 239Pu.
In the text