.. _section-3.10:

Multiple-Pin Model
------------------

.. _section-3.10.1:

Introduction
~~~~~~~~~~~~

The new multiple pin treatment was added to the code to account for
pin-to-pin variations within a subassembly. This new multiple pin
treatment can be used in at least two different ways. One approach is to
use the new treatment to compute nominal or "best estimate" variations
within a subassembly. Another approach is to compute hot channel
behavior due to postulated deviations in coolant flow rate, coolant flow
area, and pin power. When doing a whole-core analysis, one probably
would not want the hot channel temperatures used in reactivity feedback
calculation or in the core outlet temperatures that feed into the
primary loop calculations. Therefore, one would probably use parallel
treatments for the same subassembly: a nominal one-pin or multi-pin
treatment used for reactivity feedback, plus a hot channel treatment
decoupled from the reactivity feedback and from the primary loop
calculation. The hot channel treatment would be de-coupled by setting
reactivity feedback coefficients to zero and by setting the number of
subassemblies represented by this treatment to zero.

The new multiple pin treatment allows a number of coupled channels to be
used to model a single subassembly. Thus a channel can represent a part
of a subassembly instead of the whole subassembly. Peaking factors can
be mechanistically calculated by reducing coolant flow areas and flow
rates or increasing pin power levels in some channels.

As previously mentioned, the multiple pin option is currently only
available for steady-state and single phase transient calculations.

.. _section-3.10.2:

Physical Model
~~~~~~~~~~~~~~

In the new multiple pin option, the regions above and below the pin
section of a subassembly are still represented by a single channel; but
a number of channels can be used to represent the pin section. Each
channel in the pin section can represent one or more concentric rows of
pins and their associated coolant. It is also possible for channels to
represent slices of pins for a subassembly with a strong lateral power
skew or for one with a hot subassembly on one side and a cool
subassembly on the opposite side. :numref:`figure-3.10-1` illustrates one way that
a number of channels can be used to model the pin section of a
subassembly as concentric rings of pins and coolant subchannels. The
thimble flow region in this figure is a feature of the EBR-II XX09
instrumented subassembly and is not found in ordinary subassemblies.
This case concentrates on the coolant subchannels and splits the pins.
It is also possible to shift the channel boundaries half a pin and use a
pin-based representation with intact pins and split coolant subchannels.
Currently up to 56 channels can be used to represent a subassembly.

The new option accounts for coolant-to-coolant heat transfer between
adjacent channels, including the effects of both conduction and
turbulent mixing. It also accounts for subassembly-to-subassembly heat
transfer from the duct wall of a subassembly, through the interstitial
sodium, to the duct wall of a neighboring subassembly. In addition,
axial conduction in the coolant is accounted for.

:numref:`figure-3.10-2` illustrates the axial representation of the subassembly
flow, with parallel flow paths through the pin section. Each channel
used to represent the pin section of subassembly has its own separate
time-dependent flow rate, and the flow rates in all channels in a
subassembly are driven by common pressures at the inlet and outlet of
the pin section. Thus, transient flow redistribution among channels is
accounted for. The single flow rate in the regions above and below the
pin section is the sum of the pin section flow rates, so the subassembly
flow orifice sees the correct total flow rate. Even though the
coolant-to-coolant heat transfer coefficient includes a term for
turbulent mixing between coolant subchannels, one effect that the flow
calculation does not account for is mass flow between channels in the
pin section, although cross-flow at the ends of the pin section is
allowed. Therefore, if recirculation loops occur within a subassembly at
low flows, the model would calculate them; but the recirculation loops
would go to the ends of the pin section.

.. _figure-3.10-1:

..  figure:: media/image16.png
	:align: center
	:figclass: align-center
	:width: 6.20764in
	:height: 5.60764in

	|SAS| Multiple Pin Representation and Thermocouple Locations for the EBR-II XX09 Instrumented Subassembly

.. _figure-3.10-2:

..  figure:: media/image17.png
	:align: center
	:figclass: align-center
	:width: 5.33056in
	:height: 7.44583in

	|SAS| Multiple Pin Treatment of a Subassembly

In the new model the coolant in channel I can transfer heat directly to
the coolant in channel :math:`\text{I}-1` and channel :math:`\text{I}+1`. Using correlations of the
same form as those used in THI3D code [3-6] and the HOTCHAN code [3-7],
the channel-to-channel heat flow per pin per unit height from channel I
to channel :math:`\text{I}+1` is calculated as:

