# 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

Authors: O. M. Dada; O. B. Awojoyogbe; K. Boubaker; M. Agida
Journal: ISRN High Energy Physics
 (2012)

---

## 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.


---

## Body

## 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]:(3.1){1v∂∂t+Ω⋅∇+Σ(r,E,t)}ψ(r,Ω,E,t)=∫0∞dE′∫4πdΩ′Σs(r,Ω′⋅Ω,E′⟶E)ψ(r,Ω′,E′,t)+χ(E)4π∫0∞dE′∫4πdΩ′VE(E′)Σf(r,E′,t)ψ(r,Ω′,E′,t)+Q(r,Ω,E,t),
where v is the neutron speed, Σ and Σf are macroscopic cross-sections, Σs is the scattering cross-section, χ(E) is the distribution function, ψ is the neutron angular flux, E and E′ are energies, and Ω  and Ω′ are neutron directions, while Q(r,Ω,E,t) 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)ψ(r,Ω,E,t)=ψ1(r,t)ψ2(Ω,E),(3.3)ψ(r,Ω′,E′,t)=ψ1(r,t)ψ2′(Ω′,E′),(3.4)Σ(r,E,t)=Σ1(r)Σ2(E)Σ3(t),(3.5)Σf(r,E′,t)=Σf1(r)Σf2(E′)Σf3(t),(3.6)Σs(r,Ω′⋅Ω,E′⟶E)=Σs1(r)Σs2(Ω′⋅Ω,E′⟶E).
Introducing these expressions in (3.1), we write(3.7a)ψ2v∂ψ1∂t+Ωψ2⋅∇ψ1+Σ1Σ2Σ3ψ1ψ2=ψ1Σs1(r)∫0∞dE′∫4πdΩ′Σs2(Ω′⋅Ω,E′⟶E)ψ2′+χ(E)4π⋅ψ1ψ2′∫0∞dE′∫4πdΩ′VE(E′)Σf1Σf2Σf3=Q(r,Ω,E,t).VE represents the average number of neutrons per fission. Now, since r is actually a vector that represents x,y,z (or general) coordinates, we may write(3.7b)ψ2v∂ψ1∂t+Ωψ2⋅∂ψ1∂r+Σ1Σ2Σ3ψ1ψ2-ψ1ψ2′Σs1(r)∫0∞dE′∫4πdΩ′Σs2-χ(E)4π⋅ψ1ψ2′Σf1(r)Σf3(t)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′=Q(r,Ω,E,t).
If we also write(3.7c)Σs2(Ω′⋅Ω,E′⟶E)=Σsa(Ω′⋅Ω)Σsb(E′⟶E). It follows that(3.8)ψ2v∂ψ1∂t+Ωψ2∂ψ1∂r+{Σ1Σ2Σ3ψ2-ψ2′Σs1(r)∫0∞dE′∫4πdΩ′Σs2-χ(E)4π⋅ψ2′Σf1(r)Σf3(t)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′}ψ1=Q.
As far as we are concerned, the functions in E,E′,Ω,Ω′ would be taken as constant variables in (3.8). We may simplify the above problem further if we write(3.9)A(r,t)=Σ1Σ2Σ3ψ2-ψ2′Σs1(r)∫0∞dE′∫4πdΩ′Σs2-χ(E)4π⋅ψ2′Σf1(r)Σf3(t)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′.
Equation (3.8) becomes(3.10)ψ2v∂ψ1∂t+Ωψ2∂ψ1∂r+A(r,t)ψ1=Q(r,Ω,E,t).
It is assumed that Q is separable such that(3.11)Q(r,Ω,E,t)=Q1(r,t)Q2(Ω,E).
Equation (3.10) can be written in the form(3.12)ψ2Q2v∂ψ1∂t+Ωψ2Q2∂ψ1∂r+A(r,t)Q2ψ1=Q1(r,t).
If the function ψ1 is separable such that(3.13)ψ1(r,t)=R(r)G(t)
and also, if the possibility of expressing the other functions, as follows, exists, the problems become very easy to handle:(3.14)A(r,t)=A1(r)+A2(t),Q1(r,t)=ψ1Q0,Q0=Q0(r)+G0(t).
Hence,(3.15)ψ2Q2vRdGdt+Ωψ2Q2GdRdr+{A1(r)+A2(t)Q2}RG=RG(R0+G0),ψ2Q2v1GdGdt+Ωψ2Q21RdRdr+A1(r)Q2+A2(t)Q2=(R0+G0),(3.16)ψ2Q2v1GdGdt+A2(t)Q2-G0(t)=-Ωψ2Q21RdRdr-A1(r)Q2+R0(r).
Since both sides of (3.16) are independent of one another, they must be equal to a constant ε2, leading to the following equations:(3.17)dGdt=Q2vψ2{G0(t)-A2(t)Q2+ε2}G(t),dRdr=Q2Ωψ2{R0(r)+ε2-A1(r)Q2}R(r).
Using the method of integrating factor or by direct integration, the solutions to (3.17) are as follows:(3.18)G(t)=P0e(Q2v/ψ2)∫{G0(t)-A2(t)/Q2+ε2}dt,R(r)=P1e(Q2/Ωψ2)∫{R0(r)+ε2-A1(r)/Q2}dr,
where P0 and P1 are constants.Therefore, from (3.13), we could write(3.19)ψ(r,t)=P0P1e(Q2/ψ2)[v∫{G0(t)-A2(t)/Q2+ε2}dt+(1/Ω)∫{R0(r)+ε2-A1(r)/Q2}dr].
From (3.2), the expression for the flux is given as
(3.20)ψ(r,Ω,E,t)=P0P1ψ2(Ω,E)e(Q2/ψ2)[v∫{G0(t)-A2(t)/Q2+ε2}dt+(1/Ω)∫{R0(r)+ε2-A1(r)/Q2}dr].
This solution is very interesting because it leaves us with the choice of different expression for the functions ϕ2, Q2, G0, A2, A1, and R0. 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 functionA(r,t) into a sum of A1 and A2. 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)Σij(r,E,t)=Nj(r,t)σij(E),
for a given nuclide j and reaction type i, where Nj(r,t) is the nuclear atomic density and σij(E) is the microscopic cross-section.Comparing (3.21) and (3.4), we see that the function Σ2(E) is actually the microscopic cross-section. Hence, if we assume that nuclear atomic density is independent of time:(3.22)Σij=Nj(r)σij(E).
Then(3.23)Σ3(t)=1,Σ(r,E)=Σ1(r)Σ2(E).
This is true for the second macroscopic cross-section of the initial energy, that is, Σf(r,E′,t). Therefore, ifΣf3 is taken to be 1, it follows that(3.24)Σf(r,E′)=Σf1(r)Σf2(E′).
Hence, (3.9) becomes(3.25)A(r,t)=Σ1(r)Σ2(E)ψ2-ψ2′Σs1(r)∫0∞Σsb(E′⟶E)dE′∫4πΣsa(Ω′⋅Ω)dΩ′-χ(E)4πψ2′Σf1(r)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′.
It would be observed that (3.25) is independent of time and this means that A2(t)=0(3.26)A1(r)=A(r)=Σ1(r)Σ2(E)ψ2-ψ2′Σs1(r)∫0∞Σsb(E′⟶E)dE′∫4πΣsa(Ω′⋅Ω)dΩ′-χ(E)4πψ2′Σf1(r)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′.
Therefore, the expression for the angular flux becomes(3.27)ψ(r,Ω,E,t)=P0P1ψ2(Ω,E)e(Q2/ψ2)[v∫{G0(t)+ε2}dt+(1/Ω)∫{R0(r)+ε2-A1(r)/Q2}dr].


