International Journal of Nuclear EnergyVolume 2013 (2013), Article ID 903904, 12 pageshttp://dx.doi.org/10.1155/2013/903904

Research Article

Generalized and Stability Rational Functions for Dynamic Systems of Reactor Kinetics

Ahmed E. Aboanber1,2

¹Department of Mathematics, Faculty of Science and Arts, Qassim University, P.O. Box 1300 Buraidah, Saudi Arabia²Department of Mathematics, Faculty of Science, Tanta University, Tanta 31527, Egypt

Received 27 February 2013; Revised 2 June 2013; Accepted 7 June 2013

Academic Editor: Arkady Serikov

Copyright © 2013 Ahmed E. Aboanber. 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.

The base of reactor kinetics dynamic systems is a set of coupled stiff ordinary differential equations known as the point reactor kinetics equations. These equations which express the time dependence of the neutron density and the decay of the delayed neutron precursors within a reactor are first order nonlinear and essentially describe the change in neutron density within the reactor due to a change in reactivity. Outstanding the particular structure of the point kinetic matrix, a semianalytical inversion is performed and generalized for each elementary step resulting eventually in substantial time saving. Also, the factorization techniques based on using temporarily the complex plane with the analytical inversion is applied. The theory is of general validity and involves no approximations. In addition, the stability of rational function approximations is discussed and applied to the solution of the point kinetics equations of nuclear reactor with different types of reactivity. From the results of various benchmark tests with different types of reactivity insertions, the developed generalized Padé approximation (GPA) method shows high accuracy, high efficiency, and stable character of the solution.

1. Introduction

The multigroup neutron diffusion equations, which are derived from the general mathematical form of neutron transport equation, are represented generally by a system of time- and space-dependent coupled partial differential equations for which approximate solutions are sought from neutron transport equation. The point reactor kinetics equations are a useful simplification of the time-dependant neutronic neutron diffusion equations for most reactors operating near delayed critical. Included in that system are equations which describe the neutron level, reactivity, an arbitrary number of delayed neutron groups, and any other variables that enter into the reactivity equation in the case of reactivity feedback.

One of the important properties in a nuclear reactor is the reactivity, due to the fact that it is directly related to the control of the reactor. The startup process of a nuclear reactor requires that reactivity is varied in the system by lifting the control rods discontinuously. In practice, the control rods are withdrawn at time intervals such that reactivity is introduced in the reactor core linearly, to allow criticality to be reached in a slow and safe manner. For safety analysis and transient behaviour of the reactor dynamics, the neutron density and the delayed neutron precursor concentration are important parameters to be studied. Under simplified conditions, an analytical solution to the point reactor kinetic equations can be obtained, for example, such as there being no extraneous neutron source, assuming a prompt jump approximation [1–6], but they are elusive when the coefficients vary with time. Moreover, a time-dependent reactivity inserted into a point reactor is coupled multiplicity with neutron density to form a set of linear equations with time-dependent coefficients. This time dependence makes it difficult to obtain analytical solution, consuming CPU time on Pc’s, and numerical integration is usually employed. Furthermore, the number of computations needed for each additional term in the explicit method for the exponential function makes the computing time prohibitive. For this and to avoid the instabilities associated with the computational effort involved in using the explicit methods, a particular class of approximations for the exponential function, namely, rational function approximations, is required. On the other hand as is well known, a long-standing problem in point reactor kinetics is the stiffness arising from the orders of magnitude difference between prompt and delayed neutron lifetimes, which results in the requirement for very small time step increments in numerical solutions to the kinetics equations.

There has been a great deal of research that focused on eliminating the problem of stiffness. Several methods especially adapted for solving the initial value problems for stiff systems of ordinary differential equations. Among the methods are numerical integration using Simpson’s rule [7], finite element method [8], Runge-Kutta procedures [9, 10], piecewise polynomial approach [11, 12], singular perturbation method [13–15], finite difference method [16], and other methods which are based on integral equation formulations with the slowly varying factor in each integrand represented by an assumed functional form [17–20]. Most of these methods are successful in some specific problems but still suffer, more or less, from disadvantages, for example, limitations on the maximum permissible time step for ensuring computational stability. For slow transients in fast systems, even the analytic continuation method [21] should be changed in its basic procedure to achieve efficient computation, and inability to handle point reactor kinetics equations in their full generality, and so forth.

For quantitative computations, a numerical calculation is used. Furthermore, a very small step size has to be adopted for stability due to the problem of stiffness in the differential equations when applying a general numerical method. Many numerical methods have been proposed in the literature, such as stiffness confinement method [22] and power series solution [23–26]. Several other creative schemes have been implemented as well as da Nóbrega [27], Lawrence and Doring [28], and Aboanber and Nahla [29] use a Padé approximation to the solution of the point kinetics equations. Also, the Padé approximation technique is applied to the nonlinear system including temperature feedback via analytical inversion method [30], CORE [31]. Recently, a very efficient method for obtaining a highly accurate solution of reactor kinetic equation was proposed in [32–34]. The basic idea of this approach is to combine classical discretization schemes with extrapolation techniques to reduce numerical error (e.g., discretization error and truncation error).

