9. Depletion

When materials in a system are subject to irradiation over a long period of time, nuclides within the material will transmute due to nuclear reactions as well as spontaneous radioactive decay. The time-dependent process by which nuclides transmute under irradiation is known as depletion or burnup. To accurately analyze nuclear systems, it is often necessary to predict how the composition of materials will change since this change results in a corresponding change in the solution of the transport equation. The equation that governs the transmutation and decay of nuclides inside of an irradiated environment can be written as

\[\begin{aligned} \frac{dN_i(t)}{dt} = &\sum\limits_j \underbrace{\left [ \underbrace{f_{j \rightarrow i} \int_0^\infty dE \; \sigma_j (E, t) \phi(E,t)}_\text{transmutation} + \underbrace{\lambda_{j\rightarrow i}}_\text{decay} \right ] N_j(t)}_{\text{Production of nuclide }i\text{ from nuclide }j} \\ &- \underbrace{\left [\underbrace{\int_0^\infty dE \; \sigma_i (E,t) \phi(E,t)}_\text{transmutation} + \underbrace{\sum\limits_j \lambda_{i\rightarrow j}}_\text{decay} \right ] N_i(t)}_{\text{Loss of nuclide }i} \end{aligned}\]

where \(N_i\) is the density of nuclide \(i\) at time \(t\), \(\sigma_i\) is the transmutation cross section for nuclide \(i\) at energy \(E\), \(f_{j \rightarrow i}\) is the fraction of transmutation reactions in nuclide \(j\) that produce nuclide \(i\), and \(\lambda_{j \rightarrow i}\) is the decay constant for decay modes in nuclide \(j\) that produce nuclide \(i\). Note that we have not included the spatial dependence of the flux or cross sections. As one can see, the equation simply states that the rate of change of \(N_i\) is equal to the production rate minus the loss rate. Because the equation for nuclide \(i\) depends on the nuclide density for possibly many other nuclides, we have a system of first-order differential equations. To form a proper initial value problem, we also need the nuclide densities at time 0:

\[N_i(0) = N_{i,0}.\]

These equations can be written more compactly in matrix notation as

(1)\[\frac{d\mathbf{n}}{dt} = \mathbf{A}(\mathbf{n},t)\mathbf{n}, \quad \mathbf{n}(0) = \mathbf{n}_0\]

where \(\mathbf{n} \in \mathbb{R}^n\) is the nuclide density vector, \(\mathbf{A}(\mathbf{n},t) \in \mathbb{R}^{n\times n}\) is the burnup matrix containing the decay and transmutation coefficients, and \(\mathbf{n}_0\) is the initial density vector. Note that the burnup matrix depends on \(\mathbf{n}\) because the solution to the transport equation depends on the nuclide densities.

9.1. Numerical Integration

A variety of numerical methods exist for solving Eq. (1). The simplest such method, known as the “predictor” method, is to divide the overall time interval of interest \([0,t]\) into smaller timesteps over which it is assumed that the burnup matrix is constant. Let \(t \in [t_i, t_i + h]\) be one such timestep. Over the timestep, the solution to Eq. (1) can be written analytically using the matrix exponential

\[\mathbf{A}_i = \mathbf{A}(\mathbf{n}_i, t_i) \\ \mathbf{n}_{i+1} = e^{\mathbf{A}_i h} \mathbf{n}_i\]

where \(\mathbf{n}_i \equiv \mathbf{n}(t_i)\). The exponential of a matrix \(\mathbf{X}\) is defined by the power series expansion

\[e^{\mathbf{X}} = \sum\limits_{k=0}^\infty \frac{1}{k!} \left ( \mathbf{X} \right )^k\]

where \(\mathbf{X}^0 = \mathbf{I}\). A series of so-called predictor-corrector methods that use multiple stages offer improved accuracy over the predictor method. The simplest of these methods, the CE/CM algorithm, is defined as

\[\mathbf{n}_{i+1/2} = e^{\frac{h}{2}\mathbf{A}(\mathbf{n}_i, t_i)} \mathbf{n}_i \\ \mathbf{n}_{i+1} = e^{h \mathbf{A}(\mathbf{n}_{i+1/2},t_{i+1/2})} \mathbf{n}_i\]

Here, the value of \(\mathbf{n}\) at the midpoint is estimated using \(\mathbf{A}\) evaluated at the beginning of the timestep. Then, \(\mathbf{A}\) is evaluated using the densities at the midpoint and used to integrate over the entire timestep.

Our aim here is not to exhaustively describe all integration methods but rather to give a few examples that elucidate the main considerations one must take into account when choosing a method. Generally, there is a tradeoff between the accuracy of the method and its computational expense. The expense is driven almost entirely by the time to compute a transport solution, i.e., to evaluate \(\mathbf{A}\) for a given \(\mathbf{n}\). Thus, the cost of a method scales with the number of \(\mathbf{A}\) evaluations that are performed per timestep. On the other hand, methods that require more evaluations generally achieve higher accuracy. The predictor method only requires one evaluation and its error converges as \(\mathcal{O}(h)\). The CE/CM method requires two evaluations and is thus twice as expensive as the predictor method, but achieves an error of \(\mathcal{O}(h^2)\). An exhaustive description of time integration methods and their merits can be found in the thesis of Colin Josey.

