Generalized and Stability Rational Functions for Dynamic Systems of Reactor Kinetics

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: 𝑑𝑛(𝑡)(𝑡)𝑛(𝑡)+𝑑𝑡=𝐴^𝐺_𝑖=1𝜆_𝑖𝐶_𝑖(𝑡)+𝑓(𝑡),𝑑𝐶_𝑖(𝑡)𝑑𝑡=𝐵_𝑖(𝑡)𝑛(𝑡)−𝜆_𝑖𝐶_𝑖(𝑡), (1) with 𝐴(𝑡)=𝜌(𝑡)−𝛽(𝑡)Λ(𝑡),𝐵_𝑖(𝑡)=𝛽_𝑖Λ(𝑡)for𝑖=1,2,…,𝐺, (2) where 𝜌(𝑡) is the absolute reactivity defined as = (𝑘−1)/𝑘=𝑘_ex/𝑘, 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 Ψ(𝑡)=𝑛(𝑡)𝐶₁(𝑡)⋮𝐶_𝐺(𝑡)=Ψ₁(𝑡)Ψ₂(𝑡)⋮Ψ_𝐺+1(𝑡)(𝑡)=𝑓(𝑡)0⋮0,𝑆 (3) denotes a 𝐺+1 component column vector and 𝐴(𝑡)=𝜌(𝑡)−𝛽Λ𝜆₁⋯𝜆_𝐺𝛽₁Λ−𝜆₁𝑂𝛽⋮⋱_𝐺Λ𝑂−𝜆_𝐺=𝑎_𝑖𝑗 (4) denotes a (𝐺+1)×(𝐺+1) matrix, then (1) are simply Ψ^•(𝑡)=𝐴(𝑡)⋅Ψ(𝑡)+𝑆(𝑡). (5)

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 𝑑𝑒𝑑𝑡^{−𝑡𝐀(𝑡)}Ψ=𝑒^{−𝑡𝐀(𝑡)}𝑆_𝑖, (6) where the solution within each time step is sought starting from a piecewise constant approximation substituted in the interval [𝑡_𝑘,𝑡_𝑘+1], with a constant, namely, 𝐴_𝑘, which can be its value at some point in the interval (e.g., 𝐴_𝑘≈𝐴(𝑡_𝑘)𝑜𝑟𝐴_𝑘≈1/2(𝐴(𝑡_𝑘)+𝐴(𝑡_𝑘+1))) as 𝜌(𝑡)𝑡≈𝜌^𝑘+𝑡_𝑘+12=𝜌_𝑘,𝑆(𝑡)𝑡≈𝑆^𝑘+𝑡_𝑘+12=𝑆_𝑘for𝑡_𝑘≤𝑡≤𝑡_𝑘+1. (7) Integrating both sides of the equation over the 𝑖th time step from 𝑡_𝑖 to 𝑡_𝑖+1 gives 𝑒^{−𝑡𝐀(𝑡)}Ψ_𝑖+1−𝑒^{−𝑡𝐀(𝑡)}Ψ_𝑖=^𝑡𝑘+1_𝑡𝑘𝑒^{−𝑡𝐀(𝑡)𝑢}𝑆_𝑖𝑑𝑢. (8) Applying the usual constant parameters assumption, where 𝐴(𝑡) and 𝑆(𝑡) are replaced by their mean values 𝐴 and 𝑆, on the interval under consideration with ℎ=𝑡_𝑘+1−𝑡_𝑘. Equation (8) reduces after subtracting and evaluating the integral to the equation of the form Ψ_𝑖+1=𝑒^𝐴ℎΨ_𝑖+𝑒^𝐴ℎ−𝐼𝐴⁻¹𝑆_𝑖,whereℎ=𝑡_𝑖+1−𝑡_𝑖. (9)

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 𝑈_𝑖  (𝑖=1,2,…,ℓ). 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 𝑈^𝑇_𝑖𝑉_𝑖=1 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) 𝐅(𝐀)=𝑓(𝐀)+^ℓ_𝑖=1𝐹𝜔_𝑖−𝑓𝜔_𝑖𝑈_𝑖𝑉^𝑇_𝑖, (10) 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, exp(ℎ𝐴)) 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, 𝐹𝜔_𝑖−𝑓𝜔_𝑖⪻1. (11)

