Physical Investigation of on Nuclear Dynamics Equations in Nanoenergy Reactors Using a High D-O-F Variable Separation Protocol along with the Boubaker Polynomials Expansion Scheme BPES

Abstract

We present general analytical solutions to the nuclear dynamics-related Neutron Boltzmann Transport Equation inside nanoenergy reactors. Finding a particular solution to the neutron equation by making use of boundary conditions and initial conditions may be too much for the present study and reduce the generality of the solutions. Some simple assumptions have been introduced in the main system thanks to the Boubaker Polynomial Expansion Scheme, BPES, in order to make the general analytical procedure simple and adaptable for solving similar real-life problems.

1. Introduction

Several choices are possible for describing neutron behavior in a medium filled with nuclei. A neutron is a subatomic particle called a baryon having the characteristics strong force of the standard model. Thus, a quantum mechanical description seems appropriate, leading to an involved system of Schrodinger equations describing neutron motion between and within nuclei [1–7]. A neutron is also a relativistic particle with variation of its mass over time when travelling near the speed of light. Additionally, a neutron posses wave and classical particle properties simultaneously and therefore a collective description like that of Maxwell’s equations also seems appropriate. In reality, a neutron displays all of the above characteristics at one time or another. When a neutron collides with the nucleus, its strong force interacts with all of the individual nucleons [2, 4].

However, between nuclear collisions, neutrons move ballistically. Neutrons with energies above 20 Mev with speeds of more than 20% the speed of light exhibit relativistic motion, but most in a reactor are rarely above 0.17 c. The neutron wavelength is most important for ultra-low-energy neutrons mainly existing in the laboratory. Fortunately, classical neutral particle description with quantum mechanics describing collisions energies is the most appropriate for the investigation of neutron motion within a nuclear reactor. Neutron transport equation, which is also called neutron Boltzmann equation, is an equation which characterizes a relatively small number of neutrons colliding in a vast sea of nuclei [2, 4].

Mathematically, a neutron is a neutral point particle, experiencing deflection from or capture by a nucleus at the center of an atom. If the conditions are just right, the captured neutron causes a fissile nucleus to fission, producing more neutrons [2, 4]. The statistically large number of neutrons interacting in a reactor allows for a continuum-like description through averaging resulting in the linear Boltzmann equation. Also, a statistical mechanics formulation, first attempted by Boltzmann for interacting gases, provides appropriate descriptions.

2. Neutron Transport Equation

The neutron transport equation is a linearized version of the Boltzmann equation with wide applications in physics, geophysics, and astrophysics. The neutron transport equation models the transport of neutral particles in a scattering, fission, and absorption events with no self-interactions [4, 8]. It is used in radiation shielding and nanoenergy reactor calculations, as well as in radiative transfer of stellar and planetary atmospheres, and it also describes dispersion of light, the passage of γ-rays through dispersive media, and so forth [4, 8, 9].

The resolution of problems dealing with transport phenomena is the subject of several works, especially in the context of transfer multidimensional problems based on analytical and numerical approaches. One can refer, for example, to Fourier transform [10–12] and many others using the Laplace transform [13].

Chebyshev spectral methods for radiative transfer problems are also studied, for example, by Kim and Ishimaru in [14], by Kim and Moscoso in [15], and by Kadem in [16]. For more detailed study on Chebyshev spectral method and also approximations by the spectral methods, we refer the reader to monographs by Boyd [17] and Bernardi and Maday [18].

The neutron transport problem has been studied analytically [19] and recently, analytical solutions were sought using spectral methods [20]. However, most of these studies deal with the reduced neutron transport equation, in which one or two of the variables are left out (especially, the time-dependent part). Hence, we shall present an analytical procedure of solving the general neutron transport equation, while keeping all the fundamental variables.

3. Part I

3.1. Analytical Procedure