OpenMC does not rely on a single time integration method but rather has several classes that implement different algorithms. For example, the openmc.deplete.PredictorIntegrator class implements the predictor method, and the openmc.deplete.CECMIntegrator class implements the CE/CM method. A full list of the integrator classes available can be found in the documentation for the openmc.deplete module.

9.2. Matrix Exponential

As we saw in the previous section, numerically integrating Eq. (1) requires evaluating one or more matrix exponentials. OpenMC uses the Chebyshev rational approximation method (CRAM), which was introduced in a series of papers by Pusa (1, 2), to evaluate matrix exponentials. In particular, OpenMC utilizes an incomplete partial fraction (IPF) form of CRAM that provides a good balance of numerical stability and efficiency. In this representation the matrix exponential is approximated as

\[e^{\mathbf{A}t} \approx \alpha_0 \prod\limits_{\ell=1}^{k/2} \left ( \mathbf{I} + 2 \text{Re} \left ( \widetilde{\alpha}_\ell \left (\mathbf{A}t - \theta_\ell \mathbf{I} \right )^{-1} \right ) \right )\]

where \(k\) is the order of the approximation and \(\alpha_0\), \(\widetilde{\alpha}_\ell\), and \(\theta_\ell\) are coefficients that have been tabulated for orders up to \(k=48\). Rather than computing the full approximation and then multiplying it by a vector, the following algorithm is used to incrementally apply the terms within the product (note that the original description of the algorithm presented by Pusa contains a typo):

  1. \(\mathbf{n} \gets \mathbf{n_0}\)
  2. For \(\ell = 1, 2, \dots, k/2\)
    • \(\mathbf{n} \gets \mathbf{n} + 2\text{Re}(\widetilde{\alpha}_\ell (\mathbf{A}t - \theta_\ell)^{-1})\mathbf{n}\)
  3. \(\mathbf{n} \gets \alpha_0 \mathbf{n}\)

The \(k\)th order approximation for CRAM requires solving \(k/2\) sparse linear systems. OpenMC relies on functionality from scipy.sparse.linalg for solving the linear systems.

9.3. Data Considerations

In principle, solving Eq. (1) using CRAM is fairly simple: just construct the burnup matrix at various times and solve a set of sparse linear systems. However, constructing the burnup matrix itself involves not only solving the transport equation to estimate transmutation reaction rates but also a series of choices about what data to include. In OpenMC, the burnup matrix is constructed based on data inside of a depletion chain file, which includes fundamental data gathered from ENDF incident neutron, decay, and fission product yield sublibraries. For each nuclide, this file includes:

  • What transmutation reactions are possible, their Q values, and their products;
  • If a nuclide is not stable, what decay modes are possible, their branching ratios, and their products; and
  • If a nuclide is fissionable, the fission products yields at any number of incident neutron energies.

9.3.1. Transmutation Reactions

OpenMC will setup tallies in a problem based on what transmutation reactions are available in a depletion chain file, so any arbitrary number of transmutation reactions can be tracked. The pregenerated chain files that are available on https://openmc.org include the following transmutation reactions: fission, (n,\(\gamma\)), (n,2n), (n,3n), (n,4n), (n,p), and (n,\(\alpha\)).

9.3.2. Capture Branching Ratios

Some (n,\(\gamma\)) reactions may result in a product being in either the ground or a metastable state. The most well-known example is capture in Am241, which can produce either Am242 or Am242m. Because the metastable state of Am242m has a significantly longer half-life than the ground state, it is important to accurately model the branching of the capture reaction in Am241. This is complicated by the fact that the branching ratio may depend on the incident neutron energy causing capture.

OpenMC does not currently allow energy-dependent capture branching ratios. However, the depletion chain file does allow a transmutation reaction to be listed multiple times with different branching ratios resulting in different products. Spectrum-averaged capture branching ratios have been computed in LWR and SFR spectra and are available at https://openmc.org/depletion-chains.

9.3.3. Fission Product Yields

Fission product yields (FPY) are also energy-dependent in general. ENDF fission product yield sublibraries typically include yields tabulated at 2 or 3 energies. It is an open question as to what the best way to handle this energy dependence is. OpenMC includes three methods for treating the energy dependence of FPY:

  1. Use FPY data corresponding to a specified energy.
  2. Tally fission rates above and below a specified cutoff energy. Assume that all fissions below the cutoff energy correspond to thermal FPY data and all fission above the cutoff energy correspond to fast FPY data.
  3. Compute the average energy at which fission events occur and use an effective FPY by linearly interpolating between FPY provided at neighboring energies.

The method can be selected through the fission_yield_mode argument to the openmc.deplete.Operator constructor.

9.3.4. Power Normalization

The reaction rates provided OpenMC are given in units of reactions per source particle. For depletion, it is necessary to compute an absolute reaction rate in reactions per second. To do so, the reaction rates are normalized based on a specified power. A complete description of how this normalization can be performed is described in Normalization of Tally Results. Here, we simply note that the main depletion class, openmc.deplete.Operator, allows the user to choose one of two methods for estimating the heating rate, including:

  1. Using fixed Q values from a depletion chain file (useful for comparisons to other codes that use fixed Q values), or
  2. Using the heating or heating-local scores to obtain an nuclide- and energy-dependent estimate of the true heating rate.

The method for normalization can be chosen through the normalization_mode argument to the openmc.deplete.Operator class.