The article discusses the neutron flux stability in the core of a high-power sodium-cooled fast reactor (of the BN-type) without feedbacks. The importance of this problem for high-power BN-type 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 high-power BN-type reactor without feedbacks. A computational model of the spatial kinetics of a high-power BN-type reactor has been developed in the modal approximation based on the representation of an unsteady flux as a sum of orthogonal functions multiplied by time-dependent 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 high-power BN-type reactor without feedbacks is stable. Test calculations have illustrated the damping of perturbations of the power distribution for a reactor in a critical state.

BN-type reactor, neutron flux stability, spatial kinetics, conditionally critical problem, eigenfunction expansion

Numerous safety studies of sodium-cooled 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 (Poplavsky et al. 2011). In order to provide the required thermal power of the core within these constraints, its diameter for a high-power BN-type reactor should be more than 4.3 m. The height/diameter ratio for such a core exceeds 5. However, the Nuclear Safety Rules (NP-082-07 2008) require proof that there are no oscillations in the neutron flux density, which is especially important for such ‘flat’ cores.

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’ (Akcasu 1971, Coughanowr 1991).

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) (Bibikov 2011).

If the solution φ (t) of the system satisfies not only the Lyapunov stability condition, but also the condition

$\underset{t\to \infty }{\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 (Bibikov 2011).

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) (March-Leuba and Blakeman 1991, Hashimoto 1993, Miro et al. 2002). The advantage of this method is that, with a complication of the system of equations describing the dynamics of the reactor that is insignificant in comparison with the point kinetics, it makes it possible to estimate the possibility of oscillations of the power distribution in the reactor. As a rule, the analysis of the stability of power distribution is limited to the use of the first few harmonics of the conditionally critical problem. Although this amount is insufficient for a correct description of the spatial kinetics of the reactor in the general case, such an approach is justified for stability analysis, since, as a rule, the onset of instability of the power distribution leads to the excitation, primarily, of the first harmonics corresponding to the largest values of the eigenvalues K_i (Takeuchi et al. 1994).

In this work, the stability of the power distribution of a high-power BN-type sodium-cooled 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 non-stationary 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 (Bell and Glesston 1974), the eigenfunctions of the discretized neutron diffusion equation form a complete system, which makes it possible to use the eigenfunction expansion of the conditionally critical diffusion equation for the total flux, similar to (1).

Let us consider the non-stationary 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^{’}}\sum _s\left(\mathbf{r},E^{’}\to E\right)\Phi \left(\mathbf{r},E^{’}\right)d{E}^{’}$

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

(7)

(here γ₀ 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

, (8)

.

By entering parameters

ρ_k = 1 – 1/K_k;

Λ_k,_i = 〈Φ⁺_k, Φ_i/v〉/〈Φ⁺_k, χFΦ_k〉;

β_eff,_k,_i,_j = 〈Φ⁺_k, χ_d,_jβ_jFΦ_i〉/〈Φ⁺_k, χFΦ_k 〉;

β_eff,_k,_i = 〈Φ⁺_k, χ_d,βFΦ_i 〉/〈Φ⁺_k, χFΦ_k〉;

c_k,_j = 〈Φ⁺_k, χ_d,_jC_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 high-power BN-type 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 (Gulevich et al. 2019). To solve the partial eigenvalue problem in the modified module, the implicitly restarted Arnoldi method was used, implemented in the ARPACK software package (Lehoucq et al. 1998).

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 (Hindmarsh 1983) was used, the method of the backward differentiation formula was applied.

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 Δρ₀ = 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}0\&nbsp;\&nbsp;\text{\&nbsp;for\&nbsp;}t<10\&nbsp;\mathrm{s}\text{\&nbsp;or\&nbsp;}t>230\&nbsp;\mathrm{s};\\ \Delta {\rho }_0(t-10)/10\&nbsp;\&nbsp;\text{\&nbsp;for\&nbsp;}10\&nbsp;\mathrm{s}<t<20\&nbsp;\mathrm{s};\\ \Delta {\rho }_0(230-t)/10\text{\&nbsp;for\&nbsp;}10\&nbsp;\mathrm{s}<t<20\&nbsp;\mathrm{s};\\ \Delta {\rho }_0\text{\&nbsp;for\&nbsp;}20\&nbsp;\mathrm{s}<t<220\&nbsp;\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. 1) coincided with the power calculated in the point kinetics model with high accuracy (the difference did not exceed 0.14%), which indicates the correctness of the developed model and its software implementation.

Figure 1. 

Change in the reactor power in the test calculation.

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₀(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. 2). Since the disturbance introduced into the reactor is close to axisymmetric, the ratio |T_i (t)|/T₀(t) reaches the highest value (8∙10^–4) for the fifth harmonic, which is due to the fact that it has a radial character. For the first harmonic, which has an azimuthal character, the maximum value of the ratio |T_i (t)|/T₀(t) is 2,7∙10^–7.

Figure 2. 

Relative amplitudes (|T_i (t)|/T₀(t)) of the fourth, fifth, seventh and ninth harmonics of the neutron flux.

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\\ ⋮& ⋮& ⋮& ⋮\\ {\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 (Bibikov 2011).

To analyze the stability of the power distribution in the high-power BN-type 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 high-power BN reactor, this power is achieved at a reactivity of ρ₀ = –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 (Anderson et al. 1999) was involved, which implements the solution of the complete eigenvalue problem using the QR algorithm. The calculation showed that all the eigenvalues of the spatial kinetics matrix for the considered variant are real and less than zero, the maximum eigenvalue being –9.28 × 10^–3s^–1. Thus, the system of spatial kinetics equations for a high-power BN-type reactor with a source is asymptotically stable (Bibikov 2011).

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 1 contains the eigenvalues of the kinetic matrix using one term in the expansion and 11 maximum eigenvalues of the kinetic matrix using five expansion terms.

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, the spectrum of the matrix of point kinetics system and seven maximum eigenvalues of the matrix of spatial kinetics system are given in Table 2.

The maximum eigenvalue of kinetics system of a critical high-power BN-type reactor without feedbacks is zero; thus, its power distribution is Lyapunov stable (Bibikov 2011). The zero eigenvalue is associated with the zero flux harmonic corresponding to zero reactivity, and expresses the fact that the power of the critical reactor without feedbacks is not determined. Therefore, when the power of the fundamental harmonic is disturbed at the end of the transient process, the reactor is stabilized at a new power level that does not coincide with the power level before the disturbance. The disturbances of higher harmonics attenuate and return to zero values at the end of the transient process.

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. 3).

Figure 3. 

Decay of disturbances of the first and fourth higher harmonics in a high-power BN-type reactor.

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 high-power BN-type 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 high-power BN core with feedbacks.