Analytical Solution for Multi-Energy Groups of Neutron Di ﬀ usion Equations by a Residual Power Series Method

: In this paper, a multi-energy groups of a neutron di ﬀ usion equations system is analytically solved by a residual power series method. The solution is generalized to consider three di ﬀ erent geometries: slab, cylinder and sphere. Di ﬀ usion of two and four energy groups of neutrons is speciﬁcally analyzed through numerical calculation at certain boundary conditions. This study revels su ﬃ cient analytical description for radial ﬂux distribution of multi-energy groups of neutron di ﬀ usion theory as well as determination of each nuclear reactor dimension in criticality case. The generated results are compatible with other di ﬀ erent methods data. The generated results are practically e ﬃ cient for neutron reactors dimension.


Introduction
The nuclear reactor is a complex fuel system competition, with reflectors, coolants, control rods and other parts. The design and analysis of such reactors for different operating methods is part of a complex task involving several disciplines of nuclear engineering.
Among the most important is the determination of neutron flux distribution within the reactor core and finding the critical dimension and mass. This issue has received considerable attention in the field of reactor physics in the past decades.
The balance between neutron production from fission chain reaction and neutron loss due to radiative capture and neutron leakage must be always achieved. This balance is known as criticality, it can be mathematically represented by the steady-state neutron transport equation which can be simplified using Fick's law [1] to neutron diffusion equations. Because neutrons in a reactor have different velocities it's more convenient to represent their fluxes by multi-energy groups of neutron diffusion equations.
The nuclear reactor theory defines the probability of taking place at certain reaction by the concept of cross section which is defined as the effective size of a reacted nucleus. The cross section in nuclear reactions could be characterized due to fission, σ f ; absorption, σ a ; scattering; or σ γ radiative capture. The macroscopic cross section is a product of atomic density of the target (N) and the cross section (σ) : Σ = Nσ. It is appointed to the ith energy group cross section as Σi in multi-energy groups of neutrons system, and for each nuclear reaction in the ith energy group its can be appointed for example Σ f2 is the second energy group macroscopic fission cross section, and especially in the scattering from the ith energy group to jth energy group the macroscopic cross section will be in the specific form as Σsij and the total scattering from the ith energy group to the all other energy groups will be Σsij. Neutrons of multi-energy groups move with different speeds inside the reactors. Our analysis is focused on the two and four energy groups of neutrons systems due to the necessary simplification that already applied at other alternative solution like classical methods and transport theory data [2].
The present study introduces sufficient analytical procedure based on residual power series method (RPSM) [3][4][5][6][7][8][9][10][11][12] to provide general solution to multi-energy groups of neutron diffusion equations in rectangular, cylindrical, and spherical geometries. RPSM is an effective technique to solve multi-energy groups of neutron diffusion equations without discretization, perturbation or linearization. RPSM constructs an approximate analytical solution in a practical truncated series form by using the residual error concept.
The RPSM has many characteristic features [6]; Taylor's expansion is obtained from the solution; as a result, an exact solution is available when the solution is an elementary function. Moreover, solutions and all derivatives are applicable to each arbitrary point in a determined interval. The RPSM also requires small computational requirements at high resolution and less time.
Neutron diffusion equations are solved at several cases and neutron diffusion in a hemisphere with mixed boundary conditions are presented [13]. Both homotopy analysis and Adomian decomposition methods are used to find the solution of one energy group of neutron diffusion equation in hemispherical and cylindrical reactors [14]. Moreover, the homotopy perturbation method is used to generate particular solution of one energy group of neutron diffusion equations in hemispherical symmetry [15]. An alternative method is used to solve one energy group of a neutron diffusion equation in cylindrical symmetry [16]. There are other efforts to solve different cases of neutron diffusion equations [17][18][19][20][21].
The present study provides general analytic approximate solution for the multi-energy groups of neutron diffusion at any particular symmetry. The proposed RPSM does not only describe neutron diffusions at different reactor geometries but it also provides associated flux distributions. The developed analytical formalism is numerically computed, the obtained result is compared to the available benchmark data generated from classical method and transport theory (when the data is available) to show the efficient of RPSM.
The following outlines the study: in the following section, the multi-energy groups of neutron diffusion is described analytically, and the outline of the RPSM is illustrated in Section 3. Section 4 covers how RPSM is applied in solving equations. Numerical examples are given in Section 5 to illustrate the capability of the proposed method. The study ends in Section 6 with some concluding observations.

