A New Accurate Numerical Method Based on Shifted Chebyshev Series for Nuclear Reactor Dynamical Systems

A new method based on shifted Chebyshev series of the first kind is introduced to solve stiff linear/nonlinear systems of the point kinetics equations. The total time interval is divided into equal step sizes to provide approximate solutions. The approximate solutions require determination of the series coefficients at each step. These coefficients can be determined by equating the high derivatives of the Chebyshev series with those obtained by the given system. A new recurrence relation is introduced to determine the series coefficients. A special transformation is applied on the independent variable to map the classical range of the Chebyshev series from [−1, 1] to [0, h]. The method deals with the Chebyshev series as a finite difference method not as a spectral method. Stability of themethod is discussed and it has proved that themethod has an exponential rate of convergence.Themethod is applied to solve different problems of the point kinetics equations including step, ramp, and sinusoidal reactivities. Also, when the reactivity is dependent on the neutron density and step insertion with Newtonian temperature feedback reactivity and thermal hydraulics feedback are tested. Comparisons with the analytical and numerical methods confirm the validity and accuracy of the method.


Introduction
For reactor safety, it is important to indicate the neutron density by solving the point kinetics equations PKEs. These equations present an expedient tool to provide useful information on the dynamics of the reactor operation during startup or shut-down. The equations of the point kinetics can be derived from the classical transport problem by eliminating the spatial effects. The effect of elimination the spatial effects can be compensated by improving the solution of the PKEs. The PKEs with Newtonian temperature feedback constitute a coupled stiff system of nonlinear differential equations of the neutron density and the delayed neutron precursor concentrations. Difficulties appearing with solving such stiff systems arise from the disparity between the appearance of prompt and delayed neutrons. There are many literatures that presented several numerical studies to solve such equations. Recently, for example, Cai [1] used Magnus expansion to solve the nonlinear PKEs. Sérgio used the ITS2 method to solve the PKEs with adaptive time step [2] and used the same method to solve the PKEs in the integral formulation with arbitrary number of delayed neutron groups and Newtonian temperature feedback [3]. Hamada [4,5] introduced a new method based on Fourier series expansion and then used adaptively step size to solve stiff systems of the PKEs. Nahla [6] used analytical exponential model to solve the stochastic PKEs via eigenvalues and eigenvectors. The modified exponential time differencing method [7] is applied to solve the PKEs. Picca [8] provided a method based on piecewise constant approximation to solve the nonlinear PKEs. Wollmann da Silva [9] eliminated the stiffness character by splitting the corresponding matrix into a diagonal matrix plus a matrix that contains the remaining terms. Hamada [10,11] developed a solution of the PKEs based on the power series method with constant and with varying step integrations. Finally, there are two of the most accurate methods that currently exist, the converged accelerated Taylor series (CATS) [12] and the backward Euler difference scheme (BEFD) [13] methods. We compared our results with both methods that produced results with high accuracy.
In this paper, a new method based on shifted Chebyshev series of first kind (SCS) is introduced to solve stiff linear/nonlinear systems of the PKEs. The proposed method is designed to inherit the best properties of both spectral 2 Science and Technology of Nuclear Installations and finite difference methods. We start by dividing the total interval in the independent variable , over which a solution is desired [ 0 , ], into equal subintervals. The SCS method has employed finite Chebyshev series expansion to approximate the exact solutions at + 1 equidistance values of . We have transformed the traditional interval of the independent variable of the Chebyshev series from [−1, 1] to [0, ℎ] to be suitable for the iteration process. Hence, with large values of the truncated coefficients and suitable choice of ℎ, more accurate results can be obtained.
Determination of the Chebyshev coefficients can be obtained by equating both the solution derivatives of the Chebyshev series and those of the given system. This procedure produces an algebraic system used to generate a general formula for a recurrence relation. The new recurrence relation is very simple to programme and used to determine the successive coefficients of the Chebyshev series. To study the stability of the proposed method and also to determine the actual error, the method is applied to the model equation. Such application is usually used to test the performance of methods. We have obtained satisfactory results because of the excellent exponential convergence rate of the Chebyshev approximations.
To test the SCS method and to prove its validity and accuracy for application purposes, results of the SCS method are compared with those of analytical and different numerical methods used to solve different systems of the PKEs. Six different systems are studied; the first is a preliminary study for a stiff linear system, with constant reactivity, where the analytical solution is known. This case study aims to test the efficiency of the new technique in a simple case for further applications when the analytical solutions are not available. The remainder five test cases are stiff nonlinear systems. Numerical evaluation performed by the SCS method has been coded by Matlab for personal computer (Intel(R) Core (TM) i7-2630QM CPU @ 2.00 GHz).

