1. Introduction

Undoubtedly, increased accuracy of calculations on radioactive waste - the most important of which is depleted (or spent) nuclear fuel - effectively contributes to the development of nuclear research through the fuel cycle. At the top of these calculations is the decay heat of spent nuclear fuel during both short and long cooling periods after reactor shutdown to support safe transport and storage. The decay heat, which is considered one of the most important characteristics of spent nuclear fuel SNF, is affected by several factors, including the type of nuclear fuel and the initial fuel enrichment rate. Also, the fuel burnup and its cooling period directly affect the decay heat. Knowing the total decay heat along with the decay heats of components such as the decay heat caused by actinides or those caused by fission products is of great importance when handling a SNF for transportation or storage to ensure the safety of the SNF as well as to protect the surrounding environment, which helps in safety and shielding calculations [1]. Both modeling and simulating fresh nuclear fuel consumption and then modeling the decay of spent fuel during the cooling period - during operation and after the shutdown of the reactor - require the use of reliable and accurate nuclear data. So, it made sense for this kind of calculation and analysis to be done whenever newer and more accurate versions of nuclear data libraries emerged. In this research, all nuclear data was acquired from the Evaluated Nuclear Data File (ENDF) library, most recent version, ENDF/B-III.0. Examples of nuclear data obtained from this library are the neutron capture cross sections, fission yields of fissionable isotopes, half-lives or decay constants, branching ratios of different types of decay, and mean decay energies (alpha, beta, and gamma) of all radioactive isotopes in the fuel [2].

Many previous studies have published calculations and analyses of decay heat of spent nuclear fuel using nuclear data derived from the release of ENDF/B-VII [3 to 10]. All calculations for this research relied on the new ENDF/BVIII. 0 Nuclear Data Library. As this library fully observed the new IAEA standards, it also benefited from experimental and modern theoretical data obtained in the United States and Europe. Arguably, the integrated verification performance of the new ENDF/B-VIII.0 library has been improved compared to the previous ENDF/B-VII.1 library. And what concerns us in this research is that the improving fission-yields and decay data have been largely done [11].

In some recent research work, nuclear data has been acquired from NEW ENDF/B-VIII.0 Nuclear Data Library both in fresh fuel burnup calculations during the reactor’s operational period and during the decay of SNF in the cooling period after the reactor shutdown to identify nuclides contributing to decay heat [12-13]. Decay heat studies are based on predicting the isotopic abundances of radioisotopes constituting the decay heat present at the shutdown of the reactor, based on the operational history of the reactor and then following up on decay heat generation over a short and long period of time. Also, some research calculates only the decay heat of a single event for fission of different fissile isotopes (decay heat after thermal or rapid fission). These calculations, while important, do not give a complete picture of the entire nuclear fuel behavior [14-15].

The first step of calculations (fresh fuel burnup calculations) is usually done using codes based on the Mont Carlo method, the most important and accurate of which is the MCNP code [16]. The second part of the calculations, relating to the conclusion and prediction of isotopic content of spent fuel, used a wide spectrum of codes that draw nuclear data from a nuclear data library and rely on different methods of calculation. In this research, as well as relying entirely on NEW ENDF/B-VIII.0 Nuclear Data Library, MATLAB was also used to generate radionuclide isotopic abundances based on irradiation or burnup history. The reason for choosing MATLAB was because it had several built-in solvers dedicated to problem-solving that had a huge number of differential equations, especially very stiff ones [17], exactly as in isotopic decay of SNF. The utilization of programs that offer accurate solutions plays a significant role in advancing research on nuclear fuel cycles. Furthermore, it facilitates the development of suitable measures for health and safety, as well as regulatory strategies for the handling of radioactive materials. This includes achieving both short and long-term safety requirements for the cooling process after a reactor’s shutdown.

2. Methodology

To calculate the evolution of material composition and radionuclide isotopic abundances during 1,350 days expected core lifetime burnup, simulated irradiation of fresh uranium fuel enriched up to 4.5% was performed using the Monte Carlo N-Particle eXtended (MCNPX) code. The MCNPX code relies on randomized numbers in statistical samples where the program records the neutrons pathways emitted by any source and determines the ability of these neutrons to react [16]. The simulation process begins with the preparation of the input data file, which includes materials for the nuclear fuel used, as well as fuel geometrical specifications, as well as source data. After running the code, the required outputs, the final isotopic contents, are obtained after fresh fuel is burned after required duration (1,350 days divided into cycles). The MCNPX code was combined with homemade MATLAB R2023a code for decay heat calculations.

The homemade code was written to numerically solve the differential equations for the decay and accumulation of radioactive isotope chains of spent fuel (Bateman equations) [18].

The time dependent concentration of i^th nuclei of N_i(t) in SNF at any time “t” after reactor shutdown can be written as:

where C_i represents the initial concentration of i^th isotope of SNF just after shutdown of the reactor, λ_i is the decay constant of the i^th nuclei, and j represents decays of all j^th radioactive isotope that give out i^th isotope as a result of direct decay.

Now we can define N_j→i (t) as:

where, b_j→i is the branching ratios of i^th nuclide per decay of nuclide j^th, and λ_j is the decay constant of the j^th nuclei.

Now, by rewriting the previous two equations in the form of differential equations, in addition to converting all the inputs to the form of matrices that are suitable for the MATLAB program, as follows:

As a second step, MATLAB’s Build-in ODE Solver ODE15s was used. ODE15s is recommended for very stiff ordinary differential equations as in the case of radioactive decay equations for isotopic components of SNF containing a wide range of decay constants, where the derivatives of the variables vary with time by orders of magnitude. The concept of ODE15s lies in the Backward Differentiation Formulae, or BDF. This algorithm is unique in that it combines a variable order and variable step size, utilizing both a progressive formula and a method of assessing local errors that is typical of constant step size formulas. The code generated 1,800 time points (over a full life time of 300 years of cooling) and calculated the isotopic contents for all radioactive and non-radioactive isotopes of SNF at each point. To illustrate mathematically the problem of stiffness, in some ODE problems, the step size taken by the solver is forced to a level that is extremely small compared to the integration interval and these step sizes can be so small that passing a short time of solution would require millions of evaluations. This can cause the ODE solver to fail. And even if it succeeds, it will take a very long time to do so. Therefore, recent versions of MATLAB have focused on solving this type of problems with continuous improvements. Research shows how efficient ODE15s are in terms of numerical accuracy in cases where very steep gradients appear in solution [19-20].

The calculations presented in this paper were based solely on nuclear data obtained from the most recent edition of the Evaluated Nuclear Data File, ENDF/B-III. 0, as previously stated in the introductory section. These nuclear data were neutrons capture cross sections, decay constants/half-lives, branching ratios for all types of decay, average decay energies (Eα, Eβ, and Eγ) per decomposition, and fission yields.

3. Results and Discussion

The ratios of decay heat calculated in the entire cooling period (300 years) for all radioisotopes affecting decay heat are organized in the Table 1 by dividing them into actinides (Ac, Th, Pa, U, Np, Pu, Am, and Cm isotopes) and non-actinides (108 radioisotopes), or gaseous and nongaseous. The heat generated by actinides (52%) during the entire cooling period is roughly the same as that found in non-actinides (48%), while on the other hand, the gaseous radioisotopes consisting mainly of xenon (¹³³Xe, and ¹³⁵Xe) and radon (²¹⁷Rn, ²¹⁸Rn, ²¹⁹Rn, ²²⁰Rn, and ²²²Rn) isotopes contribute to the decay heat by less than 1%.

In this research, the results are drawn on a log-log scale so that we can present and analyze as many data points as possible clearly, especially when there are low and high decay heat values. Fig. 1 presents decay heat power vs cooling time of separate actinides during the entire cooling period. Although the results showed the presence of radioisotopes of 7 actinides, some of them had a very weak effect on decay heat, particularly actinium isotopes (²²⁵Ac, ²²⁷Ac, and ²²⁸Ac), which were neglected because they were infinitesimal. The effect of both uranium isotopes (9 isotopes from ²³²U to ²⁴⁰U respectively) and neptunium isotopes (²³⁶Np, ²³⁷Np, ²³⁸Np, and ²³⁹Np) was also weak and less than 0.5% of the total decay temperature. On the other hand, the main effect of the decay heat caused by actinides is the isotopes of americium (²⁴¹Am, ²⁴²Am, ²⁴³Am, and ²⁴⁴Am), whose effect is close to 28% of the total decay heat during the entire cooling period, followed by the isotopes of plutonium (9 isotopes from ²³⁶Pu to ²⁴⁴Pu respectively), which affect more than 18% of the total decay heat, and finally the isotopes of curium (7 isotopes from ²⁴²Cm to ²⁴⁸Cm respectively) contribute about 5% of the total decay heat. Also, by checking Fig. 1, we can say that the effect of both uranium and neptunium isotopes is significant at the beginning of the cooling period (beginning of storage of SNF) and quickly decreases only within a month or a few months. Also, curium isotopes begin with a significant effect on decay heat and extend to a few years of cooling (possibly 20 years). Finally, although the effect of plutonium and americium isotopes on decay heat has been low in short cooling periods, they are the main contributors on decay heat at the long-term cooling level. Fig. 2 tells us that at the beginning of the cooling period, the non-actinide isotopes were significantly more than several times the effect of actinides in decay heat, and then began to decrease at a greater rate than the actinides, each of which is equal after cooling for about 40 years, and then decreased at a very high rate.

