.. _section-14.2:
In‑Pin Fuel Motion
------------------
.. _section-14.2.1:
Overview and Assumptions
~~~~~~~~~~~~~~~~~~~~~~~~
When LMFBR fuel pins melt in an overpower accident, the interiors of the
pins melt first and form cavities containing molten fuel and fission
gas. Some of this fission gas is dissolved in the molten fuel and the
remainder is free fission gas which resides in bubbles too large to be
constrained by surface tension. The decrease in the void volume, caused
by the density decrease of the melting fuel and also of the molten fuel,
which is further heated, compresses the free fission gas and causes a
hydrostatic cavity pressurization which loads the cladding. This cavity
pressure increases further when more free gas becomes available due to
additional fuel melting and release of dissolved gas, and when fuel
vaporization occurs. This leads eventually to cladding failure. The
subsequent fuel and fission‑gas ejection from the fuel pin and the local
depressurization in the cavity caused by this ejection leads to fuel
motion inside the pin toward the failure location. If the pin failure
location for driver fuel is above the core midplane, this in‑pin fuel
motion reduces reactivity; if it is near the midplane, it can cause a
reactivity addition.
The concept of a pressurized molten pin cavity, which was originally
developed for TOP accidents, is also reasonable for LOF accidents which
lead to a power rise and fuel melting before the cladding melts.
Therefore, this cavity concept is used not only in PLUTO2 but also in
LEVITATE. However, in an LOF accident, leading only to a few times
nominal power or less, cladding motion and fuel swelling will precede
fuel melting. In this case, the fuel pin is likely to disrupt completely
into a mixture of molten fuel, gas, and solid fuel chunks once fuel
melting has begun. This is only modeled in LEVITATE.
The in‑pin fuel motion in the molten pin cavity in both PLUTO2 and
LEVITATE is treated as a one‑dimensional, compressible flow with a time
dependent and spatially variable flow cross section. The fuel/fission
gas flow is considered to be homogeneous and the fission gas is assumed
to be in thermal equilibrium with the fuel.
Special attention had to be given to the potentially large mass sinks
due to fuel and fission‑gas ejection from the molten pin cavity. Mass
sources due to fuel melt‑in which cause flow cross section changes are
also considered. The fuel which is melting into the cavity brings with
it dissolved (intragranular) and free (intergranular) gas. These two
types of gas are treated separately in the PLUTO2 and LEVITATE in‑pin
motion models. The free gas is considered not to be constrained by
surface tension and to act as an ideal gas that is at the local fuel
temperature. The dissolved gas, which is in the form of small bubbles,
is assumed to be strongly constrained by surface tension, and its
contribution to the local pressure depends on the input variable PRSFTN.
If PRSFTN ≤ 10\ :sup:`7` then the model assumes that the volume of the
small bubbles containing the dissolved gas is negligible and thus the
dissolved gas does not immediately affect the local pressure. If PRSFTN
>10\ :sup:`7`, the pressure calculation takes into account the volume of
the dissolved gas bubbles, as described in :numref:`section-14.2.6`. The dissolved
gas is released (coalesces) according to a decay type law using an input
decay constant. The amount released during each time step is added to
the free gas. The fission-gas pressure calculation takes the
compressibility of the liquid fuel into account and reduces to a
liquid‑fuel single‑phase pressure calculation for no fission gas and
zero void volume. The total pressure includes the fission‑gas pressure
and fuel vapor pressure. The fuel vapor pressure is based on the
radially averaged fuel temperature.
Incoherency of the pin failures is mainly an issue in lower ramp rate
TOP accidents. In these accidents, pin failure incoherency should be
helpful for the post‑accident coolability. If not all pins fail in the
lead assemblies during a slow TOP, some subchannels may remain open and
provide coolant paths. However, if near‑midplane failures are assumed in
slow TOPs in higher void worth cores, there may be a potential for
exacerbating the accident due to pin failure incoherency. This is
because the first pin failures will probably lead to some rapid fuel
sweepout but may also cause significant voiding. Pins failing later will
inject fuel into partially voided regions and this fuel may not
experience rapid sweepout.
There are two major reasons for the incoherency in pin failures. One is
the stochastic nature of pin failures and the other is due to
differences in the power, flow, or coolant temperatures experienced by
different pins. For example, the outer two fuel‑pin rows in a
subassembly see colder coolant temperatures due to the proximity of the
hexcan wall. Many subassemblies in an LMFBR also see considerable radial
power skews. In PLUTO2 an attempt has been made to treat the stochastic
pin failures by allowing different pin failure groups. However,
currently only one pin group can fail, and the other two groups of pins
have to remain unfailed; (see input variables NRPIl, NRPI2, and NRPI3).
The treatment of the delayed fuel expulsion and in‑pin fuel motion for
the two additional pin groups is not yet operational.
.. _section-14.2.2:
Initial Conditions for the In‑pin Calculation from DEFORM and the Pre‑Failure Pin Heat‑transfer Calculation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Before pin failure, the molten fuel/fission‑gas cavity is treated by the
DEFORM pin behavior module. In the DEFORM treatment, the initially
present central fuel hole is considered and the closing of this central
hole due to the expansion of the molten fuel is calculated. Although
DEFORM does not calculate axial fuel redistribution inside the molten
cavity, an axially uniform pressure is assumed for the entire molten
cavity and found by considering the total gas and volume available. (An
option is also available in DEFORM to use the nodal pressures for
loading the cladding. This option is relevant only for very high power
situations in which the molten cavity grows to a large size within a few
milliseconds.) When pin failure is predicted by DEFORM or because an
input melt‑fraction or failure time have been reached, PLUTO2 is
initialized from DEFORM conditions. These include the uniform cavity
pressure (or the axial pressure distribution if the special DEFORM
option was used), the dissolved and free fission‑gas densities, and the
molten cavity geometry. However, PLUTO2 will have exactly the same
cavity geometry as DEFORM only if the PLUTO2 input option *FNMELT* = 0.0
is chosen. This value means that the radial cavity boundary is defined
by the location of the fuel solidus temperature (an assumption in
DEFORM). However, it is recommended for PLUTO2 to use a higher melt
fraction *(FNMELT* > 0.5) for the definition of the location of the
cavity boundary. This is because the analysis of TREAT test E8 [14‑5,
14‑6] showed that fuel, which had not yet undergone significant melting,
did not move. This can probably be explained by a high viscosity of
partially molten fuel.
When the regular DEFORM option of a uniform cavity pressure is used,
there is an inconsistency in the DEFORM/PLUTO2 transition due to the
lack of an axial in‑pin fuel motion model in DEFORM. The nodal pressures
calculated in PLUTO2 which use the pressure‑generating free fission gas
and the nodal volume obtained from DEFORM can be different from the
average cavity pressure. Since the latter is considered more realistic
for all but very high power situations, the free fission‑gas content of
the PLUTO2 nodes is adjusted to give the DEFORM calculated average
cavity pressure. When the axially nonuniform pressure option in DEFORM
is chosen (only relevant for very high power situations), the
PLUTO2‑calculated nodal pressure can still be somewhat different from
the one calculated by DEFORM and has to be slightly adjusted. This is
because DEFORM takes into account the radial temperature, fuel density,
and porosity profiles in the molten fuel cavity, whereas PLUTO2 uses a
radially averaged temperature, fuel density and porosity in the cavity
nodes. The above‑mentioned inconsistencies in the transition from DEFORM
to PLUTO2 are not considered to be serious. To remove them would require
the inclusion of pre‑failure in‑pin fuel motion in DEFORM and the
accounting of the radial temperature profile in the molten cavity in
PLUTO2 or LEVITATE.
The initialization of the PLUTO2 and LEVITATE cavity is done in
subroutines PLINPT (**PL** UTO2 **IN**\ PU\ **T**) and PLSET (**PL**\ UTO2
**SET**\ UP). DEFORM calculated cavity dimensions and free and dissolved
fission‑gas densities are transferred to PLUTO2 in these routines. Also,
the radially averaged cavity temperatures are determined from the
temperature profiles obtained from the heat‑transfer routines, TSHTRN or
TSHTRV.
.. _section-14.2.3:
Coupling of the In‑pin Motion Calculation with the PLHTR Heat Transfer Calculation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
For the calculation of the growth of the molten fuel cavity, it is
important that a fuel temperature calculation is performed in the solid
fuel surrounding the molten fuel cavity. This temperature calculation is
performed with the PLHTR (PLUTO2 HEAT TRANSFER) routine which is a
modified version of the transient heat‑transfer subroutine TSHTRV. (The
latter is used to calculate fuel‑pin temperatures during the sodium
boiling phase of a LOF accident). In the PLHTR heat conduction
calculation, the fuel node adjacent to the molten cavity boundary
receives heat from the molten cavity in the form of a heat source. This
heat source is the total energy convected to the cavity wall from the
molten fuel flow during a heat‑transfer time step divided by the
heat‑transfer time‑step size. If a solid fuel node reaches the input
melt fraction criterion FNMELT, this node will be gradually added to the
molten cavity. The rate of fuel addition depends on the temperature
gradient between the node adjacent to the cavity wall and its
neighboring solid fuel node.
.. _section-14.2.4:
Overview of the Numerical Approach for the In‑pin Fuel Motion Calculation and Description of Subroutines PLlPIN and PL2PIN
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
There are several requirements for the solution algorithm for this one
dimensional, compressible flow problem with variable flow cross
section: (a) it had to be able to handle large mass sinks (due to fuel
ejection from the pins), (b) it had to conserve mass perfectly, and (c)
it had to run efficiently. This was achieved with a predominantly
explicit Eulerian solution method. All convection 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 is of
importance for treating the strong mass sinks.
For the in‑pin motion calculation three mass, one energy and one
momentum conservation equations are solved. The mass conservation
equations are for the fuel and the free and the 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. In the fuel‑ejecting
nodes, this new pressure and the advanced densities and energies are
used for calculating the fuel and fission‑gas ejection rates. The
pressure in the ejection nodes is then decreased in order to account for
the fuel and gas losses. The adjusted new pressures are then used in the
fuel/fission‑gas momentum equation. There is an input option in PLUTO2
to use a combination of the new and the old pressures. This option lets
the pressure, :math:`P`, be:
.. math::
:label: 14.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 [14‑21].
In these test calculations, it was also attempted to iterate on the
convective fluxes by repeating most of the in‑pin calculation for each
time step. This seemed to make the results somewhat smoother, but did
not allow a time step as long as the one allowable with the
above‑described explicit scheme.
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 with the densities, energies, and
temperatures defined at the cell centers and the velocities at the cell
edges is used. The spatial differencing uses full donor cell
differencing. Although this is not as accurate as higher order
differencing, it makes the calculation stable because it introduces a
numerical diffusion effect. It is interesting to note that this
stabilizing effect was not enough to keep two‑phase test calculations
stable for shock reflections from a rigid wall. Although these test
calculations could be stabilized by introducing an artificial viscous
pressure into the calculations [14‑21], the artificial viscous pressure
is not necessary for regular pin blowdown calculations. It could even
cause problems if the input constants, which affect the magnitude of the
artificial viscous pressure, are set too high. The two relevant input
parameters, *ClVIPR* and *C2VIPR*, should therefore be set to small
values such as :math:`10^{-3}`.
The free fission‑gas mass conservation, the fuel mass conservation, and
the fuel energy equation are solved in subroutine PLlPIN (*PL*\ UTO2
*1*\ st *PIN* ROUTINE). PLlPIN also computes the molten cavity geometry
changes due to fuel melt‑in, and the fuel and fission‑gas ejection from
the pins (see :numref:`section-14.3`). The fuel/fission‑gas momentum equation and
the dissolved fission‑gas mass equation are solved in PL2PIN (*PL*\ UTO2
*2*\ nd *PIN* ROUTINE). This routine also calculates the sonic
velocities for each node. The minimum time step found is the predicted
time step for the next cycle. However, the actual PLUTO2 time step,
which is used both in the pin and in the channel calculation, is the
smaller of the pin and channel hydrodynamics time steps. The channel
time step is usually the smaller one, and thus, dominates the time‑step
selection.
.. _section-14.2.5:
Definition of the Generalized Smear Densities for the In‑pin Calculation
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The use of generalized smear densities in PLUTO2 (and LEVITATE) has been
prompted by the many different moving and stationary components in this
problem. Its use also simplifies the differential and finite difference
equations for variable cross section flow. The pie chart in :numref:`figure-14.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: 14.2-2
\theta_{\text{k}}\left( z,t \right) = \frac{A_{\text{k}}\left( z,t \right)}{\text{AXMX}}
.. _figure-14.2-1:
.. figure:: media/image9.png
:align: center
:figclass: align-center
The Left-hand Side of This Pie Chart Illustrates the Possible Material Cross Sectional Areas in the Fuel Pin, the Whole Pie Representing an Area AXMW Which is an Input Parameter
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 PLUTO2 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: 14.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 cavities of the failed pins.
:math:`\theta_{\text{ca}}` = generalized volume fraction of the molten cavities in
all failed pins which can be calculated from :math:`A_{\text{ca}} \big/ \text{AXMX}`
:math:`\theta_{\text{fuca}}` = generalized volume fraction of the fuel in the
cavities of all failed pins
Generalized smear densities, which are always marked by a prime, are
defined as products of physical densities and generalized volume
fractions:
.. math::
:label: 14.2-4
{\rho'}_{\text{fuca}} \left( z,t \right) = \rho_{\text{fuca}}\left( T \right)\theta_{\text{fuca}}\left( z,t \right) = \frac{\rho_{\text{fuca}}\left( T \right) A_{\text{fuca}}\left( z,t \right)}{\text{AXMX}}
.. math::
:label: 14.2-5
{\rho'}_{\text{fica}} \left( z,t \right) = \rho_{\text{fica}}\left( T,P \right)\theta_{\text{fica}}\left( z,t \right) = \frac{\rho_{\text{fica}}\left( T,P \right) A_{\text{fica}}\left( z,t \right)}{\text{AXMX}}
.. math::
:label: 14.2-6
{\rho'}_{\text{fsca}} \left( z,t \right) = \rho_{\text{fsca}}\left( z,t \right)\theta_{\text{fuca}}\left( z,t \right) = \frac{\rho_{\text{fsca}}\left( z,t \right) A_{\text{fuca}}\left( z,t \right)}{\text{AXMX}}
where the subscript :math:`\text{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 :math:`\text{fuca}`. It should again be
pointed out that the A's refer to total cross section areas of all the
cavity fuel, free fission gas, and dissolved fission gas in the failed
pins of one subassembly.
The generalized source or sink term is written as:
.. math::
:label: 14.2-7
S' = \frac{S^{l}}{\text{AXMX}}
where the source or sink term :math:`S^{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 :math:`m^{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:`14.2-7`, a change in the generalized smear
density :math:`\rho'` due to source or sink :math:`S'` it just :math:`S' \Delta t`.
.. _section-14.2.6:
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) and one energy and one momentum conservation equation. The
continuity equation for the molten fuel in the pin cavity is written:
.. math::
:label: 14.2-8
\frac{\partial}{\partial \text{t}}(\rho_{\text{fuca}}A_{\text{fuca}}) = - \frac{\partial}{\partial \text{z}}\left( \rho_{\text{fuca}}A_{\text{fuca}} u_{\text{fuca}} \right) + S_{\text{fuca},\text{me}}^{l}\left( z,t \right) - S_{\text{fuca},\text{ej}}^{l}\left( z,t \right)
where the subscripts :math:`\text{me}` and :math:`\text{ej}` refer to fuel melting into and fuel
ejection from the pin cavities of all failed pins, respectively. By
dividing by AXMX and using the definitions of the generalized smear
densities and source and sink terms, :eq:`14.2-8` becomes:
.. math::
:label: 14.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},\text{me}} - {S'}_{\text{fuca},\text{ej}}
where the primed sources and sinks are per unit time and per unit smear
volume (see :eq:`14.2-7`). The integrated source term,
:math:`S'_{\text{fuca,me}} \Delta t_{\text{PL}}`, which is actually needed
in the finite difference equations of the code, is calculated from
.. math::
:label: 14.2-10a
{S'}_{\text{fuca},\text{me}}\Delta t_{\text{PL}} = \rho_{\text{fu},\text{cabd}}\Delta A_{\text{me}} \text{NRPI} \cdot \frac{\text{FNPI}}{\text{AXMX}}
here
:math:`\Delta t_{\text{PL}}` = PLUTO2 time step
:math:`\rho_{\text{fu,cabd}}` = fuel density (including porosity) adjacent to the
cavity boundary,
:math:`\Delta A_{\text{me}}` = area of fuel (including porosity) melted into the
cavity per PLUTO2 time step per pin
:math:`\text{NRPI}` = number of pins per subassembly,
:math:`\text{FNPI}` = fraction of the pins which are failed.
When the solid fuel node adjacent to the cavity has not yet exceeded an
input melt fraction value FNMELT:
.. math::
:label: 14.2-10b
{S'}_{\text{fuca,me}} \Delta t_{\text{PL}} = 0
The :math:`\Delta A_{\text{me}}` in :eq:`14.2-10a` is related to a change in the
cavity diameter by
.. math::
:label: 14.2-11a
\Delta A_{\text{me}} = \frac{\pi}{4}\Delta D_{\text{ca}}\left( 2D_{\text{ca}} + \Delta D_{\text{ca}} \right)
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
:math:`\text{FNMELT} + 0.1`. Once the melt fraction has become greater
than FNMELT, the change in cavity diameter is calculated from
.. math::
:label: 14.2-11b
\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 that has melted into the
cavity per PLUTO2 time step.
.. math::
:label: 14.2-12
\text{FNMECA} = \frac{\Delta t_{\text{PL}}\left( \frac{\left( T_{\text{cabd}}^{n + 1} - T_{\text{cabd}}^{n} \right)}{\Delta t_{\text{Ht}}} \right)}{\left( 0.1 \cdot \left( T_{\text{fu,liq}} - T_{\text{fu,sol}} \right) \right)} \
= \frac{\Delta T_{\text{PL}}}{\left( 0.1 \cdot \left( T_{\text{fu,liq}} - T_{\text{fu,sol}} \right) \right)}
where
:math:`T_{\text{cabd}}^{n}, T_{\text{cabd}}^{n + 1}` are 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},\text{liq}}, T_{\text{fu},\text{sol}}` are fuel
liquidus and solidus temperatures, respectively. The difference between
the two is the melting band width.
:math:`\Delta t_{\text{PL}}` is the PLUTO2 time step.
:math:`\Delta T_{\text{PL}}` is the temperature change of the fuel adjacent to the
cavity during a PLUTO2 time step.
:eq:`14.2-12` 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 sink term :math:`{S'}_{\text{fuca,ej}}` that is due to fuel ejection is also
the source term for the coolant channel equations. It will be described
in detail in :numref:`section-14.3`. Its integrated form is similar to
:eq:`14.2-10a`
.. math::
:label: 14.2-13
{S'}_{\text{fuca},\text{ej}} \cdot \Delta t_{\text{PL}} = {\rho'}_{\text{fuca}} \cdot \text{FN}
where
:math:`\text{FN}` = fraction of the fuel in an ejecting cavity cell that is actually
ejected during a PLUTO2 time step.
The free fission-gas mass conservation equation is written
.. math::
:label: 14.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},\text{me}}^{l}\left( z,t \right) - S_{\text{fica},\text{ej}}^{l} + S_{\text{fsca,rl}}^{l}
where the subscripts :math:`\text{fsca}` and :math:`r l` 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 and sink terms yields
.. math::
:label: 14.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},\text{me}}^{l} - S_{\text{fica,ej}}^{l} + S_{\text{fsca,rl}}^{l}
The integrated source term for free fission gas due to fuel melt-in is
similar to that in :eq:`14.2-10a`:
.. math::
:label: 14.2-16
{S'}_{\text{fica,me}} \Delta t_{\text{PL}} = \rho_{\text{fifs,cabd}} \Delta A_{\text{me}} \left( \text{NRPI} \cdot \frac{\text{FNPI}}{\text{AXMX}} \right) \cdot \text{FIFNGB}
where :math:`\Delta A_{\text{me}}` is the calculated change in cavity cross
sectional area described before (see :eq:`14.2-11a`),
:math:`\rho_{\text{fifs,cabd}}` is the density of all the gas in the fuel node
adjacent to the cavity, and FIFNGB is an input parameter determining
the fraction of all the gas entering the cavity with the molten fuel
which instantaneously becomes free gas. This is an input parameter
because DEFORM does not currently discriminate between intra-granular
and grain-boundary gas. The ratio of the grain-boundary gas to the total
gas in a node may be similar to the fraction of the gas that becomes
free upon melting into the cavity. However, it seems more likely that
only a smaller fraction will be instantaneously available because some
of the grain-boundary bubbles may not be large enough to be considered
free gas as soon as the fuel melts.
The quantity :math:`\rho_{\text{fifs,cabd}}` is calculated from:
.. math::
:label: 14.2-16a
\rho_{\text{fifs,cabd}} = \text{RETFG2}_{\text{cabd}} \frac{\left( \frac{\text{FUMS}_{\text{cabd}}}{\text{FUELM}} S_{\text{cabd}} \right)}{\text{VOLUME}_{\text{cabd}}}
where
:math:`\text{RETFG2}_{\text{cabd}}` is the total retained 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}}` is is the current fuel mass in the melting radial
fuel-pin node at the cavity boundary,
:math:`\text{FUELMS}_{\text{cabd}}` is 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}}` is the current volume of the melting radial
fuel-pin node at the cavity boundary.
The free fission-gas sink term due to fuel ejection in :eq:`14.2-15` is
described in detail in :numref:`section-14.3`. Its form is similar to that of the
sink term for fuel ejection :eq:`14.2-13`:
.. math::
:label: 14.2-17
{S'}_{\text{fica,ej}} \cdot \Delta t_{\text{PL}} = {\rho'}_{\text{fica}} \cdot \text{FN}
The term :math:`\text{FN}` used in this equation is the same as in :eq:`14.2-13`
because it is assumed in PLUTO2 that fission gas and fuel are
ejected with the same volume fractions that exist in the ejecting cavity
node.
The source term due to dissolved fission-gas release in :eq:`14.2-15` is:
.. math::
:label: 14.2-18
{S'}_{\text{fsca},\text{r}l} = {\rho'}_{\text{fsca}} \cdot \text{CIRTFS}
where CIRTFS is a release constant for dissolved fission gas which is
input and has the dimensions :math:`s^{-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: 14.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}}^{l} - S_{\text{fsca,ej}}^{l} - S_{\text{fsca,rl}}^{l}
By dividing this equation by :math:`\text{AXMX}` and using the definitions of the
generalized smear density and sources and sinks, one obtains:
.. math::
:label: 14.2-20
\frac{\partial}{\partial \text{t}}{\rho'}_{\text{fsca}} = - \frac{\partial}{\partial \text{z}}\left( {\rho'}_{\text{fica}} u_{\text{fuca}} \right) + S_{\text{fsca,me}}^{l} - S_{\text{fsca,ej}}^{l} - S_{\text{fsca,rl}}^{l}
where
.. math::
:label: 14.2-21
{S'}_{\text{fsca,me}} = {S'}_{\text{fica,me}} \frac{\left( 1 - \text{FIFNGB} \right)}{\text{FIFNGB}}
and :math:`S'_{\text{fica,me}}` and FIFNGB are defined for :eq:`14.2-16`. The
ejection term is connected to the one for the free fission-gas ejection,
which is given in :eq:`14.2-17`:
.. math::
:label: 14.2-22
{S'}_{\text{fsca,ej}} = {S'}_{\text{fica,ej}}\frac{{\rho'}_{\text{fsca}}}{{\rho'}_{\text{fica}}}
The absolute value of the sink term due to the dissolved fission-gas
release has been described in :eq:`14.2-18`.
The fuel energy conservation equation is written:
.. math::
:label: 14.2-23
\frac{\partial}{\partial \text{t}}\left( \rho_{\text{fuca}} e_{\text{fuca}} A_{\text{fuca}} \right) = & - \frac{\partial}{\partial \text{z}}\left( \rho_{\text{fuca}} e_{\text{fuca}} A_{\text{fuca}} u_{\text{fuca}} \right) + S_{\text{fuca,me}}^{l} e_{\text{fu,cabd}} \\
& - {S'}_{\text{fuca,ej}} e_{\text{fuca}} + 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} \cdot \text{FNPI}
Dividing :eq:`14.2-23` by AXMX and using the definitions of the
generalized smear densities and sources and sinks produces:
.. math::
:label: 14.2-24
\frac{\partial}{\partial \text{t}}\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} \right) = & - \frac{\partial}{\partial \text{z}}\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right) + {S'}_{\text{fuca,me}} e_{\text{fu,cabd}} \\
& - S_{\text{fuca,ej}}^{l} e_{\text{fuca}} + Q {\rho'}_{\text{fuca}} \\
& - h_{\text{fuca,cabd}} \cdot \left( T_{\text{fuca}} - T_{\text{fu,cabd}} \right) \pi D_{\text{ca}} \text{NRPI} \cdot \frac{\text{FNPI}}{\text{AXMX}}
By rewriting the left-hand side of :eq:`14.2-24` as
.. math::
:label: 14.2-24a
\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:`14.2-8`, one obtains
.. math::
:label: 14.2-25
{\rho'}_{\text{fuca}}\frac{\partial}{\partial \text{t}}e_{\text{fuca}} = & - \frac{\partial}{\partial \text{z}}\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,ej}} \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}} \text{NRPI} \cdot \frac{\text{FNPI}}{\text{AXMX}}
where :math:`Q` is the fission heat source per kg of fuel
.. math::
:label: 14.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:`14.2-25a`,
.. math::
:label: 14.2-25b
\text{FNPOHE} = \exp\left\lbrack \text{POWCOF}\left( 1 \right) + \Delta t\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 PLUTO2 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.
: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 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: 14.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`, :math:`K`.
The heat-transfer coefficient in :eq:`14.2-23` is the sum of a convective
and a conduction heat transfer term.
.. math:: h_{\text{fuca},\text{cabd}} = \left( h_{1} + h_{2} \right)
where
.. math::
:label: 14.2-27
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: 14.2-27a
\text{Nu} = \text{St} \cdot \text{Pr} \cdot \text{Re}
where
St = Stanton number = :math:`\frac{h}{\left( \rho u C_{\text{p}} \right)}`
Pr = Prandtl number = :math:`\mu \frac{C_{\text{p}}}{k}`
Re = Reynolds numbers = :math:`D \frac{\rho u}{\mu}`
The Deissler correlation [14-22, 14-23] was used for finding the
relationship between the three nondimensional numbers on the right-hand
side of :eq:`14.2-27a`. The Prandtl number for fuel is about 2.2. By using
this value in the Deissler correlation, it can be shown that
.. math::
:label: 14.2-28
\text{St} \approx 0.0158 \text{Re}^{- 0.2}
By using :eq:`14.2-8` and the definition of the Prandtl number,
.. math::
:label: 14.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
.. math::
:label: 14.2-30
\text{Re}_{\text{ca}} = \frac{\left( \frac{\left| u_{\text{fuca}} \right| D_{\text{ca}} {\rho'}_{\text{fuca}}}{\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 (the stagnant
condition is assumed in the simplified PLUTO2 node --- see :numref:`section-14.1.2`), is of the following form
.. math::
:label: 14.2-31
h_{2} = \frac{k_{\text{fu}}}{D_{\text{ca}}}4
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: 14.2-32
P_{\text{ca}} = P_{\text{fica}}\left( T_{\text{fuca}}, {\rho}_{\text{fica}} \right) + P_{\text{fvca}} \left( T_{\text{fuca}} \right)
where the fission-gas pressure is calculated form a special form of the
ideal-gas equation which takes the compressibility of the liquid fuel
into account:
.. math::
:label: 14.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 reasonable solution of the quadratic :eq:`14.2-33` is:
.. math::
:label: 14.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 reasonable result.
:eq:`14.2-34` reduces to a single-phase liquid pressure solution for no
fission gas and :math:`\theta_{\text{fuca}} > \theta_{\text{ca}}`:
.. math::
:label: 14.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 modified :eq:`14.2-33` in which the term with
:math:`K_{\text{fu}}` is dropped.
If the dissolved fission gas present in the molten fuel cavity is
assumed to occupy a negligible volume, in comparison with the volume
occupied by the free gas then the :math:`\theta_{\text{fuca}}` in :eq:`14.2-33` and
14.2-34 reflects the actual volume occupied by the molten fuel. This
assumption is justified when the ambient cavity pressure is
significantly smaller than the bubble surface tension pressure, as was
the case in all TREAT experiments analyzed with SAS4A. When this
assumption is used, the dissolved fission gas does not affect the cavity
pressure until it is released, in a gradual manner, to the free gas
field.
However, during a transient characterized by a rapid increase in the
power level, the cavity pressures can reach higher values, so that the
influence of the dissolved fission gas becomes important. In order to
address this issue the dissolved fission gas modeling in the pin cavity
was introduced as described below.
The dissolved fission gas is assumed to exist in the form of small
bubbles, with the gas inside the bubbles at the pressure
:math:`P_{\text{bubble}}`:
:math:`P_{\text{bubble}} = P_{\text{surface tension}} + P_{\text{cavity}}`.
The quantity :math:`P_{\text{surface tension}}` is an input to the code
(PRSFTN). Its value can be determined as:
.. math:: P_{\text{surface tension}} = \frac{2\sigma}{r},
where:
:math:`\sigma` = bubble surface tension,
:math:`r` = bubble radius.
The recommended value for :math:`P_{\text{surface tension}}` is :math:`4 \times 10^{7}` Pa,
corresponding to :math:`\sigma = 0.4` N/m (400 dyne/cm) and :math:`r = 2 \times 10^{-8}` m (200 Ã…).
The surface tension pressure is assumed to
remain constant during the calculations. Furthermore, in order to
simplify the calculations, the pressure :math:`P_{\text{cavity}}` used in the
calculation of the bubble pressure :math:`P_{\text{bubble}}` is the cavity
pressure at the end of the previous time step. The volume occupied by
the dissolved gas in the axial cavity cell i is then calculated using
the equation:
.. math:: V_{\text{D,i}} = \frac{M_{\text{D,i}} \cdot R \cdot T_{\text{i}}}{m \cdot P_{\text{bubble,i}}} = \frac{R}{m} \cdot \frac{M_{\text{D,i}} \cdot T_{\text{i}}}{P_{\text{bubble,i}}}
where:
:math:`M_{\text{D,I}}` = mass of dissolved fission gas in cell :math:`i`,
:math:`\frac{R}{m}` = universal gas constant divided by molar mass,
:math:`T_{\text{i}}` = gas temperature in cell :math:`i`, assumed to be equal to the
molten fuel temperature,
:math:`P_{\text{bubble,I}}` = the pressure of the gas contained in the small
bubbles, calculated above.
The fuel and dissolved gas are assumed to form a homologous mixture,
with the volume given by:
.. math:: V_{\text{fuel}}^{\text{mixture}} = V_{\text{fuel,i}} + V_{\text{D,i}},
where:
:math:`V_{\text{fuel,I}}` = molten fuel volume in cell :math:`i`,
:math:`V_{\text{D,i}}` = dissolved gas volume in cell :math:`i`.
The calculation then proceeds as before, using the volume
:math:`V_{\text{fuel,i}}^{\text{mixture}}` instead of
:math:`V_{\text{fuel,i}}` to determine the value of :math:`\theta_{\text{fuca}}`
used in :eq:`14.2-33` and :eq:`14.2-34`.
When there is little molten fuel in an axial cavity cell, i.e., the gas
occupies a large fraction of the cell volume, a perfect gas equation is
used to determine the free gas pressure. The underlying assumption is
that the compressibility of the fuel/dissolved-gas-mixture has a
negligible effect under these circumstances. However, if larger amounts
of molten fuel are present, the pressure calculation takes into account
the molten fuel compressibility. In this case, the dissolved fission gas
and molten fuel volumes are lumped together and assumed to have the same
compressibility as the molten fuel. This is a simplifying assumption,
but is conservative in the sense that the molten fuel/dissolved gas
mixture is likely to have a higher compressibility than the molten fuel
itself and thus would lead to somewhat lower cavity pressures during
rapid power excursions.
This treatment of the dissolved fission gas has also been implemented in
DEFORM-4, LEVITATE, and PLUTO2. If the input constant PRSFTN is less
than :math:`10^{7}` the effect of the dissolved fission gas on the ambient
pressure is ignored.
The momentum conservation equation for the fuel/fission-gas mixture is:
.. math::
:label: 14.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}} + \frac{\partial \text{P}_{\text{vi}}}{\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) \frac{A_{\text{ca}} \cdot F_{\text{friction}}}{\left( 2D_{\text{ca}} \right)} \\
- \left( S_{\text{fuca,ej}}^{l} + S_{\text{fica,ej}}^{l} \right) u_{\text{fuca}}
Dividing :eq:`14.2-36` by AXMX and using the definitions for the
generalized smear densities and mass sinks gives:
.. math::
:label: 14.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}} + \frac{\partial \text{P}_{\text{vi}}}{\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)} \\
- \left( {S'}_{\text{fuca,ej}} + {S'}_{\text{fica,ej}} \right) u_{\text{fuca}}
:math:`u_{\text{fuca}}` = vertically upward fuel velocity in the cavity (see
the sign of the gravity head term in :eq:`14.2-36`).
:math:`P_{\text{vi}}` = artificial viscous pressure for stabilizing test
problems involving shock wave propagation and shock reflection from a
rigid boundary. Not necessary for regular cavity blowdown problems.
If the net free fission-gas mass flux into a numerical node is negative
(i.e., for net outflow), then :math:`P_{\text{vi}} = 0`.
If the net convective free fission-gas mass flux into a numerical cell
is positive [14-27]:
.. math::
:label: 14.2-38
P_{\text{vi}} = \frac{\text{C2VIPR} \cdot 0.5 \cdot \left( {p'}_{\text{fu,K}} + {p'}_{\text{fi,K}} \right) \cdot \left( u_{\text{fuca,K+1}} - u_{\text{fica,K}} \right)^{2}}{\theta_{\text{ficaxx,K}}}
where
:math:`K` = index of the axial cavity node for which :math:`P_{\text{vi}}` is being
evaluated
:math:`\text{C2VIPR}` = input constant which determines the magnitude of the
numerical damping.
.. math::
:label: 14.2-38a
\theta_{\text{ficaxx}} = \begin{cases}
\theta_{\text{fica}} \text{ for } \theta_{\text{fica}} > \theta_{\text{ca}} \cdot \text{C1VIPR} \\
\theta_{\text{ca}} - \text{C1VIPR} \text{ for } \theta_{\text{fica}} < \theta_{\text{ca}} \cdot \text{C1VIPR} \\
\end{cases}
where
:math:`\text{C1VIPR}` = Input constant
The Moody friction factor :math:`F_{\text{friction}}` in :eq:`14.2-37` depends
on the Reynolds number
.. math::
:label: 14.2-39
\text{Re}_{\text{pi}} = \frac{\left| u_{\text{fuca}} \right| D_{\text{ca}} \cdot \left( {\rho'}_{\text{fuca}} + {\rho'}_{\text{fica}} \right)}{\left( \theta_{\text{cs}} \cdot \text{VIFULQ} \right)}
where VIFULQ is the viscosity of liquid fuel which is input. The
friction factor is
.. math::
:label: 14.2-39a
F_{\text{friction}} = \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}`.
The last term in :eq:`14.2-37` is necessary because the total momentum in
the cell decreases due to ejection and this has to be considered in a
momentum equation written in conservative form. 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-14.2.7:
Finite Difference Equations for the In-pin Motion
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
In the overview of the numerical scheme given in :numref:`section-14.2.4`, 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: 14.2-40
\frac{\left( {\rho'}_{\text{fuca,K}}^{n + 1} - {\rho'}_{\text{fuca,K}}^{n} \right)}{\Delta t_{\text{PL}}} = \
- \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: 14.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-14.2.4`,
which gives an overview of the numerical scheme. Densities are at the
cell centers and velocities at the cell edges (see schematic before the
momentum conservation :eq:`14.2-48`). The source and sink terms have
already been described in their finite difference form in the previous
section.
The convective fluxes at the lower and upper boundaries of the cavity
that are located in the end cells KK1 and KKMX, respectively, are
.. math::
:label: 14.2-41
\left( {\rho'} u \right)_{\text{fuca,KK1}} = 0
and
.. math::
:label: 14.2-42
\left( {\rho'} u \right)_{\text{fuca,KKMX+1}} = 0
The end cells are not always the same during a PLUTO2 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 had
to be done because of problems with overcompression of cells with a very
small cross section.
The finite difference form of the fuel energy :eq:`14.2-25` is:
.. math::
:label: 14.2-43
\frac{{\rho'}_{\text{fuca,K}} \left( e_{\text{fuca,K}}^{n + 1} - e_{\text{fuca,K}}^{n} \right)}{\Delta t_{\text{PL}}} = - \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} \cdot \text{FNPI}}{\text{AXMX}} \\
The fuel cavity temperature :math:`T_{\text{fuca}}` in the above equation has
to be calculated from the internal energy:
.. math::
:label: 14.2-44
& \text{For } e_{\text{fu,sol}} < e_{\text{fuca}} < e_{\text{fu,liq}}~ , \\
& 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::
:label: 14.2-44a
& \text{For } e_{\text{fuca}} > e_{\text{fu,liq}}~ , \\
& 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:`14.2-43` in calculated as:
.. math::
:label: 14.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 that are in cells :math:`\text{KK1}` and :math:`\text{KKMX}`, the
convective energy fluxes are zero.
.. math::
:label: 14.2-45a
\left( {\rho'}_{\text{fuca}} e_{\text{fuca}} u_{\text{fuca}} \right)_{\text{KK1}}^{n} = 0
and
.. math::
:label: 14.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:`14.2-40a` and :eq:`14.2-45`, and
the definitions of the energy gain and loss terms given earlier, and by
differencing the second term of :eq:`14.2-43` like the first, :eq:`14.2-43`
can be solved for :math:`e_{\text{fuca,K}}^{n + 1}`. Fuel temperatures
that are shown in the PLUTO2 output are calculated by using :eq:`14.2-44`
and :eq:`14.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: 14.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: 14.2-47
\theta_{\text{ca,bk}} = 0.5 \left( \theta_{\text{ca,K}} + \theta_{\text{ca,K-1}} \right)
where the subscript :math:`bk` indicates that these quantities are at the
lower boundaries of cell :math:`K`. This is shown on the schematic below. 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 the following
schematic.
.. list-table::
:header-rows: 0
:align: center
:widths: auto
* - :math:`\theta_{\text{bk}}`
- :math:`\theta_{\text{bk+1}}`
-
* - :math:`{\rho'}_{\text{bk}}`
- :math:`\rho_{\text{bk+1}}`
-
* - :math:`u_{\text{K}}`
- :math:`u_{\text{K+1}}`
-
* - :math:`\theta_{\text{K-1}}`
- :math:`\theta_{\text{K}}`
- :math:`\theta_{\text{K+1}}`
|
.. list-table::
:header-rows: 0
:align: center
:widths: auto
* - :math:`\rho_{\text{K-1}}`
- :math:`\rho_{\text{K}}`
- :math:`\rho_{\text{K+1}}`
* - :math:`P_{\text{K-1}}`
- :math:`P_{\text{K}}`
- :math:`P_{\text{K+1}}`
* - :math:`S_{\text{K-1}}`
- :math:`S'_{\text{K}}`
- :math:`S'_{\text{K+1}}`
|
.. list-table::
:header-rows: 0
:align: center
:widths: auto
* - :math:`z_{\text{K}}`
- :math:`z_{\text{K+1}}`
The finite difference form of the momentum conservation :eq:`14.2-37`, is
written using the definitions 14.2-46 and 14.2-47 as:
.. math::
:label: 14.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{PL}}} = - \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) + \left( P_{\text{vi,K}}^{n} - P_{\text{vi,K-1}}^{n} \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} F_{\text{friction,bk}} \big/ \left( 2D_{\text{ca}} \right) - \left( {S'}_{\text{fuca,ej,K-1}}^{n + 1/2} \right. \\ \left.
+ {S'}_{\text{fuca,ej,K}}^{n + 1/2} + {S'}_{\text{fica,ej,K-1}}^{n + 1/2} + {S'}_{\text{fica,ej,K}}^{n + 1/2} \right) \cdot \left( u_{\text{fuca,K}}^{n + 1} + u_{\text{fuca,K}}^{n} \right) \cdot 0.25
where
:math:`\varepsilon` = Input value EPCH that can be between zero and one (see
:eq:`14.2-1`),
and :math:`\Delta z` implies :math:`0.5 \left( \Delta z_{\text{K-1}} + \Delta z_{\text{K}} \right)`.
By defining
.. math::
:label: 14.2-49
{S'}_{\text{fuca,ej,bk}}^{n + 1/2} = \left( {S'}_{\text{fica,ej,K-1}}^{n + 1/2} + {S'}_{\text{fica,ej,K}}^{n + 1/2} + {S'}_{\text{fica,ej,K-1}}^{n + 1/2} + {S'}_{\text{fica,ej,K}}^{n + 1/2} \right) \cdot 0.5
and by collecting all terms with :math:`u_{\text{uica}}^{n + 1}` on the
left-hand side of the equation, one obtains
.. math::
:label: 14.2-50
u_{\text{fuca,K}}^{n + 1} \left\lbrack \frac{{\rho'}_{\text{fufi,bk}}^{n + 1}}{\Delta t_{\text{PL}}} + \left| u_{\text{fuca,K}}^{n} \right| {\rho'}_{\text{fufi,bk}}^{n} \frac{F_{\text{friction,bk}}}{\left( 2D_{\text{ca}} \right)} \right. \\ \left.
+ {S'}_{\text{fufi,ej,bk}}^{n + 1/2} \cdot 0.5 \right\rbrack = \frac{{\rho'}_{\text{fufi,bk}}^{n} u_{\text{fuca,K}}^{n}}{\Delta t_{\text{PL}}} \\
- \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 \big/ \Delta z - g {\rho'}_{\text{fufi,bk}} \\
- {S'}_{\text{fuca,ej,bk}}^{n + 1/2} u_{\text{fuca,K}}^{n} \cdot 0.5 \\
The convective momentum flux in :eq:`14.2-50` is calculated as
.. math::
:label: 14.2-51
\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: 14.2-51a
\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: 14.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: 14.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,KK+1}}^{2} & \text{if } u_{\text{fuca,KK1+1}} < 0
\end{cases}
.. math::
:label: 14.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:`14.2-50` can be solved for
:math:`u_{\text{fuca,K}}^{n + 1}` if :eq:`14.2-46`, :eq:`14.2-47`, and :eq:`14.2-51`
through :eq:`14.2-54` are used.
.. _section-14.2.8:
Time-step Determination for the In-Pin Motion
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The PLUTO2 time step :math:`\Delta t_{\text{PL}}` used in the numerical solution of
all the in-pin and channel conservation equations is restricted by the
sonic Courant conditions for both the in-pin and channel flows. The
sonic Courant condition for the channel flow and the determination of
the PLUTO2 time step :math:`\Delta t_{\text{PL}}` is given at the end of the
channel flow description in :numref:`section-14.4.6.4`. In the present section,
only the restriction imposed by the sonic Courant condition for the
in-pin flow is described.
The time-step size 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: 14.2-55
\Delta t_{\text{PL,pin}} = 0.4 \cdot \min\left\lbrack \frac{{z}_{\text{K}}}{v_{\text{sonic,K}} + |u_{\text{fuca,K}}|} \right\rbrack_{\text{K = KK1,KKMX}}
The minimum in :eq:`14.2-55` is evaluated over all the axial cells of the
molten fuel cavity. The sonic velocity is calculated from an expression
for an adiabatic, homogeneous two-component gas-liquid mixture which is
based on Eq. 27 in Reference [14-28].
.. math::
:label: 14.2-56
v_{\text{sonic}}^{2} = \gamma_{\text{fi}} P_{\text{fica}} \big/ \left\{ \left\lbrack \alpha_{\text{fica}}^{2} \cdot \rho_{\text{fica}} + \alpha_{\text{fica}} \left( 1 - \alpha_{\text{fica}} \right) \rho_{\text{fuca}} \right\rbrack \right. \\ \left.
+ \left\lbrack \left( 1 - \alpha_{\text{fica}} \right)^{2} \rho_{\text{fica}} + \alpha_{\text{fica}} \left( 1 - \alpha_{\text{fica}} \right) \rho_{\text{fica}} \right\rbrack \gamma_{\text{fi}} P_{\text{fica}} K_{\text{fu}} \right\} \\
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 PLUTO2)
: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 :math:`\text{PLUTO2}` 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.