Generalization of the Analytical Exponential Model for Homogeneous Reactor Kinetics Equations

Abstract

Mathematical form for two energy groups of three-dimensional homogeneous reactor kinetics equations and average one group of the precursor concentration of delayed neutrons is presented. This mathematical form is called “two energy groups of the point kinetics equations.” We rewrite two energy groups of the point kinetics equations in the matrix form. Generalization of the analytical exponential model (GAEM) is developed for solving two energy groups of the point kinetics equations. The GAEM is based on the eigenvalues and the corresponding eigenvectors of the coefficient matrix. The eigenvalues of the coefficient matrix are calculated numerically using visual FORTRAN code, based on Laguerre’s method, to calculate the roots of an algebraic equation with real coefficients. The eigenvectors of the coefficient matrix are calculated analytically. The results of the GAEM are compared with the traditional methods. These comparisons substantiate the accuracy of the results of the GAEM. In addition, the GAEM is faster than the traditional methods.

1. Introduction

In real nuclear reactor problems that involve safety considerations, the coefficients of the space-time neutron diffusion equations depend upon the neutron flux density, precursor concentration of delayed neutrons groups, space, time, and temperature feedback. The numerical methods and the analytical approximations for solving the space-time neutron diffusion equations in the nuclear reactor have been of interest in nuclear physics and reactor design.

Three-dimensional nodal mixed dual method was used for the nuclear reactor kinetics equations with improved quasistatic model and the semiimplicit scheme to solve the precursor equations [1]. A rigorous weight function in the quasistatic method was presented to solve the time-dependent multigroup diffusion equations for subcritical systems with external sources [2]. Parallelised Krylov methods were applied to the improved quasistatic approach in addition to the direct, implicit time difference, approach for solving space-time dependent multigroup neutron diffusion equations [3]. The nodal diffusion method was presented to solve space-time neutron kinetics using the finite element, primal and mixed hybrid nodal methods [4]. A one-step implicit method and a nodal modal method were studied for the solution of the time-dependent neutron diffusion equations which-based on a hexagonal spatial mesh [5–7]. An adiabatic kinetics model was developed for a boiling water reactor (BWR) core simulator AETNA [8]. Mathematical modeling of the space-time kinetics phenomena in advanced heavy water reactor, a 920 MW thermal, vertical pressure tube type thorium-based nuclear reactor, was presented using nodal modal method [9]. Class of Padé and cut-product approximations were applied to solve the two energy groups space-time nuclear reactor kinetics equations in two dimensions [10]. Adaptive Matrix Formation (AMF) method was presented and applied to homogeneous, symmetric heterogeneous, and nonsymmetric heterogeneous reactors in the cases of two and three dimensions [11]. The generalized Runge-Kutta method was developed for solving the multigroup, multidimensional, static, and transient neutron diffusion kinetics equations [12]. Computation accuracy and efficiency of a power series analytic method were presented for the time-space-dependent neutron diffusion equations with adiabatic heat up and doppler feedback [13]. The multigroup integro-differential equations of the neutron diffusion kinetics were presented and solved numerically in multislab geometry with the use of the progressive polynomial approximation [14].

The different kinds of Padé approximations were applied to the nonlinear point kinetics equations including temperature feedback via analytical inversion method [15, 16]. The straightforward recurrence relation of a power series was presented to solve the point kinetics equations in the presence of delayed neutrons with Newtonian feedback [17]. The stability of generalized Runge-Kutta method of order four was studied for the stiff system point kinetics equations with temperature feedback reactivity [18, 19]. The analytical exponential model was generalized to solve the point kinetics equations of beryllium and heavy water moderated reactors with step, ramp, and temperature feedback reactivities [20]. Power series solution (PWS) was applied for solving the point kinetics equations with lumped model temperature feedback [21]. Recently, new analytical solution, based on the roots of inhour equation, the eigenvalues of the coefficient matrix, was presented and applied for the point kinetics equations with the constant, ramp, and temperature feedback reactivities [22].

In this work, the two energy groups of the point kinetics equations and the average one group of delayed neutron are introduced. The generalization of the analytical exponential model is calculated using the eigenvalues and the eigenvectors of the coefficient matrix. The eigenvalues and the corresponding eigenvectors of the coefficient matrix are calculated. The results of the generalization of the analytical exponential model are compared with the results of the traditional methods.