Multi-Energy Groups of Neutron Diffusion System
The solution for multi-energy groups of neutron diffusion equations system of the following form is assumed to have a unique solution in the interval of integration: where C ii is known as a group buckling, D i is a group diffusion coefficient, and C ij is a constant connects between fluxes in different energy groups of neutrons. Respectively C ii , C ij , and D i are defined as following: The constants in Equation (2) have been defined in term of different macroscopic cross sections, the number of neutrons produced per fission for each group (ν i ) and the fraction of fission neutrons that emitted with energies in the ith group (χ i ).
This system of equations describes behavior of the neutrons in the nuclear reactors where each flux ∅ i expresses the neutron flux with specific speed. Each flux is maximum at the center of the reactor, its derivative vanishes, so the initial conditions can be written as where the fluxes ∅ i (r), i = 1, 2, . . . , n are functions of independent variable r, and h i ∈ R. Throughout it is assumed that ∅ i (r), i = 1, 2, . . . , n are analytic functions for r ≥ 0.
The system of equations in Equation (1) can be simplified by taking into consideration the basic nuclear reactor theory fact which issuer that the geometrical buckling B 2 should be equal to the material buckling for all energy groups in the criticality case, i.e., Substitution the values of B 2 leads to This system of equations shows that the ratio of each two fluxes is constant, classical solution, like Cramer's rule, determines the value of B 2 as well as finds each flux separately for any number of energy groups for any reactor geometry [22,23].

Basic Idea of the RPSM
The basic definition in addition to the basic theories of the RPSM and its applicability to different types of differential equations is given in [3][4][5][6]. To provide convenience to the reader, we will provide a review of the RPSM by illustrating the following algorithm: Step 1: Write the nth-order differential equation in the following form: subject to the initial conditions where F is a function of r, ∅, ∅ , ∅ , . . . , ∅ (n−1) and ∅ is unknown analytic function of r on a neighborhood of r 0 .
Step 2: Assume the solution has the following form: According to the initial conditions (7), the solution can be expressed in the form: Step 3: Define the kth-truncated series of ∅(r) as follows Step 4: Define the kth-residual functions as follows Step 5: Substitute the form of ∅ k (r) into Equation (11).
Step 6: Obtain the (k − n)th derivative of Res k (r) as: Step 7: Solve the following algebraic equations for k = n, n + 1, n + 2, . . ., equation by equation then obtain the values of the unknown coefficients; c n , c n+1 , c n+2 , . . . , c k , respectively.
Step 8: Collect the obtained forms of c m and ∅ m (r) for each m = 0, 1, 2, . . . , k in term of expanded series, then we obtain the kth-approximate solution to the Equations (6) and (7).
Step 9: If there are a pattern in the coefficients of the series as terms of some well-known elementary functions, then we have the exact solution ∅(r).

