4.6. Solution Methods

The SAS4A/SASSYS‑1 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 SAS4A/SASSYS‑1 neutronics solution methods have been adopted from SAS3A.

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 SAS4A/SASSYS‑1, this time interval is taken as the main time step. If conditions are known at time \(t\), then integration of Eq. 4.3-1 to time \(t + \Delta t\) gives

(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 \(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

(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 \(\phi_{1}\) and \(\phi_{2}\) are constants to be determined. To find \(\phi_{1}\) and \(\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 Eqs. 4.2-4 and 4.6-1. The integral term in Eq. 4.6-1 is evaluated analytically. This yields a single equation in the two unknowns \(\phi_{1}\) and \(\phi_{2}\). To obtain \(\phi_{1}\) and \(\phi_{2}\), this equation is evaluated for the full time step \(\Delta t\), and for the half time step \(\Delta t/2\), to give two linear equations in the two unknowns. Solution of this equation set yields \(\phi_{1}\) and \(\phi_{2}\), and thus an estimate of \(\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 \(\Delta t/2\) and \(\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 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 \(\phi\) is obtained with this extrapolated value. This solution is then used to drive the energy equation solutions in the SAS4A/SASSYS‑1 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 \(\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, \(k_{{\text{ext}}}^{n}\), for time step \(n\) eliminated the oscillations. First an error, \(e_{\text{n}}\), is calculated as the difference between the computed reactivity, \(k^{n}\), and the extrapolated reactivity:

(4.6‑2a)

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

Then the reactivity \(k'\) actually used for step \(n\) is computed as

(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 JREEXT.

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 IPOWER, and the choice of an input function or table is set by 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 NPREAT > 0, input values for the power (or reactivity) and problem times at which they apply are entered in arrays PREATB and PREATM.

When the power vs. time option is invoked (IPOWER = 1), additional power tables may be entered by setting NPOWDK to the total number of tables to be input (2, …, 5), and entering the additional power vs. time table in arrays PRETB2 and PRETM2. (The first table is always entered in PREATB and PREATM). Input integer IDKCRV then specifies the input power curve to be used in a particular channel. Note that this multiple input power table option is only available for IPOWER = 1 and NPOWDK > 1.