Fig. 1

Actinide isotopes decay heats as a function of storage time.

Fig. 2

Actinides vs Non- Actinides decay heats as a function of storage time.

As observed in Fig. 3, α-, β-, and γ- decay heats were calculated based on both isotopic abundances of radioactive isotopes during the entire cooling period and on Eα, Eβ, and Eγ derived from ENDF/B-VIII.0. By looking at Fig. 3, we can say that the main decay heats during the first year of cooling was caused by beta and gamma decay. In the cooling period of more than 40 years, the decay heat caused by alpha decay was the main contributor to total decay heat.

Fig. 3

α, β, and γ decay heats as a function of storage time.

Fig. 4 illustrates the dominant radionuclides, across 1, 10, 100, and 300 years of cooling (SNF storage). There were 154 isotopes that contributed to the production of heat at any time during the decay from the onset of cooling up to 300 years, of which 46 isotopes represented 7 actinides and 108 isotopes represented fission fragments. Also, these 154 isotopes, of which 66 were alpha emitters, 145 beta emitters, and 137 gamma emitters. As we conclude from the figure, during the four cooling periods after 1, 10, 100, and 300 years, about 10 to 13 radioisotopes contribute more than 90% of the decay heat, while the remaining isotopes (about 140 isotopes) combined contribute less than 10% of the decay heat. From the beginning of cooling of the SNF to a full year, the isotopes contributing 90% to decay heat were ranked downwards from the most contribution to the least as follows: ¹⁴⁴Pr, ¹⁰⁶Rh, ⁹⁵Nb, ¹³⁴Cs, ⁹⁵Zr, ¹⁴⁰La, ²⁴²Cm, ⁹¹Y, ¹⁰³Ru, ⁹⁰Y, ¹⁴⁴Ce, ⁸⁹Sr, and ²³⁹Np. And as we can see, in this period of cooling one isotope of actinides contributes to the decay heat, which is ²⁴²Cm. In the first 10 years of cooling, the isotopes contributing about 90% to decay heat were ranked downwards from the most contribution to the least as follows: ⁹⁰Y, ²⁴¹Am, ²³⁸Pu, ¹⁴⁴Pr, ¹⁰⁶Rh, ¹³⁴Cs, ¹³⁷Cs, ²⁴⁴Cm, ⁹⁰Sr, ⁹⁵Nb, ⁹⁵Zr, ²⁴⁰Pu, and ²⁴²Cm. The increased effect of actinides on decay heat has also been observed. In addition to ²⁴²Cm, ²⁴¹Am, ²³⁸Pu, ²⁴⁴Cm, and ²⁴⁰Pu also contributed to the decay heat. In the cooling periods of 100 and 300 years we found that the isotopes contributing to the decay heat did not change, although there were slight changes in the contribution ratios. In those cooling periods, isotopes that contributed more than 90% to the decay heat were as follows: ⁹⁰Y, ²⁴¹Am, ²³⁸Pu, ¹⁴⁴Pr, ¹⁰⁶Rh, ¹³⁴Cs, ¹³⁷Cs, ²⁴⁴Cm, ⁹⁰Sr, ⁹⁵Nb, ⁹⁵Zr, ²⁴⁰Pu, and ²⁴²Cm. Also, the effect of actinides on the decay heat increased, as ²³⁸Pu, ²⁴⁰Pu, ²⁴¹Am, ²⁴²Cm, and ²⁴⁴Cm contributed more than 38% to the total decay heat.

Fig. 4

The dominant radioisotopes contribute decay heat during 1, 10, 100, and 300 years of cooling.

4. Conclusions

From the consistent contribution of isotopes to thermal decomposition, shown in Fig. 4, it is concluded that there is no need to calculate cooling times up to 105 years in some research as that is unnecessarily time-consuming in the calculations. Fission products ¹⁴⁴Pr, ¹⁰⁶Rh, and ¹³⁴Cs are the isotopes with the greatest impact ever in intermediate storage periods (about 50% of total decay heat), while on the other side the fission product ⁹⁰Y and both actinide isotopes ²⁴¹Am and ²³⁸Pu (44.2% at 100 y cooling and 56.2% at 300 y cooling) are the largest ever in long storage periods. The effect of fission products ¹⁴⁴Pr, ¹⁰⁶Rh, and ¹³⁴Cs also persists in long cooling periods but to a lesser extent than in short periods (24.5% at 100 y cooling and 17% at 300 y cooling).

The isotopes of the actinides americium and plutonium only have an important contribution to decay heat in long cooling periods, while other actinides do not have a significant effect in any cooling periods except for the isotopes of curium, which begin to rapidly decline after the first 10 years of storage.