In this study, all particles including nuclei are in motion along with random collisions [3, 20]. Although several forms of the neutron transport equation exist, the integrodifferential formulation, arguably the most popular form in neutron transport and reactor physics applications, is presented and would be solved analytically for general applications. This equation is given as follows [2, 3, 20, 21]:𝑣1𝜕𝜕=𝑡+Ω⋅∇+Σ(𝑟,𝐸,𝑡)𝜓(𝑟,Ω,𝐸,𝑡)_∞0𝑑𝐸^_4𝜋+𝜒(𝐸)𝑑Ω^Σ_𝑠𝑟,Ω^⋅Ω,𝐸^⟶𝐸𝜓𝑟,Ω^,𝐸^,𝑡4𝜋^∞0𝑑𝐸^_4𝜋𝑑Ω^𝑉_𝐸𝐸^Σ_𝑓𝑟,𝐸^,𝑡𝜓𝑟,Ω^,𝐸^,𝑡+𝑄(𝑟,Ω,𝐸,𝑡),(3.1) where 𝑣 is the neutron speed, Σ and Σ_𝑓 are macroscopic cross-sections, Σ_𝑠 is the scattering cross-section, 𝜒(𝐸) is the distribution function, 𝜓 is the neutron angular flux, 𝐸 and 𝐸^ are energies, and Ω and Ω^ are neutron directions, while 𝑄(𝑟,Ω,𝐸,𝑡) is the source function [3].