The Point Kinetics Equations with Temperature Feedback
The system of the point kinetics equations for neutron density ( ) and −delayed neutron precursor groups ( ), = 1, 2, . . . , with temperature feedback can be described by the differential equations [14]: where ( ) is average neutron flux density, ( / 3 ), ( ) is ℎ groups of delayed neutron precursors, ( = 1, 2, . . . , ), ( is the ℎ group decay constant ( −1 ), is the ℎ group delayed fraction, = ∑ =1 , Λ is the neutron generation time ( ), is the reciprocal of the reactor heat capacity, and is the total number of delayed neutron groups. The functions ( ) and ( ) may be constants or functions of the independent variable. In general, the reactivity ( ) = ( , , ) is prescribed as a constant, linear, or nonlinear function of neutron density, temperature, or both. The previous system can be rewritten in a matrix form as Science and Technology of Nuclear Installations   3 Here, ( ), 0 ( ) are vectors of ( +2)-dimensional real space, ( ) is the square point kinetics matrix of order ( + 2) × ( + 2), and ( ) is a column vector which contains all the dependent variables.

Shifted Chebyshev Polynomials on [0, ℎ]
Assume that the standard finite Chebyshev series is defined by [15] where * ( ) is used here to express the usual Chebyshev polynomials of the first kind; the independent variable is defined on [−1, 1]; and is a parameter which refers to the number of the truncated coefficients. Since the interval in the independent variable is divided into subintervals or steps as we have mentioned in the introduction section, then for the iterative process it is quite more convenient to use the smaller range [0, ℎ] than the larger range[−1, 1]. Therefore, we may map the independent variable in [−1, 1] to the variable in [0, ℎ] by the transformation Such transformation leads to the shifted Chebyshev polynomials of the first kind ( ) of the form This transformation leads to a new type of the series of polynomials From (10), it can be concluded that The successive Chebyshev polynomials can be defined by the recurrence relation which together with 0 ( ) = 1 and 1 ( ) = 2 /ℎ−1 recursively generates all the polynomials { ( )} very efficiently. The high derivatives of (12) can be deduced by the following recurrence relations: For the approximation of ( ), we should expect to get a smaller mini-max error with ( ) in 0 ≤ ≤ ℎ than with * ( ) in the larger range −1 ≤ ≤ 1 according to the following theorem.
Theorem 1 (Lagrange-Chebyshev approximation [16]). Assume that ( ) is the Lagrange polynomial of order that is based on the Chebyshev interpolating nodes on [ , ].
Theorem 1 shows that the error depends basically on the length of the step size ℎ = − and on the number of the truncated coefficients . It should be noted that the interval [0, ℎ] minimizes the error ( ) better than the interval [−1, 1].

Derivation of the SCS Method
Let the approximate solution ( ) of system (5) be expressed by the finite shifted Chebyshev series where ( ) is defined in (10) as a function of the independent variable and the initial value of ( ) can be obtained by where the coefficients , = 1, 2, . . . , , are vectors of order ( + 2). The subscript refers to the number of the coefficients truncated from the Chebyshev series and also used to express the order of the SCS method. Assume that the interval of the independent variable over which a solution is desired, [ 0 , ], is divided into equal subintervals , ] producing the approximations ( ), = 1, 2, . . . , . Thus, the true solution ( ) can be approximated at + 1 equidistance values of , ( 0 , 1 , . . . , ) so that the step size ℎ is given by ℎ = ( − 0 )/ and hence we can define = 0 + ℎ with 0 ≤ ≤ and 0 ≤ ≤ .
The main aim is to determine the series coefficients for = 0, 1, . . . , over each subinterval. These coefficients can be expressed in terms of the solution derivative , , . . . , ( ) at the beginning of the integration step, = 0. Let us start differentiation term by term of both sides of (15) ℎ times at = 0, to get the recurrence relation: where the high derivatives of ( ) (0) appearing in the right hand side of (17) can be estimated directly from (12) and (13) at = 0. On the other hand, differentiating (5) of the given system ( − 1) ℎ times at = 0 yields Take the approximation ( ) (0) ≈ ( ) (0) by equating (17) and (18) and with some algebraic simplifications; we can get the following recurrence relations in reverse order: . . .
Equation (19) can be compressed into one recurrence relation: Now, the desired coefficients , = 1, 2, . . . , , are estimated and the remaining coefficient 0 can be determined also from (17) in terms of both the estimated coefficients and the Chebyshev polynomials as Substituting these estimated coefficients into (15) gives the required approximation ( ). The algorithm that computes such approximation is in Appendix for constant matrix .

