Spectrum behavior for the nonlinear fractional point reactor kinetics model

NUCLEAR ENERGY SCIENCE AND ENGINEERING

Spectrum behavior for the nonlinear fractional point reactor kinetics model

Ahmed E. Aboanber，

Abdallah A. Nahla，

A.A. Hemeda

Nuclear Science and Techniques

Vol.29, No.1

Article number 14

Published in print 01 Jan 2018

Available online 29 Dec 2017

DOI：10.1007/s41365-017-0350-4

46100

The nonlinear fractional point reactor kinetics equations in the presence of Newtonian temperature reactivity feedback with a multi-group of delayed neutrons, which describes the spectrum behavior of neutron density into the homogenous nuclear reactors, is developed. This system is one of the most important stiff coupled nonlinear fractional differentials for nuclear reactor dynamics. The generalization of Taylor’s formula that involves Caputo fractional derivatives is developed in an attempt to overcome the difficulty of the stiffness of the nonlinear fractional differential model. Moreover, the general fractional derivatives are calculated analytically throughout this work. Furthermore, the local and global estimated errors were analyzed, which suggest that the error quantification should take into account the possible grow in time of the error. This observation provides a motivation for going beyond more classical local-in-time concepts of error (Local Truncation Error). The neutron density response with time is analyzed for the anomalous diffusion, sub-diffusion, and supper-diffusion processes.

Nonlinear fractionalGeneralized Taylor’s formulaPoint kineticsMulti-group delayed neutronsTemperature feedback reactivity

1 Introduction

Fractional differential equations have been the focus of many studies due to their frequent appearance in various fields such as reactor kinetics, theoretical physics, fluid mechanics, mathematical biology and finance, chemistry, mathematics, and engineering [1-6]. The Caputo fractional derivative is considered recently in many mathematical physics and engineering problems, [7-14], more than the Riemann-Liouville fractional derivatives, because it allows traditional initial and boundary conditions to be included in the formulation of the problem. In this paper, we use the generalized Taylors formula (GTF) to solve fractional differential equations of nonlinear fractional point reactor kinetics equations. This method is very useful and can be applied to solve many important fractional differential equations with non constant coefficients.

The point kinetics equations model is one of the most important reduced pattern from the neutron transport equation of the nuclear reactor. This model has been the subject of countless studies and applications to understand the transient behavior of neutrons and its effects on nuclear reactor dynamics [15-18]. The presence of temperature reactivity feedback is beneficial to provide an estimate of the transient behavior of reactor power and other system variables of the reactor cores which are fairly tightly coupled. As a result of this reason, some of our previous works [19-23], have been developed to solve analytically and/or numerically the nonlinear point kinetics equations.

The linear and nonlinear point kinetics equations were generalized to fractional neutron point kinetics equations by Ref. [24-26], which were commented on and developed by Ref. [27-31]. Furthermore, the fractional point kinetics equations were discussed by [32-34]. The modified point kinetics equations using Fick’s approximation [35] were presented by [36,37]. As a result of scientific development, our aim of this work is to progress the nonlinear fractional point kinetics model. A technique based on the generalized Taylor’s series in the fractional form has been selected as the best or most appropriate to solve this nonlinear fractional model. Furthermore, the neutron density response with time is estimated and analyzed in two cases, sub-diffusion and supper-diffusion processes of nuclear reactor dynamics.

This paper is organized as follows: The generalized Taylor’s formula (GTF) is developed for solving the nonlinear fractional point reactor kinetics equations with multi-group of delayed neutrons and the adiabatic feedback model in Sect. 2. In addition, the fractional derivatives of the neutron density, the precursor concentrations of multi-group delayed neutrons, and the temperature reactivity feedback is derived. In Sect. 3, the local and global errors were analyzed. The numerical results of the proposed method for different types of reactivities and various fractional orders are discussed in Sect. 4. The conclusion, with a brief summary of the main feedback, is introduced in Sect. 5.

2 Preliminaries

In this section, we mention the following basic definitions of fractional calculus, which are used further in the present work.

Definition 1[2,3,9,10]: Cauchy’s formula for repeated integration of order n∈N is defined as

\begin{array}{l} J^{n} f (t) & = \frac{1}{(n - 1)!} \int_{t_{0}}^{t} {(t - τ)}^{n - 1} f (τ) d τ \\ = \frac{1}{Γ (n)} \int_{t_{0}}^{t} {(t - τ)}^{n - 1} f (τ) d τ . \end{array}

(1)

Definition 2[3]: The Riemann-Liouville fractional integral operator of order α>0 of a function f(t) is defined as follows:

J^{α} f (t) = \frac{1}{Γ (α)} \int_{t_{0}}^{t} {(t - τ)}^{α - 1} f (τ) d τ .

(2)

Definition 3[3]: The Riemann-Liouville derivative of order α, which n-1<α<n and n∈N, defines as

D^{α} f (t) = \frac{1}{Γ (n - α)} \frac{d^{n}}{d t^{n}} \int_{t_{0}}^{t} \frac{f (τ)}{{(t - τ)}^{α + 1 - n}} d τ .

(3)

According to the Riemann-Liouville fractional derivative (3), we have