Employing the method of variable separation (which is somehow similar to the multigroup theory), we shall make the following assumptions:𝜓𝜓𝑟,(𝑟Ω,_Ω,𝐸,𝑡)=𝜓₁(𝑟,𝑡)𝜓₂(Ω,𝐸),(3.2)Σ(𝑟,,𝐸𝐸,𝑡^,)𝑡==Σ₁𝜓₁(𝑟,𝑡)𝜓^2Ω^,𝐸^,(3.3)Σ(𝑟)Σ₂(𝐸)Σ₃(𝑡),(3.4)_𝑓Σ𝑟,𝐸^,𝑡=Σ_𝑓1(𝑟)Σ_𝑓2𝐸^Σ_𝑓3(𝑡),(3.5)_𝑠𝑟,Ω^⋅Ω,𝐸^⟶𝐸=Σ_𝑠1(𝑟)Σ_𝑠2Ω^⋅Ω,𝐸^⟶𝐸.(3.6) Introducing these expressions in (3.1), we write𝜓₂𝑣𝜕𝜕𝜓𝑡¹+Ω𝜓²⋅∇𝜓₁+Σ₁Σ₂Σ₃𝜓₁𝜓₂=𝜓₁Σ_𝑠1(𝑟)^∞0𝑑𝐸^_4𝜋𝑑Ω^Σ_𝑠2Ω^⋅Ω,𝐸^⟶𝐸𝜓^2+𝜒(𝐸)4𝜋⋅𝜓¹𝜓^2^∞0𝑑𝐸^_4𝜋𝑑Ω^𝑉_𝐸𝐸^Σ_𝑓1Σ_𝑓2Σ_𝑓3=𝑄(𝑟,Ω,𝐸,𝑡).(3.7a)𝑉_𝐸 represents the average number of neutrons per fission. Now, since 𝑟 is actually a vector that represents 𝑥,𝑦,𝑧 (or general) coordinates, we may write𝜓₂𝑣𝜕𝜕𝜓𝑡¹+Ω𝜓²⋅𝜕𝜓¹𝜕𝑟+Σ¹Σ₂Σ₃𝜓₁𝜓₂−𝜓₁𝜓^2Σ_𝑠1(𝑟)^∞0𝑑𝐸^_4𝜋𝑑Ω^Σ_𝑠2−𝜒(𝐸)4𝜋⋅𝜓¹𝜓^2Σ_𝑓1(𝑟)Σ_𝑓3(𝑡)^∞0𝑉_𝐸𝐸^Σ_𝑓2𝐸^𝑑𝐸^_4𝜋𝑑Ω^=𝑄(𝑟,Ω,𝐸,𝑡).(3.7b) If we also writeΣ_𝑠2Ω^⋅Ω,𝐸^⟶𝐸=Σ_𝑠𝑎Ω^⋅ΩΣ_𝑠𝑏𝐸^⟶𝐸.(3.7c) It follows that𝜓₂𝑣𝜕𝜕𝜓𝑡¹+Ω𝜓²𝜕𝜓₁𝜕𝑟+Σ¹Σ₂Σ₃𝜓₂−𝜓^2Σ_𝑠1(𝑟)^∞0𝑑𝐸′_4𝜋𝑑Ω^Σ_𝑠2−𝜒(𝐸)4𝜋⋅𝜓^2Σ_𝑓1(𝑟)Σ_𝑓3(𝑡)^∞0𝑉_𝐸𝐸^Σ_𝑓2𝐸^𝑑𝐸^_4𝜋𝑑Ω^𝜓¹=𝑄.(3.8) As far as we are concerned, the functions in 𝐸,𝐸^,Ω,Ω^ would be taken as constant variables in (3.8). We may simplify the above problem further if we write𝐴(𝑟,𝑡)=Σ₁Σ₂Σ₃𝜓₂−𝜓^2Σ_𝑠1(𝑟)^∞0𝑑𝐸^_4𝜋𝑑Ω^Σ_𝑠2−𝜒(𝐸)4𝜋⋅𝜓^2Σ_𝑓1(𝑟)Σ_𝑓3(𝑡)^∞0𝑉_𝐸𝐸^Σ_𝑓2𝐸^𝑑𝐸^_4𝜋𝑑Ω^.(3.9) Equation (3.8) becomes𝜓₂𝑣𝜕𝜕𝜓𝑡¹+Ω𝜓²𝜕𝜓₁𝜕𝑟+𝐴(𝑟,𝑡)𝜓¹=𝑄(𝑟,Ω,𝐸,𝑡).(3.10) It is assumed that 𝑄 is separable such that𝑄(𝑟,Ω,𝐸,𝑡)=𝑄₁(𝑟,𝑡)𝑄₂(Ω,𝐸).(3.11) Equation (3.10) can be written in the form𝜓₂𝑄₂𝑣𝜕𝜕𝜓𝑡¹+Ω𝑄𝜓²₂𝜕𝜓₁𝜕𝑟+𝐴𝑄(𝑟₂,𝑡)𝜓₁=𝑄₁(𝑟,𝑡).(3.12) If the function 𝜓₁ is separable such that𝜓₁(𝑟,𝑡)=𝑅(𝑟)𝐺(𝑡)(3.13) and also, if the possibility of expressing the other functions, as follows, exists, the problems become very easy to handle:𝐴(𝑟𝑄,𝑡)=𝐴₁(𝑟)+𝐴₂(𝑡),₁𝑄₀(𝑟,𝑡)=𝜓₁𝑄₀,=𝑄₀(𝑟)+𝐺₀(𝑡).(3.14) Hence,𝜓₂𝑄₂𝑣𝑅𝑑𝑑𝐺𝑡+Ω𝑄𝜓²₂𝐺𝑑𝑅𝑑𝑟+𝐴¹(𝑟)+𝐴₂(𝑡)𝑄₂𝑅𝐺=𝑅𝐺𝑅⁰𝜓₂+𝐺₀,𝑄₂𝑣𝐺1𝑑𝑑𝐺𝑡+Ω𝑄𝜓²₂1𝑅𝑑𝑑𝑅𝑟+𝐴¹(𝑟)𝑄₂+𝐴²(𝑡)𝑄₂𝜓=𝑅₀+𝐺₀,(3.15)₂𝑄₂𝑣𝐺1𝑑𝑑𝐺𝑡+𝐴²(𝑡)𝑄₂−𝐺₀(𝑡)=−Ω𝜓²𝑄₂1𝑅𝑑𝑑𝑅𝑟−𝐴¹(𝑟)𝑄₂+𝑅₀(𝑟).(3.16) Since both sides of (3.16) are independent of one another, they must be equal to a constant 𝜀², leading to the following equations:𝑑𝐺𝑑𝑡=𝑄²𝑣𝜓₂𝐺⁰(𝑡)−𝐴²(𝑡)𝑄₂+𝜀²𝑑𝑅𝐺(𝑡),𝑑𝑟=Ω𝑄𝜓²₂𝑅⁰(𝑟)+𝜀²−𝐴¹(𝑟)𝑄₂𝑅(𝑟).(3.17) Using the method of integrating factor or by direct integration, the solutions to (3.17) are as follows:𝑅𝐺((𝑟𝑡))==𝑃𝑃₀₁𝑒^{(𝑄2𝑣/𝜓2)∫{𝐺0(𝑡)−𝐴2(𝑡)/𝑄2+𝜀2}𝑑𝑡},𝑒^{(𝑄2/Ω𝜓2)∫{𝑅0(𝑟)+𝜀2−𝐴1(𝑟)/𝑄2}𝑑𝑟},(3.18) where 𝑃₀ and 𝑃₁ are constants.

