Neutron Density Calculation Using the Generalised Adams-Bashforth-Moulton Method

A nuclear reactor is a device to initiate, control and maintain nuclear fission reactions. The understanding of nuclear fission processes requires the study of the neutron population and the concentration of delayed neutrons, which can be modelled by the equations of point kinetics. The reactor point kinetic equations are the reduced-order model most representative in the nuclear reactor dynamics. Abstract

The calculation of reactivity depends on the neutron density which makes the point kinetics equations nonlinear, and therefore very difficult to solve analytically. The fast and delayed neutron lifetimes are of different orders of magnitude, converting the equations of point kinetics into a stiff system. Several methods have been proposed to solve the equations of point kinetics: In [1], the confinement method (SCM) is proposed to overcome stiffness, whereas other authors have used the generalised Runge-Kutta (GRK) method [2], the Padé approximations [3], the generalization of the analytical inversion method (AIM) [4], and the analytical method (AEM), which used exponential functions [5]. It has also been demonstrated that the point kinetics equations can be solved numerically using the reactivity Piecewise Constant Approximation method (PCA) [6], using a numerical algorithm called: Constant Reactivity (CORE) to calculate nuclear density [7,8] presented a numerical integral method and investigated the neutron density produced by inserting different forms of reactivity in thermal reactors with multiple groups of the neutrons using the Best Function (BBF). The Taylor series (TSM) was used to calculate the nuclear density with feedback reactivity [9], and the power series method (PWS) has been used to obtain approximate solutions with and without feedback [10][11]. [12] introduced a highly accurate algorithm, combining the Backward Euler Finite Difference method (BEFD) and [13] described a semi-analytical method to solve the equations of point kinetics with a technique called (EPCA) which iteratively corrects the error in the source term with good accuracy. The reduced form of the differential transform method (reduced DT Method) [14] was used to obtain the total neutron density. In [15] the Taylor-Lie series was applied to numerically solve the equations of point kinetics with a quadrature technique. [16] introduced the Haar Wavelet Operational Method (HWOM). In [17] an explicit analytical solution was developed from a Taylor series expansion (ITS2). [18] introduced a new method based on the Trigonometric Fourier Series (TFS) to obtain approximate solutions. In [19] the Exponential Time Differencing Method (ETDM) with a Taylor series approximation was proposed to solve the equations of point kinetics; this is a semi-analytic and self-starting method in which the point kinetic equations are integrated using an integration factor. [20] introduced a matrix form known as Treatment Theta Method (TTM) and in a more recent work the Adams-Bashforth and Adams-Moulton Methods (ABM) were used in nuclear reactivity calculations [21]. Finally, in a recent work [22], neutron density was calculated using the Magnus expansion. This new work proposes to solve the equations of point kinetics for the calculation of nuclear density with the generalised Adams-Bashforth and Adams-Moulton predictor-corrector method of order O(h 4 ), with the objective of studying various predictors using the same corrector. This type of multipass method follows a prediction and correction scheme.
The present work is structured as follows: firstly, a section considering a physical-mathematical model, which contains the theoretical framework and the description of the variables or parameters of the model; secondly, a section with the proposed method, which explains the methodology that was carried out in this research. Finally, a section with the results and conclusions.

Model formulation
The point kinetics equations describe the time behaviour of the nuclear reactor. These equations are obtained from the equation of neutron diffusion with six groups of delayed neutrons, as given by [23].
with the following initial conditions: Where P (t ) is neutron density C i (t ) is the concentration of the i-th group of delayed neutron precursors, ρ(t ) is the reactivity, Λ is the neutron generation e i time, β i is the fraction of th -th delayed neutron group β is the total fraction of delayed neutrons, (β = i β i ), λ i is the decay constant of the i-th group of delayed neutron precursors.

