Solution of neutron-transport multigroup equations system in subcritical systems *

An iteration method has been implemented to solve a neutron transport equation in a multigroup diffusion approximation. A thermoelectric generator containing plutonium dioxide, used as a source of thermal and electric power in spacecraft, was studied. Neutron yield and multigroup diffusion approximation data was used to obtain a continuous and group distribution of neutron flux density spectra in a subcritical multiplying system. Numerical multigroup approaches were employed using BNAB-78, a system of group constants, and other available evaluated nuclear data libraries (ROSFOND, BROND, BNAB, EXFOR and ENDSF). The functions of neutron distribution in the zero iteration for the system of multigroup equations were obtained by approximating an extensive list of calculated and experimental data offered by the EXFOR and ENDSF nuclear data libraries. The required neutronic functionals were obtained by solving a neutron transport equation in a 28-group diffusion approximation. The calculated data was verified. The approach used is more efficient in terms of computational efforts (the values of the neutron flux density fractions converge in the third iteration). The implemented technique can be used in nuclear and radiation safety problems.

that models the object and simulation of real neutron transport in material.
High accuracy of verified Monte Carlo codes is achieved through a great deal of computational resources consumed and is defined primarily by the error with which the cross-sections of the neutron interaction with the material nuclei were determined.
Deterministic methods based on differential and finite-difference representation of the Boltzmann equation (Bell and Glesston 1974) do not require a great deal of computational resources.Such methods simplify the main simulation phase that consists in solution of the transport equation.The major deterministic methods include the discrete ordinate method, the characteristics method, the first collision probabilities method, the PN-method, and multigroup diffusion approximation as its simplest and most widespread form.
The accuracy of the solution (Ф(r,t), Σ(r,E), k eff ) obtained using multigroup diffusion approximation is not high, but the method in question offers a reasonably good approximation and is still used in Russian and foreign programs, e.g. in the USA's MCNP and KENO and in many others.
An iteration process is implemented in the study to solve a system of multigroup diffusion equations which is used to determine the spectral and integral neutronic characteristics of a subcritical multiplying system.
The considered system is a thermoelectric generator containing 5.65 kg of PuO 2 , used as a heat and electricity source in spacecraft.