Therefore, from (3.13), we could write𝜓(𝑟,𝑡)=𝑃₀𝑃₁𝑒^{(𝑄2/𝜓2∫)[𝑣{𝐺0(𝑡)−𝐴2(𝑡)/𝑄2+𝜀2∫}𝑑𝑡+(1/Ω){𝑅0(𝑟)+𝜀2−𝐴1(𝑟)/𝑄2}𝑑𝑟]}.(3.19) From (3.2), the expression for the flux is given as 𝜓(𝑟,Ω,𝐸,𝑡)=𝑃₀𝑃₁𝜓₂(Ω,𝐸)𝑒^{(𝑄2/𝜓2∫)[𝑣{𝐺0(𝑡)−𝐴2(𝑡)/𝑄2+𝜀2∫}𝑑𝑡+(1/Ω){𝑅0(𝑟)+𝜀2−𝐴1(𝑟)/𝑄2}𝑑𝑟]}.(3.20) This solution is very interesting because it leaves us with the choice of different expression for the functions 𝜙₂, 𝑄₂, 𝐺₀, 𝐴₂, 𝐴₁, and 𝑅₀. Hence, the nature of a nuclear reactor, the intended use, and different applications may be easily imposed on the analytical expression for the neutron flux.

However, there is a change that must be overcome, that is, separating the function 𝐴(𝑟,𝑡) into a sum of 𝐴₁ and 𝐴₂. Once the possibility of obtaining such separation exists, then our solution is very easy to apply in any real life situations (nuclear reactor design, transport in porous media, fractured/random media analyses, etc.). However, we shall attempt a way around the challenge.

The macroscopic cross-section is given usually as follows:Σ_𝑖𝑗(𝑟,𝐸,𝑡)=𝑁_𝑗(𝑟,𝑡)𝜎_𝑖𝑗(𝐸),(3.21) for a given nuclide 𝑗 and reaction type 𝑖, where 𝑁_𝑗(𝑟,𝑡) is the nuclear atomic density and 𝜎_𝑖𝑗(𝐸) is the microscopic cross-section.

Comparing (3.21) and (3.4), we see that the function Σ₂(𝐸) is actually the microscopic cross-section. Hence, if we assume that nuclear atomic density is independent of time:Σ_𝑖𝑗=𝑁_𝑗(𝑟)𝜎_𝑖𝑗(𝐸).(3.22) ThenΣ(𝑟,𝐸)Σ=₃(Σ𝑡)₁=1,(𝑟)Σ₂(𝐸).(3.23) This is true for the second macroscopic cross-section of the initial energy, that is, Σ_𝑓(𝑟,𝐸^,𝑡). Therefore, ifΣ_𝑓3 is taken to be 1, it follows thatΣ_𝑓𝑟,𝐸^=Σ_𝑓1(𝑟)Σ_𝑓2𝐸^.(3.24) Hence, (3.9) becomes𝐴(𝑟,𝑡)=Σ₁(𝑟)Σ₂(𝐸)𝜓₂−𝜓^2Σ_𝑠1(𝑟)^∞0Σ_𝑠𝑏𝐸^⟶𝐸𝑑𝐸^_4𝜋Σ_𝑠𝑎Ω^⋅Ω𝑑Ω^−𝜒(𝐸)4𝜋𝜓^2Σ_𝑓1(𝑟)^∞0𝑉_𝐸𝐸^Σ_𝑓2𝐸^𝑑𝐸^_4𝜋𝑑Ω^.(3.25) It would be observed that (3.25) is independent of time and this means that 𝐴₂(𝑡)=0𝐴₁(𝑟)=𝐴(𝑟)=Σ₁(𝑟)Σ₂(𝐸)𝜓₂−𝜓^2Σ_𝑠1(𝑟)^∞0Σ_𝑠𝑏𝐸^⟶𝐸𝑑𝐸^_4𝜋Σ_𝑠𝑎Ω^⋅Ω𝑑Ω^−𝜒(𝐸)4𝜋𝜓^2Σ_𝑓1(𝑟)^∞0𝑉_𝐸𝐸^Σ_𝑓2𝐸^𝑑𝐸^_4𝜋𝑑Ω^.(3.26) Therefore, the expression for the angular flux becomes𝜓(𝑟,Ω,𝐸,𝑡)=𝑃₀𝑃₁𝜓₂(Ω,𝐸)𝑒^{(𝑄2/𝜓2∫)[𝑣{𝐺0(𝑡)+𝜀2∫}𝑑𝑡+(1/Ω){𝑅0(𝑟)+𝜀2−𝐴1(𝑟)/𝑄2}𝑑𝑟]}.(3.27)