### 3.2. Expression for the Distribution Functionχ(E)

From (3.26), we may write [2, 3, 20, 21](3.28)χ(E)4πψ2′Σf1(r)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′=Σ1(r)Σ2(E)ψ2-A(r)-ψ2′Σs1(r)∫0∞Σsb(E′⟶E)dE′∫4πΣsa(Ω′⋅Ω)dΩ′,χ(E)=4πΣ1(r)Σ2(E)ψ2-4πA(r)-4πψ2′Σs1(r)∫0∞Σsb(E′⟶E)dE′∫4πΣsa(Ω′⋅Ω)dΩ′ψ2′Σf1(r)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′.


## 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]:(3.1){1v∂∂t+Ω⋅∇+Σ(r,E,t)}ψ(r,Ω,E,t)=∫0∞dE′∫4πdΩ′Σs(r,Ω′⋅Ω,E′⟶E)ψ(r,Ω′,E′,t)+χ(E)4π∫0∞dE′∫4πdΩ′VE(E′)Σf(r,E′,t)ψ(r,Ω′,E′,t)+Q(r,Ω,E,t),
where v is the neutron speed, Σ and Σf are macroscopic cross-sections, Σs is the scattering cross-section, χ(E) is the distribution function, ψ is the neutron angular flux, E and E′ are energies, and Ω  and Ω′ are neutron directions, while Q(r,Ω,E,t) 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)ψ(r,Ω,E,t)=ψ1(r,t)ψ2(Ω,E),(3.3)ψ(r,Ω′,E′,t)=ψ1(r,t)ψ2′(Ω′,E′),(3.4)Σ(r,E,t)=Σ1(r)Σ2(E)Σ3(t),(3.5)Σf(r,E′,t)=Σf1(r)Σf2(E′)Σf3(t),(3.6)Σs(r,Ω′⋅Ω,E′⟶E)=Σs1(r)Σs2(Ω′⋅Ω,E′⟶E).
Introducing these expressions in (3.1), we write(3.7a)ψ2v∂ψ1∂t+Ωψ2⋅∇ψ1+Σ1Σ2Σ3ψ1ψ2=ψ1Σs1(r)∫0∞dE′∫4πdΩ′Σs2(Ω′⋅Ω,E′⟶E)ψ2′+χ(E)4π⋅ψ1ψ2′∫0∞dE′∫4πdΩ′VE(E′)Σf1Σf2Σf3=Q(r,Ω,E,t).VE represents the average number of neutrons per fission. Now, since r is actually a vector that represents x,y,z (or general) coordinates, we may write(3.7b)ψ2v∂ψ1∂t+Ωψ2⋅∂ψ1∂r+Σ1Σ2Σ3ψ1ψ2-ψ1ψ2′Σs1(r)∫0∞dE′∫4πdΩ′Σs2-χ(E)4π⋅ψ1ψ2′Σf1(r)Σf3(t)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′=Q(r,Ω,E,t).
If we also write(3.7c)Σs2(Ω′⋅Ω,E′⟶E)=Σsa(Ω′⋅Ω)Σsb(E′⟶E). It follows that(3.8)ψ2v∂ψ1∂t+Ωψ2∂ψ1∂r+{Σ1Σ2Σ3ψ2-ψ2′Σs1(r)∫0∞dE′∫4πdΩ′Σs2-χ(E)4π⋅ψ2′Σf1(r)Σf3(t)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′}ψ1=Q.
As far as we are concerned, the functions in E,E′,Ω,Ω′ would be taken as constant variables in (3.8). We may simplify the above problem further if we write(3.9)A(r,t)=Σ1Σ2Σ3ψ2-ψ2′Σs1(r)∫0∞dE′∫4πdΩ′Σs2-χ(E)4π⋅ψ2′Σf1(r)Σf3(t)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′.
Equation (3.8) becomes(3.10)ψ2v∂ψ1∂t+Ωψ2∂ψ1∂r+A(r,t)ψ1=Q(r,Ω,E,t).
It is assumed that Q is separable such that(3.11)Q(r,Ω,E,t)=Q1(r,t)Q2(Ω,E).
Equation (3.10) can be written in the form(3.12)ψ2Q2v∂ψ1∂t+Ωψ2Q2∂ψ1∂r+A(r,t)Q2ψ1=Q1(r,t).
If the function ψ1 is separable such that(3.13)ψ1(r,t)=R(r)G(t)
and also, if the possibility of expressing the other functions, as follows, exists, the problems become very easy to handle:(3.14)A(r,t)=A1(r)+A2(t),Q1(r,t)=ψ1Q0,Q0=Q0(r)+G0(t).
Hence,(3.15)ψ2Q2vRdGdt+Ωψ2Q2GdRdr+{A1(r)+A2(t)Q2}RG=RG(R0+G0),ψ2Q2v1GdGdt+Ωψ2Q21RdRdr+A1(r)Q2+A2(t)Q2=(R0+G0),(3.16)ψ2Q2v1GdGdt+A2(t)Q2-G0(t)=-Ωψ2Q21RdRdr-A1(r)Q2+R0(r).
Since both sides of (3.16) are independent of one another, they must be equal to a constant ε2, leading to the following equations:(3.17)dGdt=Q2vψ2{G0(t)-A2(t)Q2+ε2}G(t),dRdr=Q2Ωψ2{R0(r)+ε2-A1(r)Q2}R(r).
Using the method of integrating factor or by direct integration, the solutions to (3.17) are as follows:(3.18)G(t)=P0e(Q2v/ψ2)∫{G0(t)-A2(t)/Q2+ε2}dt,R(r)=P1e(Q2/Ωψ2)∫{R0(r)+ε2-A1(r)/Q2}dr,
where P0 and P1 are constants.Therefore, from (3.13), we could write(3.19)ψ(r,t)=P0P1e(Q2/ψ2)[v∫{G0(t)-A2(t)/Q2+ε2}dt+(1/Ω)∫{R0(r)+ε2-A1(r)/Q2}dr].
From (3.2), the expression for the flux is given as
(3.20)ψ(r,Ω,E,t)=P0P1ψ2(Ω,E)e(Q2/ψ2)[v∫{G0(t)-A2(t)/Q2+ε2}dt+(1/Ω)∫{R0(r)+ε2-A1(r)/Q2}dr].
This solution is very interesting because it leaves us with the choice of different expression for the functions ϕ2, Q2, G0, A2, A1, and R0. 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 functionA(r,t) into a sum of A1 and A2. 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)Σij(r,E,t)=Nj(r,t)σij(E),
for a given nuclide j and reaction type i, where Nj(r,t) is the nuclear atomic density and σij(E) is the microscopic cross-section.Comparing (3.21) and (3.4), we see that the function Σ2(E) is actually the microscopic cross-section. Hence, if we assume that nuclear atomic density is independent of time:(3.22)Σij=Nj(r)σij(E).
Then(3.23)Σ3(t)=1,Σ(r,E)=Σ1(r)Σ2(E).
This is true for the second macroscopic cross-section of the initial energy, that is, Σf(r,E′,t). Therefore, ifΣf3 is taken to be 1, it follows that(3.24)Σf(r,E′)=Σf1(r)Σf2(E′).
Hence, (3.9) becomes(3.25)A(r,t)=Σ1(r)Σ2(E)ψ2-ψ2′Σs1(r)∫0∞Σsb(E′⟶E)dE′∫4πΣsa(Ω′⋅Ω)dΩ′-χ(E)4πψ2′Σf1(r)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′.
It would be observed that (3.25) is independent of time and this means that A2(t)=0(3.26)A1(r)=A(r)=Σ1(r)Σ2(E)ψ2-ψ2′Σs1(r)∫0∞Σsb(E′⟶E)dE′∫4πΣsa(Ω′⋅Ω)dΩ′-χ(E)4πψ2′Σf1(r)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′.
Therefore, the expression for the angular flux becomes(3.27)ψ(r,Ω,E,t)=P0P1ψ2(Ω,E)e(Q2/ψ2)[v∫{G0(t)+ε2}dt+(1/Ω)∫{R0(r)+ε2-A1(r)/Q2}dr].


