.. _section-15.2:
In-pin Hydrodynamic Model
-------------------------
.. _section-15.2.1:
Overview of the Numerical Approach for the In-pin Fuel Motion Calculation and Description of Subroutines PN1PIN and PN2PIN
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
There are several requirements for the solution algorithm for this
one-dimensional, compressible flow problem with variable flow cross
section: (a) it has to be able to handle large mass sources (due to fuel
melting at the cavity boundary), (b) it has to conserve mass perfectly,
and (c) it has to run efficiently. This was achieved with a
predominantly explicit Eulerian solution method. All convective mass,
energy, and momentum fluxes are treated explicitly (i.e., the beginning
of time-step values are used). However, the solution sequence of the
different equations introduces a certain implicitness, which can be
important when treating the strong mass sources.
The set of conservation equations describing the in-pin fuel motion
includes the following; 3 for mass, 1 for momentum, and 1 for energy.
The mass conservation equations are for the fuel, and the free and
dissolved fission gas. The fuel and free fission-gas mass conservation
equations are solved first, followed by the fuel energy conservation
equation. The fission gas temperature change is assumed to be the same
as the fuel temperature change. From the results of the mass and energy
conservation equations, a new pressure is calculated. This is not the
true end-of-time-step pressure because the velocity changes during the
time step have not been included. However, it is a proper prediction for
the end-of-time-step value in an explicit sense. The new pressures are
then used in the fuel/fission-gas momentum equation. There is an input
option in PINACLE to use a combination of the new and the old pressures.
This option lets the pressure, P, be:
.. math::
:label: 15.2-1
P = \left( 1 - \text{EPCH} \right) \cdot \text{beginning-of-time-step pressure} \\
+ \text{EPCH} \cdot \text{advanced pressure.}
The recommended input value is :math:`\text{EPCH} = 1` because the calculation
remained stable for longer time steps in test problems involving shock
propagation and shock reflection when this input value was used [15-5].
The time-step criterion in this compressible calculation is the sonic
Courant condition. This does not require particularly small time steps
because the sonic velocity in two-phase mixtures is fairly low for the
void fractions encountered in pin blowdown calculations.
A staggered numerical grid is used with the densities, energies, and
temperatures defined at the cell centers, and the velocities defined at
the cell edges. The spatial differencing uses full donor cell
differencing. Although this is not as accurate as higher order
differencing, it makes the calculation very stable.
The free fission-gas mass conservation, the fuel mass conservation, and
the fuel energy conservation equations are solved in subroutine PN1PIN
(**P**\ I\ **N**\ ACLE **1**\ st **PIN** ROUTINE). PN1PIN also computes the
molten cavity geometry changes due to fuel melt-in. The fuel/fission-gas
momentum equation and the dissolved fission-gas mass equation are solved
in PN2PIN (**P**\ I\ **N**\ ACLE **2**\ nd **PIN** ROUTINE). This routine also
calculates the sonic velocities for each node and the maximum
hydrodynamic time step (see :numref:`section-15.3.1.3`).
.. _section-15.2.2:
Definition of the Generalized Smear Densities for the In-pin Hydrodynamic Calculation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The use of generalized smear densities in SAS4A for the PLUTO2 [15-3]
and LEVITATE [15-2] fuel motion models was prompted by the many
different moving and stationary components treated by these models.
PINACLE, which has used the in-pin fuel motion model in PLUTO2 and
LEVITATE as a starting point, has maintained the use of the smeared
densities for consistency and convenience. The use of this concept also
simplified the differential and finite difference equations for variable
cross section flow. The pie chart in :numref:`figure-15.2-1` gives an example of the
relative cross sectional areas within a subassembly or experimental loop
at a certain axial elevation
If the total area of the subassembly is AXMX, the generalized volume
fraction of a certain component :math:`k` is:
.. math::
:label: 15.2-2
\theta_{k} \left( z,t \right) = \frac{A_{\text{k}}\left( z,t \right)}{\text{AXMX}}
where :math:`A_{\text{k}}` is the cross sectional area occupied by component :math:`k`
(the latter refers, e.g., to the cavities or the moving fuel in all
failed pins within the subassembly cross section). The reference area
AXMX, which is an input quantity, can be arbitrarily chosen. This is
because the code is invariant to the choice of AXMX (i.e., as long as it
is not varied by several orders of magnitude which can lead to
differences due to truncation errors). However, the recommended value of
AXMX is the cross sectional area of a subassembly or experiment test
section (including the can wall) because the volume fractions of the
different components that appear in the PINACLE output are better
understood in this case.
The generalized volume fraction of the free fission gas and fuel vapor
in the cavity is the difference between the cavity volume fraction and
the fuel volume fraction.
.. math::
:label: 15.2-3
\theta_{\text{fica}}\left( z,t \right) = \theta_{\text{ca}}\left( z,t \right) - \theta_{\text{fuca}}\left( z,t \right)
where:
:math:`\theta_{\text{fica}}` = Generalized volume fraction of the free fission gas
and fuel vapor in the pin cavities
:math:`\theta_{\text{ca}}` = Generalized volume fraction of the molten cavities in
all pins which can be calculated from :math:`A_{\text{ca}}/\text{AXMX}`
:math:`\theta_{\text{fuca}}` = Generalized volume fraction of the fuel in the
cavities of all pins
Generalized smear densities, which are always marked by a prime, are
defined as products of physical densities and generalized volume
fractions:
.. math::
:label: 15.2-4
{\rho'}_{\text{fuca}}\left( z,t \right) = \rho_{\text{fuca}}\left( T \right) \cdot \theta_{\text{fuca}}\left( z,t \right) \
= \rho_{\text{fuca}}\left( T \right) \cdot \frac{A_{\text{fuca}}\left( z,t \right)}{\text{AXMX}}
.. math::
:label: 15.2-5
{\rho'}_{\text{fica}}\left( z,t \right) = \rho_{\text{fica}}\left( T \right) \cdot \theta_{\text{fica}}\left( z,t \right) \
= \rho_{\text{fica}}\left( T \right) \cdot \frac{A_{\text{fica}}\left( z,t \right)}{\text{AXMX}}
.. math::
:label: 15.2-6
{\rho'}_{\text{fsca}}\left( z,t \right) = \rho_{\text{fsca}}\left( T \right) \cdot \theta_{\text{fuca}}\left( z,t \right) \
= \rho_{\text{fsca}}\left( T \right) \cdot \frac{A_{\text{fuca}}\left( z,t \right)}{\text{AXMX}}
.. _figure-15.2-1:
.. figure:: media/image7.png
:align: center
:figclass: align-center
Illustration of the Generalized Volume Fraction
where the subscript **fsca** refers to the fission gas which is dissolved
in the cavity fuel. The temperature :math:`T` is the fuel temperature which
should actually be written with subscript **fuca**. It should again be
pointed out that the :math:`A`'s refer to total cross section areas of all the
cavity fuel, free fission gas, and dissolved fission gas in the pins of
one subassembly.
The generalized source or sink term is written as:
.. math::
:label: 15.2-7
S' = \frac{S^{\mathcal{l}}}{\text{AXMX}}
where the source or sink term :math:`S_{\text{l}}` represents a mass source or
sink per unit time and unit length. The primed source or sink term has
the dimension of mass per unit time and per unit smear volume. This unit
of smear volume is a m\ :sup:`3` of the cell volume :math:`\text{AXMX} \cdot \Delta z` in
which all relevant components (including components in all failed pins,
the components in all the channels and the structure) are assumed to be
uniformly smeared. From :eq:`15.2-7`, a change in the generalized smear
density, :math:`{\rho'}`, due to source or sink :math:`S'` is just :math:`S' \Delta t`.
.. _section-15.2.3:
Differential Equations for the In-Pin Fuel Motion and Description of Sink and Source Terms
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The equation set for the in-pin fuel motion includes three mass
conservation equations (for fuel, free fission gas, and dissolved
fission gas), one energy, and one momentum conservation equation. The
continuity equation for the molten fuel in the pin cavity is written:
.. math::
:label: 15.2-8
\frac{\partial}{\partial \text{t}} \left( \rho_{\text{fuca}} A_{\text{fuca}} \right) = \
- \frac{\partial}{\partial \text{z}}\left( \rho_{\text{fuca}} A_{\text{fuca}} u_{\text{fuca}} \right) +\
S_{\text{fuca,me}}^{\mathcal{l}}\left( z,t \right)
where the subscript me refers to fuel melting into the pin cavities of
all pins. By dividing by AXMX and using the definitions of the
generalized smear densities and source and sink terms, :eq:`15.2-8`
becomes:
.. math::
:label: 15.2-9
\frac{\partial}{\partial \text{t}}{\rho'}_{\text{fuca}} = \
- \frac{\partial}{\partial \text{z}}\left( {\rho'}_{\text{fuca}} u_{\text{fuca}} \right) + {S'}_{\text{fuca,me}}
where the primed source is per unit time and per unit smear volume (see
:eq:`15.2-7`). The integrated source term,
:math:`S'_{\text{fuca,me}} \Delta t_{\text{PN}}`, which is actually needed in the
finite difference equations of the code, is calculated from
.. math::
:label: 15.2-10a
{S'}_{\text{fuca,me}} \Delta t_{\text{PN}} = \frac{\rho_{\text{fu,cabd}} \Delta A_{\text{me}} \text{NRPI}}{\text{AXMX}}
where
:math:`\Delta t_{\text{PN}}` = PINACLE time step, s
:math:`\rho_{\text{fu,cabd}}` = Fuel density (including porosity) adjacent to the
cavity boundary kg/m\ :sup:`3`
:math:`\Delta A_{\text{me}}` = Area of fuel (including porosity) melted into the
cavity per PINACLE time step per pin, m\ :sup:`3`
NRPI = Number of pins per subassembly,
When the solid fuel node adjacent to the cavity has not yet exceeded an
input melt fraction value FNMELT:
.. math::
:label: 15.2‑10b
{S'}_{\text{fuca,me}} \Delta t_{\text{PN}} = 0
The :math:`\Delta A_{\text{me}}` in :eq:`15.2-10a` is related to a change in the
cavity diameter by
.. math::
:label: 15.2-11
\Delta A_{\text{me}} = \frac{\pi}{4}\Delta D_{\text{ca}} \left( 2D_{\text{ca}} + \Delta D_{\text{ca}} \right)
where
:math:`D_{\text{ca}}` = Cavity diameter, m
:math:`\Delta D_{\text{ca}}` = Change in the cavity diameter, m
In order to avoid adding the whole radial node instantaneously upon
meeting the input melt fraction criterion FNMELT, the radial node is
added gradually beginning at the time FNMELT is reached and the
addition is completed once the melt fraction has exceeded FNMELT + 0.1.
Once the melt fraction has become greater than FNMELT, the change
in cavity diameter is calculated from
.. math::
:label: 15.2-12
\Delta D_{\text{ca}} = \text{FNMECA} \cdot 2 \cdot \Delta r_{\text{node}}
where :math:`\Delta r_{\text{node}}` is the width of the heat-transfer node
adjacent to the cavity wall before it is melted into the cavity.
FNMECA is the fraction of this node width which has melted into the
cavity per PINACLE time step.
.. math::
:label: 15.2-13
\text{FNMECA} = \Delta t_{\text{PN}} \frac{\left\lbrack \frac{\left( T_{\text{cabd}}^{n + 1} - T_{\text{cabd}}^{n} \right)}{\Delta t_{\text{Ht}}} \right\rbrack}{\left\lbrack 0.1 \cdot \left( T_{\text{fu,liq}} - T_{\text{fu,sol}} \right) \right\rbrack} \
= \frac{\Delta T_{\text{PN}}}{\left\lbrack 0.1 \cdot \left( T_{\text{fu,liq}} - T_{\text{fu,sol}} \right) \right\rbrack} \\
where
:math:`T_{\text{cabd}}^{n}, T_{\text{cabd}}^{n + 1}` = The
temperatures of the fuel node adjacent to the cavity at the beginning
and at the end of the heat-transfer time step :math:`\Delta t_{\text{Ht}}`,
respectively
:math:`T_{\text{fu,liq}}, T_{\text{fu,sol}}` = Fuel
liquidus and solidus temperatures, respectively (The difference between
the two is the melting band width)
:math:`\Delta t_{\text{PN}}` = The PINACLE time step.
:math:`\Delta T_{\text{PN}}` = The temperature change of the fuel
adjacent to the cavity during a PINACLE time step
:eq:`15.2-13` implies that the whole heat-transfer node will be
melted into the cavity after its temperature has risen by 1/10 of the
melting band width beyond the input value **FNMELT**. Moreover, it is
checked whether the neighboring solid fuel node also has exceeded the
input melt fraction criterion FNMELT. If this is the case, the entire
remaining node currently melting into the cavity is added immediately to
the cavity. This situation can occur in TREAT experiments.
The free fission-gas mass conservation equation is written
.. math::
:label: 15.2-14
\frac{\partial}{\partial \text{t}}\left( \rho_{\text{fica}}A_{\text{fica}} \right) = \
- \frac{\partial}{\partial \text{z}}\left( \rho_{\text{fica}} A_{\text{fica}} u_{\text{fuca}} \right) \
+ S_{\text{fica,me}}^{l} + S_{\text{fsca,rl}}^{l}
where the subscripts **fsca** and **rl** refer to fission gas in solution
and to the release of this dissolved fission gas, respectively. Dividing
by AXMX and by using the definition of the generalized smear densities
and source terms yields
.. math::
:label: 15.2-15
\frac{\partial}{\partial \text{t}}{\rho'}_{\text{fica}} = \
- \frac{\partial}{\partial \text{z}}\left( {\rho'}_{\text{fica}} u_{\text{fuca}} \right) \
+ {S'}_{\text{fica,me}} + {S'}_{\text{fsca,rl}}
The integrated source term for free fission gas due to fuel melt-in is
similar to that in :eq:`15.2-10a`:
.. math::
:label: 15.2-16
{S'}_{\text{fica,me}} \Delta t_{\text{PN}} = \rho_{\text{figb,cabd}} \Delta A_{\text{me}} \left( \frac{\text{NRPI}}{\text{AXMX}} \right)
where
:math:`\Delta A_{\text{me}}` = Calculated change in cavity cross sectional area
from :eq:`15.2-11`
:math:`\rho_{\text{figb,cabd}}` = Density of the grain-boundary gas
in the fuel node adjacent to the cavity
The quantity :math:`\rho_{\text{figb,cabd}}` is calculated from:
.. math::
:label: 15.2-17
\rho_{\text{figb,cabd}} = \frac{\text{GNBFG} 2_{\text{cabd}} \cdot \left( \frac{\text{FUMS}_{\text{cabd}}}{\text{FUELMS}_{\text{cabd}}} \right)}{\text{VOLUME}_{\text{cabd}}}
where
:math:`\text{GNBFG2}_{\text{cabd}}` = The grain-boundary fission-gas mass in the
original radial fuel-pin node at the cavity boundary before it has
actually melted in
:math:`\text{FUMS}_{\text{cabd}}` = The current fuel mass in the melting radial
fuel-pin node at the cavity boundary
:math:`\text{FUELMS}_{\text{cabd}}` = The original fuel mass of the radial fuel-pin
node at the cavity boundary before any fuel is removed due to melt-in
:math:`\text{VOLUME}_{\text{cabd}}` = The current volume of the melting radial
fuel-pin node at the cavity boundary
The source term due to dissolved fission-gas release in :eq:`15.2-15` is:
.. math::
:label: 15.2-18
{S'}_{\text{fsca,rl}} = {\rho'}_{\text{fsca}} \cdot \text{CIRTFS} \cdot \Delta t_{\text{PN}}
where
:math:`\text{CIRTFS}` = Release time constant for dissolved fission gas which is
input and has the dimensions s\ :sup:`-1`
The same release constant is also used for the dissolved fission-gas
release in the coolant channels - see :eq:`14.4-20`. This is a relatively
simple exponential decay-type approach to treat the release of the gas
dissolved in molten fuel. However, the understanding of the mechanism of
dissolved gas release from molten fuel is still very limited.
The dissolved fission-gas mass conservation equation is:
.. math::
:label: 15.2-19
\frac{\partial}{\partial \text{t}} \left( \rho_{\text{fsca}} A_{\text{fuca}} \right) = \
- \frac{\partial}{\partial \text{z}}\left( \rho_{\text{fsca}} A_{\text{fuca}} u_{\text{fuca}} \right) \
+ {S'}_{\text{fsca,me}} - {S'}_{\text{fsca,rl}}
By dividing this equation by AXMX and using the definitions of the
generalized smear densities and sources, one obtains:
.. math::
:label: 15.2-20
\frac{\partial}{\partial \text{t}}\left( {\rho'}_{\text{fsca}} \right) = \
- \frac{\partial}{\partial \text{z}}\left( {\rho'}_{\text{fsca}} u_{\text{fuca}} \right) \
+ {S'}_{\text{fsca,me}} - {S'}_{\text{fsca,rl}}
where
.. math::
:label: 15.2-21
{S'}_{\text{fsca,me}} = \rho_{\text{fi,ig,cadb}} \cdot \Delta A_{\text{me}} \cdot \left( \frac{\text{NRPI}}{\text{AXMX}} \right)
:math:`\rho_{\text{fi,ig,cabd}}` = The density of the
intra-granular gas in the fuel node adjacent to the cavity.
The absolute value of the sink term due to the dissolved fission-gas
release has been described in :eq:`15.2-18`.
The fuel energy conservation equation is written:
.. math::
:label: 15.2-22
\frac{\partial}{\partial \text{t}} \left\lbrack \rho_{\text{fuca}} e_{\text{fuca}} A_{\text{fuca}} \right\rbrack = - \frac{\partial}{\partial \text{z}} \left\lbrack \rho_{\text{fuca}} e_{\text{fuca}} A_{\text{fuca}} u_{\text{fuca}} \right\rbrack + S_{\text{fuca,me}}^{l} e_{\text{fu,cabd}} \\
+ Q A_{\text{fuca}} \rho_{\text{fuca}} - h_{\text{fuca,cabd}} \cdot \left( T_{\text{fuca}} - T_{\text{fu,cabd}} \right) \pi D_{\text{ca}} \text{NRPI}
Dividing :eq:`15.2-22` by AXMX and using the definitions of the
generalized smear densities and sources and sinks produces:
.. math::
:label: 15.2-23
\frac{\partial}{\partial \text{t}}\left\lbrack {\rho'}_{\text{fuca}} e_{\text{fuca}} \right\rbrack = - \frac{\partial}{\partial \text{z}}{\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} + {S'}_{\text{fuca,me}} e_{\text{fu,cabd}} \\
+ Q {\rho'}_{\text{fuca}} - h_{\text{fuca,cabd}} \cdot \left( T_{\text{fuca}} - T_{\text{fu,cabd}} \right) \pi D_{\text{ca}} \frac{\text{NRPI}}{\text{AXMX}}
By rewriting the left-hand side of :eq:`15.2-23` as
.. math::
:label: 15.2-24
\frac{\partial}{\partial \text{t}}\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} \right) \
= {\rho'}_{\text{fuca}} \frac{\partial}{\partial \text{t}} e_{\text{fuca}} + e_{\text{fuca}} \
\frac{\partial}{\partial \text{t}} {\rho'}_{\text{fuca}}
and by using the mass conservation :eq:`15.2-8`, one obtains
.. math::
:label: 15.2-25
{\rho'}_{\text{fuca}} \frac{\partial}{\partial \text{t}} e_{\text{fuca}} \
= \frac{\partial}{\partial \text{t}}\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right) \
+ e_{\text{fuca}}\frac{\partial}{\partial \text{z}} \left( {\rho'}_{\text{fuca}} u_{\text{fuca}} \right) \\
+ {S'}_{\text{fuca,me}}\left( e_{\text{fu,cabd}} - e_{\text{fuca}} \right) + Q {\rho'}_{\text{fuca}} \\
- h_{\text{fuca,cabd}} \left( T_{\text{fuca}} - T_{\text{fu,cabd}} \right) \pi D_{\text{ca}} \frac{\text{NRPI}}{\text{AXMX}}
where :math:`Q` is the fission heat source per kg of fuel
.. math::
:label: 15.2-25a
Q = \text{FNPOHE} \cdot \text{POW} \cdot \text{PSHAPE} \left( K \right) \cdot \\
\left( 1 - \text{GAMSS} - \text{GAMTNC} - \text{GAMTNE} \right) \cdot \frac{F_{\text{POWER}}}{\text{FUMASS}\left( K \right)}
In :eq:`15.2-25a`,
.. math::
:label: 15.2-25b
\text{FNPOHE} = \exp\left\lbrack \text{POWCOF}\left( 1 \right) + \Delta t_{\text{x}} \left\{ \text{POWCOF}\left( 2 \right) \
+ \Delta t_{\text{x}} \text{POWCOF}\left( 3 \right) \right\} \right\rbrack
where :math:`\Delta t_{\text{x}}` is the time between the current *PINACLE* time
and the beginning of the current main (point kinetics) time step. The
POWCOF's are the coefficients that are found by fitting an
exponential function to the power levels at the end of the last three
point kinetics time step.
.. list-table::
:header-rows: 0
:align: center
:widths: auto
* - :math:`\text{POW}`
- =
- Steady-state power in the peak axial fuel-pin segment (see SAS4A input description)
* - :math:`\text{PSHAPE} \left( K \right)`
- =
- Ratio of pin power at axial node :math:`K` to :math:`\text{POW}` (see SAS4A input description)
* - :math:`\text{GAMSS, GAMTNC, GAMTNE}`
- =
- Fractions of total power for the direct heating of structure, coolant, and cladding, respectively (See SAS4A input description)
* - :math:`\text{FUMASS} \left( K \right)`
- =
- Initial total fuel mass in axial pin segment :math:`K`
* - :math:`F_{\text{POWER}}`
- =
- Power reduction factor if there is a radial power gradient in the pin (as is common in *TREAT* experiments):
.. math::
:label: 15.2-26
F_{\text{POWER}} = \frac{\sum_{\text{I=1}}^{I_{\text{cabd}}}{\text{PSHAPR}\left( I \right) \cdot \text{FUELMS} \
\left( I,K \right) \big/ \sum_{\text{I=1}}^{I_{\text{cabd}}}{\text{FUELMS}\left( I,K \right)}}}{\sum_{\text{I=1}}^{\text{NT}}{\text{PSHAPR}\
\left( I \right) \cdot \text{FUELMS} \left( I,K \right)} \big/ \sum_{\text{I=1}}^{\text{NT}}{\text{FUELMS}\left( I,K \right)}}
where
:math:`\text{NT}` = Number of radial pin nodes
:math:`I_{\text{cabd}}` = Number of radial pin nodes in the cavity
:math:`\text{PSHAPR} \left( I \right)` = Mass normalized radial power distribution in radial node
:math:`I` (See SAS4A input description)
:math:`\text{FUELMS} \left( I,K \right)` = Initial fuel mass in the radial fuel-pin node
:math:`I, K`
The heat-transfer coefficient in :eq:`15.2-22` is the sum of a convective
and a conduction heat-transfer term.
.. math::
:label: 15.2-27
h_{\text{fuca,cabd}} = \left( h_{1} + h_{2} \right)
where
.. math::
:label: 15.2-27a
h_{1} = \frac{k_{\text{fu}}}{D_{\text{ca}}} \cdot \text{St} \cdot \text{Pr} \cdot \text{Re}
This comes from the definition of the Nusselt number
.. math::
:label: 15.2-27b
\text{Nu} = \text{St} \cdot \text{Pr} \cdot \text{Re}
where
:math:`\text{St}` = Stanton number = :math:`\frac{h}{\rho u C_{\text{p}}}`
:math:`\text{Pr}` = Prandtl number in :math:`\frac{\mu C_{\text{p}}}{k}`
:math:`\text{Re}` = Reynolds numbers = :math:`\frac{D \rho u}{\mu}`
The Deissler correlation [15-7] was used for finding the relationship
between the three nondimensional numbers on the right-hand side of
:eq:`15.2-27b`. The Prandtl number for fuel is about 2.2. By using this value
in the Deissler correlation, it can be shown that
.. math::
:label: 15.2-28
\text{St} \approx 0.0158\ \text{Re}^{- 0.2}
By using :eq:`15.2-28` and the definition of the Prandtl number,
.. math::
:label: 15.2-29
h_{1} = \frac{1}{D_{\text{ca}}} \cdot \mu_{\text{fu,liq}} \cdot C_{\text{P,fu}} \cdot \text{CIA3} \cdot \text{Re}_{\text{ca}}^{0.8}
where
:math:`k_{\text{fu}}` = Conductivity of fuel which is input (CDFU)
:math:`\mu_{\text{fu,liq}}` = Liquid fuel viscosity which is input
(VIFULQ)
:math:`\text{CIA3}` = Input constant. A value of 0.0158 is recommended because of :eq:`15.2-28`
.. math::
:label: 15.2-30
Re_{\text{ca}} = \frac{\left( \left| u_{\text{fuca}} \right| D_{\text{ca}} {\rho'}_{\text{fuca}} \
\big/ \theta_{\text{ca}} \right)}{\mu_{\text{fu,liq}}}
The conduction heat-transfer coefficient, :math:`h_{2}`, which is
relevant for a low flow or a stagnant flow condition is of the following
form
.. math::
:label: 15.2-31
h_{2} = 4 \cdot \frac{k_{\text{fu}}}{D_{\text{ca}}}
The pressure calculation in the fuel-pin cavity is based on the
assumption that the fission-gas pressure and fuel-vapor pressure can be
added. The total cavity pressure is:
.. math::
:label: 15.2-32
P_{\text{ca}} = P_{\text{fica}}\left( T_{\text{fuca}}, \rho_{\text{fica}} \right) + P_{\text{fvca}} \left( T_{\text{fuca}} \right)
If the input variable :sasinp:`INAPN` is equal to 1, then the sodium
vapor pressure contribution will also be included in the pressure
calculation, accounting for the possible presence of liquid sodium in
metal fuel pins:
.. math::
:label: 15.2-32a
P_{\text{ca}} = P_{\text{fica}}\left( T_{\text{fuca}}, \rho_{\text{fica}} \right) \
+ P_{\text{fvca}} \left( T_{\text{fuca}} \right) + P_{\text{Na}}
where the sodium vapor pressure contribution is calculated as follows:
.. math::
:label: 15.2-32b
P_{\text{Na}} = P_{\text{Nv,ca}} \left( T_{\text{fuca}} \right) - \left\lbrack P_{\text{fica}}\left( T_{\text{fuca}},\
\rho_{\text{fica}} \right) + P_{\text{fvca}} \right\rbrack \cdot \text{CINAPN}
where :sasinp:`CINAPN` is an input constant with values between 0
and 1.
The fission-gas pressure in :eq:`15.2-32` and :eq:`15.2-32a` is calculated from
a special form of the ideal-gas equation which takes the compressibility
of the liquid fuel into account:
.. math::
:label: 15.2-33
P_{\text{fica}} = \frac{R_{\text{fi}} \cdot T_{\text{fuca}} {\rho'}_{\text{fica}}}{\left( \theta_{\text{ca}} \
- \theta_{\text{fuca}} + \theta_{\text{fuca}} K_{\text{fu}} \cdot P_{\text{fica}} \right)}
where
:math:`R_{\text{fi}}` = The universal gas constant divided by the molecular
weight of the mixture of fission gas and helium fill gas (:math:`R_{\text{fi}}` is
equal to an input value RGAS which should take into account the relative
amounts of krypton, xenon, and helium. The latter may be important for
near-fresh fuel)
:math:`K_{\text{fu}}` = Liquid fuel compressibility which is input (see CMFU)
The physically meaningful solution of the quadratic :eq:`15.2-33` is:
.. math::
:label: 15.2-34
P_{\text{fica}} = \frac{- \left( \theta_{\text{ca}} - \theta_{\text{fuca}} \right) \
+ \sqrt{\left( \theta_{\text{ca}} - \theta_{\text{fuca}} \right)^{2} \
+ 4 \theta_{\text{fuca}} K_{\text{fu}} \cdot R_{\text{fi}} \cdot T_{\text{fuca}}\
{\rho'}_{\text{fica}}}}{2 \theta_{\text{fuca}} K_{\text{fu}}}
There is also a second solution with a minus sign in front of the square
root which does not give a physically meaningful result.
:eq:`15.2-34` reduces to a single-phase liquid pressure solution for no
fission gas and :math:`\theta_{\text{fuca}} > \theta_{\text{ca}}`:
.. math::
:label: 15.2-35
P_{\text{fica}} = - \frac{\left( \theta_{\text{ca}} - \theta_{\text{fuca}} \right)}{\left( \theta_{\text{fuca}} - K_{\text{fu}} \right)}
which is equivalent to the definition of the fuel compressibility. For
void fractions greater than 30% the fission-gas pressure is calculated
from a simplified form of :eq:`15.2-33` in which the term with
:math:`K_{\text{fu}}` is dropped.
The momentum conservation equation for the fuel/fission-gas mixture is:
.. math::
:label: 15.2-36
\frac{\partial}{\partial \text{t}}\left( \rho_{\text{fuca}} A_{\text{fuca}} u_{\text{fuca}} + \rho_{\text{fica}} A_{\text{fica}} u_{\text{fuca}} \right) \\
= - \frac{\partial}{\partial \text{z}}\left( \rho_{\text{fuca}} A_{\text{fuca}} u_{\text{fuca}}^{2} + \rho_{\text{fica}} A_{\text{fica}} u_{\text{fuca}}^{2} \right) \\
- A_{\text{ca}} \cdot \left( \frac{\partial \text{P}_{\text{ca}}}{\partial \text{z}} \right) - g \left( \rho_{\text{fuca}} A_{\text{fuca}} + \rho_{\text{fica}} A_{\text{fica}} \right) \\
- u_{\text{fuca}}\left| u_{\text{fuca}} \right| \cdot \left( \rho_{\text{fuca}} \frac{A_{\text{fuca}}}{A_{\text{ca}}} + \rho_{\text{fica}} \frac{A_{\text{fica}}}{A_{\text{ca}}} \right) A_{\text{ca}} \cdot F_{\text{friction}} \big/ \left( 2D_{\text{ca}} \right)
Dividing :eq:`15.2-36` by AXMX and using the definitions for the
generalized smear densities and mass sinks gives:
.. math::
:label: 15.2-37
\frac{\partial}{\partial \text{t}}\left\lbrack u_{\text{fuca}}\left( {\rho'}_{\text{fuca}} + {\rho'}_{\text{fica}} \right) \right\rbrack \\
= \frac{\partial}{\partial \text{z}}\left\lbrack u_{\text{fuca}}^{2}\left( {\rho'}_{\text{fuca}} + {\rho'}_{\text{fica}} \right) \right\rbrack \\
- \theta_{\text{ca}} \left( \frac{\partial \text{P}_{\text{ca}}}{\partial \text{z}} \right) - g \left( {\rho'}_{\text{fuca}} + {\rho'}_{\text{fica}} \right) \\
- u_{\text{fuca}} \left| u_{\text{fuca}} \right| \cdot \left( {\rho'}_{\text{fuca}} + {\rho'}_{\text{fica}} \right) \frac{F_{\text{friction}}}{\left( 2D_{\text{ca}} \right)}
where
:math:`u_{\text{fuca}}` = Upward fuel velocity in the cavity (see the sign of
the gravity head term in :eq:`15.2-36`)
The Moody friction factor :math:`F_{\text{friction}}` in :eq:`15.2-37` depends on
the Reynolds number:
.. math::
:label: 15.2-38
Re_{\text{pi}} = \frac{\left| u_{\text{fuca}} \right| D_{\text{ca}} \cdot \left( {\rho'}_{\text{fuca}} \
+ {\rho'}_{\text{fica}} \right)}{\left( \theta_{\text{ca}} \cdot \text{VIFULQ} \right)}
where VIFULQ is the viscosity of liquid fuel which is input. The
friction factor is
.. math::
:label: 15.2-39
F_{\text{frication}} = \begin{cases}
\frac{64}{\text{Re}_{\text{pi}}} & \text{for } \text{Re}_{\text{pi}} < \text{CIREFU} \\
\text{CIFRFU} & \text{for } \text{Re}_{\text{pi}} > \text{CIREFU} \\
\end{cases}
where CIREFU and CIFRFU are both input and should be made consistent
in order to avoid a jump in the friction factor at :math:`\text{Re}_{\text{pi}} = \text{CIREFU}`.
There is no term accounting for the fuel melt-in because this fuel is
added to the cavity with zero axial velocity, and therefore does not
change the total momentum. However, since the generalized fuel smear
density will change in a cell with melt-in, this will lead to a velocity
decrease in such a cell.
.. _section-15.2.4:
Finite Difference Equations for the In-Pin Motion Model
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In the overview of the numerical scheme given in :numref:`section-15.2.1`, it was
pointed out that full donor cell spatial differencing and a largely
explicit time differencing are used for treating the in-pin motion. The
implicit aspect of the solution is that the mass and energy conservation
equations are solved first and then a pressure is calculated on the
basis of the mass and energy equation results. This advanced pressure is
used in the momentum conservation equation.
The finite differencing of all the mass conservation equations is the
same. The fuel mass conservation is used as an illustration.
.. math::
:label: 15.2-40
\frac{\left( {\rho'}_{\text{fuca,K}}^{n + 1} - {\rho'}_{\text{fuca,K}}^{n} \right)}{\Delta t_{\text{PN}}} = \
- \frac{\left( \left( {\rho'} u \right)_{\text{fuca,K+1}} - \left( {\rho'} u \right)_{\text{fuca,K}} \right)}{\Delta z_{\text{K}}} \
+ \sum_{\text{k}}{{S'}_{\text{K,k}}}
where
.. math::
:label: 15.2-40a
\left( {\rho'} u \right)_{\text{fuca,K}} = \begin{cases}
{\rho'}_{\text{fuca,K-1}} u_{\text{K}} & \text{for } u_{\text{K}} > 0 \\
{\rho'}_{\text{fuca,K}} u_{\text{K}} & \text{for } u_{\text{K}} < 0 \\
\end{cases}
:math:`u_{\text{K}}` = Fuel velocity at the mesh cell boundary :math:`K`.
:math:`{\rho'}_{\text{fuca,K-1}}` = Generalized fuel smear density at the
mesh cell midpoint below boundary :math:`K`
:math:`{\rho'}_{\text{fuca,K}}` = Generalized fuel smear density at the mesh
cell midpoint above boundary :math:`K`
The numerical grid used in the program was discussed in :numref:`section-15.2.1`
which gives an overview of the numerical scheme. Densities are defined
at the cell centers and velocities at the cell edges as illustrated in
:numref:`figure-15.2-2`. The source and sink terms have already been described in
their finite difference form in the previous section.
.. _figure-15.2-2:
.. figure:: media/image2-2.png
:align: center
:figclass: align-center
Numerical Grid Used in PINACLE
The convective fluxes at the lower and upper boundaries of the cavity
which are located in the end cells KK1 and KKMX, respectively, are
.. math::
:label: 15.2-41
\left( {\rho'} u \right)_{\text{fuca,KK1}} = 0
and
.. math::
:label: 15.2-42
\left( {\rho'} u \right)_{\text{fuca,KKMX+1}} = 0
The end cells are not always the same during a PINACLE run because the
molten cavity can extend axially. When a new cell is added to the
cavity, mass can only flow into this end cell if its cross section is at
least 20% of that of the neighboring cell in the molten cavity. This
restriction had to be included because of problems with overcompression
of cells with a very small cross section.
The finite difference form of the fuel energy :eq:`15.2-25` is:
.. math::
:label: 15.2-43
\frac{{\rho'}_{\text{fuca,K}} \left( e_{\text{fuca,K}}^{n + 1} - e_{\text{fuca,K}}^{n} \right)}{\Delta t_{\text{PN}}} = - \left\lbrack \left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K+1}}^{n} \right. \\ \left.
- \left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K}}^{n} \right\rbrack \big/ \Delta z_{\text{k}} + e_{\text{fuca,K}}^{n}\left\lbrack \left( {\rho'}_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K+1}}^{n} \right. \\ \left.
- \left( {\rho'}_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K}}^{n} \right\rbrack \big/ \Delta z_{\text{k}} + {S'}_{\text{fuca,me,K}}^{n} \left( e_{\text{fu,cabd,K}}^{n} - e_{\text{fuca,K}}^{n} \right) \\
+ Q_{\text{K}}^{n} {\rho'}_{\text{fuca,K}}^{n} - h_{\text{fuca,cabd,K}}^{n} \left( T_{\text{fuca,K}}^{n} - T_{\text{fu,cabd,K}}^{n} \right) \pi D_{\text{ca,K}}^{n} \frac{\text{NRPI}}{\text{AXMX}}
The fuel cavity temperature :math:`T_{\text{fuca}}` in the above equation has
to be calculated from the internal energy:
.. math:: \text{For } e_{\text{fu,sol}} < e_{\text{fuca}} < e_{\text{fu,liq}}
.. math::
:label: 15.2-44
T_{\text{fuca}} = T_{\text{fu,sol}} + \left( T_{\text{fu,liq}} - T_{\text{fu,sol}} \right) \cdot \frac{\left( e_{\text{fuca}}
- e_{\text{fu,sol}} \right)}{\left( e_{\text{fu,liq}} - e_{\text{fu,sol}} \right)}
.. math:: \text{For } e_{\text{fuca}} > e_{\text{fu,liq}}
.. math::
:label: 15.2-44a
T_{\text{fuca}} = T_{\text{fu,liq}} + \frac{\left( e_{\text{fuca}} - e_{\text{fu,liq}} \right)}{C_{\text{p,fu}}}
where :math:`C_{\text{p,fu}}` is the fuel specific heat which is
the single input value CPFU. The convective energy flux at cell
boundary :math:`K` in :eq:`15.2-43` is calculated as:
.. math::
:label: 15.2-45
\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K}} = \begin{cases}
\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} \right)_{\text{K-1}} u_{\text{fuca,K}} & \text{for } u_{\text{fuca,K}} > 0 \\
\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} \right)_{\text{K}} u_{\text{fuca,K}} & \text{for } u_{\text{fuca,K}} < 0 \\
\end{cases}
At the lower and upper cavity ends, which are in cells KK1 and KKMX, the
convective energy fluxes are zero:
.. math::
:label: 15.2-45a
\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right)_{\text{KK1}}^{n} = 0
and
.. math::
:label: 15.2-45b
\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right)_{\text{KKMX+1}}^{n} = 0
By using the convective fluxes from :eq:`15.2-40a` and :eq:`15.2-45`, and
the definitions of the energy gain and loss terms given earlier, and by
differencing the second term of :eq:`15.2-43` like the first, :eq:`15.2-43`
can be solved for :math:`e_{\text{fuca,K}}^{n + 1}`. Fuel temperatures
that are shown in the PINACLE output are calculated by using
:eq:`15.2-44` and :eq:`15.2-44a`.
For the finite difference form of the momentum equation the following
quantities have to be defined at the edges of the numerical cells: the
combined fuel/fission-gas generalized smear density and the cavity
volume fraction. These quantities become:
.. math::
:label: 15.2-46
{\rho'}_{\text{fufi,bk}} = 0.5 \left( {\rho'}_{\text{fuca,K-1}} + {\rho'}_{\text{fuca,K}} \right)
+ 0.5 \left( {\rho'}_{\text{fica,K-1}} + {\rho'}_{\text{fica,K}} \right)
.. math::
:label: 15.2-47
\theta_{\text{ca,bk}} = 0.5 \left( \theta_{\text{ca,K}} + \theta_{\text{ca,K-1}} \right)
where the subscript **bk** indicates that these quantities are at the
lower boundaries of cell :math:`K`. This is shown in :numref:`figure-15.2-2`. On the
numerical grid the velocities are also defined on the cell boundaries,
whereas the pressures, densities and volume fractions are defined at the
cell centers. This is also shown in :numref:`figure-15.2-2`. The finite difference
form of the momentum conservation, :eq:`15.2-37`, is written using
:eq:`15.2-46` and :eq:`15.2-47` as:
.. math::
:label: 15.2-48
\frac{\left( {\rho'}_{\text{fufi,bk}}^{n + 1} u_{\text{fuca,K}}^{n + 1} - {\rho'}_{\text{fufi,bk}}^{n} u_{\text{fuca,K}}^{n} \right)}{\Delta t_{\text{PN}}} = - \left\lbrack \left( u_{\text{fuca}}^{2} {\rho'}_{\text{fufi}} \right)_{\text{K}}^{n} \right. \\ \left.
- \left( u_{\text{fuca}}^{2} {\rho'}_{\text{fufi}} \right)_{\text{K-1}}^{n} \right\rbrack \big/ \Delta z - \theta_{\text{ca,bk}}\left\lbrack \left( 1 - \varepsilon \right) \cdot \left( P_{\text{ca,K}}^{n} - P_{\text{ca,K-1}}^{n} \right) \right. \\ \left.
+ \varepsilon \cdot \left( P_{\text{ca,K}}^{n + 1} - P_{\text{ca,K-1}}^{n + 1} \right) \right\rbrack \big/ \Delta z - g {\rho'}_{\text{fufi,bk}} \\
- u_{\text{fuca,K}}^{n + 1} \left| u_{\text{fuca,K}}^{n} \right| {\rho'}_{\text{fufi,bk}}^{n} \frac{F_{\text{friction,bk}}}{\left( 2D_{\text{ca}} \right)}
where
:math:`\varepsilon` = Input value EPCH that can be between zero and one (see :eq:`15.2-1`)
:math:`\Delta z = 0.5 \left( \Delta z_{\text{K-1}} + \Delta z_{\text{K}} \right)`
By collecting all terms with :math:`u_{\text{fuca}}^{n + 1}` on the
left-hand side of the equation, one obtains
.. math::
:label: 15.2-49
u_{\text{fuca,K}}^{n + 1}\left\lbrack \frac{{\rho'}_{\text{fufi,bk}}^{n + 1}}{\Delta t_{\text{PN}}} + \left| u_{\text{fuca,K}}^{n} \right| {\rho'}_{\text{fufi,bk}}^{n} \frac{F_{\text{friction,bk}}}{\left( 2D_{\text{ca}} \right)} \right\rbrack \\
= \frac{{\rho'}_{\text{fufi,bk}}^{n + 1} u_{\text{fuca,K}}^{n}}{\Delta t_{\text{PN}}} \\
- \frac{\left\lbrack \left( u_{\text{fuca}}^{2} {\rho'}_{\text{fufi}} \right)_{\text{K}}^{n} - \left( u_{\text{fuca}}^{2} {\rho'}_{\text{fufi}} \right)_{\text{K-1}}^{n} \right\rbrack}{\Delta z} \\
- \theta_{\text{ca,bk}} \left\lbrack \left( 1 - \varepsilon \right) \cdot \left( P_{\text{ca,K}}^{n} - P_{\text{ca,K-1}}^{n} \right) \right. \\ \left.
+ \varepsilon \left( P_{\text{ca,K}}^{n + 1} - P_{\text{ca,K-1}}^{n + 1} \right) + \left( P_{\text{vi,K}}^{n} - P_{\text{vi,K-1}}^{n} \right) \right\rbrack \Delta z - g {\rho'}_{\text{fufi,bk}}
The convective momentum flux in :eq:`15.2-49` is calculated as
.. math::
:label: 15.2-50
\left( u_{\text{fuca}}^{2} {\rho'}_{\text{fufi}} \right)_{\text{K}} = \begin{cases}
{\rho'}_{\text{fufi,K}} u_{\text{fuca,K}}^{2} & \text{if } \left( u_{\text{fuca,K}} + u_{\text{fuca,K+1}} \right) > 0 \\
{\rho'}_{\text{fufi,K}} u_{\text{fuca,K+1}}^{2} & \text{if } \left( u_{\text{fuca,K}} + u_{\text{fuca,K+1}} \right) < 0
\end{cases}
.. math::
:label: 15.2-51
\left( u_{\text{fuca}}^{2} {\rho'}_{\text{fufi}} \right)_{\text{K-1}} = \begin{cases}
{\rho'}_{\text{fufi,K-1}} u_{\text{fuca,K-1}}^{2} & \text{if } \left( u_{\text{fuca,K-1}} + u_{\text{fuca,K}} \right) > 0 \\
{\rho'}_{\text{fufi,K-1}} u_{\text{fuca,K}}^{2} & \text{if } \left( u_{\text{fuca,K-1}} + u_{\text{fuca,K}} \right) < 0
\end{cases}
where
.. math::
:label: 15.2-52
{\rho'}_{\text{fufi,K}} = {\rho'}_{\text{fuca,K}} + {\rho'}_{\text{fica,K}}
The momentum fluxes for the lower and upper end cells of the cavity,
which are designated by KK1 and KKMX, are:
.. math::
:label: 15.2-53
\left( u_{\text{fuca}}^{2} {\rho'}_{\text{fufi}} \right)_{\text{KK1}} = \begin{cases}
{\rho'}_{\text{fufi,KK1}} u_{\text{fuca,KK1+1}}^{2} \cdot 0.25 & \text{if } u_{\text{fuca,KK1+1}} > 0 \\
{\rho'}_{\text{fufi,KK1}} u_{\text{fuca,KK1+1}}^{2} & \text{if } u_{\text{fuca,KK1+1}} < 0
\end{cases}
.. math::
:label: 15.2-54
\left( u_{\text{fuca}}^{2} {\rho'}_{\text{fufi}} \right)_{\text{KKMX}} = \begin{cases}
{\rho'}_{\text{fufi,KKMX}} u_{\text{fuca,KKMX}}^{2} & \text{if } u_{\text{fuca,KKMX}} > 0 \\
{\rho'}_{\text{fufi,KKMX}} u_{\text{fuca,KKMX}}^{2} \cdot 0.25 & \text{if } u_{\text{fuca,KKMX}} < 0
\end{cases}
The factor 0.25 in the above convection terms comes from the assumption
of a zero velocity at the end of the cavity.
The momentum :eq:`15.2-49` can be solved for
:math:`u_{\text{fuca,K}}^{n + 1}` if :eq:`15.2-46`, :eq:`15.2-47`, and :eq:`15.2-50`
through :eq:`15.2-54` are used.
.. _section-15.2.5:
Treatment of the Fuel Ejection Above the Top of the Active Fuel Pin
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The equations presented in the preceding sections apply to all the
cavity cells before the cavity reaches the top fuel cell and the upward
axial fuel ejection is initiated. When the fuel ejection begins,
however, a different situation is created at the top of the fuel pin,
through the creation of a new cell of variable length, as illustrated in
:numref:`figure-15.2-3`. The presence of the axial fuel ejection and of the variable
length space above the fuel pin require a special treatment of the top
cell, as outlined below.
.. _section-15.2.5.1:
Initiation of the Axial Fuel Ejection
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
When PINACLE is initiated it models a bottled-up cavity, extending
axially from the cell KK1 to KKMX. Only limited fuel relocation occurs
during this period, as the pressure gradients in the cavity are very
small. For the initiation of the axial fuel motion two conditions must
be satisfied: 1.) The cavity must extend to the upper active
.. _figure-15.2-3:
.. figure:: media/image8.png
:align: center
:figclass: align-center
Geometry Used in the Calculation of the Fuel Ejection Above the Top of the Active Fuel
fuel cell, KCORE2, and 2.) the temperature of the top fuel pin surface
must be higher than the fuel solidus temperature, i.e.:
.. math::
:label: 15.2-55
T_{\text{fu,tp}} > T_{\text{fu,sol}}
The temperature of the top fuel surface :math:`T_{\text{fu,tp}}`
is calculated as follows:
- if blanket pellets are present:
.. math::
:label: 15.2-56
T_{\text{fu,tp}} = T_{\text{fu}} \left( 1,\text{KCORE2} + 1 \right) \\
+ \left( T_{\text{fu}}\left( 1,\text{KCORE2} \right) - T_{\text{fu}}\left( 1,\text{KCORE2} + 1 \right) \right) \cdot \text{CIPNTP}
- if no blanket pellets are present and the pin analyzed contains oxide
fuel no axial temperature gradients are present at the fuel pin top
and the constraint 15.2-55 is eliminated by using:
.. math::
:label: 15.2-57
T_{\text{fu,tp}} = T_{\text{fu,liq}}
- if no blanket pellets are present and the pin analyzed contains metal
fuel:
.. math::
:label: 15.2-58
T_{\text{fu,tp}} = T_{\text{Na,plenum}} + \left( T_{\text{fu}}\left( 1,\text{KCORE2} \right)
- T_{\text{Na,plenum}} \right) \cdot \text{CIPNTP}
where:
:sasinp:`CIPNTP` is an input variable. The value 0.5, used
in the TS-2 and M2 and M3 TREAT tests analyses has proven satisfactory
in predicting the time of fuel motion onset
:math:`T_{\text{fu}} \left( 1,\text{KCORE2} \right)` = The fuel temperature in the axial cell
KCORE2 and radial cell 1 (i.e. the central cell)
After the initiation of the axial fuel ejection the PINACLE calculations
are extended to the cell :math:`\text{KKMXPN}`:
.. math::
:label: 15.2-59
\text{KKMXPN} = \text{KCORE2} + 1
.. _section-15.2.5.2:
Mass Conservation Equation for the Top Cell
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
In order to minimize numerical complications due to the variable length
cell KKMXPN, the cell length used in the mass conservation equation is
:math:`\Delta z_{\text{KKMXPN}}^{\text{REF}}` the axial length of the
segment KKMXPN (see :numref:`figure-15.2-3`). The use of this fixed reference
length allows the use of :eq:`15.2-40`, unmodified, for the calculation of
the :math:`\rho_{\text{fu,ca,K}}^{n + 1}`:
.. math::
:label: 15.2-60
{\rho'}_{\text{fuca,K}}^{n + 1} = {\rho'}_{\text{fuca,K}} - \left\lbrack \left( {\rho'} u \right)_{\text{fuca,K+1}} - \left( {\rho'} u \right)_{\text{fuca,K}} \right\rbrack \\
\cdot \frac{\Delta t_{\text{PN}}}{\Delta z_{\text{KKMXPN}}^{\text{REF}}} + S_{\text{K}}
where:
:math:`K = \text{KKMXPN} = \text{KCORE2} + 1`
:math:`S_{\text{k}}` = The integrated source term :math:`{S'}_{\text{K}} \cdot \Delta t_{\text{PN}}`,
which is always zero in the :math:`\text{KCORE2} + 1` cell
:math:`\left( {\rho'} u \right)_{\text{fuca,K+1}}` = The convective mass flux at the
boundary :math:`\text{KKMXPN} + 1` is always zero
.. math::
:label: 15.2-61
\left( {\rho'} u \right)_{\text{fuca,K}} = \begin{cases}
{\rho'}_{\text{fuca,K-1}} \cdot u_{\text{K}} & \text{for } u_{\text{K}} > 0 \\
{\rho'}_{\text{fuca,K}} \cdot \frac{\Delta z_{\text{KKMXPN}}^{\text{REF}}}{\Delta z_{\text{KKMXPN}}} \cdot u_{\text{K}} \cdot \frac{\theta_{\text{K-1}}}{\theta_{\text{K}}} & \text{for } u_{\text{K}} < 0
\end{cases}
Note that the definition of the convective mass flux is more complex
when :math:`u_{\text{K}} < 0`. The smeared density :math:`{\rho'}_{\text{fuca,K}}`
is multiplied by the factor:
.. math:: \frac{\Delta z_{\text{KKMXPN}}^{\text{REF}}}{\Delta z_{\text{KKMXPN}}}
To obtain the actual smeared density in the existing cell, and the
velocity :math:`u_{\text{K}}`, which is defined at the boundary KKMXPN, but
below the boundary, is multiplied by the factor.
.. math:: \frac{\theta_{K-1}}{\theta_{K}}
To account for the abrupt area change present at the boundary :math:`\text{KKMXPN} = \text{KCORE2} + 1`.
.. _section-15.2.5.3:
The Energy Conservation Equation for the Top Cell
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Due to the use of the reference cell
:math:`\Delta z_{\text{KKMXPN}}^{\text{REF}}` conservation equation
15.2-43 remains virtually unchanged for the top cell :math:`\text{KCORE2} + 1`:
.. math::
:label: 15.2-62
e_{\text{fuca,K}}^{n + 1} = e_{\text{fuca,K}}^{n}\left\lbrack \left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K+1}}^{n} - \left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K}}^{n} \right\rbrack \\
\cdot \frac{\Delta t_{\text{PN}}}{\Delta z_{\text{KKMXPN}}^{\text{REF}} \cdot {\rho'}_{\text{fuca,K}}} + e_{\text{fuca,K}}^{n} \cdot \left\lbrack \left( {\rho'}_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K+1}}^{n} - \left( {\rho'}_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K}}^{n} \right\rbrack \\
\cdot \frac{\Delta t_{\text{PN}}}{\Delta z_{\text{KKMXPN}}^{\text{REF}} \cdot {\rho'}_{\text{fuca,K}}} + \left\lbrack {S'}_{\text{fuca,me,K}}^{n} \cdot \left( e_{\text{fu,cabd,K}}^{n} - e_{\text{fuca,K}}^{n} \right) \right. \\ \left.
+ Q_{\text{K}}^{n} {\rho'}_{\text{fuca,K}}^{n} - h_{\text{fuca,cabd,K}}^{n} \cdot \left( T_{\text{fuca,K}}^{n} - T_{\text{fu,cabd,K}}^{n} \right) \cdot \pi \cdot D_{\text{ca,K}}^{n} \right. \\ \left.
\cdot \frac{\text{NPRI}}{\text{AXMX}} \cdot \frac{\Delta z_{\text{KKMXPN}}}{\Delta z_{\text{KKMXPN}}^{\text{REF}}} \right\rbrack \cdot \frac{\Delta t_{\text{PN}}}{{\rho'}_{\text{fuca,K}}}
where
:math:`K = \text{KKMXPN} = \text{KCORE2} + 1`
:math:`{S'}_{\text{fuca,me,K}}` = The melting fuel source which is always zero in
the top cell.
:math:`\left( {\rho'}_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K}}` = The convective mass term
explained in :numref:`section-15.2.5.2`.
:math:`\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K+1}}` = The
convective energy term at boundary :math:`\text{KKMXPN}+1` and is always zero
:math:`\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right)_{\text{K}}` = The
convective energy term at boundary KKMXPN and is defined below
.. math::
:label: 15.2-63
\left( {\rho'} e u \right)_{\text{fuca,K}} = \begin{cases}
\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} \right)_{\text{K-1}} \cdot u_{\text{fuca,K}} & \text{for } u_{\text{fuca,K}} > 0 \\
{\rho'}_{\text{fuca,K}} \cdot \frac{\Delta z_{\text{KKMXPN}}^{\text{REF}}}{\Delta z_{\text{KKMXPN}}} \cdot e_{\text{fuca,K}} \cdot u_{\text{fuca,k}} \cdot \frac{\theta_{\text{K-1}}}{\theta_{\text{K}}} & \text{for } u_{\text{fuca,k}} < 0
\end{cases}
The explanation of the correction terms
:math:`\frac{\Delta z^{\text{REF}}}{\Delta z_{\text{KKMXPN}}}` and
:math:`\frac{\theta_{\text{K-1}}}{\theta_{\text{K}}}`
present in the convective energy term when
:math:`u_{\text{fuca,K}} < 0` has been given in :numref:`section-15.2.5.2`.
Note also in :eq:`15.2-62` that the term describing the heat transfer
between the molten fuel and the wall contains the correction factor:
.. math:: \frac{\Delta z_{\text{KKMXPN}}}{\Delta z_{\text{KKMXPN}}^{\text{REF}}}
to account for the fact that the heat transfer surface is determined by
the actual length of the top cell. Also for this top cell, the
temperature :math:`T_{\text{fu,cabd,K}}` is set equal to the
inner cladding temperature.
.. _section-15.2.5.4:
The Momentum Conservation Equation for the Top Cell
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The momentum conservation equation in the top cell is used to calculate
the fuel velocity :math:`u_{\text{fuca,K}}` at cell boundary
:math:`K = \text{KCORE2} + 1`. Because of the presence of the abrupt area change
and of the blanket pellet stack and/or liquid sodium it was necessary to
write a separate momentum equation for the top cell, rather than using
the momentum equation developed to calculate the new fuel velocity at
all other locations.
The geometry used to derive the momentum equation is illustrated in :numref:`figure-15.2-3`. The momentum conservation cell is shown in :numref:`figure-15.2-3`. The
momentum conservation is written first in integral form.
.. math::
:label: 15.2-64
\Delta\left\lbrack u_{\text{fuca,K}} \cdot \rho_{\text{fufi}} \cdot A_{\text{ca,K-1}} \cdot \frac{\Delta z_{\text{K-1}}}{2} + u_{\text{fuca,K}} \cdot \frac{A_{\text{K-1}}}{A_{\text{K}}} \cdot \left( \rho_{\text{fufi}} A_{\text{ca,K}} \cdot \Delta z_{\text{KKMXPN}} \right. \right. \\ \left. \left.
+ M_{\text{slug}} \cdot \text{NRPI} \right) \right\rbrack \big/ \Delta t_{\text{PN}} = - A_{\text{bk}} \cdot \left( P_{\text{PLENUM}} - P_{\text{K-1}} \right) + \left( \rho \cdot u^{2} \cdot A \right)_{\text{fuca,K}} \\
- \rho_{\text{fufi}} u_{\text{fuca,K}} \left| u_{\text{fuca,K}} \right| \cdot A_{\text{K-1}} F_{\text{FRICTION,K-1}} \cdot \frac{\Delta z_{\text{K-1}}}{4D_{\text{ca,K-1}}} \\
- \rho_{\text{fufi}} u_{\text{fuca,K}} \left| u_{\text{fuca,K}} \right| \cdot A_{\text{K}} \left( \frac{A_{\text{K-1}}}{A_{\text{K}}} \right)^{2} \cdot F_{\text{FRICTION,K}} \cdot \frac{\Delta z_{\text{KKMXPN}}}{2 D_{\text{ca,K}}}
Note that in writing :eq:`15.2-64` the blanket stack and/or liquid sodium
was assumed to move with the same velocity as the fuel in the top cell
:math:`\text{KCORE2} + 1`. More will be said about this assumption later in this chapter.
The mass of the fuel stack and or sodium slug is defined by the input
variable :sasinp:`FUSLMA`. After dividing :eq:`15.2-64` by AXMX,
multiplying by the number of pins NRPI, and using the definition of
generalized smear densities and volume fractions we obtain:
.. math::
:label: 15.2-65
\frac{\Delta z_{\text{K-1}}}{2} \cdot \Delta \left\lbrack u_{\text{fufi,K}} \cdot {\rho'}_{\text{fufi,K-1}} \right\rbrack \big/ \Delta t_{\text{PN}} + \frac{\theta_{\text{ca,K-1}}}{\theta_{\text{ca,K}}} \cdot \Delta\left\lbrack u_{\text{fuca,K}} \left( {\rho'}_{\text{fufi,K}} \cdot \Delta z_{\text{KKMXPN}}^{\text{REF}} \right. \right. \\ \left. \left.
+ \frac{M_{\text{slug}} \cdot \text{NRPI}}{\text{AXMX}} \right) \right\rbrack \big/ \Delta t_{\text{PN}} = - \theta_{\text{bk}} \cdot \left( P_{\text{PLENUM}} - P_{\text{K-1}} \right) + \left( \rho' u^{2} \right)_{\text{fufi,K}} \\
- {\rho'}_{\text{fufi,K-1}} \cdot u_{\text{fuca,K}} \cdot \left| u_{\text{fuca,K}} \right| \cdot F_{\text{FRICTION,K-1}} \cdot \frac{\Delta z_{\text{K-1}}}{4 \cdot D_{\text{ca,K-1}}} \\
- {\rho'}_{\text{fufi,K}} \cdot u_{\text{fuca,K}} \cdot \left| u_{\text{fuca,K}} \right| \cdot \left( \frac{\theta_{K-1}}{\theta_{K}} \right)^{2} \cdot F_{\text{FRICTION,K}} \cdot \frac{\Delta z_{\text{KKMXPN}}}{2 D_{\text{ca,K}}}
The convective momentum flux
:math:`\left( {\rho'} u^{2} \right)_{\text{fufi,K}}` is calculated using
:eq:`15.2-50` and the area fraction :math:`\theta_{\text{bk}}` at the abrupt area
change is calculated as follows:
.. math::
:label: 15.2-66
\theta_{\text{bk}} = \begin{cases}
\theta_{\text{K}} & \text{if } u_{\text{K}} \geq 0 \text{ i.e. expansion} \\
1.67 \frac{\theta_{\text{K}} \cdot \theta_{\text{K-1}}}{\theta_{\text{K}} + \theta_{\text{K-1}}} & \text{if } u_{\text{K}} < 0 \text{ i.e. contraction}
\end{cases}
The friction factor :math:`F_{\text{FRICTION}}` is defined by :eq:`15.2-39` and
15.2-39a. By using the identity:
.. math::
:label: 15.2-67
\Delta \left( u \rho' \right) = {\rho'}^{n + 1} \cdot \Delta u + u^{n} \cdot \Delta \rho'
:eq:`15.2-65` is now rearranged in the final form:
.. math::
:label: 15.2-68
\frac{\Delta u_{\text{fufi,K}}}{\Delta t_{\text{PN}}} \cdot {\rho'}_{\text{fufi,K-1}} \cdot \frac{\Delta z_{\text{K-1}}}{2} + {\rho'}_{\text{fufi,K}} \cdot \Delta z_{\text{KKMXPN}}^{\text{REF}} \cdot \frac{\theta_{\text{ca,K-1}}}{\theta_{\text{ca,K}}} + \frac{\theta_{\text{ca,K-1}}}{\theta_{\text{ca,K}}} \\
\cdot \frac{M_{\text{slug}} \cdot \text{NRPI}}{\text{AXMX}} \cdot C_{\text{TOP}} = - \theta_{\text{bk}} \cdot P_{\text{PLENUM}} \cdot C_{\text{TOP}} + P_{\text{K}} \left( 1 - C_{\text{TOP}} \right) - P_{\text{K-1}} \\
+ \left( \rho' u^{2} \right)_{\text{fufi,K}} - {\rho'}_{\text{fufi,K-1}} \cdot u_{\text{fuca,K}} \cdot \left| u_{\text{fuca,K}} \right| \cdot F_{\text{FRICTION,K-1}} \\
\cdot \frac{\Delta z_{\text{K-1}}}{4 \cdot D_{\text{ca,K-1}}} - {\rho'}_{\text{fufi,K}} \cdot u_{\text{fuca,K}} \cdot \left| u_{\text{fuca,K}} \right| \cdot \left( \frac{\theta_{\text{K-1}}}{\theta_{\text{K}}} \right)^{2} \cdot F_{\text{FRICTION,K}} \cdot \frac{\Delta z_{\text{KKMXPN}}}{2 D_{\text{ca,K}}} \\
:eq:`15.2-68` is used to calculate the change in velocity
:math:`\Delta u_{\text{fufi,K}}` at the top cavity boundary and thus determines the
new fuel ejection velocity. Note that :eq:`15.2-68` contains the variable
CTOP, which was not present in :eq:`15.2-66`. This variable has
initially the value 1, so the :eq:`15.2-68` is identical to :eq:`15.2-67`.
However, when the solid pellet stack reaches a rigid obstacle, i.e.,
when :math:`\Delta z_{\text{KKMXPN}} =` :sasinp:`FUSLDT`, the variable
:math:`C_{\text{TOP}}` is set to zero.
This has the effect of decoupling the slug from molten fuel motion. The
slug is assumed to remain fixed in place. It cannot move upwards because
of the rigid obstacle, such as dimples, and presumably it cannot move
downwards because of fuel freezing and crust formation in the space
underneath. When :math:`C_{\text{TOP}} = 0` the pressure difference used in the
momentum equation thus become :math:`P_{\text{K}} - P_{\text{K-1}}`
rather than the :math:`P_{\text{PLENUM}} - P_{\text{K-1}}` used
previously.
.. _section-15.2.6:
Time-step Determination for the Pre-Failure In-Pin Motion
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The PINACLE time step :math:`\Delta t_{\text{PN}}` used in the numerical solution
of all the in-pin and channel conservation equations is restricted by
the sonic Courant condition for the in-pin flow. The determination of
the PINACLE time step, :math:`\Delta t_{\text{PN}}`, is given at the end of the
channel flow description in :numref:`section-15.3.1.3`. In the present section,
only the restriction imposed by the sonic Courant condition for the
in-pin flow is described.
The time-step for the in-pin motion is computed to be a fraction, 0.4,
of the minimum time step based on the sonic Courant condition
.. math::
:label: 15.2-69
\Delta t_{\text{PN}} = 0.4 \cdot \min \left\lbrack \Delta z_{\text{K}}/\left( v_{\text{sonic,K}} \
+ \left| u_{\text{fuca,K}} \right| \right) \right\rbrack_{\text{K=KK1,KKMX}}
The minimum in :eq:`15.2-69` is evaluated over all the axial cells of the
molten fuel cavity. The sonic velocity, :math:`v_{\text{sonic}}`, is calculated
from an expression for an adiabatic, homogeneous two-component
gas-liquid mixture which is based on Eq. 27 in Reference [15-8].
.. math::
:label: 15.2-70
v_{\text{sonic}}^{2} = \frac{\gamma_{\text{fi}}P_{\text{fica}}}{\left\lbrack \alpha_{\text{fica}}^{2} \cdot \rho_{\text{fica}} + \alpha_{\text{vica}}\left( 1 - \alpha_{\text{fica}} \right) \rho_{\text{fuca}} \right\rbrack} \\
+ \left\lbrack \left( 1 - \alpha_{\text{fica}} \right)^{2} \rho_{\text{fuca}} + \alpha_{\text{fica}} \left( 1 - \alpha_{\text{vica}} \right) \rho_{\text{fica}} \right\rbrack \gamma_{\text{fi}} P_{\text{fica}} K_{\text{fu}} \\
where
:math:`\alpha_{\text{fica}} = \frac{\theta_{\text{fica}}}{\theta_{\text{ca}}}` = Void fraction in
the cavity
:math:`\gamma_{\text{fi}} = \frac{C_{\text{p,fi}}}{C_{\text{v,fi}}} = 1.4` (value assumed in PINACLE)
:math:`K_{\text{fu}} = \text{CMFU}` = Input liquid fuel compressibility
The above equation holds for adiabatic gas behavior, although the in-pin
fission-gas treatment in PINACLE is isothermal (the gas temperature is
assumed to be always equal to the fuel temperature). However, the sonic
velocity for adiabatic gas behavior is higher than that for isothermal
gas behavior, and thus, leads to a more conservative (i.e., smaller)
time step. Moreover, if a pure gas flow were treated in sections of pins
with a prefabricated central hole, the current time-step determination
would actually be necessary.