3.2. Expression for the Distribution Function 𝜒(𝐸)

From (3.26), we may write [2, 3, 20, 21]𝜒(𝐸)4𝜋𝜓^2Σ_𝑓1(𝑟)^∞0𝑉_𝐸𝐸^Σ_𝑓2𝐸^𝑑𝐸^_4𝜋𝑑Ω^=Σ₁(𝑟)Σ₂(𝐸)𝜓₂−𝐴(𝑟)−𝜓^2Σ_𝑠1(𝑟)^∞0Σ_𝑠𝑏𝐸^⟶𝐸𝑑𝐸^_4𝜋=4𝜋Σ₁Σ_𝑠𝑎Ω^⋅Ω𝑑Ω^,𝜒(𝐸)(𝑟)Σ₂(𝐸)𝜓₂−4𝜋𝐴(𝑟)−4𝜋𝜓^2Σ_𝑠1(𝑟)∫^∞0Σ_𝑠𝑏𝐸^⟶𝐸𝑑𝐸^∫_4𝜋Σ_𝑠𝑎Ω^⋅Ω𝑑Ω^𝜓^2Σ_𝑓1(𝑟)∫^∞0𝑉_𝐸𝐸^Σ_𝑓2𝐸^𝑑𝐸^∫_4𝜋𝑑Ω^.(3.28)

4. Part II

4.1. Criticality and Analytical Solutions

The neutron transport equation without delayed neutrons is given as [2, 3, 20, 21]𝑣1𝜕𝜕=𝑡+𝑄Ω_ext∇+𝜎𝑟,⃗𝐸𝜓𝑟,⃗Ω,𝐸,𝑡×𝑑𝐸𝑟,⃗_Ω,𝐸,𝑡+𝑑𝐸^′𝑑Ω^𝜎_𝑠𝑟,⃗Ω^⋅Ω,𝐸^→𝐸𝜓𝑟,⃗Ω,𝐸^,𝑡+𝜒(𝐸)𝑑Ω^𝑉_𝐸𝜎_𝑓𝑟,⃗𝐸^𝜓𝑟,⃗Ω^,𝐸^,𝑡,(4.1) where 𝑄_ext is the external sources of neutrons, 𝜎 is the microscopic cross-section, 𝑣 is the neutron speed, and 𝑉_𝐸 is the average number of neutrons per fission.

The equation assumes that all neutrons are emitted instantaneously at the time of fission. In fact, small fraction of neutrons is emitted later due to certain fission products [2, 3].

