Demonstration of Pronghorn’s Subchannel Code Modeling of Liquid-Metal Reactors and Validation in Normal Operation Conditions and Blockage Scenarios

Abstract

Pronghorn-SC is a subchannel code within the Multiphysics Object-Oriented Simulation Environment (MOOSE). Initially designed to simulate flows in water-cooled, square lattice, subchannel assemblies, Pronghorn-SC has been expanded to simulate liquid-metal-cooled flows in triangular lattices, hexagonal subchannel assemblies. For this purpose, the algorithm of Pronghorn-SC was adapted to solve the subchannel equations as they are applicable to a hexagonal wire-wrapped sodium-cooled fast reactor. Cheng–Todreas models for pressure drop and cross-flow models were adopted and a coolant heat conduction term was added. To solve these equations, an improved implicit algorithm was developed robust enough to deal with the numerical issues, associated with low flow and recirculation phenomena. To confirm the prediction capability of Pronghorn-SC, calculations and comparisons with available experimental data of 19- and 37-pin assemblies were performed, as well as other subchannel codes. Finally, a flow blockage modeling feature was added. This capability was validated for both water-cooled square sub-assemblies and sodium-cooled hexagonal sub-assemblies, using experimental data of partially and fully blocked cases.

1. Introduction

Liquid Metal Reactors (LMRs) including Sodium Fast Reactors (SFRs) and Lead Fast Reactors (LFRs) are advanced nuclear systems selected as promising designs by the Generation-IV International Forum [1]. The success of these reactor designs primarily depends on the safety characteristics and efficiency, which is tied to reactor core performance. There is a need to demonstrate that the fuel design can operate reliably at high burnup, high power densities and high coolant temperatures. Some of the limitations such as clad thermal creep, corrosion and clad wastage formation are highly temperature dependent phenomena. Hence, predicting the temperature distribution accurately, during normal operation and postulated transient scenarios, is key to demonstrating the safety of these designs.

Subchannel codes are thermal-hydraulic codes that offer an efficient compromise between Computational Fluid Dynamics (CFD) codes and system codes to simulate a nuclear reactor core. They use a quasi-3D model formulation and allow for a finer mesh than system codes, without the high computational cost of CFD. The degrees of freedom (DOF) of a typical subchannel assembly is in the order of a few hundred thousand, while in the corresponding CFD model it would require a few million. That is why a thermal-hydraulic analysis of a nuclear reactor core, performed using a subchannel code, is more efficient than using a CFD code and offers higher fidelity than using a system code. Moreover, a subchannel model can be coupled to 1D thermal-hydraulic system codes or CFD, as well as fuel performance codes, to perform efficient multi-physics/multi-scale calculations.

Various subchannel codes have been developed or adapted to LMR designs. One of the earliest subchannel analysis codes for SFR wire-wrapped sub-assemblies is SUPERENERGY-2 [2]. MATRA-LMR [3] is an adaptation of MATRA [4] to an SFR. SABRE4 [5] is a subchannel code designed for SFRs. SLTHEN [6] is a steady-state LMR thermal-hydraulic analysis code based on the ENERGY model code. COBRA-WC [7] is derived from COBRA-IV to analyze liquid-metal fast breeder reactor assembly transients. Since most of the nuclear reactors are designed based on the subchannel analysis codes, it was important to develop a subchannel analysis capability within MOOSE.

Pronghorn-SC is a submodule of the coarse-mesh thermal-hydraulics reactor core simulator Pronghorn [8], which is a MOOSE [9]-based code. MOOSE was developed to solve systems of coupled, nonlinear partial differential equations that often arise when simulating nuclear systems. The main advantage of the MOOSE framework is that it is a flexible finite element tool that leverages the multiphysics capabilities inherent in MOOSE. Neutronics, fuel performance and thermal-hydraulics form the primary set of physics that needs to be resolved in a coupled fashion.

Pronghorn-SC’s initial development is documented in [10] and was first developed for water-cooled reactors, with subchannels and bare fuel pins in a square-lattice arrangement. This current work expands Pronghorn-SC to be applicable to liquid-metal-cooled reactors, with subchannels and wire-wrapped fuel rods in triangular-lattice arrangements, such as those encountered in SFRs. Additionally, it incorporates the ability to model blockages. The ability to model blockages is an important feature of Pronghorn-SC since the existence of blockages in subchannels can cause local temperature spikes in the bulk fluid and fuel rods that can possibly endanger the fuel integrity. These temperature peaks are a result of local re-circulation phenomena around the blockage location. For this reason, along with the development of a blockage modeling capability, Pronghorn-SC’s solver had to be improved upon in order to be able to solve for flows where re-circulation/low flow is dominant.

2. Governing Equations

The subchannel thermal-hydraulic analysis is based on the conservation equations of mass, linear momentum and energy on specific control volumes. The control volumes are connected in both axial and radial directions to capture the three-dimensional effect of the reactor core. The control volumes for the case of a triangular lattice arrangement are presented in Figure 1.

Integration over control volume $V_i$ produces the subchannel equations for conservation of mass, energy and linear momentum in the axial direction. Integration over transversal volume $V_{ij}$ produces the conservation equation for linear momentum in the lateral direction. All results presented in this work are obtained using Pronghorn-SC’s collocated formulation (velocities and pressures reside at the same nodes). In addition to the collocated formulation, Pronghorn-SC supports a staggered formulation of mass and momentum in the axial direction, in which the linear momentum equation is integrated over a volume $V_i$, which is axially shifted by half a control volume height.

• Conservation of mass:

  $$\frac{\mathrm{d}{\rho }_i}{\mathrm{d}t}{V}_i+\mathsf{\Delta }\dot{m_i}+\sum _j{w}_{ij}=0,$$

  (1)

  where i is the subchannel index and j the index of the neighbor subchannel. $\mathsf{\Delta }$ refers to the difference between the inlet and outlet of the control volume in the axial direction. $\dot{m_i}[kg/sec]$ is the mass flow rate of subchannel i in the axial direction. $w_{ij}[kg/sec]$ is the diversion crossflow in the lateral direction, from subchannel i to neighboring subchannel j, resulting from local pressure differences between the two.
• Conservation of linear momentum in the axial direction:

  $$\frac{\mathrm{d}{\dot{m}}_i}{\mathrm{d}t}\mathsf{\Delta }Z+\mathsf{\Delta }\left(\frac{{\dot{m_i}}^2}{S_i{\rho }_i}\right)+\sum _j{w}_{ij}{U}^{★}=-S_i\mathsf{\Delta }{P}_i+Frictio{n}_i+Dra{g}_{ij}-g{\rho }_i{S}_i\mathsf{\Delta }Z.$$

  (2)
  On the left-hand side of Equation (2), $\frac{\mathrm{d}{\dot{m}}_i}{{\mathrm{d}}_t}$ is the time derivative term and the change of momentum in the axial direction is $\mathsf{\Delta }\left(\frac{{\dot{m_i}}^2}{S_i^z{\rho }_i}\right)$. The inertia transfer between subchannels due to diversion crossflow is ${\sum }_j{w}_{ij}{U}^{★}$, where $U^{*}$ is the axial velocity of the donor cell. $-g{\rho }_i{S}_i\mathsf{\Delta }z$ represents the gravity force and g is the acceleration of gravity. It is assumed that gravity is the only significant body force in the momentum equation. The donor cell is the cell from which crossflow flows and depends on the sign of $w_{ij}$. If it is positive, the donor cell is i and if it is negative, the donor cell is j. Henceforth, donor cell quantities will be denoted with the star (${}^{*}$) symbol. $Frictio{n}_i$ represents an irreversible pressure drop, caused by interaction of the fluid with the solid structures, including possible local form losses, due to spacers and mixing vanes. $Frictio{n}_i$ is computed from correlations adapted to the particular geometry of the assemblies. $Dra{g}_{ij}$ is caused by viscous stresses within the fluid at the interface between subchannels and is also computed from closure correlations. $\mathsf{\Delta }{P}_i$ is the axial pressure drop associated with volume $V_i$.
