Open AccessArticle

Solving a Novel System of Time-Dependent Nuclear Reactor Equations of Fractional Order

Doaa Filali

Mohammed Shqair

²,

Fatemah A. Alghamdi

Sherif Ismaeel

^4,5 and

Ahmed Hagag

^6,*

Department of Mathematical Science, College of Sciences, Princess Nourah Bint Abdulrahman University, Riyadh 84428, Saudi Arabia

College of Science, Zarqa University, Zarqa 13110, Jordan

Financial Sciences Department, Applied College, Imam Abdulrahman Bin Faisal University, Dammam 31441, Saudi Arabia

⁴

Department of Physics, Faculty of Science, Prince Sattam Bin Abdulaziz University, Al-Kharj 16278, Saudi Arabia

⁵

Department of Physics, Faculty of Science, Ain Shams University, Cairo 11566, Egypt

⁶

Department of Basic Science, Faculty of Engineering, Sinai University, Ismailia 41636, Egypt

Author to whom correspondence should be addressed.

Symmetry 2024, 16(7), 831; https://doi.org/10.3390/sym16070831

Submission received: 9 May 2024 / Revised: 15 June 2024 / Accepted: 21 June 2024 / Published: 2 July 2024

(This article belongs to the Section Mathematics)

Download

Browse Figures

Figure 1
Graphs of neutron flux <math display="inline"><semantics> <mrow> <mi>φ</mi> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> </semantics></math> in two dimensions for different values of <math display="inline"><semantics> <mrow> <mi>α</mi> </mrow> </semantics></math>. "> Figure 2
Graphs of neutron flux <math display="inline"><semantics> <mrow> <mi>φ</mi> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> </semantics></math> in three dimensions with <math display="inline"><semantics> <mrow> <msub> <mrow> <mi>φ</mi> </mrow> <mrow> <mn>0</mn> </mrow> </msub> <mfenced separators="|"> <mrow> <mi>x</mi> </mrow> </mfenced> <mo>=</mo> <msup> <mrow> <mi>x</mi> </mrow> <mrow> <mn>2</mn> </mrow> </msup> </mrow> </semantics></math> at α = 0.1, 0.2, 0.3, 0.4, 0.5, 0.90, 0.95, and 1. "> Figure 2 Cont.
Graphs of neutron flux <math display="inline"><semantics> <mrow> <mi>φ</mi> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> </semantics></math> in three dimensions with <math display="inline"><semantics> <mrow> <msub> <mrow> <mi>φ</mi> </mrow> <mrow> <mn>0</mn> </mrow> </msub> <mfenced separators="|"> <mrow> <mi>x</mi> </mrow> </mfenced> <mo>=</mo> <msup> <mrow> <mi>x</mi> </mrow> <mrow> <mn>2</mn> </mrow> </msup> </mrow> </semantics></math> at α = 0.1, 0.2, 0.3, 0.4, 0.5, 0.90, 0.95, and 1. ">

Versions Notes

Abstract

Building upon the previous research that solved neutron diffusion equations in simplified slab geometry, this study advances the field by addressing the more complex cylindrical geometry, focusing on neutron diffusion equations that are coupled with delayed neutrons in cylindrical reactors of fractional order. The method of solving used integrates the technique of residual power series (RPS) with the Laplace transform (LT) method. Anomalous neutron behavior is explained by examining the non-Gaussian scenario with various fractional parameters α. The LRPSM Laplace transform and residual power series method employed in this approach eliminates the complex difficulties. This simplicity makes the method particularly coherent with different fractional calculus applications. To validate the proposed method, numerical simulations are conducted with two different initial conditions representing distinct scenarios. The obtained results are presented in suitable tables and figures. It should be emphasized that this system is solved for the first time utilizing fractional calculus techniques. The outcomes are consistent with those achieved using the Adomian decomposition method.

Keywords:

nuclear reactor equation; Laplace transform method; residual power series method; factional differential equations; cylindrical reactors

1. Introduction

The utilization of nuclear reactor energy marks a pivotal aspect of a clean energy source. Presently, nuclear reactors contribute to 10% of the global energy demand. The inception of nuclear energy emerged from the concept of inducing fission in heavy nuclei through the bombardment of (thermal) slow neutrons.

While the neutron diffusion equation has been extensively studied in slab geometries, this research extends these concepts to the more intricate cylindrical reactor geometry [1,2,3]. Researchers in this field pursue scientific advancements by investigating various aspects. Different scenarios of nuclear reactor equations have been addressed using diverse methods. For instance, the homotopy perturbation method (HPM) has been employed to solve neutron diffusion equations for reactors exhibiting bare hemispherical and cylindrical symmetries for non-fractional systems [4,5]. Additionally, LRPSM has been used to tackle different branches of sciences [6,7,8].

The target of this work is solving the neutron diffusion equations coupled with delayed neutrons in cylindrical reactors fractionally using the Laplace residual power series method where the theoretical results are verified numerically.

Studying the nuclear reactor equation system with delayed neutrons provides a real application of fractional calculus in nuclear reactor theory. This application has historically been studied only for the slab reactor which can be considered as a simple example that is solved in Cartesian coordinates; in this work, the cylindrical reactor equations is solved in cylindrical coordinates which have an additional term in the equations of the system. This improvement of system has never been considered before.

This paper is prepared as follows. Section 2 introduces essential information about the studied phenomena. Section 3 presents the needed concepts and theories related to fractional calculus and the Laplace transforms. LRPSM details are provided in Section 4. Numerical examples of the theoretical study are considered in Section 5. The paper’s conclusion is given in Section 6.

2. The Physical Phenomena

Neutron distribution in reactors is analogous to the dissolution of a solute in physical chemistry. This procedure is expressed using a continuity equation.

\frac{1}{v_{c}} \frac{\partial φ (r, t)}{\partial t} = s (r, t) - Σ_{a} φ (r, t) - \nabla \cdot J (r, t),

(1)

In the context where

φ (r, t)

represents the neutron flux,

J (r, t)

denotes the neutron current density,

Σ_{a}

is the absorption cross-section in the macroscopic scale,

v_{c}

represents the velocity of the neutrons, and

