Print
Assessment of the neutron flux stability in a high power BN-type reactor in terms of modal spatial kinetics*
expand article infoIvan V. Tormyshev, Andrey V. Gulevich, Vladimir A. Yeliseev, Victor Yu. Stogov
‡ JSC SSC RF-IPPE n.a. A.I. Leypunsky, Obninsk, Russia
Open Access

Abstract

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.

Keywords

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

Introduction

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

limtx(t)-φ(t)=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 Ki (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.

Deriving a system of reactor kinetics equations when expanding a non-stationary flux in eigenfunctions of a conditionally critical problem

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:

Φ(r,Ω,E,t)=iTi(t)Φi(r,Ω,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)

1vΦt+(A+ΔA)Φ=χp(1-β)(F+ΔF)Φ+jλjχd,jCj

Cjt=βj(F+ΔF)Φ-λjCj (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Φ(r,E)=-D(r,E)Φ(r,E)+Σtot(r,E)Φ(r,E)-Esr,EEΦr,EdE

F is the neutron fission operator of the stationary conditionally critical problem describing the reactor critical state, defined as

FΦ(r,E)=EvΣf(r,E)Φ(r,E)dE

ΔA are ΔF are the operator disturbances that brought the reactor out of the critical state; Q is the external source; Cj 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 /Ki, (3)

iTit1vΦi+iTi(A+ΔA)Φi=χp(1-β)iTi(F+ΔF)Φijλjχd,jCj+Q,

Cjt=βj(F+ΔF)iTiΦi-λjCj. (4)

Taking into account the equalities AΦi = χFΦi /Ki 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

iTit1vΦi+iTiχKiFΦi+iTiΔAΦi=

=χiTi(F+ΔF)Φi-χiβiTi(F+ΔF)Φi+jλjχd,jCj+Q, (5)

Cjt=βj(F+ΔF)iTiΦi-λjCj

.

We multiply the first equation of the resulting system by the eigenfunctions Φ+k of the adjoint eigenvalue problem

A +Φ+k = χF+Φ+i/Kk, (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

Φk+,χFΦi=γk&nbsp;for&nbsp;k=i,0&nbsp;for&nbsp;ki;(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

iTitΦ+,k1vΦi=

=1-Kk-1TkΦ+,kχFΦk-iTiΦk+,χdβFΦi+jλjΦ+,kχd,jCj+

+Φ+,kQ-iTiΦ+,kΔAΦi+iTiΦk+,χΔFΦi-iTiΦk+,χdβΔFΦi, (8)

Φ+,kχd,jCjt=

=iTiΦk+,χd,jβjFΦi+iTiΦk+,χd,jβjΔFΦi-λjΦk+,χd,jCj.

By entering parameters

ρk = 1 – 1/Kk;

Λk,i = 〈Φ+k, Φi/v〉/〈Φ+k, χFΦk〉;

βeff,k,i,j = 〈Φ+k, χd,jβjFΦi〉/〈Φ+k, χFΦk 〉;

βeff,k,i = 〈Φ+k, χdFΦi 〉/〈Φ+k, χFΦk〉;

ck,j = 〈Φ+k, χd,jCj〉/〈Φ+k, χFΦk〉;

qk = 〈Φ+k, Q〉/〈Φ+k, χFΦk〉;

δAk,i = 〈Φ+k, ΔAΦi〉/〈Φ+k, χFΦk〉;

δFk,i = 〈Φ+k, χΔFΦi〉/〈Φ+k, χFΦk〉;

δDk,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

ΛdTdt=(ρ-β)T+jλjcj+q-δAT+δFT-jδDjT

dcjdt=βjT-λjcj (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 Tk; cj is the vector of elements ck,j; q is the matrix of elements qk; δA is the matrix of elements δAk,i; δF is the matrix of elements δFk,i; δDj is the matrix of elements δDk,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 Nd and is equal to (Nd + 1)N.

Test calculations in the modal kinetics model of a BN-type reactor

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 δDj 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 δDj were calculated for the disturbed state of the reactor. It was assumed that the elements of the matrices δA, δF and δDj 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

Δρ(t)=0&nbsp;&nbsp;&nbsp;for&nbsp;t<10&nbsp;s&nbsp;or&nbsp;t>230&nbsp;s;Δρ0(t-10)/10&nbsp;&nbsp;&nbsp;for&nbsp;10&nbsp;s<t<20&nbsp;s;Δρ0(230-t)/10&nbsp;for&nbsp;10&nbsp;s<t<20&nbsp;s;Δρ0&nbsp;for&nbsp;20&nbsp;s<t<220&nbsp;s. (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 |Ti (t)|/T0(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 |Ti (t)|/T0(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 |Ti (t)|/T0(t) is 2,7∙10–7.

Figure 2. 

Relative amplitudes (|Ti (t)|/T0(t)) of the fourth, fifth, seventh and ninth harmonics of the neutron flux.

Analysis of the stability of a system of spatial kinetics equations of a high-power BN-type reactor

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

dTdt=Λ-1(ρ-β)T+Λ-1jλjcj+Λ-1q, (11)

dcjdt=βjT-λjcj.

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

Λ-1(ρ-β)Λ-1λ1Λ-1λNdβ1-λ100βNd00-λNd. (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 ρ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 (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.

Table 1.

Spectra of matrices of the system of kinetic equations of a high-power BN-type 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 11th 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 2.

Table 2.

Spectra of matrices of the system of kinetic equations of a high-power BN-type 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 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.

Conclusion

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.

References

  • Akcasu Z (1971) Mathematical Methods in Nuclear Reactor Dynamics. Academic Press, New York, 460 pp.
  • Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J, Du Croz J, Greenbaum A, Hammarling S, McKenney A, Sorensen D (1999) LAPACK Users’ Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, 407 pp. https://doi.org/10.1137/1.9780898719604
  • Bell J, Glesston S (1974) Nuclear Reactor Theory. Atomizdat Publ., Мoscow, 494 pp. [in Russian]
  • Bibikov YuN (2011) Course in Ordinary Differential Equations. Lan Publ., St. Petersburg, 304 pp. [in Russian]
  • Coughanowr DR (1991) Process Systems Analysis and Control. McGraw-Hill, New York, 566 pp.
  • Gulevich AV, Dolgikh VP, Eliseev VA, Peregudova OO, Rozhihin EV, Semenov MYu, Stogov VYu, Tormyshev IV (2019) Calculation of the Dominant Ratio for a BN-type Reactor Core. VANT. Ser.: Yaderno-Reaktornye Konstanty [VANT. Series: Nuclear Reaction Constants] 4: 50–54. [in Russian]
  • Hindmarsh AC (1983) ODEPACK, A Systematized Collection of ODE Solvers. In: Stepleman RS et al. (Eds) Scientific Computing. North-Holland, Amsterdam, (Vol. 1 of IMACS Transactions on Scientific Computation), 55–64.
  • Lehoucq RB, Sorensen DC, Yang C (1998) ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. Society for Industrial and Applied Mathematics, Philadelphia, PA, 137 pp. https://doi.org/10.1137/1.9780898719628
  • March-Leuba J, Blakeman ED (1991) A Mechanism for Out-of-Phase Power Instabilities in Boiling Water Reactors. Nuclear Science and Engineering 107(2): 173–179. https://doi.org/10.13182/NSE91-A15730
  • Miro R, Ginestarb D, Verdu G, Hennigc D (2002) A Nodal Modal Method for the Neutron Diffusion Equation. Application to BWR Instabilities Analysis. Annals of Nuclear Energy 29: 1171–1194. https://doi.org/10.1016/S0306-4549(01)00103-7
  • NP-082-07 (2008) Federalnaya Sluzhba po Ekologicheskomu, Tekhnologicheskomu i Atomnomu Nadzoru. Nuclear Safety Rules for Reactor Installations of Nuclear Power Plants. Yadernaya i Radiatsionnaya Bezopasnost [Nuclear and Radiation Safety] 1: 52–77. [in Russian]
  • Poplavsky VM, Matveev VI, Eliseev VA, Kuznetsov IA, Volkov AV, Shvetsov YuE, Khomyakov YuS, Tsiboula AM (2011) Studies on Influence of Sodium Void Reactivity Effect on the Concept of the Core and Safety of Advanced Fast Reactor. Journal of Nuclear Science and Technology 48(4): 538–546. https://doi.org/10.1080/18811248.2011.9711731
  • Takeuchi Y, Takigawa Y, Uematsu H (1994) A Study on Boiling Water Reactor Regional Stability from the Viewpoint of Higher Harmonics. Nuclear Technology 106(3): 300–314. https://doi.org/10.13182/NT94-A34960

* Russian text published: Izvestiya vuzov. Yadernaya Energetika (ISSN 0204-3327), 2021, n. 2, pp. 39–49.