4.4. Decay Heat

A more detailed decay heat model has been developed for SAS4A/SASSYS‑1 version 5.0 that allows greater flexibility than the previous model and can more accurately represent decay heat during long term transients. The key features of the new model include extending the number of exponential terms that can be used to represent decay heat from six to 24, the inclusion of built-in decay heat parameters from the most recent ANS decay heat power standard, the ability to mix multiple user-supplied and/or built-in decay heat curves within a single region, and the inclusion of a pre-defined non-decay-heat region.

In the following sections, the ANS decay heat standard is described, along with how it is applied to the SAS4A/SASSYS‑1 decay heat model. The implementation of the new model with multiple decay heat curves per region is also discussed in this section, while a detailed description of existing and new input related to the updated decay heat model is provided in Section 4.8.

4.4.1. The ANS Decay Heat Power Standard

The American Nuclear Society (ANS) has published a standard for decay heat power in light water reactors.[4-4] For over 25 years, an extension to the standard has been planned to address fast reactor fuel cycles. In the absence of any such extension, the standard for light-water reactors is used for incorporation into the new decay heat model.

In the standard, decay heat power is defined for each of four different fissionable isotopes in terms of 23-term exponential functions:

(4.4-1)\[f(t) = \sum_{\text{n} = 1}^{23}{\alpha_{\text{n}}e^{- \lambda_{\text{n}}t}} \quad \text{(MeV/Fission-s)}\]

where \(\alpha_{\text{n}}\) is the immediate contribution (in MeV/s) of exponential term \(n\) to the decay power resulting from one fission event, \(\lambda_{\text{n}}\) is the decay constant for term \(n\), and \(t\) is time in seconds after the fission event. The standard currently defines decay heat power for thermal fission in U-235 and Pu-239, fast fission in U-238, and thermal fission in Pu-241.

The ANS standard defines a method for calculating the decay heat power, after shutdown, which results from a known reactor power history. Although not stated in the standard, the method given makes an implicit assumption that the ratio between total power and fission power is fixed when calculating the number of fissions for a given power level. This is not precisely correct, especially for fresh fuel at the start of irradiation. An examination of the parameters given for U-235 suggest that it takes over ten hours of steady-state operation to reach 90% of the equilibrium decay power, and roughly 3½ months to reach 98%. At a fixed total power, then, the fission power will depend on the current level of the decay heat power. The standard avoids the issue by defining decay heat power on a per-fission basis, leaving it up to the user of the standard to provide an appropriate value for the recoverable energy per fission.

In light of this discrepancy and in consideration of the fact that SAS4A/SASSYS‑1 performs decay heat calculations based on fission power (not total power) the methods prescribed in the ANS standard are modified so that they can be adapted for use in SAS4A/SASSYS‑1. Furthermore, the standard describes “adjustments” that can be made to the calculated decay heat power to account for neutron capture in fission products. Because these adjustments are specific to a thermal-spectrum reactor, they are not accounted for here. In addition, the standard describes a method for including contributions from U-239 and Np-239 decay heat power. Because a fast reactor generally has a significant quantity of actinides contributing to decay heat power, this adjustment is also not accounted for here. Instead, user-supplied decay heat curves can be combined with the ANS standard curves to accomplish the same task in a manner that is relevant to the problem being solved.

4.4.1.1. Energy per Fission

The total recoverable energy per fission, \(Q_{\text{t}}\), in a fissionable isotope consists of “prompt” fission energy, \(Q_{\text{f}}\), and the energy from complete decay of fission products, \(Q_{\text{d}}\):

(4.4-2)\[Q_{\text{t}} = Q_{\text{f}} + Q_{\text{d}} \quad \text{(MeV)}\]

The energy from decay of fission products can be calculated by integrating Eq. (4.4-1) over all time:

(4.4-3)\[Q_{\text{d}} = \sum_{\text{n} = 1}^{23}\frac{\alpha_{\text{n}}}{\lambda_{\text{n}}} \quad \text{(MeV)}\]

If the total recoverable energy per fission is known (or input by the user) then the prompt energy per fission can be determined as \(Q_{\text{f}} = Q_{\text{t}} - Q_{\text{d}}\). A useful term, to be used later, is the ratio of the decay energy to the prompt energy,

(4.4-4)\[\beta_{\text{d}} = \frac{Q_{\text{d}}}{Q_{\text{f}}} = \sum_{\text{n} = 1}^{23}\beta_{\text{n}} = \sum_{\text{n} = 1}^{23}\frac{\alpha_{\text{n}}}{\lambda_{\text{n}}Q_{\text{f}}}\]

Therefore, \(\beta_{\text{d}}\) is the total decay heat yield per unit fission heat. Note that for user-defined decay heat parameters, values for \(\beta_{\text{n}}\) and \(\lambda_{\text{n}}\) are defined by user input.