.. math::
    :label: 3.10-1

    Q_{\text{I,I} + 1} = \left\lbrack U_{1}\overline{k} + U_{2} \overline{C}\left( w_{1} + w_{\text{I} + 1} \right) \right\rbrack \left( T_{\text{I}} - T_{\text{I} + 1} \right)

where

:math:`\overline{k}` = average thermal conductivity of the
coolant

:math:`\overline{C}` = average specific heat

:math:`w_{\text{i}}` = coolant mass flow per pin (kg/s) in channel :math:`I`,

and

:math:`T_{\text{I}}` = coolant temperature

In this equation, :math:`U_{\text{I}}` is a geometry factor for thermal
conduction, and :math:`U_2` is a product of a turbulent-mixing
coefficient and a geometry factor for turbulent mixing between channels.
Since a |SAS| channel usually models a group of coolant
subchannels, the values used for :math:`U_1` and :math:`U_2` must
account for a combination of parallel and series heat flow paths between
the middle of channel :math:`I` and the middle of channel :math:`I + 1`.
Subassembly-to-subassembly heat transfer is handled in a somewhat
simpler manner. For heat transfer from the outer surface of the
structure in channel :math:`I` to the outer surface of the structure in
channel :math:`J`, a constant value is used for the produce of the heat
transfer coefficient and the heat transfer area per unit height.

.. _section-3.10.3:

Numerical Methods
~~~~~~~~~~~~~~~~~

Most of the transient heat transfer calculations and flow rate
calculations in SASSYS-1 use semi-implicit time differencing in order to
obtain stable and accurate solutions with reasonably long time steps.
Before the onset of coolant boiling or pin disruption, time step sizes
of a second or more are commonly used; and the code usually runs
significantly faster than real time on a Cray XMP computer.

From a numerical computation point of view, the two main tasks in adding
the multiple pin model to the code were the coolant-to-coolant heat
transfer calculation and the coolant flow rate calculation with parallel
flow paths in the pin section. Both of these calculations use
semi-implicit time differencing. In the single-pin model, coolant
temperatures for all of the radial temperature nodes in the pin,
coolant, and structure at one axial node are solved for simultaneously
in order to obtain a semi-implicit time differencing solution without
iteration. In the new multiple-pin treatment, this concept is carried
one step further. At a given axial node, temperature at all of the
radial nodes for all channels representing a subassembly are solved for
simultaneously. In the heat transfer calculations, the axial conduction
terms, which are small, are treated with explicit forward time
differencing so that axial nodes are decoupled and can be treated
separately except for the coolant convection terms. The axial coupling
due to the coolant convection terms is handled by starting at one end of
the subassembly and solving for axial nodes one at a time in the
direction of the flow. If flow has reversed in some channels but not in
others, the calculation progresses in the direction of the dominate
flow; and explicit forward differencing is used for the coolant
convection terms in the non-dominate flow direction channels. The
subassembly-to-subassembly heat fluxes are calculated with explicit
forward differencing in time, and this does impose a stability limit on
the time step size. For typical subassembly duct wall thicknesses, the
explicit subassembly-to-subassembly heat flux calculation limits the
maximum time step size to a value in the range from .25 to .5 seconds.
For the coolant flow rate calculations, the incompressible flow momentum
equations are linearized about conditions at the beginning of the time
step. Then, flow at the end of the step are calculated for all channels
in a subassembly simultaneously.

A null transient is used to obtain steady-state temperatures at the
start of the regular transient. First all coolant, pin, structure, and
reflector temperatures in all subassemblies are set to the coolant inlet
temperature. Then, the power levels and coolant flow rates are held
constant while a number of transient heat transfer time steps are made.
Since the pin thermal time constant and the coolant transit time through
a subassembly are both less than a second, the null transient results
converge rapidly if reasonably large time steps are used.

.. _section-3.10.4:

Detailed Mathematical Treatment
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The main computational parts of the new multiple pin model are the heat
transfer calculations in the pin section of a subassembly, the coolant
flow rate calculations for a subassembly, the subassembly-to-subassembly
heat transfer, and changes to the driver routines to call the
appropriate new routines at the proper times. The pin section heat
transfer calculations are done in two new subroutines: TSHTM3 calculates
temperatures in the core and axial blankets, and TSHTM2 calculates
temperatures in the gas plenum region. The existing TSHTN1 still
calculates temperatures in the reflector zones above and below the pin
section, but it was modified. TSHTN1 had to be modified slightly to pick
up mixed mean outlet temperatures from the multiple channels
representing the pin section. Also, axial conduction in the coolant was
added to TSHTN1. The new TSCLM1 routine calculates transient coolant
flow rates for a subassembly. The subassembly-to-subassembly heat fluxes
are calculated in CHCHFL.

.. _section-3.10.4.1:

Heat Transfer Calculations in the Core and Axial Blankets: Subroutine TSHTM3
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The multiple pin heat transfer calculations for the core and axial
blankets in subroutine TSHTM3 are similar to the single pin calculations
in the existing subroutine TSHTN3, as described in :numref:`section-3.3.1`. One
difference is that TSHTN3 does one time step for one channel each time
it is called, whereas TSHTM3 solves for temperatures in all pins or
channels representing a subassembly when it is called. Also, TSHTM3 adds
extra terms for coolant-to-coolant heat transfer; and TSHTM3 includes
axial conduction in the coolant.

Within a fuel pin, one-dimensional radial heat transfer is used. The
basic heat transfer equation is :eq:`3.3-1`. This part of the calculation
in TSHTM3 is identical to the corresponding calculation in TSHTN3.

For the coolant in channel I, the heat transfer calculations are
basically one-dimensional in the axial direction, with extra terms for
coolant-to-coolant heat transfer from channel :math:`I-1` and :math:`I+1`. The basic
heat transfer equation is:

.. math::
    :label: 3.10-2

    \rho c A_{\text{c}}\frac{\partial \text{T}}{\partial \text{t}} + \frac{\partial}{\partial \text{z}}\left( wcT \right) = Q_{\text{T}}A_{\text{c}}

where

:math:`A_{\text{c}}` = coolant flow area

:math:`z` = axial position

:math:`w` = coolant mass flow rate

and

:math:`Q_{\text{r}}` = total heat source per unit volume

The heat source contains a number of terms:

.. math::
    :label: 3.10-3

    Q_{\text{T}} = Q_{\text{c}} + Q_{\text{ec}} + Q_{\text{sc}} + Q_{\text{ax}} + {Q'}_{\text{I} - 1,\text{I}} - {Q'}_{\text{I,I} + 1}

where

:math:`Q_{\text{c}}` = source due to direct heating of the coolant by neutrons
and gamma rays,

:math:`Q_{\text{ec}}` = heat flow from cladding to coolant,

:math:`Q_{\text{sc}}` = heat flow from structure to coolant,

:math:`Q_{\text{ax}}` = axial conduction source, given by

.. math::
    :label: 3.10-4

    Q_{\text{ax}} = \frac{\partial}{\partial \text{z}} k \frac{\partial \text{T}}{\partial \text{z}}

and