The present work is organized as follows. In Section 2 a newly developed general expression for a special type of functions, which allows us to approximate the exponential function in an economical manner, is derived. In addition, the factorization methods based on using temporarily the complex plane with the analytical inversion are generalized and applied to obtain the general solution of the point kinetics equations. In order to ensure the validity and stability of this class of approximations, the test is applied to the rational function approximations (Padé approximations). Numerical examples of applying the method to a variety of primary problems confirm that the time step increment size can be greatly increased and that much computing time (CPU) can be saved, as compared to other conventional methods, in Section 3. Concluding remarks are summarized in Section 4.

2. Mathematical Background and Point Kinetics Model

2.1. Basic Idea

In most situations encountered in the analysis of different types of reactors, it is sufficient to model the neutronic behaviour of the reactor by a low order approximation to the formally exact neutron transport equation. Practically, analytical solutions of the coupled time-dependent transport and precursor equations for reactor kinetics problems are stiff and difficult; so approximate methods are usually used. The most widely used of these approximations is the multigroup neutron diffusion theory. In many operations and accident transients, large and fast changes in neutron distribution can occur. Therefore, a large number of recalculations of the shape function are needed to obtain an acceptable accuracy of calculations. The basic problem which it attacks is that of finding solutions to the space-time neutron group diffusion equation derived from the neutron Boltzmann transport equations for a reactor. Every approximation and stability study presented in this work is based on the point reactor kinetics equations, whose derivation and particular analytical solutions may be found in the standard reactor physics literature. In particular, Duderstadt and Hamilton [35] provide a rather intuitive approach to the equations while Henry [36] shows how they can be rigorously obtained from the time-dependant diffusion equation.

The point reactor kinetics equations are a useful simplification of the time-dependant neutronic Boltzmann problem for most reactors operating near delayed critical. The variables and parameters that appear in these equations are the result of integrating the original partial derivatives formulation of the neutron conservation problem with respect to both the spatial and angular coordinates and to all available neutron energies. Nevertheless, these particular details are not of interest for the analysis that follows. The “global” dynamic behavior of the reactor can be described adequately using the point reactor model. For the G-delayed neutron precursor group, the point reactor model is written as follows: with where is the absolute reactivity defined as = , which is generally a function of time, is the reactor power (or neutron density), is the effective multiplication factor, is the neutron generation time, is the precursor concentrations of a group and is the average decay constant, is the effective delayed neutron fraction of a group, and is the source.

Equations (1) may be compactly written in vector form. If denotes a component column vector and denotes a matrix, then (1) are simply

A wealth of more or less economical methods have been devised to solve the system of (5). If the matrix elements , in particular the reactivity, are constants, then exact solution of (5) can be easily obtained, for example, Porsching [19]; Duderstadt and Hamilton [35] and Henry [36] provide a rather intuitive approach to the equations. However, if as is usually the case (and therefore ) is a function of time, then (5) no longer has an exact solution. The main problem that has been recognized for a long time consists in the so-called stiffness of the system of (5). Stiff differential equations frequently arise in physical situation characterized by the existence of greatly differing time constants and are common in particular to such other fields as common engineering applications, circuit analysis, missile guidance system, chemical engineering, and chemical analysis. For a review of the numerical implications of this problem, see, for example, Gear’s book [37].

Multiplying both sides of (5) by the integrating factor , we obtained the exact solution where the solution within each time step is sought starting from a piecewise constant approximation substituted in the interval , with a constant, namely, , which can be its value at some point in the interval (e.g., ) as Integrating both sides of the equation over the th time step from to gives Applying the usual constant parameters assumption, where and are replaced by their mean values and , on the interval under consideration with . Equation (8) reduces after subtracting and evaluating the integral to the equation of the form

The numerical scheme proposed in the present investigation is based on this equation. One class of numerical methods for solving the kinetics equations is based upon the Taylor series approximation to the matrix series of (9). The number of computations needed for each additional term makes the computing time prohibitive. For this and to avoid the instabilities associated with the computational effort involved in using the explicit methods, (9) requires a particular class of approximations for the exponential function, namely, Padé rational approximations.

2.2. Generalities of Analytical Inversion Method

Let us first state in this section the following important theorem before proceeding with the generalized rational function. The stated theorem for the matrices and , Abaonber and Nahla [30], provides us with a mathematical formula, which permits to approximate some special functions. Let the matrix have eigenvalues and eigenvectors   . Also, let have corresponding eigenvalues , similar to its eigenvalues which are the same as those of . Thus, the following theorem is valid for any for which is bounded for all . If and with that is normalized to unity, so the following expression holds for any where represents a special type of functions (e.g., exponential functions, Trigonometric functions, and Hyperbolic functions) where the eigenvalues of both the matrices and form a biorthonormal set when properly normalized.

Equation (10) is so far a mere mathematical manipulation. It has, however, a form that permits us to approximate some special functions (e.g., the exponential function, ) in an economical manner. The important note here is that if is a good approximation for , the th term in the summation should drop and consequently it will have a very small coefficient,

Since , thus, to a high degree of accuracy, we have where the sum is over only those for which (11) does not hold.

The vectors and are easily calculated from their defining equation: or in the normalized form: where is the normalization factor, satisfying the normalization condition , given by

In this connection the analytical inversion method is applied to the point kinetic matrix A. For numerical purpose, is then approximated by some simple function of , typically a real rational function of type : with some normalization (usually or ).

