Corresponding author: Ivan V. Tormyshev ( itormyshev@ippe.ru ) Academic editor: Georgy Tikhomirov
© 2021 Ivan V. Tormyshev, Andrey V. Gulevich, Vladimir A. Yeliseev, Victor Yu. Stogov.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Tormyshev IV, Gulevich AV, Yeliseev VA, Stogov VYu (2021) Assessment of the neutron flux stability in a high power BNtype reactor in terms of modal spatial kinetics. Nuclear Energy and Technology 7(3): 201206. https://doi.org/10.3897/nucet.7.73488

The article discusses the neutron flux stability in the core of a highpower sodiumcooled fast reactor (of the BNtype) without feedbacks. The importance of this problem for highpower BNtype reactors is associated with the specific features of the layout of their cores, including a large diameter and height/diameter ratio about 5. The technique used to substantiate the stability of neutron fields is based on the analysis of the spectrum of the matrix of the system of spatial kinetics equations describing the core of a highpower BNtype reactor without feedbacks. A computational model of the spatial kinetics of a highpower BNtype reactor has been developed in the modal approximation based on the representation of an unsteady flux as a sum of orthogonal functions multiplied by timedependent amplitudes. The eigenfunctions of the conditionally critical problem are used in the diffusion approximation, which in the discrete case form a complete system. The spectrum of the matrix of the system of ordinary differential equations describing the spatial kinetics of the reactor has been calculated. It is shown that the neutron flux in the core of a highpower BNtype reactor without feedbacks is stable. Test calculations have illustrated the damping of perturbations of the power distribution for a reactor in a critical state.
BNtype reactor, neutron flux stability, spatial kinetics, conditionally critical problem, eigenfunction expansion
Numerous safety studies of sodiumcooled fast reactors have shown that, in order to exclude the possibility of destruction of the core in beyond design basis accidents of the ULOF type, the core height should be no more than 85–90 cm, and above it there should be a sodium plenum (
The concept of ‘stability’ as applied to nuclear reactors is formulated similarly to the concept of stability of any dynamic systems, i.e., the reactor state is called ‘stable’ if the deviations of the reactor parameters (including power) from the stationary values that have arisen after the introduction of a disturbance into the reactor do not increase indefinitely with time (limited response to limited impact). Otherwise, if the disturbance introduced into the reactor increases with time indefinitely, the reactor state is called ‘unstable’ (
Mathematical criteria for the stability of dynamical systems are formulated in the theory of differential equations.
The solution φ(t) of the system of differential equations
ẋ = f(t,x)
with the initial conditions x(0) = φ(0) is said to be Lyapunov stable if for any ε > 0 there is δ(ε) > 0 such that if x(0) – φ(0) < δ(ε), then x(t) – φ (t) < δ(ε) for any t ≥ 0. Thus, Lyapunov stability means that small perturbations of the initial condition do not lead to an unbounded increase in time of the perturbation of the solution x(t) (
If the solution φ (t) of the system satisfies not only the Lyapunov stability condition, but also the condition
$\underset{t\to \infty}{\mathrm{lim}}\Vert \mathbf{x}\left(t\right)\mathit{\phi}\left(t\right)\Vert =0,$
then the solution of the system is said to have the property of asymptotic stability (
One of the methods used to assess the stability of power distribution in nuclear reactors is to analyze the stability of a system of differential equations describing the kinetics of a reactor in the approximation of expansion in some system of functions (as a rule, in terms of eigenfunctions of a conditionally critical diffusion problem) (
In this work, the stability of the power distribution of a highpower BNtype sodiumcooled fast reactor without feedbacks is considered by analyzing the eigenvalues of a system of differential equations describing the reactor kinetics when the neutron flux is expanded in terms of the eigenfunctions of the conditionally critical problem. The nuclear reactor dynamics model without feedbacks is applicable to the analysis of the processes occurring in the reactor core when the reactor operates in the starting power range.
One of the ways to solve the problems of spatial kinetics is to expand the nonstationary neutron flux into a series in terms of the system of orthogonal functions:
$\Phi (\mathbf{r},\mathit{\Omega},E,t)=\sum _{i}{T}_{i}\left(t\right){\Phi}_{i}(\mathbf{r},\mathit{\Omega},E)$. (1)
As shown in (
Let us consider the nonstationary diffusion equation (coordinates r, E, t, on which the values entering into the equation depend, are omitted for brevity)
$\frac{1}{v}\frac{\partial \Phi}{\partial t}+(A+\Delta A)\Phi ={\chi}_{p}(1\beta )(F+\Delta F)\Phi +\sum _{j}{\lambda}_{j}{\chi}_{d,j}{C}_{j}$
$\frac{\partial {C}_{j}}{\partial t}={\beta}_{j}(F+\Delta F)\Phi {\lambda}_{j}{C}_{j}$ (2)
where A is the operator of diffusion, absorption and scattering of a neutron of the stationary conditionally critical problem describing the reactor critical state, defined as
$A\Phi (\mathbf{r},E)=\nabla D(\mathbf{r},E)\nabla \Phi (\mathbf{r},E)+{\Sigma}_{tot}(\mathbf{r},E)\Phi (\mathbf{r},E){\int}_{{E}^{\u2019}}\sum _{s}\left(\mathbf{r},{E}^{\u2019}\to E\right)\Phi \left(\mathbf{r},{E}^{\u2019}\right)d{E}^{\u2019}$
F is the neutron fission operator of the stationary conditionally critical problem describing the reactor critical state, defined as
$F\Phi (\mathbf{r},E)={\int}_{E}v{\Sigma}_{f}(\mathbf{r},E)\Phi (\mathbf{r},E)dE$
ΔA are ΔF are the operator disturbances that brought the reactor out of the critical state; Q is the external source; C_{j} is the concentration of delayed neutron group j precursor nuclei; v is the neutron speed; χ_{p} is the prompt neutron spectrum; β is the total effective fraction of delayed neutrons; β_{j} is the group fraction of delayed neutrons; λ_{j} is the is the decay constant of group j delayed neutrons; χ_{d},_{j} is the spectrum of group j delayed neutrons.
Let us substitute expansion (1) into system (2), choosing as functions Φ_{i} the eigenfunctions of the eigenvalue problem:
AΦ_{i} = χFΦ_{i} /K_{i}, (3)
$\sum _{i}\frac{\partial {T}_{i}}{\partial t}\frac{1}{v}{\Phi}_{i}+\sum _{i}{T}_{i}(A+\Delta A){\Phi}_{i}={\chi}_{p}(1\beta )\sum _{i}{T}_{i}(F+\Delta F){\Phi}_{i}\sum _{j}{\lambda}_{j}{\chi}_{d,j}{C}_{j}+Q$,
$\frac{\partial {C}_{j}}{\partial t}={\beta}_{j}(F+\Delta F)\sum _{i}{T}_{i}{\Phi}_{i}{\lambda}_{j}{C}_{j}$. (4)
Taking into account the equalities AΦ_{i} = χFΦ_{i} /K_{i} and χFΦ_{i} = (1 – β)χ_{p} FΦ_{i} + βχ_{d} FΦ_{i}, (χ is the total spectrum of fission neutrons, χ_{d} is the spectrum of delayed neutrons) system (4) can be transformed as
$\sum _{i}\frac{\partial {T}_{i}}{\partial t}\frac{1}{v}{\Phi}_{i}+\sum _{i}{T}_{i}\frac{\chi}{{K}_{i}}F{\Phi}_{i}+\sum _{i}{T}_{i}\Delta A{\Phi}_{i}=$
$=\chi \sum _{i}{T}_{i}(F+\Delta F){\Phi}_{i}{\chi}_{i}\beta \sum _{i}{T}_{i}(F+\Delta F){\Phi}_{i}+\sum _{j}{\lambda}_{j}{\chi}_{d,j}{C}_{j}+Q$, (5)
$\frac{\partial {C}_{j}}{\partial t}={\beta}_{j}(F+\Delta F)\sum _{i}{T}_{i}{\Phi}_{i}{\lambda}_{j}{C}_{j}$
.
We multiply the first equation of the resulting system by the eigenfunctions Φ^{+}_{k} of the adjoint eigenvalue problem
A ^{+}Φ^{+}_{k} = χF^{+}Φ^{+}_{i}/K_{k}, (6)
and the equations for delayed neutrons by χ_{d},_{j} Φ^{+}_{k}.
Taking into account the orthogonality of the eigenfunctions of the direct and adjoint problems
$\left.\right)"\; close=">">{\Phi}_{k}^{+},\chi F{\Phi}_{i}$(7)
(here γ_{0} is the fission neutrons importance, for k > 0 γ_{k} is a generalization of the concept of the fission neutrons importance for higher harmonics of the neutron flux) we obtain the system of equations
$\sum _{i}\frac{\partial {T}_{i}}{\partial t}\left.\right)"\; close=">">{\Phi}^{+}{}_{k},\frac{1}{v}{\Phi}_{i}$
$=\left(1{K}_{k}^{1}\right){T}_{k}\left.\right)"\; close=">">{\Phi}^{+}{}_{k},\chi F{\Phi}_{k}$
$+\left.\right)"\; close=">">{\Phi}^{+}{}_{k},Q$, (8)
$\frac{\partial \left.\right)"\; close=">">{\Phi}^{+}{}_{k},{\chi}_{d,j}{C}_{j}}{}\partial t$
$=\sum _{i}{T}_{i}\left.\right)"\; close=">">{\Phi}_{k}^{+},{\chi}_{d,j}{\beta}_{j}F{\Phi}_{i}$.
By entering parameters
ρ_{k} = 1 – 1/K_{k};
Λ_{k},_{i} = 〈Φ^{+}_{k}, Φ_{i}/v〉/〈Φ^{+}_{k}, χFΦ_{k}〉;
β_{eff},_{k},_{i},_{j} = 〈Φ^{+}_{k}, χ_{d},_{j}β_{j}FΦ_{i}〉/〈Φ^{+}_{k}, χFΦ_{k} 〉;
β_{eff},_{k},_{i} = 〈Φ^{+}_{k}, χ_{d},βFΦ_{i} 〉/〈Φ^{+}_{k}, χFΦ_{k}〉;
c_{k},_{j} = 〈Φ^{+}_{k}, χ_{d},_{j}C_{j}〉/〈Φ^{+}_{k}, χFΦ_{k}〉;
q_{k} = 〈Φ^{+}_{k}, Q〉/〈Φ^{+}_{k}, χFΦ_{k}〉;
δA_{k},_{i} = 〈Φ^{+}_{k}, ΔAΦ_{i}〉/〈Φ^{+}_{k}, χFΦ_{k}〉;
δF_{k},_{i} = 〈Φ^{+}_{k}, χΔFΦ_{i}〉/〈Φ^{+}_{k}, χFΦ_{k}〉;
δD_{k},_{i,j} = 〈Φ^{+}_{k}, χ_{d}β_{j}ΔFΦ_{i}〉/〈Φ^{+}_{k}, χFΦ_{k} 〉,
which are a generalization of the parameters of point kinetics, we bring the system to the form
$\mathbf{\Lambda}\frac{d\mathbf{T}}{dt}=(\mathbf{\rho}\mathbf{\beta})\mathbf{T}+\sum _{j}{\mathbf{\lambda}}_{j}{\mathbf{c}}_{j}+\mathbf{q}\delta \mathbf{A}\mathbf{T}+\delta \mathbf{FT}\sum _{j}\delta {\mathbf{D}}_{j}\mathbf{T}$
$\frac{d{\mathbf{c}}_{j}}{dt}={\mathit{\beta}}_{j}\mathbf{T}{\mathit{\lambda}}_{j}{\mathbf{c}}_{j}$ (9)
Here Λ is the matrix composed of elements Λ_{k},_{i}; ρ is the diagonal matrix, on the diagonal of which the elements ρ_{k} are located; β is the matrix composed of the elements β_{eff},_{k},_{i}; β_{j} is the matrix of elements β_{eff},_{k},_{i},_{j}; λ_{j} is the diagonal matrix with elements of λ_{j} on the diagonal; T is the vector of the amplitudes of the harmonics T_{k}; c_{j} is the vector of elements c_{k},_{j}; q is the matrix of elements q_{k}; δA is the matrix of elements δA_{k},_{i}; δF is the matrix of elements δF_{k},_{i}; δD_{j} is the matrix of elements δD_{k},_{i,j}.
The number of equations in system (9) is determined by the number of eigenfunctions used to represent the spatial dependence of the neutron flux N and the number of delayed neutron groups N_{d} and is equal to (N_{d} + 1)N.
In order to check the correctness of the equations obtained, a calculation was performed using the modal kinetics model of the transient process in a highpower BNtype reactor.
To calculate the coefficients included in system (9), and to solve this system, computer programs were developed.
The higher harmonics of the conditionally critical problem, necessary for calculating the coefficients of system (9), were calculated in the diffusion approximation using the modified neutronic calculation module of the REACTOR software package (
The parameters of system (9) were calculated by summing over the volume of the calculated region the values of the direct and adjoint flux harmonics, multiplied by the corresponding coefficients, i.e., the inverse neutron velocity, fission neutron source, delayed neutron source, external neutron source, deviations of the neutron transfer operator and the fission neutron source and delayed neutrons from a stationary state.
To solve system of equations (9), the software package for solving systems of ordinary differential equations ODEPACK (
The transient process consisted in the rise of the reactor power after a change in state, leading to the introduction of positive reactivity. In the initial state, the reactor was critical (matrices δA, δF and δD_{j} had zero elements) and operated at a power of 1 kW. The reactor was brought out of the critical state due to a 5% decrease in the sodium density in the active part of 18 central fuel assemblies within 10 seconds, the reactor parameters remained unchanged within 200 seconds, after which, within 10 seconds, the sodium density in the central fuel assemblies returned to the original value. The values of the elements of the matrices δA, δF and δD_{j} were calculated for the disturbed state of the reactor. It was assumed that the elements of the matrices δA, δF and δD_{j} when a disturbance is introduced into the reactor, change linearly with time, i.e., when a disturbance was introduced, they would linearly increase from zero values to values corresponding to the disturbed state; when the reactor returned to its initial state, they would linearly decrease to zero values.
Conditionally critical calculation of the disturbed state of the reactor showed that the introduced disturbance of the composition corresponded to an increase in reactivity Δρ_{0} = 4,25∙10^{–5}ΔK/K.
The transient process was calculated in two models: in the standard model of point kinetics with the assignment of the introduced reactivity according to the law
$\Delta \rho \left(t\right)=\left\{\begin{array}{l}0nbsp;nbsp;\text{nbsp;fornbsp;}t10nbsp;\mathrm{s}\text{nbsp;ornbsp;}t230nbsp;\mathrm{s};\\ \Delta {\rho}_{0}(t10)/10nbsp;nbsp;\text{nbsp;fornbsp;}10nbsp;\mathrm{s}t20nbsp;\mathrm{s};\\ \Delta {\rho}_{0}(230t)/10\text{nbsp;fornbsp;}10nbsp;\mathrm{s}t20nbsp;\mathrm{s};\\ \Delta {\rho}_{0}\text{nbsp;fornbsp;}20nbsp;\mathrm{s}t220nbsp;\mathrm{s}.\end{array}\right)$ (10)
and in modal kinetics model (9) using 10 eigenfunctions in expansion (1).
The reactor power calculated in the modal approximation (Fig.
The temporal behavior of all the higher harmonics of the neutron flux during the transient process is characterized by the following regularity: at the initial stage, there is a significant (by several orders of magnitude) increase in the ratio of their amplitudes to the amplitude of the fundamental harmonic T_{i} (t)/T_{0}(t). Then this ratio for all the harmonics is stabilized, which means that the reactor accelerates with a constant distribution of the neutron flux. After the reactor returns to its initial state, the higher harmonics damp, and the neutron flux distribution in the reactor is again established, corresponding to the fundamental harmonic of the conditionally critical problem (Fig.
The stability analysis was carried out for the system of equations of the spatial kinetics of an undisturbed reactor without feedbacks.
This system can be rearranged into the form
$\frac{d\mathbf{T}}{dt}={\mathbf{\Lambda}}^{1}(\mathbf{\rho}\mathbf{\beta})\mathbf{T}+{\mathbf{\Lambda}}^{1}\sum _{j}{\mathbf{\lambda}}_{j}{\mathbf{c}}_{j}+{\mathbf{\Lambda}}^{1}\mathbf{q}$, (11)
$\frac{d{\mathbf{c}}_{j}}{dt}={\mathit{\beta}}_{j}\mathbf{T}{\mathit{\lambda}}_{j}{\mathbf{c}}_{j}$.
According to the theory of ordinary differential equations, the stability of the solution to a system of linear differential equations with constant coefficients is determined by the spectrum of the system matrix.
In the case under consideration, the system matrix takes the form
$\left(\begin{array}{cccc}{\Lambda}^{1}(\rho \beta )& {\Lambda}^{1}{\lambda}_{1}& \cdots & {\Lambda}^{1}{\lambda}_{{N}_{d}}\\ {\beta}_{1}& {\lambda}_{1}& 0\cdots & 0\\ \vdots & \vdots & \vdots & \vdots \\ {\beta}_{{N}_{d}}& 0& 0\cdots & {\lambda}_{{N}_{d}}\end{array}\right)$. (12)
The absence in the spectrum of the matrix of eigenvalues with a positive real part is a criterion for the stability of the system of ordinary differential equations. If all the eigenvalues have negative real parts, the system has asymptotic stability, if among the eigenvalues there is a simple eigenvalue with zero real part, the system is Lyapunov stable (
To analyze the stability of the power distribution in the highpower BNtype reactor, parameters (9) and matrix elements (12) were calculated.
The stability analysis was carried out for two options, i.e., stationary operation at a power level close to the minimum controllable power, when the power of the subcritical reactor is maintained by its own neutron source associated with the decay of transuranic elements in the core, and for stationary operation in a critical state.
The thermal power of the reactor with the source was taken to be 1 kW. According to the estimates made, for the considered version of a highpower BN reactor, this power is achieved at a reactivity of ρ_{0} = –5,51 × 10^{–4}ΔK/K.
When we use five terms in expansion (1) and six groups of delayed neutrons, the dimension of system (10) is 35. Thus, matrix (12) has 35 eigenvalues. To solve the complete eigenvalue problem for this matrix, the LAPACK software package (
To test the effect on the eigenvalues of the system of kinetic equations of changes in the number of harmonics used to describe the spatial dependence of the neutron flux in expansion (1), the matrix spectrum of the system (11) was calculated using one term in the expansion, which corresponds to the approximation of point kinetics. Table
Spectra of matrices of the system of kinetic equations of a highpower BNtype reactor with a source
Point kinetics, s^{–1}  Spatial kinetics, s^{–1} 

–9.28278E–3  –9.28278E–3 
–1.68761E–2  –1.33580E–2 
–8.05185E–2  –1.33598E–2 
–2.01789E–1  –1.33948E–2 
–7.43827E–1  –1.33950E–2 
–2.78100E+0  –1.68761E–2 
–9.83543E+3  –2.99014E–2 
–2.99298E–2  
–3.05003E–2  
–3.05010E–2  
–8.05185E–2 
For the case of point kinetics, the number of eigenvalues is by one greater than the number of groups of delayed neutrons. In the case of expansion in harmonics, each eigenvalue of the system of point kinetics equations corresponds to five (in the considered case) eigenvalues, since the system includes the concentrations of precursor nuclei and the neutron flux expanded in the first five harmonics of the neutron flux. The eigenvalue minimum in absolute value of these five corresponds to the expansion of the flux and the concentration of precursor nuclei in terms of the fundamental harmonic and coincides with the corresponding eigenvalue of the point kinetics system. Thus, the first eigenvalue of the system written for the expansion with five terms coincides with the first eigenvalue of the system of point kinetics equations, the sixth one coincides with the second eigenvalue of the point kinetics system, the 11^{th} eigenvalue coincides with the third, etc.
As can be seen from the table, refining the description of the reactor kinetics does not lead to an increase in the maximum eigenvalue of the system of reactor kinetics equations.
For a reactor operating in a critical state at ρ_{0} = 0, the spectrum of the matrix of point kinetics system and seven maximum eigenvalues of the matrix of spatial kinetics system are given in Table
Spectra of matrices of the system of kinetic equations of a highpower BNtype reactor in a critical state
Point kinetics, s^{–1}  Spatial kinetics, s^{–1} 

0.0  0.0 
–1.5019E–2  –1.3356E–2 
–7.0303E–2  –1.3358E–2 
–1.8911E–1  –1.3395E–2 
–7.2721E–1  –1.3395E–2 
–2.7576E+0  –1.5019E–2 
–8.5260E+3  –2.9874E–2 
The maximum eigenvalue of kinetics system of a critical highpower BNtype reactor without feedbacks is zero; thus, its power distribution is Lyapunov stable (
Calculations of the decay of disturbances of the first and fourth higher harmonics of the neutron flux were performed. After a disturbance equal to 10 W is introduced, the amplitude of the harmonic decreases to a value close to zero in a time of the order of 10^{–4} s (Fig.
As shown by the calculation of the spectra of the matrices of the system of differential spatial kinetics equations for the subcritical and critical states of the highpower BNtype reactor without feedbacks, the stability criterion for the states considered is fulfilled, the neutron flux in the subcritical core of the BN reactor is asymptotically stable, the neutron flux in the critical core of the BN reactor is Lyapunov stable.
We plan to continue our work by adding feedbacks to the spatial kinetics model and substantiating the stability of the neutron flux in the highpower BN core with feedbacks.