Now, if we seek as asymptotic solutions to (4.1) in the following form:𝜓𝑟,⃗Ω,𝐸,𝑡=𝜓_𝑎𝑟,⃗Ω,𝐸𝑒^𝛼𝑡,(4.2) where the solution satisfies the boundary conditions, and if the system is source free (𝑄_ext=0), (4.1) becomes𝛼𝑣+Ω⋅∇+𝜎𝑟,⃗𝐸𝜓^𝑎𝑟,⃗Ω,𝐸=+𝜒𝑑𝐸(𝐸^′)𝑑Ω𝑑𝐸^𝜎_𝑠𝑟,⃗Ω^⋅Ω,𝐸^⟶𝐸𝜓_𝑎𝑟,⃗Ω,𝐸^𝑑Ω^𝑉_𝐸𝜎_𝑓𝑟,⃗𝐸^𝜓_𝑎𝑟,⃗Ω^,𝐸^.(4.3) We shall make the following assumptions:𝜎𝜎𝑟,⃗𝑟𝐸,⃗_𝐸=𝜎₁𝑟⃗𝜎₂(𝐸),𝜓_𝑎=𝜎_𝑓1𝑟⃗𝜎_𝑓2𝐸^,𝜓𝑟,⃗Ω,𝐸=𝜓_𝑎1𝑟⃗𝜓_𝑎2(Ω)𝜓_𝑎3(𝐸),_𝑎𝜓_𝑎𝑟,⃗Ω,𝐸^=𝜓_𝑎1𝑟⃗𝜓_𝑎2(Ω)𝜓^𝑎3𝐸^,𝜎_𝑠𝑟,⃗Ω^,𝐸^=𝜓_𝑎1𝑟⃗𝜓_𝑎2Ω^𝜓^𝑎3𝐸^,𝑟,⃗Ω^⋅Ω,𝐸^⟶𝐸=𝜎_𝑠1𝑟⃗𝜎_𝑠2Ω^⋅Ω𝜎_𝑠3𝐸^⟶𝐸.(4.4) Equation (4.3) then becomes𝛼𝑣𝜓^𝑎1𝜓_𝑎2𝜓_𝑎3+Ω⋅𝜓_𝑎3𝜓_𝑎2𝑑𝜓_𝑎1𝑑𝑟+𝜎¹𝜎₂𝜓_𝑎1𝜓_𝑎2𝜓_𝑎3=𝜎_𝑠1𝜎^𝑠2𝜓^𝑎3𝑑𝐸^𝜎^𝑠3𝜓_𝑎2𝑑Ω^+𝜒(𝐸)𝑉_𝐸𝜎_𝑓2𝜓^𝑎3𝑑𝐸^𝜎^𝑓1Ω𝑑𝜓_𝑎1𝜓^𝑎2𝑑Ω^𝜓_𝑎1,(4.5)𝑑𝑟=𝜓^𝑎1𝜎^𝑠1𝑟⃗𝜓_𝑎3(𝐸)𝜎^𝑠2𝐸^⟶𝐸𝜓^𝑎3𝐸^𝑑𝐸^𝜎^𝑠3Ω^⋅Ω𝜓_𝑎2𝑑Ω^−𝛼𝑣−𝜎¹𝑟⃗𝜎₂+𝜒(𝐸)𝑉_𝐸(𝐸)𝜓_𝑎2(Ω)𝜓_𝑎3(𝐸)𝜎^𝑓2𝐸^𝜓^𝑎3𝐸^𝑑𝐸′𝜎_𝑓1𝑟⃗𝜓^𝑎2Ω^𝑑Ω^.(4.6) In the differential equation (4.6), all variables other than functions of 𝑟 are taken to be constant. Hence, we may write𝐵𝑟⃗=𝜎^𝑠1𝑟⃗𝜓_𝑎3(𝐸)𝜎^𝑠2𝐸^⟶𝐸𝜓^𝑎3𝐸^𝑑𝐸^𝜎^𝑠3Ω^⋅Ω𝜓_𝑎2𝑑Ω^−𝛼𝑣−𝜎¹𝑟⃗𝜎₂+𝜒(𝐸)𝑉_𝐸(𝐸)𝜓_𝑎2(Ω)𝜓_𝑎3(𝐸)𝜎^𝑓2𝐸^𝜓^𝑎3𝐸^𝑑𝐸^𝜎^𝑓1𝑟⃗𝜓^𝑎2Ω^𝑑Ω^Ω𝑑𝜓_𝑎1,(4.7)𝑑𝑟=𝐵𝑟⃗𝜓^𝑎1.(4.8) If we suppose that the integral ∫^∞0(𝐵(𝑟)⃗/Ω)𝑑𝑟⃗ is convergent, and taking into account the characteristics of a given nuclear reactor with spherical symmetry (𝐵(𝑟)⃗=‖𝐵(𝑟)⃗‖(𝑟/⃗‖𝑟‖⃗)=‖𝐵(𝑟)⃗‖𝑢_𝑟⃗), within radial range [0,𝑅],𝑑‖‖𝐵𝐵(𝑟()⃗𝑟||)⃗_⃗0⃗𝑟=‖‖=𝑘₁,𝑑𝑟||||_⃗0⃗𝑟=𝐵(𝑟)⃗||_{⃗𝑟=𝑅⃗𝑢𝑟}=0,𝑑‖‖𝐵(𝑟)⃗‖‖=0,𝑑𝑟||||_{⃗𝑟=𝑅⃗𝑢𝑟}=𝑘₂,(4.9) where 𝑘₁ and 𝑘₂ are core reactor characteristic constants.

