A solution for the neutron diffusion equation in the spherical and hemispherical reactors using the residual power series

A novel analytical solution to the neutron diffusion equation is proposed in this study using the residual power series approach for both spherical and hemispherical fissile material reactors. Various boundary conditions are investigated, including zero flux on the boundary, zero flux on the extrapolated boundary distance, and the radiation boundary condition (RBC). The study also shows how two hemispheres with opposing flat faces interact. We give numerical results for the same energy neutrons diffused in pure ${}^{239}Pu$. By qualitative comparison with the homotopy perturbation method and Bessel function-based solutions, the residual power series method (RPSM) presents accurate series solutions that converge to the exact solutions, as shown in this study. Moreover, numerical results were shown to be improved by the computer implementation of the analytic solutions.

1 Introduction

Modern physics is taken into consideration after classical physics failed to explain various phenomena, and many branches of modern physics developed, including nuclear physics. To perform nuclear fission in a nuclear reactor, it is essential to have a neutron in the thermal state to break through the nucleus, as each fission produces more than one neutron, leading to a new reaction where a chain reaction is created [1].

The continuity equation, known as the time-dependent transfer equation, expresses the distribution of neutrons in the reactor, which is given as

$\frac{1}{v}\frac{\partial \varphi \left(\mathit{r},t\right)}{\partial t}=s\left(\mathit{r},t\right)-{\mathrm{\Sigma }}_a\left(\mathit{r}\right)\varphi \left(\mathit{r},t\right)-\nabla .J\left(\mathit{r},t\right),$

where $\upsilon$ is the neutron velocity, $\varphi \left(\mathit{r},t\right)$ is the neutron flux, $s\left(\mathit{r},t\right)$ is the neutron source term, ${\mathrm{\Sigma }}_a$ is an absorption macroscopic cross section, and $\mathit{r}$ is the position vector at the time $t$.

The continuity equation involves both neutron current $J\left(\mathit{r},t\right)$ and neutron flux $\varphi \left(\mathit{r},t\right)$, where the relationship between current and the flux can be shown using Fick’s law, which is given in the following mathematical equation:

$J\left(\mathit{r},t\right)=-D\nabla \varphi \left(\mathit{r},t\right),$

where $\mathrm{D}$ is the diffusion coefficient.

By applying Fick’s law, the time-dependent transport equation can be simplified to the following time-dependent neutron diffusion equation:

$\frac{1}{v}\frac{\partial \varphi \left(\mathit{r},t\right)}{\partial t}=s\left(\mathit{r},t\right)-{\mathrm{\Sigma }}_a\left(\mathit{r}\right)\varphi \left(\mathit{r},t\right)+D{\nabla }^2\varphi \left(\mathit{r},t\right),$

whose solution is considered very important in studying the behavior of neutrons in reactors [2, 3].

The stability of the nuclear reactor has been studied in detail in the past [4]. The condition of criticality in a nuclear reactor is essential to ensure its stability, as the number of neutrons in the reactor must be constant at any time. The flow in the critical state must be time-independent, which can be mathematically expressed using a time-independent diffusion equation as follows:

$D{\nabla }^2\varphi \left(\mathit{r},\theta \right)+\left({\upsilon \sum }_f-{\sum }_a\right)\varphi \left(\mathit{r},\theta \right)=0.(1)$

Various techniques, analytical and numerical, have been employed to solve the neutron diffusion equation throughout the past few decades, such as the homotopy perturbation method, the Laplace transform process, and the Galerkin method with B-splines [5–21]. To the best of our knowledge, the residual power series method (RPSM) [21–29], an analytical technique, has never been used to handle issues relating to the neutron diffusion equation.

RPSM is a relatively new technique that provides exact and approximate analytical solutions in the series form to many differential and integrodifferential equations of integer and fractional orders [21–29]. The idea of the RPSM is based on using the residual function and the derivative operator concepts in determining the coefficients of the solution, which is assumed to be in a power series form. In RPSM, it is usual to use a classical derivative operator to determine the coefficients of a power series solution if the differential equation is of integer order. However, a fractional derivative operator is used to determine the coefficients if the differential equation is of fractional order.

Later, we will rewrite Eq. 1 using spherical coordinates in spherical and hemispherical reactors, which gives rise to single equations. This work aims to adapt and modify RPSM to solve Eq. 1 after being formulated in spherical coordinates. This modification is a regeneration of RPSM and is presented for the first time in this work. This modification relies on using a fractional derivative operator, such as the Caputo fractional derivative (CFD) [30, 31] or the conformable fractional derivative [21, 32], in differential equations of the integer order. Since the equations to be solved are singular, we impose the solution in a Frobenius series form. Therefore, the new modified method can be considered an alternative method to the Frobenius method in determining the power series solution coefficients.

In fact, several types of fractional derivative operators can be used in our hypothetical modification of the RPSM, such as the CFD operator, the Riemann–Liouville fractional derivative operator [33], and the conformable fractional derivative operator. The CFD is employed in the modified method without bias or preference to determine the Frobenius series coefficients. The CFD operator is defined using the following formula:

$D_t^{\delta }f\left(t\right)=\left\{\begin{array}{c}\frac{1}{\mathrm{\Gamma }\left(m-\delta \right)}{\int }_0^t{\left(t-\zeta \right)}^{m-\delta -1}{f}^{\left(m\right)}\left(\zeta \right)d\zeta ,m-1<\delta <m,0\le \zeta <t,\\ f^{\left(m\right)}\left(t\right),\delta =m,\end{array}\right.(2)$

where $\left.\delta \in (m,m-1\right]$, $D_t^{\delta }$ is the CFD operator and $\mathrm{\Gamma }\left(\bullet \right)$ is the gamma function.

Moreover, the numerical results in this study are contrasted using pure ${}^{239}Pu$ for the first time in studying spherical and hemispherical neutron diffusion equations, whereas the previous works in history studied pure ${}^{238}U$.

The following section presents the methodology of the RPSM and its application to spheres and hemispheres. The radiation boundary condition is also covered in depth to set the stage for the computations. Section 3 presents numerical results for the flux distribution in the studied nuclear reactors and critical radius calculations. The associated technical problems encountered throughout the computations are presented, along with the corresponding solutions.

2 Theory

2.1 The basic idea of the RPSM

RPSM is a technique to find the power series solution coefficients without having a recurrence relation. Indeed, the coefficients of the power series are calculated by the concept of residual error through a sequence of algebraic equations, and in the end, a truncated series solution (approximate solution) is obtained. The major merit of RPSM is that it can be implemented to the problem directly without linearization, perturbation, or discretization and without any transformation by selecting appropriate initial conditions.

Most of the problems resolved using RPSM were to provide series solutions about regular points and did not address situations where there are singular points. In this work, we introduce a simple modification to the construction of the RPSM to be able to solve issues containing these types of points.

The basic definition in addition to the basic theories of RPSM and its applicability to different types of differential equations is given in [21–29]. For the convenience of the reader, we will introduce a review of the RPSM by clarifying the following algorithm on an ordinary differential equation:

(1) Express the $n$th-order ordinary differential equation in the following general form:

${{\left(x-x_0\right)}^n\varphi }^{\left(n\right)}\left(x\right)=N_x\left[\varphi \right],(3)$

subject to the initial conditions

${\varphi }^{\left(k\right)}\left(x_0\right)={\lambda }_k,k=\mathrm{0,1,2},\dots ,n-1,(4)$

where $\varphi$ is an unknown function of x on an open interval $I$ containing $x_0$ and $N_x$ is a nonlinear differential operator of order $j\le n-1$, and its coefficients are polynomials.

(2) If $x_0$ is a regular point for Eq. 3, assume the solution has a power series expansion:

$\varphi \left(x\right)={\displaystyle \sum _{k=0}^{\infty }}{c}_k{\left(x-x_0\right)}^k,(5)$

whereas if $x_0$ is a singular point for Eq. 3, assume the solution has the following expansion (this is the new case that the RPSM did not address):

$\varphi \left(x\right)={\displaystyle \sum _{i=0}^{\infty }}{c}_k{\left(x-x_0\right)}^{k+\delta },c_0\ne 0.(6)$

(3) Define the residual functions as follows:

$\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)={\left(x-x_0\right)}^n{\varphi }^{\left(n\right)}\left(x\right)-N_x\left[\varphi \right].(7)$

