Analysis of Fractional-order Point Reactor Kinetics Model with Adiabatic Temperature Feedback for Nuclear Reactor with Subdiffusive Neutron Transport

Abstract

In this paper, a fractional-order nonlinear model is developed for the nuclear reactor with subdiffusive neutron transport. The proposed fractional-order point reactor kinetics model is a system of three coupled, nonlinear differential equations. The model represents subprompt critical condition. The nonlinearity in the model is due to the adiabatic temperature feedback of reactivity. This model originates from the fact that neutron transport inside the reactor core is subdiffusion and should be better modeled using fractional-order differential equations. The proposed fractional-order model is analyzed for step and sinusoidal reactivity inputs. The stiff system of differential equations is solved numerically with Adams-Bashforth-Moulton method. The proposed model is stable with self-limitting power excursions. The issue of convergence of this method for the proposed model for different values of fractional derivative order is also discussed.

Appendix: Fractional Second-order Adams-Bashforth-Moulton Method

Here the fractional second-order Adams-Bashforth-Moulton (ABM) method which is used in Sect. 4 is explained in brief. The main computational steps involved in the algorithm are presented here for the equispaced grid points. For details, refer to [33, 37, 38]. It is an extension of the classical ABM method used to numerically solve the first-order ODEs. It comes in the category of the so-called PECE (Predict, Evaluate, Correct, Evaluate) type since it involves calculation of the predictor value which is in turn used to compute the corrector value. This method and its variants are popular in the field of fractional calculus and applied areas [39, 40]. The algorithm explained below is for a single fractional differential equation. However, it can be easily modified to handle a system of FDEs.

Consider the single term FDE with Caputo FD

$$\begin{aligned} _{0}D^{\alpha }_{t}y(t) = f(t,y(t)), \end{aligned}$$

(21)

where $\alpha \in \mathbb {R}^+$ and with the appropriate initial conditions:

$$\begin{aligned} D^k_{t}y(0) = y_0^{(k)}, \quad k = 0, 1, \ldots , m-1, \end{aligned}$$

(22)

where, $m=\lceil \alpha \rceil$ is the ceil function. The equivalent Volterra integral equation is

$$\begin{aligned} y(t) = \sum _{k=0}^{m-1}\frac{t^k}{k!} D^k_{t} y(0) + \frac{1}{\varGamma (\alpha )} \int _0^t (t-\tau )^{\alpha -1} f(\tau ,y(\tau )) d\tau . \end{aligned}$$

(23)

The integration limits from 0 to t imply the nonlocal structure of the fractional derivatives.

The next step is to use the product trapezoidal quadrature formula to replace the integral in (23). We approximate the following integral

$$\begin{aligned} \int _0^{t_{k+1}} (t_{k+1}-\tau )^{\alpha -1} g(\tau ) d\tau , \end{aligned}$$

(24)

as

$$\begin{aligned} \int _0^{t_{k+1}} (t_{k+1}-\tau )^{\alpha -1} g_{k+1}(\tau ) d\tau , \end{aligned}$$

(25)

where $\tilde{g}_{k+1} \equiv$ piecewise linear interpolation for g(t) with grid points at $t_j,$ $j = 0, 1, 2, \ldots , k+1$. Thus we can write the integral (25) as

$$\begin{aligned} \int _0^{t_{k+1}} (t_{k+1}-\tau )^{\alpha -1} g_{k+1}(\tau ) d\tau = \sum _{j=0}^{k+1} a_{j, k+1}g(t_j), \end{aligned}$$

(26)

for the equispaced nodes ($t_j = jh$ with some fixed step-size h). The values of $a_{j, k+1}$ are given for $j=0$ as

$$\begin{aligned} \frac{h^\alpha }{\alpha (\alpha +1)} \left( k^{\alpha +1} - (k-\alpha )(k+1)^\alpha \right) , \end{aligned}$$

for $1 \le j \le k$ as

$$\begin{aligned} \left( \frac{h^\alpha }{\alpha (\alpha +1)}\right) (d), \end{aligned}$$

where

$$\begin{aligned} d = (k-j+2)^{\alpha +1} + (k-j)^{\alpha +1} - 2(k-j+1)^{\alpha +1}, \end{aligned}$$

and for $j = k + 1$ as

$$\begin{aligned} \frac{h^\alpha }{\alpha (\alpha +1)}. \end{aligned}$$

So the corrector formula is

$$\begin{aligned} y_{k+1} =&\sum _{j=0}^{m-1}\frac{t_{k+1}^j}{j!} y_0^{(j)} \nonumber \\&+ \frac{1}{\varGamma (\alpha )} \sum _{j=0}^k a_{j, k+1} f(t_j,y_j) \\&+ \frac{1}{\varGamma (\alpha )} \left( a_{k+1,k+1} f(t_{k+1},y^P_{k+1})\right) \nonumber , \end{aligned}$$

(27)

where now the predictor $y^P_{k+1}$ is evaluated as

$$\begin{aligned} y^P_{k+1} = \sum _{j=0}^{m-1}\frac{t_{k+1}^j}{j!} y_0^{(j)} + \frac{1}{\varGamma (\alpha )} \sum _{j=0}^k b_{j, k+1} f(t_j,y_j), \end{aligned}$$

(28)

with

$$\begin{aligned} b_{j, k+1} = \frac{h^\alpha }{\alpha } \left( (k+1-j)^\alpha - (k-j)^\alpha \right) . \end{aligned}$$

(29)

For $0<\alpha <1$, the predictor and corrector expressions get modified as

$$\begin{aligned} y^P_{k+1} = y_0 + \frac{1}{\varGamma (\alpha )} \sum _{j=0}^k b_{j, k+1} f(t_j,y_j), \end{aligned}$$

(30)

and

$$\begin{aligned} y_{k+1} =&\, y_0 + \frac{1}{\varGamma (\alpha )} \sum _{j=0}^k a_{j, k+1} f(t_j,y_j) \nonumber \\&+ \frac{1}{\varGamma (\alpha )} \left( a_{k+1,k+1} f(t_{k+1},y^P_{k+1})\right) . \end{aligned}$$

(31)

As already mentioned in Sect. 4, the convergence of this algorithm deteriorates as $\alpha \rightarrow 0$. This algorithm was coded in MATLAB.