## 3.2. Expression for the Distribution Functionχ(E)

From (3.26), we may write [2, 3, 20, 21](3.28)χ(E)4πψ2′Σf1(r)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′=Σ1(r)Σ2(E)ψ2-A(r)-ψ2′Σs1(r)∫0∞Σsb(E′⟶E)dE′∫4πΣsa(Ω′⋅Ω)dΩ′,χ(E)=4πΣ1(r)Σ2(E)ψ2-4πA(r)-4πψ2′Σs1(r)∫0∞Σsb(E′⟶E)dE′∫4πΣsa(Ω′⋅Ω)dΩ′ψ2′Σf1(r)∫0∞VE(E′)Σf2(E′)dE′∫4πdΩ′.


## 4. Part II

### 4.1. Criticality and Analytical Solutions

The neutron transport equation without delayed neutrons is given as [2, 3, 20, 21](4.1){1v∂∂t+Ω∇+σ(r⃗,E)}ψ(r⃗,Ω,E,t)=Qext(r⃗,Ω,E,t)+∫dE′∫dΩ′σs(r⃗,Ω′⋅Ω,E′→E)ψ(r⃗,Ω,E′,t)+χ(E)×∫dE′∫dΩ′VEσf(r⃗,E′)ψ(r⃗,Ω′,E′,t),
where Qext is the external sources of neutrons, σ is the microscopic cross-section, v is the neutron speed, and VE 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)ψ(r⃗,Ω,E,t)=ψa(r⃗,Ω,E)eαt,
where the solution satisfies the boundary conditions, and if the system is source free (Qext=0), (4.1) becomes(4.3){αv+Ω⋅∇+σ(r⃗,E)}ψa(r⃗,Ω,E)=∫dE′∫dΩ′σs(r⃗,Ω′⋅Ω,E′⟶E)ψa(r⃗,Ω,E′)+χ(E)∫dE′∫dΩ′VEσf(r⃗,E′)ψa(r⃗,Ω′,E′).
We shall make the following assumptions:(4.4)σ(r⃗,E)=σ1(r⃗)σ2(E),σ(r⃗,E′)=σf1(r⃗)σf2(E′),ψa(r⃗,Ω,E)=ψa1(r⃗)ψa2(Ω)ψa3(E),ψa(r⃗,Ω,E′)=ψa1(r⃗)ψa2(Ω)ψa3′(E′),ψa(r⃗,Ω′,E′)=ψa1(r⃗)ψa2(Ω′)ψa3′(E′),σs(r⃗,Ω′⋅Ω,E′⟶E)=σs1(r⃗)σs2(Ω′⋅Ω)σs3(E′⟶E).
Equation (4.3) then becomes(4.5)αvψa1ψa2ψa3+Ω⋅ψa3ψa2dψa1dr+σ1σ2ψa1ψa2ψa3=σs1∫σs2ψa3′dE′∫σs3ψa2dΩ′+χ(E)∫VEσf2ψa3′dE′∫σf1ψa2′dΩ′ψa1,(4.6)Ωdψa1dr=ψa1{χ(E)VEψa2(Ω)ψa3(E)σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′∫σs3(Ω′⋅Ω)ψa2dΩ′-αv-σ1(r⃗)σ2(E)+χ(E)VEψa2(Ω)ψa3(E)∫σf2(E′)ψa3′(E′)dE′∫σf1(r⃗)ψa2′(Ω′)dΩ′σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′}.
In the differential equation (4.6), all variables other than functions of r are taken to be constant. Hence, we may write(4.7)B(r⃗)=σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′∫σs3(Ω′⋅Ω)ψa2dΩ′-αv-σ1(r⃗)σ2(E)+χ(E)VEψa2(Ω)ψa3(E)∫σf2(E′)ψa3′(E′)dE′∫σf1(r⃗)ψa2′(Ω′)dΩ′,(4.8)Ωdψa1dr=B(r⃗)ψa1.
If we suppose that the integral ∫0∞(B(r⃗)/Ω)dr⃗ is convergent, and taking into account the characteristics of a given nuclear reactor with spherical symmetry (B(r⃗)=∥B(r⃗)∥(r⃗/∥r⃗∥)=∥B(r⃗)∥u⃗r), within radial range [0,R],(4.9)B(r⃗)|r⃗=0⃗=k1,d‖B(r⃗)‖dr|r⃗=0⃗=0,B(r⃗)|r⃗=Ru⃗r=0,d‖B(r⃗)‖dr|r⃗=Ru⃗r=k2,
where k1 and k2 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:(4.10)‖B(r⃗)‖=12N0∑k=1N0λk×B4k(rRμk),
where B4k are the 4k-order Boubaker polynomials, r is the radius (r∈[0,R]), μk are B4k minimal positive roots, N0 is a prefixed integer, and λk|k=1⋯N0 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 μk|k=1⋯N0 are roots of B4k|k=1⋯N0, (4.9) is reduced to(4.11)∑k=1N0λk×B4k(rRμk)|r⃗=0⃗=∑k=1N0λk×(-2)=2k1N0,∑k=1N0λk×dB4k((r/R)μk)dr|r⃗=0⃗=∑k=1N0λk×0=0,∑k=1N0λk×B4k(rRμk)|r⃗=Ru⃗r=∑k=1N0λk×B4k(μk)|=0,∑k=1N0λk×dB4k((r/R)μk)dr|r⃗=Ru⃗r=∑k=1N0λk×dB4k(μk)dr=∑k=1N0λk×Hk=2k2N0,(4.12)withHk=dB4k((r/R)μk)dr|r⃗=Ru⃗r=(4μk[2-μk2]×∑j=1kB4j2(μk)B4(k+1)(μk)+4μk3).
The solution is then assigned to the set of pondering real coefficients λ̃k|k=1⋯N0, which minimizes the Minimum Square functional ΨN0:(4.13)ΨN0=(∑k=1N0λ̃k×(-2)-2k1N0)2+(∑k=1N0λ̃k×Hk-2k2N0)2,
which gives the following solution to (4.8):(4.14)ψa(r⃗)=P3e∫((1/2N0)∑k=1N0λ̃k×B4k((r/R)μk)/Ω)dr⃗
with P3 constant.From (4.8) and our earlier assumptions, we write
(4.15)ψa(r⃗,Ω,E)=P3ψa2(Ω)ψa3(E)∫(1/2N0)∑k=1N0λ̃k×B4k((r/R)μk)Ωdr⃗,(4.16)ψ(r⃗,Ω,E,t)=P3ψa2(Ω)ψa3(E)eαt+∫((1/2N0)∑k=1N0λ̃k×B4k((r/R)μk)/Ω)dr⃗.
Equation (4.16) is the neutron flux which defined the distribution neutron in the reactor core.


