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

Authors: Yasser Mohamed Hamada
Journal: Science and Technology of Nuclear Installations
 (2018)

---

## 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[-1,1] to [0,h]. 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.


---

## Body

## 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 variablet, over which a solution is desired [t0,tn], into n equal subintervals. The SCS method has employed finite Chebyshev series expansion to approximate the exact solutions at n+1 equidistance values of t. We have transformed the traditional interval of the independent variable of the Chebyshev series from [-1,1] to [0,h] to be suitable for the iteration process. Hence, with large values of the truncated coefficients and suitable choice of h, 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 densityN(t) and G-delayed neutron precursor groups Cit, i=1,2,…,G with temperature feedback can be described by the differential equations [14]:(1)dNtdt=ρt,N,T-βΛNt+∑i=1GλiCit(2)dCitdt=βiΛNt-λiCit(3)dTtdt=EtNt+UtTt(4)ρt=Ft,N,Twhere N(t) is average neutron flux density, (neutron/cm3), Ci(t) is Gth groups of delayed neutron precursors, (i=1,2,…,G), (atom/cm3), ρ(t) is reactivity, where 1.0$ = ρ/β, λi is the ith group decay constant(s-1), βi is the ith group delayed fraction, β=∑i=1Gβi, Λ is the neutron generation time (s), kc is the reciprocal of the reactor heat capacity, and G is the total number of delayed neutron groups.The functionsEt and Ut may be constants or functions of the independent variable. In general, the reactivity ρt=Ft,N,T 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(5)y′t=At·yt,y0=y0where (6)
y
′
t
=
d
d
t
N
t
C
1
t
C
2
t
⋮
C
G
t
T
t
,
y
t
=
N
t
C
1
t
C
2
t
⋮
C
G
t
T
t
,
y
0
=
N
0
N
0
β
1
Λ
λ
1
N
0
β
2
Λ
λ
2
⋮
N
0
β
G
Λ
λ
G
T
0
,
A
=
ρ
t
,
N
,
T
-
β
Λ
λ
1
λ
2
⋯
λ
G
0
β
1
Λ
-
λ
1
0
⋯
0
0
⋮
⋮
⋮
⋱
⋮
⋮
β
G
Λ
0
0
⋯
-
λ
G
0
E
t
0
0
⋯
0
U
tHere, yt, y0(t) are vectors of (G+2)-dimensional real space, A(t) is the square point kinetics matrix of order (G+2)×(G+2), and y(t) is a column vector which contains all the dependent variables.


## 3. Shifted Chebyshev Polynomials on[0,h]

Assume that the standard finite Chebyshev series is defined by [15](7)yυx=∑k=0υakTk∗xwhereTk∗(x) is used here to express the usual Chebyshev polynomials of the first kind; the independent variable x 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 t 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,h] than the larger range[-1,1]. Therefore, we may map the independent variable x in [-1,1] to the variable t in [0,h] by the transformation (8)x=2t-hhSuch transformation leads to the shifted Chebyshev polynomials of the first kind Tk(t) of the form(9)Tkt=Tk2t-hh,k=0,1,…,υThis transformation leads to a new type of the series of polynomials(10)T0t=1T1t=2th-1T2t=8t2h2-8th+1T3t=32t3h3-48t2h2+18th-1T4t=128t4h4-256t3h3+160t2h2-32th+1T5t=512t5h5-1280t4h4+1120t3h3-400t2h2+50th-1⋮From (10), it can be concluded that (11)Tk0=-1kTkh=1,k=0,1,2,…The successive Chebyshev polynomials can be defined by the recurrence relation(12)Tk+1t=22th-1Tkt-Tk-1twhich together with T0t=1 and T1t=2t/h-1 recursively generates all the polynomials Tk(t) very efficiently. The high derivatives of (12) can be deduced by the following recurrence relations:(13)Tk+1it=22th-1Tkit+4ihTki-1-Tk-1itFor the approximation ofy(t), we should expect to get a smaller mini–max error with Tkt in 0≤t≤h than with Tk∗x in the larger range -1≤x≤1 according to the following theorem.Theorem 1 (Lagrange-Chebyshev approximation [16]).
Assume thatyυt is the Lagrange polynomial of order υ that is based on the Chebyshev interpolating nodes on [a,b]. Ifyυt∈Cυ+1[a,b], then(14)Et=yext-yυt=2b-aυ+14υ+1υ+1!maxa<ξ<b⁡yυ+1ξholds for all ξ∈a,b and t∈a,b.Theorem1 shows that the error depends basically on the length of the step size h=b-a and on the number of the truncated coefficients υ. It should be noted that the interval [0,h] minimizes the error E(t) better than the interval [-1,1].