Whatever the degree, , of the denominator, it can always be factored into real linear and quadratic factors: For the sake of simplicity, suppose that , and is approximated by which can be obtained by solving recursively the matrix problems: with Due to the particular structure of , a semianalytical inversion (Analytical Inversion Methods (AIMs) [30]) can be performed for each elementary step resulting eventually in substantial time savings. They developed a general expression for the inverse of , where for a real and a point kinetic matrix of (4) with the reactivity appropriately averaged, the following expression is introduced: where, generally for delayed neutron group: Here, is any point kinetics-like matrix, that is, except when or or . Consequently, the previous expression is more general compared with that introduced by Aboanber and Nahla [30]. The in (18) being (generally) real numbers, the quadratic factors will have complex conjugate roots.

The expression (21) is of no great advantage by itself since we can solve directly the system of equations implied by the inverse shown with the same computational effort spending essentially in this case. However, the utility of the analytical inversion is evident when form a complex conjugate pair. In this connection the following pair of factors is considered, with , and which is a real matrix and thus has a real inverse. This expression can be expressed in the general form after carrying out some involved but straightforward algebra: where with When A is the point kinetics matrix, this general expression reduces, of course, to the formula previously derived by Aboanber and Nahla [30] for the particular case of the Padé 02 approximations where we have .

The remarkable feature of expressions (21) and (24) for the point kinetics matrix is that only the coefficient depends on the reactivity; unless the rational approximation chosen changes from one time step to the other, it is thus possible to predetermine one for all the vectors a and b and the matrix D, while only modifying at each time step, if at least the reactivity various with time. In general, the time dependence of the macroscopic cross-sections would also affect the entries of the point kinetics matrix, so that these terms should generally be regarded as time-dependent quantities; in this case, one should also reevaluate vector a, but this should probably be not done at each time step, as the is likely to change much slower than .

2.3. Generalized Padé Approximations

Numerical methods for differential equations when applied to a linear test problem generate a sequence of approximations which satisfy a difference equation. In the problem of finding an approximation to the solution of such differential equations by numerical methods, many investigators show that rational approximation to the exponential function is of great importance. Among the rational approximations to the exponential, the Padé approximations are the most common. The Padé approximation often gives better approximation of the function than truncating its Taylor series, and it may still work where the Taylor series does not converge. For these reasons Padé approximations are used extensively in computer calculations. Some famous one-step integration methods have since long been recognized as particular Padé approximations, like the forward Euler (FE), backward Euler (BE), and Crank-Nicholson (CN) schemes that are equivalent to the Padé 10, Padé 01 and Padé 11 approximations, respectively. Many classical integration schemes of the Runge-Kutta type have been shown, moreover, to lead to the Padé approximation when applied to the equation .

Given a rational function , the Padé approximation to is obtained if the remaining coefficients are chosen so that the Maclaurin expansion of agrees with the Maclaurin expansion of to terms. As such, the Padé approximations of appear as local approximations of at and are consequently rather poor approximations for large leading to restrictions on the time steps that can be taken.

In the present work, we follow the same ideas, but with different techniques; in particular, we derive and apply to the point kinetics equations global approximations of Padé in an attempt to increase the time steps to be used or, equivalently, to minimize the amount of work needed to achieve certain accuracy. For numerical purposes, is then approximated by some simple functions of , typically a real rational function of type given by (16), with usual renormalization ( or ).

In the case of exponential function of (16), the explicit expressions for the ’s and ’s are obtained. Dealing with the following representative equations for the exponential functions: the following system of algebraic equations of the unknown coefficients ’s and ’s are obtained [29]: The utility of the previous system of algebraic equations is evident, especially, when one tries to determine any type of Padé approximations of any order. Actually, the Padé approximation to is a one-point Padé approximation this point being , and we have .

2.4. Stability Considerations

The stability property is strongly problem dependent. In order to get rid of this constraint, we pay more attention to convergence and stability problems, by introducing a test problem and , where is a complex constant. Convergence guarantees us that we can approach the exact solution as closely as we want by diminishing correspondingly the time step, and the order of convergence gives us an idea of how fast this can be achieved when is reduced; these notions, especially the last, turn out to be of little practical interest when we are interested in approximating correctly the exact solution in the most economical way, that is, in a minimum number of time steps. This is true also when stability is considered, following the concept of time A-stability, a method is said to be A-stable if all numerical approximations tends to zero as the number of time steps tends to infinity when it is applied to the test differential equation with a fixed positive and any complex constant , with [37]. This concept of stability is clearly much stronger than the usual (traditional) ones and must usually replace them when dealing with stiff systems of ODEs.

Applying the rational approximation methods developed in the previous section for different type of Padé approximations to the differential equation , we obtain where and A-stability requires where is any of the rational approximations considered. If , the exact solution itself is growing exponentially and accuracy becomes the only requirement. A given integration method, (29), will be A-stable if it satisfies (30). This, in turn, will be true if and only if as a consequence of the maximum modulus theorem. Varga [38] has shown that all the Padé approximations with are absolutely stable on the real negative axis: These approximations are thus possible if one is sure that the spectrum of A is real. Truly, this is not the case for the space-time kinetics described in our case. A method is A()-stable if its stability domain is a wedge in the left half of the complex plane where the angle between the edge and the real axis is equal to . Consequently, it is preferable to use an A()-stable scheme for which a rational function is such that where then the method is A()-stable (or acceptable), where . The angle is only one of a number of parameters which have been proposed for measuring the extent of the stability region, but it is probably the best such measure. Computer calculations were performed to estimate the value of for different cases of Padé approximations as shown in Table 1. The method is A(0)-acceptable if there exists an such that is A()-acceptable and -acceptable (A-stable) if and only if is A()-acceptable for all [39].