### 4.2. TheK-Eigenvalue

TheK-eigenvalue from of the critically problem is formulated by assuming that VE, 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 VE by VE/k [2, 3, 20, 21] so that such a solution is equivalent to (4.11). Therefore, it follows that(4.17)B(r⃗)=σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′∫σs3(Ω′⋅Ω)ψa2dΩ′-αv-σ1(r⃗)σ2(E)+χ(E)VEkψa2(Ω)ψa3(E)∫σf2(E′)ψa3′(E′)dE′∫σf1(r⃗)ψa2′(Ω′)dΩ′,k=χ(E)VE∫σf2(E′)ψa3′(E′)dE′∫σf1(r⃗)ψa2′(Ω′)dΩ′ψa2ψa3{B(r⃗)-σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′∫σs3(Ω′⋅Ω)ψa2dΩ′+αv+σ1(r⃗)σ2(E)}.
We may write
(4.18)k=λ1λ2,
where(4.19)λ1=χ(E)VE∫σf2(E′)ψa3′(E′)dE′∫σf1(r⃗)ψa2′(Ω′)dΩ′,λ2=ψa2ψa3{B(r⃗)-σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′∫σs3(Ω′⋅Ω)ψa2dΩ′+αv+σ1(r⃗)σ2(E)}.