• Conservation of linear momentum in the lateral direction:

  $$\frac{\mathrm{d}{w}_{ij}}{\mathrm{d}t}{L}_{ij}+\frac{L_{ij}}{\mathsf{\Delta }Z}\mathsf{\Delta }\left(w_{ij}{\overline{U}}_{ij}\right)=-S_{ij}\mathsf{\Delta }{P}_{ij}+Frictio{n}_{ij},$$

  (3)

  where $S_{ij}=g_{ij}\mathsf{\Delta }Z$ is the surface area between subchannels i, j and $g_{ij}$ is the gap between them. Lateral pressure gradient ($\mathsf{\Delta }{P}_{ij}/L_{ij}$) across the subchannels and forced mixing between subchannels owing to mixing vanes and spacer grids is the driving force behind diversion crossflow $w_{ij}$. $L_{ij}$ is the distance between the centers of subchannels i, j. ${\overline{U}}_{ij}$ is the average axial velocity of the two subchannels. The overall friction loss term $Frictio{n}_{ij}$ models all the viscous effects and form losses associated with the transversal fluid motion through the gap.
• Conservation of enthalpy:

  $$\frac{\mathrm{d}{〈\rho h〉}_i}{{\mathrm{d}}_t}{V}_i+\mathsf{\Delta }\left({\dot{m}}_i{h}_i\right)+\sum _j{w}_{ij}{h}^{★}+h_{ij}^{\prime }=q_i^{\prime }\mathsf{\Delta }Z-\sum _j{Y}_{ij}\frac{S_{ij}{\eta }_{ij}}{L_{ij}}(T_i-T_j)+\frac{\mathsf{\Delta }\left(Y_i{S}_i{T}_i\right)}{\mathsf{\Delta }Z},$$

  (4)

  where $h_{ij}^{\prime }$ is defined as the turbulent enthalpy transfer between subchannels i, j and $q_i^{\prime }$ is the average linear power density $\left[\frac{kW}{m}\right]$ going into the control volume $V_i$ of subchannel i from the neighboring fuel rods. $T_i$ is the average fluid temperature of volume $V_i$, $Y_{ij}\left[\frac{W}{mK}\right]$ is the thermal conductivity of the coolant between the adjacent subchannels $i,j$ and ${\eta }_{ij}$ is the radial conduction factor, yielding the effective conduction length between them. The second term on the right-hand side models the radial heat conduction and the last term on the left-hand side models the axial heat conduction. Both these terms are particularly important when the coolant is liquid metal and are neglected when the coolant is water. Dissipation due to viscous stresses can be neglected and the total derivative of pressure (work of pressure) can be set to zero. Also neglected is the volumetric heat source due to moderation, since heat is mainly transferred to the fluid through the fuel rods’ surface.

2.1. Geometric Parameters

For square lattice arrangements, the definition of the subchannel geometry is straightforward. The reader ought to refer to [11] for details. However, a specific geometric treatment is implemented for triangular lattices. Geometric parameters, such as flow area, wetted perimeter and subchannel equivalent hydraulic diameter, are first defined for bare and then expanded for wire-wrapped fuel bundles. Subindices 1, 2 and 3 correspond to inner, edge and corner subchannel types, respectively. P is the pitch, W is the gap between fuel rods and duct wall, D is the fuel rod diameter, $D_w$ is the wire diameter, $\theta$ is the wire angle, H is the wire lead length, A is the flow area, $P_w$ is the wetted perimeter and $D_e$ is the equivalent hydraulic diameter. Prime superscript (${}^{\prime }$) is used to denote bare bundles throughout this section. Using this nomenclature, the geometric parameters are defined as follows:

• Bare rod, flow area and wetted perimeter:

  $$\begin{array}{cc}\hfill & A_1^{\prime }=\frac{\sqrt{3}}{4}{P}^2-\frac{\pi D^2}{8},\hfill \end{array}$$

  (5)

  $$\begin{array}{cc}\hfill & A_2^{\prime }=P(W-\frac{D}{2})-\frac{\pi D^2}{8},\hfill \end{array}$$

  (6)

  $$\begin{array}{cc}\hfill & A_3^{\prime }=\frac{{(W-\frac{D}{2})}^2}{\sqrt{3}}-\frac{\pi D^2}{24},\hfill \end{array}$$

  (7)

  $$\begin{array}{cc}\hfill & P_{w1}^{\prime }=\frac{\pi D}{2},\hfill \end{array}$$

  (8)

  $$\begin{array}{cc}\hfill & P_{w2}^{\prime }=P+\frac{\pi D}{2},\hfill \end{array}$$

  (9)

  $$\begin{array}{cc}\hfill & P_{w3}^{\prime }=\frac{\pi D}{6}+\frac{2(W-\frac{D}{2})}{\sqrt{3}},\hfill \end{array}$$

  (10)
• Wire wrapped rod, flow area and wetted perimeter:

  $$\begin{array}{cc}\hfill & A_1=A_1^{\prime }-\frac{\pi D_w^2}{8cos\left(\theta \right)},\hfill \end{array}$$

  (11)

  $$\begin{array}{cc}\hfill & A_2=A_2^{\prime }-\frac{\pi D_w^2}{8cos\left(\theta \right)},\hfill \end{array}$$

  (12)

  $$\begin{array}{cc}\hfill & A_3=A_3^{\prime }-\frac{\pi D_w^2}{24cos\left(\theta \right)},\hfill \end{array}$$

  (13)

  $$\begin{array}{cc}\hfill & P_{w1}=P_{w1}^{\prime }-\frac{\pi D_w}{2cos\left(\theta \right)},\hfill \end{array}$$

  (14)

  $$\begin{array}{cc}\hfill & P_{w2}=P_{w2}^{\prime }-\frac{\pi D_w}{2cos\left(\theta \right)},\hfill \end{array}$$

  (15)

  $$\begin{array}{cc}\hfill & P_{w3}=P_{w3}^{\prime }-\frac{\pi D_w}{6cos\left(\theta \right)},\hfill \end{array}$$

  (16)

  $$\begin{array}{cc}\hfill & cos\left(\theta \right)=\frac{H}{\sqrt{H^2+{\left(\pi (D+D_w)\right)}^2}},\hfill \end{array}$$

  (17)
• Wire projected area:

  $$\begin{array}{cc}\hfill & A_{r1}=\pi (D+D_w)\frac{D_w}{6},\hfill \end{array}$$

  (18)

  $$\begin{array}{cc}\hfill & A_{r2}=\pi (D+D_w)\frac{D_w}{4},\hfill \end{array}$$

  (19)

  $$\begin{array}{cc}\hfill & A_{r3}=\pi (D+D_w)\frac{D_w}{6},\hfill \end{array}$$

  (20)
• Hydraulic diameter:

  $$\begin{array}{cc}\hfill & D_{ei}^{\prime }=\frac{4A_i^{\prime }}{P_{wi}^{\prime }},\hfill \end{array}$$

  (21)

  $$\begin{array}{cc}\hfill & D_{ei}=\frac{4A_i}{P_{wi}},\hfill \end{array}$$

  (22)

2.2. Closure Models

In computational tools that simulate fluid motion, any phenomenon not resolved by the method itself (i.e., something that has smaller length scales than the grid) must be taken into consideration, thus modeled in some way. CFD discretizes the flow domain in high detail and uses the no-slip boundary condition on the wall surface, such that there is no need to use a friction model. On the other hand, there is a need to model the turbulent effect of the unresolved scales of motion. This is done using Reynolds Averaged Navier–Stokes or Large Eddie Simulation, or turbulent modeling. The subchannel formulation used in Pronghorn-SC considers one relatively large control volume between fuel rods and integrates/averages over that volume. This means that, similarly to CFD, closure models are required to account for the effect of averaging over the unresolved scales of motion. The required closure models in the subchannel formulation compute the axial and lateral friction and form pressure losses, turbulent mixing of momentum and enthalpy in the lateral direction and the conduction factor. Moreover, when structures are included in the subchannel formulation (e.g., fuel pins or ducts in LMFRs), heat exchange coefficient correlations are required to compute the heat transfer between the fluid and structures.

• Axial direction friction term:

  $$Frictio{n}_i=-\frac{1}{2}{K}_i\frac{{\dot{m}}_i\left|{\dot{m}}_i\right|}{S_i{\rho }_i}.$$

  (23)

  where $K_i=[\frac{f_w}{Dh{y}_i}\mathsf{\Delta }Z+k_i]$ is an overall axial loss coefficient, encompassing concentrated form losses $k_i$ and frictional losses $\frac{f_w}{Dh{y}_i}\mathsf{\Delta }Z$ due to the fluid and rod interaction. $S_i$ is the axial flow area, $f_w=4f$ is the Darcy friction factor, f is the Fanning friction factor and $Dh{y}_i$ is the hydraulic diameter.
