Prediction of the Average Decay Heat per Fission for MOX Nuclear Fuel ()
1. Introduction
Modeling the inventory of fission products after the fission burst of different fissionable nuclei as well as calculating the decay heat released due to the decay of these fission fragments has been presented in several kinds of literature. Most of these calculations are done depending on two methods exclusively for finding the inventory of fission products, either simulating statistically decomposing all fission fragments by tracking each event individually or solving the Batman equation [1] - [16]. Table 1 shows some of the programs and codes that have been used to solve problems related to the calculation of fission products inventory and decay heat released due to fission fragments.
In this paper, the second method, which depends on solving Bateman’s equation numerically, was used to model the decay of fission fragments. In this case, a specially developed program that can work with large data sets can be used to solve huge numbers of differential equations of type initial value problems in less time, for example, MATLAB software [23], MATHCAD [24], which have a number of built-in tools that make it easier to work with, access, and manipulate a large number of initial value problems differential equations [25] [26].
In this paper, the decay heat of thermal-induced fission of MOX fuel averaged over 100 fissions of its fissile nuclei (U235, Pu239, and Pu241) was calculated. MOX fuel is a mixture of plutonium oxide (from spent fuel) and uranium oxide (natural uranium), with about 3 to 5 percent by weight of plutonium [27]. The main difference between MOX fuels from uranium fuel is that MOX fuel contains more fissile material than uranium fuel. This leads to a difference in nuclear properties which in turn affects the decay heat [28].
Inventories of all fission products need to be well-known with high accuracy as possible because these inventories are required as input to the summation calculations of decay heat. Theoretical calculations of decay heat of different fissile nuclei can play a prominent role in studies of potential new nuclear fuel mixtures [29].
2. Methodology
The source code used in the solution was created and implemented using
Table 1. Some of the Important codes for fission products inventory calculations and decay heat released due to fission fragments.
MATLAB. The elapsed time for each code run was measured and found to be between 30 - 60 minutes using a computer with a 2.4 GHz dual-core i7 CPU. These elapsed times showed a hardware-dependent behavior when using computers with lower speeds. The current code is largely an evolution of the HEATKAU code [22]. Except that the code used does not have an interface, it performs the same calculations in almost the same way, but with some fundamental improvements. All entries have been updated and extracted from ENDF/B-VIII.0 Nuclear Data Library [30] rather than ENDF/B-VII Nuclear Data Library [31]. The latest version of MATLAB (MATLAB R2021b [32] ) was used along with writing the code in a direct way and improving the computational procedures to exclude zero arithmetic operations as much as possible, which made the computation time significantly reduced to approximately 20% of the time of HEATKAU. The list of fissionable nuclei that the code can handle has been updated to include all nuclei for which fission yield data is available in the ENDF/B-VIII.0 Nuclear Data Library as shown in Table 2. Finally, the calculation method has been developed to be able to calculate the average decay heat after fission burst in a mixed fuel like MOX in which more than one type of nuclide such as undergoes fission.
The algorithm of the solution is described in the flowchart in Figure 1. The code follows the following sequences:
Figure 1. The flowchart of the solution.
Table 2. Fission yield data available in ENDF/B-VIII.0 [33].
*Fission yield of Pu-239 is also available at energy of 2 MeV of incident neutron.
Choose the nuclide that will undergo fission from Table 1.
Choose the incident neutron energy based on the availability in ENDF/B-VIII.0 (Table 1), for example thermal fission or 500 keV neutron energy fission.
Extract fission fragment mass distribution (independent fission yield) from ENDF/B-VIII.0 and convert it into matrix form (Z A I) (Input file Y0).
Extract decay modes of all isotopes from ENDF/B-VIII.0.
Use decay modes of all isotopes to generate the branching matrix (Input file B).
Extract decay constants from ENDF/B-VIII.0, and generate decay constants matrix (Input file λ).
Use one of the MATLAB built-in ODE solvers to solve Bateman equation numerically, and generate fission fragments inventory (Output file N(t)).
Extract average energy released per decay of all isotopes from ENDF/B-VIII.0, and generate Alpha, Beta, and Gamma energy matrices (Input files
${E}_{\alpha},{E}_{\beta}$ and
${E}_{\gamma}$ ).
Use the summation method to calculate decay heats of a specific fission event (Output files
${H}_{\alpha},{H}_{\beta}$ and
${H}_{\gamma}$ ).
Calculate the total decay heat
${H}_{\alpha}+{H}_{\beta}+{H}_{\gamma}$ (Output file
${H}_{Total}$ ).
Enter fission ratios for every 100-fission event in the case of mixed fissions like MOX.
Calculate the average decay heat in the case of mixed fissions (Output file
${H}_{Av}$ ).
Sort the decay heat from all isotopes, in descending order, before summation at a specific time after the fission event to conclude the most contributors to decay heat.
The output results can be obtained directly in graphs, in addition to the presence of the complete digital data for the outputs saved anywhere you specify on the device.
The decay rate of each fission fragment after fission burst is simply given by an exponential law. Despite this, a large number of these nuclei are in turn unstable and are depleted after some time by some method of decay. The general solution to this correlative decay that leads to a system of differential equations which must be solved to get the fission fragment inventory was put forward by Bateman in 1910 [34]. So, the number of i^{th} nuclei of fission fragments via thermal fission of three fissile materials U^{235}, Pu^{239}, and Pu^{241}, at time t is
${N}_{i}\left(t\right)$ and can be written as:
${\left[{N}_{i}\left(t\right)={\left({Y}_{i}\right)}_{{\text{U}}^{\text{235}}}\ast {\text{e}}^{\left(-{\lambda}_{i}\ast t\right)}+\underset{1}{\overset{{M}_{{\text{U}}^{\text{235}}}}{{\displaystyle \sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{N}_{j\to i}\left(t\right),j\ne i\right]}_{{\text{U}}^{\text{235}}}$ (1)
${\left[{N}_{i}\left(t\right)={\left({Y}_{i}\right)}_{{\text{Pu}}^{\text{239}}}\ast {\text{e}}^{\left(-{\lambda}_{i}\ast t\right)}+\underset{1}{\overset{{M}_{{\text{Pu}}^{\text{239}}}}{{\displaystyle \sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{N}_{j\to i}\left(t\right),j\ne i\right]}_{{\text{Pu}}^{\text{239}}}$ (2)
${\left[{N}_{i}\left(t\right)={\left({Y}_{i}\right)}_{{\text{Pu}}^{\text{241}}}\ast {\text{e}}^{\left(-{\lambda}_{i}\ast t\right)}+\underset{1}{\overset{{M}_{{\text{Pu}}^{\text{241}}}}{{\displaystyle \sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{N}_{j\to i}\left(t\right),j\ne i\right]}_{{\text{Pu}}^{\text{241}}}$ (3)
where
${\left({Y}_{i}\right)}_{{\text{U}}^{\text{235}}}$,
${\left({Y}_{i}\right)}_{{\text{Pu}}^{\text{239}}}$, and
${\left({Y}_{i}\right)}_{{\text{Pu}}^{\text{241}}}$ are the independent fission yields of the thermal fission of the three fissile nuclei U^{235}, Pu^{239}, and Pu^{2}^{41}, respectively.
${\lambda}_{i}$ is the decay constant of the i^{th} fission product,
${M}_{{\text{U}}^{235}}$,
${M}_{{\text{Pu}}^{\text{239}}}$, and
${M}_{{\text{Pu}}^{\text{241}}}$ are the total number of fission product nuclides of the three fissile nuclei U^{235}, Pu^{239}, and Pu^{2}^{41}, respectively. In addition, the symbol j represents the decays of all j^{th} fission products that give out i^{th} fission products as daughter.
Equations (4) to (6) give the total decay heat generated by fission products
$f\left(t\right)$ as a function of cooling time t after a fission burst for the fissile nuclei U^{235}, Pu^{239}, and Pu^{2}^{41}, respectively.
$f{\left(t\right)}_{{\text{U}}^{\text{235}}}=\underset{1}{\overset{{M}_{{\text{U}}^{235}}}{{\displaystyle \sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\lambda}_{i}\ast {N}_{i}\left(t\right)\ast {E}_{i}$ (4)
$f{\left(t\right)}_{{\text{Pu}}^{\text{239}}}=\underset{1}{\overset{{M}_{{\text{Pu}}^{\text{239}}}}{{\displaystyle \sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\lambda}_{i}\ast {N}_{i}\left(t\right)\ast {E}_{i}$ (5)
$f{\left(t\right)}_{{\text{Pu}}^{\text{241}}}=\underset{1}{\overset{{M}_{{\text{Pu}}^{\text{241}}}}{{\displaystyle \sum}}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}{\lambda}_{i}\ast {N}_{i}\left(t\right)\ast {E}_{i}$ (6)
where
${E}_{i}$ is the average energy released per decay of the i^{th} isotope, and is given by
${E}_{i}={\left({E}_{\alpha}+{E}_{\beta}+{E}_{\gamma}\right)}_{i}$ (7)
where
${E}_{\alpha},{E}_{\beta}$ and
${E}_{\gamma}$ are the average energy released per heavy particles-, beta-, and gamma-decays of the i^{th} isotope, and extracted from ENDF/B-VIII.0.
The total decay heat due to thermal fissions is the sum of the decay heat from each fissile nuclide in the fuel. Therefore, the burst function as a function of time
$f\left(t\right)$ of MOX fuel can be calculated by using the weight percentages
$\%\text{\hspace{0.05em}}\text{\hspace{0.05em}}i$ of the three fissile nuclei as:
$f{\left(t\right)}_{\text{MOX}}=\frac{{\displaystyle \sum \%\text{\hspace{0.05em}}\text{\hspace{0.05em}}i\ast f{\left(t\right)}_{i}}}{100}$ (8)
To predict the average heat of decay per fission of the MOX nuclear fuel, we need to calculate the weight percentages
$\%\text{\hspace{0.05em}}\text{\hspace{0.05em}}i$ of the three fission nuclei, U^{235}, Pu^{239}, and Pu^{241}, respectively. As shown in Figure 2, it is assumed that the MOX fuel contains a mix of depleted uranium and a mix of plutonium isotopes. Also as mentioned in the figure, for 1000 kg of the MOX fuel, the weights of the three fissile nuclei U^{235}, Pu^{239}, and Pu^{241} are 1.9, 7.1, and 34.3 kg, respectively.
It is known that in the operation of reactors not designed to burn plutonium, no more than 30% of the MOX groups can be loaded along with 70% of conventional uranium fuel (UOX)) in the core of these reactors [28]. Moreover, Figure 3 shows the weight composition of 30% MOX fuel and 70% Natural fuel for a 1000 kg of metal fuel.
The weight ratios of the three fissile nuclei U^{235}, Pu^{239}, and Pu^{241} can be calculated by using Equations (9)-(11):
$\%\text{}{U}^{235}=\frac{Weight\left({U}^{235}\right)}{Weight\left({U}^{235}\right)+Weight\left(P{u}^{239}\right)+Weight\left(P{u}^{241}\right)}\ast 100$ (9)
Figure 2. Weight composition (kg) of MOX fuel for a 1000 kg of metal fuel.
Figure 3. Weight composition (kg) of 30% MOX (300 kg) and 70% Natural fuel (700 kg) for a 1000 kg of metal fuel.
$\%\text{}P{u}^{239}=\frac{Weight\left(P{u}^{239}\right)}{Weight\left({U}^{235}\right)+Weight\left(P{u}^{239}\right)+Weight\left(P{u}^{241}\right)}\ast 100$ (10)
$\%\text{}P{u}^{241}=\frac{Weight\left(P{u}^{241}\right)}{Weight\left({U}^{235}\right)+Weight\left(P{u}^{239}\right)+Weight\left(P{u}^{241}\right)}\ast 100$ (11)
So, for 1000 kg metal fuel, the calculated weights of U^{235} = 0.57 + 5.04 = 5.61 kg, Pu^{239} = 2.13 kg, and Pu^{241} = 10.29 kg, respectively.
$\%\text{}{U}^{235}=\frac{5.61}{5.61+2.13+10.29}\ast 100=31\text{\%}$ (12)
$\%\text{}P{u}^{239}=\frac{2.13}{5.61+2.13+10.29}\ast 100=12\text{\%}$ (13)
$\%\text{}P{u}^{241}=\frac{5.61}{5.61+2.13+10.29}\ast 100=57\text{\%}$ (14)
Finally, Figure 4 shows the weight % of the three fissile materials of 30% MOX and 70% natural fuel. Therefore, Equation (8) can be rewritten as:
$f{\left(t\right)}_{\text{MOX}}=\frac{31\ast f{\left(t\right)}_{{\text{U}}^{\text{235}}}+12\ast f{\left(t\right)}_{{\text{Pu}}^{\text{239}}}+57\ast f{\left(t\right)}_{{\text{Pu}}^{\text{241}}}}{100}$ (15)
Figure 4. Weight % of fissile materials of 30% MOX and 70% Natural fuel.
3. Results
The first-order differential equations (Equations (1)-(3)) of every fission product in each fissile material in MOX fuel were solved numerically to find the exact time-dependent distribution of fission product. Time-dependent distribution of fission product until cooling time t = 10^{7} sec after burst thermal fission with (0.0235 eV neutrons) for the three fissile nuclei U^{235}, Pu^{239}, and Pu^{241} were calculated numerically, and complete fission product inventories are shown in Figures 5-7 , respectively.
By using the aforementioned method (Summation Method), the burst function
$f\left(t\right)$ (decay heat) of the fission products after the fission burst were calculated for each of the three fissile nuclei in the MOX fuel U^{235}, Pu^{239}, and Pu^{241} using Equations (4)-(6), respectively. In addition, the average decay heat of MOX fuel
$f{\left(t\right)}_{\text{MOX}}$ can be calculated by using Equation (15). Moreover, Figure 8 shows a comparison of the calculated total decay heat of the three fissionable nuclei U^{235}, Pu^{239}, and Pu^{241}, with the predicted total decay of MOX fuel averaged over 100 fissions.
In addition, one of the most important aspects associated with the prediction of the decay heat is identifying the contribution of different fission product nuclides to the total heat of decay. Calculating the contribution at multiple times after the fission process is important in assessing the fission product’s decay heat and evaluating the decay data. The overall contribution of specific nuclides that contribute the total thermal energy released after thermal fission of any fissile isotope is concluded from the contribution matrix F (Equation (16)), which was calculated by integration of all specific contributions between starting time t_{1} and ending time t_{2}. Figures 9-11 showed the overall contribution of nuclides that contribute more than 2% of the total thermal energy of the decay of fission products of thermal fission of U^{235}, Pu^{239}, and Pu^{241}, respectively.
$F=\left[\begin{array}{c}{f}_{1}\\ {f}_{2}\\ \vdots \\ {f}_{M}\end{array}\right]=\left[\begin{array}{c}{\displaystyle {\int}_{{t}_{1}}^{{t}_{2}}{f}_{1}\left(t\right)\text{d}t}\\ {\displaystyle {\int}_{{t}_{1}}^{{t}_{2}}{f}_{2}\left(t\right)\text{d}t}\\ \vdots \\ {\displaystyle {\int}_{{t}_{1}}^{{t}_{2}}{f}_{M}\left(t\right)\text{d}t}\end{array}\right]$ (16)
Figure 5. Complete inventory of fission products decay of fission-due to 0.0235 eV neutron-of U-235.
Figure 6. Complete inventory of fission products decay of fission-due to 0.0235 eV neutron of Pu-239.
Figure 7. Complete inventory of fission products decay of fission-due to 0.0235 eV neutron-of Pu-241.
Figure 8. Calculated total decay heat of the three fissionable nuclei U^{235}, Pu^{239} vs the average decay heat of MOX fuel.
Figure 9. Fission fragments with contribution more than 2% on overall total decay heat of fission due to 0.0235 eV neutron of U-235.
Figure 10. Fission fragments with contribution more than 2% on overall total decay heat of fission due to 0.0235 eV neutron of Pu-239.
Figure 11. Fission fragments with contribution more than 2% on overall total decay heat of fission-due to 0.0235 eV neutron-of Pu-241.
4. Conclusions
The above-mentioned code flowchart shown in Figure 1, gave good capability to estimate the decay heat of different fuels using the fission yields and the nuclear decay data, which are available from ENDF nuclear data libraries. Since these libraries are subjected to comprehensive screening procedures and fruitful verifications before releasing to the scientific community. Therefore, it is crucial to use the latest version of ENDF to obtain a good decay heat estimation. In the present study, ENDF/B-VIII has been used instead of ENDF/B-VII.
The present code is a powerful tool to help researchers to implement the heat of nuclear decay after a fission explosion, which is extremely important for safety and design studies related to spent fuel management, transportation or storage. The maximum values for the highest decay heats always set the rules in the design goals of safety margins. The code can be used to calculate fission fragments inventory and decay heats of all available fission yield in ENDF/B-VIII.0. In addition, the capabilities of the present code have been verified by calculating the decay heat due to the decay of fission fragments of 31 fissionable nuclei listed in (Table 2).
The decay heat of U^{235}, Pu^{239}, and Pu^{241} due to thermal induced neutron fission is calculated using the present code. Moreover, the decay heat for the MOX fuel using the equivalent content (weighted values) of its fissile nuclei is calculated. The results show that for the case of thermal fission, the average decay heat per fission calculated for the MOX fuel is lower than the decay heat per thermal fission of U^{235} and Pu^{239}, but higher than that per thermal fission of Pu^{241}. Moreover, the most contributors at different cooling times after fission burst (10, 100, 1000, 10,000, 100,000 seconds) to the decay heat have been identified for U^{235}, Pu^{239} and Pu^{241}.