Stability of the SCS Method
To investigate the stability of the proposed method, we consider the model problem that is generally used to test the performance of various numerical methods with the exact solution ( ) = 0 ; is negative real. During the iteration process, the approximate solution depends basically on the step size at each step. To guarantee convergence of the Chebyshev series to the solution from which its coefficients were computed, it is essential to place additional conditions on the step sizeℎ. Such study of the stability of the model equation can be easily transferred to study the stability of a system of differential equations [17]. Now, we apply the SCS method of order five (as a simple case to be generalized later) to the model equation over the interval [0, ℎ]. Estimation of the high derivatives of the test equation (22) As explained in Section 3, by taking the approximation (17) and (23) Solving the algebraic system (24) gives the coefficient values , = 1, . . . , 5 The coefficient 0 can be estimated from (21). Taking into account (11), the required solution ( 1 ) at 1 = 0 + ℎ, 0 = 0, can be obtained from (15).
where ( ℎ) is the stability function It is clear that the stability function represents a good approximation to the exponential function ℎ (exact solution) only if the combination of ℎ is small in particular, if | ℎ| < 1 or ℎ < 1/ . Hence, we can write ℎ = (ℎ) + O(ℎ 6 ). Since the terms omitted from the Chebyshev series constitute the truncation error, the absolute local truncation error associated with the approximate solution is = | ( 1 ) − ( 1 )| = O(ℎ 6 ) so that this error is of the order of ℎ 6 . As we take ℎ small, as the local truncation error can be reduced, for a single iteration. In addition, adding more terms to the truncated series should produce results with comparable accuracy and also can increase the value of ℎ. A procedure for stepping from one value of to another, that is from 0 to 0 + ℎ, can be applied for stepping from to +1 = + ℎ. Thus, (27) can be generalized to obtain the solution ( +1 ) Solving the difference equation (29) gives For large values of , we can take the approximation ( ) ≈ . For negative value of with | ( )| < 1, we can conclude that ( ) → 0 as → ∞. Therefore, the proposed method is stable. The proposed method has been applied for the model problem of = −5 and ℎ = 0.1. In this case, Figure 1 shows the relative errors associated with the method for different orders = 5, 7, 10, 12, 15, and 20. The figure indicates that, as the order decreases, the relative error increases.

Results and Discussion
To verify the accuracy and the efficiency of the proposed method presented in this study for solving the point kinetics equations, we compare the SCS results with those obtained analytically and numerically. The first test case is a preliminary test case of a step reactivity insertion which aims to test the accuracy of the proposed method where the analytical solution is known. The second and third test cases are linear systems of ramp and sinusoidal reactivities, respectively. The fourth test case is a nonlinear system where the reactivity is dependent on the neutron density. The fifth test case is a nonlinear system of Newtonian temperature feedback reactivity. Finally, the SCS method is applied to solve a complex system of thermal hydraulic feedback.

Test Case 1 (Step Reactivity in Thermal Reactor).
This problem is a preliminary test case of a linear system of thermal reactor of six delayed precursor groups, without feedback. Such study aims to testing the accuracy of the SCS results where the analytical solution is available. The kinetic parameters are shown in Table 1 [18].
In this model the step reactivity makes the system linear so that the system should be asymptotically stable at the origin if and only if 0 < 0; see [19]. The exact solutions used in Table 2 are estimated via the eigenvalues and the corresponding eigenvectors of the matrix [20]. In Table 2