Analytical Solution by the RPSM
RPSM is employed to find sufficient series solution for the multi-energy groups of neutrons time-independent neutron diffusion system of Equation (1) for three nuclear reactor essential geometries, namely spherical, cylindrical and slab reactors which will be studied respectively.
The RPSM express the solution of Equation (14) subject to the initial conditions Equation (3) as a power series expansion about the initial point r = 0. The series solution is supposed to take the form ∅ i (r) = ∞ m=0 a i,m r m , i = 1, 2, . . . , n.
Repeat same procedure for k = 4, k = 5 and k = 6, values of a i4 , a i5 and a i6 are given in the following formulas: Define the recurrence relation: Therefore, the solution of the system (14) subject to the initial conditions (3) can be arranged in the form In a special case, if T i,2k = 1, then It should be noted that the simplification of Equation (24) coincides with a one energy group of a neutron diffusion system.
Using the same methodology of RPSM which used in treating the spherical reactor in the last subsection we reach to: Substitute Equation (15) into Equation (26), we have: To obtain the second approximate solutions, we place k = 2 in Equation (27), differentiate both sides of the equation with respect to r and substitute r = 0 in the resulting equation, we conclude d dr The algebraic equations, d dr Res 2 i (0) = 0, i = 1, 2, . . . , n, give the following values for a i2 : Thus, the second truncated series solutions (the second approximate solutions) for Equation (25) can be written as Similarly, to find the third approximate solution, put k = 3 in Equation (27) and solve the algebraic equations, d 2 dr 2 Res 2 i (0) = 0, i = 1, 2, . . . , n, that give us a i3 = 0, i = 1, 2, . . . , n. Repeat same procedure for k = 4, k = 5 and k = 6, the values of a i4 , a i5 and a i6 will become available as follows: Using of the recurrence relation: . ., will help us write the series coefficients obtained as follows: Thus, the solution of the cylindrical reactor system (25) with the conditions (3) can be expressed in the following form: In a special case, if T i,2k = 1, then where J 0 (r) is a Bessel's function of the first kind. It is noted that the solution given in of Equation (33), is the same as the solution of one energy group of a neutron diffusion system.

Multi-Energy Groups Slab Reactor
The final example is the slab reactor, the multi-energy groups of time-independent neutron diffusion system of a slab reactor can be expressed as: . . . ∅ n (r) + (C n,1 ∅ 1 (r) + C n,2 ∅ 2 (r) + C n,3 ∅ 3 (r) + · · · + C n,n ∅ n (r)) = 0. According to Equation (11), the kth-residual function of the Equation (34) will be as follows: Substitute Equation (15) into Equation (35), we have: For, k = 2, differentiate both sides of Equation (36) with respect to r and substitute r = 0, to obtain d dr The solution of the algebraic equations, d dr Res 2 i (0) = 0, i = 1, 2, . . . , n gives the following values for a i2 : Thus, the second truncated series solutions for Equation (35) can be written as Like previous cases, the third approximate solution of the system in Equation (38) is the same as the second approximate solution since the fourth coefficient of the series (15), a i3 , is zero.
Repeat same procedure for k = 4, k = 5, and k = 6, values of a i4 , a i5 , and a i6 are given in the following formulas: Again, by the recurrence relation: , m = 2, 4, 6, . . ., we can arrange the coefficients of the series as follows: Therefore, the solution of the slab reactor system (34) with the conditions (3) can be expressed in the following form: In a special case, if T i,2k = 1, then ∅ i (r) = cos r which coincides with one energy group of a neutron diffusion system as well.

Practical Numerical Results
Generating numerical results requires specifying boundary conditions. The flux is assumed to vanish at the surface of the reactor according zero flux boundary condition. Indeed, extrapolated boundary condition (EBC) assumes the flux vanishes at small distance beyond the surface. Couple boundary conditions are considered in such a way diffusion flux vanishes at the critical dimension. As in the same procedure that used in previous works like [20,21], one of the fluxes is normalized (to one) at the core of the reactor, and the other fluxes will be calculated depending on Equation (5) as ratios of the normalized flux, for two energy groups of neutrons the ratio between thermal and fast fluxes is given in the data source [2] which convenient with Equation (5) calculations.
Following numerical examples results are compare with the classical solutions, some of them are also compared with transport theory data in case transport theory data is available.
It should be mentioned that, as any series solution, increasing the order of approximation increases the accuracy of the solution. Therefore, to obtain accurate numerical solutions, all numerical values and graphs are calculated for the 10th approximation of the power series solution obtained in Equations (23), (32), and (41) and it is compared with the 10th approximation of the series solution obtained by the classical method.