For solving (4.8), the Boubaker Polynomials Expansion Scheme, BPES, [22–35] is proposed. This scheme is applied through setting the following expression:‖‖𝐵𝑟⃗‖‖=12𝑁_0𝑁⁰_𝑘=1𝜆_𝑘×𝐵_4𝑘𝑟𝑅𝜇^𝑘,(4.10) where 𝐵_4𝑘 are the 4𝑘-order Boubaker polynomials, 𝑟 is the radius (𝑟∈[0,𝑅]), 𝜇_𝑘 are 𝐵_4𝑘 minimal positive roots, 𝑁₀ is a prefixed integer, and 𝜆_𝑘|_{𝑘=1⋯𝑁0} are unknown pondering real coefficients.

The main advantage of this step lies in (4.10), which ensures verifying the four boundary conditions in (4.9), at the earliest stage of resolution protocol. In fact, due to the properties of the Boubaker polynomials [25–33], and since 𝜇_𝑘|_{𝑘=1⋯𝑁0} are roots of 𝐵_4𝑘|_{𝑘=1⋯𝑁0}, (4.9) is reduced to^𝑁0_𝑘=1𝜆_𝑘×𝐵_4𝑘𝑟𝑅𝜇^𝑘|||_⃗0⃗𝑟==^𝑁0_𝑘=1𝜆_𝑘×(−2)=2𝑘₁𝑁₀,_𝑁0_𝑘=1𝜆_𝑘×𝑑𝐵^4𝑘(𝑟/𝑅)𝜇_𝑘𝑑𝑟||||_⃗0⃗𝑟==^𝑁0_𝑘=1𝜆_𝑘×0=0,_𝑁0_𝑘=1𝜆_𝑘×𝐵_4𝑘𝑟𝑅𝜇^𝑘|||_{⃗𝑟=𝑅⃗𝑢𝑟}=^𝑁0_𝑘=1𝜆_𝑘×𝐵_4𝑘𝜇_𝑘||=0,_𝑁0_𝑘=1𝜆_𝑘×𝑑𝐵^4𝑘(𝑟/𝑅)𝜇_𝑘𝑑𝑟||||_{⃗𝑟=𝑅⃗𝑢𝑟}=^𝑁0_𝑘=1𝜆_𝑘×𝑑𝐵^4𝑘𝜇_𝑘𝑑𝑟=^𝑁0_𝑘=1𝜆_𝑘×𝐻_𝑘=2𝑘₂𝑁₀,(4.11)with𝐻_𝑘=𝑑𝐵^4𝑘(𝑟/𝑅)𝜇_𝑘𝑑𝑟||||_{⃗𝑟=𝑅⃗𝑢𝑟}=4𝜇^𝑘2−𝜇^2𝑘×∑^𝑘𝑗=1𝐵^24𝑗𝜇_𝑘𝐵_4(𝑘+1)𝜇_𝑘+4𝜇^3𝑘.(4.12) The solution is then assigned to the set of pondering real coefficients 𝜆̃_𝑘|_{𝑘=1⋯𝑁0}, which minimizes the Minimum Square functional Ψ_𝑁0:Ψ_𝑁0=⎛⎜⎜⎝^𝑁0_𝑘=1𝜆̃_𝑘×(−2)−2𝑘₁𝑁₀⎞⎟⎟⎠²+⎛⎜⎜⎝^𝑁0_𝑘=1𝜆̃_𝑘×𝐻_𝑘−2𝑘₂𝑁₀⎞⎟⎟⎠²,(4.13) which gives the following solution to (4.8):𝜓_𝑎𝑟⃗=𝑃₃𝑒^{∫((1/2𝑁0)∑𝑁0𝑘=1}𝜆̃_{𝑘×𝐵4𝑘((𝑟/𝑅)𝜇𝑘)/Ω)𝑑⃗𝑟}(4.14) with 𝑃₃ constant.

From (4.8) and our earlier assumptions, we write 𝜓_𝑎𝑟,⃗Ω,𝐸=𝑃₃𝜓_𝑎2(Ω)𝜓_𝑎3(𝐸)1/2𝑁⁰∑^𝑁0_𝑘=1𝜆̃_𝑘×𝐵_4𝑘(𝑟/𝑅)𝜇_𝑘𝜓𝑟,⃗Ω,𝐸,𝑡=𝑃₃Ω𝑑𝑟,⃗(4.15)𝜓_𝑎2(Ω)𝜓_𝑎3(𝐸)𝑒^{∫𝛼𝑡+((1/2𝑁0)∑𝑁0𝑘=1}𝜆̃_{𝑘×𝐵4𝑘((𝑟/𝑅)𝜇𝑘)/Ω)𝑑⃗𝑟}.(4.16) Equation (4.16) is the neutron flux which defined the distribution neutron in the reactor core.