Indeed, $\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)=0$, for every $x\in I$ and so $\frac{d^j}{d{x}^j}\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)=0,j=\mathrm{0,1,2},\dots$ at every interior point of $I$.

(4) Substitute the expansions (5) and (6) of $\varphi \left(x\right)$ into Eq. 7 to obtain, respectively, the following expansions:

$\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)={\displaystyle \sum _{k=n}^{\infty }}\frac{c_{k}k!}{\left(k-n\right)!}{\left(x-x_0\right)}^k-N_x\left[{\displaystyle \sum _{k=0}^{\infty }}{c}_k{\left(x-x_0\right)}^k\right],(8)$

$\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{c_k\mathrm{\Gamma }\left(k+\delta +1\right)}{\mathrm{\Gamma }\left(k+\delta -n+1\right)}{\left(x-x_0\right)}^{k+\delta }-N_x\left[{\displaystyle \sum _{k=0}^{\infty }}{c}_k{\left(x-x_0\right)}^k\right],(9)$

where $\mathrm{\Gamma }$ is the gamma function.

(5) In the singular point case, we need to find the value of $\delta \ge 0$. Therefore, we find the CFD of order $\delta$ for the residual function in Eq. 9 as follows:

$D_{x_0}^{\delta }\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)={\displaystyle \sum _{k=0}^{\infty }}\frac{c_k{\left(\mathrm{\Gamma }\left(k+\delta +1\right)\right)}^2}{\mathrm{\Gamma }\left(k+\delta -n+1\right)k!}{\left(x-x_0\right)}^k-D_{x_0}^{\delta }{N}_x\left[{\displaystyle \sum _{k=0}^{\infty }}{c}_k{\left(x-x_0\right)}^k\right],(10)$

where $D_{x_0}^{\delta }$ is the CFD of order $\delta \ge 0$, which has the following two properties: 1) $D_{x_0}^{\delta }\left(a\right)=0,a$ is a constant 2) $D_{x_0}^{\delta }{\left(x-x_0\right)}^{\gamma }=\frac{\mathrm{\Gamma }\left(\gamma +1\right)}{\mathrm{\Gamma }\left(\gamma -\delta +1\right)}{\left(x-x_0\right)}^{\gamma -\delta }$.

(6) Since the CFD of a constant is 0, the value of $\delta$ is the solution of the following indicial equation:

$D_{x_0}^{\delta }\mathrm{R}\mathrm{e}\mathrm{s}\left(x_0\right)=0.(11)$

Since $\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)=0$ and $D_{x_0}^{\delta }\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)=0$ for all $x$ in the interior of $I$, the following equations are true:

$\frac{d^i}{d{x}^i}\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)=0,i=\mathrm{0,1,2},\dots ,(12)$

$\frac{d^i}{d{x}^i}{D}_{x_0}^{\delta }Res\left(x\right)=0,i=\mathrm{0,1,2},\dots .(13)$

(7) If $x_0$ is a regular point for Eq. 3, then $c_k={\lambda }_k,k=\mathrm{0,1,2},\dots ,n-1$, and $c_k$ for $k=n,n+1,n+2,\dots ,$ can be determined by solving the following algebraic equations, sequentially:

$\frac{d^{k-n}}{d{x}^{k-n}}Res\left(x_0\right)=0,k=n,n+1,n+2,\dots ,(14)$

whereas if $x_0$ is a singular point for Eq. 3, then $c_0\ne 0$ is arbitrary and can be determined later by the given initial conditions, and $c_k$ for $k=\mathrm{1,2,3},\dots$ can be specified by solving the following algebraic equations, sequentially:

$\frac{d^k}{d{x}^k}{D}_{x_0}^{\delta }Res\left(x_0\right)=0,k=\mathrm{1,2,3},\dots .(15)$

(8) Replace the values obtained for $c_k$, $k=\mathrm{0,1,2},\dots m$ with their place in the expanded series of $\varphi \left(x\right)$ to obtain the $m$th approximate solution to Eq. 3.

(9) If there is a pattern in the coefficients of the series in terms of some well-known elementary functions, then we have the exact solution $\varphi \left(x\right)$.

2.2 Bare sphere

First, the application of RPSM in spherical reactors is studied. We consider first the rather simple example of a bare spherical reactor of radius $a$. The spherical symmetry of the system implies that the flux is a function of $r$ only. The time-independent diffusion Eq. 1 can thus be written as

${\nabla }^2\varphi \left(r\right)+B^2\varphi \left(r\right)=0,(16)$

where

$B^2=\frac{{\upsilon \sum }_f-{\sum }_a}{D}(17)$

is the buckling of the reactor. Substituting the Laplacian in spherical coordinates and using $x=Br$, we obtain

$x^2\frac{d^2\varphi \left(x\right)}{d{x}^2}+2x\frac{d\varphi \left(x\right)}{dx}+x^2\varphi \left(x\right)=0.(18)$

To solve Eq. 18 using RPSM, we assume that the solution has a fractional power series representation at $x=0$ as follows:

$\varphi \left(x\right)={\displaystyle \sum _{n=0}^{\infty }}{c}_n{x}^{n+\delta },\hspace{1em}{c}_0\ne 0.(19)$

Define the residual function of Eq. 18 as

$\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)=x^2\frac{d^2\varphi \left(x\right)}{d{x}^2}+2x\frac{d\varphi \left(x\right)}{dx}+x^2\varphi \left(x\right).(20)$

Substitute Eq. 19 into the residual function (20), we obtain the following series:

$\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)={\displaystyle \sum _{n=0}^{\infty }}\left(n+\delta \right)\left(n+\delta -1\right)c_n{x}^{n+\delta }+{\displaystyle \sum _{n=0}^{\infty }}2\left(n+\delta \right)c_n{x}^{n+\delta }+{\displaystyle \sum _{n=0}^{\infty }}{c}_n{x}^{n+\delta +2}.(21)$

According to Eq. 10, the CFD of the residual function in Eq. 21 will be as