s (r, t)

is the head of the neutron’s idiom, the relationship between

φ (r, t)

and

J (r, t)

can be expressed through Fick’s law.

J (r, t) = - D \nabla φ (r, t)

(2)

The neutrons are diffused according to the following equation which can be stated when D is designated as the neutron diffusion coefficient:

\frac{1}{v_{c}} \frac{\partial φ (r, t)}{\partial t} = s (r, t) - Σ_{a} φ (r, t) + D \nabla^{2} φ (r, t)

(3)

The source term

s (r, t)

can be written in nuclear reactor physics as

ν Σ_{f} φ (r, t)

where

ν

is the number of neutrons produced per fission and

Σ_{f}

is the fission cross-section in the macroscopic scale.

Now, it is important to introduce buckling

B^{2}

which is an essential nuclear reactor concept defined as

\frac{ν Σ_{f} φ (r, t) - Σ_{a}}{D}

This signifies the widely recognized equation for time-dependent neutron diffusion in nuclear reactors.

Taking a step further, the fractional analysis of the time-dependent form of this equation adds more detail by examining the physical impact of delayed neutrons. The investigation of a system that involves coupled neutron diffusion equations with delayed neutrons in fractional form is central to this study, which is important in nuclear reactor dynamics.

Delayed neutrons are crucial due to their role in influencing the operational lifespan of neutrons in reactor environments.

The effective lifetime of neutrons encompasses the immediate lifetime and an extended period marked by precursor elements. This additional delay far exceeds the prompt neutron lifetime, thereby prolonging the reactor’s time constant, which facilitates improved control.

Incorporating the effect of delayed neutrons into the

s (r, t)

term in Equation (3). The source term will be

s (r, t) = γ \sum_{f} φ (r, t) + λ ω (r, t)

(4)

A group of time-variant equations functions within a single-dimensional framework.

\frac{1}{v_{c}} \frac{\partial}{\partial t} φ (x, t) = D \frac{\partial^{2}}{\partial x^{2}} φ (x, t) + (- \sum_{a} + γ \sum_{f}) φ (x, t) + λ ω (x, t),

(5)

\frac{\partial}{\partial t} ω (x, t) = β ν \sum_{f} φ (x, t) - λ ω (x, t),

(6)

where the following

φ (x, 0) = φ_{0} (x), ω (x, 0) = \frac{β γ \sum_{f}}{λ} φ_{0} (x),

(7)

are needed conditions.

Here, the belated neutron intensity is symbolized by

φ (x, t)

, the average neutron yield per fission is

γ

, the decay constant is

λ

, and the proportion of delayed fission neutrons is

β

The term

ν Σ_{f} φ (r, t) - Σ_{a}

is interpreted as the overall cross-section

Σ

The dynamic equation highlights how neutron distribution changes over space and time within reactor settings.

Fractional calculus becomes notably important in the time-varying diffusion equations, especially when considering the delayed neutrons.

This aspect is termed anomalous diffusion, denoting a departure from the Gaussian distribution characteristic of regular diffusion.

Irregular cases of study are described by a non-linear, asymptotic association, especially one that follows a power-law pattern, between time and mean-squared displacement.

〈x^{2} (t)〉 ~ t^{α}