Table 1: Coefficients and explicit Padé approximations with their stability for the exponential function. .

Approximation
Stability

0

1
A-Stable

0

−1

0

0

0

0
A-Stable

0

0

0

1
A(59.62)-Stable

0

0

0

1
A-Stable

0

−1

0

0

1
A-Stable

0

0

1
A(82.81)-Stable

0

0

1
A-Stable

The scheme of rational approximations described in Table 1 showed that the Padé , , and approximations lead to A-stable schemes for any , while this is not any more true for Padé (, ) approximation. For example, the Padé (0, 3) is A()-stable where . As will be shown in the next section, many approximations which are proved to be accurate in most applications are not A-stable.

3. Numerical Results for Selected Benchmark Test Cases

In order to verify the developed method proposed in this study for the prediction of the neutron density distribution, the Padé approximation method was used for the numerical solution of the point kinetics equations. The nuclear parameters used in the validation of formulation proposed in this study which are to be used for the proposed benchmarks against the formulation proposed by [32–34] and other reference literatures are reported in Table 2. Several authors claim to present benchmark solutions where they only report 2 to 3 accurate digits. Several reactor configurations are considered, including thermal and fast reactors and multiprecursor materials with varying delayed emission fractions. All benchmark problems correspond to transient cases starting from a critical state, with equilibrium between neutron and precursor concentration. The initial conditions and should be given to complete the formulation of the problem. In the absence of an external neutron source, the equilibrium condition for a steady-state density is given by . In the following, the accuracy of the developed Padé method is benchmarked against the extremely accurate EPCA (Ganapol, 2013) and CATS methods [40], which can produce 9-digit or more reference values. The comparison with other results are also reported to highlight the performance of developed Padé technique.

Table 2: Delayed-neutron parameters of typical reactors used in the benchmark test cases.

Fast reactor
Thermal reactor (I)
Thermal reactor (II)
Thermal reactor (III)

1
0.0129
16.720
0.0127
28.500
0.0127
26.600
0.0124
21.500

2
0.0311
123.20
0.0317
159.75
0.0317
149.10
0.0305
142.40

3
0.1340
95.040
0.1150
141.00
0.1150
131.60
0.1110
127.40

4
0.3310
144.30
0.3110
305.25
0.3110
284.90
0.3010
256.80

5
1.2600
45.340
1.4000
96.000
1.4000
89.600
1.1300
74.800

6
3.2100
15.400
3.8700
19.500
3.8700
18.200
3.0000
27.300

(pcm)  sec
(pcm)  sec
(pcm)  sec
(pcm) sec

Benchmark Series 1: Step Reactivity Insertion. In order to compare the different methods a series of benchmarks is considered. The first example simulates a step-reactivity insertion without feedback. The Padé approximation method results are compared to well-established benchmarks. The first benchmark considered is for a step reactivity insertion of $−0.5 and $1 into a thermal reactor. The neutron generation time is 5 × 10⁻⁴ s and the kinetics parameters are given in Table 2. Table 3 presents the results corresponding to this case for both negative and positive reactivity insertion in a thermal reactor configuration, thermal I. The average step size used in the developed FORTRAN code based on the present generalized Padé method (GPA) is 0.1 s. The converged densities agree to all of the digits of the recent reference solution EPCA and CATS. The approximately zero relative percentage error between GPA, EPCA, and GPA indicates that the relative error of the neutron density is much lower than the predefined error tolerance with the advantage of small time step size and of course CPU time. This simple preliminary test shows that the presented methods are consistent with the both EPCA and CATS solutions, when rounded, as well as with the analytical solution to the 9th digit.

Table 3: Neutron density of benchmark 1: step of $−0.5 and $1 into a thermal reactor.

Time (sec)
Analytical (A)
CATS
Padé (P)
Percentage error (A-P)/A

0.1
6.9892522557 (−1)
6.989252256 (−1)
6.989252557 (−1)
−8.375 (−6)

1.0
6.0705356561 (−1)
6.070535656 (−1)
6.070535946 (−1)
−2.744 (−6)

5.0
4.8255302992 (−1)
4.825530299 (−1)
4.825530607 (−1)
−5.317 (−6)

10.0
3.9607769071 (−1)
3.960776907 (−1)
3.960777198 (−1)
−6,890 (−6)

100
7.1582854438 (−2)
7.158285444 (−2)
7.158286234 (−2)
−1.194 (−5)

200
1.6934791658 (−2)
1.693479166 (−2)
1.693479374 (−2)
−1.276 (−5)

0.1
2.5157661414 (00)
2.515766141 (00)
2.515766093 (00)
2.003 (−6)

0.5
1.0362533810 (01)
1.036253381 (01)
1.036253327 (01)
2.921 (−6)

1.0
3.2183540945 (01)
3.218354095 (01)
3.218353783 (01)
7.816 (−6)

10.0
3.2469788980 (09)
3.246978898 (09)
3.246975705 (09)
9.581 (−5)

