Solution of the one-speed neutron transport equation in anisotropic scattering medium with a Pomraning-Eddington approximation approach

In this work, a combined modified K-integral approach and Pomraning-Eddington approximation to solve the neutron transport equation for general anisotropic scattering law is presented. This combined approach leads to the establishment of general expressions of the even e(x, μ) and odd o(x, μ) functions of the angular variable μ. Both functions are slowly varying in x. Furthermore, these expressions can be reduced to those reported in literature. In addition, the parameters such as D p = ∫ − 1 1 μ 2 p e ( x , μ ) d μ and d p = ∫ − 1 1 μ 2 p + 1 o ( x , μ ) d μ with p = 1, 2, … are obtained without any integral calculus and the variables we used α1, α2, r1 and r2 verify the relation α 1 r 2 r 1 − α 2 r 1 r 2 = 0 . As application, diffusion lengths are calculated for linear scattering and the neutron heat fluxes are evaluated graphically for four-parameter scattering phase function.

D e x d ,

Introduction
In modern physics, neutron transport being an exciting topic and has been at the center of many studies. Its technological applications, firstly limited to nuclear reactor span in other fields of sciences such as radiative transfer, and biomedical domain just to name a few. The study of neutron transport has a long history and concerns the understanding of the complex mechanisms involve with the migration of neutrons in bulk media. These mechanisms start immediately after a neutron is released into a diffusing medium, this neutron can collide several times with the atomic nuclei of the medium, or may also result in absorption (where some neutrons disappear) or fission (where other neutrons reappear) or simply a change in their speed and direction (scattering).
In a slab medium those complex mechanisms are well modeled by one-speed neutron transport (Boltzmann) equation (NTE). Though obtaining one speed neutron equation remains straightforward, its solution is complicated. As a result, numerical and approximation methods are often applied. A large collection of numerical methods were developed, this include: the discrete ordinate method (DOM) [1], finite element method (FEM) [2], finite volume method (FVM) [3], monte carlo method (MCM) [4], spectral element method [5], discrete spherical harmonics method [6] an so on. More informations on this can be found in [7]. These method are versatile, but susceptible from numerical error due to finite step size [8].
A closed form analytical solution of the NTE have been obtained for some limited number of simple anisotropic scattering expressions. Some solution methods and classic references are reported in [9][10][11][12] . These solutions are useful in a variety of applications, such as: conducting sensitivity analyses to investigate the effects of various parameters, serving as screening models or benchmark solutions for more complex neutron-nucleus interaction processes that cannot be solved exactly, and for validating more comprehensive numerical solutions of the NTE [13].
These analytical solutions can be obtained by using the Pomraning-Eddington approximation (PEA) method [14][15][16]. The objective of the latter is that it separates variable of the NTE into angular and spatial part. However, until now the PEA approximation has been used only with one or two terms of the scattering phase function (expanded as an infinite series in terms of the Legendre polynomials), in order to accurately represent complex interacting phenomenon between neutron and the atom of the propagating medium, more terms of that series may be needed to represent scattering phase function accurately. According to [17] the larger the number of terms, the greater the accuracy of the phase function.
In his present form the PEA approximation can hardly handle a phase function with a great number of terms, like the four-parameter phase function of [18](this is due to the fact that diffusion coefficient D e x d , ) are embedded in the angular part (e(x, μ) and o(x, μ)), henceforth their expressions lead to transcendental equation which is difficult to resolve analytically when the number of term increases. Henceforth, it is necessary to extend the Pomraning-Eddington method to any kind of anisotropic scattering [19]. reported the above mentioned solution, giving only the generalized expression of the angular part of the NTE equation. This work aims at extending that of [19] by solving the NTE with generalized anisotropic scattering by (PEA) approach, compared to [19] the novelties introduced are: • Using a modified K-integrals approach combined with PEA approximation, to derive the solution of the spatial and the angular part of the NTE equation, • The diffusion coefficients such as D e x d , ò m m m = -+ ( ) with p=1, 2, K are obtained without any integral calculus for the first time but with the well-known Chandrasekhar Polynomial h n (λ).
The K-integral methods was developed by [20]. According to [21] there are two special weighted integrals of the one-dimensional transport equation, which are valid regardless of boundary conditions. The solution can be obtained based on the following steps: (a) Insert the PEA in the NTE; (b) Obtain even and odd functions using a modified K-integrals approach; (c) Obtain two differential equations; (d) Solve the resultant system of equations.
In order to attain this objectives, the work is organized as follows: in section 2, the mathematical description of the model is given and comparison of the PEAA model with existing ones is also made. Section 3 is devoted to the application of the model for calculation of diffusion length for linear scattering, four and six terms expansion of the scattering phase function and evaluation of diffusions parameters d 1 and D 1 via Chandrasekhar first order polynomial. The overview is founded in section 4 as conclusion.

