Print
Solution of neutron-transport multigroup equations system in subcritical systems*
expand article infoIgor V. Shamanin, Sergey V. Bedenko, Vladimir N. Nesterov, Igor O. Lutsik, Anatoly A. Prets
‡ National Research Tomsk Polytechnic University, Tomsk, Russia
Open Access

Abstract

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.

Keywords

Subcritical system, neutron distribution function, neutron transport, multigroup diffusion approximation.

Research status

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), KENO-3D (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 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), keff) 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 PuO2, 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:

(1)

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 ith group; u(i) is the neutron velocity corresponding to the average energy of the ith group; F(i), F(k) are the neutron flux densities in respective groups; Σa(i) is the macroscopic neutron absorption cross-section for the ith group; ΣRik, ΣRki are the macroscopic cross-sections of the neutron transfer from the ith group to the underlying kth group (from the overlying kth group to the considered ith group) respectively; εfp(i) is the fraction of prompt neutrons in the spectrum χfp (E) falling within the ith 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 kth group; Σf(k) is the macroscopic fission cross-section for neutrons in the kth group; εfd(i),n is the fraction of delayed neutrons in the ith group; λjn is the decay constant for the jth group of delayed neutrons of the nuclide n; Cjn (r,t) is the concentration of the delayed neutron precursor nuclei; Ssf 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); Ssf 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 (Abagyan et al. 1981)

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∙ν*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.25a2b.

The yield of spontaneous fission neutrons is

Ssf = ν*sf Rsf λN = ν*sf Rsf (ln 2 /T1/2) N,

where ν*sf is the average number of spontaneous fission neutrons; Rsf is the spontaneous fission branching ratio; λ is the decay constant; N = N0×exp(–λt) is the current inventory of the spontaneously fissile nuclide; and T1/2 is the half-life period.

The spectral distribution of spontaneous fission neutrons is

Ssf (E) = ν*sf Rsf (ln 2 /T1/2) N×χsf (E).

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

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 Ef per spontaneous fission product nucleon and is equal to b = 4Ef/a2 = 3.04/a2 (Vlaskin et al. 2015).

The average energy of the spontaneous fission neutron spectrum (Vlaskin et al. 2015, Wilson et al. 2002) is

E * = (3/2)T + Ef = 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 (Bulanenko 1980, Vukolov and Chukreev 1987) is

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; nj is the concentration of the j-type material nuclei; σαnj 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 σαnj 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)-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 (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

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 238PuO2. 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 PuO2 with a mass of 5.65 kg. The PuO2-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∙103 1/s/g of PuO2; 1.18∙104 1/s/g of PuO2 thanks to the (α,n)-reaction on oxygen, and 9.21∙103 1/s/g of PuO2 thanks to reactions on impurities. The total yield of neutrons for the considered composition is ~2.35∙104 1/s/g of PuO2.

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.

Figure 1. 

Spectral distribution of PuO2 neutrons: 1 – spontaneous fission; 2 – desired, (α,n)-reaction (a = 0.85, b = 2.50); 3 – combined desired, {csf (E) + can (E)}; 4 – combined (Taherzadeh and Gingo Peter 1972); 5 – combined (approximation of the EXFOR and ENDSF library data).

The normalized distribution of the 238Pu spontaneous fission neutron spectrum (see Fig. 2) is

Figure 2. 

Normalized distribution of the continuous (a) and the 28-group spectrum (b) of neutrons: 1 – spontaneous fission; 2 – (α,n)-reactions; 3 – total vector.

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 238PuO2 (α,n)-reactions (see Fig. 2) is

where а = 0.78; b = 2.5. The average spectrum energy is E* = 2.5 MeV.

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)

DΔF(r,t) – ΣaF(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) – ΣaF(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 (Golovatskiy et al. 2010):

(2)

The first term in (2), describing the neutron leakage, is determined from the relation

D ( i )ΔF(i) = – D(i)B2(i)F(i), (3)

where B2(i) is the geometrical parameter for the ith 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), ki, (4)

where j is the number of the iteration, beginning with the first one, by transforming system (2), with regard for (3), as follows

(5)

By expressing the value of the flux density in the ith group from (5), we get

(6)

The system of equations is transformed as follows:

(6a)

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 ith group.

To start the iteration process with the zero iteration, the flux density for neutrons in the ith 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 238PuO2 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./(cm2∙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).

Figure 3. 

Spectrum of the neutron flux density for the sphere of 238PuO2 with a mass of 5.65 kg: 1 – spontaneous fission, continuous and 28-group; 2 – (α,n)-reactions; 3, 4 – total spectrum.

Figure 4. 

Spectrum of the neutron flux density at a distance of 50 cm from the PuO2 sphere: 1 – spontaneous fission (authors); 2 − (α,n)-reactions (authors); 3 – total spectrum (authors); 4 – total spectrum of PuO2 (Monte Carlo method (Taherzadeh and Gingo Peter 1972)).