## 4.1. Criticality and Analytical Solutions

The neutron transport equation without delayed neutrons is given as [2, 3, 20, 21](4.1){1v∂∂t+Ω∇+σ(r⃗,E)}ψ(r⃗,Ω,E,t)=Qext(r⃗,Ω,E,t)+∫dE′∫dΩ′σs(r⃗,Ω′⋅Ω,E′→E)ψ(r⃗,Ω,E′,t)+χ(E)×∫dE′∫dΩ′VEσf(r⃗,E′)ψ(r⃗,Ω′,E′,t),
where Qext is the external sources of neutrons, σ is the microscopic cross-section, v is the neutron speed, and VE 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)ψ(r⃗,Ω,E,t)=ψa(r⃗,Ω,E)eαt,
where the solution satisfies the boundary conditions, and if the system is source free (Qext=0), (4.1) becomes(4.3){αv+Ω⋅∇+σ(r⃗,E)}ψa(r⃗,Ω,E)=∫dE′∫dΩ′σs(r⃗,Ω′⋅Ω,E′⟶E)ψa(r⃗,Ω,E′)+χ(E)∫dE′∫dΩ′VEσf(r⃗,E′)ψa(r⃗,Ω′,E′).
We shall make the following assumptions:(4.4)σ(r⃗,E)=σ1(r⃗)σ2(E),σ(r⃗,E′)=σf1(r⃗)σf2(E′),ψa(r⃗,Ω,E)=ψa1(r⃗)ψa2(Ω)ψa3(E),ψa(r⃗,Ω,E′)=ψa1(r⃗)ψa2(Ω)ψa3′(E′),ψa(r⃗,Ω′,E′)=ψa1(r⃗)ψa2(Ω′)ψa3′(E′),σs(r⃗,Ω′⋅Ω,E′⟶E)=σs1(r⃗)σs2(Ω′⋅Ω)σs3(E′⟶E).
Equation (4.3) then becomes(4.5)αvψa1ψa2ψa3+Ω⋅ψa3ψa2dψa1dr+σ1σ2ψa1ψa2ψa3=σs1∫σs2ψa3′dE′∫σs3ψa2dΩ′+χ(E)∫VEσf2ψa3′dE′∫σf1ψa2′dΩ′ψa1,(4.6)Ωdψa1dr=ψa1{χ(E)VEψa2(Ω)ψa3(E)σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′∫σs3(Ω′⋅Ω)ψa2dΩ′-αv-σ1(r⃗)σ2(E)+χ(E)VEψa2(Ω)ψa3(E)∫σf2(E′)ψa3′(E′)dE′∫σf1(r⃗)ψa2′(Ω′)dΩ′σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′}.
In the differential equation (4.6), all variables other than functions of r are taken to be constant. Hence, we may write(4.7)B(r⃗)=σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′∫σs3(Ω′⋅Ω)ψa2dΩ′-αv-σ1(r⃗)σ2(E)+χ(E)VEψa2(Ω)ψa3(E)∫σf2(E′)ψa3′(E′)dE′∫σf1(r⃗)ψa2′(Ω′)dΩ′,(4.8)Ωdψa1dr=B(r⃗)ψa1.
If we suppose that the integral ∫0∞(B(r⃗)/Ω)dr⃗ is convergent, and taking into account the characteristics of a given nuclear reactor with spherical symmetry (B(r⃗)=∥B(r⃗)∥(r⃗/∥r⃗∥)=∥B(r⃗)∥u⃗r), within radial range [0,R],(4.9)B(r⃗)|r⃗=0⃗=k1,d‖B(r⃗)‖dr|r⃗=0⃗=0,B(r⃗)|r⃗=Ru⃗r=0,d‖B(r⃗)‖dr|r⃗=Ru⃗r=k2,
where k1 and k2 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:(4.10)‖B(r⃗)‖=12N0∑k=1N0λk×B4k(rRμk),
where B4k are the 4k-order Boubaker polynomials, r is the radius (r∈[0,R]), μk are B4k minimal positive roots, N0 is a prefixed integer, and λk|k=1⋯N0 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 μk|k=1⋯N0 are roots of B4k|k=1⋯N0, (4.9) is reduced to(4.11)∑k=1N0λk×B4k(rRμk)|r⃗=0⃗=∑k=1N0λk×(-2)=2k1N0,∑k=1N0λk×dB4k((r/R)μk)dr|r⃗=0⃗=∑k=1N0λk×0=0,∑k=1N0λk×B4k(rRμk)|r⃗=Ru⃗r=∑k=1N0λk×B4k(μk)|=0,∑k=1N0λk×dB4k((r/R)μk)dr|r⃗=Ru⃗r=∑k=1N0λk×dB4k(μk)dr=∑k=1N0λk×Hk=2k2N0,(4.12)withHk=dB4k((r/R)μk)dr|r⃗=Ru⃗r=(4μk[2-μk2]×∑j=1kB4j2(μk)B4(k+1)(μk)+4μk3).
The solution is then assigned to the set of pondering real coefficients λ̃k|k=1⋯N0, which minimizes the Minimum Square functional ΨN0:(4.13)ΨN0=(∑k=1N0λ̃k×(-2)-2k1N0)2+(∑k=1N0λ̃k×Hk-2k2N0)2,
which gives the following solution to (4.8):(4.14)ψa(r⃗)=P3e∫((1/2N0)∑k=1N0λ̃k×B4k((r/R)μk)/Ω)dr⃗
with P3 constant.From (4.8) and our earlier assumptions, we write
(4.15)ψa(r⃗,Ω,E)=P3ψa2(Ω)ψa3(E)∫(1/2N0)∑k=1N0λ̃k×B4k((r/R)μk)Ωdr⃗,(4.16)ψ(r⃗,Ω,E,t)=P3ψa2(Ω)ψa3(E)eαt+∫((1/2N0)∑k=1N0λ̃k×B4k((r/R)μk)/Ω)dr⃗.
Equation (4.16) is the neutron flux which defined the distribution neutron in the reactor core.