Mathematical description
The stationary one-speed neutron or particle transport equation in homogeneous anisotropic scattering medium with internal energy source can be written as [22] I x x Here x is the optical depth space variable, m¢ and μ are the direction of neutron velocity before and after scattering, respectively, I(x, μ) is the neutron density distribution, ω o is the mean number of secondary neutrons per collision, Q(x) is the internal sources of energy, f i , i d r are the externally incident flux, diffusive reflectivity, respectively on the left and right boundaries (i=1, 2). P , m m¢ ( ) is the scattering phase function considered as [19] where P n (μ) is the Legendre Polynomial and a n are the expansion coefficients with a 1 o = . Following the Pomraning-Eddington approximation the solution of the one-speed neutron transport equation is written as [14] I f ( ) and f 1 (x) are the neutron radiant energy and net neutron flux, respectively, defined as is odd in μ both functions are slowly varying in x and normalized such that:

Modified K-integrals approach
This work starts by substituting equation (4) in equation (1), to obtain The derivation of the K-integrals in the context of neutron thermalisation by Kladnik and Kuscer [20], the reader is referred to the excellent work of [21] for it's mathematical foundation. Inspired by the K-integrals method [21], it was assumed that the angular part of the second expression on the left and the fourth expression on the left are linked by the following couple system where α 1 and α 2 are constants to be determined laters. These are two couple integral equation for o(x, μ) and e(x, μ). It should be noted that the classical K-integrals approach can be recovered by setting α 1 =λ 2 (eigenvalue) and α 2 =1. Henceforth, equations (8a) and (8b) generalize the K-integrals approach. Every odd and even function of μ can be expressed as: where g 1 (x, μ) and g 2 (x, μ) will be determined later.
In order to take into account the coupling between e(x, μ) and o(x, μ), it can be considered that g 1 (x, μ) and g 2 (x, μ) are linearly proportional to the same function g(x, μ) as: , , , , where the parameters r 1 and r 2 will be determined later. Inserting equations (9) and (10) in equations (8a) and (8b), considering equations (11a) and (11b), it can be verified that g(x, μ) satisfies the integral equation: In order to eliminate g(x,−μ) in equation (12) it is assumed that r r r r 0, 13 It's important to mention that the relation equation (13) is not only a simplification consideration (in fact the system equations (8a) and (8b) is only solvable under the condition equations (11a), (11b) and equation (13), from this simple relation the eigenvalue λ can be also determined by the monoenergetic neutrons. It's also noticed that this relation does not affect the generality of the problem under consideration. Using the relation equation (13) in equation (12) the following relation is obtained where the eigenvalue λ is given by Equation (14) is a well known classic eigenvalue problem [19]. h n (λ) is defined by Multiplying equation (14) by P n (μ) using equation (22) and integrating in μ from −1 to 1, the following recurrence relation is obtained   (24) and (18) are equal, likewise equations (25) and (19). Those equalities completely determine α 1 , α 2 , r 1 , r 2 , D p and d p . Hence, the formulation of the method is complete.