Two Energy Groups of Neutrons Numerical Example
The two energy groups of a neutron diffusion equations system is discussed in this section at three cases; spherical reactor, infinite cylindrical reactor and slab reactor such diffusion system satisfies the following form: subject to the initial conditions Equation (42) consider only couple speeds of neutrons: fast and thermal, this simplification is previous considered at existing transport data [2].
So, the fast and thermal fluxes for spherical reactor according to RPSM formula are given by and the fluxes of the cylindrical reactor are given by whereas for slab reactor is The derived analytical formalism is computationally determined not only to verify the theory but also to compare numerical results with a classical solution [22][23][24][25]. For this purpose, the Mathematics 9 software package is implemented to generate numerical solutions. The solution is obtained numerically for cross sections related to interactions of fast and thermal neutron diffusing in 93% enriched Uranium for geometries slab, sphere and cylinder geometries. Table 1 obtained from [2] is correspondingly used. Fast Energy Group Thermal Energy Group According this data, values of C ij , i, j = 1, 2 are determined as shown in Table 2. Table 2. The values of the coefficients C i j , i, j = 1, 2 are calculated from Equation (3).
Following subsections show numerical results for each reactor geometry.

Spherical Reactor
The spherical reactor is studied by finding its critical radius of a 93% enriched Uranium and distribution of both fast and thermal fluxes as well as their sum, the total flux.
To calculate the spherical reactor critical radius a c , Equation (44) is used, at ZF and EBC boundary conditions, respectively. The generated data are listed in Table 3, other numerical data that generated by classical methods and transport theory is illustrated for comparison. This reactor critical radius (using RPSM) reproduced that of classical Method and it is in good approximation with it and with our benchmark (Transport Theory).
The fluxes distribution in the spherical reactor is reported in Table 4 and plotted in Figure 1. Table 4 provides the normalized thermal, fast and total flux values across the system for different values of r/a c .   Table 4 and Figure 1 gives the values of the sub fluxes and their total, it is obvious that all fluxes converge at the same point.
The sufficient of RPSM solution for two energy groups of neutrons system is evident, the method can definitely reproduce canonical results for such geometry.

Infinite Cylindrical Reactor
RPSM is used in this section to calculate the reactor critical radius of 93% enriched Uranium cylindrical reactor, Equation (45) is implemented with ZF and EBC boundary conditions, the results are listed in Table 5.   Table 4 and Figure 1 gives the values of the sub fluxes and their total, it is obvious that all fluxes converge at the same point.
The sufficient of RPSM solution for two energy groups of neutrons system is evident, the method can definitely reproduce canonical results for such geometry.

Infinite Cylindrical Reactor
RPSM is used in this section to calculate the reactor critical radius a c of 93% enriched Uranium cylindrical reactor, Equation (45) is implemented with ZF and EBC boundary conditions, the results are listed in Table 5. The flux distribution of such cylindrical reactor is reported in Table 5 at which fluxes values across the system for different values of r/a c are listed. The RPSM values are compared with classical diffusion calculations [22][23][24][25]. Transport results for both nuclear reactor critical radius and fluxes distribution are unfortunately not available to compare with. Figure 2 illustrates graphically the fluxes distribution in this fissile system, while the numerical values of them is tabulated in Table 6, the maximum value of the flux is at the axis of the cylinder (r = 0), it decreases towards the surface as implied by the symmetry of the system. Applying RPSM is also sufficient for two energy groups of a neutron diffusion system at cylindrical symmetry as shown in Figure 2 and Table 6.

Slab Reactor
The last two energy groups of neutrons computational analysis in this study is an infinite slab geometry. To calculate the nuclear reactor critical dimension of a 93% enriched Uranium slab, Equation (46) computationally implemented at ZF and EBC boundary conditions, the exact nuclear reactor critical dimension results of both RPSM and Classical Method compared with transport theory data is illustrated in Table 7.   Applying RPSM is also sufficient for two energy groups of a neutron diffusion system at cylindrical symmetry as shown in Figure 2 and Table 6.