• Lateral direction friction term:

  $$Frictio{n}_{ij}=-\frac{1}{2}{g}_{ij}\mathsf{\Delta }Z{K}_{ij}{\rho }_{}\left|u_{ij}\right|u_{ij}=-\frac{1}{2}{K}_{ij}\frac{w_{ij}\left|w_{ij}\right|}{S_{ij}{\rho }^{★}}.$$

  (24)

  where $K_{ij}$ is an overall loss coefficient encompassing lateral concentrated form and friction losses, $S_{ij}$ is the lateral flow area between subchannel i and subchannel j and ${\rho }^{*}$ is the donor cell density.
• Friction factor:
  There are three friction models implemented in Pronghorn-SC. For bare square lattices, the default model is a standard flow-dependent friction coefficient formulation [12].

  $$f_i\to \left\{\begin{array}{cc}\frac{1}{64},\hfill & R{e}_i<1\hfill \\ \frac{64}{R{e}_i},\hfill & 1\le R{e}_i<5000\hfill \\ 0.316R{e}_i^{-0.25},\hfill & 5,000\le R{e}_i<30,000\hfill \\ 0.184R{e}_i^{-0.20},\hfill & 30,000\le R{e}_i\hfill \end{array}\right.$$

  (25)
  The other two models [11,13,14] for bare and wire-wrapped assemblies, both square and hexagonal, are presented below. The friction factor $f_i^{\prime }$ for a bare assembly is given by:

  $$\begin{array}{c}\hfill C_{fi}^{\prime }=a+b_1(P/D-1)+b_2{(P/D-1)}^2,\end{array}$$

  (26)

  $$\begin{array}{c}\hfill f_i^{\prime }=\frac{C_{fi}^{\prime }}{R{e}_i^n}.\end{array}$$

  (27)
  Specifically for triangular arrays, the flow regime is modeled as a function of the ratio of pitch over diameter (P/D) for laminar and turbulent flows as follows:

  $$\begin{array}{c}\hfill R{e}_L=320\times {10}^{P/D-1},\end{array}$$

  (28)

  $$\begin{array}{c}\hfill R{e}_T=10,000\times {10}^{0.7(P/D-1)},\end{array}$$

  (29)

  where $R{e}_L$ is the critical Reynolds number for laminar flow and $R{e}_T$ is the critical Reynolds number for turbulent flow. In between the two, a transition regime is obtained. $n=1$ is for laminar flow and $n=0.18$ is for turbulent flow. For the transition regime, the following equation is used:

  $$\begin{array}{c}\hfill f_{iTR}=f_{iL}(1-{\psi }^{\gamma })+f_{iT}{\psi }^{\gamma },\psi =\frac{log\left(Re\right)-logR{e}_L}{log\left(R{e}_T\right)-logR{e}_L},\gamma =\frac{1}{3},\end{array}$$

  (30)
  When Equation (26) is used for edge or corner channels, $P/D$ is replaced by $W/D$. The effect of $P/D$ (or $W/D$) is separated into two regions: $1.0\le P/D\le 1.1$ and $1.1\le P/D\le 1.5$.

  Table 1. Coefficients for bare rod subchannel friction factor constants in triangular array.

  Flow Regime	Subchannel	$1.0\le \mathit{P}/\mathit{D}\le 1.1$			$1.1\le \mathit{P}/\mathit{D}\le 1.5$
  		$\mathit{a}$	${\mathit{b}}_{\mathbf{1}}$	${\mathit{b}}_{\mathbf{2}}$	$\mathit{a}$	${\mathit{b}}_{\mathbf{1}}$	${\mathit{b}}_{\mathbf{2}}$
  Laminar	Interior	26.00	888.2	−3334	62.97	216.9	−190.2
  	Edge	26.18	554.5	−1480	44.40	256.7	−276.6
  	Corner	26.98	1636	−10050	87.26	38.59	−55.12
  Turbulent	Interior	0.09378	1.398	−8.664	0.1458	0.03632	−0.03333
  	Edge	0.09377	0.8732	−3.341	0.1430	0.04199	−0.04428
  	Corner	0.1004	1.625	−11.85	0.1499	0.006706	−0.009567

  and

  Table 2. Coefficients for bare rod subchannel friction factor constants in square array.

  Flow Regime	Subchannel	$1.0\le \mathit{P}/\mathit{D}\le 1.1$			$1.1\le \mathit{P}/\mathit{D}\le 1.5$
  		$\mathit{a}$	${\mathit{b}}_{\mathbf{1}}$	${\mathit{b}}_{\mathbf{2}}$	$\mathit{a}$	${\mathit{b}}_{\mathbf{1}}$	${\mathit{b}}_{\mathbf{2}}$
  Laminar	Interior	26.37	374.2	−493.9	35.55	263.7	−190.2
  	Edge	26.18	554.5	−1480	44.40	256.7	−267.6
  	Corner	28.62	715.9	−2807	58.83	160.7	−203.5
  Turbulent	Interior	0.09423	0.5806	−1.239	0.1339	0.09059	−0.09926
  	Edge	0.09377	0.8732	−3.341	0.1430	0.04199	−0.04428
  	Corner	0.09755	1.127	−6.304	0.1452	0.02681	−0.03411

  present the parameters a, $b_1$ and $b_2$ for the bare hexagonal and square assemblies, respectively. These parameters have been fitted to the results from Rehme [15] by Cheng and Todreas [11] with polynomials for each subchannel type.
  For wire-wrapped rods, the friction factors are derived from their bare-rod counterparts as a function of geometry, wire angle, wire drag and sweeping constants, as follows:

  $$\begin{array}{c}\hfill f_1=\frac{1}{R{e}_1^m}\left[C_{f1}^{\prime }\left(\frac{P_{w1}^{\prime }}{P_{w1}}\right)+W_d\left(\frac{3A_{r1}}{A_1^{\prime }}\right)\left(\frac{D_{e1}}{H}\right){\left(\frac{D_{e1}}{D_w}\right)}^m\right],\end{array}$$

  (31)

  $$\begin{array}{c}\hfill f_2=\frac{C_{f2}^{\prime }}{R{e}_2^m}{\left[1+W_s\left(\frac{A_{r2}}{A_2^{\prime }}\right)ta{n}^2\left(\theta \right)\right]}^{\frac{3-m}{2}},\end{array}$$

  (32)

  $$\begin{array}{c}\hfill f_3=\frac{C_{f3}^{\prime }}{R{e}_3^m}{\left[1+W_s\left(\frac{A_{r3}}{A_3^{\prime }}\right)ta{n}^2\left(\theta \right)\right]}^{\frac{3-m}{2}},\end{array}$$

  (33)

  where m is fitted to the specific geometry, $W_d,W_s$ are the wire drag and sweeping coefficients defined according to the flow regimes as follows:

  $$\begin{array}{c}\hfill W_{dT}=\left[19.56+98.71\left(\frac{D_w}{D}\right)+303.47{\left(\frac{D_w}{D}\right)}^2\right]{\left(\frac{H}{D}\right)}^{-0.541},\end{array}$$

  (34)

  $$\begin{array}{c}\hfill W_{dL}=1.4W_{dT},\end{array}$$

  (35)

  $$\begin{array}{c}\hfill W_{sL}=W_{sT}=-11.0log\left(\frac{H}{D}\right)+19.0.\end{array}$$

  (36)
• Turbulent momentum transfer:

  $$Dra{g}_{ij}=-C_T\sum _j{w}_{ij}^{\prime }\mathsf{\Delta }{U}_{ij}=-C_T\sum _j{w}_{ij}^{\prime }\left[\frac{\dot{m_i}}{{\rho }_i{S}_i}-\frac{\dot{m_j}}{{\rho }_j{S}_j}\right].$$

  (37)

  where $C_T$ is a turbulent modeling parameter, which has been calibrated for bare square lattice assemblies [10].
• Turbulent enthalpy transfer:

  $$h_{ij}^{\prime }=\sum _j{w}_{ij}^{\prime }\mathsf{\Delta }{h}_{ij}=\sum _j{w}_{ij}^{\prime }\left[h_i-h_j\right].$$

  (38)
• Turbulent crossflow:

  $$w_{ij}^{\prime }=\beta S_{ij}\overline{G},\frac{\mathrm{d}{w}_{ij}^{\prime }}{\mathrm{d}z}=\frac{w_{ij}^{\prime }}{\mathsf{\Delta }Z}=\beta g_{ij}\overline{G}.$$

  (39)

  where $\beta$ is the turbulent mixing parameter or thermal diffusion coefficient, $\overline{G}$ is the average mass flux of the adjacent subchannels. $\beta$ is the tuning parameter for the mixing model. Physically, it is a nondimensional coefficient that represents the ratio of the lateral mass flux due to mixing to the axial mass flux. In single-phase flow, no net turbulent mass exchange occurs. Both momentum and energy are exchanged between subchannels and their rates of exchange are characterized in terms of hypothetical turbulent interchange flow rates ($w_{ij}^{{}^{\prime }H},w_{ij}^{{}^{\prime }M}$) [11] for enthalpy and momentum, respectively. Simplifying for single-phase flows, the rates of turbulent exchange for energy and momentum are assumed to be equal:

  $$w_{ij}^{\prime }=w_{ij}^{{}^{\prime }H}=w_{ij}^{{}^{\prime }M}.$$

  (40)
  For bare square assemblies, $\beta$ typically takes a value between 0.05 and 0.15 and can be calibrated for specific geometries. For wire wrapped assemblies, according to the Cheng–Todreas model [16], turbulent mixing is modeled separately for inner and peripheral mixing as a function of geometry and wire angle. The turbulent mixing coefficient between the inner subchannels is given as follows:

  $$\begin{array}{c}\hfill \beta =C_m{\left(\frac{A_{r1}}{A_1}\right)}^{0.5}tan\left(\theta \right),\end{array}$$

  (41)

  and in the corners and edge subchannels the turbulent mixing coefficient is given by:

  $$\begin{array}{c}\hfill \beta =C_s{\left(\frac{A_{r2}}{A_2}\right)}^{0.5}tan\left(\theta \right),\end{array}$$

  (42)

  where the $C_m$ and $C_s$ coefficients are defined as follows for the laminar and turbulent regimes:

  $$\begin{array}{c}\hfill C_{mL}=0.077{\left(\frac{P-D}{D}\right)}^{-0.5},C_{sL}=0.033{\left(\frac{H}{D}\right)}^{0.3},\end{array}$$

  (43)

  $$\begin{array}{c}\hfill C_{mT}=0.14{\left(\frac{P-D}{D}\right)}^{-0.5},C_{sT}=0.75{\left(\frac{H}{D}\right)}^{0.3}.\end{array}$$

  (44)
  For the transition flow regime, a log-blending model is used, which is defined as follows:

  $$\begin{array}{c}\hfill C_{sTR}=C_{sT}+(C_{sT}-C_{sL}){\psi }^{\gamma },\end{array}$$

  (45)

  $$\begin{array}{c}\hfill C_{mTR}=C_{mT}+(C_{mT}-C_{mL}){\psi }^{\gamma },\end{array}$$

  (46)

  $$\begin{array}{c}\hfill \psi =\frac{log\left(Re\right)-log\left(R{e}_L\right)}{log\left(R{e}_T\right)-log\left(R{e}_L\right)},\gamma =\frac{2}{3}.\end{array}$$

  (47)
• Radial conduction factor:
  The radial conduction factor is defined after the Cheng–Todreas model [14] as follows:

  $$\begin{array}{c}\hfill {\eta }_{ij}=\eta =0.66\left(\frac{P_{ij}}{D_i}\right){\left(\frac{S_{ij}}{D_i}\right)}^{-0.3}\end{array}$$

  (48)
• Equations of state:
  These are a set of equations relating the properties of the fluid, such as density and viscosity, to enthalpy, temperature and pressure:

  –

  $T=T(P,h)$.

  –

  $\rho =\rho (T,P)$.

  –

  $\mu =\mu (T,P)$

2.3. Discretization

The collocated discretization of the variables is presented in Figure 2. $i,j$ are the subchannel indexes. $ij$ is the index of the gap between subchannels $i,j$. k is the index in the axial direction.

The governing equations are discretized as follows:

• Conservation of mass:

  $$\dot{m_{i,k}}-\dot{m_{i,k-1}}=-\sum _j{w}_{ij,k}-\frac{{\rho }_{i,k}^{n+1}{V}_{i,k}-{\rho }_{i,k}^n{V}_{i,k}}{\mathsf{\Delta }t}$$

  (49)
  The above equation can be written in matrix form as follows:

  $$\left[\begin{array}{cccc}1& 0& ...& 0\\ -1& 1& ...& 0\\ :& :& ...& :\\ 0& ...& -1& 1\end{array}\right]\times \left[\begin{array}{c}\dot{m_{0,1}}\\ \dot{m_{0,2}}\\ :\\ \dot{m_{i,k-1}}\\ \dot{m_{i,k}}\end{array}\right]=\left[\begin{array}{c}\dot{m_{0,0}}-\sum _j{w}_{0j,1}-\frac{{\rho }_{0,1}^{n+1}{V}_{0,1}-{\rho }_{0,1}^n{V}_{0,1}}{\mathsf{\Delta }t}\\ -\sum _j{w}_{0j,2}-\frac{{\rho }_{0,2}^{n+1}{V}_{0,2}-{\rho }_{0,2}^n{V}_{0,2}}{\mathsf{\Delta }t}\\ :\\ -\sum _j{w}_{ij,k}-\frac{{\rho }_{i,k}^{n+1}{V}_{i,k}-{\rho }_{i,k}^n{V}_{i,k}}{\mathsf{\Delta }t}\end{array}\right]$$

  (50)

  which is equivalent to:

  $${\mathit{M}}_{\mathit{mm}}\overrightarrow{\dot{m}}=\overrightarrow{b_m}-{\mathit{M}}_{\mathit{mw}}\overrightarrow{w}$$

  (51)
• Conservation of linear momentum in the axial direction:

  $$\begin{array}{c}\mathsf{\Delta }{P}_{i,k}=P_{i,k-1}-P_{i,k}=\frac{1}{S_{i,k-1}}[\frac{{\dot{m}}_{i,k}^{n+1}-{\dot{m}}_{i,k}^n}{\mathsf{\Delta }t}\mathsf{\Delta }Z+\hfill \\ \frac{{\dot{m}}_{i,k}^2}{S_{i,k}{\rho }_{i,k}}-\frac{{\dot{m}}_{i,k-1}^2}{S_{i,k-1}{\rho }_{i,k-1}}\\ +\sum _j{w}_{ij,k}{U}^{★}+C_T\sum _j{w}_{ij,k}^{\prime }\left[\frac{\dot{m_{i,k}}}{{\rho }_{i,k-1}{S}_{i,k}}-\frac{\dot{m_{j,k}}}{{\rho }_{j,k-1}{S}_{j,k}}\right]\\ \hfill +\frac{1}{2}{K}_i\frac{{\dot{m}}_{i,k}\left|{\dot{m}}_{i,k}\right|}{S_{i,k}{\rho }_{i,k}}+-g{\rho }_{i,k}{S}_{i,k}\mathsf{\Delta }Z]\end{array}$$

  (52)
  The above equation, can be written in matrix form as follows:

  $$\begin{array}{c}\hfill {\mathit{M}}_{\mathit{pm}}(\overrightarrow{w},\overrightarrow{\dot{m}})\overrightarrow{\dot{m}}=\mathit{S}\overrightarrow{\mathsf{\Delta }P}+\overrightarrow{b_P}\end{array}$$

  (53)

  $$\begin{array}{c}\hfill \mathit{S}\overrightarrow{\mathsf{\Delta }P}=-{\mathit{M}}_{\mathit{pp}}\overrightarrow{P}\end{array}$$

  (54)

  where the matrix ${\mathit{M}}_{\mathit{pm}}$ is calculated using the lagged values of the unknown variables $\overrightarrow{w},\overrightarrow{\dot{m}}$.
• Conservation of linear momentum in the lateral direction:

  $$\begin{array}{c}\hfill 2S_{ij,k}{L}_{ij}{\rho }^{*}\frac{w_{ij,k}^{n+1}-w_{ij,k}^n}{\mathsf{\Delta }t}+\frac{S_{ij,k}{\rho }^{*}{L}_{ij}}{\mathsf{\Delta }Z}\left(\frac{{\dot{m}}_{i,k}}{S_{i,k-1}{\rho }_{i,k-1}}+\frac{{\dot{m}}_{j,k}}{S_{j,k-1}{\rho }_{j,k-1}}\right)w_{ij,k}\\ \hfill +\frac{S_{ij,k}{\rho }^{*}{L}_{ij}}{\mathsf{\Delta }Z}\left(\frac{{\dot{m}}_{i,k-1}}{S_{i,k-1}{\rho }_{i,k-1}}+\frac{{\dot{m}}_{j,k-1}}{S_{j,k-1}{\rho }_{j,k-1}}\right)w_{ij,k-1}\\ \hfill +K_{ij}{w}_{ij,k}\left|w_{ik,k}\right|-2S_{ij,k}^2{\rho }^{*}\left[P_{i,k-1}-P_{j,k-1}\right]=0\end{array}$$

  (55)
  The above equation can be written in matrix form as follows:

  $$\begin{array}{c}\hfill {\mathit{M}}_{\mathit{wp}}\overrightarrow{P}+{\mathit{M}}_{\mathit{ww}}(\overrightarrow{\dot{m}},\overrightarrow{w})\overrightarrow{w}=\overrightarrow{b_w}\end{array}$$

  (56)

  where the matrix ${\mathit{M}}_{\mathit{ww}}$ is calculated using the lagged values of the unknown variables $\overrightarrow{w},\overrightarrow{\dot{m}}$.
• Conservation of enthalpy:

  $$\begin{array}{c}\hfill \frac{{\rho }_{i,k}^{n+1}{h}_{i,k}^{n+1}-{\rho }_{i,k}^n{h}_{i,k}^n}{\mathsf{\Delta }t}{V}_{i,k}+{\dot{m}}_{i,k}{h}_{i,k}-{\dot{m}}_{i,k-1}{h}_{i,k-1}+\sum _j{w}_{ij,k}{h}_k^{★}\\ \hfill +\sum _j{w}_{ij,k}^{\prime }\left[h_{i,k-1}-h_{j,k-1}\right]={〈q^{\prime }〉}_{i,k}\mathsf{\Delta }{z}_k-\sum _j{Y}_{ij,k}\frac{S_{ij,k}{\eta }_{ij,k}}{L_{ij,k}}(T_{i,k}-T_{j,k})\\ \hfill +\frac{Y_{i,k}{S}_{i,k}{T}_{i,k}-Y_{i,k-1}{S}_{i,k-1}{T}_{i,k-1}}{\mathsf{\Delta }{z}_k}.\end{array}$$

  (57)
  The above equation can be written in matrix form as follows:

  $$\begin{array}{c}\hfill {\mathit{M}}_{\mathit{hh}}(\overrightarrow{\dot{m}},\overrightarrow{w})\overrightarrow{h}=\overrightarrow{b_h}\end{array}$$

  (58)

where the matrix ${\mathit{M}}_{\mathit{hh}}$ is calculated using the lagged values of the unknown variables $\overrightarrow{w},\overrightarrow{\dot{m}}$. It is possible to combine all the governing equations in one big system that can be written in matrix form as follows:

$$\left[\begin{array}{cccc}{\mathit{M}}_{\mathit{mm}}& 0& {\mathit{M}}_{\mathit{mw}}& 0\\ {\mathit{M}}_{\mathit{pm}}& {\mathit{M}}_{\mathit{pp}}& 0& 0\\ 0& {\mathit{M}}_{\mathit{wp}}& {\mathit{M}}_{\mathit{ww}}& 0\\ 0& 0& 0& {\mathit{M}}_{\mathit{hh}}\end{array}\right]\times \left[\begin{array}{c}\overrightarrow{\dot{m}}\\ \overrightarrow{P}\\ \overrightarrow{w}\\ \overrightarrow{h}\end{array}\right]=\left[\begin{array}{c}\overrightarrow{b_m}\\ \overrightarrow{b_p}\\ \overrightarrow{b_w}\\ \overrightarrow{b_h}\end{array}\right]$$

(59)

2.4. Interpolation Schemes

The standard interpolation scheme for a dummy variable $\varphi$ in the axial direction reads as follows:

$${\varphi }_{i+1/2}={\alpha }_{i+1/2}{\varphi }_i+(1-{\alpha }_{i+1/2}){\varphi }_{i+1}$$

(60)

where $i+1/2$ denotes the position at the cell interface, i is the bottom cell, $i+1$ is the top cell and ${\alpha }_{i+1/2}$ is the interpolation coefficient. Pronghorn-SC supports four different interpolation schemes:

• Standard upwind:

  $${\alpha }_{i+1/2}=1.0,$$

  (61)
• Standard downwind:

  $${\alpha }_{i+1/2}=0.0,$$

  (62)
• Central difference:

  $${\alpha }_{i+1/2}=\frac{1}{2},$$

  (63)
• Exponential scheme:

  $${\alpha }_{i+1/2}=\frac{(Pe-1)e^{Pe}+1}{Pe(e^{Pe}-1)},$$

  (64)

where Pe is the Peclet number in the axial direction, defined as $Pe=\frac{\mathsf{\Delta }\left(\frac{{\dot{m_i}}^2}{S_i^z{\rho }_i}\right)}{Frictio{n}_i}$.

3. Algorithm

For the purposes of implementing a single-phase subchannel code within MOOSE, a hybrid numerical method of solving the subchannel equations has been developed [10]. The hybrid scheme utilized in this algorithm allows the division of the assembly into $NB$ blocks of $NC\le N$ axial cells per block, as shown in Figure 3. Note that $NB\times NC=N$. The blocks are solved from bottom to top. If $NC=1$ and $NB=N$, the standard, fast-running marching scheme is obtained. At the other extreme, if $NC=N$ and $NB=1$, a monolithic solve for the whole assembly is obtained, which is suggested for highly recirculating flows. In between these two limits, one obtains relatively fast-running schemes that are able to deal with recirculation phenomena, provided they are limited within the block size. The recirculation size would have to be determined by the user beforehand in order to optimize the block size and block number. In general though, Pronghorn-SC is fast enough, such that in most cases, one block performs adequately. The parameters $NB$ and N are defined by the user in the input file. Then, $NC$ is automatically computed, adapting the number of cells so that the total number is N.

The core of the solution method anchors on the calculation of crossflows ${\overrightarrow{w}}_{ij}$. When the crossflows are known, the axial mass flow rates and pressure drops are directly defined by Equations (1) and (2). Hence, the focus of Pronghorn-SC is to first solve for each $w_{ij}$. The residual function used for this nonlinear solve is based on the lateral momentum equation, Equation (65). A Jacobian Free Newton–Krylov type Method is used to solve this equation. The workhorse of this solve is the nonlinear equation solvers found in the Portable, Extensible Toolkit for Scientific Computation [17].

$$f\left(w_{ij}\right)=\frac{\mathrm{d}{w}_{ij}}{\mathrm{d}t}{L}_{ij}+\frac{L_{ij}}{\mathsf{\Delta }z}\mathsf{\Delta }\left(w_{ij}\overline{U}\right)-S_{ij}\mathsf{\Delta }{P}_{ij}+\frac{1}{2}{K}_{ij}\frac{w_{ij}\left|w_{ij}\right|}{{\rho }^{*}}=0.$$

(65)

The residual function constructed via Pronghorn-SC calculates the nonlinear residual $f\left(w_{ij}\right)$ while also updating the other main flow variables: mass flow ${\dot{m}}_i$, turbulent crossflow $w_{ij}^{\prime }$, pressure drop $\mathsf{\Delta }{P}_i$ and pressure $P_i$ through Equations (1), (39), and (2), using the current $w_{ij}$ as needed. This means that the flow variables are tightly coupled with the solution of the crossflows. Similarly as in COBRA [18], the variable $P_i$ is defined as the local pressure minus the exit pressure, $P_i\left(z\right)-P_{exit}$, so at the exit $P_i$ is zero. Each block is solved sequentially from inlet to outlet. The mass flow at the outlet of the previous block and the pressure at the inlet of the next block provide the needed boundary conditions. Each block solution calculates all the unknowns that exist in that block.

The hybrid solution algorithm is shown in Figure 4. Once the main flow variables converge in a block, the enthalpy conservation equation, Equation (4), is solved and enthalpy $\left(h\right)$ is retrieved. In the case where no heat is added to the fluid, the user has the option to either define power as zero or not perform the enthalpy calculation. Using enthalpy, pressure and the equations of state, the temperature $T_i$ and fluid properties, such as density ${\rho }_i$ and viscosity ${\mu }_i$, are calculated. After the fluid properties are updated, the procedure is repeated until the temperature solution has converged. Once the temperature solution converges, the solution moves on to the next block. Once the temperature solution converges in all blocks, pressure is checked for convergence. If there is no convergence, the procedure is repeated, starting again from the first block, until pressure converges.

Algorithm Variations

There are three variations of the algorithm presented above: explicit (default), implicit segregated and implicit monolithic. The last two implicit methods are novel additions to Pronghorn-SC pertaining to this current work, specifically developed to make the code more robust and allow the calculation of recirculation flows.

• Explicit.
  This is the default algorithm, where the unknown flow variables are calculated in an explicit manner through their governing equations. In the explicit version, Equations (49), (52), (55) and (57) are used to solve for the flow variables: $\dot{m_{i,k}},\mathsf{\Delta }{P}_{i,k},P_{i,k-1},w_{ij,k},h_{i,k}$.
• Implicit Segregated.
  In this case, the governing equations are recast in matrix form and the flow variables are calculated by solving the corresponding system. Otherwise, the solution algorithm is the same as in the default method. In the implicit segregated version, Equations (51), (53), (54), (56) and (58) are used to solve for the flow variables: $\overrightarrow{\dot{m}},\overrightarrow{P},\overrightarrow{w_{ij}},\overrightarrow{h}$.
• Implicit Monolithic.
  In this case, the conservation equations are recast in matrix form and combined into a single system, Equation (59). The user can decide whether or not they will include the enthalpy calculation in the matrix formulation. The flow variables are calculated by solving that big system to retrieve: $\overrightarrow{\dot{m}},\overrightarrow{P},\overrightarrow{w_{ij}},\overrightarrow{h}$. The solution algorithm is the same as in the default method, but the solver used in this version is a fixed point iteration instead of a Newton method.

4. Validation

Pronghorn-SC has been validated in multiple benchmarks for water-cooled assemblies and one-directional flows [10]. The effort detailed in this paper aims to expand on this validation work to include sodium-cooled hexagonal assemblies in normal operating conditions and also blockage scenarios for sodium-cooled hexagonal assemblies and water-cooled square assemblies. In summary, this section presents a validation of Pronghorn-SC on four test cases: (i) validation on the Fuel Failure Mockup (FFM) experiment built at Oak Ridge National Laboratory (ORNL) (called ORNL’s 19-pin benchmark henceforward), (ii) validation on Toshiba’s 39-pin benchmark, (iii) validation on the $7\times 7$ sleeve blockage assembly built at Pacific Northwest Laboratory (PNL) and (iv) validation at ORNL’s Thermal-Hydraulic Out-of-Reactor Safety (THORS) facility in a 19-pin sodium-cooled bundle with a central blockage of six channels.

4.1. Oak Ridge National Laboratory’s 19-pin Benchmark

The FFM experiment was built at ORNL for studying the thermal-hydraulic flow characteristics in SFR assemblies [19]. The Test Series 2 results were used to validate Pronghorn-SC. In this series, the fuel rods of the FFM experiment were heated by 19 identical electric cartridges. The configuration and external heat fluxes of these cartridges match typical values expected for SFRs. The measurements in Test Series 2 focused on the distribution of temperatures at the exit of the fuel assembly, duct walls and rod bundle. The nature of the heat source and the lack of neighboring assemblies make the distribution of temperatures at the rod bundle and duct atypical for LMFRs. In contrast, the distribution of temperatures at the exit of the fuel assembly is indicative of the heating of the coolant and flow mixing in the fuel bundle, which is expected to be representative of an actual SFR. Therefore, the validation work focused on the distribution of temperatures at the exit of the fuel assembly.

The design parameters for the Test Series 2 of the FFM experiment are presented in

Table 3. Design and operational parameters for ORNL’s 19-pin benchmark.

Experiment Parameter (Unit)	Value
Number of pins (—)	19
Rod pitch (cm)	0.726
Rod diameter (cm)	0.584
Wire wrap diameter (cm)	0.142
Wire wrap axial pitch (cm)	30.48
Flat-to-flat duct distance (cm)	3.41
Inlet length (cm)	40.64
Heated length (cm)	53.34
Outlet length (cm)	7.62
Outlet pressure (Pa)	1.01 $\times {10}^5$
Inlet temperature (K)	588.5
Power profile (—)	Uniform

. Pressure is assumed to be constant at the outlet of the assembly and temperature and mass flux are assumed to be constant at the inlet. The fuel bundle is divided between inlet, heated and outlet sections along the rod in increasing elevation. A nonzero linear heat rate is only assigned to the heated part of the rod, while no power is imposed at the inlet and outlet sections.

The numbering of rods and subchannels in the fuel assembly is presented in Figure 5. The left panel of this figure shows the rod position along with the assigned index, while the subchannel centers are indicated with red dots. On the right panel, the rods are removed and the red dots indicate the centers of the subchannels. Due to hexagonal symmetry, the temperature distribution was measured over the subchannels that approximately lie on the diagonal line that connects the opposed vertices in the duct. The orientation of the line connects the southwest vertex to the northeast one. This lines includes Subchannels 37, 36, 20, 10, 4, 1, 14 and 28.

The key phenomena dominating the temperature profile in the outlet of the domain is the competing effect between heat convection and conduction in the coolant. An example of the axial and lateral mass flow rates and the temperature and viscosity fields obtained for a high-flow-rate configuration in ORNL’s 19-pin benchmark is presented in Figure 6. The flow rapidly develops along the assembly. However, there is a significantly larger flow in the outer gaps of the rod bundle. This results in outer rods and channels that are colder than the inner ones. The temperature difference between outer and inner subchannels increases with increasing inlet mass flow rate. There is a small competing effect, due to the viscosity in the center of the channels being smaller due to heating, but this effect is of second order, compared with the flow driven by pressure drop. However, with lower inlet mass flow rates, heat conduction in the sodium coolant starts to dominate over heat convection and so the temperature profiles become flatter at the exit. In summary, the temperature distribution measured at the assembly outlet can be regulated by the balance between convection and conduction, which is experimentally regulated by changing the axial mass flux and power at the rod bundles.

Several different combinations of power in the fuel rods and flow rates have been tested in the FFM experiments. However, in the single-phase flow-rate cases, the temperature profile at the exit is simply regulated by the physical balance between convection and conduction. Therefore, it is sufficient for this validation exercise to select only three cases: a high flow rate, a low flow rate and a medium flow rate. As the flow rate decreases, conduction dominates convection and the temperature profiles at the exit become more uniform. On the other hand, as the flow rate increases, convection dominates conduction and the temperature profiles at the exit become more pronounced. Details on the selected cases are provided in

Table 4. Validation cases selected in the ORNL benchmark.

Naming	Run ID	Rod Power	Flow Rate	Reynolds Number
		(W/cm)	(m${}^3$/s)
High flow rate	022472-hf	318.2	3.47 $\times {10}^{-3}$	6.72 $\times {10}^4$
Medium flow rate	020372	30.8	3.15 $\times {10}^{-4}$	7.35 $\times {10}^3$
Low flow rate	022472-lf	4.9	4.67 $\times {10}^{-5}$	9.05 $\times {10}^2$

.

Pronghorn-SC is a relatively new effort, compared to other subchannel codes that have been developed, improved and validated over many years. As such, it is appropriate to compare Pronghorn-SC, against the results obtained by the SUBAC [20] and MATRA-LMR [3] codes. These codes were selected because they have previously presented publicly available results for the ORNL 19-pin benchmark and hence, to the authors’ knowledge, are the most developed subchannel codes applicable to LMFRs. This comparison allows researchers to have a metric of the accuracy that can be reasonably expected from Pronghorn-SC for the ORNL 19-pin benchmark and in sodium-cooled, wire-wrapped, hexagonal assemblies in general. The results obtained by the Pronghorn-SC simulation are compared against the experimental measurements the SUBAC and MATRA-LMR codes in Figure 7. As discussed before, the temperature distribution at the exit of the assembly has a parabolic profile, with a temperature peak in the central subchannels because average mass flow is directed in the outer subchannels. For lower mass flow rates, the temperature profile becomes flatter and the relative error between the codes and the experimental result increases, though for every case, Pronghorn-SC’s prediction is closer to the experimental results. The relative error in the interior subchannels is more sensitive to the change in mass flow rate than that in the outer subchannels. This is because the power density in the inner regions is higher than the outer ones, meaning that small changes in the mass flow can lead to big changes in the coolant temperature in that region.

4.2. Toshiba’s 37-pin Benchmark

The ORNL 19-pin benchmark validated the performance of Pronghorn-SC in a small assembly with relatively large mass flow rates. Under these conditions, the effect of radial heat conduction is limited due to the size of the assembly and the high flow rate. In addition, thermal buoyancy has little effect on the velocity profile. The Toshiba 37-pin benchmark extended the validation range of Pronghorn-SC [21]. This benchmark is based on liquid-sodium experiments conducted by the Toshiba Corporation Nuclear Engineering Laboratory in Japan. It consists of a larger assembly than the ORNL 19-pin benchmark, with one more outer ring of heated rods. The specific power per rod is smaller than the ORNL 19-pin benchmark, but the rods have a slightly larger diameter and a larger axial heated length, which increases the influence of thermal buoyancy on the flow profile.

The characteristics of Toshiba’s benchmark are provided in

Table 5. Design and operational parameters for Toshiba’s 37-pin benchmark.

Experiment Parameter (Unit)	Value
Number of pins (—)	37
Rod pitch (cm)	0.787
Rod diameter (cm)	0.650
Wire wrap diameter (cm)	0.132
Wire wrap axial pitch (cm)	30.70
Flat-to-flat duct distance (cm)	5.04
Inlet length (cm)	—
Heated length (cm)	93.0
Outlet length (cm)	—
Outlet pressure (Pa)	1.01 $\times {10}^5$
Inlet temperature (K)	Experiment dependent
Power profile (—)	Chopped cosine (peaking factor 1.21)

. As in the FFM experiment, the rods are electrically heated. However, contrary to FFM, the resistances in the electrically heated rods are adapted to reproduce a chopped cosine power distribution in the axial direction. All heating rods are assumed to have the same power distribution. The cross section of the fuel assembly is presented in Figure 8. As in the FFM experiments, the quantity of interest is the temperature distribution at the assembly outlet. Due to symmetry, it is enough to analyze the temperature distributions over a symmetry line. This line involves, from south to north, Subchannels 72, 49, 32, 20, 10, 4, 3, 2, 1, 7, 14, 26, 39, and 58.

Similar to the previous 19-pin benchmark, three flow configurations are selected for the validation exercise, which are described in

Table 6. Validation cases selected in Toshiba’s 37-pin benchmark.

Naming	Run ID	Rod Power	Inlet	Flow Rate	Reynolds
		(W/cm)	Temperature (K)	(m${}^3$/s)	Number
High flow rate	B37P02	15.57	484.3	1.48 $\times {10}^{-3}$	1.12 $\times {10}^4$
Medium flow rate	0C37P06	11.92	476.5	3.34 $\times {10}^{-4}$	2.81 $\times {10}^3$
Low flow rate	E37P13	3.89	479.4	1.07 $\times {10}^{-4}$	7.39 $\times {10}^2$

. Note that the high flow-rate case presents a significantly smaller flow rate than the ORNL 19-pin benchmark. Therefore, the temperature profiles for the high flow-rate case will be flatter in the present benchmark.

An example flow distribution for the high flow-rate case is depicted in Figure 9. As expected, a flatter temperature profile is obtained in the bulk of the fuel assembly, when compared to the ORNL 19-pin case. However, in this experiment, the ratio of gap distance to pitch is larger than the 19-pin benchmark case. This produces a significantly larger mass flow in the outer subchannels, as observed in Figure 9a. As a result, the outer subchannels are significantly colder than the center ones. Thus, the expected temperature distribution is a flat distribution in the central region of the assembly, with sharp drops next to the wrapper.

It can be observed in Figure 9a,b that, due to the significant difference between inlet mass flow rates at the outer and center subchannels, there is a considerable flow development length at the entry of the fuel assembly. Inlet velocity conditions were unclear in the experiment report [21], so a uniform mass flow at the inlet was assumed. If the assumption of uniform inlet flow rates turns out to be incorrect, a small deterioration of accuracy of the predicted outlet temperature can be expected. However, the flow field fully develops before the outlet of the assembly, which suggests that a possible error in the inlet conditions will have little effect over the temperature distribution at the outlet of the fuel assembly.

The results obtained for the high, medium and low flow-rate validation cases are presented in Figure 10. As in ORNL’s 19-pin case, Pronghorn-SC is compared with the SUBAC code [20]. SUBAC is, to the authors’ knowledge, the subchannel code for wire-wrapped SFRs, with publicly available results that presented the best agreement to the current benchmark. As observed in Figure 10, for the high mass flow rate case, the Pronghorn-SC predicts results closer to the experimental results than SUBAC. However, when comparing the results predicted for the medium and low flow-rate cases in Figure 10b,c, respectively, Pronghorn-SC over-predicts the temperature distribution. Further analysis determined that the more pronounced distribution of temperatures predicted by Pronghorn-SC towards the center of the assembly may be the result of an overestimation of the momentum mixing rates, which would produce larger than expected flows in the outer channels. As was the case in the previous benchmark, the relative error in the interior subchannels is more sensitive to the change of mass flow rate. This again, is because the power density in the inner regions is higher than the outer ones, meaning that small changes in the mass flow can lead to big changes in the coolant temperature in that region.

4.3. Pacific Northwest Laboratory’s $7\times 7$ Sleeve Blockage Benchmark

PNL’s $7\times 7$ sleeve blockage facility was designed to investigate the turbulent flow phenomena near postulated sleeve blockages in a model nuclear fuel rod bundle. The sleeve blockages were characteristic of fuel-clad “swelling” or “ballooning”, which could occur during loss-of-coolant accidents in pressurized-water reactors [22]. The experimental parameters are presented in

Table 7. Design and operational parameters for PNL’s $7\times 7$ sleeve blockage benchmark.

Experiment Parameter (Unit)	Value
Number of pins (—)	$7\times 7$
Rod pitch (cm)	1.36906
Rod diameter (cm)	0.99568
Length (cm)	144.78
Outlet pressure (Pa)	101325
Inlet temperature (K)	302.594
Reynolds number (—)	$2.9\times {10}^4$
Inlet velocity (m/sec)	1.73736
Power profile (—)	Uniform zero power
Grid spacer location (cm)	40.64, 142.24
Grid spacer loss coefficient (—)	1.14, 1.14
Sleeve blockage location (cm)	64.135 (midway between the spacers)

.

Sleeve blockages (three inches in length) were positioned on the center nine rods of the bundle. Area reductions, of 70 and 90%, were obtained in the center four subchannels of the bundle. The 70 and 90% blockages corresponded to area reductions of 35 and 45% in the subchannels adjacent to the sides of the cluster and 17 and 22% in the subchannels next to the corners of the blockage, respectively. These area reductions were not intended to define those expected during loss-of-coolant accidents but were chosen to provide a severe test case to verify subchannel computer programs. Axial components of local mean velocity and intensity of turbulence were measured, using a one-velocity component 1aser Doppler anemometer.

The 70% and 90% blockage was chosen to validate Pronghorn-SC performance. Pronghorn-SC models the blockage by decreasing the surface area of the affected subchannels accordingly. Since the subchannel formulation is based on the concept of the hydraulic diameter, reducing the surface area affects the system of equations in multiple ways. Most significantly through the $Re$ number and the friction model, pressure drop calculation. Restricting the flow area increases the pressure drop and causes flow to diverge to the adjacent subchannels. Furthermore, the user has the option to apply a concentrated form loss coefficient on the affected subchannels at the corresponding axial cell. This will have an effect similar to area reduction. Pronghorn-SC was run with 28 axial cells for the 70% blockage case and 84 axial cells for the 90% blockage case. A user-set local form loss coefficient at the blockage, $K_{bl}=0.3$ and $K_{bl}=0.9$, was also applied for the two cases, respectively, which was axially distributed among the blocked cells. These values were fitted to produce the best agreement. The default modeling parameters $C_T=2.6,\beta =0.006$ were used. In addition to the subchannel code, a CFD simulation (Star-CCM+) of the experiment was made with 10.5 million cells, for the 70% blockage case. The results presented in Figure 11 and Figure 12 showcase the relative velocity of a center subchannel across the length of the assembly.

Pronghorn-SC utilized the implicit monolithic solver, specifically developed to deal with recirculation, as it was the only one that managed to robustly solve the problem. Predicted subchannel average velocities agreed well with measured values for both cases. Pronghorn-SC’s predicted flow of a central subchannel over-predicts the mixing length downstream of the blockage and is quicker to reduce upstream of the blockage. One possible explanation of this behavior has to do with the nature of Pronghorn-SC’s calculation. Averaged quantities over relatively large volumes are expected to be slower to adapt to local rapid changes. This could also indicate that the inter-channel mixing is underestimated.

For this reason, the subchannel simulation was run again, this time with a larger turbulent mixing parameter of $\beta =0.06,36axialcells,K_{bl}=1$ and the result is presented in Figure 13. The simulation with the adjusted mixing effect agrees much better with the experimental results, especially downstream of the blockage. This suggests that the calibrated default value of $\beta$ is not general enough to adequately model scenarios where a blockage augments the mixing effects in the wake.

At the exit region of the blockage, the experimental velocity profile obtained with the 70% blockage exhibits a jetting characteristic that was not measured in the 90% blockage case. According to the authors of the experimental analysis [22], jetting may not have been detected with the 90% blockage because the measuring volume could not be positioned as close to the blockage axial center line as was possible with the 70% blockage. Though it is also probable that no jetting was present due to flow choking. Pronghorn-SC overestimates the jetting effect in both cases.

It should be noted that the CFD simulation took about 3 h to converge, while Pronghorn-SC took about 3 s. Considering that, along with the agreement of the Pronghorn-SC data with the experimental data, one can say that the Pronghorn-SC is a useful engineering tool for modeling blockage scenarios in square water-cooled assemblies.

4.4. Thermal-Hydraulic Out-of-Reactor Safety Six-Channel Center Blockage Benchmark

THORS bundle 3A also simulates the Fast Flux Test Facility and Clinch River Breeder Reactor configurations. Nineteen electrically heated pins are contained inside a round duct, which has unheated dummy pins along the duct wall. The central six channels (1, 2, 3, 4, 5 and 6) are blocked by a non-heat-generating 6.35-mm-thick stainless-steel plate [23]. The bundle cross section is shown in Figure 14. The circles with the crosses indicate the position of thermocouples at the assembly exit. Pronghorn-SC modeled the THORS bundle 3A blockage with a 90% area reduction on the affected subchannels. The Pronghorn-SC model’s geometry and subchannel and rod index notation is shown in Figure 15. The experimental parameters are presented in

Table 8. Design and operational parameters for THORS six-channel blockage benchmark.

Experiment Parameter (Unit)	Value
Number of pins (—)	19
Rod pitch (cm)	0.726
Rod diameter (cm)	0.584
Wire wrap diameter (cm)	0.142
Wire wrap axial pitch (cm)	30.48
Flat-to-flat duct distance (cm)	3.41
Inlet length (cm)	30.48
Heated length (cm)	53.34
Outlet length (cm)	7.62
Blockage location (cm)	68.58
Outlet pressure (Pa)	2.0 $\times {10}^5$
Inlet temperature (K)	714.261
Inlet flow rate (m${}^3$/s)	3.41 $\times {10}^{-3}$
Power profile (—)	Uniform
Pin power (kW/m)	33

. The Pronghorn-SC model cross section in this benchmark is identical to the one used in Section 4.1. Similar studies have been performed for a wire-wrapped-rod bundles with lead-bismuth eutectic (LBE) coolant [24,25,26].

Run 101 was chosen to validate Pronghorn-SC performance. The THORS experiment measured the temperatures at the exits of selected subchannels. Figure 16 presents the exit temperature distribution, expressed as $T-T_{in}$ for the experimental case of 33 kW/m per pin and 100% flow at 54 gpm (Run 101), along with the Pronghorn-SC prediction. Due to the approximation of the circular experimental test section with a hexagonal Pronghorn-SC model, there is a subchannel index correspondence between the two geometries as follows: 43(37), 42(36), 17(20), 16(10), 3(4), 6(1), 8(14) and 28(28). Where the number outside the parentheses refers to the Pronghorn-SC model and the number inside the parentheses refers to the experimental facility (Figure 14).

Predicted subchannel average temperatures agreed relatively well with measured experimental values for Run 101 with a six-channel center blockage with a bigger error in Subchannel 17(20). Nevertheless, Pronghorn-SC consistently over-predicted the exit temperature for all subchannels, which can be attributed to Pronghorn-SC’s smaller cross-sectional area. It should be noted that Pronghorn-SC calculates surface averages while the experimental results are measured at the subchannel centers. As such, it is expected that Pronghorn-SC results will be a bit higher than the experimental values, since the location of the measurements is away from the heated walls in the center of the subchannels. The discrepancy in Subchannel 17(20) might very well be attributed to the location of the thermocouples and the approximate relationship between the model and actual experiment geometry. For the center subchannels where the Pronghorn-SC model geometry is more representative, the agreement is better. The poorer agreement in the exterior subchannels may be due to steeper temperature gradients in that region since the Pronghorn-SC code calculates average channel temperatures, whereas the thermocouples might be in a subchannel temperature gradient.

A CFD model was developed to further evaluate Pronghorn-SC’s performance. The CFD model had about 1 million cells and utilized an implicit unsteady transient solver. Segregated fluid and energy solvers, $k-\omega$ turbulence modeling and the default polyhedral STAR-CCM+ mesher were used. Figure 17 presents the CFD simulation results on a 2D plane around the blockage location. Furthermore, the axial profiles of massflow and temperature are plotted for a center subchannel along the stream-wise direction in Figure 18. Massflow is forced around the blockage, which causes flow to be reduced in the axial direction. The blockage causes a recirculation region to be formed downstream, which can been seen in CFD results. Due to the axial flow being reduced around the area of the blockage, a temperature peak is observed. On the other hand, downstream of the blockage, recirculation causes a cooling effect, as massflow rushes back in the central region from the outer cooler subchannels, causing the temperature to drop back down. The axial profile of the average temperature of the center subchannel agrees well with the CFD calculation of the center-line temperature. Before the blockage, the average value is a bit higher than the center value. After the blockage, enhanced mixing causes the values to overlap. Using the Pronghorn-SC temperature profile, a broad estimation of the recirculation length can be made by measuring the distance between the end of the heating pick and the end of the blockage. The result is 1.765 inches, which is consistent with the experimentally reported 2 inches [23].

Finally, it should be noted that the effect of the simulated blockage, in both cases presented in this study, depends on the axial discretization, flow area reduction and user-defined local form loss coefficient, along with the turbulent modeling parameters. Due to the nonlinear nature of the friction pressure drop calculation, the effect of these parameters is not straightforward and special care must be taken by the user to properly simulate the blockage effects and produce consistent results.

5. Conclusions

This work details the expansion of the Pronghorn-SC code. This expansion allows the code to calculate the flows of a liquid-metal coolant in hexagonal wire-wrapped assemblies, such as SFRs. Moreover, the expanded Pronghorn-SC code gained the ability to model blockage thanks to a new implicit solver integrated into the code. The improvements to the Pronghorn-SC code can be summarized as follows:

• Added liquid sodium fluid properties.
• Added radial and axial heat conduction modeling to the enthalpy conservation equation for liquid-metal coolant.
• Added effect of the wire-wrap in the pressure drop calculation, turbulent mixing and sweep enthalpy flow.
• Created an implicit solver that is more numerically robust and allows solutions in low flow and re-circulation conditions.
• Added blockage modeling for square and hexagonal assemblies.

Furthermore, the validation work already done for Pronghorn-SC in water-cooled square assemblies, has been expanded to include: sodium-cooled hexagonal assemblies in normal operating conditions and also blockage scenarios for sodium-cooled hexagonal assemblies and water-cooled square assemblies. This extends the applicability of Pronghorn-SC from LWRs to SFRs. This validation work suggests that Pronghorn-SC performs better in high flow conditions as opposed to low flow conditions, as seen in Section 4.1 and Section 4.2. Moreover, enhanced turbulent phenomena, specifically noticeable in the cases of flow blockages, required that modeling parameters had to be tuned beyond their default values, to match the experimental results, seen in Section 4.3.

In conclusion, Pronghorn-SC is a novel subchannel code with a robust numerical solver that can simulate steady-state and transient phenomena in a wide range of problems. The code performs (prediction of temperature field) similarly to other established subchannel codes tailored for hexagonal assemblies in liquid-metal-cooled reactors, such as SUBAC and MATRA-LMR. However, it can be natively coupled to other MOOSE-based applications, such as radiation transport, heat conduction, thermomechanics and computational fluid dynamics, to perform multiphysics simulations of liquid-metal-cooled fuel assemblies. The solver is robust enough to handle nonstandard flow conditions, such as flow reversal and blockage scenarios.

Future work will include adding the ability to model lead-bismuth reactors, improving the numerical robustness of the solvers, developing an adapted user interface to deal with unconventional fuel or control assemblies and additional verification and validation efforts.