2. Mathematical Model

The two energy groups neutron diffusion equations governing the dynamic groups diffusion neutron flux and average one group of the precursor concentration of delayed neutrons behavior are written as follows [23–25]: 1𝑣₁𝜕Φ𝜕𝑡¹(𝑟,𝑡)=𝐷₁∇²Φ₁(𝑟,𝑡)−Σ_𝑎1+Σ_𝑠12Φ₁(𝑟,𝑡)+𝜈Σ_𝑓11−𝛽Φ₁(𝑟,𝑡)+𝜈Σ_𝑓21−𝛽Φ₂()(),1𝑟,𝑡+𝜆𝐶𝑟,𝑡𝑣₂𝜕Φ𝜕𝑡²(𝑟,𝑡)=𝐷₂∇²Φ₂(𝑟,𝑡)−Σ_𝑎2Φ₂(𝑟,𝑡)+Σ_𝑠12Φ₁(),𝜕𝑟,𝑡𝐶()𝜕𝑡𝑟,𝑡=𝛽𝜈Σ_𝑓1Φ₁(𝑟,𝑡)+𝜈Σ_𝑓2Φ₂(𝑟,𝑡)−𝜆𝐶(𝑟,𝑡), (2.1) where Φ₁(𝑟,𝑡) and Φ₂(𝑟,𝑡) are fast and thermal neutron flux, 𝐶(𝑟,𝑡) is average one group of the precursor concentration of delayed neutrons, 𝐷₁ and 𝐷₂ are fast and thermal diffusion coefficients, Σ_𝑎1 and Σ_𝑎2 are fast and thermal absorption cross-sections, Σ_𝑓1 and Σ_𝑓2 are fast and thermal fission cross-sections, Σ_𝑠12 is the scattering cross-section from fast-to-thermal neutron, 𝜈 is the neutron fission, 𝑣₁ and 𝑣₂ are fast and thermal neutron speed, 𝜆 is decay constant of average one group of delayed neutrons, and 𝛽 is the fraction of delayed neutrons.

The neutron flux and the precursor concentration of delayed neutrons can be written as separable functions of space and time as Φ₁(𝑟,𝑡)=𝑣₁𝑛₁(𝑡)Ψ(𝑟),Φ₂(𝑟,𝑡)=𝑣₂𝑛₂(𝑡)Ψ(𝑟),𝐶(𝑟,𝑡)=𝑐(𝑡)Ψ(𝑟), (2.2) where 𝑛₁(𝑡) and 𝑛₂(𝑡) are the fast and thermal neutron density, 𝑐(𝑡) is the precursor concentration density of delayed neutrons and Ψ(𝑟) is the fundamental function which can be determined from the following: ∇²Ψ(𝑟)+𝐵²Ψ(𝑟)=0, (2.3) where 𝐵² is the material buckling.

Substituting (2.2) and (2.3) into (2.1) yields 𝑑𝑛₁(𝑡)𝑑𝑡=−𝐷₁𝐵²𝑣₁𝑛₁(𝑡)−Σ_𝑎1+Σ_𝑠12𝑣₁𝑛₁(𝑡)+𝜈Σ_𝑓11−𝛽𝑣₁𝑛₁(𝑡)+𝜈Σ_𝑓21−𝛽𝑣₂𝑛₂(𝑡)+𝜆𝑐(𝑡),𝑑𝑛₂(𝑡)𝑑𝑡=−𝐷₂𝐵²𝑣₂𝑛₂(𝑡)−Σ_𝑎2𝑣₂𝑛₂(𝑡)+Σ_𝑠12𝑣₁𝑛₁(𝑡),𝑑𝑐(𝑡)𝑑𝑡=𝛽𝜈Σ_𝑓1𝑣₁𝑛₁(𝑡)+𝜈Σ_𝑓2𝑣₂𝑛₂(𝑡)−𝜆𝑐(𝑡). (2.4) Let us consider that 𝑙₁=1/𝑣₁𝜈Σ_𝑓1 and 𝑙₂=1/𝑣₂𝜈Σ_𝑓2 are fast and thermal generation time of neutrons, 𝐿²₁=𝐷₁/Σ_𝑎1 and 𝐿²₂=𝐷₂/Σ_𝑎2 are fast and thermal diffusion length, 𝑘₁=𝜈Σ_𝑓1/Σ_𝑎1[1+𝐿²₁𝐵²] and 𝑘₂=𝜈Σ_𝑓2/Σ_𝑎2[1+𝐿²₂𝐵²] are fast and thermal multiplication factor, 𝜌₁=(𝑘₁−1)/𝑘₁ and 𝜌₂=(𝑘₂−1)/𝑘₂ are fast and thermal reactivities, 𝜅=𝑣₁Σ_𝑠12, 𝜇₁=𝛽/𝑙₁ and 𝜇₂=𝛽/𝑙₂.