Slab Reactor
The last two energy groups of neutrons computational analysis in this study is an infinite slab geometry. To calculate the nuclear reactor critical dimension a c of a 93% enriched Uranium slab, Equation (46) computationally implemented at ZF and EBC boundary conditions, the exact nuclear reactor critical dimension results of both RPSM and Classical Method compared with transport theory data is illustrated in Table 7. The fluxes distribution in the slab reactor is reported in Table 8 at which fluxes values across the system for different values of r/a c are listed, the transport theory fluxes data is available to compare with, the good agreement with it can be considered as one step forward.  Figure 3 illustrates the flux distribution of a slab fissile material, the RPSM values are compared to those obtained from classical diffusion calculations as well as from transport theory. The fluxes distribution in the slab reactor is reported in Table 8 at which fluxes values across the system for different values of / are listed, the transport theory fluxes data is available to compare with, the good agreement with it can be considered as one step forward.   Figure 3 illustrates the flux distribution of a slab fissile material, the RPSM values are compared to those obtained from classical diffusion calculations as well as from transport theory.

Multi-Energy Groups of Neutrons Numerical Example
Four energy groups of a neutron diffusion equations case is discussed as one step forward, it is represented by the following system: subject to the initial conditions:

Multi-Energy Groups of Neutrons Numerical Example
Four energy groups of a neutron diffusion equations case is discussed as one step forward, it is represented by the following system: subject to the initial conditions: The spherical reactor, infinite cylindrical reactor and slab reactor are considered respectively, in addition to assuming four different speeds of neutrons.
So, the solution of spherical reactor for four energy groups of neutrons according to RPSM formula are given by and the fluxes of the cylindrical reactor are given by where that for slab reactor is The derived analytical formalism is computationally determined not only to verify the theory, but also to compare numerical results classical solution [22][23][24][25], Mathematics 9 software package is implemented to generate numerical solutions. The solution is obtained numerically for cross sections related to interactions of four energy groups of neutrons for slab, sphere and cylinder geometries. The following data, Table 9, obtained from [25], is used correspondingly.
Unfortunately, the transport theory data is not available to compare with for all critical reactors dimensions nor for their fluxes values.
The following subsections illustrate numerical results for three kinds of reactors geometry.

Spherical Reactor
The spherical reactor is studied by finding its critical radius of the reactors and distribution of four fluxes as well as their sum.
The critical radius a c for spherical reactor is calculated at ZF and EBC boundary conditions, respectively. The generated data are listed in Table 11. The values of the fluxes in the reactor is given in Table 12 and Figure 4. Table 12 provides four fluxes and their total flux values across the system, it is clear (in both Table 12 and Figure 4) that all fluxes decreases when the reactor radius increases and vanishes at the critical radius. Unfortunately, the transport theory data is not available to compare with for all critical reactors dimensions nor for their fluxes values.
The following subsections illustrate numerical results for three kinds of reactors geometry.

Spherical Reactor
The spherical reactor is studied by finding its critical radius of the reactors and distribution of four fluxes as well as their sum.
The critical radius a for spherical reactor is calculated at ZF and EBC boundary conditions, respectively. The generated data are listed in Table 11. The values of the fluxes in the reactor is given in Table 12 and Figure 4. Table 12 provides four fluxes and their total flux values across the system, it is clear (in both Table 12 and Figure 4) that all fluxes decreases when the reactor radius increases and vanishes at the critical radius. RPSM results (as well as for classical method) for multi-energy groups of neutrons system in spherical reactor assure and improve the calculations of that of two energy groups of neutrons. RPSM results (as well as for classical method) for multi-energy groups of neutrons system in spherical reactor assure and improve the calculations of that of two energy groups of neutrons.