$D_{x_0}^{\delta }\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)={\displaystyle \sum _{n=0}^{\infty }}\left(n+\delta \right)\left(n+\delta -1\right)\frac{\mathrm{\Gamma }\left(n+\delta +1\right)}{n!}{c}_n{x}^n+{\displaystyle \sum _{n=0}^{\infty }}2\left(n+\delta \right)\frac{\mathrm{\Gamma }\left(n+\delta +1\right)}{n!}{c}_n{x}^n+{\displaystyle \sum _{n=0}^{\infty }}\frac{\mathrm{\Gamma }\left(n+\delta +3\right)}{\left(n+2\right)!}{c}_n{x}^{n+2}.(22)$

The solution of the indicial equations $D_{x_0}^{\delta }\mathrm{R}\mathrm{e}\mathrm{s}\left(0\right)=0$ gives $\delta =0$. Thus, Eq. 22 becomes

$D_{x_0}^{\delta }\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)=\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)={\displaystyle \sum _{n=0}^{\infty }}n\left(n-1\right)c_n{x}^n+{\displaystyle \sum _{n=0}^{\infty }}2n{c}_n{x}^n+{\displaystyle \sum _{n=0}^{\infty }}{c}_n{x}^{n+2}.(23)$

Thus, the $j$th derivative of the residual function $\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)$ has the following general form:

$\frac{d^j\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)}{d{x}^j}=\left\{\begin{array}{l}{\displaystyle \sum _{n=1}^{\infty }}{n}^2\left(n-1\right)c_n{x}^{n-1}+{\displaystyle \sum _{n=1}^{\infty }}2n^2{c}_n{x}^{n-1}+{\displaystyle \sum _{n=0}^{\infty }}\left(n+2\right)c_n{x}^{n+1},j=1,\\ {\displaystyle \sum _{n=2}^{\infty }}{n}^2{\left(n-1\right)}^2{c}_n{x}^{n-2}+{\displaystyle \sum _{n=2}^{\infty }}2n^2\left(n-1\right)c_n{x}^{n-2}+{\displaystyle \sum _{n=0}^{\infty }}\left(n+2\right)\left(n+1\right)c_n{x}^n,j=2,\\ {\displaystyle \sum _{n=j}^{\infty }}\frac{n\left(n-1\right)n!}{\left(n-j\right)!}{c}_n{x}^{n-j}+{\displaystyle \sum _{n=j}^{\infty }}\frac{2n\left(n!\right)}{\left(n-j\right)!}{c}_n{x}^{n-j}+{\displaystyle \sum _{n=j-2}^{\infty }}\frac{\left(n+2\right)!c_n}{\left(n-j+2\right)!}{x}^{n-j+2},j=\mathrm{3,4},\dots .\end{array}\right.(24)$

According to step (9) in Subsection 2.1, the solution of the algebraic equations

${\left.\frac{d^j\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)}{d{x}^j}\right|}_{x=0}=0,j=\mathrm{1,2},\dots (25)$

can be summarized in Table 1.

TABLE 1

TABLE 1. Some values of the coefficients of the expansion (19).

If we continue in the same manner, then we conclude that the general form of the coefficients in series (19) will be in the following form:

$c_n=\left\{\begin{array}{l}\frac{{\left(-1\right)}^{\frac{n}{2}}{c}_0}{\left(n+1\right)!},n=\mathrm{0,2,4,6},\dots ,\\ 0,n=\mathrm{1,3,5},\dots .\end{array}\right.(26)$

Therefore, the solution for Eq. 18 has the following exact solution:

$\varphi \left(x\right)={\displaystyle \sum _{n=0}^{\infty }}\frac{{\left(-1\right)}^n{c}_0}{\left(2n+1\right)!}{x}^{2n}=\frac{c_0}{x}\sin x.(27)$

Consequently, the final solution for Eq. 16, in a closed form, is given by

$\varphi \left(Br\right)=\frac{c_0\mathit{sin}\left(Br\right)}{Br}.(28)$

This flux must satisfy the relevant boundary conditions.

2.3 Hemispherical reactor

In this subsection, we will consider the hemispherical reactor. This case is more complicated because the flux is a function in two variables $r$ and $\theta$. Therefore, the time-independent diffusion Eq. 1 can be expressed, by writing the Laplacian in the spherical coordinates and using $\mu =\cos \theta$, as follows:

$r^2\frac{{\partial }^2\varphi \left(r,\mu \right)}{\partial r^2}+2r\frac{\partial \varphi \left(r,\mu \right)}{\partial r}+\frac{\partial }{\partial \mu }\left(\left(1-{\mu }^2\right)\frac{\partial \varphi \left(r,\mu \right)}{\partial \mu }\right)+B^2\varphi \left(r,\mu \right)=0,(29)$

where the buckling $B^2$ is again given by Eq. 17.

2.3.1 Solution by the separation of variables and residual power series method

Applying the separation of variables, we obtain

$\varphi \left(r,\mu \right)=R\left(r\right)\psi \left(\mu \right).(30)$

We obtain the eigenvalue equations as follows:

$r^2\frac{d^{2}R\left(r\right)}{d{r}^2}+2r\frac{dR\left(r\right)}{dr}+\left(B^2{r}^2-m\left(m+1\right)\right)R\left(r\right)=0,(31)$

$\frac{d}{d\mu }\left(\left(1-{\mu }^2\right)\frac{d\psi \left(\mu \right)}{d\mu }\right)+m\left(m+1\right)\psi \left(\mu \right)=0.(32)$

To solve the radial part using RPSM, we first, for simplicity, consider $x=Br$ and rewrite Eq. 31 as

$x^2\frac{d^{2}R\left(x\right)}{d{x}^2}+2x\frac{dR\left(x\right)}{dx}+\left(x^2-m\left(m+1\right)\right)R\left(x\right)=0,(33)$

and we assume the solution has a power series representation at $x=0$,

$R\left(x\right)={\displaystyle \sum _{n=0}^{\infty }}{c}_n{x}^{n+\delta }.(34)$

The residual function of Eq. 33 has the following form:

$Res\left(x\right)=x^2\frac{d^{2}R\left(x\right)}{d{x}^2}+2x\frac{dR\left(x\right)}{dx}+x^{2}R\left(x\right)-m\left(m+1\right)R\left(x\right).(35)$

If we substitute Eq. 34 into the residual function (35), then we obtain the following series function:

$Res\left(x\right)={\displaystyle \sum _{n=0}^{\infty }}{\left(n+\delta \right)\left(n+\delta -1\right)c}_n{x}^{n+\delta }+{\displaystyle \sum _{n=0}^{\infty }}2\left(n+\delta \right)c_n{x}^{n+\delta }+{\displaystyle \sum _{n=2}^{\infty }}{c}_{n-2}{x}^{n+\delta }-m\left(m+1\right){\displaystyle \sum _{n=0}^{\infty }}{c}_n{x}^{n+\delta },(36)$

and so, the CFD of the residual function of order $\delta$ in Eq. 36 has the following expansion:

$D_{x_0}^{\delta }\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)={\displaystyle \sum _{n=0}^{\infty }}{\left(n+\delta \right)\left(n+\delta -1\right)\frac{\mathrm{\Gamma }\left(n+\delta +1\right)}{n!}c}_n{x}^n+{\displaystyle \sum _{n=0}^{\infty }}2\left(n+\delta \right)c_n\frac{\mathrm{\Gamma }\left(n+\delta +1\right)}{n!}{x}^n+{\displaystyle \sum _{n=2}^{\infty }}\frac{\mathrm{\Gamma }\left(n+\delta +1\right)}{n!}{c}_{n-2}{x}^n-m\left(m+1\right){\displaystyle \sum _{n=0}^{\infty }}\frac{\mathrm{\Gamma }\left(n+\delta +1\right)}{n!}{c}_n{x}^n.(37)$