Since exp(ℎ𝜔_𝑖)≈𝑓(ℎ𝜔_𝑖), thus, to a high degree of accuracy, we have ()()()+expℎ𝐀≅𝑔ℎ𝐀=𝑓ℎ𝐀^/_𝐤expℎ𝜔_𝑘−𝑓ℎ𝜔_𝑘𝑈_𝑘𝑉^𝑇_𝑘, (12) where the sum ∑^/_𝑘 is over only those 𝑘 for which (11) does not hold.

The vectors 𝐔_𝑘 and 𝐕_𝑘 are easily calculated from their defining equation: {𝜔𝐈−𝐀}𝐔_𝑘=0,so𝐔_𝑘1𝜇=col¹𝜆₁+𝜔_𝑘⋯𝜇_𝐺𝜆_𝐺+𝜔_𝑘,𝜔𝐈−𝐀^𝑇𝑉_𝑘=0,so𝐕_𝑘1𝜆=col¹𝜆₁+𝜔_𝑘⋯𝜆_𝐺𝜆_𝐺+𝜔_𝑘, (13) or in the normalized form: 𝐕_𝑘=𝜈_𝑘1𝜆col¹𝜆₁+𝜔_𝑘⋯𝜆_𝐺𝜆_𝐺+𝜔_𝑘, (14) where 𝜈_𝑘 is the normalization factor, satisfying the normalization condition 𝑈^𝑇_𝑖𝑉_𝑖=1, given by 𝜈_𝑘=1+^𝐺_𝑖=1𝜇_𝑖𝜆_𝑖(𝜆_𝑖+𝜔_𝑘)²⁻¹≺1. (15)

In this connection the analytical inversion method is applied to the point kinetic matrix A. For numerical purpose, exp(ℎ𝐴) is then approximated by some simple function of (ℎ𝐴), typically a real rational function of type (𝑝,𝑞): 𝑅_𝑝𝑞ℎ𝐴=𝑃_𝑝ℎ𝐴𝑄_𝑞ℎ𝐴=^𝑝_𝑘=0𝑎_𝑘ℎ𝐴^𝑘𝑞_ℓ=0𝑏_ℓℎ𝐴^ℓ−1, (16) with some normalization (usually 𝑎₀ or 𝑏₀=1).