Neutronic characteristics of the PuO2 fuel composition.

Parameter Code
Diffusion approximation (Vlaskin et al. 2015, Voropayev et al. 1979) ANISN (Taherzadeh and Gingo Peter 1972) Diffusion approximation (Vlaskin et al. 2015, Voropayev et al. 1979) SCALE-4.3, KENO-V.a (ENDF/B-V) MCNP4b (JENDL-3.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−keff) 1.64 1.55 0 0 0

It is shown in (Taherzadeh and Gingo Peter 1972) that, at a distance of 50 cm and more, the PuO2 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/R2:

F(r,E) = F(Er2/(R + r)2,

where F(E) is the neutron flux density spectrum; r is the radius of the PuO2 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 PuO2 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 continuous distribution (the spectra were obtained by Monte Carlo method (Taherzadeh and Gingo Peter 1972)). The same shift can be also observed in Fig. 2. Besides, the fuel composition in (Taherzadeh and Gingo Peter 1972) has a more complex geometry than that being modeled and includes components of the PuO2-containing capsule structures. In the region of energies above 4 MeV, the spectrum χsf (E) is higher than used in (Taherzadeh and Gingo Peter 1972), and the contributor to curve 3 is the displaced group spectrum χαn(E).

A range of group constants was obtained in the iteration process organization which can be used to estimate keff. The expression defining keff has the form (Abagyan et al. 1981)

The constants were averaged based on relations (Abagyan et al. 1981)

Table 1 presents a comparison of the found values keff 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), keff) 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.

References

  • Abagyan LP, Bazazyants NO, Nikolayev MN, Tsibulya AM (1981) Group constants for the calculation of reactors and protection. Atomizdat Publ., Moscow, 139 pp. [in Russian]
  • Arkhipov VA, Gorshkov GV, Grebenskii BS (1972) Neutron radiation of 238PuO2 containing different amounts of 18O. Soviet Atomic Energy 32(4): 347–348.
  • Bell D, Glesston S (1974) Theory of nuclear reactors. Atomizdat Publ., Moscow, 489 pp. [in Russian]
  • Bulanenko VI (1980) Neutron yield of (α,n) reaction on oxygen. Soviet Atomic Energy 47(1): 531–534.
  • Golovatskiy AV, Nesterov VN, Shamanin IV (2010) Organization of the iterative process in the numerical reconstruction of the neutron spectrum in a multiplying system with a graphite retarder. Russian Physics Journal 53(11): 10–14 [in Russian]
  • Gorshkov VA (1962) Neutron Yield form the (α,n) Reaction in Be, B, C, O, F, Mg, Al, Si, and Granite Irradiation with Polonium a-particles. Soviet Atomic Energy 13: 123–135.
  • Hansen LF, Anderson JD, McClure JW (1967) The (α,n) cross section on 17O and 18O between 5 and 12,5 MeV. Nuclear Physics A 98(1): 23–32.
  • Heaton R, Lee H, Skensved P, Robertson BC (1989) Neutron Production from Thick-Target (α,n) Reactions. Nuclear Instruments and Methods in Physics Research A 276(3): 529–538. https://doi.org/10.1016/0168-9002(89)90579-2
  • Shamanin IV, Bulanenko VI, Bedenko SV (2010) Neutron Radiation Field of the Irradiated Ceramic Nuclear Fuel of Different Types. Izvestiya vuzov. Yadernaya Energetika 2: 97–103 [in Russian]
  • Vlaskin GN, Khomyakov YS, Bulanenko VI (2015) Neutron Yield of the Reaction (α,n) on Thick Targets Comprised of Light Elements. Atomic Energy 117(5): 357–365. https://doi.org/10.1007/s10512-015-9933-5
  • Voropayev AI, Van’kov AA, Vozyakov VV (1979) Group neutron cross sections for fission and radiation capture of transactinides. VANT. Ser. Yadernye konstanty 3(34): 34–60. [in Russian]
  • Vukolov VA, Chukreev FE (1987) Neutron yield for chemical compounds of actinides. Soviet Atomic Energy. 62(4): 271–276.
  • West D, Sherwood AC (1982) Measurements of Thick-Target (α,n) Yields from Light Elements. Annals of Nuclear Energy. 1982, v. 9, 551–577.
  • Wilson WB, Perry RT, Charlton WS (2002) SOURCES 4С: A Сode for Calculating (α,n), Spontaneous Fission, and Delayed Neutron Sources and Spectra. Los Alamos National Laboratory report LA UR-02-1839.
  • Ziegler JF, Ziegler MD, Biersack JP (2010) SRIM – The Stopping and Ranges of Ions in Matter. Nuclear Instruments and Methods in Physics Research B 268: 1818–1823. https://doi.org/10.1016/j.nimb.2010.02.091

* Russian text published: Izvestiya vuzov. Yadernaya Energetika (ISSN 0204-3327), 2017, n.4, pp. 38-49.