Infinite Cylindrical Reactor
Now, RPSM is used in this section to calculate the infinite cylindrical reactor critical radius a c , Equation (50) is implemented with ZF and EBC boundary conditions, the results are listed in Table 13. The critical radius calculations using RPSM is confirmed using classical method calculation (there is no transport theory data to compare with as in the case of two energy groups).
The flux distribution of such cylindrical reactor is reported in Table 14 at which fluxes values across the system for different values of r/a c are listed. The RPSM values are compared with classical diffusion calculations [22][23][24][25]. Figure 5 illustrates fluxes distribution in this fissile system, the maximum value of the flux is at the axis of the cylinder (r = 0), it decreases towards the surface as implied by the symmetry of the system.

Slab Reactor
Computational analysis of four energy groups of neutrons at an infinite slab geometry is also performed. To calculate the critical dimension a in slab reactors, Equation (51) implemented at ZF and EBC boundary conditions, the results are listed in Table 15.

Slab Reactor
Computational analysis of four energy groups of neutrons at an infinite slab geometry is also performed. To calculate the critical dimension a c in slab reactors, Equation (51) implemented at ZF and EBC boundary conditions, the results are listed in Table 15. The four energy groups of neutron fluxes values in the slab reactor is reported in Table 16 and all fluxes are plotted in Figure 6; this figure illustrates the flux distribution of a slab fissile material. The RPSM values are compared to those obtained from classical diffusion calculations.

Error Analysis
Since the transport theory data for four energy groups of neutrons is not available for comparison, we compared the results obtained in RPSM with the results obtained in the classical method only. In fact, both RPSM and classical methods provide us with series solutions. It is noted that the general terms in both series are different, but both solutions give us the same numerical results at the same order of approximation. Therefore, we will make another comparison between the results obtained by RPSM and classical method by calculating the residual errors (Tables 17 and  18). For the sake of brevity, we will limit the comparison to one case, the spherical reactor for the four energy groups of neutron diffusion.

Error Analysis
Since the transport theory data for four energy groups of neutrons is not available for comparison, we compared the results obtained in RPSM with the results obtained in the classical method only. In fact, both RPSM and classical methods provide us with series solutions. It is noted that the general terms in both series are different, but both solutions give us the same numerical results at the same order of approximation. Therefore, we will make another comparison between the results obtained by RPSM and classical method by calculating the residual errors (Tables 17 and 18). For the sake of brevity, we will limit the comparison to one case, the spherical reactor for the four energy groups of neutron diffusion.  The residual errors for the system in Equation (47) is defined as follows: Res i (r) = r∅ i (r) + 2∅ i (r) + r(C i,1 ∅ 1 (r) + C i,2 ∅ 2 (r) + C i,3 ∅ 3 (r) + C i,4 ∅ 4 (r)) , i = 1, 2, 3, 4, (52) and the kth-residual error is defined as: Res k i (r) = r d 2 dr 2 ∅ k i + 2 d dr ∅ k i (r) + r C i,1 ∅ k 1 (r) + C i,2 ∅ k 2 (r) + C i,3 ∅ k 3 (r) + C i,4 ∅ k 4 (r) , i = 1, 2, 3, 4, (53) where ∅ k i is the kth-approximation of ∅ i , i = 1, 2, 3, 4.

Conclusions
RPSM is applied in solving multi-energy groups of a neutron diffusion equation for different reactor geometries. This method successfully analyzes both two and four energy groups of neutron diffusion systems numerically. The flux distributions of multi-energy groups of a neutron diffusion system are presented and the determination of each nuclear reactor critical radius is illustrated. The calculated data is compatible with other compared methods such as the classical method and the available transport theory data. The efficiency of RPSM appears in introducing it as a new approximate method which can reproduce nuclear reactor fluxes directly, without using additional methods, the numerical results can be reached easily too. Moreover, RPSM, like other analytical methods, provided a series solution for multiple energy groups of neutron diffusion equations and gave us the same results obtained in classical and homotopy perturbation methods at the same order of approximation because we deal with linear equations. Despite all this, the advantage of this method is the ease, and speed of the solution in RPSM. We expect that the RPSM application on nonlinear equation models will prove to be more accurate than other analytical methods.