100
2.5964846465 (89)
2.596484647 (89)
2.596458993 (89)
9.880 (−4)

Benchmark Series 2: Ramp Reactivity Insertion. The second type of reactivity triggering a transient is a ramp insertion. The transient is studied here in both a thermal (thermal II) and a fast reactor (fast I) configurations, for two values of the ramp. The reactivity driving function is for thermal reactor and for fast reactor. The neutron density at five time points between 0 and 10 s is computed and shown in Table 4. The reference value was provided by Hermite polynomial method with the fixed step size of 0.0001 s, Li et al. [20]. An average step size of 0.001 s was obtained with the developed (GPA) method. By examining the relative error given in Table 4, the developed method gives accurate values consistently. All relative errors are much lower than the predefined error tolerance.

Table 4: Neutron density of benchmark 2: Ramp reactivity insertion for thermal reactor and for fast reactor.

(a)  

Time (sec)
MAPLE
CORE (C)
PCA
EPCA
Padé (P)
Percentage errors (C-P)/C

2
1.3382
1.3382 (00)
1.338200434 (00)
1.338200050 (00)
1.338196072 (00)
2.935 (−4)

4
2.2284
2.2285 (00)
2.228442526 (00)
2.228441897 (00)
2.228435296 (00)
2.903 (−3)

6
5.5820
5.5822 (00)
5.582053977 (00)
5.582052450 (00)
5.582036367 (00)
2.931 (−3)

8
4.2783 (01)
4.2790 (01)
4.278630638 (01)
4.278629574 (01)
4.278618740 (01)
9.910 (−3)

10
4.5116 (05)
4.5391 (05)
4.511636225 (05)
4.511636240 (05)
4.511647526 (05)
−6.048 (−1)

(b)  

Time (sec)
CATS
PCA
EPCA (E)
Padé (P)
Percentage errors (E-P)/E

0.01
1.010097111 (00)
1.010098578 (00)
1.010097112 (00)
1.009610584 (00)
−4.817 (−2)

0.1
1.113320112 (00)
1.113321641 (00)
1.113320113 (00)
1.112733108 (00)
−5.273 (−2)

0.2
1.260559925 (00)
1.260561549 (00)
1.260559925 (00)
1.259817235 (00)
−5.892 (−2)

0.5
2.136409107 (00)
2.136411382 (00)
2.136409108 (00)
2.134468802 (00)
−9.082 (−2)

1.0
1.207814197 (03)
1.207814360 (03)
1.207814198 (03)
1.203402181 (03)
−3.653 (−1)

1.1
3.257593355 (99)
3.257552723 (99)
3.257593356 (99)
3.196193006 (99)
1.882 (00)

1.15
1.028975360 (219)
1.028957612 (219)
1.028975360 (219)
1.009691369 (219)
1.874 (00)

A fast ramp of $1/s produces the neutron density given in Table 4. The stiffness is relatively strong. The neutron density and delayed neutron precursor density exponentially increase after 0.5 sec. The relatively small average step size of 0.001 sec of this case also shows that the step size can be automatically determined based on the stiffness of the problem. The relative error that is lower than the predefined tolerance in Table 4 also proves that the GPA method gives the exact results as the mathematical expectation.

For the ramp variation of reactivity as mentioned by Ganapol et al., 2013, even in the case of linear kinetics the PCA is inaccurate due to the variation of the reactivity within the discretization step, the problem which is not included in the GPA method. Indeed the numerical error due to the constant approximation can be reduced by choosing a time step sufficiently short as compared with reactivity variation and neutron lifetime. In the EPCA, this problem is reduced through coupling with the submesh that ensures an update of the reactivity within each time step. Results in Table 4 compare PCA and EPCA as well as CATS with the GPA benchmark solutions. While PCA gets the first four to at most seven digits, GPA is exact to at worst the 5th digit, which agrees very well with EPCA. Also the reference values reported in Table 4, showing that its accuracy is at best 4 digits of the true benchmark result. In addition, the GPA average step size of 0.01 sec is smaller than the GRK average step size of 0.02 sec [10]. By comparing the GPA CPU time with PWS CPU time of 3.02 sec [24], the GPA method is about a hundred times faster than the PWS method.

Benchmark Series 3: Zigzag Reactivity Ramp. A zigzag ramp reactivity is inserted into a thermal reactor taken from [27], where the GPA method is applied to compute the response of a reactor core to a zigzag ramp input of reactivity. Consider The physical parameters in this example described in Table 2, Thermal Reactor I. The results obtained with the adopted GPA method are compared against those obtained by PCA and EPCA methods as well as with Chao and Attard [22] method (SCM). These computations are done with a fixed time step size of 0.01 s. According to all zero relative errors in Table 5, the GPA method provides accurate results. The errors associated with the four methods show the consistency of the GPA method. The neutron density profiles versus time are shown in Figure 1. The zigzag reactivity ramp brings fluctuation to the neutron density profile whose curve shape follows the reactivity driving function, but the relative magnitude of the response of the delayed neutron precursor density to the zigzag reactivity ramp is relatively small, Figure 2.

Table 5: Neutron density of benchmark 3: Zigzag reactivity rampfor thermal reactor.

Time (sec)
CATS
PCA
EPCA
da Nóbrega
Padé

0.5
1.721422422 (00)
1.721116849 (00)
1.721422424 (00)
1.721422
1.721580296 (00)