Equations (2.4) are rewritten in the following form: 𝑑𝑛₁(𝑡)=𝜌𝑑𝑡𝑙¹₁−𝜇₁−𝜅𝑛₁(𝑡)+1𝑙₂−𝜇₂𝑛₂(𝑡)+𝜆𝑐(𝑡),𝑑𝑛₂(𝑡)𝑑𝑡=𝜅𝑛₁(𝑡)+𝜌₂−1𝑙₂𝑛₂(𝑡),𝑑𝑐(𝑡)𝑑𝑡=𝜇₁𝑛₁(𝑡)+𝜇₂𝑛₂(𝑡)−𝜆𝑐(𝑡). (2.5) Aboanber [26] presented two-energy group two-point kinetic model of reflected reactor. So, (2.5) can be called “two energy groups of the point kinetics equations”.

The initial values of this system, 𝑑𝑛₁(0)/𝑑𝑡=𝑑𝑛₂(0)/𝑑𝑡=𝑑𝑐(0)/𝑑𝑡=0, take the following form 𝑛₁(0)=1−𝜌₂𝑙₂𝜅𝑛₂(0),𝑐(0)=𝜇₁𝜆𝑛₁(0)+𝜇₂𝜆𝑛₂(0), (2.6) and the material buckling takes the form: 𝐵²Σ=−^𝑎1+Σ_𝑠12−𝜈Σ_𝑓12𝐷₁−Σ_𝑎22𝐷₂+Σ_𝑎1+Σ_𝑠12−𝜈Σ_𝑓12𝐷₁+Σ_𝑎22𝐷₂²−Σ_𝑎2Σ_𝑎1+Σ_𝑠12−𝜈Σ_𝑓1−Σ_𝑠12𝜈Σ_𝑓2𝐷₁𝐷₂. (2.7)

3. Generalization of the Analytical Exponential Model

Let us rewrite (2.5) in matrix form as follows 𝑑Ψ(𝑡)(𝑡)𝑑𝑡=𝐀Ψ, (3.1) where Ψ(𝑡)=𝑛₁(𝑡)𝑛₂(𝑡)𝑐(𝑡)𝜇,𝐀=−𝛼𝜁𝜆𝜅−𝜂0₁𝜇₂,−𝜆𝛼=𝜇₁𝜌+𝜅−¹𝑙₁1,𝜁=𝑙₂−𝜇₂,𝜂=1−𝜌₂𝑙₂. (3.2) Using the exponential integrator [27, 28], we can determine the general solution of (3.1) as follows: Ψ(𝑡)=𝐵exp^𝑡_𝑡0𝐀𝑑𝑡, (3.3) where 𝐵 is the integral constant.

This solution includes the exponential function of the matrix 𝐀 which can be calculated using generalization of the analytical exponential model [20, 29] as Ψ(𝑡)=³_𝑘=1𝐵_𝑘exp𝜔_𝑘𝑡𝐔_𝑘, (3.4) where 𝐵_𝑘 are the constants, 𝜔_𝑘 are the eigenvalues of the matrix 𝐀 and 𝐔_𝑘 are the corresponding eigenvectors of the matrix 𝐀.