4.2. The 𝐾-Eigenvalue

The 𝐾-eigenvalue from of the critically problem is formulated by assuming that 𝑉_𝐸, the average number of neutrons per fission, can be adjusted to obtain a time-independent solution such as the expression of (4.16). Hence, we replace 𝑉_𝐸 by 𝑉_𝐸/𝑘 [2, 3, 20, 21] so that such a solution is equivalent to (4.11). Therefore, it follows that𝐵𝑟⃗=𝜎^𝑠1𝑟⃗𝜓_𝑎3(𝐸)𝜎^𝑠2𝐸^′⟶𝐸𝜓^𝑎3𝐸^𝑑𝐸^𝜎^𝑠3Ω^⋅Ω𝜓_𝑎2𝑑Ω′−𝛼𝑣−𝜎¹𝑟⃗𝜎₂+𝜒(𝐸)𝑉_𝐸(𝐸)𝑘𝜓_𝑎2(Ω)𝜓_𝑎3(𝐸)𝜎^𝑓2𝐸^𝜓^𝑎3𝐸^𝑑𝐸^𝜎^𝑓1𝑟⃗𝜓^𝑎2Ω^𝑑Ω^𝑘=𝜒(𝐸)𝑉_𝐸,∫𝜎_𝑓2𝐸^𝜓^𝑎3𝐸^𝑑𝐸^∫𝜎_𝑓1𝑟⃗𝜓^𝑎2Ω^𝑑Ω^𝜓_𝑎2𝜓_𝑎3𝐵𝑟⃗−𝜎^𝑠1𝑟⃗𝜓_𝑎3(𝐸)∫𝜎^𝑠2𝐸^⟶𝐸𝜓^𝑎3𝐸^𝑑𝐸^∫𝜎_𝑠3Ω^⋅Ω𝜓_𝑎2𝑑Ω^+𝛼𝑣+𝜎¹𝑟⃗𝜎₂((𝐸4.)17). We may write 𝑘=𝜆¹𝜆₂,(4.18) where𝜆₁=𝜒(𝐸)𝑉_𝐸𝜎^𝑓2𝐸^𝜓^𝑎3𝐸^𝑑𝐸^𝜎^𝑓1𝑟⃗𝜓^𝑎2Ω^𝑑Ω^𝜆₂,=𝜓_𝑎2𝜓_𝑎3𝐵𝑟⃗−𝜎^𝑠1𝑟⃗𝜓_𝑎3(𝐸)𝜎^𝑠2𝐸^⟶𝐸𝜓^𝑎3𝐸^𝑑𝐸^𝜎^𝑠3Ω^⋅Ω𝜓_𝑎2𝑑Ω^+𝛼𝑣+𝜎¹𝑟⃗𝜎₂(𝐸()4.19.)

5. Conclusion

The analytical solution of the neutron transport equation applying the method of separation of variable has been presented. It is very interesting to note from (4.18) that the determination of the criticality of a nuclear nanoreactor depends on how term 𝜆₁ weighs compared to 𝜆₂. Three interesting cases can be illustrated as follows [2, 3, 20, 21].(1)If 𝜆₁>𝜆₂, the system is supercritical.(2)If 𝜆₁=𝜆₂, the system is critical.(3)If 𝜆₁<𝜆₂, the system is subcritical.

Although working with these analytical solutions may not be an easy task, it is very crucial to note that the advantages are quite enormous.(1)For instance, we could now see how each parameters involved in nuclear reactor modeling interact with one another and we can now understand the in-depth physics better, giving us the opportunity to manipulate the parameters within natural and experimental restrictions.(2)The solutions have a lot of undefined functions and hence, a nuclear scientist is at liberty of defining such functions as the field of application of neutron transport requires. This may be very helpful in applications related to the Boltzmann equation.(3)We now have analytical expressions that may be used to access the reliability of the computational/numerical tools, which are used in reactor modeling presently. Algorithms designed from any of these expressions are definitely going to be very fast compared to Monte Carlo analysis.(4)Boundary conditions are now easier to use because of the separation of the main functions in the neutron transport equation.

These advantages will be explored one at a time in our next investigation.

Acknowledgment

The authors wish to thank Dr. Teresa Kulikowska for making some important materials available to us in the course of this study.