\begin{array}{l} \lim_{α \to n} D^{α} f (t) = f^{(n)} (t) = \frac{d^{n} f (t)}{d t^{n}}, \forall n \in N \\ D^{α} c = \frac{c}{Γ (1 - α)} t^{- α} \neq 0, c is aconstant . \\ D^{α} (D^{γ} f (t)) = D^{γ} (D^{α} f (t)) = D^{α + γ} f (t) . \end{array}

(4)

Definition 4[1]: The Caputo fractional derivative of order α, which m-1< α< m and m∈N, is defined as follows:

D_{t}^{α} f (t) = \frac{1}{Γ (m - α)} \int_{t_{0}}^{t} \frac{f^{(m)} (τ)}{{(t - τ)}^{α + 1 - m}} d τ .

(5)

Based on the Caputo definition (5), the following properties for the Caputo operator hold

\begin{array}{l} \lim_{α \to m} D_{t}^{α} f (t) = f^{(m)} (t), \\ \lim_{α \to m - 1} D_{t}^{α} f (t) = f^{(m - 1)} (t) - f^{(m - 1)} (t_{0}), \\ D_{t}^{α} c = 0 c is aconstant . \\ D_{t}^{α} (D_{t}^{γ} f (t)) \neq D_{t}^{γ} (D_{t}^{α} f (t)) \neq D_{t}^{α + γ} f (t) (In general) \end{array}

(6)

Let us assume that the fractional order, α, tends to one, i.e. 1-ϵ<α< 1+ϵ and ϵ→ 0, and $\frac{d f (t)}{d t} |_{t = t_{0}} = 0$ . Then, it is easy to prove that the following properties for the Caputo operator hold [9,10]