The α parameter is instrumental in categorizing the fractional nuclear reactor equation for the diffused neutrons, where super-diffusion (

α > 1

) and sub-diffusion (

α < 1

) are the divisions of irregular diffusion (

α \neq 1

When

α = 1

, the diffusion process adheres to a Gaussian model, represented mathematically by the integer-order neutron diffusion equation.

When neutrons are confined, their movement becomes slower, leading to what is termed sub-diffusion. Conversely, fast neutron movement occurs when they travel extensively in one direction without encountering collisions, a phenomenon referred to as super-diffusion [9,10].

Numerous studies have investigated the non-integer-order diffusion equation, particularly focusing on delayed neutrons. This includes research on the matrix representation of the fractional two-group in-reactor dynamics for slab reactors when the external source is applied [11], addressing neutron diffusion equations where the delayed neutrons are present in slab reactors [12], exploring a solution of kinetic nuclear reactor equation with a single delayed precursor concentration in for the one-dimensional slab reactors [13]. Within this paradigm, the fractional form of diffusion equations for neutrons with a single delayed neutron group may be articulated as follows:

\frac{1}{v_{c}} D_{t}^{α} φ (x, t) = D \nabla^{2} φ (x, t) + (γ \sum_{f} - \sum_{a}) φ (x, t) + λ ω (x, t),

(8)

D_{t}^{α} ω (x, t) = β γ \sum_{f} φ (x, t) - λ ω (x, t),

(9)

φ (x, 0) = φ_{0} (x), ω (x, 0) = \frac{β γ \sum_{f}}{λ} φ_{0} (x),

(10)

A wide array of computational and theoretical techniques including LTM, HPM, RPSM, and other methods are formulated to tackle fractional differential equations (FDEs) [6,7,14,15,16,17,18,19,20,21,22,23,24,25,26].

These techniques frequently face constraints, such as the need for extensive and intricate calculations.

To overcome these challenges, a connection to a transform operator is often necessary. In a previous work [24], the authors introduced a novel approach by combining RPSM with LTM to provide practical and series solutions for not only linear FDEs but also nonlinear FDEs.

Historically, all cases studying nuclear reactor equations have considered slab reactors in Cartesian geometry, representing the simplest case. This work, however, marks a significant advancement by exploring a more complex geometry, namely the cylindrical geometry, for cylindrical reactors.

Contrary to RPSM, LRPSM does not depend on fractional derivation for determining series coefficients but rather employs the concept of limits.

Consequently, LRPSM requires fewer calculations for the coefficient generation of RPSM. This technique is distinguished by its speed, low computer memory usage, and resistance to codes’ round-off errors. Furthermore, quicker convergence is achieved since LRPSM uses a set of equations involving several variables to obtain the power series coefficients.

This study utilizes LRPSM to effectively resolve the interconnected fractional system of prompt and delayed neutrons.

3. Preliminaries

The definition of fractional derivatives using integrals makes them non-local operators. As a result, a time-fractional derivative displays a memory effect by encapsulating associations about the function from previous times. These derivatives provide a critical improvement for a more accurate and thorough explanation of dynamic and complex system behavior as they account for non-local and historical dispersed impacts.

This section introduces the Laplace transform and fractional calculus within the framework of Caputo’s method. The development of the Laplace residual power series (LRPS) method, which is employed to solve fractional neutron diffusion equations involving a single delayed neutron group, is largely dependent on these components.

The basic concepts and ideas associated with these mathematical tools are presented.

Definition 1.

The Caputo operator of order

α

for

φ (χ, t)

is given by [6]

D_{t}^{α} φ (χ, t) = J_{t}^{n - α} \partial_{t}^{n} φ (χ, t), n - 1 < α \leq n, χ \in I . t > 0,

(11)

D_{t}^{α}

represents the Caputo derivative operator with a fractional order of

α

n \in N

, and

I

denotes an interval. Additionally,

J_{t}^{β}

is defined as the Riemann–Liouville integral operator with a fractional time order of β.

J_{t}^{β} φ (χ, t) = \{\begin{matrix} \frac{1}{Γ (β)} \int_{0}^{t} {(t - τ)}^{β - 1} φ (χ, τ) d τ, β > 0, t > τ \geq 0 \\ φ (χ, t), β = 0 \end{matrix}

(12)

Essential characteristics of the fractional derivative

D_{t}^{α}

will be introduced in the next lemma. Additional features can be found in Reference [26].

Lemma 1.

For

n - 1 < α \leq n

ε > - 1, t \geq 0

, and

μ \in R

, we have: [24]

\begin{matrix} (i) D_{t}^{α} t^{ε} = \frac{Γ (ε + 1)}{Γ (ε + 1 - α)} t^{ε - α} . \\ \begin{matrix} {(i i) D}_{t}^{α} μ = 0 \end{matrix} . \\ {(i i i) D}_{t}^{α} J_{t}^{α} φ (χ, t) = φ (χ, t) . \\ {(i v) J}_{t}^{α} D_{t}^{α} φ (χ, t) = φ (χ, t) - \sum_{j = 0}^{n - 1} \partial_{t}^{j} φ (χ, 0^{+}) \frac{t^{j}}{j!} . \end{matrix}

(13)

Definition 2.

The following definition applies the Laplace transform of

φ (χ, t)

[6]:

Φ (χ, s) = L [φ (χ, t)] ≔ \int_{0}^{\infty} e^{- s t} φ (χ, t) d t, s > δ,

(14)

Then, the inverse Laplace transform of

Ψ (χ, s)

will be applied

φ (χ, t) = L^{- 1} [Φ (χ, s)] ≔ \int_{c - i \infty}^{c + i \infty} e^{s t} Φ (χ, s) d s, c = R e (s) > c_{0} .

(15)

In this context,

c_{0}

is located within the right half-plane where the Laplace integral absolutely converges.

The next lemma delineates the fundamental properties of the LT and the Caputo sense of fractional derivative.

Lemma 2.

φ (χ, t)

is a piecewise continuous function on the interval

I \times [0, \infty)

and has exponential orders δ, and if

Φ (χ, s)

is the Laplace transform of

φ (χ, t)

, then [24]

\begin{matrix} (i) \lim_{s \to \infty} s Φ (χ, s) = φ (χ, 0), x \in I . \\ (i i) L [J_{t}^{α} φ (χ, t)] = s^{α - 1} Φ (χ, s), α > 0 . \\ (i i i) L [D_{t}^{α} φ (χ, t)] = s^{α} Φ (χ, s) - \sum_{k = 0}^{n - 1} s^{α - k - 1} \partial_{t}^{k} φ (χ, 0), n - 1 < α < n . \\ (i v) L [D_{t}^{m α} φ (χ, t)] = s^{m α} Φ (χ, s) - \sum_{k = 0}^{m - 1} s^{(m - k) α - 1} D_{t}^{k α} φ (χ, 0), 0 < α < 1, \end{matrix}

(16)

where

D_{t}^{m α} = D_{t}^{α} . D_{t}^{α} \dots D_{t}^{α}

(

m

-times).

The next result we provide is the basis for developing an LRPS solution for the partial differential equations (PDEs). It is a unique fractional expansion.

Theorem 1.

Assume

φ (χ, t)

is a piecewise continuous function on

I \times [0, \infty)

with exponential order δ. Consider that the function

Φ (χ, s),

which is the Laplace transform

L [φ (χ, t)]

, exhibits the following fractional expansion [24]:

Φ (χ, s) = \sum_{n = 0}^{\infty} \frac{y_{n} (χ)}{s^{n α + 1}}, 0 < α \leq 1, χ \in I, s > δ .

(17)

Then,

y_{n} (χ) = D_{t}^{n α} φ (χ, 0)

The next theorem describes the necessary and sufficient conditions for the series in expansion (17) to converge.

Theorem 2.

Assume

φ (χ, t)

is a piecewise continuous function on the interval

I \times [0, \infty),

and it has an exponential order of δ. Additionally,

Φ (χ, s)

is the Laplace transform of

φ (χ, t)

. The fractional expansion in Theorem 1 can serve as a representation of it [24].

It must be noted that the creators of this method used it and other methods in history to solve a nuclear reactor system of equations [6,7].

4. The Nuclear Reactor Equations Analytical Solution

Utilizing the Laplace residual power series method (LRPSM) to derive a series solution for the fractional neutron diffusion equations is the target in this stage.

\begin{matrix} \frac{1}{v_{c}} D_{t}^{α} φ (x, t) = D \nabla^{2} φ (x, t) + (γ \sum_{f} - \sum_{a}) φ (x, t) + \frac{1}{x} \frac{\partial φ (x, t)}{\partial x} + λ ω (x, t), \\ D_{t}^{α} ω (x, t) = β γ \sum_{f} φ (x, t) - λ ω (x, t), \end{matrix}

(18)

given the specified initial conditions

φ (x, 0) = φ_{0} (x), ω (x, 0) = \frac{β ν \sum_{f}}{λ} φ_{0} (x),

(19)

Rewriting the coupled equations in (18) as

B = β γ \sum_{f}

and

\sum = γ \sum_{f} - \sum_{a}

, we obtain the following:

\begin{matrix} D_{t}^{α} φ (x, t) = v_{c} D φ_{x x} (x, t) + v_{c} \sum φ (x, t) + \frac{v_{c}}{x} φ_{x} (x, t) + v_{c} λ ω (x, t), \\ D_{t}^{α} ω (x, t) = B φ (x, t) - λ ω (x, t), \end{matrix}

(20)

The following nuclear physics initial conditions will be used

φ (x, 0) = φ_{0} (x), ω (x, 0) = \frac{B}{λ} φ_{0} (x),

(21)

When the coupled Equation (20) are subjected to the Laplace transform, we have

\begin{matrix} L [D_{t}^{α} φ (x, t)] = v_{c} D L [φ_{x x} (x, t)] + v_{c} \sum L [φ (x, t)] + \frac{v_{c}}{x} L [φ_{x} (x, t)] + v_{c} λ L [ω (x, t)], \\ L [D_{t}^{α} ω (x, t)] = B [φ (x, t)] - λ L [ω (x, t)] . \end{matrix}

(22)

Using the initial conditions in Equation (21) and Lemma 2, we may formulate Equation (22) in the following form

\begin{matrix} s^{α} Φ (x, s) - s^{α - 1} φ_{0} (x) = v_{c} D Φ_{xx} (x, s) + v_{c} \sum Φ (x, s) + \frac{v_{c}}{x} Φ_{x} (x, s) + v_{c} λ Ω (x, s), \\ s^{α} Ω (x, s) - s^{α - 1} \frac{B}{λ} φ_{0} (x) = B Φ (x, s) - λ Ω (x, s), \end{matrix}

(23)

where

Φ (x, s) = L [φ (x, t)],

Ω (x, s) = L [ω (x, t)] .

We can rewrite Equation (23) as

\begin{matrix} Φ (x, s) - \frac{1}{s} φ_{0} (x) - \frac{v_{c} D}{s^{α}} Φ_{xx} (x, s) - \frac{v_{c} \sum}{s^{α}} Φ (x, s) - \frac{v_{c}}{x s^{α}} Φ_{x} (x, s) - \frac{v_{c} λ}{s^{α}} Ω (x, s) = 0, \\ Ω (x, s) - \frac{B}{λ s} φ_{0} (x) - \frac{B}{s^{α}} Φ (x, s) + \frac{λ}{s^{α}} Ω (x, s) = 0 . \end{matrix}

(24)

The coupled Equation (24) forms a system of linear partial differential equations which involve derivatives with respect to the variable

x

. According to the LRPS approach, the series solution for the system (24) will be as follows:

\begin{matrix} Φ (x, s) = \sum_{n = 0}^{\infty} \frac{y_{n} (x)}{s^{n α + 1}}, x \in I, s > δ \geq 0, \\ Ω (x, s) = \sum_{n = 0}^{\infty} \frac{z_{n} (x)}{s^{n α + 1}}, x \in I, s > δ \geq 0 . \end{matrix}

(25)

The kth truncated series of

Φ (x, s), Ω (x, s)

can be represented as follows using Lemma 2:

\begin{matrix} Φ_{k} (x, s) = \frac{φ_{0} (x)}{s} + \sum_{n = 1}^{k} \frac{y_{n} (x)}{s^{n α + 1}}, x \in I, s > δ \geq 0, \\ Ω_{k} (x, s) = \frac{B}{λ s} φ_{0} (x) + \sum_{n = 1}^{k} \frac{z_{n} (x)}{s^{n α + 1}}, x \in I, s > δ \geq 0 . \end{matrix}

(26)

The next step is to obtain the anonymous coefficients of Equation (26) by establishing the Laplace residual functions of Equation (24)

\begin{matrix} L R e s (Φ (x, s)) = Φ (x, s) - \frac{1}{s} φ_{0} (x) - \frac{v_{c} D}{s^{α}} Φ_{xx} (x, s) - \frac{v_{c} \sum}{s^{α}} Φ (x, s) - \frac{v_{c}}{x s^{α}} Φ_{x} (x, s) - \frac{v_{c} λ}{s^{α}} Ω (x, s) \\ L R e s (C (x, s)) = C (x, s) - \frac{B}{λ s} φ_{0} (x) - \frac{B}{s^{α}} Φ (x, s) + \frac{λ}{s^{α}} Ω (x, s) . \end{matrix}

(27)

Moreover, the kth Laplace residual serves as

\begin{matrix} {L R e s}_{k} (Φ (x, s)) = Φ_{k} (x, s) - \frac{1}{s} φ_{0} (x) - \frac{v_{c} D}{s^{α}} {Φ_{k}}_{xx} (x, s) - \frac{v_{c} \sum}{s^{α}} Φ_{k} (x, s) - \frac{v_{c}}{x s^{α}} {Φ_{k}}_{x} (x, s) - \frac{v_{c} λ}{s^{α}} Ω_{k} (x, s) \\ {L R e s}_{k} (Ω (x, s)) = Ω_{k} (x, s) - \frac{B}{λ s} φ_{0} (x) - \frac{B}{s^{α}} Φ_{k} (x, s) + \frac{λ}{s^{α}} Ω_{k} (x, s) . \end{matrix}

(28)

Since

L R e s (Φ (x, s)) = 0

and

L R e s (Ω (x, s)) = 0

, we find that

s^{k α + 1} L R e s (Φ (x, s)) = 0, s^{k α + 1} L R e s (Ω (x, s)) = 0

. Therefore,

\lim_{s \to \infty} (s^{k α + 1} {L R e s}_{k} (Φ (x, s))) = 0, \lim_{s \to \infty} (s^{k α + 1} {L R e s}_{k} (Ω (x, s))) = 0 f o r k = 0,1, 2, \dots .

(29)

Then, to obtain

y_{1} (x)

and

z_{1} (x)

in Equation (26), we substitute in the first Laplace residual functions as

Φ_{1} (x, s) = \frac{φ_{0} (x)}{s} + \frac{y_{1} (x)}{s^{α + 1}}

and

Ω_{1} (x, s) = \frac{B}{λ s} φ_{0} (x) + \frac{z_{1} (x)}{s^{α + 1}}

, and we obtain

\begin{matrix} {L R e s}_{1} (Φ (x, s)) = \frac{y_{1} (x)}{s^{α + 1}} - \frac{v_{c} D}{s^{α}} [\frac{φ_{0}^{″} (x)}{s} + \frac{y_{1}^{″} (x)}{s^{α + 1}}] - \frac{v_{c} \sum}{s^{α}} [\frac{φ_{0} (x)}{s} + \frac{y_{1} (x)}{s^{α + 1}}] - \frac{v_{c}}{x s^{α}} [\frac{φ_{0}^{'} (x)}{s} + \frac{y_{1}^{'} (x)}{s^{α + 1}}] - \frac{v_{c} λ}{s^{α}} [\frac{B}{λ s} φ_{0} (x) + \frac{z_{1} (x)}{s^{α + 1}},] \\ {L R e s}_{1} (Ω (x, s)) = \frac{z_{1} (x)}{s^{α + 1}} - \frac{B}{s^{α}} [\frac{φ_{0} (x)}{s} + \frac{y_{1} (x)}{s^{α + 1}}] + \frac{λ}{s^{α}} [\frac{B}{λ s} φ_{0} (x) + \frac{z_{1} (x)}{s^{α + 1}}] . \end{matrix}

(30)

Now, by solving

\lim_{s \to \infty} s^{α + 1} {L R e s}_{1} (Φ (x, s)) = 0

\lim_{s \to \infty} s^{α + 1} {L R e s}_{1} (C (x, s)) = 0,

we find that

y_{1} (x) = v_{c} D φ_{0}^{″} (x) + v_{c} \sum φ_{0} (x) + \frac{v_{c}}{x} φ_{0}^{'} (x) + v_{c} B ω_{0} (x),

(31)

z_{1} (x) = 0 .

(32)

To obtain

y_{2} (x)

and

z_{2} (x)

in Equation (26), we substitute in the second Laplace residual functions

Φ_{2} (x, s) = \frac{φ_{0} (x)}{s} + \frac{y_{1} (x)}{s^{α + 1}} + \frac{y_{2} (x)}{s^{2 α + 1}}

and

Ω_{2} (x, s) = \frac{B}{λ s} φ_{0} (x) + \frac{z_{2} (x)}{s^{2 α + 1}}

, and we obtain

\begin{matrix} {L R e s}_{2} (Φ (x, s)) = \frac{y_{1} (x)}{s^{α + 1}} + \frac{y_{2} (x)}{s^{2 α + 1}} - \frac{v_{c} D}{s^{α}} [\frac{φ_{0}^{″} (x)}{s} + \frac{y_{1}^{″} (x)}{s^{α + 1}} + \frac{y_{2}^{″} (x)}{s^{2 α + 1}}] - \frac{v_{c} \sum}{s^{α}} [\frac{φ_{0} (x)}{s} + \frac{y_{1} (x)}{s^{α + 1}} + \frac{y_{2} (x)}{s^{2 α + 1}}] - \frac{v_{c}}{x s^{α}} [\frac{φ_{0}^{'} (x)}{s} + \frac{y_{1}^{'} (x)}{s^{α + 1}} + \frac{y_{2}^{'} (x)}{s^{2 α + 1}}] - \frac{v_{c} λ}{s^{α}} [\frac{B}{λ s} φ_{0} (x) + \frac{z_{2} (x)}{s^{2 α + 1}}], \\ {L R e s}_{2} (Ω (x, s)) = \frac{z_{2} (x)}{s^{2 α + 1}} - \frac{B}{s^{α}} [\frac{φ_{0} (x)}{s} + \frac{y_{1} (x)}{s^{α + 1}} + \frac{y_{2} (x)}{s^{2 α + 1}}] + \frac{λ}{s^{α}} [\frac{B}{λ s} φ_{0} (x) + \frac{z_{2} (x)}{s^{2 α + 1}}] . \end{matrix}

(33)

Now, by solving

\lim_{s \to \infty} s^{2 α + 1} {L R e s}_{2} (Φ (x, s)) = 0

\lim_{s \to \infty} s^{2 α + 1} {L R e s}_{2} (Ω (x, s)) = 0,

we find that

y_{2} (x) = v_{c} D y_{1}^{″} (x) + v_{c} \sum y_{1} (x) + \frac{v_{c}}{x} y_{1}^{'} (x),

(34)

z_{2} (x) = B y_{1} (x) .

(35)

Then, we have the following:

\begin{matrix} {L R e s}_{3} (Φ (x, s)) = \frac{y_{1} (x)}{s^{α + 1}} + \frac{y_{2} (x)}{s^{2 α + 1}} + \frac{y_{3} (x)}{s^{3 α + 1}} - \frac{v_{c} D}{s^{α}} [\frac{φ_{0}^{″} (x)}{s} + \frac{y_{1}^{″} (x)}{s^{α + 1}} + \frac{y_{2}^{″} (x)}{s^{2 α + 1}} + \frac{y_{3}^{″} (x)}{s^{3 α + 1}}] - \frac{v_{c} \sum}{s^{α}} [\frac{φ_{0}^{'} (x)}{s} + \frac{y_{1}^{'} (x)}{s^{α + 1}} + \frac{y_{2}^{'} (x)}{s^{2 α + 1}} + \frac{y_{3}^{'} (x)}{s^{3 α + 1}}] - \frac{v_{c} λ}{s^{α}} [\frac{B}{λ s} φ_{0} (x) + \frac{z_{2} (x)}{s^{2 α + 1}} + \frac{z_{3} (x)}{s^{3 α + 1}}], \\ {L R e s}_{3} (Ω (x, s)) = \frac{z_{2} (x)}{s^{2 α + 1}} + \frac{z_{3} (x)}{s^{3 α + 1}} - \frac{B}{s^{α}} [\frac{φ_{0} (x)}{s} + \frac{y_{1} (x)}{s^{α + 1}} + \frac{y_{2} (x)}{s^{2 α + 1}} + \frac{y_{3} (x)}{s^{3 α + 1}}] + \frac{λ}{s^{α}} [\frac{B}{λ s} φ_{0} (x) + \frac{z_{2} (x)}{s^{2 α + 1}} + \frac{z_{3} (x)}{s^{3 α + 1}}] . \end{matrix}

(36)

Now, by solving

\lim_{s \to \infty} s^{3 α + 1} {L R e s}_{3} (Φ (x, s)) = 0

\lim_{s \to \infty} s^{3 α + 1} {L R e s}_{3} (Ω (x, s)) = 0,

we find that

y_{3} (x) = v_{c} D y_{2}^{″} (x) + v_{c} \sum y_{2} (x) + \frac{v_{c}}{x} y_{2}^{'} (x) + v_{c} λ z_{2} (x),

(37)

z_{3} (x) = B y_{2} (x) - λ z_{2} (x) .

(38)

Proceeding similarly, we replace the kth truncated series

Φ_{k} (x, s)

Ω_{k} (x, s)

with the kth Laplace residual function

{L R e s}_{k} (Φ (x, s)) {, L R e s}_{k} (Ω (x, s))

. Subsequently, we multiply the resulting equations by

s^{k α + 1}

and consider the limit as

s \to \infty

. The recurrence relation shown below can be used to obtain

y_{k + 1} (x)

z_{k + 1} (x)

for

k \geq 2

y_{n} (x) = v_{c} D y_{n - 1}^{″} (x) + v_{c} \sum y_{n - 1} (x) + \frac{v_{c}}{x} y_{n - 1}^{'} (x) + v_{c} λ z_{n - 1} (x),

(39)

z_{n} (x) = B y_{n - 1} (x) - λ z_{n - 1} (x), n = 1,2, 3, \dots

(40)

where

z_{0} (x) = \frac{B}{λ} φ_{0} (x)

The system of Equation (24) in accordance with what was introduced is

\begin{matrix} Φ (x, s) = \frac{φ_{0} (x)}{s} + \sum_{n = 1}^{\infty} \frac{y_{n} (x)}{s^{n α + 1}}, x \in I, s > δ \geq 0 . \\ Ω (x, s) = \frac{B}{λ s} φ_{0} (x) + \sum_{n = 1}^{\infty} \frac{z_{n} (x)}{s^{n α + 1}}, x \in I, s > δ \geq 0 . \end{matrix}

(41)

Here, the inverse Laplace transform is applied in Equation (41) enables us to derive the series solution for the system outlined in Equation (18) as follows:

\begin{matrix} φ (x, t) = φ_{0} (x) + \sum_{n = 1}^{\infty} \frac{y_{n} (x) t^{n α}}{Γ (n α + 1)}, t \geq 0, x \in I . \\ ω (x, t) = \frac{B}{λ} φ_{0} (x) + \sum_{n = 2}^{\infty} \frac{z_{n} (x) t^{n α}}{Γ (n α + 1)}, t \geq 0, x \in I, \end{matrix}

(42)

where

\begin{matrix} y_{n} (x) = v_{c} D y_{n - 1}^{″} (x) + v_{c} \sum y_{n - 1} (x) + \frac{v_{c}}{x} y_{n - 1}^{'} (x) + v_{c} λ z_{n - 1} (x), \\ z_{n} (x) = B y_{n - 1} (x) - λ z_{n - 1} (x), n = 1,2, 3, \dots \end{matrix}

(43)

where

z_{0} (x) = \frac{B}{λ} φ_{0} (x)

5. Numerical Results and Discussion

The proposed theoretical framework is verified by tackling the fractional time-dependent nuclear reactor equations and neutron diffusion equation and incorporating delayed neutrons through LRPSM. The results obtained from this approach will be compared with numerical results derived from the neutron diffusion cross-section data as referenced, as in [10]. The necessary nuclear reactor constants are provided in Table 1.

It is important to emphasize that these data are specified in Section 1.

The functions

y_{n} (x)

and

z_{n} (x)

were calculated using codes implemented in MATHEMATICA software 12. Once these values are determined, the data necessary for constructing the tables and figures are obtained.

In this study, we will examine two numerical cases depending on the initial condition.

The first example considers the initial condition

φ (x, 0) = 1

. The distribution of the neutron in the cylindrical reactor, represented by the neutron flux, is tabulated in Table 1 and illustrated in Figure 1.

The flux was computed at various values of time, considering both a non-Gaussian distribution, (

α < 1

, sub-diffusion), and Gaussian distribution, (

α = 1

, integer order case). In addition to that, Figure 1 and Table 2 show that the neutron flux increases with time.

The second case study considers the initial condition

φ (x, 0) = x^{2}

where the neutron flux for different values of time under this condition is illustrated in Figure 2.

Presently, a varied initial value problem

φ (x, 0) = x^{2}

is employed in examining the following scenario, and the outcomes are elucidated in Figure 2.

It is evident that the behavior of the flux becomes closer to the non- fractional case as

α

approaches to unity. This indicates that the anomalous behavior decreases which means the effect of the delayed neutrons diminishes and eventually vanishes.

The nonfractional case (α = 0.1) for the cylindrical reactor is presented in the literature [5].

It is important to note that, in both examined numerical cases [9,10,12], the time-dependent flux exhibits an increase over time.

The significance of this work lies in the investigation of the behavior of delayed neutron flux, which affects irregular diffusion. The solution of the time-dependent fractional neutron diffusion equations, as shown in Equation (42), provides a mathematical expression for this flux behavior.

6. Conclusions

The fractional neutron diffusion equations, including belated neutrons, were effectively studied utilizing the Laplace residual power series approach. We were able to produce an analytical fractional series solution using this method.

The soundness of the method’s theoretical framework is verified through two distinct numerical case studies in nuclear physics. Flux calculations for both cases are derived, shedding light on the phenomenon of anomalous diffusion,

α < 1

(sub-diffusion). The investigation involves the consideration of neutron fluxes that exhibit an increase over time. Using LRPSM in fractional calculus, the novel notion of irregular diffusion, which characterizes diffusion patterns that are not Gaussian, has been thoroughly described. The application of LRPSM in cylindrical nuclear reactors can improve the study more complicated nuclear reactor geometries and the investigation of different scientific problems.

Author Contributions

Conceptualization, M.S.; methodology, A.H.; software, A.H.; validation S.I., D.F. and F.A.A.; formal analysis, A.H.; investigation, F.A.A.; resources, M.S.; data curation, F.A.A.; writing—original draft preparation, D.F. and F.A.A.; writing—review and editing, S.I.; visualization, S.I.; supervision, F.A.A.; project administration, A.H.; funding acquisition, F.A.A. All authors have read and agreed to the published version of the manuscript.

Funding

Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2024R174), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia.

Data Availability Statement

Data are contained within the article.

Acknowledgments

Thanks for Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2024R174), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia.

Conflicts of Interest

The authors declare no conflicts of interest.

References

Stacey, W.M. Nuclear Reactor Physics; John Wiley & Sons: Boston, MA, USA, 2001. [Google Scholar]
Duderstadt, J.J.; Hamilton, L.J. Nuclear Reactor Analysis; John Wiley & Sons: Boston, MA, USA, 1976. [Google Scholar]
Lamarsh, J.R. Introduction to Nuclear Engineering, 2nd ed.; Addison-Wesley: Boston, MA, USA, 1983. [Google Scholar]
Khasawneh, K.; Dababneh, S.; Odibat, Z. A solution of the neutron diffusion equation in hemispherical symmetry using the homotopy perturbation method. Ann. Nucl. Energy 2009, 36, 1711–1717. [Google Scholar] [CrossRef]
Dababneh, S.; Khasawneh, K.; Odibat, Z. An alternative solution of the neutron diffusion equation in cylindrical symmetry. Ann. Nucl. Energy 2010, 38, 1140–1143. [Google Scholar] [CrossRef]
El-Ajou, A.; Shqair, M.; Ghabar, I.; Burqan, A.; Saadeh, R. A solution for the neutron diffusion equation in the spherical and hemispherical reactors using the residual power series. Front. Phys. 2023, 11, 1229142. [Google Scholar] [CrossRef]
Alquran, M.; Ali, M.; Alsukhour, M.; Jaradat, I. Promoted residual power series technique with Laplace transform to solve some time-fractional problems arising in physics. Results Phys. 2020, 19, 103667. [Google Scholar] [CrossRef]
Alaroud, M. Application of Laplace residual power series method for approximate solutions of fractional IVP’s. Alex. Eng. J. 2022, 61, 1585–1595. [Google Scholar] [CrossRef]
Vyawahare, V.; Nataraj, P. Fractional-Order Modeling of Nuclear Reactor: From Subdiffusive Neutron Transport to Control-Oriented Models; Springer Nature: Singapore, 2018. [Google Scholar]
Ray, S. Fractional Calculus with Applications for Nuclear Reactor Dynamics; CRC Press: Broken, NW, USA, 2016. [Google Scholar]
Aboanber, A.; Nahla, A.; Aljawazneh, S. Fractional two energy groups matrix representation for nuclear reactor dynamics with an external source. Ann. Nucl. Energy 2021, 153, 108062. [Google Scholar] [CrossRef]
Sardar, T.; Ray, S.; Bera, R. The solution of coupled fractional neutron diffusion equations with delayed neutrons. Int. J. Nuclear Energy Sci. Technol. 2010, 5, 105–133. [Google Scholar] [CrossRef]
Khaled, S. Exact solution of the one-dimensional neutron diffusion kinetic equation with one delayed precursor concentration in Cartesian geometry. AIMS Math. 2022, 7, 12364–12373. [Google Scholar] [CrossRef]
Sarhan, A.; Burqan, A.; Saadeh, R.; Al-Zhour, Z. Analytical Solutions of the Nonlinear Time-Fractional Coupled Boussinesq-Burger Equations Using Laplace Residual Power Series Technique. Fractal Fract. 2022, 6, 631. [Google Scholar] [CrossRef]
Kazem, S. Exact solution of some linear fractional differential equations by Laplace transform. Int. J. Nonlinear Sci. 2013, 16, 3–11. [Google Scholar]
Yavuz, M.; Özdemir, N. Numerical inverse Laplace homotopy technique for fractional heat equations. Therm. Sci. 2018, 22, 185–194. [Google Scholar] [CrossRef]
Wang, Q. Numerical solutions for fractional KdV–Burgers equation by Adomian decomposition method. Appl. Math. Comput. 2006, 182, 1048–1055. [Google Scholar] [CrossRef]
Cao, Y.; Nikan, O.; Avazzadeh, Z. A localized meshless technique for solving 2D nonlinear integro-differential equation with multi-term kernels. Appl. Numer. Math. 2023, 183, 140–156. [Google Scholar] [CrossRef]
Akram, T.; Abbas, M.; Ali, A. A numerical study on time fractional Fisher equation using an extended cubic B-spline approximation. J. Math. Comput. Sci. 2021, 22, 85–96. [Google Scholar] [CrossRef]
Oderinu, R.A.; Owolabi, J.A.; Taiwo, M. Approximate solutions of linear time-fractional differential equations. J. Math. Comput. Sci. 2023, 29, 60–72. [Google Scholar] [CrossRef]
Ali, A.; Abbas, M.; Akram, T. New group iterative schemes for solving the two-dimensional anomalous fractional sub-diffusion equation. J. Math. Comput. Sci. 2021, 22, 119–127. [Google Scholar] [CrossRef]
Liao, H.; Zhang, D.; Zhao, Y. Sharp error estimate of a nonuniform L1 formula for time-fractional reaction subdiffusion equations. SIAM J. Numer. Anal. 2018, 56, 1122–1133. [Google Scholar] [CrossRef]
Lyu, P.; Vong, S. A symmetric fractional-order reduction method for direct nonuniform approximations of semilinear diffusion-wave equations. J. Sci. Comput. 2022, 93, 34. [Google Scholar] [CrossRef]
Eriqat, T.; El-Ajou, A.; Moa’ath, N.O.; Al-Zhour, Z.; Momani, S. A new attractive analytic approach for solutions of linear and nonlinear neutral fractional pantograph equations. Chaos Solitons Fractals 2020, 138, 109957. [Google Scholar] [CrossRef]
Batiha, I.; Shqair, M.; Abu-Sei’leek, M.H.; Alshorm, S.; Abd El-Azeem, S.; Momani, S. Two Energy Groups Neutron Diffusion Model in Spherical Reactors. Results Nonlinear Anal. 2024, 7, 160–173. [Google Scholar]
Burqan, A.; Shqair, M.; El-Ajou, A.; Ismaeel, S.M.E.; AlZhour, Z. Analytical solutions to the coupled fractional neutron diffusion equations with delayed neutrons system using Laplace transform method. AIMS Math. 2023, 8, 19297–19312. [Google Scholar] [CrossRef]

Figure 1. Graphs of neutron flux

φ (x, t)

in two dimensions for different values of

α

Figure 1. Graphs of neutron flux

φ (x, t)

in two dimensions for different values of

α

Figure 2. Graphs of neutron flux

φ (x, t)

in three dimensions with

φ_{0} (x) = x^{2}

at α = 0.1, 0.2, 0.3, 0.4, 0.5, 0.90, 0.95, and 1.

Figure 2. Graphs of neutron flux

φ (x, t)

in three dimensions with

φ_{0} (x) = x^{2}

at α = 0.1, 0.2, 0.3, 0.4, 0.5, 0.90, 0.95, and 1.

Table 1. The neutron diffusion cross-section data.

$v_{c}$	$B$	$D$	$λ$	$Σ$
$220,000 c m / s$	$0.000735 {c m}^{- 1}$	$0.356$	$0.08 s^{- 1}$	$0.005 {c m}^{- 2}$ .

Table 2. The flux simulations for various values of

α

Table 2. The flux simulations for various values of

α

t	$α = 0.1$	$α = 0.2$	$α = 0.3$	$α = 0.4$	$α = 0.5$	$α = 0.6$	$α = 0.7$	$α = 0.8$	$α = 0.9$	$α = 1$
0.00010	1.0001 × 10¹⁴	2.3998 × 10¹²	5.1358 × 10¹⁰	1.0172 × 10⁹	1.9570 × 10⁷	409,582	12,938.5	1069.01	218.945	67.1205
0.00039	1.7233 × 10¹⁴	7.1214 × 10¹²	2.6187 × 10¹¹	8.8447 × 10⁹	2.8255 × 10⁸	8.93672 × 10⁶	311,117	15,528	1642	366.761
0.00068	2.1524 × 10¹⁴	1.1106 × 10¹³	5.0966 × 10¹¹	2.1440 × 10¹⁰	8.4808 × 10⁸	3.26117 × 10⁷	1.30286 × 10⁶	63,492.3	5186.97	887.538
0.00097	2.4809 × 10¹⁴	1.4753 × 10¹³	7.7999 × 10¹¹	3.7766 × 10¹⁰	1.7145 × 10⁹	7.50461 × 10⁷	3.32891 × 10⁶	167,752	12,286.8	1733.47
0.00126	2.7544 × 10¹⁴	1.8184 × 10¹³	1.0671 × 10¹¹	5.7313 × 10¹⁰	2.8813 × 10⁹	1.38975 × 10⁸	6.69906 × 10⁶	352,590	24,648.1	3033.51
0.00155	2.9923 × 10¹⁴	2.1459 × 10¹³	1.3677 × 10¹²	7.9754 × 10¹⁰	4.3478 × 10⁹	2.26674 × 10⁸	1.17040 × 10⁷	643,464	44,222.3	4941.57
0.00184	3.2047 × 10¹⁴	2.4613 × 10¹³	1.6798 × 10¹²	1.0485 × 10¹¹	6.1139 × 10⁹	3.40119 × 10⁸	1.86209 × 10⁷	1.0668 × 10⁶	73,184.9	7636.5
0.00213	3.3979 × 10¹⁴	2.7669 × 10¹³	2.0019 × 10¹²	1.3244 × 10¹¹	8.1791 × 10⁹	4.81079 × 10⁸	2.77163 × 10⁷	1.64988 × 10⁶	113,920	11,322.1
0.00242	3.5758 × 10¹⁴	3.0642 × 10¹³	2.3329 × 10¹²	1.6236 × 10¹¹	1.0543 × 10¹⁰	6.51161 × 10⁸	3.92481 × 10⁷	2.42072 × 10⁶	169,007	16,227.1
0.00271	3.7414 × 10¹⁴	3.3544 × 10¹³	2.6718 × 10¹²	1.9451 × 10¹¹	1.3206 × 10¹⁰	8.51847 × 10⁸	5.34675 × 10⁷	3.40801 × 10⁶	241,213	22,605.1
0.00300	3.8966 × 10¹⁴	3.6384 × 10¹³	3.0181 × 10¹²	2.2879 × 10¹¹	1.6168 × 10¹⁰	1.08452 × 10⁹	7.06193 × 10⁷	4.64106 × 10⁶	333,482	30,734.9

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Filali, D.; Shqair, M.; Alghamdi, F.A.; Ismaeel, S.; Hagag, A. Solving a Novel System of Time-Dependent Nuclear Reactor Equations of Fractional Order. Symmetry 2024, 16, 831. https://doi.org/10.3390/sym16070831

AMA Style

Filali D, Shqair M, Alghamdi FA, Ismaeel S, Hagag A. Solving a Novel System of Time-Dependent Nuclear Reactor Equations of Fractional Order. Symmetry. 2024; 16(7):831. https://doi.org/10.3390/sym16070831

Chicago/Turabian Style

Filali, Doaa, Mohammed Shqair, Fatemah A. Alghamdi, Sherif Ismaeel, and Ahmed Hagag. 2024. "Solving a Novel System of Time-Dependent Nuclear Reactor Equations of Fractional Order" Symmetry 16, no. 7: 831. https://doi.org/10.3390/sym16070831

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Solving a Novel System of Time-Dependent Nuclear Reactor Equations of Fractional Order

Abstract

1. Introduction

2. The Physical Phenomena

3. Preliminaries

4. The Nuclear Reactor Equations Analytical Solution

5. Numerical Results and Discussion

6. Conclusions

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI