Science and Technology of Nuclear Installations

Volume 2018, Article ID 7105245, 15 pages

https://doi.org/10.1155/2018/7105245

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

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

Correspondence should be addressed to Yasser Mohamed Hamada; moc.liamg@adamah.m.ressay

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 aswhere Here, , are vectors of -dimensional real space, is the square point kinetics matrix of order , and is a column vector which contains all the dependent variables.

#### 3. Shifted Chebyshev Polynomials on

Assume that the standard finite Chebyshev series is defined by [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 formThis transformation leads to a new type of the series of polynomialsFrom (10), it can be concluded that The successive Chebyshev polynomials can be defined by the recurrence relationwhich together with and recursively generates all the polynomials very efficiently. The high derivatives of (12) can be deduced by the following recurrence relations:

For the approximation of , we should expect to get a smaller mini–max error with in than with in the larger range according to the following theorem.

Theorem 1 (Lagrange-Chebyshev approximation [16]). *Assume that is the Lagrange polynomial of order that is based on the Chebyshev interpolating nodes on . If, thenholds for all and .*

Theorem 1 shows that the error depends basically on the length of the step size and on the number of the truncated coefficients . It should be noted that the interval minimizes the error better than the interval .

#### 4. Derivation of the SCS Method

Let the approximate solution of system (5) be expressed by the finite shifted Chebyshev series where is defined in (10) as a function of the independent variable and the initial value of can be obtained bywhere the coefficients , , are vectors of order . The subscript refers to the number of the coefficients truncated from the Chebyshev series and also used to express the order of the SCS method. Assume that the interval of the independent variable over which a solution is desired, , is divided into equal subintervals producing the approximations , . Thus, the true solution can be approximated at equidistance values of , so that the step size is given by and hence we can define with and .

The main aim is to determine the series coefficients for over each subinterval. These coefficients can be expressed in terms of the solution derivative at the beginning of the integration step, . Let us start differentiation term by term of both sides of (15) times at , to get the recurrence relation: where the high derivatives of appearing in the right hand side of (17) can be estimated directly from (12) and (13) at . On the other hand, differentiating (5) of the given system times at yields

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 methodswith the exact solution ; is negative real. During the iteration process, the approximate solution depends basically on the step size at each step. To guarantee convergence of the Chebyshev series to the solution from which its coefficients were computed, it is essential to place additional conditions on the step size. Such study of the stability of the model equation can be easily transferred to study the stability of a system of differential equations [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.