3.9. Steady-State and Single-Phase Transient Hydraulics
3.9.1. Introduction
The core assembly hydraulics treatment in SAS4A/SASSYS‑1 includes the calculation of coolant flow rates and pressure distributions within each core channel. Coolant flow rates and pressures are calculated in a number of different places in the code. They are used in the steady-state thermal hydraulics initialization; the pre-voiding thermal hydraulics module calculates coolant flows and pressures; and after the onset of voiding, they are calculated in the boiling module. This section describes only the steady-state and pre-voiding calculations. Chapter 12 describes the hydraulics calculations after the onset of voiding. Chapter 14 and Chapter 16 describe the hydraulics calculations after pin disruption.
The coolant flow provides the heat removal from the fuel pins, and so the main reason for calculating coolant flow rates before the onset of voiding is to provide information for the fuel-pin temperature calculations. Core channel flow rates are also needed for the PRIMAR-4 primary loop thermal hydraulics calculations. The pressure distribution within a channel is calculated mainly to obtain the pressure-dependent coolant saturation temperature used to determine the onset of boiling.
For the steady-state initialization the user specifies the initial coolant flow rate for each channel and the outlet plenum pressure. The code then calculates the pressure distribution in each channel, starting from the outlet and working down to the inlet. The channel with the largest steady-state pressure drop is used to determine the steady-state inlet plenum pressure; and the inlet orifice coefficients in all other channels are adjusted so that all channels have the same total steady-state pressure drop.
In the transient hydraulic calculations, both the coolant flow rates and the pressure distributions are calculated for each channel. Pressure boundary conditions are used for the transient coolant flow rate calculations. The driving pressures for these calculations are the inlet and outlet coolant plenum pressures supplied by the PRIMAR-4 module. Since all channels use the same inlet and outlet pressures, flow redistribution between channels is automatically accounted for as temperatures and flows change. Extrapolated coolant temperatures are used to evaluate coolant properties in the transient hydraulics calculations, since coolant flow rates for a time step are calculated before temperatures are calculated.
Figure 3.1.2 shows the logic flow in subroutine TSCL0, the driver for pre-voiding core channel thermal hydraulics. As indicated previously, a multi-level time step approach is used for the transient calculations in SAS4A/SASSYS‑1. In general, the routines in a transient module start with all quantities knows at the beginning of a time step, and a pass through a module results in calculating the values of relevant parameters of the module for one channel at the end of the time step. Different time-step sizes can be used for different phenomena. In particular, the pre-voiding coolant hydraulics time step can be smaller than, but not larger than, the heat-transfer time step or the PRIMAR time step. One pass through TSCL0 calculates one coolant time step for one channel. If the coolant time step finishes a heat-transfer step, then the temperature calculations for the heat-transfer step are also carried out.
Figure 3.9.1 shows the logic flow in subroutine TSCNV1. This subroutine and the routines called by it carry out the pre-viding coolant flow and pressure calculations. Subroutine TSCNV1 also checks on whether to switch to the boiling model.
As indicated in Figure 3.9.2, the pre-voiding transient hydraulics module interacts with PRIMAR and TSHTRN. PRIMAR supplies the inlet and outlet pressures that drive the coolant flows in the core channels. In return, TSCL0 supplies PRIMAR with the current channel flow rates, as well as hydraulic parameters to use in estimating future flow rates for the next PRIMAR time step. TSCL0 supplies the coolant flow rates to TSHTRN, and TSHTRN computes the coolant temperatures used by TSCL0.
Figure 3.9.1 Subroutine TSCNV1, Pre-Boiling Coolant Flow Rates and Pressure Distribution
Figure 3.9.2 Interactions Between Pre-Voiding Transient Hydraulics and Other Modules
3.9.2. Basic Equations
Before the onset of voiding, the coolant is treated as incompressible, and the basic equation used for liquid coolant flow in non-voided channels is
with
\(w\) = \(w(t)\) - independent of \(z\)
\(w\) = coolant flow rate (kg/s) \(= \rho_{\text{c}} v A_{\text{c}}\)
\(A_{\text{c}}\) = coolant flow area
\(p\) = pressure
\(\rho_{\text{c}}\) = density
\(v\) = coolant velocity (m/s)
\(z\) = axial position
\(t\) = time
\(f\) = friction factor, which is approximated differently depending on the chosen model option.
When IFRFAC = 0:
When IFRFAC = 1:
When IFRFAC = 2:
where \(A_{\text{fr}}\) (AFR), \(b_{\text{fr}}\) (BFR), \(C_{\text{1}}\) (C1TRAN),
\(C_{\text{2}}\) (C2TRAN), \(C_{\text{3}}\) (C3TRAN), and \(A_{\text{fL}}\) (AFLAM) are
user-supplied correlation coefficients and the correlations are in the form given by [3-16].
\(\text{Re} = \frac{D_{\text{h}} \left| w \right|}{\mu A_{\text{c}}}\) = Reynolds number
\(Re_{\text{L}}\) = Maximum Reynolds number for flow in the laminar regime (RELAM)
\(Re_{\text{TR}}\) = Minimum Reynolds number for flow in the turbulent regime (RETRAN)
\(\mu\) = viscosity
\(D_{\text{h}}\) = hydraulic diameter
\(\left( \frac{\partial \text{p}}{\partial \text{z}} \right)_{\text{K}}\) = orifice pressure drop
\(g\) = acceleration of gravity
\(\rho\) = density
The axial node structure shown in Figure 3.2.3 is used for the coolant. Figure 3.9.3 shows which variables are defined at node boundaries and which are averages or integrals over the length of a node. With this mesh structure, integrating Eq. (3.9-1) over the length of the channel gives
where
and
\(\left( \frac{\Delta z_{\text{i}}}{A} \right)_{\text{b}}\) and \(\left( \frac{\Delta z_{\text{i}}}{A} \right)_{\text{t}}\) are effective inertial terms at the bottom and top of the subassembly.
\(p_{\text{b}}\) = pressure at bottom of channel, and \(p_{\text{t}}\) = pressure at top of channel. \(A'_{fL}\) and \(A'_{fr}\) are defined as:
Figure 3.9.3 Definition of Variables on the Axial Coolant Mesh
3.9.3. Flow Orifices
As mentioned in the section on the SAS4A/SASSYS‑1 channel treatment in
Section 3.2, the channel is divided axially into a number of zones. One
zone represents the fuel-pin region. Other zones represent regions above
and below the fuel pins. At the bottom of each zone, flow direction-dependent
orifice coefficients can be specified by the user. Additionally,
orifice coefficients can be used at the top of the upper zone.
These axial locations are hereafter called zone interfaces.
The orifice coefficients at zone interfaces can represent orifice blocks,
subassembly inlet and exit losses, and the pin support grid at the bottom of the pins.
Users may define two different types of orifice coefficients at zone interfaces,
a constant average coefficient: \(K_{\text{or}} \left( \text{JC} \right)\)
and a Re-dependent interface coefficient, \(K_{i}\left( Re_i \right)\).
For each type of orifice, the user supplies one orifice coefficient for positive flow and
another for negative flow (anisotropic values).
The anisotropic constant average coefficients are defined using XKORI.
The anisotropic Re-dependent coefficients are defined using ZoneLossCoefTableID
and a Zone Interface Loss Coefficient Tables, described in Section 3.9.3.1.
In addition, in the fuel-pin region, constant average orifice coefficients representing equally-spaced grid spacers can be specified.
These orifice coefficients are used as the values of \(K_{\text{or}} \left( \text{JC} \right)\) in Eq. (3.9-18) for the appropriate axial nodes.
Grid spacer loss coefficients are defined using NGRDSP and XKORGD.
Spacer grid loss coefficients are assumed to be the same in the forward and reverse directions.
The effect of each constant average orifice is distributed evenly over an axial node, unlike the Re-dependent interface coefficients which are concentrated at a point.
3.9.3.1. Anisotropic Re-dependent Loss Coefficients
Anisotropic Re-dependent loss coefficients can be defined at zone interfaces using Interface Loss Coefficient Tables. Zone interfaces are defined as the inlet of a channel axial zone and the outlet of a channel. By default, all zone interfaces are assumed to have a loss coefficient of zero. Only zone interfaces identified in an Interface Loss Coefficient Table will overwrite the default interface loss coefficient. The pressure drop at zone interface \(i\) caused by the loss coefficient is defined as
where
\({variable}_i\) is the variable at interface \(i\), and \(P_{i+}\)/\(P_{i-}\) are the pressure downstream and upstream the interface, respectively. In the core channels, the pressure reported at a channel zone interface is the upstream pressure, \(P(JC_i) = P_{i-}\). At the channel outlet, the reported pressure is the downstream pressure, \(P(MZC) = P_{(NZONE+1)+}\).
Two functional forms are currently available to capture the dependence of the loss coefficient on the Reynolds number:
Two sets of coefficients may be defined for a given zone interface; one set for forward flow and one set for reverse flow. Additionally, an interface area, \(A_i\), and hydraulic diameter, \(Dh_i\), may be defined for the evaluation of the interface pressure drop and interface Re. If an interface area and interface hydraulic diameter are not provided or are defined as zero, the area and hydraulic diameter of the node upstream or downstream the interface are used. If the upstream node has a smaller or equal area than the downstream zone, the upstream node’s area and hydraulic diameter are used for the interface. If the upstream node has a bigger area than the downstream zone, the downstream node’s area and hydraulic diameter are used for the interface.
At very low flow rates, the derivative of the loss coefficients may approach infinity. This introduces numerical difficulty when SAS attempts to reverse flow through a channel. In order to improve the robustness of the semi-implicit solution scheme, an optional lower limit has been introduced for Re-dependent loss coefficients, \(ReL\). When the Reynold number drops below \(ReL\), Eq. (3.9-24) is evaluated at \(ReL\). Additionally, the derivative of \(K_i(Re_i)\) is assumed to be zero.
Column Label |
Description |
Range |
|---|---|---|
iZone |
Index of the channel zone where the coefficient will be defined at the inlet. iZone = 1 is the inlet of the channel. iZone = NZONE + 1 is the outlet of the channel. |
1 \(\leq\) iZone \(\leq\) NZONE + 1 |
Kind |
Kind of functional form for calculating the anisotropic loss coefficient. See Eq. (3.9-24) |
1 or 2 |
C1f, C2f, C3f |
Function coefficients for forward flow |
See Eq. (3.9-25) |
C1b, C2b, C3b |
Function coefficients for backward flow |
See Eq. (3.9-25) |
A (optional) |
Area used to evaluate the Reynolds number and minor pressure loss at the interface. Optional area and hydraulic diameter must be provided as a pair. |
\(A \geq 0\) |
Dh (optional) |
Hydraulic Diameter used to evaluate the Reynolds number at the interface. Optional area and hydraulic diameter must be provided as a pair. |
\(Dh \geq 0\) |
ReL (optional) |
Lower Reynolds number limit for evaluating the functional form. For Reynolds numbers below the lower limit, the functional form is evaluated at ReL and the derivative is assumed to be 0.0. |
\(ReL \geq 0\) |
Example Channel Anisotropic Re-dependent Loss coefficient table input:
TABLE <id> CHAN1
iZone Kind C1f C2f C3f C1b C2b C3b A Dh ReL
1 1 0.5 100000 -1 1.5 2 0 2.2 6.5 1.0
3 2 1 10 -0.000001 2 0.0001 0.00001 0.021 0.08 0.0
4 1 1 100000 -1 2.2 20000 -1 0.002 0.003183099 10.0
END
3.9.4. Finite Difference Equations - Coolant Flow Rates
If \(w_{1} = w\left( t \right),\ w_{2} = w\left( t + \Delta t \right),\ \Delta w = w_{2} - w_{1}\), then we approximate \(\Delta w\) as
and
in which \(\theta_2\) is the degree of implicitness. For small time steps a semi-implicit calculation is used, with \(\theta_1\) = \(\theta_2 = .5~\). For a fully implicit calculation for large time steps, \(\theta_1 = 0\) and \(\theta_2 = 1.0~\). The value used for \(\theta_2\) is
where
and \(\tau\) is the time constant for flow rate changes given by
The definition of \(d_1\) is given by Eq. (3.9-37) below. This correlation for \(\theta_2\) is discussed in Section 3.19.1. The additional approximations are made that
These equations are combined to give
where \(p_{\text{t}1}\) and \(p_{\text{t}2}\) are the pressures at the top of the subassembly at the beginning and end of the time step. Similarly, \(p_{\text{b}1}\) and \(p_{\text{b}2}\) are the pressures at the bottom of the subassembly at the beginning and end of the time step, respectively.
Solving for the coolant flow rate change, \(\Delta w\), gives
where
The temperatures used to evaluate \(\rho_{\text{c}}\) and \(\overline{\mu}\) in Eq. (3.9-12), Eq. (3.9-14), Eq. (3.9-16), Eq. (3.9-18) and Eq. (3.9-20) are extrapolated coolant temperatures, extrapolated to the end of the coolant step.
3.9.5. Coolant Pressures
After the flow rate has been calculated, the coolant pressures in a channel at the end of a time step are calculated. First, the pressure at node \(\text{MZC}\), the last node at the top of the subassembly, is calculated as
where
\(z_{\text{plu}}\) = outlet plenum elevation
\(\rho_{\text{cout}}\) = density at outlet
\(p_{\text{x}2}\) = coolant outlet plenum pressure at the end of the time step
Note that \(p_{\text{x}2}\), the outlet plenum pressure supplied by the primary loop module, is defined at an elevation \(z_{\text{plu}}\), the elevation of the upper plenum, whereas \(p\left( \text{MZC} \right)\) is calculated at \(z\left( \text{MZC} \right)\), so the gravity head term occurs in Eq. (3.9-39).
After \(p\left( \text{MZC} \right)\) has been calculated, the other pressures are calculated, starting at node \(\text{MZC}\)-1 working down, using
and
Note \(\text{JC}_i = \text{MZC}-1\) for the channel outlet, \(i=\text{NZONE}+1\).
In Eq. (3.9-38) and Eq. (3.9-40), the value used for \(\frac{\partial w}{\partial t}\) is calculated at the end of the time step for the transient calculation as
The pressure, \(p_{\text{b}2}\) at the bottom of the subassembly at the end of the time step is calculated as
\(z_{\text{pll}}\) = plenum location. Note that \(p_{\text{in}2}\), the inlet plenum pressure supplied by the primary loop module, is defined at an elevation \(z_{\text{pll}}\), the elevation of the lower plenum.
The pressure distribution in the steady state is calculated the same as in the transient, except that \(\frac{\partial w}{\partial t}\) is zero in the steady-state calculation.
At the end of the steady state calculation, SAS will adjust the values of XKORI such that the pressure drop
across each channel is consistent, given the user provided flow rates.
These adjustments are reported in the standard output file.
The magnitude of these adjustments can be limited using ADJKTOL.
If ADJKTOL is set to a positive value and SAS determines an adjustment to XKORI that is larger than the
defined tolerance is required, SAS will report an error and the simulation will be terminated.