ABSTRACT
A new approach to solve the inverse equation of point kinetics is presented. The neutron density history can be considered as a series of infinite Bernoulli numbers. In order to simplify the problem, the proposed method uses an approximation which considers only two derivatives of the neutron population density. The reactivity is calculated for different forms of neutron population density and for different time steps. Due to its high precision, the results obtained suggest that the proposed method can be used in a real-time digital reactivity meter.
1. Introduction
In order to safely control a nuclear power station, it is necessary to know the value of the reactivity with the highest precision, this value is known from the population density of neutrons that depends on the fission interaction of different materials [Citation1].
A great variety of scientific works have been carried out around the calculation of the reactivity, each of which proposes a different solution method, in the quest for both a better approximation and low computational cost. The reported works [Citation2–Citation7] determine the reactivity and the neutron source, taking into account sub-critical conditions for the neutron population density. Another method, that determines the reactivity without the dependence on the history of neutron population density was proposed [Citation8], where the reactivity term was modified to solely depend on the derivatives of the neutron population. Later on, a method that uses the discrete version of the Laplace transform – which was implemented in a non-recursive way in what is known as the FIR filter – was proposed [Citation9]. This method allows the discretisation of the integral of the history of neutron population density in a convolution sum, but with a significant computational cost. Another method that uses a Lagrangian polynomial form [Citation10] with an accuracy of was proposed, but with the drawback that if the order of the polynomial is increased, the computational cost is also increased. Subsequently, a method based on the Euler–Maclaurin formula was published where the first Bernoulli number was used as an approximation [Citation11], this approximation increased the accuracy in the calculation of nuclear reactivity. Recently, a paper with a very good accuracy in the calculation of reactivity was reported, it is based on a matrix formulation for the calculation of reactivity [Citation12]. For this work, a new approach is proposed for the integral of the history of neutron population density by means of an infinite sum containing the Bernoulli numbers.
This work’s motivation mainly points out to state that there are currently two methods with very good accuracy for the calculation of reactivity and they were named above, such as the method of derivatives (DM) [Citation8] and the matrix formulation method (MF) [Citation12]. The latter method demonstrated superiority in accuracy over the former, because it has no restrictions on the shape of neutron density, but has a higher computational cost since it requires multiplication of matrices. However, the matrix formulation method (MF) is accurate and fast enough to be implemented in real-time digital reactivity meters.
The objective of this work is to provide an alternative method to the matrix formulation (MF) method that has equal or greater accuracy in the calculation of reactivity, although the proposed method uses restrictions in the form of neutron density as in the derivative method (DM). To validate the proposed method, different numerical experiments with different forms of neutron density distributed in eight cases were performed.
This article is organized in the following way: in Section 2 the inverse equation of point kinetics is deduced. In Section 2.1 the reactivity is calculated using the first two Bernoulli numbers. In Section 3 the proposed method for reactivity with the generalisation of the Bernoulli numbers is developed; in Section 4 the results are shown, and finally, in Section 5 the conclusions are given.
2. Inverse equation of point kinetics
The equations that model the dynamics of a nuclear reactor are known as the equations of point kinetics, which can be obtained from the neutron diffusion equation. These equations describe the time evolution of the neutron distribution and the concentration of delayed neutron precursors [Citation1]:
With initial conditions:
where is the concentration of the ith group of delayed neutron precursors, is the density of the neutron population, is the reactivity, is the neutron generation time, is the ith fraction of delayed neutrons, is the total effective fraction of delayed neutrons and is the decay constant of the ith group of delayed neutron precursors.
It is possible to solve EquationEquation (2)(2) (2) using the integrating factor method and making use of the initial conditions (3) and (4), thus obtaining an expression for the concentration of precursors of the ith group of delayed neutrons:
Now, isolating in EquationEquation (1)(1) (1) we get an expression for the reactivity:
By replacing EquationEquation (5)(5) (5) into EquationEquation (6)(6) (6) we get:
EquationEquation (7)(7) (7) represents the inverse equation of point kinetics and it is the basis for the construction of digital on-line reactivity meters [Citation13]. This equation is not used directly in the calculation of reactivity because it presents difficulties due to the computational cost in real-time measurements.
2.1. Euler–Maclaurin identity for the calculation of reactivity
In this section, two approximations for the integral in EquationEquation (7)(7) (7) for the calculation of reactivity using the first two Bernoulli numbers are presented, which will allow us to observe a repetitive behavior that helps us obtain a general formula that includes the infinite Bernoulli numbers.
In order to solve the integral that appears in EquationEquation (7)(7) (7) , the Euler–Maclaurin formula [Citation11] is used:
where the are the Bernoulli numbers shown in .
It is possible to simplify the integrand of EquationEquation (7)(7) (7) by doing the following:
EquationEquation (9)(9) (9) represents the continuous response of the system to a unit impulse function [Citation14]. With this response and the neutron population density in the integrand in EquationEquation (7)(7) (7) , it can be obtained:
The discrete version of EquationEquation (10)(10) (10) can be written in the form:
making in EquationEquation (8)(8) (8) and making use of EquationEquation (11)(11) (11) , it is possible to show that EquationEquation (7)(7) (7) can be written as [Citation11]:
where time is discretised as
Being T the time step in the reactivity calculation, , that is n is a non-negative integer. EquationEquation (13)(13) (13) allows to go from a continuous time to a discrete one.
It is possible to find an expression for the second Bernoulli number, by making in EquationEquation (8)(8) (8) . EquationEquations (9)(9) (9) and (Equation11(11) (11) ) are derived three times to obtain:
Making in EquationEquation (14)(14) (14) we get:
Making in EquationEquation (14)(14) (14) we get:
Subtracting EquationEquation (15)(15) (15) from EquationEquation (16)(16) (16) we get,
By replacing EquationEquation (17)(17) (17) into EquationEquation (8)(8) (8) we get:
By replacing EquationEquation (18)(18) (18) into EquationEquation (7)(7) (7) we get:
where
EquationEquation (19)(19) (19) represents the reactivity with two Bernoulli numbers. In the next section, we present the proposed method to find a general formula for the calculation of reactivity that can reproduce particular cases with one or two Bernoulli numbers.
3. Proposed method
A relationship with the Pascal’s triangle can be seen for each one of the terms in EquationEquation (19)(19) (19) . In order to find a method for calculating the reactivity with any Bernoulli number, the binomial theorem can be used to obtain a given order for the derivatives of the system’s response to a unit impulse function. It is possible to express EquationEquations (12)(12) (12) and (Equation19(19) (19) ) in the following way:
The summation of binomial terms in EquationEquation (21)(21) (21) can be expanded in the following form:
EquationEquations (22)(22) (22) and (Equation23(23) (23) ) contain the odd and even terms of the derivatives for the system’s response to the unit impulse and for the neutron population density. By replacing EquationEquations (22)(22) (22) and (Equation23(23) (23) ) into EquationEquation (21)(21) (21) we get:
From EquationEquation (9)(9) (9) it is possible to obtain the derivatives of the system’s response to a unit impulse in the following way:
For the nth derivative of EquationEquation (25)(25) (25) , we have:
By replacing the discrete version of EquationEquation (26)(26) (26) in EquationEquation (24)(24) (24) we get:
The infinite derivatives of the neutron population density that appear in EquationEquation (27)(27) (27) can be reduced by making the following assumption [Citation8]:
where a is given by:
The validity of the approximations given by EquationEquations (28(28) (28) –Equation30(30) (30) ) is met for 7 forms of neutron density, these are: , , , , , and , where and are constants.
By replacing EquationEquations (28)(28) (28) and (Equation29(29) (29) ) into EquationEquation (27)(27) (27) we get:
By taking common factor and in EquationEquation (31)(31) (31) we get:
By the definition of the binomial theorem we have:
It is possible to write the summation in EquationEquation (33)(33) (33) as a sum of odd and even terms:
Equating EquationEquations (33)(33) (33) and (Equation34(34) (34) ) we get:
It is possible to prove the convergence for the second sum in EquationEquation (35)(35) (35) , that is:
By replacing EquationEquation (36)(36) (36) into EquationEquation (35)(35) (35) , it is possible to obtain the convergence for the summation:
By replacing EquationEquations (36)(36) (36) and (Equation37(37) (37) ) into EquationEquation (32)(32) (32) we get the reactivity in the form:
By multiplying term by term in EquationEquation (38)(38) (38) and simplifying, we obtain:
It is possible to simplify EquationEquation (39)(39) (39) by doing:
When replacing EquationEquations (40)(40) (40) and (Equation41(41) (41) ) into EquationEquation (39)(39) (39) it can be verified that the summations over represent a Laurent series of the form [Citation15]:
By replacing EquationEquation (42)(42) (42) into EquationEquation (39)(39) (39) we get an expression for reactivity in the form:
EquationEquation (43)(43) (43) looks complicated, but, it can be written as:
where
EquationEquation (43)(43) (43) represents the reactivity with the approximation of the infinite numbers of Bernoulli and it will be the equation used as the proposed method. In the following section, the results of the most representative numerical experiments are presented. EquationEquation (7)(7) (7) represents the reference method. EquationEquation (12)(12) (12) represents the Euler–Maclaurin method [Citation11]; This is the calculation of reactivity using Bernoulli’s first number and EquationEquation (43)(43) (43) is the proposed method: the Bernoulli number generalisation method (BNGM).
4. Results
In this section, the results of the calculation of the reactivity with the generalisation of the Bernoulli numbers are presented; the computational simulations were made for different forms of neutron population density, as well as for different time steps. The values for the constants considered in this work are the decay constants , the ith fraction of delayed neutrons , the neutron generation time and the total effective fraction of delayed neutrons . The reference method is given by the analytical solution of EquationEquation (7)(7) (7) and is the basis for comparison as to know the accuracy of each method. To avoid division by zero, due to the term ““ in EquationEquation (43)(43) (43) , it is necessary to recourse to the floating-point relative accuracy in the Fortran software: Epsilon (eps). Actually, for numerical experiments using eps, we take eps. We take the eps value in Fortran as equal to . Epsilon (X) returns the smallest number E of the same kind as X such that 1 + E > 1. The difference in reactivity is the error calculated as the absolute value: difference of reactivity = numerical experiment – reference method. The numerical experiment is realized using the proposed method – given by EquationEquation (43)(43) (43) – and the reference method is given by EquationEquation (7)(7) (7) .
The different numerical experiments performed in this work are presented by case studies, cases 1–5 meet the conditions of the derivative method (DM) imposed on the forms of neutron density given by EquationEquations (28(28) (28) –Equation30(30) (30) ). These are linear functions and exponential functions primarily where the value of ω is varied. Cases 6–8 were chosen in such a way that they do not generally meet the conditions imposed on neutron density and were reported in the matrix formulation method (MF) [Citation12].
Case study 1: As the first experiment, we considered the neutron density of the form with ω = 0.00243 with different time steps for the calculation of the reactivity, in a time interval [0,1000] s. In it is possible to observe that the proposed method of the generalisation of Bernoulli numbers (BNGM) has a better approximation than the Euler–Maclaurin method [Citation11]. For large time steps, it can be seen that the method of generalisation of Bernoulli numbers converges to a value which is compared with the Matrix Formulation (MF) method [Citation12].
Case study 2: As the second and third numerical experiments, the same exponential form of neutron density was taken as in the first experiment, but with different values of ω. and show the excellent approximation of the method using the generalisation of the Bernoulli numbers for the different time steps and different time intervals, when compared to the five-point Lagrange method [Citation10]. When compared with the Matrix Formulation method [Citation12], it is evident that the Bernoulli’s generalisation method presents an excellent approximation for time step . It is also observed that the implementation of the eps has no effect on the effectiveness of the method since the errors still have a good approximation.
Case study 3: In this case, the fourth numerical experiment was performed. shows that the method of generalisation of Bernoulli numbers continues to present an excellent approximation with respect to the Euler–Maclaurin [Citation11] for all time steps. With respect to the Matrix Formulation method [Citation12], Bernoulli’s generalisation method surpasses it by two orders of magnitude at and remains with excellent approximation up to .
Case study 4: – present the results obtained from the different numerical experiments from the fifth to the seventh, for the values , and , respectively. It can be observed that the method of generalisation of Bernoulli numbers presents a better approximation for all time steps simulated compared to the matrix formulation method [Citation12]. By comparing the generalisation of the Bernoulli numbers with the Euler–Maclaurin method [Citation11], one can observe once again how a better precision is achieved.
When a time step of T = 0.1 s -being this time step a common value used for the calculation of the reactivity value given in digital meters- is taken, differences on the results were of the order between 10−12 and 10−13 as shown in Tables 3–8 for the proposed method, while the matrix formulation method displays a slightly higher differences which fell between 10−11 and 10−12 magnitude orders.
Case study 5: shows the eighth experiment, the results of the errors for a linear neutron density form. The generalisation method diverges due to the term in EquationEquation (43)(43) (43) since a becomes zero; it is for this reason that the computational arrangement of Fortran (eps) is used; it is evident that even with the eps, the generalisation method is still superior in accuracy to the Matrix Formulation method [Citation12]. Computational simulations were done for the time steps to , since these values are not reported in the literature.
shows a part of the total reactivity curve in pcm to enable the appreciation of the difference between the Matrix Formulation methods [Citation12] and the generalisation of the Bernoulli numbers, where the reference method is obtained by solving EquationEquation (7)(7) (7) for a neutron population density form with a time step . It is possible to see that the matrix formulation method moves away from the solution given by the reference method.
The cases presented so far include numerical experiments that meet the approximations given by EquationEquations (28(28) (28) –Equation30(30) (30) ), where the matrix formulation method (FM) has very good approximations. In general terms, case studies 6–8 did not meet the conditions given by EquationEquations (28(28) (28) –Equation30(30) (30) ).
Case study 6: In the following ninth numerical experiment, the power form of the neutron population density is changed to the form . The results are shown in for a time step for the calculation of the reactivity. It is observed that the method of generalisation of the Bernoulli numbers for this form of neutron density has a difference in the calculation of reactivity lower than the method of the derivatives (DM) [Citation8] and the method of Matrix Formulation (MF) [Citation12]. For the different values of d, the method of generalisation of Bernoulli numbers continues to have a better approximation. It can be observed from this numerical experiment how the proposed method has better approximations when compared to the derivative method (DM), although the two methods used the approximations given by EquationEquations (28(28) (28) –Equation30(30) (30) ). In the proposed method represented by EquationEquation (43)(43) (43) , the dominant term for calculating reactivity for the case study is the first term given by , which has no restriction on the shape of the neutron density.
shows the reactivity curve in pcm for a form of neutron population density where the methods of Matrix Formulation [Citation12], derivatives [Citation8] and the generalisation of Bernoulli numbers are compared; The reference method is obtained by solving EquationEquation (7)(7) (7) and the calculation in the simulation is done up to a final time of with the purpose of better showing the shape of the reactivity curve. The low accuracy of the derivatives method is evidenced with respect to the good precision of the Matrix Formulation method and the excellent precision of the Bernoulli number generalisation method.
Case study 7: The tenth and eleventh numerical experiments were dealt in this case. and show the excellent accuracy of the method of generalisation of Bernoulli numbers for a hyperbolic form of neutron population density. It can be observed that the proposed method overwhelmingly exceeds the Matrix Formulation method, which until now had the best approximations in the calculation of reactivity.
shows a comparison between the proposed method and the Matrix Formulation method for a neutron population density of the form making use of eps. Because of the high precision of the proposed method and the Matrix Formulation method, the axis ranges are reduced to be able to observe the differences, so that a part of the total reactivity curve is shown. These results show the accuracy of the method using the generalisation of the Bernoulli numbers.
Case study 8: The last numerical experiment is presented in this case. shows the difference in the calculation of the reactivity for a neutron population density of the form . It is noted that the method of generalisation of Bernoulli numbers exceeds the Matrix Formulation method for all time steps.
The results presented in the studied cases 6–8, which do not meet the conditions (28–30), are of high precision using the proposed method, reaching several orders of magnitude when compared to the derivative methods (DM) and matrix formulation (MF).
shows the reactivity curve for a neutron population density of the form . Because of the high precision of the proposed method and the Matrix Formulation method with respect to the reference method, the axis ranges are narrowed, so that a part of the total reactivity curve is displayed. It is evident that the method of generalisation of Bernoulli numbers improves upon the reference method given by EquationEquation (7)(7) (7) ; it is observed that the method of Matrix Formulation exceeds the reference method and the method of generalisation of the Bernoulli numbers due to the precision of the methods.
5. Conclusions
In this work, a generalisation for the calculation of reactivity using the infinite numbers of Bernoulli was presented, this expression is simplified, firstly, by an approximation in the derivatives of the neutron density, secondly, using a Laurent’s series. The results presented are highly accurate in the calculation of reactivity when compared to those reported in the literature. It is evident from the different numerical experiments for different forms of neutron population density that the proposed method has a low computational cost since it is highly accurate regardless of the time step. Due to its high accuracy, it is recommended the proposed method to be implemented in a real-time digital reactivity meter.
Acknowledgments
The authors thank the research seed of Computational Physics, the research group in Applied Physics FIASUR, and the academic and financial support of the Universidad Surcolombiana.
Disclosure statement
No potential conflict of interest was reported by the authors.
References
- Weston MS. Nuclear reactor physics. 3rd ed. Weinheim (Germany): Wiley-VCH; 2018.
- Shimazu Y, Nakano Y, Tahara Y, et al. Development of a compact digital reactivity meter and reactor physics data processor. Nucl Technol. 1987;77:247–254.
- Ansari SA. Development of on-line reactivity meter for nuclear reactors. IEEE Trans Nucl Sci. 1991;38:946–952.
- Binney SE, Bakir AIM. Design and development of a personal computer based reactivity meter for a nuclear reactor. Nucl Technol. 1989;85:12–21.
- Hoogenboom JE. Van Der Sluijs AR. Neutron source strength determination for on-line reactivity measurements. Ann Nucl Energy. 1988;15:553–559.
- Tamura S. Signal fluctuation and neutron source in inverse kinetics method for reactivity measurement in the sub-critical domain. J Nucl Sci Technol. 2003;40:153–157.
- Kitano A, Itagaki M, Narita M. Memorial-index-based inverse kinetics method for continuous measurements of reactivity and source strength. J Nucl Sci Technol. 2000;37:53–59.
- Suescún DD, Martinez SA, Da Silva CF. Formulation for the calculation of reactivity without nuclear power history. J Nucl Sci Technol. 2007;44:1149–1155.
- Suescún DD, Martinez SA, Da Silva CF. Calculation of reactivity using a finite impulse response filter. Ann Nucl Energy. 2008;35:472–477.
- Hessam M, Vosoughi N. On-line reactivity calculation using Lagrange method. Ann Nucl Energy. 2013;62:463–467.
- Suescún DD, Rodriguez JA, Figueroa JH. Reactivity calculation using the Euler-Maclaurin formula. Ann Nucl Energy. 2013;53:104–108.
- Suescún DD, Cabrera CE, Lozano PJH. Matrix formulation for the calculation of nuclear reactivity. Ann Nucl Energy. 2018;116:137–142.
- Duderstadt JJ, Hamilton LJ. Nuclear reactor analysis. New York (NY): Wiley; 1976.
- Haykin S, Veen BV. Signal and system. New York (NY): Wiley; 1999.
- Kwok YK. Applied complex variable for scientist and engineers. 2nd ed. New York (NY): Cambridge; 2010.