Substituting $x=0$ in Eq. 37 gives the indicial equation of series (34), whose solution is $\delta =m$. Hence, the residual function in Eq. 37 can be rewritten as

$D_{x_0}^m\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)={\displaystyle \sum _{n=0}^{\infty }}{\left(n+m\right)\left(n+m-1\right)\frac{\left(n+m\right)!}{n!}c}_n{x}^n+{\displaystyle \sum _{n=0}^{\infty }}2\left(n+m\right)\frac{\left(n+m\right)!}{n!}{c}_n{x}^n+{\displaystyle \sum _{n=2}^{\infty }}\frac{\left(n+m\right)!}{n!}{c}_{n-2}{x}^n-m\left(m+1\right){\displaystyle \sum _{n=0}^{\infty }}\frac{\left(n+m\right)!}{n!}{c}_n{x}^n.(38)$

Following the same manner as in Section 2.2, if we apply the sequential operator $\frac{d^j}{d{x}^i}$, $j=\mathrm{1,2,3},\dots$ on both sides of Eq. 38, we obtain

$\frac{d^j{D}_{x_0}^m\mathrm{R}\mathrm{e}\mathrm{s}\left(x\right)}{d{x}^j}=\left\{\begin{array}{l}{\displaystyle \sum _{n=1}^{\infty }}\frac{\left(n+m\right)!}{\left(n-1\right)!}\left(n\left(2m+n+1\right)\right)c_n{x}^{n-1}+{\displaystyle \sum _{n=2}^{\infty }}\frac{\left(n+m\right)!}{\left(n-1\right)!}{c}_{n-2}{x}^{n-1},j=1,\\ {\displaystyle \sum _{n=2}^{\infty }}\frac{\left(n+m\right)!}{\left(n-2\right)!}\left[n\left(2m+n+1\right)c_n+c_{n-2}\right]x^{n-2},j=2,\\ {\displaystyle \sum _{n=j}^{\infty }}\frac{\left(n+m\right)!}{\left(n-j\right)!}\left[n\left(2m+n+1\right)c_n+c_{n-2}\right]x^{n-j},j=\mathrm{3,4},\dots .\end{array}\right.(39)$

Consequently, the first few coefficients of the series in Eq. 34 are given in Table 2.

TABLE 2

TABLE 2. Some values of the coefficients of the expansion (34).

It is clear that there is a pattern in the coefficients of our series. So, the $n$th coefficient, $c_n$, has a general form as follows:

$c_n=\left\{\begin{array}{l}\frac{{\left(-1\right)}^{\frac{n}{2}}{c}_0}{2\left(4\right)\dots \left(n\right)\left(2m+3\right)\left(2m+5\right)\dots \left(2m+n+1\right)},n=\mathrm{0,2,4,6},\dots ,\\ 0,n=\mathrm{1,3,5},\dots ,\end{array}\right.(40)$

which can be reformulated as follows:

$c_n=c_0\frac{{\left(-1\right)}^n\left(2m+1\right)!\left(m+n\right)!}{m!n!\left(2m+2n+1\right)!},m\in \left\{\mathrm{0,1,2,3},\dots \right\},n=\mathrm{0,1,2,3},\dots .(41)$

Therefore, the eigensolutions to the Sturm–Liouville equation in Eq. 31 appear in the following series:

$R_m\left(Br\right)={\displaystyle \sum _{n=0}^{\infty }}{C}_{mn}\frac{{\left(-1\right)}^n\left(2m+1\right)!\left(m+n\right)!}{m!n!\left(2m+2n+1\right)!}{\left(Br\right)}^{2n+m}.(42)$

Now, to solve the angular part, which is a Legendre equation, we construct the residual function of Eq. 32 as

$Res\left(\mu \right)=\frac{d^2\psi \left(\mu \right)}{d{\mu }^2}-{\mu }^2\frac{d^2\psi \left(\mu \right)}{d{\mu }^2}-2\mu \frac{d\psi \left(\mu \right)}{d\mu }+m\left(m+1\right)\psi \left(\mu \right),(43)$

and since $\mu =0$ is a regular point of Eq. 32, we assume the solution of the following expansion:

$\psi \left(\mu \right)={\displaystyle \sum _{n=0}^{\infty }}{c}_n{\mu }^n.(44)$

Substitute Eq. 44 into Eq. 43, and we obtain

$Res\left(\mu \right)={\displaystyle \sum _{n=2}^{\infty }}n\left(n-1\right)c_n{\mu }^{n-2}-{\displaystyle \sum _{n=2}^{\infty }}n\left(n-1\right)c_n{\mu }^n-{\displaystyle \sum _{n=1}^{\infty }}2n{c}_n{\mu }^n+m\left(m+1\right){\displaystyle \sum _{n=0}^{\infty }}{c}_n{\mu }^n.(45)$

Applying the differential operator $\frac{d^j}{d{x}^i}$, $j=\mathrm{0,1,2,3},\dots$ on both sides of Eq. 45 gives

$\frac{d^{j}Res\left(x\right)}{d{x}^j}=\left\{\begin{array}{l}{\displaystyle \sum _{n=0}^{\infty }}\left[\left(n+2\right)\left(n+1\right)c_{n+2}+m\left(m+1\right)c_n\right]{\mu }^n-{\displaystyle \sum _{n=2}^{\infty }}n\left(n-1\right)c_n{\mu }^n-{\displaystyle \sum _{n=1}^{\infty }}2n{c}_n{\mu }^n,j=0,\\ {\displaystyle \sum _{n=1}^{\infty }}n\left[\left(n+2\right)\left(n+1\right)c_{n+2}+\left(m\left(m+1\right)-2n\right)c_n\right]{\mu }^{n-1}-{\displaystyle \sum _{n=2}^{\infty }}{n}^2\left(n-1\right)c_n{\mu }^{n-1},j=1,\\ {\displaystyle \sum _{n=2}^{\infty }}n\left(n-1\right)\left[\left(n+2\right)\left(n+1\right)c_{n+2}+\left(m-n\right)\left(m+n+1\right)c_n\right]{\mu }^{n-2},j=2,\\ {\displaystyle \sum _{n=j}^{\infty }}\frac{n!}{\left(n-j\right)!}\left[\left(n+2\right)\left(n+1\right)c_{n+2}+\left(m-n\right)\left(m+n+1\right)c_n\right]{\mu }^{n-j},j=\mathrm{3,4},\dots .\end{array}\right.(46)$

Thus, the first few coefficients of the series in Eq. 44 are given in Table 3.

TABLE 3

TABLE 3. Some values of the coefficients of expansion (44).

Note that there are two linearly independent solutions one of them is an infinite series and the other is a polynomial of degree m. If m is an even number, then the solution that contains terms of even exponents will be the polynomial whereas the solution that contains terms of odd exponents will be infinite series and vice versa if m is an odd number. Anyway, we can keep track of the pattern in the coefficients and write the solutions as follows:

${\psi }_{m0}\left(\mu \right)=c_{m0}\left(1+\frac{{\left(\frac{m}{2}!\right)}^2}{m!}{\displaystyle \sum _{n=1}^{\frac{m}{2}}}\frac{{\left(-1\right)}^n\left(m+2n\right)!{\mu }^{2n}}{\left(\frac{m}{2}-n\right)!\left(\frac{m}{2}+n\right)!\left(2n\right)!}\right),m=\mathrm{2,4,6},\dots ,{\psi }_{m1}\left(\mu \right)=c_{m1}\left(x+\frac{{\left(\frac{m-1}{2}!\right)}^2}{m!}\sum _{n=1}^{\frac{m-1}{2}}\frac{{\left(-1\right)}^n\left(m+2n\right)!{\mu }^{2n+1}}{\left(\frac{m-1}{2}-n\right)!\left(\frac{m-1}{2}+n\right)!\left(2n+1\right)!}\right),m=\mathrm{1,3,5},\dots .(47)$

Legendre was able to combine both solutions in a single format, regardless of whether the number $m$ was odd or even as

${\psi }_m\left(\mu \right)=c_m{P}_m\left(\mu \right)=c_m\left(\frac{1}{2^m}{\displaystyle \sum _{n=0}^{\left[m/2\right]}}\frac{{\left(-1\right)}^n\left(2m-2n\right)!{\mu }^{m-2n}}{n!\left(m-n\right)!\left(m-2n\right)!}\right).(48)$

By using Rodrigues’s formula for the Legendre polynomials, the solution can be expressed as follows:

${\psi }_m\left(\mu \right)=c_m{P}_m\left(\mu \right)=c_m\frac{1}{2^{m}m!}\frac{d^m}{d{\mu }^m}{\left({\mu }^2-1\right)}^m.(49)$

Therefore, the general solution for Eq. 29 can be expressed as

$\varphi \left(r,\mu \right)={\displaystyle \sum _{m=0}^{\infty }}{\varphi }_m\left(r,\mu \right)={\displaystyle \sum _{m=0}^{\infty }}{A}_m{R}_m\left(r\right)P_m\left(\mu \right).(50)$

2.4 The radiation boundary condition

Although the choice of suitable boundary conditions in solving the neutron diffusion equation is at the core of the nuclear reactor physicist mission, here, we will apply three boundary conditions, namely, zero flux boundary condition (ZFBC), extrapolated boundary condition (EBC), and radiation boundary condition (RBC). These boundary conditions will be applied on both spherical and hemispherical reactors.

ZFBC cures about the point at the surface of the reactor ${\mathrm{a}}_{\mathrm{c}}$ that the flux is 0, and for more reality, the flux will vanish at a point less than that point which is called extrapolated distance ${\mathrm{a}}_{\mathrm{c}}+2D$. To be more accurate, especially for the hemisphere, the mixed boundary condition gives the closest results when comparing to the benchmark; this boundary condition is called RBC $\left[n.\nabla \varphi \left(r,\mu \right)=g\left(r\right)\varphi \left(r,\mu \right)\right].$

The use of the application of RBC on the solution of the hemispherical reactor using RPSM based on the separation of variable techniques can be performed as follows:

$n.\nabla \varphi \left(r,\mu \right)=g\left(r\right)\varphi \left(r,\mu \right),(51)$

where $g\left(r\right)$ varies over the surface, and this RBC is reported in previous works [34, 35].

The RBC on the curved surface where $r=a$ will be

${\frac{\partial \varphi \left(r,\mu \right)}{\partial x}}_{r=a}=-g_1\varphi \left(a,\mu \right).(52)$

However, RBC on the flat surface will be

${\frac{\partial \varphi \left(r,\mu \right)}{\partial x}}_{\mu =0}=-g_2\varphi \left(r,0\right),(53)$

where $g_1=\frac{1}{2D},g_2=\frac{g\left(a/b\right)}{2D},$ where the values of g change from 0 to 1, and the reactor geometry changes from spherical ($g=$ 0) to hemispherical ($g=$ 1) geometry as mentioned in [34, 35].

Applying RBC on the flat surface in Eq. 50 yields

$\frac{1}{r}\sum A_m{R}_m\left(Br\right)P_m^{\prime }\left(0\right)=g_2\sum A_m{R}_m\left(Br\right)P_m\left(0\right).(54)$

Using the following recurrence relationships of $R_n\left(Br\right)$:

$\frac{R_{2m+1}\left(Br\right)}{r}=\frac{{BR}_{2m+2}\left(Br\right)}{\left(4m+5\right)\left(4m+3\right)}+B{R}_{2m}\left(Br\right),\frac{m+1}{Br}{R}_m\left(Br\right)+\frac{R_m^{\prime }\left(Br\right)}{B}=\left(2m+1\right)R_{m-1}\left(Br\right),(55)$

and some Legendre polynomials properties, we obtain the following relationship:

$B\left(2m+1\right)A_{2m+1}-\frac{2mB}{\left(4m+1\right)\left(4m-1\right)}{A}_{2m-1}=\frac{g_2}{2m+1}{A}_{2n}.(56)$

We can use the previous result to write the odd term amplitude in terms of the even term amplitude as

$A_{2m+1}={\sum }_{k=0}^n{\gamma }_{mk}{A}_{2k},(57)$

where

${\gamma }_{mk}=\frac{\left(2m\right)‼\left(2k-1\right)!!\left(4k+1\right)‼}{\left(2k\right)!!\left(2m+1\right)‼\left(4m+1\right)‼}\frac{g_2}{B},(58)$

as reported in [ 36].

Similarly, applying the RBC on the curved surface, we obtain

${\displaystyle \sum _{m=0}^{\infty }}{A}_m{R}^{\prime }\left(Ba\right)P_m\left(\mu \right)=-g_1{\displaystyle \sum _{m=0}^{\infty }}{A}_{m}R\left(Ba\right)P_m\left(\mu \right),(59)$

or it can be written as

${\displaystyle \sum _{m=0}^{\infty }}{A}_m{f}_m{P}_m\left(\mu \right)=0,(60)$

where $f_m$ $\mathrm{i}\mathrm{s}\mathrm{g}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{n}\mathrm{b}\mathrm{y}$

$f_m=\left(\frac{m}{a}+g_1\right)R_m\left(Ba\right)-\frac{B}{2m+3}{R}_{m+1}\left(Ba\right).(61)$

Following [35], we obtain

${\displaystyle \sum _{m=0}^{\infty }}{\widehat{M}}_{mk}{\mathrm{A}}_{2\mathrm{m}}=0,\mathrm{k}=\mathrm{0,1,2,3}\dots ,(62)$

where

${\widehat{M}}_{mk}={\beta }_{\left(2m,k\right)}{f}_{2m}+{\displaystyle \sum _{l=m}^{\infty }}{\beta }_{\left(2l+1,k\right)}{f}_{2l+1}{\gamma }_{l,m},(63)$

and

${\beta }_{mk}=\left(2k+1\right){\int }_0^1{P}_m\left(\mu \right)P_k\left(2\mu -1\right)d\mu .(64)$

The alternating method of ${\beta }_{mk}$ is

${\beta }_{mk}=P_m\left({\mu }_k\right),(65)$

where ${\mu }_k=\cos \left(k\pi /2N\right)$, which is used in these calculations.

3 Numerical results

The numerical data of neutrons diffusing in pure ${}^{239}Pu$ are taken from [36], which are shown in Table 4.

TABLE 4

TABLE 4. Pure ${}^{239}Pu$ cross sections.

The important diffusion constant and buckling are

$D=\frac{1}{3\left({\mathrm{\Sigma }}_f+{\mathrm{\Sigma }}_s+{\mathrm{\Sigma }}_{\gamma }\right)}=1.021241834cm,B=\sqrt{\frac{ʋ{\mathrm{\Sigma }}_f-\left({\mathrm{\Sigma }}_f+{\mathrm{\Sigma }}_{\mathrm{\gamma }}\right)}{D}}=0.357553294{cm}^{-1}.$

Depending on these data, to calculate the flux, applying different boundary conditions, to find the critical radii, and to study the flux distribution for both spherical and hemispherical reactors,, this work has been compared with benchmark calculations [34].

3.1 Critical radius calculations

The critical radius $a_c$ is the dimension of the reactor determined as the flux vanished. The determination of the critical radius depends on the boundary conditions, and the simplest way of calculating it is by using ZFBC, where the flux vanishes. In real reactors, an extrapolated distance must be considered ($a_c$ +2D) at which the value of the flux is 0; this boundary condition is the EBC. The use of ZFBC and EBC is suitable for sphere, cylinder, and slab reactors where there is high symmetry, but for more complicated reactors like hemispherical and cone reactors, we need to use the RBC.

The first case we will study is the spherical reactor. Here, the critical radius using ZFBC and EBC is considered, and the results are compared with RBC. A $C^{++}$ computer code based on Eq. 28 has been written to find the first root of Eq. 28. The value of $r$ at which the flux is 0 is $a_{c,ZF}=8.7864$ cm. Since the flux is assumed to converge at $a_c+2D$ using EBC, then $a_{c,EBC}=8.7864$. These results are tabulated in Table 5.

TABLE 5

TABLE 5. Critical radius $a_c$ (in cm) for pure ${}^{239}Pu$ for a sphere and hemisphere using different boundary conditions: ZFBC, EBC, and RBC.

After taking the bare sphere, the hemispherical reactor is considered, and the neutron flux using ZFBC will vanish on the flat surface (µ = 0). Using the Legendre polynomial properties, the even amplitudes must be equal to 0, and on the curved surface, we need to find the first zero of the solution $R_n\left(Br\right)$ of the radial part. The unique available case is for $m=1$. So Eq. 50 will reduce to the following form:

$\varphi \left(r,\mu \right)=A_1{R}_1\left(Br\right){\psi }_1\left(\mu \right)={\displaystyle \sum _{m=0}^{\infty }}\frac{{c_m\left(-1\right)}^n\left(3\right)!\left(m+1\right){\left(Br\right)}^{1+2m}}{\left(3+2m\right)!}\mu .(66)$

The critical radii in pure ${}^{239}Pu$ that are calculated using the $C^{++}$ code depending on Eq. 59 are tabulated in Table 5.

Unlike the spherical reactor, applying RBC is essential for the hemispherical reactor.

The calculation of the critical dimension using RBC is illustrated in the following equation, and the infinite summation in Eq. 62 will be cut off to the value N, where the dimension of the matrix is finite, and rewritten as

${\displaystyle \sum _{m=0}^N}{\widehat{M}}_{mk}{A}_{2m}=0,(67)$

where

${\widehat{M}}_{mk}={\beta }_{\left(2m,k\right)}{f}_{2m}+{\displaystyle \sum _{l=n}^N}{\beta }_{\left(2l+1,k\right)}{f}_{2l+1}{\gamma }_{l,m.}(68)$

This simplification is an important numerical calculation to make cut-off limit value N. Then, the flux in Eq. 50 will be simplified to

$\varphi \left(r,\mu \right)={\displaystyle \sum _{m=0}^{2N+1}}{\varphi }_m\left(r,\mu \right)={\displaystyle \sum _{m=0}^{2N+1}}{A}_m{R}_m\left(r\right){\psi }_m\left(\mu \right).(69)$

A third computer code has been used to calculate the critical radius using RBC. The code calculates the determinant of ${\widehat{M}}_{nk}$. The value of a affects the determinant through the parameters $f_n$ (Eq. 61). As $a$ increases, the value of the determinant will change its sign. On the other hand, the calculated $a_c$ depends on the value of the interaction parameter g through ${\gamma }_{nk}$ (Eq. 58). Hence, the results for the sphere and the infinitely separated hemispheres are obtained by using $g=0$ and $g=1$, respectively.

The determination of the critical dimension depends on the cut-off value N, where the number of terms will be $2N+2$, and the value of the critical dimension with the number of terms is given in Table 6.

TABLE 6

TABLE 6. Convergence of the critical radius $a_c$ (in cm) for the hemisphere using RBC.

The effect of the number of terms is clear; the value of the critical radius is converged by increasing the number of terms. Now, the effect of changing the interaction parameter $g$ on determining the critical radius is clarified in Table 7, where the cut-off value $N=10$ and the number of terms (22) are considered.

TABLE 7

TABLE 7. Variation in the critical radius $a_c$ with the interaction parameter $g$. $N=10$ (22 terms).

The reactor is a sphere when the interaction parameter $g=0$, whereas it is a hemisphere when the interaction parameter $g=1$, where the values of the critical radii are shown in Table 4 and Table 5, respectively.

3.2 Flux distribution

Now, the flux destruction inside the hemispherical reactor will be considered. It is straightforward that the symmetry of the spherical reactor makes the flux distribution a simple case study in [35], so the flux in Eq. 69 will be written as

$\varphi \left(r,\mu \right)={\displaystyle \sum _{m=0}^{2N+1}}\left[A_{2m}{R}_{2m}\left(Br\right){\psi }_{2m}\left(\mu \right)+A_{2m+1}{R}_{2m+1}\left(Br\right){\psi }_{2m+1}\left(\mu \right)\right].(70)$

For computational calculation purpose, this equation can be expressed as

$\overline{\varphi }\left(r,\mu \right)={\displaystyle \sum _{m=0}^N}{A}_{2m}\left[R_{2m}\left(Br\right){P\psi }_{2m}\left(\mu \right)+{\displaystyle \sum _{l=m}^N}{\gamma }_{l,n}{R}_{2m+1}\left(Br\right){\psi }_{2m+1}\left(\mu \right)\right].(71)$

The amplitudes $A_{2m}$ in Eq. 62 are needed to be written in the form

${\displaystyle \sum _{n=0}^N}{\widehat{M}}_{nk}{A}_{2m}=-{\widehat{M}}_{0k}{A}_0,(72)$

where $A_1$ has been normalized to unity. So A₀ $=\frac{1}{{\gamma }_{00}}$, and the flux values are normalized to the averaged flux.

$\overline{\varphi }=\frac{{\int }_0^a{\int }_0^1{r}^2\varphi \left(r,\mu \right)d\mu dr}{{\int }_0^a{\int }_0^1{r}^{2}d\mu dr}.(73)$

The value of the average flux is

$\overline{\varphi }=\frac{3}{a^3}{\int }_0^a{r}^{2}dr{\int }_0^1\varphi \left(r,\mu \right)d\mu .(74)$

Based on these data, the averaged flux calculated using the $C^{++}$ code was determined to be 0.756965. The variation in the flux at $\theta =\pi /4$ with $r$ changes from zero, and the critical radius is given in Table 8 for different values of N (different number of terms).

TABLE 8

TABLE 8. Flux variation across the critical hemisphere at $\theta =\pi /4$

The convergence of the flux values is clearly given in Table 8.

As shown in Figure 1, the fast convergence of the flux on the curved surface is obtained with the increasing number of terms (for $N=\mathrm{1,2,3,4},$ and 30).

FIGURE 1

FIGURE 1. Angular distribution over the curved surface.

It is important to clarify, from Figure 1, that the convergence is stable and acceptable, and when the number of terms is 30, the angular flux, which expresses the behavior of neutrons when they are diffusing in pure ${}^{239}Pu$, is expected which is the same as in using the homotopy perturbation method in solving the neutron diffusion equation where neutrons diffuse in the ${}^{235}U$ hemispherical reactor [35].

After studying the effect of the number of terms on the flux, it is taken into consideration at a wide range of angle values in Figure 2.

FIGURE 2

FIGURE 2. Flux distribution across the core of a bare hemisphere at different angles $N=10\left(22\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{s}\right)$

As shown in Figure 2, it is noted that the radial flux does not decrease when $r$ increases and the maximum value of the flux shifted with the increase in the angle. One more time, the neutron distribution is in the same manner as in previous work [35].

To clarify the flux behavior in more detail, the flux values for different angles and radii are tabulated in Table 9 and depicted in Figure 3.

TABLE 9

TABLE 9. Flux variation across the critical hemisphere at different values of $\theta$. $N=10$ (22 terms).

FIGURE 3

FIGURE 3. Flux distribution across the core of a bare hemisphere as a function of $r$ and $\theta$. $N=10(22$ terms).

Furthermore, a three-dimensional plot of the flux distribution across the core of a bare hemisphere as a function of $r$ and $\theta$ is shown in Figure 3.

The flux in Table 6 is represented graphically in Figure 3. The data in Figure 1 and Figure 2 are merged in the same figure, where the full neutron behavior is represented; however, this three-dimensional flux representation is not figured out previously for the hemispherical reactor.

Although the neutrons diffusing in pure ${}^{239}Pu$ in this paper are studied using RPSM for the first time and the neutrons diffusing in ${}^{235}U$ were previously studied using the homotopy perturbation method [35, 36], it is clear that their behavior is the same.

4 Conclusion

In this article, we introduced the application of RPSM in solving the neutron diffusion equation in hemispherical symmetry. The method is efficient and can present the solution in the form of a rapid series that converge to the exact solution. Two important cases are discussed: the spherical reactor and hemispherical reactor. The results are presented in tables and illustrative figures using $C++$ codes. To show the validity of the presented method, we made qualitative comparisons between the obtained results using RPSM and those computed using the homotopy perturbation method. The used method and its computational implementation provide superior results when applied to more complex and reflected reactor geometries for one-group problems. It is recommended, however, for future work to further investigate the methodology for two- and multi-group problems.

In future work, we attend to solve more partial differential equations with physical applications, using the proposed method that showed its accuracy in handling similar problems, and further investigation of the convergence-produced series should be the goal of future specialized research [37–41].

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Krane KS. Introductory nuclear physics. John Wiley & Sons (1991).

Google Scholar

2. Duderstadt JJ, Hamilton LJ. Nuclear reactor analysis. Boston, MA, USA: John Wiley & Sons (1976).

Google Scholar

3. Lamarsh JR. Introduction to nuclear engineering. 2nd ed. Boston, MA, USA: Addison-Wesley (1983).

Google Scholar

4. Kochenov AS. The stability of a nuclear power plant. J Nucl Energ A. Reactor Sci (1960) 12(4):194–200. doi:10.1016/0368-3265(60)90103-1

CrossRef Full Text | Google Scholar

5. Dababneh S, Khasawneh K, Odibat Z. An alternative solution of the neutron diffusion equation in cylindrical symmetry. Ann Nucl Energ (2011) 38(5):1140–3. doi:10.1016/j.anucene.2010.12.011

CrossRef Full Text | Google Scholar

6. Nairat M, Shqair M, Alhalholy T. Cylindrically symmetric fractional Helmholtz equation. Appl Math E-notes (2019) 19:708–17.

Google Scholar

7. Shqair M. Developing a new approaching technique of homotopy perturbation method to solve two-group reflected cylindrical reactor. Results Phys (2019) 12:1880–7. doi:10.1016/j.rinp.2019.01.063

CrossRef Full Text | Google Scholar

8. Shqair M, Farrag EA, Al-Smadi M. Solving multi-group reflected spherical reactor system of equations using the homotopy perturbation method. Mathematics (2022) 10(10):1784. doi:10.3390/math10101784

CrossRef Full Text | Google Scholar

9. Ala’yed O, Saadeh R, Qazza A. Numerical solution for the system of Lane-Emden type equations using cubic B-spline method arising in engineering. AIMS Maths (2023) 8(6):14747–66. doi:10.3934/math.2023754

CrossRef Full Text | Google Scholar

10. Edwan R, Saadeh R, Hadid S, Al-Smadi M, Momani S. Solving time-space-fractional Cauchy problem with constant coefficients by finite-difference method. Comput Math Appl (2020) 25–46.

CrossRef Full Text | Google Scholar

11. Shqair M, Ghabar I, Burqan A. Using Laplace residual power series method in solving coupled fractional neutron diffusion equations with delayed neutrons system. Fractal Fract (2023) 7:219. doi:10.3390/fractalfract7030219

CrossRef Full Text | Google Scholar

12. Aboanber AE, Nahla AA, Aljawazneh SM. Fractional two energy groups matrix representation for nuclear reactor dynamics with an external source. Ann Nucl Energ (2021) 153:108062. doi:10.1016/j.anucene.2020.108062

CrossRef Full Text | Google Scholar

13. Sardar T, Saha Ray S, Bera RK, Biswas BB, Das S. The solution of coupled fractional neutron diffusion equations with delayed neutrons. Int J Nucl Energ Sci Technol (2010) 5(2):105–13. doi:10.1504/ijnest.2010.030552

CrossRef Full Text | Google Scholar

14. El-Nabulsi RA. Fractal neutrons diffusion equation: Uniformization of heat and fuel burn-up in nuclear reactor. Nucl Eng Des (2021) 380:111312. doi:10.1016/j.nucengdes.2021.111312

CrossRef Full Text | Google Scholar

15. El-Nabulsi RA. Nonlocal effects to neutron diffusion equation in a nuclear reactor. J Comput Theor Transport (2020) 49(6):267–81. doi:10.1080/23324309.2020.1816551

CrossRef Full Text | Google Scholar

16. El-Nabulsi RA. Neutrons diffusion variable coefficient advection in nuclear reactors. Int J Adv Nucl Reactor Des Technol (2021) 3:102–7. doi:10.1016/j.jandt.2021.06.005

CrossRef Full Text | Google Scholar

17. Patra A, Shone TT, Mishra BB. Natural decomposition approximation solution for first order nonlinear differential equations. Int J Eng Technol (2018) 7(4):442–5. doi:10.14419/ijet.v7i4.5.20202

CrossRef Full Text | Google Scholar

18. Nahla AA, Al-Malki FA, Rokaya M. Numerical techniques for the neutron diffusion equations in the nuclear reactors. Adv Stud Theor Phys (2012) 6(14):649–64.

Google Scholar

19. Tumelero F, Lapa CM, Bodmann BE, Vilhena MT. Analytical representation of the solution of the space kinetic diffusion equation in a one-dimensional and homogeneous domain. Braz J Radiat Sci (2019) 7(2B):7. doi:10.15392/bjrs.v7i2b.389

CrossRef Full Text | Google Scholar

20. Khaled SM. Exact solution of the one-dimensional neutron diffusion kinetic equation with one delayed precursor concentration in Cartesian geometry. AIMS Math (2022) 7:12364–73. doi:10.3934/math.2022686

CrossRef Full Text | Google Scholar

21. El-Ajou A, Al-Smadi M, Oqielat MN, Momani S, Hadid S. Smooth expansion to solve high-order linear conformable fractional PDEs via residual power series method: Applications to physical and engineering equations. Ain Shams Eng J (2020) 11(4):1243–54. doi:10.1016/j.asej.2020.03.016

CrossRef Full Text | Google Scholar

22. El-Ajou A. Taylor’s expansion for fractional matrix functions: Theory and applications. J Math Comput Sci (2020) 21(1):1–17. doi:10.22436/jmcs.021.01.01

CrossRef Full Text | Google Scholar

23. Oqielat MN, Eriqat T, Al-Zhour Z, Ogilat MN, El-Ajou A, Hashim I. Construction of fractional series solutions to nonlinear fractional reaction–diffusion for bacteria growth model via Laplace residual power series method. Int J Dyn Control (2023) 11(2):520–7. doi:10.1007/s40435-022-01001-8

CrossRef Full Text | Google Scholar

24. Oqielat MN, Eriqat T, Ogilat O, El-Ajou A, Alhazmi SE, Al-Omari S. Laplace-residual power series method for solving time-fractional reaction–diffusion model. Fractal and Fractional (2023) 7(4):309. doi:10.3390/fractalfract7040309

CrossRef Full Text | Google Scholar

25. Saadeh R, Burqan A, El-Ajou A. Reliable solutions to fractional Lane-Emden equations via Laplace transform and residual error function. Alexandria Eng J (2022) 61(12):10551–62. doi:10.1016/j.aej.2022.04.004

CrossRef Full Text | Google Scholar

26. Saadeh R. A reliable algorithm for solving system of multi-pantograph equations. WSEAS Trans Math (2022) 21:792–800. doi:10.37394/23206.2022.21.91

CrossRef Full Text | Google Scholar

27. Khresat H, El-Ajou A, Al-Omari S, Alhazmi SE, Oqielat MA. Exact and approximate solutions for linear and nonlinear partial differential equations via Laplace residual power series method. Axioms (2023) 12(7):694. doi:10.3390/axioms12070694

CrossRef Full Text | Google Scholar

28. Korpinar Z, Baleanu D, İnç M, Almohsen B. Some applications of the least squares-residual power series method for fractional generalized long wave equations. J Ocean Eng Sci (2021). doi:10.1016/j.joes.2021.09.001

CrossRef Full Text | Google Scholar

29. Salah E, Qazza A, Saadeh R, El-Ajou A. A hybrid analytical technique for solving multi-dimensional time-fractional Navier-Stokes system. AIMS Maths (2023) 8(1):1713–36. doi:10.3934/math.2023088

CrossRef Full Text | Google Scholar

30. El-Ajou A, Al-Zhour Z. A vector series solution for a class of hyperbolic system of Caputo-time-fractional partial differential equations with variable coefficients. Front Phys (2021) 9:276. doi:10.3389/fphy.2021.525250

CrossRef Full Text | Google Scholar

31. Burqan A, El-Ajou A, Saadeh R, Al-Smadi M. A new efficient technique using Laplace transforms and smooth expansions to construct a series solution to the time-fractional Navier-Stokes equations. Alexandria Eng J (2022) 61(2):1069–77. doi:10.1016/j.aej.2021.07.020

CrossRef Full Text | Google Scholar

32. El-Ajou A. A modification to the conformable fractional calculus with some applications. Alexandria Eng J (2020) 59(4):2239–49. doi:10.1016/j.aej.2020.02.003

CrossRef Full Text | Google Scholar

33. El-Nabulsi RA, Torres DF. Necessary optimality conditions for fractional action-like integrals of variational calculus with Riemann–Liouville derivatives of order (α, β). Math Methods Appl Sci (2007) 30(15):1931–9. doi:10.1002/mma.879

CrossRef Full Text | Google Scholar

34. Cassell JS, Williams MM. A solution of the neutron diffusion equation for a hemisphere with mixed boundary conditions. Ann Nucl Energ (2004) 31(17):1987–2004. doi:10.1016/j.anucene.2004.04.008

CrossRef Full Text | Google Scholar

35. Khasawneh K, Dababneh S, Odibat Z. A solution of the neutron diffusion equation in hemispherical symmetry using the homotopy perturbation method. Ann Nucl Energ (2009) 36(11-12):1711–7. doi:10.1016/j.anucene.2009.09.001

CrossRef Full Text | Google Scholar

36. Sood A, Forster RA, Parsons DK. Analytical benchmark test set for criticality code verification. Prog Nucl Energ (2003) 42(1):55–106. doi:10.1016/s0149-1970(02)00098-7

CrossRef Full Text | Google Scholar

37. Yang L, Su H, Zhong C, Meng Z, Luo H, Li X, et al. Hyperspectral image classification using wavelet transform-based smooth ordering. Int J Wavelets, Multiresolution Inf Process (2019) 17(06):1950050. doi:10.1142/s0219691319500504

CrossRef Full Text | Google Scholar

38. Guariglia E. Harmonic Sierpinski gasket and applications. Entropy (2018) 20(9):714. doi:10.3390/e20090714

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zheng X, Tang YY, Zhou J. A framework of adaptive multiscale wavelet decomposition for signals on undirected graphs. IEEE Trans Signal Process (2019) 67(7):1696–711. doi:10.1109/tsp.2019.2896246

CrossRef Full Text | Google Scholar

40. Guariglia E, Silvestrov S. Fractional-Wavelet Analysis of Positive definite Distributions and Wavelets on D'(C) D′(C). In Engineering mathematics II: Algebraic, stochastic and analysis structures for networks, data classification and optimization 2016. Springer International Publishing (2017). p. 337–53.

CrossRef Full Text | Google Scholar

41. Guido RC, Pedroso F, Contreras RC, Rodrigues LC, Guariglia E, Neto JS. Introducing the Discrete Path Transform (DPT) and its applications in signal analysis, artefact removal, and spoken word recognition. Digital Signal Process. (2021) 117:103158. doi:10.1016/j.dsp.2021.103158

CrossRef Full Text | Google Scholar