## 4.2. TheK-Eigenvalue

TheK-eigenvalue from of the critically problem is formulated by assuming that VE, 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 VE by VE/k [2, 3, 20, 21] so that such a solution is equivalent to (4.11). Therefore, it follows that(4.17)B(r⃗)=σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′∫σs3(Ω′⋅Ω)ψa2dΩ′-αv-σ1(r⃗)σ2(E)+χ(E)VEkψa2(Ω)ψa3(E)∫σf2(E′)ψa3′(E′)dE′∫σf1(r⃗)ψa2′(Ω′)dΩ′,k=χ(E)VE∫σf2(E′)ψa3′(E′)dE′∫σf1(r⃗)ψa2′(Ω′)dΩ′ψa2ψa3{B(r⃗)-σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′∫σs3(Ω′⋅Ω)ψa2dΩ′+αv+σ1(r⃗)σ2(E)}.
We may write
(4.18)k=λ1λ2,
where(4.19)λ1=χ(E)VE∫σf2(E′)ψa3′(E′)dE′∫σf1(r⃗)ψa2′(Ω′)dΩ′,λ2=ψa2ψa3{B(r⃗)-σs1(r⃗)ψa3(E)∫σs2(E′⟶E)ψa3′(E′)dE′∫σs3(Ω′⋅Ω)ψa2dΩ′+αv+σ1(r⃗)σ2(E)}.


## 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 λ1 weighs compared to λ2. Three interesting cases can be illustrated as follows [2, 3, 20, 21].(1)
Ifλ1>λ2, the system is supercritical.(2)
Ifλ1=λ2, the system is critical.(3)
Ifλ1<λ2, 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.



---

*Source: 725318-2012-03-05.xml*