Test Case 2 (Ramp Reactivity).
A ramp reactivity input, ( ) = 0.1 , is considered with kinetic parameters listed in Table 1. Table 3 shows a comparison between our tenthorder SCS method, stiffness confinement method SCM, −weighting method [21], generalized Runge-Kutta method GRK [22], the analytical inversion method AIM [23], the piecewise constant approximation PCA [24], the numerical algorithm CORE [25], the better basis function BBF, Hermit polynomial methods [26], the generalized analytical exponential method, GAEM and Pade' approximation with treatment of inhour equation root [27], the efficient technique ET [28], the power series method PWS [10], the converged accelerated Taylor series CATS [12], the accurate solution   [29], the enhancement of the piecewise constant approximation EPCA [8], the fundamental backward Euler difference scheme BEFD [13], the integral formulation Taylor series expansions ITS2 [2], the Haar wavelet operational method HWOM [30], the modified exponential time differencing ETD [7], the trigonometric Fourier-series solutions TFS [6], and the treatment theta method TTM [31]. Thus, the proposed method is compared with twenty of different numerical methods in Table 3. The estimated values of the neutron density are rounded to 9 decimal places for 0 ≤ ≤ 11 consuming computational time 0.16 s and number of steps of 1100. The table shows that the best agreement with our SCS results occurs with the reference solutions of CATS, BEFD, and TFS methods for the all 9 digits accuracy.

Test Case 3 (Sinusoidal Reactivity).
In this case, we will examine a sinusoidal reactivity of the form [8,13]: where is the half period and 0 = 8 /(8 + ). In this case, we consider one neutron energy group with the following parameters: = 0.077 −1 , = 0.0079. Calculations are performed using fixed value of the half period, = 50 and 0 = 0.005333. Two values of neutron generation times are considered: the first is Λ = 10 −6 for fast reactor I and the second is Λ = 10 −8 for fast reactor II. In Table 4(a), the results of our SCS method, using step size ℎ = 0.0001, are compared with those of CATS and EPCA methods [8]. Up to = , the results of the SCS and CATS are completely coinciding while EPCA method differ in the last digit at = 50 . Table 4(b) for fast reactor II shows a comparison between the SCS and BEFD [28] results for one full cycle (0 ≤ ≤ 2 ). Table 4(b), all the 9 digits of CATS and the 12 digits of the SCS results are fully compatible. The SCS method for this case of fast reactor II requires a very small step size (ℎ = 1 − 6) consuming about 9.5E-5 s/step. Figures 2(a) and 2(b) show plots of the neutron densities for one full cycle for the two cases of fast reactors I and II. The figures also show the time to the peak and the neutron density at . It is found that, for fast reactor I, the time to the first peak is = 39.1073 and ( ) = 61.515290571while for fast reactor II, = 39.10721 and ( ) = 61.530146020. It should be mentioned that the value of ( ) of the fast reactor II is the same as the value that was provided in [28].

