Science and Technology of Nuclear Installations
Volume 2018 (2018), Article ID 7105245, 15 pages
https://doi.org/10.1155/2018/7105245

Research Article

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

Yasser Mohamed Hamada

Department of Basic Science, Faculty of Computers & Informatics, Suez Canal University, Ismailia 41522, Egypt

Correspondence should be addressed to Yasser Mohamed Hamada; yasser.m.hamada@gmail.com

Received 11 May 2018; Revised 27 August 2018; Accepted 2 September 2018; Published 1 October 2018

Academic Editor: Alejandro Clausse

Copyright © 2018 Yasser Mohamed Hamada. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

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.

1. 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 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 [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 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 [14]:

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 as

where

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 [15]

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 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 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 [16]). Assume that is the Lagrange polynomial of order that is based on the Chebyshev interpolating nodes on . If, then

holds 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 by

where 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

Take the approximation 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 , , 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 methods

with 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 [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 . Estimation of the high derivatives of the test equation (22) at gives

As explained in Section 3, by taking the approximation of (17) and (23), with some algebraic simplifications, we get the algebraic system

Solving the algebraic system (24) gives the coefficient values ,

The coefficient can be estimated from (21). Taking into account (11), the required solution at , , can be obtained from (15).

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.

Figure 1: Relative errors associated with different method orders , and 20 for the model problem using and .

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 [18]. 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., .

Table 1: Six precursor group parameters for the test cases considered.

Test Cases 1, 2
Test Case 4
Test Case 5
Test Case 6

1
0.0127
2.6600E-4
0.0124
2.150E-4
0.0124
2.10E-4
0.0124
2.20E-4

2
0.0317
1.4910E-3
0.0305
1.424E-3
0.0305
1.41E-3
0.0305
1.42E-3

3
0.1150
1.3160E-3
0.1110
1.274E-3
0.1110
1.27E-3
0.111
1.27E-3

4
0.3110
2.8490E-3
0.3010
2.568E-3
0.3010
2.55E-3
0.301
2.57E-3

5
1.4000
8.9600E-4
1.1400
7.480E-4
1.1300
7.40E-4
1.14
7.50E-4

6
3.8700
1.8200E-4
3.0100
2.730E-4
3.0000
2.70E-4
3.01
2.70E-4

,
,
,
,

,

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 [19]. The exact solutions used in Table 2 are estimated via the eigenvalues and the corresponding eigenvectors of the matrix [20]. 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 .

Table 2: Comparison of the exact and SCS results for test case 1 of step reactivity insertions.

SCS method results
Theoretical results

0.003
0.20
1.851268287511E+00
1.851268287511E+00

0.40
1.947593411379E+00
1.947593411379E+00

0.60
2.037922055925E+00
2.037922055925E+00

0.80
2.124831636423E+00
2.124831636423E+00

1.00
2.209840456983E+00
2.209840456983E+00

5.00
4.111945003689E+00
4.111945003689E+00

10.0
8.019199973227E+00
8.019199973227E+00

20.0
2.829739978101E+01
2.829739978101E+01

30.0
9.774957522848E+01
9.774957522848E+01

40.0
3.364929691090E+02
3.364929691090E+02

50.0
1.157602990957E+03
1.157602990957E+03

100
5.573863220601E+05
5.573863220602E+05

0.007
0.01
4.508858486352E+00
4.508858486352E+00

0.05
1.961041160367E+01
1.961041160367E+01

0.10
4.522304458456E+01
4.522304458456E+01

0.20
1.597257698629E+02
1.597257698629E+02

0.30
5.188964648308E+02
5.188964648308E+02

0.40
1.667287692551E+03
1.667287692551E+03

0.50
5.345887612045E+03
5.345887612045E+03

1.00
1.805731634228E+06
1.805731634228E+06

5.00
3.053491380444E+26
3.053491380444E+26

10.0
5.888108082190E+51
5.888108082190E+51

50.0
1.125662952084E+254
1.125662952084E+254

0.008
0.01
6.202853575094E+00
6.202853575094E+00

0.05
9.564079408383E+01
9.564079408383E+01

0.10
1.410421803591E+03
1.410421803591E+03

0.20
2.779627082462E+05
2.779627082462E+05

0.30
5.460343932944E+07
5.460343932944E+07

0.40
1.072625143566E+10
1.072625143566E+10

0.50
2.107055258349E+12
2.107055258349E+12

1.00
6.163333749908E+23
6.163333749908E+23

5.00
3.303191632501E+115
3.303191632502E+115

10.0
1.514716564142E+230
1.514716564142E+230

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 [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 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.

Table 3: Comparison between the SCS method of order ten and various numerical methods for ramp reactivity (moderately fast reactor).

Methods

Time (s)

0.00
2.00
4.00
6.00
8.00
10.0
11.0

- weighting
1.0E-4
1.0000E+00
1.3382E+00
2.2283E+00
5.5815E+00
4.2781E+01
4.5109E+05
1.7919E+16

GRK
5.6E-2
1.0000E+00
1.3382E+00
2.2284E+00
5.5820E+00
4.2786E+01
4.5115E+05
1.7904E+16

AIM
1.0E-3
1.0000E+00
1.3382E+00
2.2284E+00
5.5820E+00
4.2786E+01
4.5116E+05
1.7916E+16

PCA
1.0E-2
1.0000E+00
1.3382E+00
2.2278E+00
5.5802E+00
4.2772E+01
—
—

CORE
1.0E-4
1.0000E+00
1.3382E+00
2.2285E+00
5.5822E+00
4.2790E+01
4.5144E+05
1.7957E+16

BBF
1.0E-1
1.0000E+00
1.338200E+00
2.228442E+00
5.582052E+00
4.278629E+01
4.504121E+05
—

Hermite Polynomial
1.0E-4
1.0000E+00
1.338200E+00
2.228440E+00
5.582052E+00
4.278630E+01
4.511630E+05
—

GAEM
1.0E-1
1.0000E+00
1.338196E+00
2.228435E+00
5.582036E+00
4.278619E+01
4.511648E+05
—

Pade’
1.0E-3
1.0000E+00
1.338198E+00
2.228436E+00
5.582037E+00
4.278634E+01
4.511685E+05
—

ET
1.0E-3
1.0000E+00
1.3382E+00
2.2284E+00
5.5820E+00
4.2786E+01
4.5116E+05
—

PWS
0.01
1.0000E+00
1.338200050E+00
2.228441896E+00
5.582052448E+00
4.278629573E+01
4.511636239E+05
1.792213607E+16

CATS

1.0000E+00
1.338200050E+00
2.228441897E+00
5.582052449E+00
4.278629573E+01
4.511636239E+05
1.792213607E+16

Accurate Solution
0.001
1.0000E+00
1.3382E+00
2.2284E+00
5.582E+00
4.2786E+01
4.5117E+05
—

EPCA
2.0E-3
1.0000E+00
1.338200050E+00
2.228441897E+00
5.582052450E+00
4.278629574E+01
4.511636240E+05
1.792213608E+16

BEFD
Variable
1.0000E+00
1.338200050E+00
2.228441897E+00
5.582052449E+00
4.278629573E+01
4.511636239E+05
1.792213607E+16

ITS2
Variable
1.0000E+00
1.338200050E+00
2.228441897E+00
5.582052449E+00
4.278629573E+01
4.511636239E+05
1.792213608E+16

HWOM
—
1.0000E+00
1.338200011E+00
2.228441300E+00
5.582043255E+00
4.278584501E+01
4.522953180E+05
1.828233653E+16

TTM
1.0E-5
1.0000E+00
1.338200050E+00
2.228441896E+00
5.582052447E+00
4.278629572E+01
4.511636238E+05
1.792213607E+16

TFS
1.0E-2
1.0000E+00
1.338200050E+00
2.228441897E+00
5.582052449E+00
4.278629573E+01
4.511636239E+05
1.792213607E+16

Modified ETD
1.0E-2
1.0000E+00
1.338200E+00
2.228442E+00
5.582052E+00
4.278630E+01
4.511636E+05
—

SCS
1.0E-2
1.0000E+00
1.338200050E+00
2.228441897E+00
5.582052449E+00
4.278629573E+01
4.511636239E+05
1.792213607E+16

is variable, the total number of iteration was 372, and the smallest time calculation per step was 0.0156 s.

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 [8]. Up to , the results of the SCS and CATS are completely coinciding while EPCA method differ in the last digit at .

Table 4

(a) Results of the SCS, EPCA, and CATS methods for sinusoidal reactivity of fast reactor I

SCS
CATS
EPCA

10
2.065311114165E+00
2.065311114E+00
2.065311114E+00

20
8.852831171498E+00
8.852831171E+00
8.852831171E+00

30
4.063328341702E+01
4.063328342E+01
4.063328342E+01

40
6.134134599119E+01
6.134134599E+01
6.134134599E+01

50
4.609559630821E+01
4.609559631E+01
4.609559630E+01

(b) Results of the SCS and BEFD methods for sinusoidal reactivity of fast reactor II

SCS
BEFD

10
2.065383519141E+00
2.065383519E+00

20
8.854133920594E+00
8.854133921E+00

30
4.064354221749E+01
4.064354222E+01

40
6.135607516744E+01
6.135607517E+01

50
4.610628770211E+01
4.610628770E+01

60
2.912634839580E+01
2.912634840E+01

70
1.895177042303E+01
1.895177042E+01

80
1.393829211020E+01
1.393829211E+01

90
1.253353405688E+01
1.253353406E+01

100
1.544816513660E+01
1.544816514E+01

Table 4(b) for fast reactor II shows a comparison between the SCS and BEFD [28] results for one full cycle .

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 [28].

(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

Figure 2

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 [4], CATS [12], and BEFD [13]; 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.

Table 5: Comparison between the SCS method versus the highly accurate numerical methods: BEFD, CATS, and TFS methods for the case when the reactivity is a function of neutron flux.

($)

BEFD
CATS
TFS method

SCS method

0.1
1.0817840726
1.0817840726
1.0817840726
1.0817840726

0.5
1.1435623844
1.1435623844
1.1435623844
1.1435623844

1.0
1.1676688085
1.1676688085
1.1676688085
1.1676688085

2.0
1.2079218113
1.2079218113
1.2079218113
1.2079218113

5.0
1.3199424643
1.3199424643
1.3199424643
1.3199424643

10
1.5296094953
1.5296094953
1.5296094953
1.5296094953

15
1.8179563592
1.8179563592
1.8179563592
1.8179563592

0.1
0.93345679785
0.93345679785
93345679785
0.93345679785

0.5
0.90728689207
0.90728689207
90728689207
0.90728689207

1.0
0.89815154947
0.89815154947
89815154947
0.89815154947

2.0
0.88404456939
0.88404456939
88404456940
0.88404456939

5.0
0.85296776460
0.85296776460
85296776460
0.85296776460

10
0.81577146903
0.81577146903
81577146903
0.81577146903

15
0.78747716481
0.78747716481
78747716481
0.78747716481

0.90s
0.078s
8.45s
1.12s

Total calculation times in seconds.

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.

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 [12] and 9 digits BEFD [13] methods show that our results are completely agree with them. For $, the method of EPCA [8] is included to be compared with SCS, BEFD, and CATS methods. Results of EPCA using agree to only 6 figures for .

Table 6

(a) Neutron density for step reactivity feedback: , with = 1$

= 1$

SCS
CATS
BEFD

10
1.320385964321E+02
1.3203859643E+02
1.320385964E+02

20
5.169986094674E+01
5.1699860947E+01
5.169986095E+01

30
2.817468536129E+01
2.8174685361E+01
2.817468536E+01

40
1.814632999521E+01
1.8146329995E+01
1.814633000E+01

50
1.277957703240E+01
1.2779577032E+01
1.277957703E+01

60
9.474932500982E+00
9.4749325010E+00
9.474932501E+00

70
7.244477493764E+00
7.2444774937E+00
7.244477494E+00

80
5.646289700168E+00
5.6462897001E+00
5.646289700E+00

90
4.456834254635E+00
4.4568342546E+00
4.456834255E+00

100
3.550102765997E+00
3.5501027660E+00
3.550102766E+00

(b) Neutron density for step reactivity feedback:  , with = 1.5$

= 1.5$

SCS
CATS
BEFD

10
1.079116832006E+02
1.0791168320E+02
1.079116832E+02

20
4.160428127533E+01
4.1604281275E+01
4.160428128E+01

30
2.329893149652E+01
2.3298931496E+01
2.329893150E+01

40
1.530342749492E+01
1.5303427495E+01
1.530342749E+01

50
1.089014314534E+01
1.0890143145E+01
1.089014315E+01

60
8.101031858998E+00
8.1010318590E+00
8.101031859E+00

70
6.182690459164E+00
6.1826904592E+00
6.182690459E+00

80
4.793307819639E+00
4.7933078196E+00
4.793307820E+00

90
3.755614629341E+00
3.7556146293E+00
3.755614629E+00

100
2.966074951555E+00
2.9660749516E+00
2.966074952E+00

(c) Neutron density for step reactivity feedback: , with = 2$

= 2$

SCS
CATS
BEFD
EPCA

10
1.033808535461E+02
1.0338085355E+02
1.033808535E+02
1.03380957E+02

20
3.913886902261E+01
3.9138869026E+01
3.913886903E+01
3.91388504E+01

30
2.200377720638E+01
2.2003777206E+01
2.200377721E+01
2.20037934E+01

40
1.449367192746E+01
1.4493671927E+01
1.449367193E+01
1.44936820E+01

50
1.031861108271E+01
1.0318611083E+01
1.031861108E+01
1.03186180E+01

60
7.663319203219E+00
7.6633192032E+00
7.663319203E+00
7.66331816E+00

70
5.829395377947E+00
5.8293953779E+00
5.829395378E+00
5.82939902E+00

80
4.499427072767E+00
4.4994270728E+00
4.499427073E+00
4.49942978E+00

90
3.507422662628E+00
3.5074226626E+00
3.507422663E+00
3.50742470E+00

100
2.755126886406E+00
2.7551268864E+00
2.755126886E+00
2.75512843E+00

Figure 3: Peak and time to the peak of the flux for step feedback with ρ = 1$, 1.2$, 1.5$, and ρ = 2$.

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 [32]

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/cm³). The reactivity is given by

where (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 [13] and the simple kinetics thermal hydraulics SKINATH (with ) methods [32].

Table 7

(a) Reactor power (watts): comparison between the SCS, BEFD, and SKINATH results

Methods

SKINATH
BEFD
SCS

0.000
1.000E-02
1.000000000E-02
1.00000000000E-02

1.000
1.391E-02
1.390664400E-02
1.39066439992E-02

10.00
1.033E-01
1.033074569E-01
1.03307456949E-01

100.0
7.946E-01
7.946137382E-01
7.94613738201E-01

250.0
2.118E-03
2.118680532E-03
2.11868053231E-03

500.0
9.819E-01
9.818854372E-01
9.81885437265E-01

750.0
4.033E+00
4.033472004E+00
4.03347200416E+00

1000
7.680E+00
7.680950791E+00
7.68095079127E+00

1500
—
1.502421669E+01
1.50242166936E+01

2000
—
1.358851928E+01
1.35885192813E+01

2500
—
1.348250322E+01
1.34825032231E+01

3000
—
1.359755709E+01
1.35975570943E+01

4000
—
1.357214030E+01
1.35721402958E+01

5000
—
—
1.35731363004E+01

(b) Reactivity ($): comparison between the SCS, BEFD, and SKINATH results

Methods

SKINATH
BEFD
SCS

0.000
4.300E-02
4.300000000E-02
4.30000000000E-02

1.000
4.230E-02
4.299982607E-02
4.29998260709E-02

10.00
4.299E-02
4.299416065E-02
4.29941606543E-02

100.0
-2.950E-02
-2.949530409E-02
-2.94953040871E-02

250.0
6.388E-03
6.388467468E-03
6.38846746738E-03

500.0
-4.826E-03
-4.825173107E-03
-4.82517310699E-03

750.0
-1.858E-03
-1.857286712E-03
-1.85728671222E-03

1000
1.379E-03
1.379148486E-03
1.37914848596E-03

1500
—
4.962590743E-04
4.96259074364E-04

2000
—
-2.104995975E-04
-2.10499597521E-04

2500
—
3.094249800E-05
3.09424979793E-05

3000
—
3.819718029E-06
3.81971800210E-06

4000
—
5.947391521E-07
5.94739152658E-07

5000
—
—
-3.71867594560E-08

(c) Temperature (°C): comparison between the SCS, BEFD, and SKINATH results

Methods

SKINATH
BEFD
SCS

0.000
2.000E+01
2.000000000E+01
2.00000000000E+01

1.000
2.000E+01
2.000005684E+01
2.00000568396E+01

10.00
2.000E+01
2.000190828E+01
2.00019082829E+01

100.0
4.369E+01
4.369127585E+01
4.36912758455E+01

250.0
3.196E+01
3.196455312E+01
3.19645531152E+01

500.0
3.563E+01
3.562914154E+01
3.56291415382E+01

750.0
3.466E+01
3.465924402E+01
3.46592440236E+01

1000
3.360E+01
3.360158546E+01
3.36015854621E+01

1500
—
3.389011141E+01
3.38901114136E+01

2000
—
3.412107830E+01
3.41210782998E+01

2500
—
3.404217565E+01
3.40421756543E+01

3000
—
3.405103931E+01
3.40510393078E+01

4000
—
3.405209322E+01
3.40520932225E+01

5000
—
—
3.40522997342E+01

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.

Figure 4: Power versus time for thermal hydraulic feedback problem.

Figure 5: Reactivity versus time for thermal hydraulic feedback problem.

Figure 6: Temperature versus time for thermal hydraulic feedback problem.

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].

7. 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 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.

Appendix

See Algorithm 1.

Input:  % Insert the initial values of the vectors and ,

% Define both the step size and the number of iteration ;

% Define the Chebyshev polynomials at the endpoints of the interval , that is and

for acceding to Eq. (11)

% Define at for and acceding to

Eq. (13).

for

;

%Define the successive derivatives of the solution vector according to Eq. (18).

for

end

% If we refer to the derivative of the Chebyshev coefficient by , as an

example, the recurrence relation (19) can be used to estimate the successive values of the

Chebyshev coefficient vectors with a reverse order as follows:

;

;

;

;

;

;

end

% Equation (21) is used to estimate the remaining coefficient .

;

for  

;

end

% Substitute by the estimated coefficients into Eq. (15) to get the required approximation at

the end of the integration step.

;

for  

;

end

;

end

Algorithm 1: The Chebyshev series method of order with constant matrix .

Data Availability

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.

Acknowledgments

The author would like to thank Professor Barry Ganapol who provided him with some accurate results of the CATS method.

References

Y. Cai, X. Peng, Q. Li, and K. Wang, “A numerical solution to the nonlinear point kinetics equations using Magnus expansion,” Annals of Nuclear Energy, vol. 89, pp. 84–89, 2016.

S. Q. Bogado Leite, D. A. P. Palma, M. T. De Vilhena, and B. E. J. Bodmann, “Analytical representation of the solution of the point reactor kinetics equations with adaptive time step,” Progress in Nuclear Energy, vol. 70, pp. 112–118, 2014.

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.

Y. M. Hamada, “Trigonometric Fourier-series solutions of the point reactor kinetics equations,” Nuclear Engineering and Design, vol. 281, pp. 142–153, 2015.

Y. M. Hamada, “Generalized trigonometric Fourier-series method with automatic time step control for non-linear point kinetics equations,” Computational & Applied Mathematics, vol. 36, no. 3, pp. 1–30, 2017.

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.

M. M. A. Razak, K. Devan, and T. Sathiyasheela, “The modified exponential time differencing (ETD) method for solving the reactor point kinetics equations,” Annals of Nuclear Energy, vol. 76, pp. 193–199, 2015.

P. Picca, R. Furfaro, and B. D. Ganapol, “A highly accurate technique for the solution of the non-linear point kinetics equations,” Annals of Nuclear Energy, vol. 58, pp. 43–53, 2013.

M. Wollmann Da Silva, S. Bogado Leite, M. T. Vilhena, and B. E. J. Bodmann, “On an analytical representation for the solution of the neutron point kinetics equation free of stiffness,” Annals of Nuclear Energy, vol. 71, pp. 97–102, 2014.

Y. M. Hamada, “Generalized power series method with step size control for neutron kinetics equations,” Nuclear Engineering and Design, vol. 241, no. 8, pp. 3032–3041, 2011.

Y. M. Hamada, “Confirmation of accuracy of generalized power series method for the solution of point kinetics equations with feedback,” Annals of Nuclear Energy, vol. 55, pp. 184–193, 2013.

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.

B. D. Ganapol, “A highly accurate algorithm for the solution of the point kinetics equations,” Annals of Nuclear Energy, vol. 62, pp. 564–571, 2013.

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.

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.

Y. M. Hamada, “Liapunov's stability on autonomous nuclear dynamical systems,” Progress in Nuclear Energy, vol. 73, pp. 11–20, 2014.

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.

Y.-A. Chao and A. Attard, “Resolution of the stiffness problem of reactor kinetics,” Nuclear Science and Engineering, vol. 90, no. 1, pp. 40–46, 1985.

J. Sanchez, “On the numerical solution of the point reactor kinetics equations by generalized Runge-Kutta methods,” Nuclear Science and Engineering, vol. 103, no. 1, pp. 94–99, 1989.

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.

M. Kinard and E. J. Allen, “Efficient numerical solution of the point kinetics equations in nuclear reactor dynamics,” Annals of Nuclear Energy, vol. 31, no. 9, pp. 1039–1051, 2004.

B. Quintero-Leyva, “CORE: A numerical algorithm to solve the point kinetics equations,” Annals of Nuclear Energy, vol. 35, no. 11, pp. 2136–2138, 2008.

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.

A. A. Nahla, “Analytical solution to solve the point reactor kinetics equations,” Nuclear Engineering and Design, vol. 240, no. 6, pp. 1622–1629, 2010.

A. A. Nahla, “An efficient technique for the point reactor kinetics equations with Newtonian temperature feedback effects,” Annals of Nuclear Energy, vol. 38, no. 12, pp. 2810–2817, 2011.

S. Yamoah, E. H. K. Akaho, and B. J. B. Nyarko, “An accurate solution of point reactor neutron kinetics equations of multi-group of delayed neutrons,” Annals of Nuclear Energy, vol. 54, pp. 104–108, 2013.

A. Patra and S. Saha Ray, “A numerical approach based on Haar wavelet operational method to solve neutron point kinetics equation involving imposed reactivity insertions,” Annals of Nuclear Energy, vol. 68, pp. 112–117, 2014.

A. A. Nahla, “Numerical treatment for the point reactor kinetics equations using theta method, eigenvalues and eigenvectors,” Progress in Nuclear Energy, vol. 85, pp. 756–763, 2015.

H. L. Dodds, Jr. and R. M. Westfall, “SKINATH - a computer program for solving the reactor point kinetics equations with simple thermal-hydraulic feedback,” Tech. Rep. ORNL/CSD/TM-210, 1984.

Y. M. Hamada and M. G. Brikaa, “Nonstandard finite difference schemes for numerical solution of the fractional neutron point kinetics equations,” Annals of Nuclear Energy, vol. 102, pp. 359–367, 2017.

Y. M. Hamada, “Modified fractional neutron point kinetics equations for finite and infinite medium of bar reactor core,” Annals of Nuclear Energy, vol. 106, pp. 118–126, 2017.

Y. M. Hamada, “Solution of the fractional neutron point kinetics equations considering time derivative of the reactivity,” Progress in Nuclear Energy, vol. 98, pp. 153–166, 2017.