Whatever the degree, 𝑞, of the denominator, it can always be factored into real linear and quadratic factors: ^𝑞_ℓ=0𝑏_ℓ𝑥^ℓ=^𝑠1_𝑖=1𝐿_𝑖(𝑥)^𝑠2_𝑗=1𝑞_𝑖(𝑥),𝑠₁+2𝑠₂=𝑞. (17) For the sake of simplicity, suppose that 𝑆=0, and Ψ_𝑖+1 is approximated by Ψ_𝑖+1=^𝑝_𝑘=0𝑎_𝑘ℎ𝐴^𝑘𝑞_ℓ=0𝑏_ℓℎ𝐴^ℓ−1Ψ_𝑖 (18) which can be obtained by solving recursively the matrix problems: 𝑞_𝑗ℎ𝐴Φ_𝑗=Φ_𝑗−1,𝑗=1,…,𝑠₂,𝐿_𝑖ℎ𝐴Φ_𝑠2+𝑖=Φ_{𝑠2+𝑖−1},𝑖=1,…,𝑠₂ (19) with Φ_𝑖=^𝑝_𝑘=0𝑎_𝑘ℎ𝐴^𝑘Ψ_𝑖,Ψ_𝑖+1=Φ_𝑠2+𝑠1. (20) 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: [𝐈−𝜀𝐀]⁻¹=𝛾𝐚𝐛^𝑇+𝐃, (21) where, generally for 𝐺 delayed neutron group: 𝛾=1−𝜀𝑎₁₁+^𝐺+1_𝑖=2𝜀²𝑎_1𝑖𝑎_𝑖11−𝜀𝑎_𝑖𝑖⁻¹,𝐚=1,𝜀𝑎₂₁1−𝜀𝑎₂₂,𝜀𝑎₃₁1−𝜀𝑎₃₃,…,𝜀𝑎_𝐺+1,11−𝜀𝑎_{𝐺+1,𝐺+1}^𝑇,𝐛=1,𝜀𝑎₁₂1−𝜀𝑎₂₂,𝜀𝑎₁₃1−𝜀𝑎₃₃,…,𝜀𝑎_1,𝐺+11−𝜀𝑎_{𝐺+1,𝐺+1}^𝑇,1𝐃=diag0,1−𝜀𝑎₂₂,11−𝜀𝑎₃₃1,…,1−𝜀𝑎_{𝐺+1,𝐺+1}, (22) Here, 𝐀=[𝐚_𝑖𝑗] is any point kinetics-like matrix, that is, 𝑎_𝑖𝑗=0 except when 𝑖=1 or 𝑗=1 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 [𝐈−𝜀𝐀]⁻¹[𝐈−𝜀𝐀]⁻¹=𝐈−2𝜀_𝑅𝐀+|𝜀|²𝐀²⁻¹ (23) 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: [𝐈−𝜀𝐀]⁻¹[𝐈−𝜀𝐀]⁻¹=Re𝜀𝛾𝛾²𝐚₁𝐛^𝑇₂+𝐚₂𝐛^𝑇₁+Im𝜀𝛾𝛾²𝜀_𝐼⁻¹𝐚₁𝐛^𝑇₁−𝜀_𝐼𝐚₂𝐛^𝑇₂+𝐃, (24) where 𝛾=𝛾_𝑅+𝑖𝛾_𝐼,𝛾_𝑅=1−𝜀_𝑅𝑎₁₁+^𝐺+1_𝑖=2𝜁_𝑖𝑎_1𝑖𝑎_𝑖1𝜀²_𝑅−𝜀²_𝐼−𝜀_𝑅𝑎_𝑖𝑖|𝜀|²,𝛾_𝐼=𝜀_𝐼𝑎₁₁+^𝐺+1_𝑖=2𝜁_𝑖𝑎_1𝑖𝑎_𝑖12𝜀_𝑅𝜀_𝐼−𝜀_𝐼𝑎_𝑖𝑖|𝜀|²,𝑎₁=1,𝜀_𝑅−𝑎₂₂|𝜀|²𝑎₂₁𝜁₂𝜀,…,_𝑅−𝑎_{𝐺+1,𝐺+1}|𝜀|²𝑎_𝐺+1,1𝜁_𝐺+1^𝑇,𝑎₂=0,𝑎₂₁𝜁₂,…,𝑎_𝐺+1,1𝜁_𝐺+1^𝑇,𝑏₁=1,𝜀_𝑅−𝑎₂₂|𝜀|²𝑎₁₂𝜁₂𝜀,…,_𝑅−𝑎_{𝐺+1,𝐺+1}|𝜀|²𝑎_1,𝐺+1𝜁_𝐺+1^𝑇,𝑏₂=0,𝑎₁₂𝜁₂,…,𝑎_1,𝐺+1𝜁_𝐺+1^𝑇,𝐷=diag0,𝜁₂,…,𝜁_𝐺+1, (25) with 𝜁_𝑖=1−2𝜀_𝑅𝑎_𝑖𝑖+|𝜀|²𝑎²_𝑖𝑖⁻¹,𝑖=2,…,𝐺+1. (26) 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 𝜀_𝑅=𝜀_𝐼=ℎ/2.

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 exp(ℎ𝐴) is obtained if the remaining (𝑝+𝑞+1) coefficients are chosen so that the Maclaurin expansion of 𝑅_𝑝𝑞(𝑥) agrees with the Maclaurin expansion of exp(𝑥) to (𝑝+𝑞+1) terms. As such, the Padé approximations of exp(𝑥) appear as local approximations of exp(𝑥) at 𝑥=0 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, exp(ℎ𝐴) is then approximated by some simple functions of (ℎ𝐴), typically a real rational function of type (𝑝,𝑞) given by (16), with usual renormalization (𝑎₀ or 𝑏₀=1).

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: exp(ℎ𝐀)=^∞_𝑘=0(ℎ𝐀)^𝑘,𝑅𝑘!_𝑝𝑞ℎ𝐴=^𝑝_𝑘=0𝑎_𝑘ℎ𝐴^𝑘𝑞_ℓ=0𝑏_ℓℎ𝐴^ℓ−1, (27) the following system of algebraic equations of the unknown coefficients 𝑎_𝑘’s and 𝑏_ℓ’s are obtained [29]: ^𝑘_𝑗=0𝑐_{𝑛−𝑠−𝑗}𝑏_𝑗𝑐=0,𝑠=0,1,…,𝑛−𝑚−1,_𝑖=0if𝑖<0,𝑏₀𝑎=1,𝑛=𝑙+𝑚,_𝑟=^𝑟_𝑗=0𝑐_𝑟−𝑗𝑏_𝑗𝑏𝑟=0,1,…,𝑚,_𝑗=0if𝑗>𝑘. (28) 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 exp(𝑥) is a one-point Padé approximation this point being 𝑥=0, and we have 𝑅^(𝑛)_𝑝𝑞(0)=exp^(𝑛)(0),𝑛=0,1,…,𝑝+𝑞.

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 Ψ(0)=Ψ₀, 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 Re(𝜆)<0 [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 Ψ_𝑛+1=𝑅_𝑝𝑞(ℎ𝜆)Ψ_𝑛, (29) where ℎ=Δ𝑡 and A-stability requires 𝑅_𝑝𝑞(ℎ𝜆)≺1,∀𝜆withRe(𝜆)≺0,∀ℎ≻0, (30) where 𝑅_𝑝𝑞(ℎ𝜆) is any of the rational (𝑝,𝑞) approximations considered. If Re(𝜆)≻0, 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 𝑅_𝑝𝑞((𝜆)𝑅ℎ𝜆)≤1onRe=0,_𝑝𝑞(ℎ𝜆)analyticinRe(𝜆)≺0, (31) 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: 𝑅_𝑝𝑞(ℎ𝜆)≤1,∀𝑝,𝑞≥𝑝,ℎ𝜆realnegative. (32) 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 𝑅_𝑝𝑞(ℎ𝜆)≺1∀ℎ𝜆∈𝑆(𝛼), (33) where 𝑆(𝛼)=ℎ𝜆∈𝐶;Re(ℎ𝜆)≺0,0≤arg(−ℎ𝜆)≺𝛼, (34) then the method is A(𝛼)-stable (or acceptable), where 𝛼∈(0,𝜋/2). 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 𝛼∈(0,𝜋/2) such that 𝑅_𝑝𝑞(ℎ𝜆) is A(𝛼)-acceptable and 𝐴(𝜋/2)-acceptable (A-stable) if and only if 𝑅_𝑝𝑞(ℎ𝜆) is A(𝛼)-acceptable for all 𝛼∈(0,𝜋/2) [39].

Table 1 

Coefficients and explicit Padé approximations with their stability for the exponential function. exp(ℎ𝐀)=∑^𝑝_𝑖=1𝑎_𝑖(ℎ𝐀)^𝑖/∑^𝑞_𝑗=1𝑏_𝑗(ℎ𝐀)^𝑗.

The scheme of rational approximations described in Table 1 showed that the Padé (𝑛,𝑛), (𝑛,𝑛+1), and (𝑛,𝑛+2) approximations lead to A-stable schemes for any 𝑛, while this is not any more true for Padé (𝑛, 𝑛+3) approximation. For example, the Padé (0, 3) is A(𝛼)-stable where 𝛼≤88.23. 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 𝑛(𝑡=0) and 𝑐_𝑖(𝑡=0) 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 𝑐_𝑖(0)=(𝛽_𝑖/𝜆_𝑖Λ)𝑛₀. 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.

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.

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 𝜌=0.0007𝑡 for thermal reactor and 𝜌=0.0044𝑡 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 𝜌=0.0007𝑡 for thermal reactor and 𝜌=0.0044𝑡 for fast reactor.

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 𝜌(𝑡)=()()0.0075𝑡0≤𝑡≤0.5−0.0075𝑡−0.5+0.003750.5≤𝑡≤10.0075𝑡−11≤𝑡≤1.50.003751.5≤𝑡 (35) 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.

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: 𝜌(𝑡)=𝜌₀+𝜇sin(𝜋/180)𝑡, 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 sin(𝜋/180)𝑡 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 𝛼 (𝛼>0) when a small step reactivity 𝜌₀ (𝜌₀<𝛽) is inserted. When we consider the temperature feedback, the real reactor reactivity is given by 𝜌=𝜌₀−𝛼𝑇(𝑡), (36) 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: 𝑑𝑇(𝑡)𝑑𝑡=𝐾_𝑐𝑁(𝑡), (37) where 𝐾_𝑐 is the reciprocal of thermal capacity of reactor. Combining (36) and (37) results in 𝑑𝜌(𝑡)𝑑𝑡=−𝛼𝐾_𝑐𝑛(𝑡). (38) Generally, 𝜌(𝑡)=𝑎(𝑡)+𝑏(𝑡)^𝑡₀𝑁𝑡^′𝑑𝑡^′, (39) 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.

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.