.. _section-4.4:

Decay Heat
----------

A more detailed decay heat model has been developed for |SAS|
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 |SAS| 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 :numref:`section-4.8`.

.. _section-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:

.. math::
	:label: eq-ANS-decay

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

where :math:`\alpha_{\text{n}}` is the immediate contribution (in MeV/s) of
exponential term :math:`n` to the decay power resulting from one fission
event, :math:`\lambda_{\text{n}}` is the decay constant for term :math:`n`, and :math:`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
|SAS| 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 |SAS|.
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.

.. _section-4.4.1.1:

Energy per Fission
^^^^^^^^^^^^^^^^^^

The total *recoverable* energy per fission, :math:`Q_{\text{t}}`, in a
fissionable isotope consists of "prompt" fission energy, :math:`Q_{\text{f}}`,
and the energy from complete decay of fission products, :math:`Q_{\text{d}}`:

.. math::
	:label: eq-ANS-total

	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:`eq-ANS-decay` over all time:

.. math::
	:label: eq-4.4-1

	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
:math:`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,

.. math::
	:label: eq-ANS-beta

	\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, :math:`\beta_{\text{d}}` is the total decay heat yield per unit fission
heat. Note that for user-defined decay heat parameters, values for
:math:`\beta_{\text{n}}` and :math:`\lambda_{\text{n}}` are defined by user input.

.. _section-4.4.1.2:

Number of Fissions
^^^^^^^^^^^^^^^^^^

:eq:`eq-ANS-decay` gives the decay heat power as a function of time after a single
fission event. The fission rate at time :math:`t'` can be written as

.. math::
	:label: eq-4.4-2

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

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

.. math::
	:label: eq-ANS-fission

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

.. _section-4.4.1.3:

Decay Power
^^^^^^^^^^^

The relative decay power at time :math:`t > t'` resulting from the
fissions represented by :eq:`eq-ANS-decay` can be calculated by combining
:eq:`eq-ANS-fission` with :eq:`eq-ANS-decay`:

.. math::
	:label: eq-4.4-3

	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 :math:`P_{\text{d}}` is the relative decay power normalized to the
nominal reactor power, :math:`P_0`. The total decay power at time :math:`t`
from all fissions prior to time :math:`t` can be determined by integrating the
above equation:

.. math::
	:label: eq-ANS-DecayPower

	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:`eq-ANS-DecayPower` 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, :math:`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, |SAS| evaluates
decay heat in terms of the fission power and the correct form of
:eq:`eq-ANS-DecayPower` is more convenient.

.. _section-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 :math:`x_{\text{i}}` represents the
fraction of the total fission power from isotope :math:`i`, then the
normalized fission power from isotope :math:`i` is written as

.. math::
	:label: eq-4.4-4

	P_{\text{fi}} = x_{\text{i}}P_{\text{f}}

and the decay power from isotope :math:`i` is

.. math::
	:label: eq-4.4-5

	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, :math:`f_{\text{i}} \left( t \right)` represents the decay heat parameters for
isotope :math:`i`.

In practice (as described in the next section on "Implementation") the
decay heat power from isotope :math:`i` is calculated using the normalized
fission power only, as shown in :eq:`eq-ANS-DecayPower`, such that

.. math::
	:label: eq-4.4-6

	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

.. math::
	:label: eq-ANS-PowerSum

	P_{\text{d}}\left( t \right) = \sum_{\text{i}}{x_{\text{i}}P_{\text{di}}'\left( t \right)}

.. _section-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
:math:`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 :numref:`section-4.8`. In the
following subsections, only a single decay heat curve is considered.
Using the :math:`x_{\text{i}}` provided in input, total (normalized) decay heat
for each region can be calculated as described by :eq:`eq-ANS-PowerSum`.

.. _section-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:`eq-ANS-DecayPower`. To simplify notation, only a single term from :eq:`eq-ANS-decay` is
considered, and a term for decay heat (i.e. energy) is introduced, where
the decay heat, :math:`H_{\text{n}}`, is defined by
:math:`P_{\text{d}} \left( t \right)` = :math:`\Sigma \ \lambda_{\text{n}} \ H_{\text{n}} \left( t \right)`. At time
step :math:`k`, the time varies from :math:`t_{\text{k}}` to :math:`t_{\text{k}} + \tau`.
Term :math:`n` of :eq:`eq-ANS-DecayPower` (in terms of decay heat) can then be written as

.. math::
	:label: eq-4.4-7

	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, :math:`H_{\text{n}} \left( t_{\text{k}} \right)`. To write the
second integral in terms of the time within a time step, let
:math:`t = t' - t_{\text{k}}`. Then the above equation can be written
as

.. math::
	:label: eq-4.4-8

	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 :math:`k`
is represented by a second-order polynomial:

.. math::
	:label: eq-4.4-9

	P_{\text{f}}^{k}\left( t_{\text{k}} + t \right) = P_{0} + P_{1}t + P_{2}t^{2}

where :math:`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:

.. math::
	:label: eq-ANS-H

	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

.. math::
	:label: eq-ANS-I

	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)

.. _section-4.4.2.2:

Steady-State Initialization
^^^^^^^^^^^^^^^^^^^^^^^^^^^

Prior to commencing with the transient calculations, |SAS|
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, :math:`k`, of duration :math:`\tau_{\text{k}}`
can be written as

.. math::
	:label: eq-4.4-10

	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

.. math::
	:label: eq-ANS-H-total

	{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:`eq-ANS-H-total`, two assumptions are made. First, the
fraction of fission power contributed by each fissionable isotope is
assumed to be fixed by the values of :math:`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:`eq-ANS-H-total` can be solved as before:


.. math::
	:label: eq-ANS-H-const

	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:`eq-ANS-H` with
:math:`P_0 = P_{\text{f}}` and :math:`P_1` and :math:`P_2` set to
zero. By multiplying by :math:`\lambda_{\text{n}}` and summing the above equation
over all decay terms, the decay *power* at the end of the time step is

.. math::
	:label: eq-4.4-11

	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
:math:`x_{\text{i}}`, as shown in :eq:`eq-ANS-PowerSum`, is implied.

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

.. math::
	:label: eq-ANS-PowerInit

	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:`eq-ANS-PowerInit` into :eq:`eq-ANS-H-const`, the decay heat can
be determined at the end of this step in the histogram. For an
infinitely-long power interval, :eq:`eq-ANS-PowerInit` simplifies to

.. math::
	:label: eq-4.4-12

	P_{\text{f}} = P_{\text{t}}\frac{1}{1 + \beta_{\text{d}}}

and the decay heat from :eq:`eq-ANS-H-const` is

.. math::
	:label: eq-4.4-13

	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.