# Generalized and Stability Rational Functions for Dynamic Systems of Reactor Kinetics

Authors: Ahmed E. Aboanber
Journal: International Journal of Nuclear Energy
 (2013)

---

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


---

## Body

## 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 Section2 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 A(α) 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)dn(t)dt=A(t)n(t)+∑i=1GλiCi(t)+f(t),dCi(t)dt=Bi(t)n(t)-λiCi(t),
with
(2)A(t)=ρ(t)-β(t)Λ(t),Bi(t)=βiΛ(t)fori=1,2,…,G,
where ρ(t) is the absolute reactivity defined as = (k-1)/k=kex/k, which is generally a function of time, n(t) is the reactor power (or neutron density), k is the effective multiplication factor, Λ(t) is the neutron generation time, Ci(t) is the precursor concentrations of a group i and λi is the average decay constant, βi is the effective delayed neutron fraction of a group, and f(t) is the source.Equations (1) may be compactly written in vector form. If
(3)Ψ(t)=[n(t)C1(t)⋮CG(t)]=[Ψ1(t)Ψ2(t)⋮ΨG+1(t)],S(t)=[f(t)0⋮0]
denotes a G+1 component column vector and
(4)A(t)=[ρ(t)-βΛλ1⋯λGβ1Λ-λ1O⋮⋱βGΛO-λG]=[aij]
denotes a (G+1)×(G+1) matrix, then (1) are simply
(5)Ψ•(t)=A(t)·Ψ(t)+S(t).A wealth of more or less economical methods have been devised to solve the system of (5). If the matrix elements [aij], 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 ρ(t) (and therefore A(t)) 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 e-tA(t), we obtained the exact solution
(6)ddt(e-tA(t)Ψ¯)=e-tA(t)Si,
where the solution within each time step is sought starting from a piecewise constant approximation substituted in the interval [tk,tk+1], with a constant, namely, Ak, which can be its value at some point in the interval (e.g., Ak≈A(tk)orAk≈1/2(A(tk)+A(tk+1))) as
(7)ρ(t)≈ρ(tk+tk+12)=ρk,S(t)≈S(tk+tk+12)=Skfortk≤t≤tk+1.
Integrating both sides of the equation over the ith time step from ti to ti+1 gives
(8)e-tA(t)Ψ¯i+1-e-tA(t)Ψ¯i=∫tktk+1e-tA(t)uSidu.
Applying the usual constant parameters assumption, where A(t) and S(t) are replaced by their mean values A- and S-, on the interval under consideration with h=tk+1-tk. Equation (8) reduces after subtracting and evaluating the integral to the equation of the form
(9)Ψ¯i+1=eA-hΨ¯i+[eA-h-I](A-)-1Si,whereh=ti+1-ti.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 matricesA and AT, Abaonber and Nahla [30], provides us with a mathematical formula, which permits to approximate some special functions. Let the matrix A have eigenvalues ωi and eigenvectors Ui  (i=1,2,…,ℓ). Also, let AT have corresponding eigenvalues Vi, similar to its eigenvalues ωi which are the same as those of A. Thus, the following theorem is valid for any f(A) for which f(ωi) is bounded for all i. If AUi=ωiUi and ATVi=ωiVi with UiTVi=1 that is normalized to unity, so the following expression holds for any F(A) where F(A) represents a special type of functions (e.g., exponential functions, Trigonometric functions, and Hyperbolic functions)
(10)F(A)=f(A)+∑i=1ℓ[F(ωi)-f(ωi)]UiViT,
where the eigenvalues of both the matrices A and AT 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(hA-)) in an economical manner. The important note here is that if f(ωi) is a good approximation for F(ωi), the ith term in the summation should drop and consequently it will have a very small coefficient,
(11)[F(ωi)-f(ωi)]⪻1.Sinceexp(hωi)≈f(hωi), thus, to a high degree of accuracy, we have
(12)exp(hA)≅g(hA)=f(hA)+∑k/[exp(hωk)-f(hωk)]UkVkT,
where the sum ∑k/ is over only those k for which (11) does not hold.The vectorsUk and Vk are easily calculated from their defining equation:
(13){ωI-A}Uk=0,soUk=col[1μ1λ1+ωk⋯μGλG+ωk],{ωI-AT}Vk=0,soVk=col[1λ1(λ1+ωk)⋯λG(λG+ωk)],
or in the normalized form:
(14)Vk=νkcol[1λ1(λ1+ωk)⋯λG(λG+ωk)],
where νk is the normalization factor, satisfying the normalization condition UiTVi=1, given by
(15)νk=[1+∑i=1Gμiλi(λi+ωk)2]-1≺1.In this connection the analytical inversion method is applied to the point kinetic matrixA. For numerical purpose, exp(hA-) is then approximated by some simple function of (hA-), typically a real rational function of type (p,q):
(16)Rpq(hA-)=Pp(hA-)Qq(hA-)=(∑k=0pak(hA-)k)(∑ℓ=0qbℓ(hA-)ℓ)-1,
with some normalization (usually a0 or b0=1).Whatever the degree,q, of the denominator, it can always be factored into real linear and quadratic factors:
(17)∑ℓ=0qbℓxℓ=∏i=1s1Li(x)∏j=1s2qi(x),s1+2s2=q.
For the sake of simplicity, suppose that S-=0, and Ψi+1 is approximated by
(18)Ψi+1=[∑k=0pak(hA-)k][∑ℓ=0qbℓ(hA-)ℓ]-1Ψi
which can be obtained by solving recursively the matrix problems:
(19)qj(hA-)Φj=Φj-1,j=1,…,s2,Li(hA-)Φs2+i=Φs2+i-1,i=1,…,s2
with
(20)Φi=[∑k=0pak(hA-)k]Ψi,Ψi+1=Φs2+s1.
Due to the particular structure of A-, 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 [I-εA], where for a real ε and a point kinetic matrix A(t) of (4) with the reactivity ρ(t) appropriately averaged, the following expression is introduced:
(21)[I-εA]-1=γabT+D,
where, generally for G delayed neutron group:
(22)γ=[1-εa11+∑i=2G+1ε2a1iai11-εaii]-1,a=[1,εa211-εa22,εa311-εa33,…,εaG+1,11-εaG+1,G+1]T,b=[1,εa121-εa22,εa131-εa33,…,εa1,G+11-εaG+1,G+1]T,D=diag[0,11-εa22,11-εa33,…,11-εaG+1,G+1],
Here, A=[aij] is any point kinetics-like matrix, that is, aij=0 except when i=1 or j=1 or i=j. Consequently, the previous expression is more general compared with that introduced by Aboanber and Nahla [30]. The bℓ in (18) being (generally) real numbers, the quadratic factors qi 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 εn form a complex conjugate pair. In this connection the following pair of factors
(23)[I-εA]-1[I-ε-A]-1=[I-2εRA+|ε|2A2]-1
is considered, with ε=εR+iεI, and ε-=εR-iεI 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:
(24)[I-εA]-1[I-ε-A]-1=Re(εγ)|γ2|(a1b2T+a2b1T)+Im(εγ)|γ2|(εI-1a1b1T-εIa2b2T)+D,
where
(25)γ=γR+iγI,γR=1-εRa11+∑i=2G+1ζia1iai1(εR2-εI2-εRaii|ε|2),γI=εIa11+∑i=2G+1ζia1iai1(2εRεI-εIaii|ε|2),a1=[1,(εR-a22|ε|2)a21ζ2,…,(εR-aG+1,G+1|ε|2)aG+1,1ζG+1]T,a2=[0,a21ζ2,…,aG+1,1ζG+1]T,b1=[1,(εR-a22|ε|2)a12ζ2,…,(εR-aG+1,G+1|ε|2)a1,G+1ζG+1]T,b2=[0,a12ζ2,…,a1,G+1ζG+1]T,D=diag[0,ζ2,…,ζG+1],
with
(26)ζi=[1-2εRaii+|ε|2aii2]-1,i=2,…,G+1.
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 εR=εI=h/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 μi=βi/Λ 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 μi 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 equationdΨ/dt=λΨ.Given a rational functionRpq(x), the Padé (p+q) approximation to exp(hA) is obtained if the remaining (p+q+1) coefficients are chosen so that the Maclaurin expansion of Rpq(x) agrees with the Maclaurin expansion of exp(x) to (p+q+1) terms. As such, the Padé approximations of exp(x) appear as local approximations of exp(x) at x=0 and are consequently rather poor approximations for large x 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(hA) is then approximated by some simple functions of (hA), typically a real rational function of type (p,q) given by (16), with usual renormalization (a0 or b0=1).In the case of exponential function of (16), the explicit expressions for the ak’s and bℓ’s are obtained. Dealing with the following representative equations for the exponential functions:
(27)exp(hA)=∑k=0∞(hA)kk!,Rpq(hA-)=(∑k=0pak(hA-)k)(∑ℓ=0qbℓ(hA-)ℓ)-1,
the following system of algebraic equations of the unknown coefficients ak’s and bℓ’s are obtained [29]:
(28)∑j=0kcn-s-jbj=0,s=0,1,…,n-m-1,ci=0ifi<0,b0=1,n=l+m,ar=∑j=0rcr-jbjr=0,1,…,m,bj=0ifj>k.
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é (p,q) approximation to exp(x) is a one-point Padé approximation this point being x=0, and we have Rpq(n)(0)=exp(n)(0),n=0,1,…,p+q.