## 4. Derivation of the SCS Method

Let the approximate solutionyυt of system (5) be expressed by the finite shifted Chebyshev series (15)yυt=∑i=0υaiTitwhere Ti(t) is defined in (10) as a function of the independent variable t and the initial value of yυt can be obtained by(16)yυ0=a0-a1+a2-⋯+aυwhere the coefficients ai, i=1,2,…,υ, are vectors of order (G+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 t over which a solution is desired, t0,tf, is divided into n equal subintervals t0,t1,…,ti,ti+1,…,[tn-1,tn] producing the approximations yυti, i=1,2,…,n. Thus, the true solution y(t) can be approximated at n+1 equidistance values of t, t0,t1,…,tn so that the step size h is given by h=tn-t0/n and hence we can define ti=t0+ih with 0≤i≤n and t0≤ti≤tn.The main aim is to determine the series coefficientsai for i=0,1,…,υ over each subinterval. These coefficients can be expressed in terms of the solution derivative yυ′,yυ′′,…,yυ(υ) at the beginning of the integration step, t=0. Let us start differentiation term by term of both sides of (15)  υth times at t=0, to get the recurrence relation: (17)yυk0=∑i=0υaiTik0,k=1,2,…,υwhere the high derivatives of Tik(0) appearing in the right hand side of (17) can be estimated directly from (12) and (13) at t=0. On the other hand, differentiating (5) of the given system (υ-1)th times at t=0 yields(18)yk0=∑i=0k-1k-1iAi0yk-i-10,k=1,2,…,υTake the approximationyk(0)≈yυk(0) by equating (17) and (18) and with some algebraic simplifications; we can get the following recurrence relations in reverse order:(19)aυ=1Tυυ0yυ0aυ-1=1Tυ-1υ-10yυ-10-aυTυυ-10aυ-2=1Tυ-2υ-20yυ-20-aυ-1Tυ-1υ-20-aυTυυ-20⋮a1=1T110y10-a2T210-a3T310-⋯-aυTυ10Equation (19) can be compressed into one recurrence relation: (20)ai=1Tii0yi0-∑l=i+1l≤υυalTli0,i=1,2,…,υNow, the desired coefficientsai, i=1,2,…,υ, are estimated and the remaining coefficient a0 can be determined also from (17) in terms of both the estimated coefficients ai and the Chebyshev polynomials as(21)a0=y0-∑i=1υaiTi0Substituting these estimated coefficients into (15) gives the required approximation yυ(t). The algorithm that computes such approximation is in Appendix for constant matrix A.


## 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(22)y′t=λyt,y0=y0with the exact solution yt=y0eλt; λ 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 sizeh. 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,h]. Estimation of the high derivatives of the test equation (22) at t=0 gives(23)yk0=λky0,k=1,2,…,υAs explained in Section3, by taking the approximation yk(0)≈yυk(0) of (17) and (23), with some algebraic simplifications, we get the algebraic system(24)2-818-3250016-96320-80000192-1536672000030720-30720000061440a1a2a3a4a5=λhy0λ2h2y0λ3h3y0λ4h4y0λ5h5y0Solving the algebraic system (24) gives the coefficient values ai, i=1,…,5(25)a5=λh561440y0a4=λh43072y0+λh56144y0a3=λh3192y0+λh4384y0+3λh54096y0a2=λh216y0+λh332y0+7λh4768y0+λh5512y0a1=λh2y0+λh24y0+5λh364y0+7λh4384y0+7λh52048y0The coefficienta0 can be estimated from (21). Taking into account (11), the required solution yυt1 at t1=t0+h, t0=0, can be obtained from (15).(26)yυt1=a0+a1+a2+⋯+a5Substitute the estimated coefficients into (26) to get (27)yυh=y0+2a1+2a3+2a5=1+λh1!+λh22!+λh33!+λh44!+λh55!y0=φλhy0whereφλh is the stability function (28)φλh=1+λh1!+λh22!+λh33!+λh44!+λh55!It is clear that the stability function represents a good approximation to the exponential functioneλh (exact solution) only if the combination of λh is small in particular, if λh<1 or h<1/λ. Hence, we can write eλh=yυh+O(h6). Since the terms omitted from the Chebyshev series constitute the truncation error, the absolute local truncation error associated with the approximate solution is E=yt1-yυt1=O(h6) so that this error is of the order of h6. As we take h 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 h. A procedure for stepping from one value of  t to another, that is from t0 to t0+h, can be applied for stepping from ti to ti+1=ti+h. Thus, (27) can be generalized to obtain the solution yυti+1(29)yυti+1=φλhyυtiSolving the difference equation (29) gives(30)yυti=φλihyυt0=φλtiy0For large values ofυ, we can take the approximation φλti≈eλti. For negative value of  λ with φλti<1, we can conclude that yυti→0 as i→∞. Therefore, the proposed method is stable. The proposed method has been applied for the model problem of λ=-5 and h=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.Figure 1
Relative errors associated with different method ordersν=5,7,10,12,15, and 20 for the model problem using λ=-5 and h=0.1.


## 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 Table1 [18]. Step reactivity insertion, ρt=ρ0 for t≥0, is considered. Because the reactivity is constant, the derivatives of the point kinetics matrix A(t) appearing in (18) will be vanished; i.e., A′0=A′′0=⋯=0.Table 1
Six precursor group parameters for the test cases considered.
Test Cases 1, 2 Test Case 4 Test Case 5 Test Case 6 i λ i  ( s - 1 ) β i λ i  ( s - 1 ) β i λ i  ( s - 1 ) β i λ i  ( s - 1 ) β i 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 β = 0.007, β = 0.0065, β = 0.00645, β = 0.0065, Λ = 0.00002  s Λ = 0.0005  s Λ = 0.00005  s Λ = 0.00005  s α = 5 × 10 - 5  K - 1, k c = 0.05  k / M W sIn 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 A [20]. In Table 2, three different values of step reactivity insertions, ρ=0.003,0.007, and 0.008 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 υ=50 and constant step size h=0.001. For the prompt subcritical with ρ=0.003, calculations are performed up to 100s and for the prompt critical with ρ=0.007, calculations are performed up to 50s while for the prompt supercritical with ρ=0.008, calculations are performed up to 10s. 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 ρ=0.003 and ρ=0.008.Table 2
Comparison of the exact and SCS results for test case 1 of step reactivity insertions.
ρ t ( s ) 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,ρt=0.1βt, 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 0≤t≤11s 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 for0.1βt ramp reactivity (moderately fast reactor).
Methods h 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 V a r i a b l e a 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 a h 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]:(31)ρt=ρ0sin⁡πtTwhere T is the half period and ρ0=8β/8+λT. In this case, we consider one neutron energy group with the following parameters: λ=0.077s-1, β=0.0079. Calculations are performed using fixed value of the half period, T=50s 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 h=0.0001, are compared with those of CATS and EPCA methods [8]. Up to t=T, the results of the SCS and CATS are completely coinciding while EPCA method differ in the last digit at t=50s.Table 4  (a)
Results of the SCS, EPCA, and CATS methods for sinusoidal reactivity of fast reactor IΛ=10-6
t ( s ) 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Λ=10-8
t ( s ) 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+01Table4(b) for fast reactor II shows a comparison between the SCS and BEFD [28] results for one full cycle (0≤t≤2T).As shown in Table4(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 (h=1E-6) consuming about 9.5E-5 s/step.Figures2(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 tp and the neutron density at tp. It is found that, for fast reactor I, the time to the first peak is tp=39.1073s and ntp=61.515290571while for fast reactor II, tp=39.10721s and ntp=61.530146020. It should be mentioned that the value of ntp of the fast reactor II is the same as the value that was provided in [28].Figure 2 (a)
The time to the neutron peak and the neutron density at that peak for sinusoidal reactivity for fast reactor I withΛ=10-6 (b)
The time to the neutron peak and the neutron density at that peak for sinusoidal reactivity for fast reactor II withΛ=10-8


### 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ρt=±0.1βN(t). The parameters are listed in Table 1. Calculations are performed up to t=15s 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 (h=0.05) [4], CATS [12], and BEFD [13]; and the results obtained by the SCS method of order 20 (h=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.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.
ρ($) t ( s ) BEFD CATS TFS methodh=0.05 SCS methodh=0.1 +0.1βN(t) 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βN(t) 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  T C a  0.90s 0.078s 8.45s 1.12s aTotal 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.(32)ρt=ρ0-αTt-T0(33)dTtdt=kcNtwhere T(t) is the temperature of the reactor, T0 is the initial temperature of the reactor, α is the temperature coefficient of reactivity, and kc 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 ρ0($): 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 h=0.005 consuming CPU times about 1E-3s 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 ρ0=2$, the method of EPCA [8] is included to be compared with SCS, BEFD, and CATS methods. Results of EPCA using h=0.005=100Λ agree to only 6 figures for t>20s.Table 6  (a)
Neutron density for step reactivity feedback:ρt=ρ0-α[Tt-T0], with ρ0 = 1$
t ( s ) ρ 0 = 1$ SCS(h=0.01) 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:ρt=ρ0-α[Tt-T0], with ρ0 = 1.5$
t ( s ) ρ 0 = 1.5$ SCS(h=0.01) 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:ρt=ρ0-α[Tt-T0], with ρ0 = 2$
t ( s ) ρ 0 = 2$ SCS(h=0.01) 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+00Figure 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](34)dTtdt=1CpεNt-AHD0.75Tt-Tc1.25Tt-0.25where T(t) is the reactor temperature, Cp=13006.193586Erg°C is the heat capacity of the reactor, H=0.23m and D=0.2m are the reactor height and diameter, respectively, Tc=20°C is the coolant temperature, A=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 N(t) is the neutron density (neutrons/cm3). The reactivity is given by(35)ρt=ρ0+βTTt-Tcwhere ρt (dollars) is the system reactivity which consists of step term (initial reactivity, ρ0=0.043β) plus a linear feedback component. Here, t is in seconds and the reactivity feedback coefficients βT=-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.Tables7(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 h=0.01s) methods [32].Table 7  (a)
Reactor power (watts): comparison between the SCS, BEFD, and SKINATH results
t ( m i n ) 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
t ( m i n ) 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
t ( m i n ) 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+01As 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 att=500min where they differ in the last digit. Also, there is a difference between the SCS results of the reactivity values at t=3000and4000min, 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].