1
1.211127415 (00)
1.210794797 (00)
1.211127414 (00)
1.211127
1.211290195 (00)

1.5
1.892226140 (00)
1.891835965 (00)
1.892226142 (00)
1.892226
1.892442369 (00)

2.0
2.521600530 (00)
2.521590673 (00)
2.521600530 (00)
2.521601
2.521977814 (00)

10
1.204710535 (01)
1.204687181 (01)
1.204710536 (01)
12.04711
1.205221563 (01)

100
6.815556888 (07)
6.815428434 (07)
6.815556889 (07)
—
6.839867192 (07)

Figure 1: Variations of neutron density for zigzag reactivity at different initial values.

Figure 2: Variations of the precursors density for zigzag reactivity for six groups of delayed neutrons.

Benchmark Series 4: Sinusoidal Reactivity. For further investigation of the developed method, more calculations were made. A sinusoidal reactivity insertion transient is simulated in Figure 3 for a thermal configuration (Thermal III). In this case the reactivity has the following form: , where the notation has its conventional meaning. The equations of reactivity combined with the system of point kinetics will give a crude simulation of a slow self-limiting excursion, where is the constant part of the excess reactivity, and is a given function characterizing the time dependence of the reactivity. The parameter is a positive number that represents the magnitude of the variable part of the excess reactivity in dollars. This parameter is assumed to be sufficiently small compared to unity. That is a real assumption since, except in the accidental case, the excess reactivity is always less than one dollar. The wave form oscillations of the flux is such that a modified sinusoidal curve is obtained for different values of the parameter in Figure 3 with reactivity oscillation versus time.

Figure 3: Variations of neutron density for oscillatory reactivity at different initial values.

Benchmark Series 5: Feedback Reactivity. The presence of temperature feedback is useful in providing an estimate of the transient behaviour of a reactor power. The previous benchmark series of reactivity inputs discussed the neutron kinetics (i.e., the response of neutron density in a nuclear reactor to an external reactivity input) under the implicit assumption that the level of the neutron density does not affect the properties of the system that determine the neutron kinetics, most notably the reactivity. This is the situation when the neutron density is sufficiently small that the fission heat does not affect the temperature of the system (i.e., zero power). However, in an operating nuclear reactor the neutron density is large enough that any change in fission heating resulting from change in neutron density will produce changes in temperature, which in turn will produce changes in reactivity, or in reactivity feedback.

It is assumed that the reactor has a negative temperature coefficient of reactivity () when a small step reactivity () is inserted. When we consider the temperature feedback, the real reactor reactivity is given by where is the temperature increment of the reactor. After the reactivity is inserted into the reactor, the power responds quickly and the adiabatic model can be used for the calculation of reactor temperature [41, 42]. Then the derivative of temperature to time can be given as follows: where is the reciprocal of thermal capacity of reactor. Combining (36) and (37) results in Generally, where represented the impressed reactivity variation (generally polynomial in ), and is the “shutdown coefficient” of the reactor system. It is usually taken as a negative constant of magnitude ranging from 10⁻¹³ cm³/sec for slow system to 10⁻⁷ cm³/sec for fast metal systems. The results of comparison for this case are shown in Table 6.

Table 6: Neutron density of benchmark 5: compensated response to ramp function reactivity changes with prompt neutron generation time, Λ = 5 * 10⁻⁴ sec, and shutdown coefficients are in the range of 10⁻¹¹ cm³/sec to 10⁻¹³ cm³/sec.

Method
First peak
Time to the first peak
First asymptotic neutron density

0.1
1 × 10⁻¹¹
PWS
2.4238380 × 10¹¹
0.224000
1.0066510 × 10¹⁰

GRK4-A
2.4649574 × 10¹¹
0.22626
1.0066499 × 10¹⁰

GRK4-T
2.4652371 × 10¹¹
0.22626
1.0066500 × 10¹⁰

GRK4-H
2.4654628 × 10¹¹
0.22626
1.0066501 × 10¹⁰

Padé 11
2.4238380 × 10¹¹
0.22400
1.0060053 × 10¹⁰

1 × 10⁻¹³
PWS
2.8666680 × 10¹³
0.23890
1.0066520 × 10¹²

GRK4-A
2.9580990 × 10¹³
0.24066
1.0066480 × 10¹²

GRK4-T
2.9580157 × 10¹³
0.24066
1.0066480 × 10¹²

GRK4-H
2.9577276 × 10¹³
0.24066
1.0066482 × 10¹²

Padé 11
2.8929486 × 10¹³
0.23800
1.0066439 × 10¹²

0.01
1 × 10⁻¹¹
PWS
2.0064900 × 10¹⁰
1.10617
1.0263960 × 10⁹

GRK4-A
2.0232131 × 10¹⁰
1.10736
1.0246410 × 10⁹

GRK4-T
2.0232020 × 10¹⁰
1.10736
1.0246424 × 10⁹

GRK4-H
2.0232065 × 10¹⁰
1.10736
1.0246394 × 10⁹

Padé 11
2.0120302 × 10¹⁰
1.10600
1.0255595 × 10⁹

1 × 10⁻¹³
PWS
2.4831920 × 10¹²
1.15528
9.6705500 × 10¹⁰

GRK4-A
2.5061997 × 10¹²
1.15668
9.6374374 × 10¹⁰

