Corresponding author: Sergey V. Bedenko ( bedenko@tpu.ru ) Academic editor: Yury Korovin
© 2018 Igor V. Shamanin, Sergey V. Bedenko, Vladimir N. Nesterov, Igor O. Lutsik, Anatoly A. Prets.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Shamanin IV, Bedenko SV, Nesterov VN, Lutsik IO, Prets AA (2018) Solution of neutrontransport multigroup equations system in subcritical systems. Nuclear Energy and Technology 4(1): 7985. https://doi.org/10.3897/nucet.4.29837

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 BNAB78, 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 28group 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.
Subcritical system, neutron distribution function, neutron transport, multigroup diffusion approximation.
Neutron transport calculation methods are used in reactor physics, for process monitoring during nuclear fuel fabrication and reprocessing, to determine the irradiated fuel burnup, in biological shielding design, and to improve the radiation measurement procedures in the nuclear material accounting and control system.
There are a lot of dedicated programs and codes that implement the Monte Carlo method based on selection of neutron interaction probabilities. The internationally recognized program packages include MMKKENO and MMKC (IPPE, Russia), MCU (Kurchatov Institute, Russia), MCNP (USA), KENO3D (USA), MONK (Great Britain) and others. These programs make it possible to obtain highly accurate results thanks to a 3D geometry 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 crosssections of the neutron interaction with the material nuclei were determined.
Deterministic methods based on differential and finitedifference representation of the Boltzmann equation (
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.
We shall write the Boltzmann equation (
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; F^{(}^{i}^{)}, F^{(}^{k}^{)} are the neutron flux densities in respective groups; Σ_{a}^{(}^{i}^{)} is the macroscopic neutron absorption crosssection for the i^{th} group; Σ_{R}^{i}^{→}^{k}, Σ_{R}^{k}^{→}^{i} are the macroscopic crosssections 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 ε^{(}^{i}^{)} = ʃχ_{fp} (E)dE); ν_{f}^{(}^{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 crosssection 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}^{(}^{i}^{)} 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:
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 (
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 BNAB78 and BNAB93 systems of constants, a = 0.965(0.8 + 0.083∙ν^{*}_{f}) and b = 2.245/(0.8+0.083∙ ν^{*}_{f})^{2}. 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
S_{sf} = ν^{*}_{sf} R_{sf} λN = ν^{*}_{sf} R_{sf} (ln 2 /T_{1/2}) N,
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 halflife period.
The spectral distribution of spontaneous fission neutrons is
S_{sf} (E) = ν^{*}_{sf} R_{sf} (ln 2 /T_{1/2}) N×χ_{sf} (E).
Spontaneous and induced fission neutrons are specific in being timecorrelated unlike the (α,n)neutrons and, besides, can be approximated by one and the same function
where a = Т relates to ν^{*}_{sf} through a = 0.48 + 0.2∙(1 + ν^{*}_{sf})^{1/2}. 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} (
The average energy of the spontaneous fission neutron spectrum (
E ^{*} = (3/2)T + E_{f} = 1.48 + 0.3×(1 + ν^{*}_{sf})^{1/2}.
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 (
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 jtype material nuclei; σ^{αn}_{j} is the integral macroscopic crosssection 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 (
The dependence σ^{αn}_{j} was obtained based on data in (
The spectral distribution of (α,n)neutrons is
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 (
where the constants а and b depend on the type of the nuclei i and j.
The investigations conducted in (
The reference data used for the spectra acquisition were the results obtained in (
To be able to compare correctly the calculated data we have obtained with (
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 neutrons for the considered composition is ~2.35∙10^{4} 1/s/g of PuO_{2}.
Figs
The normalized distribution of the ^{238}Pu spontaneous fission neutron spectrum (see Fig.
where с = 0.36; а = 0.84; and b = 3.63. The average spectrum energy is E^{*} = 2.02 MeV, and the average number of neutrons is ν^{*}_{sf} = 2.21.
The normalized distribution of the neutron spectrum for the ^{238}PuO_{2} (α,n)reactions (see Fig.
where а = 0.78; b = 2.5. The average spectrum energy is E^{*} = 2.5 MeV.
In subcritical multiplying systems, the equation written in a diffusion approximation for neutrons in all energy groups (
DΔF(r,t) – Σ_{a}F(r,t) – ν_{f} Σ_{f} F(r,t) < 0,
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
DΔF(r,t) – Σ_{a}F(r,t) – ν_{f} Σ_{f} F(r,t) + S (r,t) = 0.
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 (
(2)
The first term in (2), describing the neutron leakage, is determined from the relation
D ^{(} ^{i} ^{)}ΔF^{(}^{i}^{)} = – D^{(}^{i}^{)}B^{2(}^{i}^{)}F^{(}^{i}^{)}, (3)
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
F^{(}^{i}^{)}_{j} = f (F^{(–1)}_{j}_{–1}, F^{(0)}_{j}_{–1}, ... , F^{(}^{k}^{)}_{j}_{–1}, ... , F^{(26)}_{j}_{–1}), k ≠ i, (4)
where j is the number of the iteration, beginning with the first one, by transforming system (2), with regard for (3), as follows
By expressing the value of the flux density in the i^{th} group from (5), we get
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
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
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 BNAB78 group constants (
The result of the calculation for the 28group spectrum Ф(E), neutr./(cm^{2}∙s) is shown in Figs
Neutronic characteristics of the PuO_{2} fuel composition.
Parameter  Code  
Diffusion approximation ( 
ANISN ( 
Diffusion approximation ( 
SCALE4.3, KENOV.a (ENDF/BV)  MCNP4b (JENDL3.2)  
Mass, kg  5.65  5.65  24.30  25.42  24.97 
k _{eff}  0.39  0.35  1.0  1.0  1.0 
1/(1−k_{eff})  1.64  1.55  0  0  0 
It is shown in (
F(r,E) = F(E)×r^{2}/(R + 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.
A range of group constants was obtained in the iteration process organization which can be used to estimate k_{eff}. The expression defining k_{eff} has the form (
The constants were averaged based on relations (
Table
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 (
The spectra S (E) and χ(E) used in the calculations can be calculated in dedicated codes such as SOURCE4C (
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 28group 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.