Design and analysis of a numerical method for fractional neutron diffusion equation with delayed neutrons

Abstract

The main purpose of this work is to construct and analyze an efficient numerical scheme for solving the fractional neutron diffusion equation with delayed neutrons, which describes neutron transport in a nuclear reactor. The L1 approximation is used for discretization of time derivative and finite difference method is used for discretization of space derivative. The stability and convergence analysis of the proposed method are studied. The method is shown to be second-order convergent in space and $(2-2\alpha )$-th order convergent in time, where α is the order of fractional derivative. Numerical experiments are carried out to demonstrate the performance of the method and theoretical analysis. The effects of fractional order derivative, relaxation time and radioactive decay constant on the neutron flux behaviour are investigated. Moreover, the CPU time of the present method is provided.

Introduction

The neutron population distribution in a nuclear reactor is described by transport equation. One of the simplest approximations to neutron transport is the approximation given by the neutron diffusion theory. The diffusion equation describes the behaviour of a large amount of neutrons when individual properties of separate neutron trajectories in matter are “averaged”. The one-group neutron diffusion equation in one-dimensional space when delayed neutrons are averaged by one group of delayed neutrons is given by (see [33], [1], [31]):

$$\frac{1}{\nu }\frac{\partial \mathrm{\Phi }(x,t)}{\partial t}=D\frac{{\partial }^2\mathrm{\Phi }(x,t)}{\partial x^2}+(\gamma {\mathrm{\Sigma }}_f-{\mathrm{\Sigma }}_a)\mathrm{\Phi }(x,t)+\lambda C(x,t),(x,t)\in (0,X)\times (0,T),$$

$$\frac{\partial C(x,t)}{\partial t}=\beta \gamma {\mathrm{\Sigma }}_f\mathrm{\Phi }(x,t)-\lambda C(x,t),$$

with initial condition

$$\mathrm{\Phi }(x,0)={\mathrm{\Phi }}^0(x)$$

and boundary conditions

$$\mathrm{\Phi }(0,t)=0,\mathrm{\Phi }(X,t)=0.$$

Here, ν is the neutron velocity, $\mathrm{\Phi }(x,t)$ is the neutron flux, $C(x,t)$ is the concentration neutron precursors, D is the neutron diffusion coefficient, γ is the average number of neutrons produced per fission, ${\mathrm{\Sigma }}_f$ is the macroscopic fission cross-section, ${\mathrm{\Sigma }}_a$ is the macroscopic absorption cross-section, λ is the radioactive decay constant and β is the fraction of delayed neutrons. In [32], Sardar et al. presented a fractional neutron diffusion model with delayed neutrons of the following form

$$\frac{1}{\nu }\frac{{\partial }^{\alpha }\mathrm{\Phi }(x,t)}{\partial t^{\alpha }}=D\frac{{\partial }^2\mathrm{\Phi }(x,t)}{\partial x^2}+(\gamma {\mathrm{\Sigma }}_f-{\mathrm{\Sigma }}_a)\mathrm{\Phi }(x,t)+\lambda C(x,t),0<\alpha <\frac{1}{2},(x,t)\in (0,X)\times (0,T),$$

$$\frac{{\partial }^{\alpha }C(x,t)}{\partial t^{\alpha }}=\beta \gamma {\mathrm{\Sigma }}_f\mathrm{\Phi }(x,t)-\lambda C(x,t),$$

with initial condition

$$\mathrm{\Phi }(x,0)={\mathrm{\Phi }}^0(x)$$

and boundary conditions

$$\mathrm{\Phi }(0,t)=0,\mathrm{\Phi }(X,t)=0.$$

It is worth pointing out that the above time-fractional model is not dimensionally correct. In this work, we present a dimensionally correct fractional neutron diffusion model with delayed neutrons of the following form