Computational model
We shall write the Boltzmann equation (Bell and Glesston 1974) in a 28-group diffusion approximation: XXX where i is the number of the neutron energy group for which the equation has been written; k is the energy group number; D (i) is the neutron diffusion constant for the i th group; u (i) is the neutron velocity corresponding to the average energy of the i th group; k) are the neutron flux densities in respective groups; Σ a (i) is the macroscopic neutron absorption cross-section for the i th group; Σ R i→k , Σ R k→i are the macroscopic cross-sections of the neutron transfer from the i th group to the underlying k th group (from the overlying k th group to the considered i th group) respectively; ε fp (i) is the fraction of prompt neutrons in the spectrum χ fp (E) falling within the i th group (for all spectra ε (k) is the average number of prompt neutrons per the induced fission event for neutrons in the k th group; Σ f (k) is the macroscopic fission cross-section for neutrons in the k th group; ε fd (i),n is the fraction of delayed neutrons in the i th group; λ j n is the decay constant for the j th group of delayed neutrons of the nuclide n; C j n (r,t) is the concentration of the delayed neutron precursor nuclei; S sf is the integral yield of spontaneous fission neutrons; ε sf (i) is the fraction of neutrons in the spectrum χ sf (E); ε αn is the fraction of neutrons in the spectrum χ αn (E); S sf is the integral yield of spontaneous fission neutrons; S αn is the integral yield of neutrons for the (α,n) reaction; and Q (i) (r,t) is the intensity of the external neutron source.Equation (1) needs to be supplemented with a boundary condition and the initial condition: XXX To organize the iteration process, it is required to define the spectra of induced fission neutrons, spontaneous fission neutrons and the (α,n)-reaction neutrons: χ f (E) (ε (i) = ∫χ(E)dE), χ sf (E) and χ αn (E) respectively.And the selection of the neutron distribution function and the spectrum S(E,t,r) = Q(r,t)•χ(E) is an important simulation phase that makes it possible to improve considerably the accuracy of the solutions sought.
The normalized spectrum of induced fission neutrons χ f (E) has been calculated in an assumption that the neutron yield can be approximated by Watt's spectral function (Abagyan et al. 1981) XXX 3 2 exp( / 4) ( ) exp( / ) sh .
The parameters a and b depend on the type of the fissionable nuclide and the energy of the neutrons that have caused the fission and are clearly associated with the energy and fissile isotope averaged quantity ν * f .In the BNAB-78 and BNAB-93 systems of constants, a = 0.965(0.8+ 0.083 The parameters are used to determine the average fission neutron spectrum energy E * = 1.5a + 0.25a 2 b.
The yield of spontaneous fission neutrons is where ν * sf is the average number of spontaneous fission neutrons; R sf is the spontaneous fission branching ratio; λ is the decay constant; N = N 0 ×exp(-λt) is the current inventory of the spontaneously fissile nuclide; and T 1/2 is the half-life period.
The spectral distribution of spontaneous fission neutrons is Spontaneous and induced fission neutrons are specific in being time-correlated unlike the (α,n)-neutrons and, besides, can be approximated by one and the same function XXX The parameter b is defined by the temperature Т and the average kinetic energy E f per spontaneous fission product nucleon and is equal to b = 4E f /a 2 = 3.04/a 2 (Vlaskin et al. 2015).
The average energy of the spontaneous fission neutron spectrum (Vlaskin et al. 2015, Wilson et al. 2002) is The latter relation makes it possible to obtain the average energy of the spectrum χ sf (E) equal to the experimental value for the most of fissile nuclides.
The yield of neutrons in the (α,n)-reaction (Bulanenko 1980, Vukolov andChukreev 1987) where ω i = (1 -ω j ) is the fraction of the α-emitter (material i); ω j is the fraction of the material on which the reaction runs (material j); q α i is the yield of α-particles from the material of the type i; n j is the concentration of the j-type material nuclei; σ αn j is the integral macroscopic cross-section of the (α,n)-reaction; ε ij (E) = (-dE/dx) is the stopping power of the α-particle; E α is the average spectrum energy for α-particles; and B αn is the threshold of the (α,n)-reaction.
The data on the mass stopping power ε ij (E) were calculated using the SRIM (The Stopping and Range of Ions in Matter) code (Ziegler et al. 2010).
The dependence σ αn j was obtained based on data in (Hansen et al. 1967, Bair and Gomez del Campo 1979, West and Sherwood 1982, Heaton et al. 1989, Murata and Shibata 2002) recommended for solving such type of problems, and the compilation and the analysis were performed in (Vlaskin et al. 2015, Wilson et al. 2002, Bulanenko 1980, Shamanin et al. 2010).
The spectral distribution of (α,n The normalized distribution of the neutron spectrum χ αn (E) was obtained through approximation of experimental (EXFOR − Experimental Nuclear Reaction Data) and evaluated (ENDSF − Evaluated Nuclear Structure Data File) data (Gorshkov 1962, Herold 1968, Anderson and Neff 1969, Taherzadeh 1971, Taherzadeh and Gingo Peter 1972, Arkhipov et al. 1972, Bair and Butler 1973, Anderson 1980) by the function of the type XXX where the constants а and b depend on the type of the nuclei i and j.
The investigations conducted in (Gorshkov 1962, Herold 1968, Anderson and Neff 1969, Taherzadeh 1971, Taherzadeh and Gingo Peter 1972, Arkhipov et al. 1972, Bair and Butler 1973, Anderson 1980) had the purpose of obtaining highly accurate values of the neutron yields and spectra for 238 PuO 2 .The measurement error does not exceed 5%, and the compilation and the evaluation were performed in (Vlaskin et al. 2015, Bulanenko 1980, Shamanin et al. 2010).Full texts and the results of the experiments are available in the EXFOR and ENDSF nuclear data libraries.
The reference data used for the spectra acquisition were the results obtained in (Taherzadeh and Gingo Peter 1972) which present accurate characteristics of the neutron spectra for the fuel composition PuO 2 with a mass of 5.65 kg.The PuO 2 -containing thermoelectric generator studied in (Taherzadeh and Gingo Peter 1972) is used as a heat and electricity source in space satellites.The measurements were performed to design the radiation protection for the spacecraft's main components.
To be able to compare correctly the calculated data we have obtained with (Taherzadeh and Gingo Peter 1972) and to have the desired distribution as the result, the isotopic composition of Pu (Pu, %: 238 -80.27; 239 -15.87; 240 -3.022; 241 -0.643, 242 -0.132) and the presence of light impurity elements (Li, Be, B, C, F, Na, Al, Si and others) were taken into account.
Therefore, the yield of neutrons due to the spontaneous fission of Pu isotopes is equal to approximately 2.48•10 3 1/s/g of PuO 2 ; 1.18•10 4 1/s/g of PuO 2 thanks to the (α,n)-reaction on oxygen, and 9.21•10 3 1/s/g of PuO 2 thanks to reactions on impurities.The total yield of neu-

Solution of a system of multigroup neutron diffusion equations by iteration method
In subcritical multiplying systems, the equation written in a diffusion approximation for neutrons in all energy groups (Bell and Glesston 1974) has the form (Golovatskiy et al. 2010) and, hence, the value of the neutron flux density Ф(r,t) in such system is expected to become zero.The value F(r,t) in subcritical systems is known not to become zero but has a particular value characterized by neutron sources S(r,t) and by the neutronic parameters of the system.
With regard for the sources of neutrons, the above equation will become stationary Therefore, equation ( 1) in a stationary form, with no regard for delayed neutrons and in the absence of an external source, can be written as follows (Golovatskiy et al. 2010 The first term in (2), describing the neutron leakage, is determined from the relation where B 2(i) is the geometrical parameter for the i th group.
Organizing the iteration process requires the system of equations to be composed as follows where j is the number of the iteration, beginning with the first one, by transforming system (2), with regard for (3), as follows XXX XXX XXX By expressing the value of the flux density in the i th group from (5), we get XXX The system of equations is transformed as follows: ; . . . . . . . . . . . . . . . . . .
All parameters in system (6a) are known, except the flux densities for the previous iteration (k)  j-1 and, therefore, the sum

∑
This sum defines the number of the neutrons produced in the current generation during the fission of nuclei by all neutrons in the previous generation, except the neutrons in the i th group.
To start the iteration process with the zero iteration, the flux density for neutrons in the i th group is determined from (2), with regard for (3), by the relation XXX where the number of the neutrons produced in the current generation during the fission of nuclei by all neutrons in the previous generation is defined as equal to unity, that is XXX After the neutron flux spectrum in the zero iteration is determined using system (6), the iteration process is implemented.The calculation was performed for a subcritical spherical system of 238 PuO 2 with a mass of 5.65 kg.
The system of the BNAB-78 group constants (Abagyan et al. 1981) and recommendations in (Voropayev et al. 1979) were used, which made it possible to compare correctly the results obtained with the data presented in (Taherzadeh and Gingo Peter 1972).

Calculation results and discussions
The result of the calculation for the 28-group spectrum Ф(E), neutr./(cm 2 •s) is shown in Figs. 3 and 4. We shall note that the values of the neutron flux density fractions converge in the third iteration, with the data obtained being fairly accurate (see Fig. 4 and Table 1).
It is shown in (Taherzadeh and Gingo Peter 1972) that, at a distance of 50 cm and more, the PuO 2 fuel composition with a mass of 5.65 kg can be viewed as a point source, and the neutron flux decay takes place by the law ~1/R 2 : where F(E) is the neutron flux density spectrum; r is the radius of the PuO 2 sphere; and R is the distance from the sphere surface to the detection point.
Fig. 3 presents the flux density spectrum for the neutrons produced by the PuO 2 point source.It can be seen from the figure that curves 3 (authors) and 4 (Taherzadeh and Gingo Peter 1972) have a shift.The explanation is that curve 3 is a group distribution, and curve 4 is a con- .
Table 1 presents a comparison of the found values k eff and the mass with those obtained using ANISN (Taherzadeh and Gingo Peter 1972), SCALE-4.3(KENO-V.a)and MCNP4b.

Conclusion
The continuous and group distributions of spectra and the system subcriticality value have been obtained using neutron yield data and a multigroup diffusion approximation.The results fairly accurately agree with the data in (Gorshkov 1962, Herold 1968, Anderson and Neff 1969, Taherzadeh 1971, Taherzadeh and Gingo Peter 1972, Arkhipov et al. 1972, Bair and Butler 1973, Anderson 1980) and with the results obtained in the codes ANISN (Taherzadeh and Gingo Peter 1972), MCNP (JENDL-3.2.) and SCALE-4.3 (KENO-V.a,ENDF/B-V).
The spectra S(E) and χ(E) used in the calculations can be calculated in dedicated codes such as SOURCE-4C (Vlaskin et al. 2015) and NEDIS-2m (Wilson et al. 2002) which allow finding the intensity and the spectrum of neutrons produced in the (α,n) reactions and during spontaneous fission and preparing a file with output data in such form as is convenient for solving the neutron transport equation using such codes as MCU, MCNP and SCALE.
The spectra S(E) and χ(E) were obtained by approximating an extensive list of calculated and experimental data offered by the EXFOR and ENDSF nuclear data libraries, and the required functionals (Ф(r,E), k eff ) were found as part of solving the neutron transport equation in a 28-group diffusion approximation.
The approach used is more efficient in terms of computational efforts and the nuclear data bank storage cost and can be used to solve applied nuclear and radiation safety problems.
Figs. 1 and 2 show the sought-after distributions S(E) and spectra χ(E) used to solve the transport equation in a 28-group diffusion approximation.The normalized distribution of the 238 Pu spontaneous fission neutron spectrum (see Fig.2) is XXX

Table 1 .
Neutronic characteristics of the PuO 2 fuel composition.