### 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)=Ψ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 h 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 h 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
(29)Ψn+1=Rpq(hλ)Ψn,
where h=Δt and A-stability requires
(30)|Rpq(hλ)|≺1,∀λwithRe(λ)≺0,∀h≻0,
where Rpq(hλ) is any of the rational (p,q) 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
(31)|Rpq(hλ)|≤1onRe(λ)=0,Rpq(hλ)analyticinRe(λ)≺0,
as a consequence of the maximum modulus theorem. Varga [38] has shown that all the Padé approximations with q≥p are absolutely stable on the real negative axis:
(32)|Rpq(hλ)|≤1,∀p,q≥p,hλrealnegative.
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 Rpq(hλ) is such that
(33)|Rpq(hλ)|≺1∀hλ∈S(α),
where
(34)S(α)={hλ∈C;Re(hλ)≺0,0≤|arg(-hλ)|≺α},
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 Rpq(hλ) is A(α)-acceptable and A(π/2)-acceptable (A-stable) if and only if Rpq(hλ) is A(α)-acceptable for all α∈(0,π/2) [39].Table 1
Coefficients and explicit Padé approximations with their stability for the exponential function.exp(hA)=∑i=1pai(hA)i/∑j=1qbj(hA)j.
Approximation
Stability
P
01
a
0
0
b
0
1
A-Stable
a
1
0
b
1
−1
a
2
0
b
2
0
P
11
a
0
0
b
0
0
A-Stable
a
1
1
/
2
b
1
1
/
2
a
2
0
b
2
0
P
11
n
a
0
0
b
0
1
A(59.62)-Stable
a
1
1
/
2
b
1
-
1
/
4
a
2
0
b
2
0
P
02
a
0
0
b
0
1
A-Stable
a
1
0
b
1
−1
a
2
0
b
2
1
/
2
P
12
a
0
0
b
0
1
A-Stable
a
1
1
/
3
b
1
-
2
/
3
a
2
0
b
2
1
/
6
P
12
n
a
0
0
b
0
1
A(82.81)-Stable
a
1
1
/
3
b
1
1
/
3
a
2
0
b
2
1
/
12
P
22
a
0
0
b
0
1
A-Stable
a
1
1
/
2
b
1
-
1
/
2
a
2
1
/
12
b
2
1
/
12The scheme of rational approximations described in Table1 showed that the Padé (n,n), (n,n+1), and (n,n+2) approximations lead to A-stable schemes for any n, while this is not any more true for Padé (n, n+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.


## 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)dn(t)dt=A(t)n(t)+∑i=1GλiCi(t)+f(t),dCi(t)dt=Bi(t)n(t)-λiCi(t),
with
(2)A(t)=ρ(t)-β(t)Λ(t),Bi(t)=βiΛ(t)fori=1,2,…,G,
where ρ(t) is the absolute reactivity defined as = (k-1)/k=kex/k, which is generally a function of time, n(t) is the reactor power (or neutron density), k is the effective multiplication factor, Λ(t) is the neutron generation time, Ci(t) is the precursor concentrations of a group i and λi is the average decay constant, βi is the effective delayed neutron fraction of a group, and f(t) is the source.Equations (1) may be compactly written in vector form. If
(3)Ψ(t)=[n(t)C1(t)⋮CG(t)]=[Ψ1(t)Ψ2(t)⋮ΨG+1(t)],S(t)=[f(t)0⋮0]
denotes a G+1 component column vector and
(4)A(t)=[ρ(t)-βΛλ1⋯λGβ1Λ-λ1O⋮⋱βGΛO-λG]=[aij]
denotes a (G+1)×(G+1) matrix, then (1) are simply
(5)Ψ•(t)=A(t)·Ψ(t)+S(t).A wealth of more or less economical methods have been devised to solve the system of (5). If the matrix elements [aij], 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 ρ(t) (and therefore A(t)) 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 e-tA(t), we obtained the exact solution
(6)ddt(e-tA(t)Ψ¯)=e-tA(t)Si,
where the solution within each time step is sought starting from a piecewise constant approximation substituted in the interval [tk,tk+1], with a constant, namely, Ak, which can be its value at some point in the interval (e.g., Ak≈A(tk)orAk≈1/2(A(tk)+A(tk+1))) as
(7)ρ(t)≈ρ(tk+tk+12)=ρk,S(t)≈S(tk+tk+12)=Skfortk≤t≤tk+1.
Integrating both sides of the equation over the ith time step from ti to ti+1 gives
(8)e-tA(t)Ψ¯i+1-e-tA(t)Ψ¯i=∫tktk+1e-tA(t)uSidu.
Applying the usual constant parameters assumption, where A(t) and S(t) are replaced by their mean values A- and S-, on the interval under consideration with h=tk+1-tk. Equation (8) reduces after subtracting and evaluating the integral to the equation of the form
(9)Ψ¯i+1=eA-hΨ¯i+[eA-h-I](A-)-1Si,whereh=ti+1-ti.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 matricesA and AT, Abaonber and Nahla [30], provides us with a mathematical formula, which permits to approximate some special functions. Let the matrix A have eigenvalues ωi and eigenvectors Ui  (i=1,2,…,ℓ). Also, let AT have corresponding eigenvalues Vi, similar to its eigenvalues ωi which are the same as those of A. Thus, the following theorem is valid for any f(A) for which f(ωi) is bounded for all i. If AUi=ωiUi and ATVi=ωiVi with UiTVi=1 that is normalized to unity, so the following expression holds for any F(A) where F(A) represents a special type of functions (e.g., exponential functions, Trigonometric functions, and Hyperbolic functions)
(10)F(A)=f(A)+∑i=1ℓ[F(ωi)-f(ωi)]UiViT,
where the eigenvalues of both the matrices A and AT 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(hA-)) in an economical manner. The important note here is that if f(ωi) is a good approximation for F(ωi), the ith term in the summation should drop and consequently it will have a very small coefficient,
(11)[F(ωi)-f(ωi)]⪻1.Sinceexp(hωi)≈f(hωi), thus, to a high degree of accuracy, we have
(12)exp(hA)≅g(hA)=f(hA)+∑k/[exp(hωk)-f(hωk)]UkVkT,
where the sum ∑k/ is over only those k for which (11) does not hold.The vectorsUk and Vk are easily calculated from their defining equation:
(13){ωI-A}Uk=0,soUk=col[1μ1λ1+ωk⋯μGλG+ωk],{ωI-AT}Vk=0,soVk=col[1λ1(λ1+ωk)⋯λG(λG+ωk)],
or in the normalized form:
(14)Vk=νkcol[1λ1(λ1+ωk)⋯λG(λG+ωk)],
where νk is the normalization factor, satisfying the normalization condition UiTVi=1, given by
(15)νk=[1+∑i=1Gμiλi(λi+ωk)2]-1≺1.In this connection the analytical inversion method is applied to the point kinetic matrixA. For numerical purpose, exp(hA-) is then approximated by some simple function of (hA-), typically a real rational function of type (p,q):
(16)Rpq(hA-)=Pp(hA-)Qq(hA-)=(∑k=0pak(hA-)k)(∑ℓ=0qbℓ(hA-)ℓ)-1,
with some normalization (usually a0 or b0=1).Whatever the degree,q, of the denominator, it can always be factored into real linear and quadratic factors:
(17)∑ℓ=0qbℓxℓ=∏i=1s1Li(x)∏j=1s2qi(x),s1+2s2=q.
For the sake of simplicity, suppose that S-=0, and Ψi+1 is approximated by
(18)Ψi+1=[∑k=0pak(hA-)k][∑ℓ=0qbℓ(hA-)ℓ]-1Ψi
which can be obtained by solving recursively the matrix problems:
(19)qj(hA-)Φj=Φj-1,j=1,…,s2,Li(hA-)Φs2+i=Φs2+i-1,i=1,…,s2
with
(20)Φi=[∑k=0pak(hA-)k]Ψi,Ψi+1=Φs2+s1.
Due to the particular structure of A-, 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 [I-εA], where for a real ε and a point kinetic matrix A(t) of (4) with the reactivity ρ(t) appropriately averaged, the following expression is introduced:
(21)[I-εA]-1=γabT+D,
where, generally for G delayed neutron group:
(22)γ=[1-εa11+∑i=2G+1ε2a1iai11-εaii]-1,a=[1,εa211-εa22,εa311-εa33,…,εaG+1,11-εaG+1,G+1]T,b=[1,εa121-εa22,εa131-εa33,…,εa1,G+11-εaG+1,G+1]T,D=diag[0,11-εa22,11-εa33,…,11-εaG+1,G+1],
Here, A=[aij] is any point kinetics-like matrix, that is, aij=0 except when i=1 or j=1 or i=j. Consequently, the previous expression is more general compared with that introduced by Aboanber and Nahla [30]. The bℓ in (18) being (generally) real numbers, the quadratic factors qi 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 εn form a complex conjugate pair. In this connection the following pair of factors
(23)[I-εA]-1[I-ε-A]-1=[I-2εRA+|ε|2A2]-1
is considered, with ε=εR+iεI, and ε-=εR-iεI 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:
(24)[I-εA]-1[I-ε-A]-1=Re(εγ)|γ2|(a1b2T+a2b1T)+Im(εγ)|γ2|(εI-1a1b1T-εIa2b2T)+D,
where
(25)γ=γR+iγI,γR=1-εRa11+∑i=2G+1ζia1iai1(εR2-εI2-εRaii|ε|2),γI=εIa11+∑i=2G+1ζia1iai1(2εRεI-εIaii|ε|2),a1=[1,(εR-a22|ε|2)a21ζ2,…,(εR-aG+1,G+1|ε|2)aG+1,1ζG+1]T,a2=[0,a21ζ2,…,aG+1,1ζG+1]T,b1=[1,(εR-a22|ε|2)a12ζ2,…,(εR-aG+1,G+1|ε|2)a1,G+1ζG+1]T,b2=[0,a12ζ2,…,a1,G+1ζG+1]T,D=diag[0,ζ2,…,ζG+1],
with
(26)ζi=[1-2εRaii+|ε|2aii2]-1,i=2,…,G+1.
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 εR=εI=h/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 μi=βi/Λ 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 μi 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 equationdΨ/dt=λΨ.Given a rational functionRpq(x), the Padé (p+q) approximation to exp(hA) is obtained if the remaining (p+q+1) coefficients are chosen so that the Maclaurin expansion of Rpq(x) agrees with the Maclaurin expansion of exp(x) to (p+q+1) terms. As such, the Padé approximations of exp(x) appear as local approximations of exp(x) at x=0 and are consequently rather poor approximations for large x 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(hA) is then approximated by some simple functions of (hA), typically a real rational function of type (p,q) given by (16), with usual renormalization (a0 or b0=1).In the case of exponential function of (16), the explicit expressions for the ak’s and bℓ’s are obtained. Dealing with the following representative equations for the exponential functions:
(27)exp(hA)=∑k=0∞(hA)kk!,Rpq(hA-)=(∑k=0pak(hA-)k)(∑ℓ=0qbℓ(hA-)ℓ)-1,
the following system of algebraic equations of the unknown coefficients ak’s and bℓ’s are obtained [29]:
(28)∑j=0kcn-s-jbj=0,s=0,1,…,n-m-1,ci=0ifi<0,b0=1,n=l+m,ar=∑j=0rcr-jbjr=0,1,…,m,bj=0ifj>k.
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é (p,q) approximation to exp(x) is a one-point Padé approximation this point being x=0, and we have Rpq(n)(0)=exp(n)(0),n=0,1,…,p+q.