$${\mathrm{\tau }}^{\alpha }{\mathrm{\Sigma }}_a\frac{{\partial }^{\alpha }\mathrm{\Phi }(x,t)}{\partial t^{\alpha }}=D\frac{{\partial }^2\mathrm{\Phi }(x,t)}{\partial x^2}+[(1-\beta )\gamma {\mathrm{\Sigma }}_f-{\mathrm{\Sigma }}_a]\mathrm{\Phi }(x,t)+\lambda C(x,t),0<\alpha <\frac{1}{2},(x,t)\in (0,X)\times (0,T),$$

$${\mathrm{\tau }}^{\alpha -1}\frac{{\partial }^{\alpha }C(x,t)}{\partial t^{\alpha }}=\beta \gamma {\mathrm{\Sigma }}_f\mathrm{\Phi }(x,t)-\lambda C(x,t),$$

subject to the initial conditions

$$\mathrm{\Phi }(x,0)={\mathrm{\Phi }}^0(x),C(x,0)=\frac{\beta \gamma {\mathrm{\Sigma }}_f}{\lambda }{\mathrm{\Phi }}^0(x),\frac{{\partial }^{\alpha }\mathrm{\Phi }(x,0)}{\partial t^{\alpha }}=0,\frac{{\partial }^{\alpha }C(x,0)}{\partial t^{\alpha }}=0$$

and boundary conditions

$$\mathrm{\Phi }(0,t)=0,\mathrm{\Phi }(X,t)=0,C(0,t)=0,C(X,t)=0,$$

where α is the anomalous diffusion order and ${\mathrm{\tau }}^{\alpha }$ is the fractional relaxation time which is given by

$${\mathrm{\tau }}^{\alpha }={\left(\frac{3D}{\nu }\right)}^{\alpha }.$$

The derivation of this fractional neutron diffusion model is given in the next section. Fractional differential equations have been found to provide a more realistic representation for anomalous diffusion occurring in nature and the theory of complex systems. For instance, models based on fractional-order differential structures have been used in studying diffusion systems [11], [24], [3], [22], [19], financial systems [28], [34], [18], [23], [7], biological systems [27], [26], [12], [20], [21], electrical systems [17], [35], [2] and heat conduction problem [4]. Fractional order modelling of nuclear reactor has been presented in [15], [13], [14], [16], [30], [10], [29]. In many cases, analytical solutions to fractional differential equations are typically difficult or impossible to be obtained. In fact, analytical solution to the fractional neutron diffusion equation with delayed neutrons defined by equations (9)-(10) is not known. We note that no numerical technique has been proposed in literature to solve the problem (9)-(12). The authors of [32] used Adomian decomposition method to obtain series solution of fractional neutron diffusion model described by (5)-(8). In [31], Sapagovas and Vileiniskis used method of summary approximation for numerical solution of two-dimensional classical neutron diffusion equation with one group of delayed neutron.

The aim of this work is to develop an efficient numerical method for solving the fractional neutron diffusion model (9)-(12) on a uniform grid. In this method, we first transform the coupled fractional differential equations (9)-(10) into a fractional differential equation and then develop a numerical method to tackle the resulting fractional differential equation. More specifically, the time-fractional derivative in the resulting equation is described in the sense of Caputo. The L1 approximation is used in time direction to obtain a backward Euler-type method and finite difference method (FDM) is used for discretization of space derivative. The stability and convergence of the proposed method are investigated. The method is second order accurate with respect to space variable and $(2-2\alpha )$-th order accurate with respect to time variable. Numerical experiments are carried out to demonstrate the performance of the method and to validate the theoretical results. The effects of the order of fractional derivative α, relaxation time τ and radioactive decay constant λ on the behaviour of neutron flux are investigated. The CPU time of the proposed method is provided.

The rest of the paper is organized as follows: In section 2, we derive the fractional neutron diffusion model as described by equations (9)-(10). In section 3, we construct a numerical method based on finite difference method to approximate the solution of the problem considered. Stability of the proposed method is discussed in section 4. Section 5 is devoted to the convergence analysis of the method. In section 6, numerical experiment is carried out to illustrate the applicability of the present numerical scheme. Some concluding remarks are given in section 7.