Using the initial condition Ψ(0)=Ψ₀ to determine the constants 𝐵_𝑘=𝐕^𝑇_𝑘Ψ₀, (3.4) becomes as follows Ψ(𝑡)=³_𝑘=1exp𝜔_𝑘𝑡𝐔_𝑘𝐕^𝑇_𝑘Ψ₀, (3.5) where, 𝐕_𝑘 are the eigenvectors of matrix 𝐀^𝑇 which corresponding to eigenvalues 𝜔_𝑘 and satisfies the normalization condition 𝐔^𝑇_𝑘𝐕_𝑘=1. The eigenvectors 𝐔_𝑘 and 𝐕_𝑘 of the coefficient matrix 𝐀 and 𝐀^𝑇, respectively, are calculated analytically in the next subsections.

Equation (3.5) gives good results only for a constant matrix 𝐀. But for variable matrix 𝐀, the solution of the equation (3.1) can be takes the form Ψ𝑡_𝑚+1=³_𝑘=1expℎ𝜔_𝑘𝐔_𝑘𝐕^𝑇_𝑘Ψ𝑡_𝑚, (3.6) where, the values of the eigenvalues 𝜔_𝑘 and the corresponding eigenvectors 𝐔_𝑘 and 𝐕_𝑘 are calculated at time 𝑡_𝑚+(ℎ/2), and ℎ=𝑡_𝑚+1−𝑡_𝑚 is the time step interval [𝑡_𝑚,𝑡_𝑚+1].

Equations (3.5) and (3.6) represent the solutions of the two energy groups of point kinetics equation with average one group of delayed neutrons. These solutions represent a generalization of the analytical exponential model for step and variable coefficient matrix 𝐀, respectively.

3.1. Eigenvalues of the Matrix A

The eigenvalues of the matrix 𝐀 are determined from the following algebraic equation: 𝜔³+𝜔²𝛼+𝜂+𝜆+𝜔𝛼𝜂+𝛼𝜆+𝜂𝜆−𝜅𝜁−𝜆𝜇₁+𝛼𝜂𝜆−𝜅𝜁𝜆−𝜅𝜆𝜇₂−𝜂𝜆𝜇₁=0. (3.7) The equation (3.7) is called “inhour equation.” This equation has a 3rd degree algebraic equation. To determine the roots of (3.7), we use the visual FORTRAN code, based on Laguerre’s numerical method, to calculate the roots of an algebraic equation with real coefficients.

3.2. Eigenvectors of the Matrix A

The eigenvectors of the matrix 𝐀 are determined analytically as follows: 𝐔_𝑘=1𝜅𝜔_𝑘𝜇+𝜂₁𝜔_𝑘+𝜂+𝜅𝜇₂𝜔_𝑘+𝜆𝜔_𝑘+𝜂. (3.8)

3.3. Eigenvectors of the Matrix A^T

The eigenvectors of the matrix 𝐀^𝑇 are determined analytically as follows: 𝐕_𝑘=11+𝜅𝜁𝜔_𝑘+𝜆+𝜅𝜆𝜇₂/𝜔_𝑘+𝜆𝜔_𝑘+𝜂²+𝜆𝜇₁𝜔_𝑘+𝜂+𝜅𝜆𝜇₂/𝜔_𝑘+𝜆²𝜔_𝑘×1𝜁𝜔+𝜂_𝑘+𝜆+𝜆𝜇₂𝜔_𝑘+𝜆𝜔_𝑘𝜆+𝜂𝜔_𝑘,+𝜆 (3.9) where (1+(𝜅𝜁(𝜔_𝑘+𝜆)+𝜅𝜆𝜇₂)/(𝜔_𝑘+𝜆)(𝜔_𝑘+𝜂)²+(𝜆𝜇₁(𝜔_𝑘+𝜂)+𝜅𝜆𝜇₂)/(𝜔_𝑘+𝜆)²(𝜔_𝑘+𝜂))⁻¹ is the normalization factor.

4. Results and Discussions

To check the accuracy of the results of the mathematical model of two energy groups of the homogeneous reactor kinetics equations, it is applied to the three-dimensional homogeneous reactor. This reactor is a bare homogeneous cube of side length 200 cm on each side with two neutron groups and one precursor group. The parameters of the three-dimensional homogeneous reactor are taken from [11] and shown in Table 1. The mathematical model of two energy groups of the point kinetics equations are solved by the generalization of the analytical exponential model (GAEM). Two different types of problems; step and ramp reactivities, are applied.

Table 1 

Parameters of three-dimensional homogeneous reactor [11].