## 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 Table1 [18]. Step reactivity insertion, ρt=ρ0 for t≥0, is considered. Because the reactivity is constant, the derivatives of the point kinetics matrix A(t) appearing in (18) will be vanished; i.e., A′0=A′′0=⋯=0.Table 1
Six precursor group parameters for the test cases considered.
Test Cases 1, 2 Test Case 4 Test Case 5 Test Case 6 i λ i  ( s - 1 ) β i λ i  ( s - 1 ) β i λ i  ( s - 1 ) β i λ i  ( s - 1 ) β i 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 β = 0.007, β = 0.0065, β = 0.00645, β = 0.0065, Λ = 0.00002  s Λ = 0.0005  s Λ = 0.00005  s Λ = 0.00005  s α = 5 × 10 - 5  K - 1, k c = 0.05  k / M W sIn 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 A [20]. In Table 2, three different values of step reactivity insertions, ρ=0.003,0.007, and 0.008 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 υ=50 and constant step size h=0.001. For the prompt subcritical with ρ=0.003, calculations are performed up to 100s and for the prompt critical with ρ=0.007, calculations are performed up to 50s while for the prompt supercritical with ρ=0.008, calculations are performed up to 10s. 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 ρ=0.003 and ρ=0.008.Table 2
Comparison of the exact and SCS results for test case 1 of step reactivity insertions.
ρ t ( s ) 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,ρt=0.1βt, 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 0≤t≤11s 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 for0.1βt ramp reactivity (moderately fast reactor).
Methods h 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 V a r i a b l e a 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 a h 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]:(31)ρt=ρ0sin⁡πtTwhere T is the half period and ρ0=8β/8+λT. In this case, we consider one neutron energy group with the following parameters: λ=0.077s-1, β=0.0079. Calculations are performed using fixed value of the half period, T=50s 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 h=0.0001, are compared with those of CATS and EPCA methods [8]. Up to t=T, the results of the SCS and CATS are completely coinciding while EPCA method differ in the last digit at t=50s.Table 4  (a)
Results of the SCS, EPCA, and CATS methods for sinusoidal reactivity of fast reactor IΛ=10-6
t ( s ) 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Λ=10-8
t ( s ) 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+01Table4(b) for fast reactor II shows a comparison between the SCS and BEFD [28] results for one full cycle (0≤t≤2T).As shown in Table4(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 (h=1E-6) consuming about 9.5E-5 s/step.Figures2(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 tp and the neutron density at tp. It is found that, for fast reactor I, the time to the first peak is tp=39.1073s and ntp=61.515290571while for fast reactor II, tp=39.10721s and ntp=61.530146020. It should be mentioned that the value of ntp of the fast reactor II is the same as the value that was provided in [28].Figure 2 (a)
The time to the neutron peak and the neutron density at that peak for sinusoidal reactivity for fast reactor I withΛ=10-6 (b)
The time to the neutron peak and the neutron density at that peak for sinusoidal reactivity for fast reactor II withΛ=10-8


## 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ρt=±0.1βN(t). The parameters are listed in Table 1. Calculations are performed up to t=15s 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 (h=0.05) [4], CATS [12], and BEFD [13]; and the results obtained by the SCS method of order 20 (h=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.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.
ρ($) t ( s ) BEFD CATS TFS methodh=0.05 SCS methodh=0.1 +0.1βN(t) 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βN(t) 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  T C a  0.90s 0.078s 8.45s 1.12s aTotal 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.(32)ρt=ρ0-αTt-T0(33)dTtdt=kcNtwhere T(t) is the temperature of the reactor, T0 is the initial temperature of the reactor, α is the temperature coefficient of reactivity, and kc 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 ρ0($): 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 h=0.005 consuming CPU times about 1E-3s 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 ρ0=2$, the method of EPCA [8] is included to be compared with SCS, BEFD, and CATS methods. Results of EPCA using h=0.005=100Λ agree to only 6 figures for t>20s.Table 6  (a)
Neutron density for step reactivity feedback:ρt=ρ0-α[Tt-T0], with ρ0 = 1$
t ( s ) ρ 0 = 1$ SCS(h=0.01) 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:ρt=ρ0-α[Tt-T0], with ρ0 = 1.5$
t ( s ) ρ 0 = 1.5$ SCS(h=0.01) 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:ρt=ρ0-α[Tt-T0], with ρ0 = 2$
t ( s ) ρ 0 = 2$ SCS(h=0.01) 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+00Figure 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](34)dTtdt=1CpεNt-AHD0.75Tt-Tc1.25Tt-0.25where T(t) is the reactor temperature, Cp=13006.193586Erg°C is the heat capacity of the reactor, H=0.23m and D=0.2m are the reactor height and diameter, respectively, Tc=20°C is the coolant temperature, A=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 N(t) is the neutron density (neutrons/cm3). The reactivity is given by(35)ρt=ρ0+βTTt-Tcwhere ρt (dollars) is the system reactivity which consists of step term (initial reactivity, ρ0=0.043β) plus a linear feedback component. Here, t is in seconds and the reactivity feedback coefficients βT=-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.Tables7(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 h=0.01s) methods [32].Table 7  (a)
Reactor power (watts): comparison between the SCS, BEFD, and SKINATH results
t ( m i n ) 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
t ( m i n ) 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
t ( m i n ) 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+01As 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 att=500min where they differ in the last digit. Also, there is a difference between the SCS results of the reactivity values at t=3000and4000min, 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[-1,1] to [0,h]. 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.



---

*Source: 7105245-2018-10-01.xml*