.. _section-4.6:

Solution Methods
----------------

The |SAS| point kinetics equations are solved using a technique
model by Kaganove [4-13] and extended by Fuller [4-14]. The decay-heat
model is solved according to the method developed by Woodruff for the
SAS3A computer code [4-1]. In general, the |SAS| neutronics
solution methods have been adopted from SAS3A.

.. _section-4.6.1:

Point Kinetics
~~~~~~~~~~~~~~

The solution of the point kinetics model equations is carried out
assuming that (1) the reactivity has been specified, and (2) the number
of delayed-neutron precursors is limited to six or less. This latter
assumption requires the interpretation of the precursor yields and decay
constants as average values, typical for a particular blend of isotopes.

Given these assumptions, the first step in the solution of the point
kinetics equations is the integration of the delayed-neutron precursor
energy production balance equations given by :eq:`4.3-1` over a time
interval. In |SAS|, this time interval is taken as the main
time step. If conditions are known at time :math:`t`, then integration
of :eq:`4.3-1` to time :math:`t + \Delta t` gives

.. math::
    :label: 4.6-1

	C_{\text{i}}\left( t + \Delta t \right) = C_{\text{i}}\left( t \right)e^{- \lambda_{\text{i}}\Delta t} + \frac{\beta_{\text{i}}}{\Lambda}e^{- \lambda_{\text{i}}\left( t + \Delta t \right)}\int_{t}^{t + \Delta t}{\phi\left( t' \right)e^{- \lambda_{\text{i}}t'}\text{dt}'}

The expressions given by :eq:`4.6-1` for the :math:`C_{\text{i}}` at the advanced
time are substituted directly into :eq:`4.2-4`. Next, the time variation
of the fission power amplitude over the time step is assumed to be

.. math::
    :label: 4.6-2

	\phi\left( t + \Delta t \right) = \phi\left( t \right) + \phi_{1}\Delta t + \phi_{2}\left( \Delta t \right)^{2}

where :math:`\phi_{1}` and :math:`\phi_{2}` are constants to be
determined. To find :math:`\phi_{1}` and :math:`\phi_{2}`, :eq:`4.6-2` is
first differentiated and substituted for in :eq:`4.2-4`. In addition,
:eq:`4.6-2` is substituted directly for in :eq:`4.2-4` and :eq:`4.6-1`. The integral
term in :eq:`4.6-1` is evaluated analytically. This yields a single
equation in the two unknowns :math:`\phi_{1}` and :math:`\phi_{2}`. To
obtain :math:`\phi_{1}` and :math:`\phi_{2}`, this equation is evaluated
for the full time step :math:`\Delta t`, and for the half time step
:math:`\Delta t/2`, to give two linear equations in the two unknowns.
Solution of this equation set yields :math:`\phi_{1}` and
:math:`\phi_{2}`, and thus an estimate of
:math:`\phi\left( t + \Delta t \right)`.

To control the precision of this approximate solution method, an
internal time-step control algorithm has been implemented. In addition
to the full time-step solution, a half time step solution is also
obtained (i.e. evaluation of :eq:`4.2-4` at :math:`\Delta t/2` and :math:`\Delta t/4`). The half-step
solution is extrapolated to the end of the time step and compared with
the full step solution. If the fractional difference is within a
user-specified tolerance, the full-step solution is accepted. If not,
the internal time step is halved and the entire solution process is
repeated. This procedure is carried out until the solution is advanced
to the end of the main time step. The half/full-step fission power
amplitude convergence precision is entered as :sasinp:`EPSPOW`.

On each main time step, the point kinetics solution is obtained twice.
At the beginning of the time step, the net reactivity is linearly
extrapolated in time to the end of the time step, and a solution for
:math:`\phi` is obtained with this extrapolated value. This solution is
then used to drive the energy equation solutions in the |SAS|
models through the time step. When all channels have been advanced to
the end of the main time step, the reactivities are calculated and the
net reactivity is used to solve once more for :math:`\phi`. This second
solution is taken as the final value for the time step, and solutions
for the delayed-neutron precursor and decay-heat equations are based on
this value.

For long, slow transients, an oscillation or instability in the solution
for reactivity and power has been observed when time steps larger than
0.25 seconds are used. Modifications to the code have been made to
eliminate this instability or oscillation. With the modified code, time
step sizes of 1 second can be used in the early part of the transient,
when flows and temperatures are changing rapidly, and 5 second step
sizes can be used after the first few hundred seconds of the transient.
The oscillations or instabilities come about because the reactivity used
to compute the power level for a time step is an extrapolated
reactivity, based on the last two previous computed steps. The
extrapolated reactivity is used to compute a power level that is then
used to compute fuel pin temperatures. After the temperatures are
computed, the reactivity is computed. The computed reactivity may not be
the same as the extrapolated reactivity. If the extrapolated reactivity
is high, it will lead to higher fuel temperatures, which will lead to a
lower computed reactivity. For the next step the extrapolated reactivity
will be low, and the computed reactivity will be high. For large time
steps, this leads to oscillations. The fuel expansion and Doppler
reactivity oscillate up-down-up-down from step to step.

Adding a correction term to the extrapolated reactivity,
:math:`k_{{\text{ext}}}^{n}`, for time step :math:`n` eliminated
the oscillations. First an error, :math:`e_{\text{n}}`, is calculated as the
difference between the computed reactivity, :math:`k^{n}`, and the
extrapolated reactivity:

.. math::
    :label: 4.6-2a

	e_{\text{n}} = k^{n} - k_{{\text{ext}}}^{n}

Then the reactivity :math:`k'` actually used for step :math:`n` is
computed as

.. math::
    :label: 4.6-2b

	k' = k_{{\text{ext}}}^{n} + \left| e_{\text{n} - 1} \right|{\text{sign}}\left( e_{\text{n} - 2} \right)

The correction term tends to cancel out the oscillations. The
application of this correction is controlled by input variable :sasinp:`JREEXT`.

.. _section-4.6.2:

Specified Power History
~~~~~~~~~~~~~~~~~~~~~~~

The point kinetics and decay-heat models may be overridden through the
specification on input of a user-supplied power history. This option is
intended to be used for analyses in which the power history is known,
such as those for in-pile experiments. It is also possible to specify an
external, driving reactivity that is summed with the internally
calculated feedback reactivity. This option may employ either a
user-specified subroutine function named PREA, or a user-specified table
of input power (or reactivity) values for input problem time values. The
choice of specifying power or external reactivity is controlled with
input variable :sasinp:`IPOWER`, and the choice of an
input function or table is set by :sasinp:`NPREAT`,
which gives the number of entries in the power (or reactivity) table and
defaults to zero, which indicates usage of the input function PREA.
Should a PREA not be supplied for this case, steady-state power or
reactivity conditions are assumed.

For :sasinp:`NPREAT` > 0, input values for the power (or reactivity) and problem
times at which they apply are entered in arrays :sasinp:`PREATB` and :sasinp:`PREATM`.

When the power vs. time option is invoked (:sasinp:`IPOWER` = 1), additional power
tables may be entered by setting :sasinp:`NPOWDK` to
the total number of tables to be input (2, …, 5), and entering the
additional power vs. time table in arrays :sasinp:`PRETB2` and :sasinp:`PRETM2`.
(The first table is always entered in :sasinp:`PREATB` and :sasinp:`PREATM`). Input integer
:sasinp:`IDKCRV` then specifies the input power curve to be used in a particular
channel. Total reactor power reported in the output represents only channels assigned to the first input power curve
and does not include power from channels assigned to input power curves 2-5.
Note that this multiple input power table option is only available for :sasinp:`IPOWER` = 1 and :sasinp:`NPOWDK` > 1.