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 to . The method deals with the Chebyshev series as a finite difference method not as a spectral method. Stability of the method is discussed and it has proved that the method has an exponential rate of convergence. The method 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.
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 start-up 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  used Magnus expansion to solve the nonlinear PKEs. Sérgio used the ITS2 method to solve the PKEs with adaptive time step  and used the same method to solve the PKEs in the integral formulation with arbitrary number of delayed neutron groups and Newtonian temperature feedback . 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  used analytical exponential model to solve the stochastic PKEs via eigenvalues and eigenvectors. The modified exponential time differencing method  is applied to solve the PKEs. Picca  provided a method based on piecewise constant approximation to solve the nonlinear PKEs. Wollmann da Silva  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)  and the backward Euler difference scheme (BEFD)  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 and finite difference methods. We start by dividing the total interval in the independent variable , over which a solution is desired , into equal subintervals. The SCS method has employed finite Chebyshev series expansion to approximate the exact solutions at equidistance values of . We have transformed the traditional interval of the independent variable of the Chebyshev series from to 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).
2. The Point Kinetics Equations with Temperature Feedback
The system of the point kinetics equations for neutron density and delayed neutron precursor groups , with temperature feedback can be described by the differential equations :where is average neutron flux density, , is groups of delayed neutron precursors, , , is reactivity, where 1.0$ = , is the group decay constant), is the group delayed fraction, , 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 aswhere Here, , are vectors of -dimensional real space, is the square point kinetics matrix of order , and is a column vector which contains all the dependent variables.
3. Shifted Chebyshev Polynomials on
Assume that the standard finite Chebyshev series is defined by 
where is used here to express the usual Chebyshev polynomials of the first kind; the independent variable is defined on ; 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 than the larger range. Therefore, we may map the independent variable in to the variable in by the transformation Such transformation leads to the shifted Chebyshev polynomials of the first kind of the formThis transformation leads to a new type of the series of polynomialsFrom (10), it can be concluded that The successive Chebyshev polynomials can be defined by the recurrence relationwhich together with and 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 than with in the larger range according to the following theorem.
Theorem 1 (Lagrange-Chebyshev approximation ). Assume that is the Lagrange polynomial of order that is based on the Chebyshev interpolating nodes on . If, thenholds for all and .
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 minimizes the error better than the interval .
4. 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 bywhere the coefficients , , are vectors of order . 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, , is divided into equal subintervals producing the approximations , . Thus, the true solution can be approximated at equidistance values of , so that the step size is given by and hence we can define with and .
The main aim is to determine the series coefficients for over each subinterval. These coefficients can be expressed in terms of the solution derivative at the beginning of the integration step, . Let us start differentiation term by term of both sides of (15) times at , to get the recurrence relation: where the high derivatives of appearing in the right hand side of (17) can be estimated directly from (12) and (13) at . On the other hand, differentiating (5) of the given system times at yields
Equation (19) can be compressed into one recurrence relation:
Now, the desired coefficients , , are estimated and the remaining coefficient 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 .
5. 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 methodswith the exact solution ; 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 .
Now, we apply the SCS method of order five (as a simple case to be generalized later) to the model equation over the interval . Estimation of the high derivatives of the test equation (22) at gives
Solving the algebraic system (24) gives the coefficient values ,
Substitute the estimated coefficients into (26) to get
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 or . Hence, we can write . Since the terms omitted from the Chebyshev series constitute the truncation error, the absolute local truncation error associated with the approximate solution is so that this error is of the order of . 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 to , can be applied for stepping from to . Thus, (27) can be generalized to obtain the solution
Solving the difference equation (29) gives
For large values of , we can take the approximation . For negative value of with , we can conclude that as . Therefore, the proposed method is stable. The proposed method has been applied for the model problem of and . In this case, Figure 1 shows the relative errors associated with the method for different orders , and 20. The figure indicates that, as the order decreases, the relative error increases.
6. 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.
6.1. 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 . Step reactivity insertion, for , is considered. Because the reactivity is constant, the derivatives of the point kinetics matrix appearing in (18) will be vanished; i.e., .
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 ; see . The exact solutions used in Table 2 are estimated via the eigenvalues and the corresponding eigenvectors of the matrix . In Table 2, three different values of step reactivity insertions, , and for prompt subcritical, prompt critical, and prompt supercritical, respectively, have been compared. The table shows a comparison between the SCS and the analytical solutions. The SCS method is applied with high order and constant step size . For the prompt subcritical with , calculations are performed up to and for the prompt critical with , calculations are performed up to while for the prompt supercritical with , calculations are performed up to . The average computational times were 5E-4s/step. It can be seen that the SCS and the analytical results are identical to the 12 digits for all reactivity types except only two values at and .
6.2. Test Case 2 (Ramp Reactivity)
A ramp reactivity input, , is considered with kinetic parameters listed in Table 1. Table 3 shows a comparison between our tenth-order SCS method, stiffness confinement method SCM, weighting method , generalized Runge-Kutta method GRK , the analytical inversion method AIM , the piecewise constant approximation PCA , the numerical algorithm CORE , the better basis function BBF, Hermit polynomial methods , the generalized analytical exponential method, GAEM and Pade' approximation with treatment of inhour equation root , the efficient technique ET , the power series method PWS , the converged accelerated Taylor series CATS , the accurate solution , the enhancement of the piecewise constant approximation EPCA , the fundamental backward Euler difference scheme BEFD , the integral formulation Taylor series expansions ITS2 , the Haar wavelet operational method HWOM , the modified exponential time differencing ETD , the trigonometric Fourier- series solutions TFS , and the treatment theta method TTM . 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 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.
6.3. 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 . In this case, we consider one neutron energy group with the following parameters: , . Calculations are performed using fixed value of the half period, and . Two values of neutron generation times are considered: the first is for fast reactor I and the second is for fast reactor II. In Table 4(a), the results of our SCS method, using step size , are compared with those of CATS and EPCA methods . Up to , the results of the SCS and CATS are completely coinciding while EPCA method differ in the last digit at .
As shown in 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 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 and while for fast reactor II, and . It should be mentioned that the value of of the fast reactor II is the same as the value that was provided in .
(a) The time to the neutron peak and the neutron density at that peak for sinusoidal reactivity for fast reactor I with
(b) The time to the neutron peak and the neutron density at that peak for sinusoidal reactivity for fast reactor II with
6.4. 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 . The parameters are listed in Table 1. Calculations are performed up to 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 , CATS , and BEFD ; and the results obtained by the SCS method of order 20 , 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. . Indeed, the calculation times obtained by both CATS and BEFD methods are still superior to the SCS method.
6.5. Test Case 5 (Temperature Feedback Reactivity)
To check the accuracy and efficiency of the new technique, more comparisons are made with other conventional methods. In this case, we have applied the SCS method to solve the point kinetics equations in the presence of Newtonian temperature feedback reactivity where the reactivity is dependent on the temperature.where is the temperature of the reactor, 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 ($): 1.0, 1.2, 1.5, and 2. Tables 6(a), 6(b), and 6(c) show the SCS results of the neutron densities for the step reactivities of 1.0, 1.5, and 2$, to 12 digits. The SCS method is applied with order 30 and consuming CPU times about per step. Comparisons with the 10 digits of CATS  and 9 digits BEFD  methods show that our results are completely agree with them. For $, the method of EPCA  is included to be compared with SCS, BEFD, and CATS methods. Results of EPCA using agree to only 6 figures for .
6.6. Test Case 6 (Thermal Hydraulics Feedback)
This problem is a stiff nonlinear system with simple thermal hydraulic feedback. The SCS method is applied to solve the point kinetics equations (1) and (2) coupled with the thermal hydraulic feedback equation where is the reactor temperature, is the heat capacity of the reactor, and are the reactor height and diameter, respectively, is the coolant temperature, is a physical constant for air at standard temperature and pressure, is the fraction of the energy from fission deposited as heat in the system, and is the neutron density (neutrons/cm3). The reactivity is given bywhere (dollars) is the system reactivity which consists of step term (initial reactivity, ) plus a linear feedback component. Here, is in seconds and the reactivity feedback coefficients 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 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  and the simple kinetics thermal hydraulics SKINATH (with ) methods .
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 where they differ in the last digit. Also, there is a difference between the SCS results of the reactivity values at , 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].
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 to . 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 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.
See Algorithm 1.
No data were used to support this study. All parameters and constants used for the test cases considered are available in the reference section. The algorithms used in this manuscript are coded by the authors themselves using Matlab 8 for personal computer.
Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this paper.
The author would like to thank Professor Barry Ganapol who provided him with some accurate results of the CATS method.
Q. Sérgio and L. Bogado, “Solution of the point kinetics equations with temperature feedback by ITS2 method,” Prog. Nucl. Energy, vol. 91, pp. 240–249, 2016.View at: Google Scholar
A. A. Nahla and A. M. Edress, “Analytical exponential model for stochastic point kinetics equations via eigenvalues and eigenvectors,” Journal of Nuclear Science and Technology, pp. 20–27, 2016.View at: Google Scholar
B. D. Ganapol et al., The Solution of the Point Kinetics Equation: A Converged Accelerated Taylor Series (CATS), PHYSOR 2012, American Nuclear Society, LaGrange Park, Knoxville, TN, USA, 2012.
D. L. Hetrick, Dynamics of Nuclear Reactors, American Nuclear Society, JBC, Illinois, IL, USA, 1993.
J. C. Mason and D. C. Handscomb, Chebyshev Polynomials, CRC Press, Boca Raton, FL, USA, 2003.View at: MathSciNet
G. M. M. Phillips and J. P. Taylor, Theory and Application of Numerical Analysis, Elsevier Science & Technology Books, 2nd edition, 1996.
E. Kendall, W. H. Atkinson, and E. S. David, Numerical Solution of Ordinary Differential Equations, John Wiley & Sons, Inc, Hoboken, NJ, USA, 2009.
J. Nobrega, “A new solution of the point kinetics equations,” Nuclear Science and Engineering, vol. 46, pp. 366–375, 1971.View at: Google Scholar
Y. M. Hamada, “Liapunov's stability on autonomous nuclear dynamical systems,” Progress in Nuclear Energy, vol. 73, pp. 11–20, 2014.View at: Google Scholar
E. S. Lawrence, J. I. Arnold, and H. F. Stephen, Elementary Linear Algebra: A Matrix Approach, Pearson Education, Inc, Upper Saddle River, NJ, USA, 2nd edition, 2008.
A. E. Aboanber and A. A. Nahla, “Solution of the point kinetics equations in the presence of Newtonian temperature feedback by Padé approximations via the analytical inversion method,” Journal of Physics A: Mathematical and General, vol. 35, no. 45, pp. 9609–9627, 2002.View at: Publisher Site | Google Scholar | MathSciNet
L. Haofeng, W. Chen, and Q. Zhu, “A new integral for solving the point reactor neutron kinetics equations,” Annals of Nuclear Energy, vol. 36, pp. 427–432, 2009.View at: Google Scholar