4.4.1.2. Number of Fissions

Eq. (4.4-1) gives the decay heat power as a function of time after a single fission event. The fission rate at time \(t'\) can be written as

(4.4-5)\[F\left( t' \right) = P_{0}\frac{P_{\text{f}}\left( t' \right)}{Q_{\text{f}}K} \quad \text{(Fissions/s)}\]

where \(P_0\) is the nominal reactor power (Watts), \({\text{f}}\) is the relative fission power normalized to the nominal reactor power, and \(K = 1.602177 \times 10^{-13}\) J/MeV. The total number of fission events between time \(t'\) and \(t' + \text{dt}\) is then

(4.4-6)\[F\left( t' \right)\text{dt} = P_{0}\frac{P_{\text{f}}\left( t' \right)}{Q_{\text{f}}K}\text{dt}' \quad \text{(Fissions)}\]

4.4.1.3. Decay Power

The relative decay power at time \(t > t'\) resulting from the fissions represented by Eq. (4.4-1) can be calculated by combining Eq. (4.4-6) with Eq. (4.4-1):

(4.4-7)\[P_{\text{d}}\left( t \right) = \frac{P_{\text{f}}\left( t' \right)}{Q_{\text{f}}}f\left( t - t' \right)\text{dt}' \quad \left( t > t' \right)\]

where \(P_{\text{d}}\) is the relative decay power normalized to the nominal reactor power, \(P_0\). The total decay power at time \(t\) from all fissions prior to time \(t\) can be determined by integrating the above equation:

(4.4-8)\[P_{\text{d}}\left( t \right) = \int_{- \infty}^{t}{\frac{P_{\text{f}}\left( t' \right)}{Q_{\text{f}}}f\left( t - t' \right)\text{dt}'}\]

Note that Eq. (4.4-8) is similar in form to Equation (3) given in the ANS standard, with the exception that the ratio of total power to total recoverable energy, \(P_{\text{t}} / Q_{\text{t}}\), is used in the ANS standard. The equation in the standard inaccurately calculates the fission rate in fresh fuel shortly after startup. While this is likely to be a very small error in most calculations, SAS4A/SASSYS‑1 evaluates decay heat in terms of the fission power and the correct form of Eq. (4.4-8) is more convenient.

4.4.1.4. Contributions from Multiple Isotopes

In the preceding sections, a single fissionable isotope is implied. The determination of decay heat power contributions from multiple fissionable isotopes is straightforward. If \(x_{\text{i}}\) represents the fraction of the total fission power from isotope \(i\), then the normalized fission power from isotope \(i\) is written as

(4.4-9)\[P_{\text{fi}} = x_{\text{i}}P_{\text{f}}\]

and the decay power from isotope \(i\) is

(4.4-10)\[P_{\text{di}}\left( t \right) = x_{\text{i}}\int_{- \infty}^{t}{\frac{P_{\text{f}}\left( t' \right)}{Q_{\text{f}}}f\left( t - t' \right)\text{dt}'}\]

Here, \(f_{\text{i}} \left( t \right)\) represents the decay heat parameters for isotope \(i\).

In practice (as described in the next section on “Implementation”) the decay heat power from isotope \(i\) is calculated using the normalized fission power only, as shown in Eq. (4.4-8), such that

(4.4-11)\[P_{\text{di}}'\left( t \right) = \int_{- \infty}^{t}{\frac{P_{\text{f}}\left( t' \right)}{Q_{\text{f}}}f\left( t - t' \right)\text{dt}'}\]

The fractional fission powers are only applied when total (or integrated) decay heat power is calculated. Therefore, total decay heat power is defined as

(4.4-12)\[P_{\text{d}}\left( t \right) = \sum_{\text{i}}{x_{\text{i}}P_{\text{di}}'\left( t \right)}\]

4.4.2. Implementation

Prior to version 5.0, a channel is assigned to a single decay heat curve by the parameter IDKCRV. Since multiple curves may now be assigned to a given channel, the parameter IDKCRV now refers to which decay heat region a channel belongs. The relationship between decay heat curves and decay heat regions is defined by new user input in the form of two input matrices, DKFRAC and DKANSI, that contains the various \(x_{\text{i}}\) for each region. With this, multiple curves can be assigned (with different weights) to a single region, and a single curve can be assigned to multiple regions. To support older input decks, the matrix defaults to an identity matrix where curve 1 is assigned to region 1, etc. In this case, curves and regions have the same meaning. A description of new and existing input is given in Section 4.8. In the following subsections, only a single decay heat curve is considered. Using the \(x_{\text{i}}\) provided in input, total (normalized) decay heat for each region can be calculated as described by Eq. (4.4-12).

4.4.2.1. Transient Calculations

At the beginning of a transient calculation, the initial decay heat power and fission power are known as a result of the steady-state initialization (see below). From that point on, fission power is determined based on the point (or spatial) kinetics equations. Given the normalized fission power, the normalized decay power can be solved using Eq. (4.4-8). To simplify notation, only a single term from Eq. (4.4-1) is considered, and a term for decay heat (i.e. energy) is introduced, where the decay heat, \(H_{\text{n}}\), is defined by \(P_{\text{d}} \left( t \right)\) = \(\Sigma \ \lambda_{\text{n}} \ H_{\text{n}} \left( t \right)\). At time step \(k\), the time varies from \(t_{\text{k}}\) to \(t_{\text{k}} + \tau\). Term \(n\) of Eq. (4.4-8) (in terms of decay heat) can then be written as

(4.4-13)\[H_{\text{n}}\left( t_{\text{k}} + \tau \right) = e^{- \lambda_{\text{n}}\tau}\int_{-\infty}^{t_{\text{k}}}{P_{\text{f}}\left( t' \right){\beta_{\text{n}}e}^{- \lambda_{\text{n}}\left( t_{\text{k}} - t' \right)}\text{dt}'} + \int_{t_{\text{k}}}^{t_{\text{k}} + \tau}{P_{\text{f}}\left( t' \right){\beta_{\text{n}}e}^{- \lambda_{\text{n}}\left( t_{\text{k}} + \tau - t' \right)}\text{dt}'}\]

Note that the first integral is just the (known) decay heat at the beginning of the time step, \(H_{\text{n}} \left( t_{\text{k}} \right)\). To write the second integral in terms of the time within a time step, let \(t = t' - t_{\text{k}}\). Then the above equation can be written as

(4.4-14)\[H_{\text{n}}\left( t_{\text{k}} + \tau \right) = H_{\text{n}}\left( t_{\text{k}} \right)e^{- \lambda_{\text{n}}\tau} + \beta_{\text{n}}e^{- \lambda_{\text{n}}\tau}\int_{0}^{\tau}{P_{\text{f}}^{k}\left( t_{\text{k}} + t \right)e^{\lambda_{\text{n}}t}\text{dt}}\]

From the point kinetics solution, the fission power over time step \(k\) is represented by a second-order polynomial:

(4.4-15)\[P_{\text{f}}^{k}\left( t_{\text{k}} + t \right) = P_{0} + P_{1}t + P_{2}t^{2}\]

where \(P_0 = P_{\text{f}} \left( t_{\text{k}} \right)\). Using this expansion, the decay heat at the end of the point-kinetics time step can be solved:

(4.4-16)\[H_{\text{n}}\left( t_{\text{k}} + \tau \right) = H_{\text{n}}\left( t_{\text{k}} \right)e^{- \lambda_{\text{n}}\tau} + \frac{\beta_{\text{n}}}{\lambda_{\text{n}}}\left\lbrack P_{0}I_{0}\left( \tau \right) + \frac{P_{1}I_{1}\left( \tau \right)}{\lambda_{\text{n}}} + \frac{P_{2}I_{2}\left( \tau \right)}{\lambda_{\text{n}}^{2}} \right\rbrack\]

where

(4.4-17)\[\begin{split}I_{0}\left( \tau \right) &= 1 - e^{- \lambda_{\text{n}}\tau} \\ I_{1}\left( \tau \right) &= \lambda_{\text{n}}\tau - I_{0}\left( \tau \right) \\ I_{2}\left( \tau \right) &= \left( \lambda_{\text{n}}\tau \right)^{2} - 2I_{1}\left( \tau \right)\end{split}\]

4.4.2.2. Steady-State Initialization

Prior to commencing with the transient calculations, SAS4A/SASSYS‑1 performs a steady-state initialization that includes the determination of initial decay heat. The user has the option of entering a reactor power history as a histogram of relative reactor power. Assuming no initial decay heat is present (i.e. fresh fuel with no fission products) the initial fission power will be equal to the power level of the first histogram. However, as decay heat builds, fission power will decrease to maintain constant power.

Like the decay heat calculations during a transient, the decay heat during the constant-power interval, \(k\), of duration \(\tau_{\text{k}}\) can be written as

(4.4-18)\[H_{\text{n}}\left( t_{\text{k}} + \tau_{\text{k}} \right) = H_{\text{n}}\left( t_{\text{k}} \right)e^{- \lambda_{\text{n}}\tau_{\text{k}}} + \beta_{\text{n}}e^{- \lambda_{\text{n}}\tau_{\text{k}}}\int_{0}^{\tau_{\text{k}}}{P_{\text{f}}^{\text{k}}\left( t_{\text{k}} + t \right)e^{\lambda_{\text{n}}t}\text{dt}}\]

However, only the total power is known over the interval, not the fission power. The above equation can be rewritten as

(4.4-19)\[{H_{\text{n}}\left( t_{\text{k}} + \tau_{\text{k}} \right) = H_{\text{n}}\left( t_{\text{k}} \right)e^{- \lambda_{\text{n}}\tau_{\text{k}}}}{+ \beta_{\text{n}}e^{- \lambda_{\text{n}}\tau_{\text{k}}}\int_{0}^{\tau_{\text{k}}}{\left\lbrack P_{\text{t}}^{k} - \sum_{\text{n} = 1}^{{\text{NDKGRP}}}{\lambda_{\text{n}}H_{\text{n}}\left( t_{\text{k}} + t \right)} \right\rbrack e^{\lambda_{\text{n}}t}\text{dt}}}\]

For a single fissionable isotope, this represents a coupled set of NDKGRP integral equations. When multiple fissionable isotopes contribute to a single decay heat region in the calculation, the number of coupled equations can increase dramatically.

To decouple the terms in Eq. (4.4-19), two assumptions are made. First, the fraction of fission power contributed by each fissionable isotope is assumed to be fixed by the values of \(x_{\text{i}}\). Second, the fission power is assumed to be a constant such that the total power at the end of the constant-power interval matches the value supplied by the user. The first assumption only applies to decay heat regions that use more than one decay heat curve and may not be valid where there is significant depletion or breeding of fissionable isotopes during the course of steady-state initialization.

The second assumption was chosen for a number of reasons. First, by matching power at the end of the initiating power interval, it maintains continuity between the steady-state initialization and the transient calculation. Second, it can be shown (by example) to introduce at most only a few hundredths of a percent error in most cases, and up to a few tenths of a percent error in very unusual situations. Third, it correctly predicts the two bounding cases of fresh fuel (no decay heat) and infinitely-long, steady-state equilibrium.

By assuming a constant (although initially unknown) fission power, the integral in Eq. (4.4-19) can be solved as before:

(4.4-20)\[H_{\text{n}}\left( t_{\text{k}} + \tau \right) = H_{\text{n}}\left( t_{\text{k}} \right)e^{- \lambda_{\text{n}}\tau} + P_{\text{f}}\frac{\beta_{\text{n}}}{\lambda_{\text{n}}}\left( 1 - e^{- \lambda_{\text{n}}t} \right)\]

Note that this is the same as Eq. (4.4-16) with \(P_0 = P_{\text{f}}\) and \(P_1\) and \(P_2\) set to zero. By multiplying by \(\lambda_{\text{n}}\) and summing the above equation over all decay terms, the decay power at the end of the time step is

(4.4-21)\[P_{\text{d}}\left( t_{\text{k}} + \tau \right) = \sum_{\text{n}}{{\lambda_{\text{n}}H}_{\text{n}}\left( t_{\text{k}} \right)e^{- \lambda_{\text{n}}\tau}} + P_{\text{f}}\sum_{\text{n}}{\beta_{\text{n}}\left( 1 - e^{- \lambda_{\text{n}}t} \right)}\]

For simplicity, the additional sum over isotopes with the appropriate \(x_{\text{i}}\), as shown in Eq. (4.4-12), is implied.

Given the user-supplied total power and that \(P_{\text{t}} = P_{\text{f}} + P_{\text{d}}\), the fission power can now be solved:

(4.4-22)\[P_{\text{f}} = \frac{P_{\text{t}} - \sum_{\text{n}}{{\lambda_{\text{n}}H}_{\text{n}}\left( t_{\text{k}} \right)e^{- \lambda_{\text{n}}\tau}}}{1 + \sum_{\text{n}}{\beta_{\text{n}}\left( 1 - e^{- \lambda_{\text{n}}t} \right)}}\]

Again, the sum over isotopes is implied. If the preexisting decay power at the end of the time step (the second term in the numerator) is greater than the user-supplied total power, then the fission power for that histogram is set to zero by the code. This will be the case, for example, if the user specifies a zero-power shutdown period. Otherwise, the numerator represents the fission power if no new decay power is generated during the time step. The denominator then adjusts that value so that the fission power plus accumulated decay power will be equal to the total power at the end of the time step.

Finally, by substituting Eq. (4.4-22) into Eq. (4.4-20), the decay heat can be determined at the end of this step in the histogram. For an infinitely-long power interval, Eq. (4.4-22) simplifies to

(4.4-23)\[P_{\text{f}} = P_{\text{t}}\frac{1}{1 + \beta_{\text{d}}}\]

and the decay heat from Eq. (4.4-20) is

(4.4-24)\[H_{\text{n}}\left( t_{\text{k}} + \tau \right) = P_{\text{f}}\frac{\beta_{\text{n}}}{\lambda_{\text{n}}}\]

These two expressions match the analytical solution exactly.