Nuclear Reactivity Calculation with Reduction of Fluctuations
Corresponding email: daniel.suescun@usco.edu.co
Published at : 28 Jul 2023
Volume : IJtech
Vol 14, No 5 (2023)
DOI : https://doi.org/10.14716/ijtech.v14i5.5970
Suescun-Diaz, D., Ule-Duque, G., Chala-Casanova, J.A., 2023. Nuclear Reactivity Calculation with Reduction of Fluctuations. International Journal of Technology. Volume 14(5), pp. 1113-1122
| Daniel Suescun-Diaz | Department of Exact and Natural Sciences, Surcolombiana University, Neiva, 410001, Colombia |
| Geraldyne Ule-Duque | Department of Exact and Natural Sciences, Surcolombiana University, Neiva, 410001, Colombia |
| Jesus Antonio Chala-Casanova | Department of Exact and Natural Sciences, Surcolombiana University, Neiva, 410001, Colombia |
This article presents a study
based on a series of numerical experiments. It demonstrates the possibility of
reducing fluctuations in the calculation of reactivity using the second
Bernoulli number based on the approximation of the Euler-Maclaurin formula.
This approach requires knowledge of the first three derivatives, which are
implemented progressively. The fluctuations are assumed to occur around an
average value of the neutron density with a Gaussian distribution. Jitter
reduction is performed with a first-order delayed low-pass filter for different
forms of neutron population density, with different time steps and with
different filter constants. The numerical results show that the method can be
used as a digital reactivity meter.
Inverse point kinetics equation; Neutron population density; Numerical experiment; reactivity; Second Bernoulli number
The increasing world population has led to a higher
demand for electrical energy in recent years. As a critical component of
supporting a country's economic growth (Saroji
et al., 2022), the electric power
system plays a crucial role in meeting this demand. An alternative path would
be to provide a solution where home automation systems can reduce electricity
consumption (Rabbani and Foo, 2022). Another viable option is urban wind energy which is one of the new renewable
ways of producing electricity; however, researchers have not studied very much
in this field (Krasniqi, Dimitrieska, and Lajqi, 2022). Nuclear energy is another
viable option obtained through nuclear reactors. However, it is necessary to
know the reactivity parameters with good accuracy in nuclear reactors to
operate a nuclear power plant more safely. Therefore, one of the main functions
of a nuclear power plant is to control and monitor reactivity (Hyvarinen et al., 2022).
In recent decades, different experimental and
computational methods have been developed to estimate the reactivity value in
the core of a nuclear reactor. Studies have been carried out in a BAEC TRIGA
Mark- II research reactor to analyze the effects of reactivity insertion, as
well as in a prototype fast breeder reactor Monju (Hossain
et al., 2022; Ohgama et al., 2022). The results are validated by a deterministic model given by the Inhour
equation and the Monte Carlo method. The in-hour equation is also used
in estimating reactivity in experiments conducted in the light water reactor at
the VENUS-II experimental facility (Jiang et al., 2022).
Solving the inverse point kinetics equation provides a mathematical model that allows the calculation of reactivity. This method is commonly employed in computer-based simulations and facilitates the development of digital reactivity meters. To apply this model, the neutron density is required as an input, which can be measured using devices like ionization chambers (Vasilenko, Pankin, and Skvortsov, 2019). Studies that take this perspective are being conducted in the context of the decommissioning of the Fukushima Daiichi nuclear power plant. This includes a preliminary analysis that identifies the range that applies to the MIK code, such as ramped reactivity insertion (Fukuda, Nishiyama, and Obara, 2021).
A method has been developed to correct reactivity
values by accounting for changes in both the neutron flux function and detector
efficiency (Zhitarev et al., 2021). Based on an analysis of the inverse point kinetics equation, the
influence of the background current on the measured reactivity is analyzed and
a method for iterative calculation of reactivity is introduced (Huo et al., 2019). An
accurate numerical solution for the inverse point kinetics equation is given
using the discrete Fourier transform (Suescun-Diaz, Lozano-Parada,
and Rasero-Causil, 2019). The extended Kalman filter technique (Bhatt et
al., 2013) and the wavelet-based
multiscale extended Kalman filter technique have also been proposed for
reactivity estimation (Patel Mukhopadhyay,
and Tiwari, 2018).
However, reactivity meters based on the inverse point kinetics equation have
sufficient capabilities to accurately estimate reactivity (Shimazu, 2014).
Due to gamma radiation, electronic system noise, and
environmental radiation, there is considerable noise in the electrical current
during reactivity measurement by external detectors, which leads to significant
reactivity error, especially at low powers (Huo et
al., 2019). Under these noise conditions,
it is necessary to apply a filtering tool to smooth or reduce the fluctuations
in the measurements.
In the present work, we study the approximation to
solve the integral in the inverse point kinetic equation using the Euler-Maclaurin formula (Arfken, Weber, and Harris, 2013), which
provides a powerful connection between integrals and sums, considering the approximation of the second
Bernoulli number, with the combination using the first-order delayed low-pass
filter to reduce fluctuations in the reactivity calculation. The results
indicate that it is an alternative method for reactivity calculation with good
approximation and can be used to design a digital reactivity meter.
2.1. Theoretical
Considerations
The point kinetics equations are a set of differential equations that
describe the time evolution of the expected values of the neutron density and
the concentration of delayed neutron precursor groups in the core of a nuclear
reactor. These equations accurately describe the reactor core dynamics and
correspond to the time component of the neutron diffusion equation under the
assumption of an isotropic and homogeneous medium (Stacey,
2018). They are described as follows:
where, P(t) is the neutron density, Ci(t) is
the concentration of the i-th group of delayed neutron precursor, 
is the reactivity, 
is the neutron generation time,
is the i-th fraction of delayed neutrons, 
is the total effective
fraction of delayed neutrons, 
is the decay constant of the
i-th group of delayed neutron precursors.
Solving for from equation (1) leads to the following reactivity equation:
The unknown term in equation (5) is the concentration of precursors Ci(t), which can be found by solving equation (2) by an integrating factor or by
the Laplace transform -considering equations (3) and (4)- will lead to the
following expression:
By replacing equation (6) with equation (5), a new equation for reactivity
is obtained, which needs all the values of the neutron density to be known:
Equation (7) is the so-called inverse point kinetic
equation and allows the calculation of the nuclear reactivity, which provides
information on the behavior of a reactor core. This equation has been a model
that has been applied in the design of digital reactivity meters that
contribute significantly to the control and safe operation of a nuclear
reactor. However, its dependence on all the values of the neutron density in a
non-recursive way causes a high computational cost, which makes it difficult to
calculate the reactivity in real-time. To reduce the computational cost, the
following section presents a method that discretizes the integral containing
the dependence on the neutron density by using the Euler-Maclaurin formula with
an approximation of two Bernoulli numbers.
2.2. Proposed Method
To achieve greater accuracy in
reactivity calculations while minimizing computational costs, it is necessary
to discretize the integral term in equation (7). This is accomplished by
applying the Euler-Maclaurin formula as follows (Kwok,
2010):
where the term Bk represents the
Bernoulli numbers.
Substituting equation (8) into
equation (7), reactivity with the approximation of the first Bernoulli number B1=1/6
is obtained (Suescun-Diaz, Ule-Duque,
and Pena-Lara, 2020). To find a descriptive expression for reactivity with the approximation
of the second Bernoulli number B2=1/30, substitute k=2
into equation (8), taking derivative three times leads to:
Being 
the time step in
the reactivity computation, n indicates the discrete-time, and its
relation to the continuous time is 
is the
system response to a unit impulse function (Haykin,
2014) defined here as 
Replacing equation (9) into equation (7), the
following expression for reactivity is obtained:
Equation (10) represents the reactivity with the
approximation of the first two Bernoulli numbers, being 
the correction of the first Bernoulli number represented in equation (12), and 
the correction with the second Bernoulli number represented in the equation
(13).
For the calculation of the first, second and third
derivatives, the progressive derivatives are implemented (Mathews and Fink, 2005) as given
by equations (14-16):
where
is the average
value of an arbitrary deterministic signal, 
is a Gaussian
noise, 
is the standard
deviation and 
will be the new
signal containing fluctuations. To reduce the fluctuations in the calculation
of reactivity, a low-pass filter with first-order delay is considered here (Tamura, 2003; Kitano et al., 2000), which is given by the expression:
Where 
is the filtering constant.
There are 422 nuclear reactors in operation and
another 56 under construction with 377 872 MWe and 58 584 MWe total net
installed capacity, respectively. The reactivity value is critical for ensuring the
safe operation of nuclear reactors. Therefore, the objective of the proposed
method in this work is to numerically solve the inverse point kinetics equation
using equation (10). The low-pass filter given by equation (18) is proposed to
reduce the fluctuations of an input signal associated with neutron population
density measurements.
The simulations were implemented using the MATLAB
computational tool. The physical constants used in this work are due to the
interaction of neutrons with the combustible material 235U.
These constants are the decay constant = [0.0127, 0.0317,
0.115, 0.311, 1.4, 3.87]s-1, the delayed neutron fraction = [0.000266, 0.001491, 0.001316, 0.002849, 0.000896, 0.000182], the total
delayed neutron fraction 
and the instantaneous neutron
generation time =2×105s (Ganapol,
2013). For noise simulation, a seed
generating random numbers of 231-1 is used. Some results of
the different numerical experiments for calculating reactivity are presented in
the next section. Consider different forms of neutron density, time steps and
filter constants.
The
physical parameters in this results section are considered as above for a
thermal reactor with 235U fuel elements. The most outstanding
results obtained for different numerical experiments are presented using the
proposed method given by equation (10) and denoted by 
with different values of 
obtained from the inhour
equation (Kinard and Allen, 2004). Finally, the neutron
population of form P(t)=a+bt3
is considered with different
values of b. The reference method is obtained by the analytical solution
of equation (7), which allows an analysis of the accuracy of the method.
It is necessary to know
the reactivity with high precision; however, in practice, the neutron
population density contains noise, which has a Gaussian distribution that
produces fluctuations or uncertainties in the reactivity calculation. To reduce
fluctuations, the first-order delay low-pass filter given by equation (18) with
a filtering constant of 
and 
is used. In all numerical experiments, the time step
varies in the interval [0.001, 0.1] s, and the standard
deviation is 
For the different derivatives of neutron
density represented in equations (11-13), the progressive derivatives are taken
(Mathews and Fink, 2005), represented by equations (14-16).
Table 1 shows the maximum
differences in reactivity in pcm (parts per hundred thousand) for a
neutron density of the form 
with a value of 
obtained from the inhour
equation that provides a reactivity of about 50 pcm. It is
possible to observe that for this reactivity value, the reduction of the
fluctuations (RF) is effective for any time step, obtaining a reduction above 50%
for any case, although the most significant occurs when a constant filter 
and a time step 
are used, reaching a reduction of
84.06%. This significant reduction indicates that the uncertainty in the
reactivity value decreases. In other words, it increases the precision in the
calculation of the reactivity that would achieve greater control of the
reactor. This RF is calculated as follow:
This RF is calculated as shown in equation (19).
To
validate the results obtained with the proposed method for the exponential form
of neutron density, a time-step of 
is considered for the EM2F case.
The results show a maximum difference of 0.88 pcm. In a study by Suescun-Díaz, Lozano-Parada, and Rasero-Causil (2019) under the
same conditions, a maximum difference of 2.06 pcm was reported using the
discrete Fourier transform. These results clearly demonstrate that the EM2F
method yields the highest level of reduction.
Table
1 Maximum difference in reactivity
In order to be able to make decisions in the operation
of a nuclear reactor using control rods, the reduction of fluctuations is
presented in Figures 1 and 2, the reactivity for a form of neutron density with 
and without low-pass filter (EM2) and with a low-pass filter (EM2F) at a filter
constant of ,
respectively. It can be observed that when a first-order low-pass filter is
applied, the fluctuations are effectively reduced.
Figure 1 Reactivity
for a neutron density , without low-pass filter
Figure 2 Reactivity
for a neutron density with with low-pass filter
Table 2 considers the neutron density of the form , , which gives a reactivity of about 20 pcm.
It is evident that for time steps 
the reduction of
fluctuations with the proposed method is significant. The achieved reduction
ranges from 56.64 % to 92.17 %. In this numerical experiment, it
was found that the fluctuations are reduced even for a time step , reaching a minimum reduction of 56.64 %.
Table 2 Maximum difference in reactivity
Table 3 shows the maximum differences in reactivity
caused by the neutron density of the form with a value of 
where the reactivity is about 70 pcm.
In this case, it can be observed that for a constant filter ,
the reduction of fluctuations is very good for time steps of .
Increasing the value of the filter constant up to , we find
that the best reduction is obtained for , with the
efficiency decreasing with increasing time step, reaching a reduction level of 11.15%
for 
the results found
show that it is not possible to reduce the fluctuations.
Table 3 Maximum difference in reactivity
In the following numerical experiment, a neutron
density of the form 
is taken. For this value of 
a reactivity value of 140 pcm is produced. In
Table 4 can be observed that the efficiency of the method decreases when 
either with 
It is noted
that the greatest reduction occurs at 
with a reduction percentage of 77.98 %, compared to the case in which
the filter constant is , which represents a reduction of
31.41 % for the same time step.
Table 4 Maximum difference in reactivity
In another type of numerical experiment, Table 5
presents the neutron density for the form 
with 
It is possible to observe that
for the different values of b, the value of the maximum difference of
the proposed method EM2 with noise remains constant. When applying the
low pass filter with a filtering constant of , the value of
the maximum difference is constant for the first two values of b, with a
reduction of 78.37% in both cases. When the filtering constant is
increased to , it is evident that the maximum differences
remain constant at 0.49 pcm for any b value, thus producing constant
reductions of 80% in the fluctuations presented in the calculation of
reactivity. This constant value in the maximum
difference is due to the attenuation caused by using the low-pass filter, as
cited in (Suescun-Díaz,
Lozano-Parada,
and Rasero-Causil, 2019).
Table 5 Maximum difference
in reactivity
Figures 3-4 show the reactivity for a neutron density
of the form , with a = 1, b = (0.0127)5/9
and a filtering constant of . It is possible to observe the
effective reduction of the fluctuations that agree with the data shown in Table
5.
Figure 3 Reactivity for a neutron density of the form , without low-pass
filter
Figure 4 Reactivity
for a neutron density of the form with low-pass filter
The Euler-Maclaurin method is presented here using the
approximation of the second Bernoulli number to solve the integral of the
inverse point kinetics equation that depends on the neutron density. This
approximation uses the first-order delay low-pass filter to reduce the
fluctuations in the reactivity calculation. The results of the different
numerical experiments show that the proposed method can be considered for the
exponential and cubic forms of the neutron density. The method achieves a
better reduction of fluctuations in numerical experiments considering the shape
of the exponential neutron density when the time-step is of the order of 
achieving reductions of more than 70%. When the neutron
density was changed to a cubic form, it was observed that the reductions were
almost constant and reached 80%. The limitations of the proposed method
for small values in the standard deviation of 
should be improved by
using another method for filtering and is a future work that can be studied
using digital signal processing. The results obtained indicate that the
proposed method can be an alternative to be implemented in a digital reactivity
meter when there is noise in the neutron population density.
Arfken, G.B., Weber,
H.J., Harris, F.E, 2013. Mathematical Methods for Physicists A Comprehensive
Guide. 7th Edition. UK: Elsevier
Bhatt, T.U.,
Shimjith, S.R., Tiwari, A.P., Singh, K.P., Singh, S.K., Singh, K., Patil, R.K.,
2013. Estimation of Sub-Criticality Using Extended Kalman Filtering Technique. Annals
of Nuclear Energy, Volume 60, pp. 98–105
Fukuda, K.,
Nishiyama, J., Obara, T., 2021. Supercritical
Transient Analysis For Ramp Reactivity Insertion Using Multiregion Integral
Kinetics Code. Nuclear Science and Engineering, Volume 195(5), pp.
453–463
Ganapol, B.D., 2013.
A Highly Accurate Algorithm for The Solution of The Point Kinetics Equations. Annals
of Nuclear Energy, Volume 62, pp. 564–571
Haykin, S.S., 2014. Adaptive
Filter Theory. 5th Edition. London, UK: Pearson
Hossain, N., Soner,
M.A., Prodhan, M.M., Sahadath, M.H., Kabir, K.A., 2022. Experimental Analysis
of STEP Reactivity Insertion Effect on Reactor Power, Fuel Temperature and
Reactor Period in BAEC TRIGA Research Reactor. Annals of Nuclear Energy,
Volume 165, p. 108665
Huo, X., Fan, Z., Xu,
L., Chen, X., Hu, Y., Yu, H., 2019. A New and Efficient Method to Measure
Reactivity in a Nuclear Reactor. Annals of Nuclear Energy, Volume 133,
pp. 455–457
Hyvarinen, J.,
Vihavainen, J., Ylonen, M., Valkonen, J., 2022. An Overall Safety Concept for
Nuclear Power Plants. Annals of Nuclear Energy, Volume 178, p. 109353
Jiang, W., Zhu,
Q.-F., Zhou, Q., Ma, F., Zhang, L., Liu, Y., Li, J.-Y., Ge, H.-L., Yu, R.,
Meng, H.-Y., Wang, D.-W., Chen, L., 2022. Reactivity Worth Measurement of The
Lead Target On Venus-II Light Water Reactor and Validation of Evaluated Nuclear
Data. Annals of Nuclear Energy, Volume 165, p. 108779
Kinard, M., Allen,
E.J., 2004. Efficient Numerical Solution of The Point Kinetics Equations In
Nuclear Reactor Dynamics. Annals of Nuclear Energy, Volume 31(9), pp.
1039–1051
Kitano, A., Itagaki,
M., Narita, M., 2000. Memorial-Index-Based Inverse Kinetics Method for
Continuous Measurement of Reactivity and Source Strength. Journal of Nuclear
Science and Technology, Volume 37(1), pp. 53–59
Krasniqi, G.,
Dimitrieska, C., Lajqi, S., 2022. Wind Energy Potential in Urban Area: Case
study Prishtina. International Journal of Technology, Volume 13(3), pp.
458–472
Kwok, Y., 2010. Applied
Complex Variables for Scientists and Engineers. 2nd Edition. UK:
Cambridge University Press
Mathews, J.H., Fink, K.D., 2005. Numerical Methods Using Matlab. 4th
Edition. New Delhi, India: Pearson
Ohgama, K.,
Takegoshi, A., Katagiri, H., Hazama, T., 2022. Evaluation Of Fuel Reactivity
Worth Measurement in The Prototype Fast Reactor Monju. Nuclear Technology,
Volume 208(10), pp. 1619–1633
Patel, S.B.,
Mukhopadhyay, S., Tiwari, A.P., 2018. Estimation Of Reactivity and Delayed
Neutron Precursors’ Concentrations Using a Multiscale Extended Kalman Filter. Annals
of Nuclear Energy, Volume 111, pp. 666–675
Rabbani, N.A., Foo,
Y.-L., 2022. Home Automation to Reduce Energy Consumption. International
Journal of Technology, Volume 13(6), pp. 1251–1260
Saroji, G., Berawi,
M.A., Sari, M., Madyaningarum, N., Socaningrum, J.F., Susantono, B., Woodhead,
R., 2022. Optimizing the Development of Power Generation to Increase the
Utilization of Renewable Energy Sources. International Journal of Technology,
Volume 13(7), pp. 1422–1431
Shimazu, Y., 2014. A
Simple Procedure to Estimate Reactivity With Good Noise Filtering
Characteristics. Annals of Nuclear Energy, Volume 73, pp. 392–397
Stacey, W.M., 2018. Nuclear
Reactor Physics. 3rd Edition. Atlanta, USA: Wiley-VCH
Suescun-Diaz, D.,
Lozano-Parada, J.H., Rasero-Causil, D.A., 2019. Novel Fluctuation Reduction Procedure for Nuclear Reactivity
Calculations Based on The Discrete Fourier Transform Method. Journal of Nuclear Science and Technology, Volume 56, pp. 608–616
Suescun-Diaz, D., Ule-Duque, G., Pena-Lara, D., 2020. Reduction of Reactivity
Fluctuation with the Euler-Maclaurin Method. ARPN Journal of Engineering and
Applied Sciences, Volume 15, pp. 96–103
Tamura, S., 2003. Signal Fluctuation and Neutron
Source in Inverse Kinetics Method For Reactivity Measurement in The
Sub-Critical Domain. Journal of Nuclear Science and Technology, Volume 40, pp. 153–157
Vasilenko, V.A., Pankin, A.M., Skvortsov, K. V., 2019.
Reactivity Calibrator. Atomic Energy, Volume 125(3), pp. 157–161
Zhitarev, V.E., Kachanov, V.M., Sergevnin, A.Y.,
Sutemieva, N.E., 2021. Correction of the Reactivity Measured by Oruk Technique.
Physics of Atomic Nuclei, Volume 84(8), pp. 1405–1412