4.1. Step Reactivity

The thermal neutron flux at the center point of the two energy groups of the homogeneous reactor kinetics equations with a positive step reactivity, ΔΣ_𝑎2=−0.0000369, is calculated in Table 2. The results of the GAEM for the two energy groups of the point kinetics equations have an accuracy comparable with the results of the adaptive matrix formation (AMF) [11], 3DKIN [30] and the exact solution for the space-time nuclear reactor kinetics equations.

Table 2 

Thermal neutron flux at the center point of three-dimensional homogeneous reactor with a positive step reactivity.

The fast and thermal neutron flux of the two energy groups of the point kinetics equations with a negative step reactivity, ΔΣ_𝑎2=+0.0000369, is calculated in Table 3. This table substantiates the accuracy of the results of the GAEM compared with the results of the AMF [11].

Table 3 

Neutron flux at the center point of three-dimensional homogeneous reactor with a negative step reactivity.

Figure 1. shows the fast and thermal neutron flux at three values for positive step reactivities, ΔΣ_𝑎2=−0.0000369,−2(0.0000369) and −3(0.0000369). And, the fast and thermal neutron flux at three types of negative step reactivities, ΔΣ_𝑎2=+0.0000369,+2(0.0000369) and +3(0.0000369), are shown in Figure 2.

Figure 1 

Fast and thermal neutron flux at a positive step reactivity.

Figure 2 

Fast and thermal neutron flux at a negative step reactivity.

The calculation times, CPU times, of the GAEM to the AMF method are 0.43 s: 0.79 s at the same time step using personal desktop computer (processor: Core 2 Duo, 2.93 GHz). This means that the GAEM for two energy groups of the point kinetics equations is faster than the AMF for the space-time nuclear reactor kinetics equations.

4.2. Ramp Reactivity

To substantiate the accuracy and the efficiency of the generalization of the analytical exponential model for solving two energy groups of the point kinetics equations with linear ramp reactivity, it is applied to three dimensional homogeneous reactor. The fast and thermal neutron flux for a positive ramp reactivity, ΔΣ_𝑎2=−(0.0000369/0.4)𝑡, are calculated in Table 4. This table is demonstrates the comparison between the results of the GAEM and AMF.

Table 4 

Neutron flux at the center point of three dimensional homogeneous reactor with a positive ramp reactivity.

The fast and thermal neutron flux of the two energy groups of the point kinetics equations with a negative ramp reactivity, ΔΣ_𝑎2=+(0.0000369/0.4)𝑡, is calculated in Table 5.

Table 5 

Neutron flux at the center point of three dimensional homogeneous reactor with a negative ramp reactivity.

The fast and thermal neutron flux at different types of positive ramp reactivities, ΔΣ_𝑎2=−(0.0000369𝑡),−2(0.0000369𝑡) and −3(0.0000369𝑡), is shown in Figure 3. And also, the fast and thermal neutron flux at different types of negative ramp reactivities, ΔΣ_𝑎2=+(0.0000369𝑡),+2(0.0000369𝑡) and +3(0.0000369𝑡), are shown in Figure 4.

Figure 3 

Fast and thermal neutron flux at a positive ramp reactivity.

Figure 4 

Fast and thermal neutron flux at a negative ramp reactivity.

For ramp reactivity, The calculation times, CPU times, of the GAEM to the AMF method are 0.45 s: 0.81 s. This means, the GAEM is faster than the AMF.

5. Conclusions

Two energy groups of the point kinetics equations are represented a mathematical form for three dimensional homogeneous reactor kinetics equations and average one group of precursor delayed neutrons. The mathematical form was rewritten in the matrix form. Generalization of the analytical exponential model (GAEM) was developed to solve the new mathematical matrix form for three dimensional homogeneous reactor. This GAEM based on the eigenvalues and the corresponding eigenvectors of the coefficient matrix.

The results of GAEM were compared with the AMF, 3DKIN, and the exact solution for three dimensional homogeneous reactor. This comparison substantiated the accuracy of the results of the GAEM for different types of reactivities. In addition, the GAEM for the mathematical form, two energy groups of the point kinetics equations, is faster than any other traditional methods. So, this mathematical form may be one of the most important models for the neutron diffusion equations in the nuclear reactors.