Test Case 4 (The Reactivity Is a Function of the Neutron Density).
In this case, we have considered the reactivity as a function of the neutron density. The reactivity function is of the form ( ) = ±0.1 ( ). The parameters are listed in Table 1. Calculations are performed up to = 15 for both positive and negative reactivity. The accuracy of the results is tested to 10 decimal places. Table 5 shows the results obtained by three of the highly accurate methods; the TFS method of order 10 (ℎ = 0.05) [4], CATS [12], and BEFD [13]; and the results obtained by the SCS method of order 20 (ℎ = 0.1), for positive and negative reactivities. The table indicates that the results of the SCS method and those of CATS and DEFD are identical to the ten digits.
It should be noted that all calculations in this literature are performed by MATLAP on a 2.0 GHz Intel Core i7-2630QM computer while the calculations obtained by the CATS and BEFD methods are performed on a VAIO 2.4 GHz dual core laptop as indicated in Ref. [13]. Indeed, the calculation times obtained by both CATS and BEFD methods are still superior to the SCS method.
where ( ) is the temperature of the reactor, 0 is the initial temperature of the reactor, is the temperature coefficient of reactivity, and is the reciprocal of the thermal capacity of the reactor. The parameters are listed in Table 1. In Figure 3, we have determined the times to peak neutron density and the corresponding peak densities for the step reactivities   where ( ) is the reactor temperature, = 13006. 193586 ∘ C is the heat capacity of the reactor, = 0.23 and = 0.2 are the reactor height and diameter, respectively, = 20 ∘ C is the coolant temperature, = 17.52 is a physical constant for air at standard temperature and pressure, = 1 is the fraction of the energy from fission deposited as heat in the system, and ( ) is the neutron density (neutrons/cm 3 ). The reactivity is given by where ( ) (dollars) is the system reactivity which consists of step term (initial reactivity, 0 = 0.043 ) plus a linear feedback component. Here, is in seconds and the reactivity feedback coefficients = −0.00306 / ∘ C has units of dollars/ ∘ C. The remaining parameters describing this problem are listed in Table 1. Calculations are performed up to 5000 min using the SCS method of order 25 with constant step size of 0.001 consuming computational time 5.2E-3s per step.
Tables 7(a), 7(b), and 7(c) show a comparison between the results of SCS method and the results obtained by both the backward Euler finite difference BEFD [13] and the simple kinetics thermal hydraulics SKINATH (with ℎ = 0.01 ) methods [32].
As shown in the three tables, the results of the SCS are presented with 11 digits to be compared with 9 digits of the BEFD results. Comparisons indicate that the SCS results are identical to all the 9 digits of the BEFD method except one value of the power at = 500 where they differ in the last digit. Also, there is a difference between the SCS results of the reactivity values at = 3000 and 4000 , in the last two digits. On the other hand, additional SCS results for this problem out to a final time of 5000 minutes are presented in Tables 7(a), 7(b), and 7(c) and Figures 4, 5, and 6 for the power, reactivity, and temperature, respectively. Due to the presence of feedback, more than one peak of the neutron  density, reactivity, and temperature are arising. The three figures exhibit a characteristic damped oscillatory behaviour to equilibrium levels.
Future Work. Although the present method provided accurate results in comparison with the two highly accurate methods CATS and BEFD, the last two methods are still superior because the smallness of their CPU times. The technique of using a criterion for choosing an adaptively step size length was applied effectively in a previous work [5,33]. So, our future work is to apply the technique of varying step sizes for the SCS method based on the solution behaviours. Such technique is expected to reduce both calculation times and the number of integration steps. On the other hand, this method can be applied with modifications to solve the modified fractional model of the point kinetics equations [34,35].

Conclusions
A new method based on shifted Chebyshev series is introduced to solve different systems of the point kinetics equations. Special transformation is used to map the original interval of the Chebyshev polynomials from [−1, 1] to [0, ℎ].
It has been proved that such transformation is useful and suitable for the iteration process. By differentiating both the shifted Chebyshev series and the given system th times, we end up with an algebraic system of equations. Solving such a system produced a new recurrence relation used to determine the series coefficients. One advantage of the SCS method is its simplicity whereas the Chebyshev coefficients can be easily determined using a simple recurrence relation. This recurrence relation is very simple to program since it is containing nothing more than the previous estimated coefficients. Therefore, the method is very fast consuming small computational times. Stability of the proposed method is discussed and excellent convergence rate of the Chebyshev approximations is proved. It is proved that the proposed method has exponential rate of convergence. Also, it is proved that as the orders of the method increase, the errors decrease. The present method has been applied to systems which contain six different types of reactivities. The SCS results are compared with the analytical solutions and with some of the conventional methods that are used to solve such systems. The most accurate numerical methods to solve the point kinetics equations are those of CATS and BEFD of Ganapol and TFS of Hamada. The results obtained by the SCS method are in excellent agreement with the results obtained by these methods. The success of the SCS method depends on whether the solutions could be approximated with acceptable accuracy. It is proved that the SCS method succeeded to compute more than 10 correct digits of the solutions. Based on the obtained results in the tables and figures, we conclude that the new method has produced highly accurate results when it has been applied to several linear/nonlinear problems of the point kinetics equations. The accuracy and fast calculations promise a possible applicability of the new method to solve more complex transient problems.