GRK4-T
2.5061921 × 10¹²
1.15668
9.6374817 × 10¹⁰

GRK4-H
2.5061852 × 10¹²
1.15668
9.6374301 × 10¹⁰

Padé 11
2.4909973 × 10¹²
1.15500
9.6629450 × 10¹⁰

0.001
1 × 10⁻¹¹
PWS
1.2652650 × 10⁹
7.48864
1.9212660 × 10⁸

GRK4-A
1.2760915 × 10⁹
7.48728
1.9207147 × 10⁸

GRK4-T
1.2761036 × 10⁹
7.48728
1.9207122 × 10⁸

GRK4-H
1.2760943 × 10⁹
7.48728
1.9206984 × 10⁸

Padé 11
1.2740050 × 10⁹
7.48864
1.9282197 × 10⁸

1 × 10⁻¹³
PWS
1.7068560 × 10¹¹
7.68362
2.0046470 × 10¹⁰

GRK4-A
1.7245128 × 10¹¹
7.68204
2.0040163 × 10¹⁰

GRK4-T
1.7245159 × 10¹¹
7.68204
2.0040077 × 10¹⁰

GRK4-H
1.7245002 × 10¹¹
7.68204
2.0039905 × 10¹⁰

Padé 11
1.7209210 × 10¹¹
7.68263
2.0067995 × 10¹⁰

4. Conclusions

The point reactor kinetics equations are a useful simplification of the time-dependant neutronic Boltzmann problem that allow studying the stability of a reactor core using dynamical systems theory. From a strictly mathematical point of view, the point kinetics of reactor dynamic systems are of great importance for understanding the time-dependent behaviour of the neutron density. In nuclear reactor the response to either a planned or unplanned change in the reactor conditions is of great importance to the safe and reliable operation of the reactor. It is therefore important to understand the response of the neutron density and how it relates to the speed of lifting control rods (reactivity variations). The generalized Padé rational approximation (GPA) is used to approximate the exponential function in the point kinetics equation solution. General expression for the coefficients is developed for any type of approximations. In addition, the roots of the inhour formula are used as the eigenvalues of the point kinetics matrix. A semi-analytical inversion and factorization methods are developed and generalized for each elementary step resulting eventually in substantial time saving. It has been shown that GPA can exhibit a significant computational advantage by reducing the CPU time over more conventional approximations, and that it can be applied to a wide range of important problems. By comparing the GPA CPU time with PWS CPU time of 3.02 sec, the GPA method is about a hundred times faster than the PWS method as mentioned in benchmark. If larger time steps are considered, as reported in the presented tables, the GPA would still be quite useful over fairly time steps, especially at the extreme values of the inserted reactivity. The methods presented here should, however, retain much of their interest when applied to all cases where, after a small step-like transient, the reactivity reaches new stationary values. The results confirm the theoretical analysis and indicate the range of applicability of the presented methods. Their great advantage is that all we have said earlier remains valid for full space-time kinetics. In addition, the numerical tests demonstrated that high accuracy could be achieved in strong step reactivity insertion. The low step size was required to compute the ramp reactivity to acquire the same accuracy.

Acknowledgment

This work is supported by the Scientific Research Deanship, Qassim University Kingdom of Saudi Arabia, Research Centre, Faculty of Science and Arts.

References

A. E. Aboanber, “Analytical solution of the point kinetics equations by exponential mode analysis,” Progress in Nuclear Energy, vol. 42, no. 2, pp. 179–197, 2003.

A. E. Aboanber, “An efficient analytical form for the period-reactivity relation of beryllium and heavy-water moderated reactors,” Nuclear Engineering and Design, vol. 224, no. 3, pp. 279–292, 2003.

F. Zhang, W.-Z. Chen, and X.-W. Gui, “Analytic method study of point-reactor kinetic equation when cold start-up,” Annals of Nuclear Energy, vol. 35, no. 4, pp. 746–749, 2008.

W. Z. Chen, B. Kuang, L. F. Guo, Z. Y. Chen, and B. Zhu, “New analysis of prompt supercritical process with temperature feedback,” Nuclear Engineering and Design, vol. 236, no. 12, pp. 1326–1329, 2006.

W. Z. Chen, L. F. Guo, B. Zhu, and H. F. Li, “Accuracy of analytical methods for obtaining supercritical transients with temperature feedback,” Progress in Nuclear Energy, vol. 49, no. 4, pp. 290–302, 2007.

H. F. Li, W. Z. Chen, F. Zhang, and L. Luo, “Approximate solutions of point kinetics equations with one delayed neutron group and temperature feedback during delayed supercritical process,” Annals of Nuclear Energy, vol. 34, no. 6, pp. 521–526, 2007.

G. R. Keepin and C. W. Cox, “General solution of the reactor kinetic equations,” Nuclear Science and Engineering, vol. 8, pp. 670–690, 1960.

C. M. Kang and K. F. Hansen, “Finite element methods for reactor analysis,” Nuclear Science and Engineering, vol. 51, no. 4, pp. 456–495, 1973.

J. C. Allerd and D. S. Carter, “Kinetics of homogeneous power reactors of the LAPRE Type,” Nuclear Science and Engineering, vol. 3, pp. 482–503, 1958.

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.

C. M. Kang, “Piecewise polynomial approximations for the point kinetics equations,” Transactions of the American Nuclear Society, vol. 14, pp. 201–202, 1971.