## 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)=Ψ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 h 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 h 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
(29)Ψn+1=Rpq(hλ)Ψn,
where h=Δt and A-stability requires
(30)|Rpq(hλ)|≺1,∀λwithRe(λ)≺0,∀h≻0,
where Rpq(hλ) is any of the rational (p,q) 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
(31)|Rpq(hλ)|≤1onRe(λ)=0,Rpq(hλ)analyticinRe(λ)≺0,
as a consequence of the maximum modulus theorem. Varga [38] has shown that all the Padé approximations with q≥p are absolutely stable on the real negative axis:
(32)|Rpq(hλ)|≤1,∀p,q≥p,hλrealnegative.
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 Rpq(hλ) is such that
(33)|Rpq(hλ)|≺1∀hλ∈S(α),
where
(34)S(α)={hλ∈C;Re(hλ)≺0,0≤|arg(-hλ)|≺α},
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 Rpq(hλ) is A(α)-acceptable and A(π/2)-acceptable (A-stable) if and only if Rpq(hλ) is A(α)-acceptable for all α∈(0,π/2) [39].Table 1
Coefficients and explicit Padé approximations with their stability for the exponential function.exp(hA)=∑i=1pai(hA)i/∑j=1qbj(hA)j.
Approximation
Stability
P
01
a
0
0
b
0
1
A-Stable
a
1
0
b
1
−1
a
2
0
b
2
0
P
11
a
0
0
b
0
0
A-Stable
a
1
1
/
2
b
1
1
/
2
a
2
0
b
2
0
P
11
n
a
0
0
b
0
1
A(59.62)-Stable
a
1
1
/
2
b
1
-
1
/
4
a
2
0
b
2
0
P
02
a
0
0
b
0
1
A-Stable
a
1
0
b
1
−1
a
2
0
b
2
1
/
2
P
12
a
0
0
b
0
1
A-Stable
a
1
1
/
3
b
1
-
2
/
3
a
2
0
b
2
1
/
6
P
12
n
a
0
0
b
0
1
A(82.81)-Stable
a
1
1
/
3
b
1
1
/
3
a
2
0
b
2
1
/
12
P
22
a
0
0
b
0
1
A-Stable
a
1
1
/
2
b
1
-
1
/
2
a
2
1
/
12
b
2
1
/
12The scheme of rational approximations described in Table1 showed that the Padé (n,n), (n,n+1), and (n,n+2) approximations lead to A-stable schemes for any n, while this is not any more true for Padé (n, n+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 n(t=0) and ci(t=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 n0 is given by ci(0)=(βi/λiΛ)n0. 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.
i
Fast reactor
Thermal reactor (I)
Thermal reactor (II)
Thermal reactor (III)
λ
i
(
sec
-
1
)
β
i
(
pcm
)
λ
i
(
sec
-
1
)
β
i
(
pcm
)
λ
i
(
sec
-
1
)
β
i
(
pcm
)
λ
i
(
sec
-
1
)
β
i
(
pcm
)
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
β
=
440 (pcm) Λ=1.0(-7) sec
β
=
750 (pcm) Λ=5.0(-4) sec
β
=
700 (pcm) Λ=2.0(-5) sec
β
=
645 (pcm)Λ=5.0(-5) secBenchmark 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−4 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 ρ=0.0007t for thermal reactor and ρ=0.0044t 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.0007t for thermal reactor and ρ=0.0044t 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 Table4. 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 Table4 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
(35)ρ(t)={0.0075t0≤t≤0.5-0.0075(t-0.5)+0.003750.5≤t≤10.0075(t-1)1≤t≤1.50.003751.5≤t
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: ρ(t)=ρ0+μsin(π/180)t, 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 ρ0 is the constant part of the excess reactivity, and sin(π/180)t 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 ρ0 (ρ0<β) is inserted. When we consider the temperature feedback, the real reactor reactivity is given by
(36)ρ=ρ0-αT(t),
where T is the temperature increment of the reactor. After the reactivity ρ0 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)dT(t)dt=KcN(t),
where Kc is the reciprocal of thermal capacity of reactor. Combining (36) and (37) results in
(38)dρ(t)dt=-αKcn(t).
Generally,
(39)ρ(t)=a(t)+b(t)∫0tN(t′)dt′,
where a(t) represented the impressed reactivity variation (generally polynomial in t), and b(t) is the “shutdown coefficient” of the reactor system. It is usually taken as a negative constant of magnitude ranging from 10−13 cm3/sec for slow system to 10−7 cm3/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−4 sec, and shutdown coefficients are in the range of 10−11 cm3/sec to 10−13 cm3/sec.
a
b
Method
First peak
Time to the first peak
First asymptotic neutron density
0.1
1 × 10−11
PWS
2.4238380 × 1011
0.224000
1.0066510 × 1010
GRK4-A
2.4649574 × 1011
0.22626
1.0066499 × 1010
GRK4-T
2.4652371 × 1011
0.22626
1.0066500 × 1010
GRK4-H
2.4654628 × 1011
0.22626
1.0066501 × 1010
Padé 11
2.4238380 × 1011
0.22400
1.0060053 × 1010
1 × 10−13
PWS
2.8666680 × 1013
0.23890
1.0066520 × 1012
GRK4-A
2.9580990 × 1013
0.24066
1.0066480 × 1012
GRK4-T
2.9580157 × 1013
0.24066
1.0066480 × 1012
GRK4-H
2.9577276 × 1013
0.24066
1.0066482 × 1012
Padé 11
2.8929486 × 1013
0.23800
1.0066439 × 1012
0.01
1 × 10−11
PWS
2.0064900 × 1010
1.10617
1.0263960 × 109
GRK4-A
2.0232131 × 1010
1.10736
1.0246410 × 109
GRK4-T
2.0232020 × 1010
1.10736
1.0246424 × 109
GRK4-H
2.0232065 × 1010
1.10736
1.0246394 × 109
Padé 11
2.0120302 × 1010
1.10600
1.0255595 × 109
1 × 10−13
PWS
2.4831920 × 1012
1.15528
9.6705500 × 1010
GRK4-A
2.5061997 × 1012
1.15668
9.6374374 × 1010
GRK4-T
2.5061921 × 1012
1.15668
9.6374817 × 1010
GRK4-H
2.5061852 × 1012
1.15668
9.6374301 × 1010
Padé 11
2.4909973 × 1012
1.15500
9.6629450 × 1010
0.001
1 × 10−11
PWS
1.2652650 × 109
7.48864
1.9212660 × 108
GRK4-A
1.2760915 × 109
7.48728
1.9207147 × 108
GRK4-T
1.2761036 × 109
7.48728
1.9207122 × 108
GRK4-H
1.2760943 × 109
7.48728
1.9206984 × 108
Padé 11
1.2740050 × 109
7.48864
1.9282197 × 108
1 × 10−13
PWS
1.7068560 × 1011
7.68362
2.0046470 × 1010
GRK4-A
1.7245128 × 1011
7.68204
2.0040163 × 1010
GRK4-T
1.7245159 × 1011
7.68204
2.0040077 × 1010
GRK4-H
1.7245002 × 1011
7.68204
2.0039905 × 1010
Padé 11
1.7209210 × 1011
7.68263
2.0067995 × 1010


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



---

*Source: 903904-2013-08-13.xml*