Adams-Bashforth and Adams-Moulton Methods
The Adams-Bashforth-Moulton multistep method (ABM) solves the initial value problem of the form y = f (t k , y k ) using information from points previous to y k to find y k+1 , generally, these points are obtained with the Runge-Kutta method of order O(h 4 ). The method of indeterminate coefficients is used to calculate the solution of the function f (t k , y k ), the differential equation with indeterminate coefficients for the generalised fourth-order predictor of the Adams-Bashforth method can be represented by the following equation: The coefficients α i , β i for i = 1, ..., 4 are calculated by carrying out a Taylor series expansion around h on each term y k−1 ,y k−2 ,y k−3 ,y k−1 ,y k−2 y k−3 , of Eq. (5).
We have that Replacing Eqs. (6, 7) into Eq. (5) and equating to the Taylor series for Y k+1 , that is, making i = −1 in Eq. (6), we obtain a set of algebraic equations By solving Eqs. (8) we obtain the coefficients/predictors. Table 1 shows these values, which are then replaced into Eq. (5), thus, different predictor formulas can be obtained. For the error estimation, it must be taken into account that Similarly, the corrector of the Adams-Moulton method can be obtained or deduced from the fundamental theorem of calculus, and the Lagrange's interpolator polynomial is given by: The truncation error of Eq. (9) is: The method's convergence criterion is obtained by making successive adjustments to the corrector in the following way: where, y p , y c are values in the predictor-corrector formulas, y c c is the first correction in the corrector After some algebraic manipulation, Eq. (11) can be written as: After a new correction, the result will be: Thus, by adding all the corrections, a geometric series with ratio r arises In order for the geometric series as given by Eq. (14) to converge, we must have r < 1. Rearranging for h gives the convergence criteria as The case predictor proposed (P3), is selected from Table 1, obtaining: The predictor error of Eq. (16) is given by: In order to obtain better results, we can find the corrector's modifier and the predictor in the following way Replacing Eq. (10) and Eq. (17) into Eq. (18), we get: The term that contains the derivative can be expressed as Multiplying Eq. (20) by 19 720 , and making use of Eq. (10) we get By replacing Eq. (21) into Eq. (18) we get Adding − p k+1 to both sides of Eq. (22) we get Eq. (23) can be used to improve the value of the predictor and to estimate the error both in the predictor and the corrector, without having to calculate the fifth derivative. A slow change in the difference between the predictor and corrector from one step to the next is assumed. Then P k+1 is replaced by P k , y k+1 by y k in Eq. (23) to obtain the following modified formula for the predictor: Similarly, the modified formula for the corrector is:

Results and discussion
The Adams-Bashforth-Moulton method is applied to solve the equations of point kinetics with one or six groups of delayed neutrons. Three cases of reactivity functions are presented: step, ramp and sinusoidal. Each case will be studied separately and the results will be compared to other methods available in the literature. Step Reactivity For the neutron density calculation the following parameters for a thermal reactor are considered: Table 2 and Table 3 show different values for neutron density with reactivities of 300 pcm = 0.43 ($) and 700 pcm (parts per hundred thousand) = 1 ($) reported in the literature [17], these are values for critical and subcritical reactors. The dollar symbol 1 ($) is equivalent to 1 β = 0.007 . The step size in the proposed ABM method is h = 10 −3 s , unfortunately, the authors do not show the time step used in their reported methods. Fig. 1 shows how stiffness forces an abrupt change in the neutron density for a reactivity of 300 pcm = 0.43 ($). This is due to the big differences between Λ = 5 × 10 −4 s and 1/λ i .    The ABM method can be applied to a thermal reactor with the values obtained in Yun Cai et al.  Table 4 shows the results obtained with the ABM method and those reported in the literature for a reactivity ρ = 1 ($).

Ramp reactivity
In this section, the behaviour of the neutron density is determined for ramp reactivities of two types: moderately fast and very fast, which correspond to thermal and fast reactors, respectively.
In the case of moderately fast ramp reactivities, the parameters for the thermal reactor are obtained from  Table 5 and Table 6 show the calculation of neutron density given a reactivity of the form ρ(t ) = 0.0007*t, in this case ρ = 0.1($)/s. The results using the proposed ABM method are compared with the values reported in references [17] and [22]. The proposed method uses the step size h = 10 −3 s again. It is observed that the ABM method is very accurate when compared to other methods in the literature. For this type of reactivity, the neutron density increases exponentially in such a way that for a time t = 9 s there are approximately 487 neutrons / cm 3 , as can be observed in Fig. 2  In the case of very fast ramp reactivities, the ABM method can be applied to a fast reactor with the parameters obtained in Picca [13] using the EPCA method, and from Yun Cai et al.    Table 7 shows the results obtained with the ABM method and those reported in the literature for a reactivity ρ = 1 ($)/s.

Sinusoidal Reactivity
In order to validate the ABM approach, the results obtained for the sinusoidal reactivity are compared with the results obtained in [13] and [22] using the ME3, PCA / ME2 and EPCA methods with nonlinear reactivity insertion for a fast reactor with the following parameters: where λ 1 = 0.077 s −1 , β = β 1 = 0.0079, Λ = 1 × 10 −6 s . The results are shown in Table 8.
Another validation for the approximation of the ABM method consists in comparing it with the results obtained in [19] using the ETD and TSM methods with nonlinear reactivity insertion for a thermal reactor with the following parameters: λ 1 = 0.0127 s −1 , λ 2 = 0.0317 s −1 , λ 3 Table 9. This form of inserted reactivity produces an oscillation in the neutron population density, as can be seen in Fig. 3.
The results presented using the proposed method of ABM of order 4 can be summarised in the following way: using a time step h = 0.001s, for thermal reactors with Step Reactivity, at least nine significant figures could be obtained.  Table 6. Neutron density obtained with a reactivity of ρ = 0.43 ($).   Table 7. Neutron density for reactivity ρ = 1 ($).  Table 9. Neutron density for a thermal reactor with sinusoidal reactivity.

Conclusions
In this paper, the generalised predictor of the Adams-Bashforth-Moulton method (ABM) was presented. This is a method based upon prediction and correction, that guarantees a higher precision when numerically solving a system of differential equations, in this case applied to the point kinetics equations. The results were presented for case P3 of