J. P. Hennart, “Piecewise polynomial approximations for nuclear reactor point and space kinetics,” Nuclear Science and Engineering, vol. 64, no. 4, pp. 875–901, 1977.

W. L. Hendry and G. I. Bell, “An analysis of the time-dependent neutron transport equation with delayed neutrons by the method of matched asymptotic expansions,” Nuclear Science and Engineering, vol. 35, pp. 240–248, 1969.

R. Goldstein and L. M. Shotkin, “Use of the prompt-jump approximation in fast reactor kinetics,” Nuclear Science and Engineering, vol. 38, no. 2, pp. 94–103, 1969.

T. Blenski, A. Gadomski, and J. Mika, “Higher order prompt-jump approximation in reactor kinetics,” Nuclear Science and Engineering, vol. 66, no. 3, pp. 277–283, 1978.

H. D. Brown, “A general treatment of flux transients,” Nuclear Science and Engineering, vol. 2, pp. 687–693, 1957.

F. T. Adler, “Reactor kinetics: integral equation formulation,” Journal of Nuclear Energy Parts A/B, vol. 15, no. 2-3, pp. 81–85, 1961.

K. F. Hansen, B. V. Koen, and W. W. Little, “Stable numerical solutions of the reactor kinetics equations,” Nuclear Science and Engineering, vol. 22, pp. 51–59, 1965.

T. A. Porsching, “Numerical solution of the reactor kinetics equations by approximate exponentials,” Nuclear Science and Engineering, vol. 25, pp. 183–188, 1966.

H. Li, W. Chen, L. Luo, and Q. Zhu, “A new integral method for solving the point reactor neutron kinetics equations,” Annals of Nuclear Energy, vol. 36, no. 4, pp. 427–432, 2009.

J. C. Vigil, “Solution of the reactor kinetics equations by analytic continuation,” Nuclear Science and Engineering, vol. 29, pp. 392–401, 1967.

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

J. Basken and J. D. Lewins, “Power series solutions of the reactor kinetics equations,” Nuclear Science and Engineering, vol. 122, no. 3, pp. 407–416, 1996.

A. E. Aboanber and Y. M. Hamada, “PWS: an efficient code system for solving space-independent nuclear reactor dynamics,” Annals of Nuclear Energy, vol. 29, no. 18, pp. 2159–2172, 2002.

A. E. Aboanber and Y. M. Hamada, “Power series solution (PWS) of nuclear reactor dynamics with newtonian temperature feedback,” Annals of Nuclear Energy, vol. 30, no. 10, pp. 1111–1122, 2003.

T. Sathiyasheela, “Power series solution method for solving point kinetics equations with lumped model temperature and feedback,” Annals of Nuclear Energy, vol. 36, no. 2, pp. 246–250, 2009.

J. A. W. da Nóbrega, “A new solution of the point kinetics equations,” Nuclear Science and Engineering, vol. 46, pp. 366–370, 1971.

R. D. Lawrence and J. J. Doring, “A smoothing and extrapolation method for point kinetics,” Transactions of the American Nuclear Society, vol. 24, pp. 199–201, 1976.

A. E. Aboanber and A. A. Nahla, “On pade' approximations to the exponential function and application to the point kinetics equations,” Progress in Nuclear Energy, vol. 44, no. 4, pp. 347–368, 2004.

A. E. Aboanber and A. A. Nahla, “Generalization of the analytical inversion method for the solution of the point kinetics equations,” Journal of Physics A, vol. 35, no. 14, pp. 3245–3263, 2002.

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.

B. D. Ganapol, “The refined way to solve the reactor point kinetics equations for imposed reactivity insertions,” Nuclear Technology and Radiation Protection, vol. 24, no. 3, pp. 157–166, 2009.

B. D. Ganapol and P. Picca, “A highly accurate benchmark for reactor point kinetics with feedback,” in Proceedings of the 17th Pacific Basin Nuclear Conference, Cancún, México, October 2010.

B. D. Ganapol, “A highly accurate algorithm for the solution of the point kinetics equations,” Annals of Nuclear Energy, 2013.

J. J. Duderstadt and L. J. Hamilton, Nuclear Reactor Analysis, John Wiley & Sons, New York, NY, USA, 1976.

A. F. Henry, Nuclear Reactor Analysis, The MIT Press, Cambridge, Mass, USA, 1975.

C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, Englewood Cliffs, NJ, USA, 1971.

R. S. Varga, “On higher order stable implicit methods for solving parabolic partial differential equations,” Journal of Mathematical Physics, vol. 40, pp. 220–231, 1961.

S. P. Nørsett, “C-Polynomials for rational approximation to the exponential function,” Numerische Mathematik, vol. 25, no. 1, pp. 39–56, 1975.

B. D. Ganapol, P. Picca, A. Previti, and D. Mostacci, The Solution of the Point Kinetics Equations Via Acceleration Taylor Series (CATS), Knoxville, Tenn, USA, 2012.

D. L. Hetrick, Dynamics of Nuclear Reactor, American Nuclear Society, La Grange Park, Ill, USA, 1993.

A. E. Aboanber, “Stability of generalized Runge-Kutta methods for stiff kinetics coupled differential equations,” Journal of Physics A, vol. 39, no. 8, pp. 1859–1876, 2006.