\begin{array}{l} \lim_{α \to 1^{-}} D_{t}^{α} f (t) = f^{'} (t), from case 0 < α < 1 \\ \lim_{α \to 1^{+}} D_{t}^{α} f (t) = f^{'} (t), from case 1 < α < 2 \\ \Rightarrow \lim_{α \to 1} D_{t}^{α} f (t) = f^{'} (t) . \end{array}

(7)

To suit our case, the fractional order α is chosen so that the fractional derivative laws are approximate to the ordinary derivatives laws. This means that, for a small positive value of parameter ϵ, the fractional derivatives satisfy approximately the following properties

\begin{array}{l} D_{t}^{α} (D_{t}^{γ} f (t)) ≃ D_{t}^{γ} (D_{t}^{α} f (t)) ≃ D_{t}^{α + γ} f (t) \\ D_{t}^{n α} (f (t) g (t)) ≃ \sum_{k = 0}^{n} (\begin{matrix} n \\ k \end{matrix}) (D_{t}^{(n - k) α} f (t)) (D_{t}^{k α} g (t)) . \end{array}

(8)

3 Generalized mathematical model based on Taylor’s formula

The nonlinear point reactor kinetics equations with multi-group of delayed neutrons and temperature reactivity feedback, [38,39], can be generalized to nonlinear fractional differential equations. This system of equation was stretched out by a perturbation in the procedure of a temperature feedback, where the disorder might be comprehended as an adjustment in the nuclear framework design, with an outcome of a modified specific heat flow that incites an adjustment in the temperature. Since the source of heat production is the nuclear processes, the variation of the thermal change is related linearly to the neutron density in its simplest form.

D_{t}^{α} N (t) = (\frac{ρ (t) - β}{l}) N (t) + \sum_{i = 1}^{I} λ_{i} C_{i} (t)

(9)

D_{t}^{α} C_{i} (t) = \frac{β_{i}}{l} N (t) - λ_{i} C_{i} (t), i = 1,2,3, \dots, I

(10)

ρ (t) = γ (t) - σ [T (t) - T_{0}]

(11)

D_{t}^{α} T (t) = K_{c} N (t)

(12)

where N(t) is the neutron density, ρ(t) is the reactivity as a function of time t, Ci(t) is the precursor concentrations of i-group of delayed neutrons, $D_{t}^{α}$ is the Caputo fractional derivative operator of order α, $β = \sum_{i = 1}^{I} β_{i}$ is the total fraction of delayed neutrons, βi is the fraction of i-group of delayed neutrons, λi is the decay constant of i-group of delayed neutrons, l is the prompt neutron generation time, and I is the total number of delayed neutron groups. T(t) is the temperature of the reactor, T₀ is the initial temperature of the reactor, σ is the temperature coefficient of reactivity, γ(t) represents the external reactivity, and Kc is the reciprocal of the thermal capacity of reactor. The proportionality constant, K_c, signifies a parameter for the influence of the change of heat flow on the rate of temperature change. Thus, the linear relation between temperature change rate and neutron density is supplied with a feedback mechanism.

Notice that, according to dimensional metrology, the dimension of the the prompt neutron generation time, l, is (sα) while the dimension of the the decay constant of i-group of delayed neutrons, λi, is (s-α) from Eq. (9). Using the dimension of these parameters in equation (10), we find that the fractional order, α, of the neutron density should be equal to the fractional order of the precursor concentrations of i-group of delayed neutrons. So, it is not correct to use different fractional derivatives in this system. Moreover, In order to preserve the physical meaning of the system parameters and variables, the values of the fractional order α should be close to one. For example, $\frac{d N (t)}{d t}$ means the rate of time change of the neutron density, but $D_{t}^{0.5} N (t)$ doesn’t represent the same physical meaning of $\frac{d N (t)}{d t}$ .

Differentiating Eq. (11) with respect to time, t, and substituting Eq. (12) yield

D_{t}^{α} ρ (t) = D_{t}^{α} γ (t) - σ K_{c} N (t) .

(13)

Equations (1), (2), and (5) represent the system of the nonlinear fractional differential equations. In this paper, we describe a method that is surprisingly simple while maintaining a desired degree of accuracy. The generalized Taylor’s formula [40-44] is suggested to solve the above system of stiff coupled differential equation in its fractional form as follows:

N (t_{m + 1}) = \sum_{n = 0}^{\infty} \frac{h^{n α} D_{t}^{n α} N (t_{m})}{Γ (n α + 1)} = N (t_{m}) + \frac{h^{α} D_{t}^{α} N (t_{m})}{Γ (α + 1)} + \frac{h^{2 α} D_{t}^{2 α} N (t_{m})}{Γ (2 α + 1)} + \dots,

(14)

C_{i} (t_{m + 1}) = \sum_{n = 0}^{\infty} \frac{h^{n α} D_{t}^{n α} C_{i} (t_{m})}{Γ (n α + 1)} = C_{i} (t_{m}) + \frac{h^{α} D_{t}^{α} C_{i} (t_{m})}{Γ (α + 1)} + \frac{h^{2 α} D_{t}^{2 α} C_{i} (t_{m})}{Γ (2 α + 1)} + = 1,2,3, \dots, I,

(15)

ρ (t_{m + 1}) = \sum_{n = 0}^{\infty} \frac{h^{n α} D_{t}^{n α} ρ (t_{m})}{Γ (n α + 1)} = ρ (t_{m}) + \frac{h^{α} D_{t}^{α} ρ (t_{m})}{Γ (α + 1)} + \frac{h^{2 α} D_{t}^{2 α} ρ (t_{m})}{Γ (2 α + 1)} + \dots,

(16)

where h is the length of time interval [tm,t_m+1].

The fractional derivatives of order α for the neutron density, the precursor concentrations of the i-group of delayed neutrons, and the reactivity are

D_{t}^{α} N (t_{m}) = (\frac{ρ (t_{m}) - β}{l}) N (t_{m}) + \sum_{i = 1}^{I} λ_{i} C_{i} (t_{m}),

(17)

D_{t}^{α} C_{i} (t_{m}) = \frac{β_{i}}{l} N (t_{m}) - λ_{i} C_{i} (t_{m}), i = 1,2,3, \dots, I,

(18)

D_{t}^{α} ρ (t_{m}) = D_{t}^{α} γ (t_{m}) - σ K_{c} N (t_{m}) .

(19)

The fractional derivatives of order 2α for the neutron density, the precursor concentrations of the i-group of delayed neutrons, and the reactivity are given by:

D_{t}^{2 α} N (t_{m}) ≃ - \frac{β}{l} D_{t}^{α} N (t_{m}) + \sum_{i = 1}^{I} λ_{i} D_{t}^{α} C_{i} (t_{m}) + \frac{1}{l} (ρ (t_{m}) D_{t}^{α} N (t_{m}) + N (t_{m}) D_{t}^{α} ρ (t_{m})),

(20)

D_{t}^{2 α} C_{i} (t_{m}) ≃ \frac{β_{i}}{l} D_{t}^{α} N (t_{m}) - λ_{i} D_{t}^{α} C_{i} (t_{m}), i = 1,2,3, \dots, I,

(21)

D_{t}^{2 α} ρ (t_{m}) ≃ D_{t}^{2 α} γ (t_{m}) - σ K_{c} D_{t}^{α} N (t_{m}) .

(22)

Similarly, the fractional derivatives of order nα for the neutron density, the precursor concentrations of the i-group of delayed neutrons, and the reactivity can be written as follows

\begin{array}{l} D_{t}^{n α} N (t_{m}) ≃ - \frac{β}{l} D_{t}^{(n - 1) α} N (t_{m}) + \sum_{i = 1}^{I} λ_{i} D_{t}^{(n - 1) α} C_{i} (t_{m}) \\ + \frac{1}{l} \sum_{j = 0}^{n - 1} (\begin{matrix} n - 1 \\ j \end{matrix}) (D_{t}^{j α} ρ (t_{m})) (D_{t}^{(n - j - 1) α} N (t_{m})), \end{array}

(23)

D_{t}^{n α} C_{i} (t_{m}) ≃ \frac{β_{i}}{l} D_{t}^{(n - 1) α} N (t_{m}) - λ D_{t}^{(n - 1) α} C_{i} (t_{m}), i = 1,2,3, \dots, I,

(24)

D_{t}^{n α} ρ (t_{m}) ≃ D_{t}^{n α} γ (t_{m}) - σ K_{c} D_{t}^{(n - 1) α} N (t_{m}) .

(25)

Substituting Eqs. (17) –(25) into Eqs. (14)–(16) yields the numerical solution of the nonlinear fractional point nuclear reactor kinetics equations in the presence of Newtonian temperature feedback reactivity.

4 Local and Global Estimate Error

The exact solution of the neutron density, N(t_m+1), at time t_m+1 can be written as:

N (t_{m + 1}) = N (t_{m}) + h^{α} F (t_{m}, N (t_{m}), h^{α}) + T E_{n} (ξ),

(26)

where $F (t_{m}, N (t_{m}), h^{α}) = \sum_{k = 1}^{n} \frac{h^{(k - 1) α} D_{t}^{k α} N (t_{m})}{Γ (k α + 1)}$ , and TEn(ξ) is the local truncation error given by:

T E_{n} (ξ) = \frac{h^{(n + 1) α} D_{t}^{(n + 1) α} N (ξ)}{Γ ((n + 1) α + 1)} \leq K h^{(n + 1) α}, ξ \in [t_{m}, t_{m + 1}],

(27)

where K is a constant, which is the maximum value of $\frac{D_{t}^{(n + 1) α} N (ξ)}{Γ ((n + 1) α + 1)}$ between tm and t_m+1.

The approximate solution using the generalized Taylor’s formula of order n is Nn(t_m+1), which is defined by equation (14). Such formula can be rewritten as:

N_{n} (t_{m + 1}) = N_{n} (t_{m}) + h^{α} F (t_{m}, N_{n} (t_{m}), h^{α}) .

(28)

From Eqs. (26) and (28):

| N (t_{m + 1}) - N_{n} (t_{m + 1}) | \leq | N (t_{m}) - N_{n} (t_{m}) | + h^{α} | F (t_{m}, N (t_{m}), h^{α}) - F (t_{m}, N_{n} (t_{m}), h^{α}) | + | T E_{n} (ξ) | .

(29)

The function F(tm,N(tm),hα) is Lipschitz continuous, i.e.

| F (t_{m}, N (t_{m}), h^{α}) - F (t_{m}, N_{n} (t_{m}), h^{α}) | \leq L | N (t_{m}) - N_{n} (t_{m}) |,

(30)

where L is a positive constant value.

Under this condition, the global error, ϵn(t_m+1), takes the following form [45]:

ε_{n} (t_{m + 1}) = | N (t_{m + 1}) - N_{n} (t_{m + 1}) | \leq (1 + h^{α} L) | N (t_{m}) - N_{n} (t_{m}) | + K h^{(n + 1) α} = (1 + h^{α} L) ε_{n} (t_{m}) + K h^{(n + 1) α},

(31)

and consequently, we get:

ε_{n} (t_{m}) \leq (1 + h^{α} L) ε_{n} (t_{m - 1}) + K h^{(n + 1) α} .

(32)

Now, setting η=(1+hα L) and χ= K h^(n+1)α, from equation (32) we get,

\begin{array}{l} ε_{n} (t_{m}) \leq η ε_{n} (t_{m - 1}) + χ \leq η (η ε_{n} (t_{m - 2}) + χ) + χ = η^{2} ε_{n} (t_{m - 2}) + χ (1 + η) \\ \leq η^{2} (μ ε_{n} (t_{m - 3}) + χ) + χ (1 + η) = η^{3} ε_{n} (t_{m - 3}) + χ (1 + η + η^{2}) \\ \leq \dots \leq η^{m} ε_{n} (t_{0}) + χ (1 + η + η^{2} + \dots + η^{m - 1}) . \end{array}

(33)

The initial global error is equal to zero, ϵn(t₀)=0, and then, from Eq.(33), we obtain

\begin{array}{l} ε_{n} (t_{m}) \leq χ (1 + η + η^{2} + \dots + η^{m - 1}) & = χ \frac{1 - η^{m}}{1 - η} \\ = χ \frac{{(1 + h^{α} L)}^{m} - 1}{h^{α} L} . \end{array}

(34)

Or, in an alternative form,

\begin{array}{l} ε_{n} (t_{m}) \leq \frac{K h^{(n + 1) α}}{h^{α} L} (e^{m h^{α} L} - 1) & = \frac{K h^{n α}}{L} (e^{m h^{α} L} - 1) \\ \approx O (h^{(n + 1) α}) . \end{array}

(35)

Comparing this with the local truncation error (LTE), we see that the global error (GE) has the same order as the local error with a different coefficient in the estimates. In the present calculations, the LTE and GE is tabulated according to the following formula: Percent truncation error is

P T R = | \frac{N_{n + 1} (t_{m}) - N_{n} (t_{m})}{N_{n + 1} (t_{m})} | \times 100 % .

(36)

Global truncation error is

G T E = \sum_{k = 0}^{m} | N_{n + 1} (t_{k}) - N_{n} (t_{k}) | .

(37)

5 Numerical Results and Discussions

The validity of the generalized Taylor’s formula for solving the nonlinear fractional point reactor kinetics equations in the presence of Newtonian temperature reactivity feedback and six-groups of delayed neutrons is verified for different types of reactivity. Two cases are considered for the external reactivity, such as step and ramp reactivities. The neutron density is calculated for an U²³⁵ nuclear reactor with various external reactivity. The parameters of this reactor [19] are λi: 0.0124, 0.0305, 0.111, 0.301, 1.13, 3.0 (s-α), β: 0.00021, 0.00141, 0.00127, 0.00255, 0.00074, 0.00027, β=0.00645, l=5.0× 10^-5 (sα), the temperature coefficient of reactivity is σ=5×10^-5 (K^-1), and the reciprocal of the thermal capacity of the reactor is K_c=0.05 (K/MW sα). The initial conditions are N(0)=1 (neutron) and $C_{i} (0) = \frac{β_{i} N (0)}{Λ λ_{i}}$ .

The neutron density is mainly affected by two parameters, the reactivity of the Newtonian temperature feedback and the fractional order, α, as shown in Fig. 1. The first parameter causes a decrease in the neutron density in the absence of an external neutron source, which also agrees with the quicker decreasing of the fractional order (i.e. α< 1). The pattern of this figure also shows a slow decrease for the neutron density with increasing of the fractional order α>1. Moreover, the rate of change of the neutron density is affected also by the fractional order.

Figure 1:

(Color online) Neutron density without external reactivity (γ=0)

5.1 Step external reactivity

The validity and accuracy of the nonlinear fractional model and its proposed solution technique is demonstrated through several different computational examples of the U²³⁵ nuclear reactor for step external reactivity. The neutron density with various times is shown in Figs. 2, 3, 4, 5 and 6 using various values of the step external reactivity as γ=0.25β, 0.5β, 0.75β, 1.0β and 1.25β respectively. The figures show the effect of the positive external reactivity, which have an influence on the neutron density at different fractional orders of α=0.98, 0.99, 1.0, 1.01, and 1.02. The neutron density increases with time to the maximum value, followed by a slow decrease until it is vanished due to the effect of the Newtonian reactivity temperature feedback. Moreover, the effects of sub-diffusion and super-diffusion processes can be classified as follows: the sub-diffusion processes, α< 1, generate a greater number of fission processes, which means the neutron density increases quicker. While the super-diffusion processes, α>1, generate a small number of fission processes, which means slower increases of the neutron density.

Figure 2:

(Color online) Neutron density for step external reactivity (γ=0.25β)

Figure 3:

(Color online) Neutron density for step external reactivity (γ=0.5β)

Figure 4:

(Color online) Neutron density for step external reactivity (γ=0.75β)

Figure 5:

(Color online) Neutron density for step external reactivity (γ=1.0β)

Figure 6:

(Color online) Neutron density for step external reactivity (γ=1.25β)

5.2 Ramp external reactivity

The presented example is based on the same previous parameters of the U²³⁵ nuclear reactor except for the external reactivity and the nonlinear coefficient, σ. In this case, the external reactivity is a function of time such as γ(t)=0.01 tα, 0.01tα, or 0.001tα. Also, the nonlinear coefficient takes the value of σ=2×10^-10 and 2×10^-12.

The neutron density at various times using the generalized Taylor’s formula is calculated at distinct values of the ramp external reactivity function, ρ_ex(t)=0.01tα, and ρ_ex(t)=0.1tα. Figures 7, 8, 9, 10, 11 and 12 show the behavior of the increasing neutron density with time up to the maximum value of the neutron density due to the increase of external reactivity in different fractional order, α=0.98, 1.0, and 1.02. Thereafter, the neutron density decreases due to the effect of the Newtonian temperature reactivity feedback. After that, the external reactivity increasing is equivalent to the feedback effect of Newtonian temperature. Therefore, the neutron density is almost stable. Furthermore, note that the sub-diffusion processes, α< 1, generate a greater number of fission processes, which means that the neutron density increases quicker. On the other hand, the super-diffusion processes, α>1, generates a smaller number of fission processes, which means that the neutron density increases slower.

Figure 7:

(Color online) Neutron density for ramp external reactivity γ(t)=0.01 tα and σ=2×10^-10

Figure 8:

(Color online) Neutron density for ramp external reactivity γ(t)=0.01 tα and σ=2×10^-12

Figure 9:

(Color online) Neutron density for ramp external reactivity γ(t)=0.01 tα and σ=2×10^-10

Figure 10:

(Color online) Neutron density for ramp external reactivity γ(t)=0.01 tα and σ=2×10^-12

Figure 11:

(Color online) Neutron density for ramp external reactivity γ(t)=0.003 tα and σ=2×10^-10

Figure 12:

(Color online) Neutron density for ramp external reactivity γ(t)=0.003 tα and σ=2×10^-12

Tables 1 and 2 show the approximate percent truncation error and global error of the neutron density for step external reactivity insertion at ρ_ex=0.25 $a n d < i > ρ < / i > < s u b > e x < / s u b > = 1.0$ . Similar behavior for local truncation and global errors is presented in Tables 3 and 4 for ramp external reactivity insertion γ(t)=0.001 tα, =2.0E-10, and γ(t)=0.01 tα, and =2.0E-12. In addition, further studies were carried out in Tables 1 and 2 for step external reactivity insertion ρ_ex=0.25 $a n d < i > ρ < / i > < s u b > e x < / s u b > = 1.0$ and, while Tables 3 and 4 for ramp external reactivity insertion determine if improvements in accuracy could be obtained by taking more terms in the series. Excellent results are obtained at the six terms in the series expansion for all representative cases under consideration. In all cases, the generalized Taylor’s formula (GTF) that involves the Caputo fractional derivative technique appears to be accurate and generates results that are consistent with other recently published methods such as PCA [46], EPCA [47], CATS [48], and ITS2 [49], even only by considering a first few order expansion, Table 5 for step reactivity insertion and Table 6 for ramp reactivity insertion.

Approximate percent truncation error and global error of the neutron density for step external reactivity insertion p ex =0.25$

α	t (s)	n=2	n=3	n=4	n=5	n=6	GTE
0.98	50.0	1.55×10^-6	3.36×10^-8	7.03×10^-9	4.37×10^-10	2.75×10^-10	3.18×10^-4
	100.0	3.74×10^-7	1.52×10^-8	2.97×10^-9	9.07×10^-11	1.73×10^-10	3.19×10^-4
	150.0	1.34×10^-6	2.65×10^-8	5.28×10^-9	2.35×10^-10	2.64×10^-10	3.19×10^-4
	200.0	1.45×10^-6	2.95×10^-8	6.03×10^-9	3.49×10^-10	2.73×10^-10	3.19×10^-4
	250.0	1.32×10^-6	3.05×10^-8	6.35×10^-9	4.23×10^-10	2.62×10^-10	3.19×10^-4
1.0	50.0	1.09×10^-6	2.78×10^-8	5.04×10^-9	3.22×10^-10	1.25×10^-10	1.40×10^-4
	100.0	3.40×10^-8	7.56×10^-9	1.38×10^-9	8.80×10^-11	3.41×10^-11	1.40×10^-4
	150.0	9.38×10^-7	1.72×10^-8	3.12×10^-9	1.99×10^-10	7.69×10^-11	1.40×10^-4
	200.0	1.12×10^-6	1.93×10^-8	3.49×10^-9	2.24×10^-10	8.60×10^-11	1.40×10^-4
	250.0	1.08×10^-6	1.99×10^-8	3.59×10^-9	2.30×10^-10	8.85×10^-11	1.40×10^-4
1.02	50.0	7.19×10^-7	2.15×10^-8	3.36×10^-9	2.09×10^-10	4.70×10^-11	6.15×10^-5
	100.0	2.99×10^-7	2.46×10^-9	4.14×10^-10	4.12×10^-11	4.69×10^-12	6.15×10^-5
	150.0	6.07×10^-7	1.08×10^-8	1.75×10^-9	1.48×10^-10	1.29×10^-11	6.15×10^-5
	200.0	8.46×10^-7	1.26×10^-8	2.01×10^-9	1.62×10^-10	1.65×10^-11	6.15×10^-5
	250.0	8.60×10^-7	1.30×10^-8	2.03×10^-9	1.53×10^-10	1.90×10^-11	6.15×10^-5

Approximate percent truncation error and global error of the neutron density for step external reactivity insertion p ex = 1.0$

α	t (s)	n=2	n=3	n=4	n=5	n=6	GTE
0.98	1.0	1.21×10^-2	1.52×10^-4	2.19×10^-7	4.80×10^-8	1.05×10^-9	1.41×10^-6
	2.0	7.62×10^-3	8.93×10^-5	5.63×10^-7	2.32×10^-8	1.00×10^-9	1.63×10^-6
	3.0	5.47×10^-3	6.78×10^-5	5.11×10^-7	3.69×10^-8	3.72×10^-9	2.60×10^-6
	4.0	4.44×10^-3	5.62×10^-5	4.92×10^-7	4.95×10^-8	6.93×10^-9	4.36×10^-6
	5.0	3.84×10^-3	4.93×10^-5	4.83×10^-7	5.99×10^-8	1.02×10^-8	6.66×10^-6
1.0	1.0	5.44×10^-3	5.64×10^-7	2.47×10^-6	3.65×10^-8	8.03×10^-10	6.52×10^-7
	2.0	6.53×10^-3	6.57×10^-5	2.74×10^-7	5.21×10^-9	9.27×10^-11	7.55×10^-7
	3.0	4.63×10^-3	5.03×10^-5	1.98×10^-7	3.89×10^-9	6.82×10^-11	7.92×10^-7
	4.0	3.73×10^-3	4.16×10^-5	1.60×10^-7	3.20×10^-9	5.55×10^-11	8.15×10^-7
	5.0	3.21×10^-3	3.63×10^-5	1.39×10^-7	2.78×10^-9	4.80×10^-11	8.30×10^-7
1.02	1.0	1.20×10^-3	1.66×10^-4	4.58×10^-7	4.77×10^-8	6.68×10^-10	2.55×10^-7
	2.0	5.61×10^-3	4.77×10^-5	1.52×10^-7	1.40×10^-9	1.17×10^-10	3.72×10^-7
	3.0	3.93×10^-3	3.71×10^-5	8.45×10^-8	1.70×10^-9	3.41×10^-10	4.79×10^-7
	4.0	3.15×10^-3	3.06×10^-5	5.17×10^-8	4.19×10^-9	6.38×10^-10	6.67×10^-7
	5.0	2.69×10^-3	2.66×10^-5	3.22×10^-8	6.18×10^-9	9.55×10^-10	9.21×10^-7

Approximate percent truncation error and global error of the neutron density for ramp external reactivity insertion γ(t) = 0.001tα and σ= 2.0×10^-10

α	t (s)	n=4	n=5	n=6	n=7	GTE
0.98	1.0	9.63×10^-12	3.00×10^-13	5.29×10^-14	1.76×10^-14	3.73×10^-14
	2.0	9.01×10^-12	2.47×10^-13	3.70×10^-14	0.0	3.73×10^-14
	3.0	8.81×10^-12	2.54×10^-13	2.82×10^-14	1.41×10^-14	1.13×10^-3
	4.0	6.43×10^-14	2.79×10^-13	6.43×10^-14	4.29×10^-14	1.35×10^-3
	5.0	5.08×10^-10	2.99×10^-13	3.99×10^-14	0.0	1.41×10^-3
	6.0	2.16×10^-7	6.66×10^-10	1.17×10^-12	0.0	1.44×10^-3
1.0	1.0	3.40×10^-12	1.46×10^-13	5.49×10^-14	0.0	3.26×10^-14
	2.0	2.82×10^-12	1.10×10^-13	2.75×10^-14	0.0	3.26×10^-14
	3.0	2.54×10^-12	1.09×10^-13	3.63×10^-14	0.0	2.66×10^-4
	4.0	2.37×10^-12	1.15×10^-13	3.83×10^-14	0.0	4.07×10^-4
	5.0	1.18×10^-11	1.13×10^-13	5.02×10^-14	0.0	4.59×10^-4
	6.0	8.60×10^-10	4.08×10^-13	1.63×10^-14	0.0	5.11×10^-4
1.02	1.0	1.17×10^-12	7.54×10^-14	1.89×10^-14	0.0	0.0
	2.0	7.63×10^-13	5.99×10^-14	4.49×10^-14	0.0	0.0
	3.0	6.10×10^-13	8.71×10^-14	4.35×10^-14	0.0	0.0
	4.0	4.71×10^-13	5.54×10^-14	2.77×10^-14	0.0	0.0
	5.0	9.79×10^-14	8.39×10^-14	5.59×10^-14	0.0	0.0
	6.0	1.47×10^-11	8.67×10^-14	5.20×10^-14	0.0	0.0

Approximate percent truncation error and global error of the neutron density for ramp external reactivity insertion γ(t) = 0.001t^αand σ= 2.0×10^-12

α	t (s)	n=6	n=7	n=8	n=9	n=10	GTE
0.98	0.5	3.43×10^-13	3.67×10^-14	0.0	0.0	0.0	0.0
	1.0	3.41×10^-8	7.99×10^-9	6.58×10^-12	1.12×10^-11	2.28×10^-13	5.95×10^-2
	1.5	5.00×10^-9	3.28×10^-11	1.65×10^-13	1.24×10^-13	0.0	1.69×10^-1
	2.0	1.27×10^-9	2.22×10^-12	1.48×10^-14	0.0	0.0	1.76×10^-1
	2.5	3.41×10^-11	1.19×10^-12	2.89×10^-14	1.45×10^-14	0.0	1.79×10^-1
	3.0	3.13×10^-11	6.28×10^-13	1.46×10^-14	1.46×10^-14	2.92×10^-14	1.84×10^-1
1.0	0.5	3.33×10^-13	6.24×10^-14	0.0	0.0	0.0	0.0
	1.0	1.97×10^-8	1.97×10^-10	1.61×10^-12	0.0	0.0	0.0
	1.5	1.31×10^-9	1.21×10^-11	7.52×10^-14	7.52×10^-14	0.0	4.86×10^-2
	2.0	1.24×10^-9	4.67×10^-12	6.58×10^-14	1.64×10^-14	0.0	5.45×10^-2
	2.5	1.67×10^-10	2.88×10^-14	1.44×10^-14	1.44×10^-14	0.0	5.64×10^-2
	3.0	9.34×10^-13	3.59×10^-13	1.44×10^-14	0.0	1.44×10^-14	5.72×10^-2
1.02	0.5	8.75×10^-14	0.0	0.0	0.0	0.0	0.0
	1.0	4.81×10^-10	2.92×10^-12	0.0	0.0	0.0	0.0
	1.5	2.12×10^-10	3.55×10^-12	1.55×10^-13	5.64×10^-14	0.0	0.0
	2.0	5.71×10^-10	2.85×10^-12	1.21×10^-13	4.04×10^-14	0.0	0.0
	2.5	2.48×10^-10	3.91×10^-13	3.01×10^-14	1.51×10^-14	0.0	0.0
	3.0	3.13×10^-11	8.53×10^-14	2.84×10^-14	1.42×10^-14	0.0	0.0

Neutron density for step external reactivity insertion p_ex=1.0$ and fractional order α=1.0

t (s)	GTF(n=6, h=0.01)	PCA	EPCA	CATS	ITS2
0.0	1.0	1.0	1.0	1.0	1.0
10.0	132.03859643	132.038654	132.038597	132.038596	132.03859643
20.0	51.699860947	51.6998767	51.6998611	51.6998609	51.699860947
30.0	28.174685361	28.1746922	28.1746854	28.1746854	28.174685361
40.0	18.146329995	18.1463339	18.1463300	18.1463300	18.146329995
50.0	12.779577032	12.7795795	12.7795770	12.7795770	12.779577032
60.0	9.4749325010	9.47493414	9.47493251	9.47493250	9.4749325010
70.0	7.2444774938	7.24447861	7.24447750	7.24447749	7.2444774937
80.0	5.6462897002	5.64629045	5.64628971	5.64628970	5.6462897002
90.0	4.4568342546	4.45683475	4.45683426	4.45683425	4.4568342546
100.0	3.5501027660	3.55010308	3.55010277	3.55010277	3.5501027660

Neutron density for ramp external reactivity insertion σ=2.0×10^-12and fractional order α=1.0

γ(t)	t (s)	GTF(n=10, h=0.001)	EPCA	CATS	ITS2
0.1 t	0.1	2.4733658251×10,⁺¹	2.473365830×10,⁺¹	2.473365825×10,⁺¹	2.4733658251×10,⁺¹
	0.5	1.5433617863×10⁺¹²	1.543361754×10⁺¹²	1.543361786×10⁺¹²	1.5433617863×10⁺¹²
	5.0	1.0029740921×10⁺¹²	1.002974092×10⁺¹²	1.002974092×10⁺¹²	1.0029740921×10⁺¹²
	7.5	1.0017984372×10⁺¹²	1.001798437×10⁺¹²	1.001798437×10⁺¹²	1.0017984372×10⁺¹²
	10.0	1.0011886207×10⁺¹²	1.001188621×10⁺¹²	1.001188621×10⁺¹²	1.0011886207×10⁺¹²
0.01 t	0.1	1.1672108379	1.167210838	1.167210838	1.1672108379
	0.5	4.2699528644	4.269952865	4.269952864	4.2699528644
	5.0	1.0338896655×10⁺¹¹	1.033889665×10⁺¹¹	1.033889665×10⁺¹¹	1.0338896655×10⁺¹¹
	7.5	1.0194999125×10⁺¹¹	1.019499913×10⁺¹¹	1.019499913×10⁺¹¹	1.0194999125×10⁺¹¹
	10.0	1.0124348832×10⁺¹¹	1.012434883×10⁺¹¹	1.012434883×10⁺¹¹	1.0124348832×10⁺¹¹
0.003 t	0.1	1.0453716665	1.045371667	1.045371666	1.0453716665
	0.5	1.3246619862	1.324661986	1.324661986	1.3246619862
	5.0	3.2156761131×10⁺¹⁰	3.215676113×10⁺¹⁰	3.215676113×10⁺¹⁰	3.2156761131×10⁺¹⁰
	7.5	3.2102051821×10⁺¹⁰	3.210205182×10⁺¹⁰	3.210205182×10⁺¹⁰	3.2102051821×10⁺¹⁰
	10.0	3.1456146867×10⁺¹⁰	3.145614687×10⁺¹⁰	3.145614687×10⁺¹⁰	3.1456146867×10⁺¹⁰

6 Conclusion

The objective of this study was to estimate the GTF technique as applied to point reactor kinetics problems with temperature feedback of reactivity for the multi-group of the precursor delayed neutrons. The method is based on a generalized Taylor’s formula expansions in the fractional form of the neutron density and reactivity functions to obtain explicit analytical solutions to the reactor point kinetics equations in the integral formulation. Variable time steps are used by the method to alleviate stiffness arising from the orders of magnitude with six delayed neutron groups. Furthermore, the transient neutron density of the nonlinear fractional model using the proposed technique was computed for different cases of external reactivities and various fractional orders. The obtained results were found in excellent agreement at the fractional order one with some of the conventional methods, in all cases under consideration. It can be seen the GTF is in excellent agreement with EPCA and, particularly, with CATS and ITS2, therefore also confirming the extreme accuracy of the method for these test cases. Results by GTF are believed accurate to all digits. The attached figures show the pattern of the neutron density, which increase with time to the maximum value due to the external reactivity followed by decreasing with time due to the effect of the Newtonian temperature feedback reactivity. In addition, for a positive reactivity, we note that: the sub-diffusion processes, α< 1, generate a greater number of fission processes and the super-diffusion processes, α>1, generate a smaller number of fission processes. On the other hand, for a negative reactivity, the neutron density decreases quicker with decreasing the fractional order, α< 1 and the neutron density decreases slower for increasing the fractional order, α>1. In the future work, the GTF will be applied to the two energy groups for point kinetics equations with emphasis on clarifying the assumption used in applying the method, like the values of α close to one are considered, and the same power is considered for the neutron and precursor equations.

References

[1]

M. Caputo, Elasticità Dissipazione, (Zanichelli, Bologna, Italy, 1969).

[2]

K.B. Oldham, J. Spanier, The Fractional Calculus, (Academic Press, New York, 1974).

[3]

I. Podlubny, Fractional Differential Equations, (Academic Press, San Diego, CA 1999).

[4]

K. Diethelm, J. Ford,

Analysis of fractional differential equations

, J. Math. Anal. Appl. 265, 229-248 (2002). doi: 10.1006/jmaa.2000.7194