Comparison with existing models
Equations (18) and (19) are the generalization of the equations (24) and (25) respectively derived by [19] in their work, that formulated the extended Eddington approximation with anisotropic scattering. For example the expression of the even and odd function derived by [19] can be deduced, respectively from equation (18) by considering r 1 =2C e and from equation (19) by considering r 2 =2C o , where C e and C o are constants of the paper [19]. Furthermore, in our approach the parameters such as D p and d p are directly obtained whereas in other approach they are not easily obtained and long integrals calculus and the parameters α 1 , α 2 , r 1 , r 2 verify the simplified mathematical relation equation (13).

Examples of applications
In order to show the advantage of the current method we will apply it to four-parameters phase function. As mentioned in the introduction in this work it, should be noted that the PEA can hardly handle this kind of phase function. Consider equation (1), where the neutron-nucleus scattering can be governed by the following binomial scattering law [24] P L cos 1 2 1 cos , 28 where the scattering angle Θ is defined by [22] cos 1 1 cos 29 where W¢ and Ω are unit vectors that define the direction of propagation before and after the scattering event, f is the angle between the projection of Ω and W¢ on a plane. It is known that equation (28) can also be expanded in terms of Legendre polynomials as: where l=1, 2, .., L with a 0 =1. It's well known that the problem equation (1) with general boundary conditions equations (2a) and (2b) can be linked to the corresponding source-free problem [26] x In order to solve equation (33) with it boundaries condition equations (33a) and (33b), the PEA method should be followed by assuming that: In the first step, the expression of e(x, μ) and o(x, μ) are straightforward and deduced by considering a a a a , , , 0 o 1 2 3 ¹ and a n =0 for n 0, 1, 2, 3 ¹ in equations (18) and (  The eigenvalue λ is obtained by reducing the series equation (20) to four terms and using the explicit The term of the right side of this equation is composed of four main terms. The first being the isotropic scattering and the other terms are correction due to forward or backward scattering. In the same way, by considering a a a a , , , 0 o 1 2 3 ¹ and a n =0 for n 0, 1, 2, 3 ¹ in the expression of f 1 (x, μ) and f 2 (x, μ) equations (26) and (27) (44) and (45) into equations (24) and (25)

Substituting equations
The equalities of equation ( In the first step the expression of e(x, μ) and o(x, μ) are straightforward deduced by considering a a a a a a , , , , , 0 o 1 2 3 4 5 ¹ and a n =0 for n 0, 1, 2, 3, 4, 5 ¹ in equations (18) and ( -An important result of this work is the determination of diffusion coefficient as function of Chandrasekhar polynomial For instance Thompson scattering have a o =1, a 1 =0, a 2 =1/2 and a n =0 for n 3  using equation (50) and equation (23) it's found that the diffusion coefficient D 1 is given by The formula equation (64) is the same result obtained by Professor Pomraning GC [19] in his work. Secondly with pure triplet scattering (a a a a 1, 0, 0, 0 o 1 2 3 = = = ¹ and a n =0 for n 4  ) using again equation (50) and equation (23), the diffusion coefficient d 1 is given by The formula equation (65) is the same result obtained by Sallah et al [16] by using integral calculus. Hence the current method.  To determine the unknown constants c 1 and c 2 , a weight function W(μ) is introduced to force the boundary conditions equations (16a), (16b) to be fulfilled as [15,27]

Numerical calculations
The diffusion length LL 1 = l was calculated for linear anisotropic scattering and compared with the result obtaind by Ozturk and Anli (2012) see table 1, we can see that our results are comparable with those obtained by P 1 approximation, this indicate that our model behave like the well-known P 1 approximation at lower order scattering.
In order to calculate the neutron partial flux J+ and J− at the medium boundaries, we need to give the expression of the internal energy source function