.. math::
    :label: 3.10-5

    {Q'}_{\text{I,I} + 1} = \frac{Q_{\text{I,I} + 1}}{A_{\text{c}}}

where :math:`{Q'}_{\text{I,I} + 1}` is given by equation 3.10-1. Note that
:eq:`3.10-2` is the same as :eq:`3.3-5`, except that additional source terms are
included in the multiple pin model.

Finite differencing in space and time is used to solve these heat
transfer equations. :numref:`figure-3.2-3` shows the axial and radial mesh used
for a single channel. In the multiple pin model the pin section mesh is
repeated for each channel, but the axial reflector zones are only used
in the first channel used to represent a subassembly. At each axial node
in the core and axial blankets a number of radial nodes are used for the
fuel, three radial nodes are used for the cladding, one node is used in
the coolant, and two nodes are used in the structure. All channels,
representing a subassembly must use the same axial mesh. Also, all
subassemblies connected by subassembly-to-subassembly heat transfer must
use the same axial mesh.

For a given time step, :eq:`3.3-1` and :eq:`3.10-2` become linear finite
difference equations. The temperatures are known at the beginning of the
time step, and the temperature at the end of the step are the unknowns
to be solved for. Semi-implicit time differencing is used for all terms
except the axial conduction terms :math:`Q_{\text{ax}}`, and the
coolant-to-coolant terms :math:`Q_{\text{I} - 1,\text{I}}` and :math:`{Q'}_{\text{I,I} + 1}`.
The :math:`Q_{\text{ax}}` axial conduction term is calculated explicitly based
on temperatures at the beginning of the time step. The
coolant-to-coolant terms are calculated fully implicitly, using the
coolant flow rates and temperatures at the end of the step.

After finite differencing for one axial node, one obtains N simultaneous
liner equations in N unknowns, where the unknowns are the temperatures
at the end of the time step for all radial nodes in all channels
representing the subassembly. These equations are solved by Gaussian
elimination. The equations are basically tri-diagonal, with extra
non-tri-diagonal terms for coolant-to-coolant heat transfer between
channels. No full N by N matrix is ever set up by TSHTM3, and a general
full matrix solution package is not used. Instead, only non-zero terms
are computed and stored; and the Gaussian elimination solution was
written specifically for this set of equations. The result is that the
storage requirements vary linearly, rather than quadratically, with the
maximum allowable value for N; and the computation time varies linearly,
rather than quadratically, with the actual value of N.

.. _section-3.10.4.2:

Heat Transfer Calculations in the Gas Plenum Region: Subroutine TSHTM2
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Subroutine TSHTM2 does the heat transfer calculations for one time step
for all axial nodes in the gas plenum region for all channels
representing a subassembly when it is called. TSHTM2 is similar to
TSHTM3, but TSHTM2 deals with fewer radial nodes. As shown in :numref:`figure-3.2-5`, in the plenum region there are still two radial nodes in the
structure and one in the coolant; but there is only one in the cladding;
and there are no fuel nodes. Instead of fuel temperatures there is one
gas temperature common to all of the axial nodes in a pin. As in TSHTM3,
for a given axial node TSHTM2 solves simultaneously for temperatures at
all radial modes in all channels representing a subassembly.

One special problem in TSHTM2 is the gas temperature, which is common to
all axial nodes in the gas plenum region of a pin. Axial nodes are
handled one at a time, rather than simultaneously, whereas a
semi-implicit time differencing treatment would require solving for all
axial nodes simultaneously. The new multiple-pin TSHTM2 routine handles
this problem in a simpler manner than the old single-pin TSHTN2. The way
that the gas temperature is handled is TSHTM2 is to calculate a separate
gas temperature for each axial node. At the beginning of the time step,
the gas temperatures at all axial nodes are set to one common value.
Then separate values are calculated for each axial node at the end of
the step. Finally, an average of the separate values is calculated for
the one common value at the ed of the step.

.. _section-3.10.4.3:

Coolant Flow Rates: Subassembly TSCLM1
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The coolant for rate calculation for a subassembly in the multiple-pin
routine TSCLM1 is much more complex than the corresponding calculation
in the single pin routine TSCNV1. As shown in :numref:`figure-3.10-2`, the
multiple pin calculation involves multiple parallel flow paths in the
pin section, in series with single flow paths in the reflector regions.
Therefore, the multiple-pin routine TSCLM1 was written from scratch,
rather than starting form the single pin TSCNV1 routine.

:numref:`figure-3.10-3` illustrates the main variables used in TSCLM1.

Incompressible flow is used for these calculations, so conservation of
mass gives

.. math::
    :label: 3.10-6

    N_{1}\ w_{\text{r}} = \sum_{\text{k}}{w_{\text{pk}} N_{\text{k}}}

where

:math:`w_{\text{r}}` = coolant mass flow rate per pin in the reflector zones,

:math:`w_{\text{pk}}` = pin section coolant mass flow rate per pin in channel
k,

:math:`N_{\text{k}}` = number of pins in channel k,

and

:math:`I` = channel number of the first channel representing the subassembly.

The momentum equation for the pin section is

.. math::
    :label: 3.10-7

    \frac{L_{\text{p}}}{A_{\text{cpk}}} \frac{d\text{w}_{\text{pk}}}{\text{dt}} = p_{\text{b}} - p_{\text{c}} - \Delta p_{\text{pk}}

.. _figure-3.10-3:

..  figure:: media/image18.png
	:align: center
	:figclass: align-center
	:width: 6.25000in
	:height: 6.98125in

	Subassembly Coolant Flows

where

:math:`L_{\text{p}}` = length of pins,

:math:`A_{\text{cpk}}` = coolant flow area in channel k in the pin section,

:math:`t` = time,

:math:`p_{\text{b}}` = coolant pressure at the bottom of the pins,

:math:`p_{\text{c}}` = coolant pressure at the top of the pins,

.. math::
    :label: 3.10-8

    \Delta p_{\text{pk}} = \text{ pressure loss } = \sum_{\text{jv} = \text{JCPNBT}}^{\text{JCPNTM}}{\Delta p_{\text{kjc}}}

:math:`jc` = axial coolant node,

:math:`\text{JCPNBT}` = first axial node in the pin section

:math:`\text{JCPNTM}` = last axial node in the pin section,

and

:math:`\Delta p_{\text{kjc}}` = pressure loss in axial node jc of
channel :math:`k`.

The pressure loss contains a number of terms:

.. math::
    :label: 3.10-9

    \Delta p_{\text{kjc}} = \Delta p_{\text{frkjc}} + \Delta p_{\text{grkjc}} + \Delta p_{\text{orkjc}} + \Delta p_{\text{acckjc}}

where

.. math::
    :label: 3.10-10

    \Delta p_{\text{frkjc}} = \text{ friction loss } = \frac{f_{\text{kjc}}\Delta z_{\text{jc}}w_{\text{pk}}\left| w_{\text{pk}} \right|}{2{\overline{\rho}}_{\text{ckjc}}\ A_{\text{cpk}}\ D_{\text{hpk}}}

.. math::
    :label: 3.10-11

    \Delta p_{\text{grpjc}} = \text{ gravity head } = g {\overline{p}}_{\text{ckjc}} \Delta z_{\text{jc}}

.. math::
    :label: 3.10-12

    \Delta p_{\text{orfkjc}} = \text{ orifice loss } = \frac{K_{\text{orkjc}} w_{\text{pk}} \left| w_{\text{ph}} \right|}{2{\overline{\rho}}_{\text{ckjc}}A_{\text{cpk}}^{2}}

.. math::
    :label: 3.10-13

    \Delta p_{\text{acckjc}} = \text{ acceleration term } = \frac{w_{\text{pk}}^{2}}{A_{\text{cpm}}^{2}}\left( \frac{1}{\rho_{\text{ckjc} + 1}} - \frac{1}{\rho_{\text{ckjc}}} \right)

:math:`\Delta z_{jc}` = axial node length,

:math:`f_{\text{kjc}}` = friction factor in channel k at node :math:`jc`,

:math:`\rho_{\text{ckjc}}` = coolant density at bottom of node :math:`jc` in
channel :math:`k`,

:math:`{\overline{\rho}}_{\text{ckjc}}` = average coolant
density in node :math:`jc` in channel :math:`k`,

:math:`D_{\text{hpk}}` = Hydraulic diameter,

and

:math:`K_{\text{orkjc}}` = orifice coefficient in node :math:`jc` of channel
:math:`k`

The friction factor is calculated as in :eq:`eq-3.9-3a`.

using

.. math::
    :label: 3.10-14

    \text{Re}= \text{ Reynolds number } = \frac{D_{\text{hpk}}\left| w_{\text{pk}} \right|}{\mu_{\text{kjc}}\ A_{\text{cpk}}}

and

:math:`\mu_{\text{kjc}}` = viscosity of the coolant.

After semi-implicit finite differencing in time and linearization of the
pressure drop terms, applying :eq:`3.10-7` to a time step from :math:`t` to :math:`t + \Delta t` gives

.. math::
    :label: 3.10-15

    \frac{L_{\text{p}}}{A_{\text{cpk}}} \frac{\Delta w_{\text{pk}}}{ \Delta t } = p_{\text{b}}\left( t \right) + \theta_{2}\Delta p_{\text{b}} - p_{\text{c}}\left( t \right) - \theta_{2}\Delta p_{\text{c}} - \Delta p_{\text{pk}}\left( t \right) - \theta_{2} \Delta w_{\text{pk}}K_{\text{wpk}}

where

.. math::
    :label: 3.10-16

    \Delta w_{\text{pk}} = w_{\text{pk}}\left( t + \Delta t \right) - w_{\text{pk}}\left( t \right)

.. math::
    :label: 3.10-17

    \Delta p_{\text{b}} = p_{\text{b}}\left( t + \Delta t \right) - p_{\text{b}}\left( t \right)

.. math::
    :label: 3.10-18

    \Delta p_{\text{c}} = p_{\text{c}}\left( t + \Delta t \right) - p_{\text{c}}\left( t \right)

:math:`\theta_{2}` = degree of implicitness

.. math::
    :label: 3.10-19

    \theta_{1} = 1 - \theta_{2}

.. math::
    :label: 3.10-20

    K_{\text{wpk}} = \sum_{\text{jc} = \text{JCPNBT}}^{\text{JCPNTM}}{K_{\text{wpkjc}}}

.. math::
    :label: 3.10-21

    K_{\text{wpkjc}} = \frac{\partial\Delta \text{p}_{\text{kjc}}}{\partial \text{w}_{\text{pk}}} = K_{\text{frkjc}} + K_{\text{wcrkjc}} + K_{\text{acckjc}}

.. math::
    :label: 3.10-22a

    K_{\text{frkjc}} = \begin{cases}
    \left( 2 + b_{\text{fr}} \right)\frac{\Delta p_{\text{frkjc}}\left( t \right)}{w_{\text{pk}}\left( t \right)} & \text{if }R_{\text{ckjc}}\  \geq \ R_{\text{et}} \\
    \frac{\Delta p_{\text{frkjc}}\left( t \right)}{w_{\text{pk}}\left( t \right)} & \text{if }R_{\text{ckjc}}\  < \ R_{\text{et}} \\
    \end{cases}

or

.. math::
    :label: 3.10-22b

    K_{\text{frkjc}} = \left\lbrack \left( 2 + b_{\text{fr}} \right)A_{\text{fr}} A \left( \text{RE} \right)^{b_{\text{fr}}} + \frac{A_{\text{fL}}}{\text{RE}} \right\rbrack\frac{\Delta z_{\text{jc}}\left| w_{\text{pk}}\left( t \right) \right|}{2p_{\text{ckjc}}D_{\text{hpk}}A_{\text{cpk}}}

.. math::
    :label: 3.10-23

    K_{\text{wpkjc}} = \frac{2\Delta p_{\text{orfkjc}}\left( t \right)}{w_{\text{pk}}\left( t \right)}

and

.. math::
    :label: 3.10-24

    k_{\text{acckjc}} = \frac{2\Delta p_{\text{acckjc}}\left( t \right)}{w_{\text{pk}}\left( t \right)}

:eq:`3.10-15` can then be written as

.. math::
    :label: 3.10-25

    \Delta p_{\text{b}} - \Delta p_{\text{c}} = d_{0\text{pk}} + d_{0\text{pk}} + d_{1\text{pk}}\Delta w_{\text{pk}}

where

.. math::
    :label: 3.10-26

    d_{0\text{pk}} = \frac{p_{\text{c}}\left( t \right) - p_{\text{b}}\left( t \right) + \Delta p_{\text{pk}}\left( t \right)}{\theta_{2}}

and

.. math::
    :label: 3.10-27

    d_{1\text{pk}} = \frac{\frac{L_{\text{p}}}{A_{\text{cpk}}} + \theta_{2} \Delta t K_{\text{wpk}}}{\theta_{2} \Delta t }

Applying a similar process to the upper and lower reflector regions
gives:

.. math::
    :label: 3.10-28

    \Delta p_{\text{b}} = D_{0\text{l}} + D_{1\text{l}} \sum_{\text{k}} \Delta w_{\text{pk}}

and

.. math::
    :label: 3.10-29

    \Delta p_{\text{c}} = D_{0\text{u}} + D_{1\text{u}} \sum_{\text{k}} \Delta w_{\text{pk}}

where

.. math::
    :label: 3.10-30

    D_{0\text{l}} = \frac{\theta_{1}p_{\text{a}}\left( t \right) + \theta_{2} p_{\text{a}}\left( t + \Delta t \right) - p_{\text{b}}\left( t \right) - \Delta p_{\text{l}}\left( t \right)}{\theta_{2}}

.. math::
    :label: 3.10-31

    D_{1\text{l}} = \frac{- \sum_{\text{KA} = 1}^{\text{KZPIN} - 1}{\frac{L_{\text{kz}}}{A_{\text{cKZ}}}} - \theta_{2} \Delta t K_{w\text{l}}}{\theta_{2} \Delta t }

.. math::
    :label: 3.10-32

    D_{0\text{u}} = \frac{\theta_{1}p_{\text{d}}\left( t \right) + \theta_{2}\ p_{\text{d}}\left( t + \Delta t \right) - p_{\text{c}}\left( t \right) + \Delta p_{\text{u}}\left( t \right)}{\theta_{2}}

.. math::
    :label: 3.10-33

    D_{1\text{u}} = \frac{\sum_{\text{KZ} - \text{KZPIN} + 1}^{\text{KZM}}{\frac{\text{L}_{\text{KZ}}}{A_{\text{cKZ}}}} + \theta_{2} \Delta t K_{\text{wu}}}{\theta_{2} \Delta t }

:math:`L_{\text{KZ}}` = length of zone :math:`kz`

:math:`A_{\text{cKZ}}` = coolant flow area in zone :math:`kz`

:math:`\text{KZPIN}` = axial zone number of the pin section

:math:`\text{KZM}` = last zone number

.. math::
    :label: 3.10-34

    \Delta p_{\text{l}} = \sum_{\text{jc} = 1}^{\text{JCPNBT} - 1}{ \Delta p_{\text{rjc}}}

.. math::
    :label: 3.10-35

    \Delta p_{\text{u}} = \sum_{\text{jc} = \text{JCPNTP}}^{\text{MZCM}1} \Delta p_{\text{rjc}}

:math:`\Delta p_{\text{rjc}}` = pressure loss in node :math:`jc` of a reflector zone

:math:`\text{JCPNBT}` = first axial node in pin section

:math:`\text{JCPNTP}` = first axial node above the pin section

:math:`\text{MZCM}1` = last axial node

.. math::
    :label: 3.10-36

    K_{\text{wl}} = \frac{\partial\Delta p_{\text{l}}}{\partial w_{\text{r}}}

and

.. math::
    :label: 3.10-37

    K_{\text{wu}} = \frac{\partial\Delta p_{\text{u}}}{\partial w_{\text{u}}}

Then, combining :eq:`3.10-25` and :eq:`3.10-28` gives

.. math::
    :label: 3.10-38

    \Delta p_{\text{b}} = B_{\text{o}} + B_{1} \Delta p_{\text{c}}

where

.. math::
    :label: 3.10-39

    B_{\text{o}} = \frac{D_{\text{ol}} -  D_{1\text{l}} S_{1}}{1 -  D_{1\text{f}} S_{\text{o}}}~,

.. math::
    :label: 3.10-40

    B_{1} = \frac{- D_{1\text{l}}S_{\text{o}}}{1 - D_{1\text{l}}S_{\text{o}}}

.. math::
    :label: 3.10-41

    S_{\text{o}} = \sum_{\text{k}}{ \frac{1}{d_{1\text{pk}}}}~ ,

and

.. math::
    :label: 3.10-42

    S_{1} = \sum_{\text{k}}{ \frac{d_{\text{opk}}}{d_{1\text{pk}}}}~ ,

Similarly, combining :eq:`3.10-25` and :eq:`3.10-29` gives

.. math::
    :label: 3.10-43

    \Delta p_{\text{c}} = C_{\text{o}} + C_{1} \Delta p_{\text{b}}

where

.. math::
    :label: 3.10-44

    C_{\text{o}} = \frac{D_{\text{ou}} - D_{1\text{u}} S_{1}}{1 +  D_{1\text{u}} S_{\text{o}}}~ ,

and

.. math::
    :label: 3.10-45

    C_{1} = \frac{D_{1\text{u}} S_{\text{o}}}{1 +  D_{1\text{u}} S_{\text{o}}}~ .

Combining :eq:`3.10-38` and :eq:`3.10-43` gives

.. math::
    :label: 3.10-46

    \Delta p_{\text{c}} = \frac{C_{\text{o}} + C_{1}B_{\text{o}}}{1 -  C_{1} B_{1}}

Then :math:`\Delta p_{\text{b}}` can be obtained from :eq:`3.10-38` and
:math:`\Delta w_{\text{bk}}` can be obtained from :eq:`3.10-25` for each channel.

After :math:`p_{\text{b}}`, :math:`p_{\text{c}}`, and the coolant flow rates have been
calculated for the end of the time step it is possible to calculate the
pressures. The pressure is calculated at the axial node boundaries.
:math:`p_{\text{kjc}}` is the pressure at the bottom (inlet end) of node :math:`jc`
in channel :math:`k`. First, the nodal pressure loss at :math:`t + \Delta t` is
calculated as

.. math::
    :label: 3.10-47

    \Delta p_{\text{pkjc}}\left( t +  \Delta t  \right) = \Delta p_{\text{kjc}}\left( t \right) + K_{\text{wpkjc}} \Delta w_{\text{pk}}

and :math:`\Delta p_{\text{pk}} (t + \Delta t)` is calculated from :eq:`3.10-8`.
Then :eq:`3.10-7` is used to obtain
:math:`\frac{\text{dw}_{\text{pk}}\left( t + \Delta t \right)}{\text{dt}}`.

Integrating the momentum equation over one axial node gives

.. math::
    :label: 3.10-48

    \frac{\Delta z_{\text{j}}}{A_{\text{cpk}}} \frac{\text{dw}_{\text{pk}}}{\text{dt}} = p_{\text{pkjc}} -  p_{\text{pkjc} + 1} - \Delta p_{\text{kjc}}

Starting by setting

.. math::
    :label: 3.10-49

    p_{\text{pkjJCPNTP}} = P_{\text{c}}

the code marches down the channel, using :eq:`3.10-48` to obtain
:math:`p_{\text{pkjc}}` after :math:`p_{\text{pkjc}+1}` has been calculated.
A similar procedure is used for calculating the pressures in the upper
and lower reflector zones.

The equation used to compute the degree of implicitness as a function of
time step size in TSCMV1 is

.. math::
    :label: 3.10-50

    \theta_{2} = \frac{a + bx + x^{2}}{2a + cx + x^{2}}

where

.. math::
    :label: 3.10-51

    x=\Delta t/\tau

:math:`\tau` = a time constant

:math:`a` = 6.12992

:math:`b` = 2.66054

and

:math:`c` = 3.56284

The basis for this expression is given in :numref:`section-A3.1`. For a single
channel treatment, the time constant, :math:`\tau` would be

.. math::
    :label: 3.10-52

    \tau = \frac{\sum_{\text{KZ} = 1}^{\text{KZM}}{ \frac{L_{\text{KZ}}}{A_{\text{cKZ}}} + \left( \frac{L_{\text{i}}}{A} \right)_{\text{i}} +  \left( \frac{L_{\text{i}}}{A} \right)_{\text{x}}}}{K_{\text{wpk}}\ K_{\text{w}1} + K_{\text{wu}}}

where

:math:`\left( \frac{L_{\text{i}}}{A} \right)_{\text{i}}` = extra coolant inertia term
at the subassembly inlet to account for inertia of the coolant in the
inlet plenum,

and

:math:`\left( \frac{L_{\text{i}}}{A} \right)_{\text{x}}` = same for the subassembly
outlet.

Since a simultaneous solution of flows in all channels of an assembly
plus the lower and upper reflector zones is required, the overall time
constant is calculated as

.. math::
    :label: 3.10-53

    \tau = \frac{ \frac{L_{\text{p}}}{\sum_{\text{k}}{A_{\text{cpk}}}} + \sum_{\text{kZ} \neq \text{KZPIN}}{\frac{L_{\text{KZ}}}{A_{\text{cKZ}}} + \left( \frac{L_{\text{i}}}{A} \right)_{\text{i}} +  \left( \frac{L_{\text{i}}}{A} \right)_{\text{x}}} }{K_{\text{w}1} +  K_{\text{wu}} +  \frac{1}{\sum_{k}{\frac{1}{k_{\text{wpk}}}}}}

.. _section-3.10.5:

Relationship Between Single Pin and Multiple Pin Models
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The single pin treatment and the multiple pin treatment coexist in the
|SAS| code; it is possible to use single pin treatments for
some subassemblies and multiple pin treatments for other subassemblies
in the same problem. There are now two separate sets of subroutines in
the code for computing single phase thermal hydraulics: the single pin
routines and multiple pin routines. Since the multiple pin model can
handle the case of a single pin per subassembly, the single pin single
phase thermal hydraulics subroutines are now redundant; and at some
point in the future they will probably be removed from the code.

As previously mentioned, currently the multiple pin model is only
available for steady state and single phase transient calculations; no
multiple pin treatment is available in the code for coolant boiling, pin
disruption, or relocation of fuel or cladding. In the future, the
multiple pin treatment will probably be extended to coolant boiling and
to relocation of fuel and cladding.