1. Introduction
The revolution of using energy generated from nuclear reactors can be considered as an essential part of a clean renewable source of energy; today, nuclear reactors cover 10% of the use of the global needed energy.
The inspiring idea is that the creation of nuclear energy is released when heavy nuclei undergo fission while they are bombarded by (thermal) slow neutrons. The neutron diffusion equation, which is simplified from the neutron transport equation when using Fick’s law, describes the behavior of neutrons in the nuclear fission reactor [
1,
2,
3]. Researchers have investigated their scientific ambitions by studying different issues in this field, where different cases of neutron diffusion equations have been solved by various techniques. The homotopy perturbation method (HPM) has been used in solving bare hemispherical symmetry and cylindrical symmetry one group neutron diffusion equation reactors [
4,
5]; the residual power series method was used to solve multi-energy groups of neutron diffusion equations [
6,
7]; then the HPM is used to solve reflected reactor one group different geometries neutron diffusion equations; finally, solving a multi-group reflected spherical reactor system of equations using the HPM is performed [
8].
The distribution of neutrons in nuclear reactors is constructed in the same way as dissolving the solute in the chemical solutions, which can be represented mathematically by the continuity equation:
where
is the neutron flux;
is neutron current density;
is the macroscopic absorption cross-section;
is neutron velocity; and
is a source of neutrons term. Using Fick’s law, which gives the relation between neutron flux
and the current density
where
D is the neutron diffusion coefficient, then the neutron diffusion equation can be expressed as
This is the well-known nuclear reactor time-dependent neutron diffusion equation.
One step forward, the fractional study of the time-dependent neutron diffusion equation is to make the situation more detailed by studying the physical effect of delayed neutrons. These delayed neutrons have an essential role in the reactor. Investigating a system of coupled fractional neutron diffusion equations with delayed neutrons was the target of this work.
The importance of delayed neutrons in reactors is the effective lifetime of such neutrons (which are given by their prompt lifetime, and additional delayed time), which is characterized by the β-decay of their precursor, which is considerably longer than the prompt neutron lifetime. Hence, delayed neutrons substantially increase the time constant of a reactor so that effective control is possible.
The contribution of the effect of the delayed neutron is included in Equation (3) in the source term
, which can be written as
In this study, the general three-dimension vector will be simplified to the one-dimension vector to be more suitable.
The system of time-dependent neutron diffusion equations is known as the neutron point kinetic model; this is in one dimension:
with initial conditions:
where the delayed neutron density is
, and the constants γ, λ, and β are the average number of neutrons produced per fission, the radioactive decay constant, and the fraction of the delayed fission neutrons, respectively, where
can be considered as the total macroscopic cross-section
.
The kinetic equation describes the dynamics of space-time neutron distribution in the reactors. These equations are divided into point kinetic equations and space kinetic equations; the time change in neutron density is described by the kinetic point equations, which estimate the reactor power response; the fractional calculus governs an important part of the time-dependent diffusion equation in the presence of delayed neutrons, which is known as anomalous diffusion; the anomalous diffusion is called by the variation from the Gaussian distribution of the diffusion.
Anomalous diffusion is characterized by the asymptotical nonlinear (especially power-law type) dependence of the mean-squared displacement on the time, which is .
The parameter α feeds back in the determination of the type of fractional neutron diffusion equation, where the anomalous diffusion ( is classified into sub diffusion () and super diffusion ().
In the case of α = 1: the diffusion is Gaussian, mathematically representing the non-fractional neutron diffusion equation.
The trapping of neutrons in the medium causes slower neutron movement considered as sub diffusion, while the fast movement of neutrons when they move in one direction for a very long time without any collision is characterized by super diffusion [
9,
10].
The fractional neutron diffusion equation with a delayed neutron system has been studied in many works where the fractional two energy group matrix representation for nuclear reactor dynamics with an external source is considered [
11], while the coupled fractional neutron diffusion equations with delayed neutrons have been solved [
12], and an exact solution in the case of the one-dimensional neutron diffusion kinetic equation with one delayed precursor concentration in Cartesian geometry was studied [
13].
Here, the fractional neutron diffusion equations with one delayed neutron group can be written as follows:
Numerous numerical and analytical techniques including Laplace transforms, the homotopy method, the Pade-variational iteration method, the finite difference method, the Adomian decomposition method, and others [
14,
15,
16,
17,
18,
19,
20,
21,
22,
23,
24] have been developed to solve fractional differential equations. However, these methods have limitations such as the need for long and complex computations. Therefore, they need to be connected to a transform operator. The authors of [
25] combined the residual power series method (RPSM) with the Laplace transform to introduce actual and series solutions of the linear and nonlinear FDEs. The new method, known as the Laplace-residual power series method (LRPSM), was adopted to construct series solutions of a variety of fractional deferential equations (FDEs). LRPSM does not rely on fractional derivation to determine the coefficients of the series as in RPSM, but depends on the concept of the limit, so few calculations generate the coefficients compared to the residual power series method. The current technique is quick, requires little computer memory, and is not influenced by computational round-off errors. Furthermore, this technique computes the coefficients of the power series using a chain of equations with more than one variable, indicating that the present method has a rapid convergence.
This paper solved a system of coupled fractional neutron diffusion equations with delayed neutrons efficiently by applying LRPSM. The accuracy and effectiveness of the method were clarified by displaying numerical examples and comparing the solutions with the results of some approved approaches.
The rest of paper is structured as follows. The concepts and theorems in fractional calculus and Laplace transforms that are required to reach our results are introduced in
Section 2. The LRPSM algorithm is described in
Section 3.
Section 4 introduces the numerical examples and compares the solutions with the results of some approved approaches. In
Section 5, we finally present our conclusions.
2. Fractional Expansion in Laplace Space
The fractional derivatives are non-local operators because they are defined using integrals. Therefore, the time-fractional derivative contains information about the function at earlier points, thus it possesses a memory effect. Such derivatives consider the history and non-local distributed effects, which are essential for better and more accurate descriptions and an understanding of complex and dynamic system behavior.
The fractional calculus in Caputo’s sense and the Laplace transform, which are necessary to create the LRPS solution for the fractional neutron diffusion equations with one delayed neutron group, are introduced in this part with some basic concepts and ideas.
Definition 1. The time-fractional Caputo derivative of order of sense is defined aswhere denotes the time-fractional Caputo derivative operator of order ;
; is an interval; and is the time-fractional Riemann–Liouville integral operator of order that is defined as: In the following lemma, we introduce certain properties of the operator
that we will require for the rest of the work in this paper. The reader should see [
26] for more features.
Lemma 1. For,, and, we have: Definition 2. Letbe a continuous function onand of exponential order. Then, the Laplace transform of the functionis denoted and defined as follows:whereas the inverse Laplace transform of the functionis defined as follows:wherelies in the right half plane of the absolute convergence of the Laplace integral. In the next lemma, we introduce several crucial characteristics related to the Laplace transform and the fractional derivative in the Caputo sense.
Lemma 2. Let be a continuous function onand of exponential orders, and . Then,where (-times) Now, as shown in the following result, we add a new fractional expansion that serves as the foundation for creating a LRPS solution for the partial differential equations (PDEs).
Theorem 1. Let be continuous on and of exponential order . Suppose that the functionhas the following fractional expansion:then .
The following theorem shows the conditions that guarantee the convergence of the series in the expansion (17).
Theorem 2. Let be a continuous function on and of exponential order and can be represented as the fractional expansion in Theorem 1. If on where , then the remainder of the fractional expansion as in Equation (17) satisfies the following inequality:.
3. Exact Solution of the Fractional Neutron Diffusion Equations with One Delayed Neutron Group
The main objective of this section is to implement the LRPSM to construct a series solution to the fractional neutron diffusion equations with one delayed neutron group
with the initial conditions
Rewriting the coupled equations in (18) by replacing
and
we obtain
with the initial conditions
Operating the Laplace transform on the coupled equations in (20), we obtain
According to Lemma 2 part (iii) and using the initial conditions in Equation (21), we can write equations in (22) in the following form:
where
Equation (23) can be written as
The coupled equations in Equation (24) represent a linear system of partial differential equations that contain derivatives relative to the independent variable
. According to the Laplace residual power series method, we assume that the series solution of the system in Equation (24) has the form
Using the fact in Lemma 2 part (i), of the
kth truncated series of
can be expressed as
In the next step, we define the Laplace residual functions of the coupled equations in (24) in order to find the unknown coefficients of the series in Equation (26):
and the
kth Laplace-residual functions as
Since
and
, we have
. Therefore,
To find
and
in Equation (26), we substituted
and
in the first Laplace-residual functions to obtain
Next, by solving
,
one can obtain
To find
and
in Equation (26), we substituted
and
in the second Laplace-residual functions to obtain
Next, by solving
,
one can obtain
In the third Laplace, the residual functions to obtain
By solving
,
one obtains
If we continue in the same manner, we can substitute the k
th truncated series
,
into the k
th Laplace-residual function
. Then, we can multiply the resulting equations by
and taking the limit as
,
,
for
can be obtained by the following recurrence relation:
where
.
According to what was introduced, the series solution of system (24) is
Therefore, the series solution of system (18) can be obtained by transforming the solution in (41) into the original space by using the inverse Laplace transform. Therefore, the Laplace residual power series solution of system (18) is given by:
where
4. Numerical Results and Discussion
To verify the driven theory, we solved the time-dependent fractional neutron diffusion equations with delayed neutrons using the Laplace residual power series method, and the results were tested with the numerical values of the following nuclear reactor cross-section data [
10]
It must be noted that these cross-sections are defined in
Section 1.
The values of were obtained using MATHEMATICA in 11th Gen Intel(R) Core (TM) i7-1165G7 2.80 GHz with a 16 GB RAM computer. Once are calculated, the data to construct the tables and figures are achieved.
The numerical results of neutron flux at different values of
are shown in
Table 1 for
.
The flux was calculated at different times, and the anomalous diffusion,
(sub diffusion) was considered. Additionally, it is clear from
Table 1 that the neutron flux increased with time, as shown in
Figure 1.
For more clarification, the neutron fluxes
are plotted in
Figure 1 with the initial condition
.
Now, the different initial value problem
was used in the study of the following case, where the results are clarified in
Figure 2.
The plotted profile of the neutron fluxes for different values of
(
0.3, 0.4 and
) in both cases
and
, as presented in
Table 1 and
Figure 1 and
Figure 2, shows that the results were in agreement with the Adomian decomposition method [
10].
As observed from the tables and figures, the Laplace residual power series method used to solve the time-dependent fractional neutron diffusion equations with delayed neutrons provided us with more details about the neutron flux in non-Gaussian diffusion for different selected times.
The importance of this work was manifested by uncovering the behavior of the flux in the presence of the delayed neutrons that affected them (this phenomenon is known as anomalous diffusion). This flux can be expressed mathematically by the solution of the time-dependent fractional neutron diffusion equations, as shown in Equation (42).
It must be noted that the time-dependent flux increased with time in both the studied numerical cases [
9,
10,
12].